旁听了今天的上机课，收获良多。
方阵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的地方就会是非常非常非常小的数，这是由于计算机用浮点数运算得到的效果。
放假。
