旁听了今天的上机课，收获良多。

方阵A求逆，先做LU分解。A的逆等于U的逆乘于L的逆，L的逆就利用下三角矩阵求逆算法进行求解，U的逆可以这样求：先将U转置成下三角矩阵，再像对L求逆一样对U的转置求逆，再将得到的结果转置过来，得到的就是U的逆。
因此，关键是下三角矩阵的求逆。

1.下三角矩阵求逆算法

我利用的公式计算公式如下：

对角元素.png
对角元素以下的元素.png

我的代码如下：

def triInverse(matA):
    '''
    @author：zengwei
    输入：
        matA:一个等待求逆的下三角矩阵,大小为n*n，并且希望里面的元素值是浮点数
    输出：
        matInv:matA的逆矩阵
    '''
    numRows = matA.shape[1]
    matL = matA.copy()
    matInv = np.zeros((numRows,numRows))  
    for row in np.arange(0,numRows):
        matInv[row,row] = 1/matL[row,row]
        for k in np.arange(row-1,-1,-1):
            matInv[row,k] = -(np.dot(matInv[row,k+1:row+1],matL[k+1:row+1,k]))/matL[k,k]
    return matInv

同样，生成一个矩阵来测试一下上面的函数，并得到输出。

import numpy as np
test_mat = np.array([[1,0,0,0],[2,3,0,0],[4,5,6,0],[7,8,9,10]],dtype=float)
'''
test_mat = array([[ 1.,  0.,  0.,  0.],
                  [ 2.,  3.,  0.,  0.],
                  [ 4.,  5.,  6.,  0.],
                  [ 7.,  8.,  9., 10.]])
'''
Inv = triInverse(test_mat)
'''
Inv = array([[ 1.        ,  0.        ,  0.        ,  0.        ],
             [-0.66666667,  0.33333333,  0.        ,  0.        ],
             [-0.11111111, -0.27777778,  0.16666667,  0.        ],
             [-0.06666667, -0.01666667, -0.15      ,  0.1       ]])
'''

显然，我不知道解出来对不对，但是我可以用这个逆矩阵Inv和测试矩阵test_mat相乘看看是不是单位矩阵来判断。

np.dot(test_mat,Inv)
'''
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])
'''

是单位矩阵，说明下三角矩阵求逆成功。接下来，利用上面的函数来进行矩阵的求逆。

2.矩阵求逆

首先，先贴出我LU分解的函数：

def getLandU(A):
    '''
    @author：zengwei
    输入：
    A:系数矩阵；array格式，且数值类型必须为浮点数，否则会出现截断误差。
    输出：
    LU分解中的L和U矩阵
    '''
    matA = A.copy()
    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矩阵

然后我们利用getLandU函数和triInverse函数来写矩阵求解函数。如下：

def matInverse(A):
    '''
    @author：zengwei
    输入：
        A:想要求逆的系数矩阵，n*n，希望里面的数值是浮点数
    输出：
        A的逆矩阵
    '''
    L,U = getLandU(A)             #LU分解
    L_inv = triInverse(L)         #L求逆
    U_inv = triInverse(U.T).T     #U求逆
    
    return np.dot(U_inv,L_inv)

下面，我们生成一个随机矩阵来测试，并验证求得的逆矩阵和原矩阵相乘是否为单位矩阵。

TestMat = np.float64(np.random.randint(0,10,(4,4)))
'''
TestMat = array([[3., 7., 4., 0.],
                 [8., 3., 9., 3.],
                 [5., 4., 2., 4.],
                 [7., 0., 7., 6.]])
'''
TestMatInv = matInverse(TestMat)
isI = np.int64(np.dot(TestMatInv,TestMat))  # 验证
'''
isI = array([[1, 0, 0, 0],
             [0, 1, 0, 0],
             [0, 0, 1, 0],
             [0, 0, 0, 0]], dtype=int64)
'''

可见，逆矩阵乘与原矩阵是单位矩阵，这里用了一下np.int64()函数，是为了让结果好看。否则原本是0的地方就会是非常非常非常小的数，这是由于计算机用浮点数运算得到的效果。

放假。