你有你的計(jì)劃，世界另有計(jì)劃。
列主元消去法是對高斯順序消去法在主元為0或較小值時(shí)的修正；LU分解本質(zhì)上就是高斯消去法；追趕法適用于三對角矩陣；Cholesky分解是可以在對稱正定矩陣情況下比列主元消去法速度提高大約一倍，數(shù)據(jù)存儲(chǔ)內(nèi)存減少一半。
高斯順序消去和列主元消去法的代碼已經(jīng)寫過。今天寫LU分解、追趕法和Cholesky分解。主要參考今天數(shù)學(xué)課上的PPT課件。
1.LU分解
原理不多說，可以參考相關(guān)資料。直接上代碼：
def getLandU(matA):
'''
輸入：
A:系數(shù)矩陣；array格式，且數(shù)值類型必須為浮點(diǎn)數(shù)，否則會(huì)出現(xiàn)截?cái)嗾`差。
輸出：
LU分解中的L和U矩陣
'''
if matA[0,0]==0:
print('不符合要求！需要換行操作。')
else:
numRows = matA.shape[0]
matL = np.identity(numRows) # 準(zhǔn)備給L矩陣用
for row in np.arange(0,numRows-1): # 前n-1行，row代表第幾行
for k in range(row+1,numRows): # 在第row行的情況下，遍歷剩下的行
matL[k,row] = matA[k,row]/matA[row,row] # 獲得L矩陣
if matA[k,row] != 0:
matA[k,:]=matA[k,:] - (matA[k,row]/matA[row,row])*matA[row,:]
return matL,matA # matL為L矩陣，matA變?yōu)閁矩陣
輸入一個(gè)系數(shù)矩陣A，輸出直接得到L矩陣和U矩陣。拿今天上課中例題2.3.1來測試一下。
A = np.array([[2,1,-2],[4,0,-1],[0,3,2]],dtype=float)
b = np.array([[1],[7],[-1]],dtype=float)
'''
A = array([[ 2. , 1. , -2. ],
[ 0. , -2. , 3. ],
[ 0. , 0. , 6.5]])
'''
L,U = getLandU(A)
得到如下結(jié)果：
'''
L = array([[ 1. , 0. , 0. ],
[ 2. , 1. , 0. ],
[ 0. , -1.5, 1. ]])
U = array([[ 2. , 1. , -2. ],
[ 0. , -2. , 3. ],
[ 0. , 0. , 6.5]])
'''
跟課件上一樣，說明效果不錯(cuò)。下面寫LDU分解，只需要在LU分解的基礎(chǔ)上做一些修改即可。以下是代碼，就比上面多了三四行。
def getL_D_U(matA):
'''
輸入：
A:系數(shù)矩陣；array格式，且數(shù)值類型必須為浮點(diǎn)數(shù)，否則會(huì)出現(xiàn)截?cái)嗾`差。
輸出：
LDU分解中的L、D和U矩陣
'''
if matA[0,0]==0:
print('不符合要求！需要換行操作。')
else:
numRows = matA.shape[0]
matL = np.identity(numRows) # 準(zhǔn)備給L矩陣用
for row in np.arange(0,numRows-1): # 前n-1行，row代表第幾行
for k in range(row+1,numRows): # 在第row行的情況下，遍歷剩下的行
matL[k,row] = matA[k,row]/matA[row,row] # 獲得L矩陣
if matA[k,row] != 0:
matA[k,:]=matA[k,:] - (matA[k,row]/matA[row,row])*matA[row,:]
matD = np.identity(numRows)
for i in range(numRows):
matD[i,i] = matA[i,i]
matA[i,:] = matA[i,:]/matA[i,i]
return matL,matD,matA # matL為L矩陣，matD為D矩陣，matA變?yōu)閁矩陣
同樣測試一下：
A = np.array([[2,1,-2],[4,0,-1],[0,3,2]],dtype=float)
L,D,U = getL_D_U(A)
'''
L = array([[ 1. , 0. , 0. ],
[ 2. , 1. , 0. ],
[ 0. , -1.5, 1. ]])
D = array([[ 2. , 0. , 0. ],
[ 0. , -2. , 0. ],
[ 0. , 0. , 6.5]])
U = array([[ 1. , 0.5, -1. ],
[-0. , 1. , -1.5],
[ 0. , 0. , 1. ]])
'''
可見，結(jié)果是正確的。
2.追趕法
追趕法使用于三對角矩陣。實(shí)際上，對于解三對角矩陣，直接用上面的LU分解再求解已經(jīng)沒什么問題了，那么為什么還要弄一個(gè)追趕法呢？這是因?yàn)槿绻龑蔷仃囀且环N很普遍的矩陣形式，那么就有必要編寫一套專門求解三對角矩陣LU分解的算法，這樣可以提高計(jì)算效率，還可以降低計(jì)算機(jī)內(nèi)存的存儲(chǔ)空間。
追趕法本質(zhì)是LU分解的一種特例，這里利用其定義來編寫函數(shù)。
以下是我的代碼。你會(huì)發(fā)現(xiàn)比上面的LU分解少了一層循環(huán)，速度必然提升。不過代碼里面的索引取值比較繞，可能要多花一點(diǎn)時(shí)間搞明白。
def catchUp(matA):
'''
只適用于三對角矩陣！
輸入：
A:系數(shù)矩陣；array格式，且數(shù)值類型必須為浮點(diǎn)數(shù)，否則會(huì)出現(xiàn)截?cái)嗾`差。
輸出：
LU分解中的L和U矩陣
'''
numRows = matA.shape[0]
matL = np.identity(numRows) # 準(zhǔn)備給L矩陣用
matU = np.zeros((numRows,numRows)) # 準(zhǔn)備給U矩陣用
matU[0,0] = matA[0,0] # 初始
for row in np.arange(0,numRows-1): # 前n-1行，row代表第幾行
matL[row + 1,row] = matA[row + 1,row]/matU[row,row] # L的遞推
matU[row+1,row+1] = matA[row+1,row+1]-matL[row + 1,row]*matA[row,row+1] # U的遞推,對角部分
matU[row,row+1] = matA[row,row+1]
return matL,matU # matL為L矩陣，matD為D矩陣，matA變?yōu)閁矩陣
這里用《數(shù)值方法》51頁的例題2.3.2來做測試。
A = np.array([[2,-1,0,0,0],[-1,2,-1,0,0],[0,-1,2,-1,0],[0,0,-1,2,-1],[0,0,0,-1,2]],dtype=float)
'''
A = array([[ 2., -1., 0., 0., 0.],
[-1., 2., -1., 0., 0.],
[ 0., -1., 2., -1., 0.],
[ 0., 0., -1., 2., -1.],
[ 0., 0., 0., -1., 2.]])
'''
L,U = catchUp(A)
得到的結(jié)果如下：
'''
L = array([[ 1. , 0. , 0. , 0. , 0. ],
[-0.5 , 1. , 0. , 0. , 0. ],
[ 0. , -0.66666667, 1. , 0. , 0. ],
[ 0. , 0. , -0.75 , 1. , 0. ],
[ 0. , 0. , 0. , -0.8 , 1. ]])
U = array([[ 2. , -1. , 0. , 0. , 0. ],
[ 0. , 1.5 , -1. , 0. , 0. ],
[ 0. , 0. , 1.33333333, -1. , 0. ],
[ 0. , 0. , 0. , 1.25 , -1. ],
[ 0. , 0. , 0. , 0. , 1.2 ]])
'''
可以看到，得到的結(jié)果和書上一致。
3.Cholesky分解
Cholesky分解是在矩陣式對稱正定矩陣的情況下，只需要分解出L矩陣就可以實(shí)現(xiàn)，L乘于L的轉(zhuǎn)置即可以等于系數(shù)矩陣A。
這里參考如下的計(jì)算公式：
以下是我的代碼，直接分解出L矩陣。
def Cholesky(matA):
'''
只適對稱正定矩陣！
輸入：
A:系數(shù)矩陣；array格式，且數(shù)值類型必須為浮點(diǎn)數(shù)，否則會(huì)出現(xiàn)截?cái)嗾`差。
輸出：
L矩陣
'''
numRows = matA.shape[0]
matL = np.zeros((numRows,numRows))
for row in np.arange(0,numRows):
matL[row,row] = (matA[row,row] - np.sum([s**2 for s in matL[row,0:row]]))**0.5 # 公式2.3.17
for k in np.arange(row+1,numRows):
matL[k,row] = (matA[k,row] - np.sum(np.dot(matL[k,0:row],matL[k-1,0:row])))/matL[row,row] # 公式2.3.18
return matL # matL所求
用書上例題2.3.3來測試一下。
A = np.array([[4,-1,1],[-1,4.25,2.75],[1,2.75,3.5]],dtype=float)
'''
A = array([[ 4. , -1. , 1. ],
[-1. , 4.25, 2.75],
[ 1. , 2.75, 3.5 ]])
'''
L = Cholesky(A)
得到結(jié)果如下。并且用L乘于L的轉(zhuǎn)置也得到了原來的系數(shù)矩陣A。
'''
L = array([[ 2. , 0. , 0. ],
[-0.5, 2. , 0. ],
[ 0.5, 1.5, 1. ]])
'''
reback = np.dot(L,L.T)
'''
reback = array([[ 4. , -1. , 1. ],
[-1. , 4.25, 2.75],
[ 1. , 2.75, 3.5 ]])
'''
結(jié)果與書上一致。
今日課上所講的三種分解全部寫完。今天所講的三種分解，本質(zhì)都是為了減少運(yùn)算步驟，加快計(jì)算機(jī)運(yùn)行線性方程組求解的速度。如果再考慮數(shù)據(jù)的存儲(chǔ)的話，還可以減少內(nèi)存的占用。
