你有你的计划，世界另有计划。

列主元消去法是对高斯顺序消去法在主元为0或较小值时的修正；LU分解本质上就是高斯消去法；追赶法适用于三对角矩阵；Cholesky分解是可以在对称正定矩阵情况下比列主元消去法速度提高大约一倍，数据存储内存减少一半。
高斯顺序消去和列主元消去法的代码已经写过。今天写LU分解、追赶法和Cholesky分解。主要参考今天数学课上的PPT课件。

1.LU分解

原理不多说，可以参考相关资料。直接上代码：

def getLandU(matA):
    '''
    输入：
    A:系数矩阵；array格式，且数值类型必须为浮点数，否则会出现截断误差。
    输出：
    LU分解中的L和U矩阵
    '''
    if matA[0,0]==0: 
        print('不符合要求！需要换行操作。')
    else:
        numRows = matA.shape[0] 
        matL = np.identity(numRows)         # 准备给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变为U矩阵

输入一个系数矩阵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)

得到如下结果：

'''
L = array([[ 1. ,  0. ,  0. ],
           [ 2. ,  1. ,  0. ],
           [ 0. , -1.5,  1. ]])
U = array([[ 2. ,  1. , -2. ],
           [ 0. , -2. ,  3. ],
           [ 0. ,  0. ,  6.5]])    
'''

跟课件上一样，说明效果不错。下面写LDU分解，只需要在LU分解的基础上做一些修改即可。以下是代码，就比上面多了三四行。

def getL_D_U(matA):
    '''
    输入：
    A:系数矩阵；array格式，且数值类型必须为浮点数，否则会出现截断误差。
    输出：
    LDU分解中的L、D和U矩阵
    '''
    if matA[0,0]==0: 
        print('不符合要求！需要换行操作。')
    else:
        numRows = matA.shape[0] 
        matL = np.identity(numRows)         # 准备给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变为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. ]])
'''

可见，结果是正确的。

2.追赶法

追赶法使用于三对角矩阵。实际上，对于解三对角矩阵，直接用上面的LU分解再求解已经没什么问题了，那么为什么还要弄一个追赶法呢？这是因为如果三对角矩阵是一种很普遍的矩阵形式，那么就有必要编写一套专门求解三对角矩阵LU分解的算法，这样可以提高计算效率，还可以降低计算机内存的存储空间。
追赶法本质是LU分解的一种特例，这里利用其定义来编写函数。

追赶法.png

以下是我的代码。你会发现比上面的LU分解少了一层循环，速度必然提升。不过代码里面的索引取值比较绕，可能要多花一点时间搞明白。

def catchUp(matA):
    '''
    只适用于三对角矩阵！
    输入：
    A:系数矩阵；array格式，且数值类型必须为浮点数，否则会出现截断误差。
    输出：
    LU分解中的L和U矩阵
    '''
    numRows = matA.shape[0] 
    matL = np.identity(numRows)         # 准备给L矩阵用
    matU = np.zeros((numRows,numRows))  # 准备给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变为U矩阵

这里用《数值方法》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)

得到的结果如下：

'''
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       ]])
'''

可以看到，得到的结果和书上一致。

3.Cholesky分解

Cholesky分解是在矩阵式对称正定矩阵的情况下，只需要分解出L矩阵就可以实现，L乘于L的转置即可以等于系数矩阵A。
这里参考如下的计算公式：

Cholesky.png

以下是我的代码，直接分解出L矩阵。

def Cholesky(matA):
    '''
    只适对称正定矩阵！
    输入：
    A:系数矩阵；array格式，且数值类型必须为浮点数，否则会出现截断误差。
    输出：
    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)

得到结果如下。并且用L乘于L的转置也得到了原来的系数矩阵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 ]])
'''

结果与书上一致。

今日课上所讲的三种分解全部写完。今天所讲的三种分解，本质都是为了减少运算步骤，加快计算机运行线性方程组求解的速度。如果再考虑数据的存储的话，还可以减少内存的占用。