python矩阵求逆函数_09-30:Python矩阵求逆

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

方阵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的地方就会是非常非常非常小的数，这是由于计算机用浮点数运算得到的效果。

放假。