没有 Numpy 的矩阵求逆
Posted
技术标签:
【中文标题】没有 Numpy 的矩阵求逆【英文标题】：Matrix inversion without Numpy 【发布时间】：2015-11-13 20:21:23 【问题描述】：我想在不使用 numpy.linalg.inv 的情况下反转矩阵。
原因是我正在使用 Numba 来加速代码，但不支持 numpy.linalg.inv，所以我想知道是否可以使用“经典”Python 代码反转矩阵。
使用 numpy.linalg.inv 的示例代码如下所示：
import numpy as np
M = np.array([[1,0,0],[0,1,0],[0,0,1]])
Minv = np.linalg.inv(M)
【问题讨论】：
可能不会。没有python“内置”可以为您执行此操作，并且自己编程矩阵求逆绝非易事（例如，请参阅en.wikipedia.org/wiki/… 以获取可能不完整的方法列表）。我也不知道任何numpy-python 独立线性代数包...
如果您只想反转 3x3 矩阵，您可以查找公式 here。 （您最好指定要反转的矩阵的维度和类型。在您的示例中，您使用最简单的单位矩阵。它们是真实的吗？并且是常规的？）
准确地说是一个4x4实矩阵
【参考方案1】：
这是一个更优雅和可扩展的解决方案，imo。它适用于任何 nxn 矩阵，您可能会发现用于其他方法。请注意，getMatrixInverse(m) 将数组数组作为输入。请随时提出任何问题。
def transposeMatrix(m):
return map(list,zip(*m))
def getMatrixMinor(m,i,j):
return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])]
def getMatrixDeternminant(m):
#base case for 2x2 matrix
if len(m) == 2:
return m[0][0]*m[1][1]-m[0][1]*m[1][0]
determinant = 0
for c in range(len(m)):
determinant += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c))
return determinant
def getMatrixInverse(m):
determinant = getMatrixDeternminant(m)
#special case for 2x2 matrix:
if len(m) == 2:
return [[m[1][1]/determinant, -1*m[0][1]/determinant],
[-1*m[1][0]/determinant, m[0][0]/determinant]]
#find matrix of cofactors
cofactors = []
for r in range(len(m)):
cofactorRow = []
for c in range(len(m)):
minor = getMatrixMinor(m,r,c)
cofactorRow.append(((-1)**(r+c)) * getMatrixDeternminant(minor))
cofactors.append(cofactorRow)
cofactors = transposeMatrix(cofactors)
for r in range(len(cofactors)):
for c in range(len(cofactors)):
cofactors[r][c] = cofactors[r][c]/determinant
return cofactors
【讨论】：
这非常有效。根据要求，应该是公认的答案。唯一需要的小改动是在#base case for 2x2 matrix。您需要显式转换为浮点数。
如果矩阵不是正方形，转置函数会出错，我们可以简单地找到列表的转置：zip(*theArray) 取自：***.com/questions/4937491/matrix-transpose-in-python
@MohanadKaleia 你是对的，谢谢。虽然非方阵没有逆矩阵，但我确实声称我的答案是由可重复使用的部分组成的，所以我已根据您的建议修复了转置函数。
@stackPusher 这太棒了。我希望我可以多次投票
如果你使用的是python3，那么你需要将transposeMatrix定义为list(map(list,zip(*m)))而不是map(list,zip(*m))【参考方案2】：
至少截至 2018 年 7 月 16 日，Numba 具有快速矩阵逆。 （您可以看到它们是如何重载标准 NumPy 逆运算和其他操作here。）
这是我的基准测试结果：
import numpy as np
from scipy import linalg as sla
from scipy import linalg as nla
import numba
def gen_ex(d0):
x = np.random.randn(d0,d0)
return x.T + x
@numba.jit
def inv_nla_jit(A):
return np.linalg.inv(A)
@numba.jit
def inv_sla_jit(A):
return sla.inv(A)
对于小型矩阵，它特别快：
ex1 = gen_ex(4)
%timeit inv_nla_jit(ex1) # NumPy + Numba
%timeit inv_sla_jit(ex1) # SciPy + Numba
%timeit nla.inv(ex1) # NumPy
%timeit sla.inv(ex1) # SciPy
[输出]
2.54 µs ± 467 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
67.3 µs ± 9.18 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
63.5 µs ± 7.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
56.6 µs ± 5.03 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
请注意，加速仅适用于 NumPy 逆，而不适用于 SciPy（如预期的那样）。
稍微大一点的矩阵：
ex2 = gen_ex(40)
%timeit inv_nla_jit(ex2) # NumPy + Numba
%timeit inv_sla_jit(ex2) # SciPy + Numba
%timeit nla.inv(ex2) # NumPy
%timeit sla.inv(ex2) # SciPy
[输出]
131 µs ± 12.9 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
278 µs ± 26.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
231 µs ± 24.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
189 µs ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
所以这里仍有加速，但 SciPy 正在迎头赶上。
【讨论】：
【参考方案3】：这是另一种方式，使用高斯消除代替：
def eliminate(r1, r2, col, target=0):
fac = (r2[col]-target) / r1[col]
for i in range(len(r2)):
r2[i] -= fac * r1[i]
def gauss(a):
for i in range(len(a)):
if a[i][i] == 0:
for j in range(i+1, len(a)):
if a[i][j] != 0:
a[i], a[j] = a[j], a[i]
break
else:
raise ValueError("Matrix is not invertible")
for j in range(i+1, len(a)):
eliminate(a[i], a[j], i)
for i in range(len(a)-1, -1, -1):
for j in range(i-1, -1, -1):
eliminate(a[i], a[j], i)
for i in range(len(a)):
eliminate(a[i], a[i], i, target=1)
return a
def inverse(a):
tmp = [[] for _ in a]
for i,row in enumerate(a):
assert len(row) == len(a)
tmp[i].extend(row + [0]*i + [1] + [0]*(len(a)-i-1))
gauss(tmp)
ret = []
for i in range(len(tmp)):
ret.append(tmp[i][len(tmp[i])//2:])
return ret
【讨论】：
我需要这种技术来求解马尔可夫链。 哈！这也是我做这个的原因 foobar 挑战？ ? 是的，你明白了！【参考方案4】：对于 4 x 4 矩阵，使用数学公式可能就可以了，您可以使用谷歌搜索“4 x 4 矩阵求逆公式”找到该公式。例如这里（我不能保证它的准确性）：
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
一般来说，对一般矩阵求逆并不适合胆小的人。您必须了解所有数学上困难的情况，并知道为什么它们不适用于您的使用，并在向您提供数学上的病态输入时捕获它们（或者返回低准确度的结果或已知的数字垃圾）只要您实际上并没有最终除以零或溢出 MAXFLOAT ，这在您的用例中并不重要......您可能会使用异常处理程序捕获并呈现为“错误：矩阵是奇异的或非常接近”）。
作为程序员，使用数值数学专家编写的库代码通常会更好，除非您愿意花时间了解您正在解决的特定问题的物理和数学性质，并在自己的专家中成为自己的数学专家字段。
【讨论】：
【参考方案5】：没有 numpy 的 3x3 逆矩阵 [python3]
import pprint
def inverse_3X3_matrix():
I_Q_list = [[0, 1, 1],
[2, 3, -1],
[-1, 2, 1]]
det_ = I_Q_list[0][0] * (
(I_Q_list[1][1] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][1])) - \
I_Q_list[0][1] * (
(I_Q_list[1][0] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][0])) + \
I_Q_list[0][2] * (
(I_Q_list[1][0] * I_Q_list[2][1]) - (I_Q_list[1][1] * I_Q_list[2][0]))
co_fctr_1 = [(I_Q_list[1][1] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][1]),
-((I_Q_list[1][0] * I_Q_list[2][2]) - (I_Q_list[1][2] * I_Q_list[2][0])),
(I_Q_list[1][0] * I_Q_list[2][1]) - (I_Q_list[1][1] * I_Q_list[2][0])]
co_fctr_2 = [-((I_Q_list[0][1] * I_Q_list[2][2]) - (I_Q_list[0][2] * I_Q_list[2][1])),
(I_Q_list[0][0] * I_Q_list[2][2]) - (I_Q_list[0][2] * I_Q_list[2][0]),
-((I_Q_list[0][0] * I_Q_list[2][1]) - (I_Q_list[0][1] * I_Q_list[2][0]))]
co_fctr_3 = [(I_Q_list[0][1] * I_Q_list[1][2]) - (I_Q_list[0][2] * I_Q_list[1][1]),
-((I_Q_list[0][0] * I_Q_list[1][2]) - (I_Q_list[0][2] * I_Q_list[1][0])),
(I_Q_list[0][0] * I_Q_list[1][1]) - (I_Q_list[0][1] * I_Q_list[1][0])]
inv_list = [[1 / det_ * (co_fctr_1[0]), 1 / det_ * (co_fctr_2[0]), 1 / det_ * (co_fctr_3[0])],
[1 / det_ * (co_fctr_1[1]), 1 / det_ * (co_fctr_2[1]), 1 / det_ * (co_fctr_3[1])],
[1 / det_ * (co_fctr_1[2]), 1 / det_ * (co_fctr_2[2]), 1 / det_ * (co_fctr_3[2])]]
pprint.pprint(inv_list)
inverse_3X3_matrix()
【讨论】：
【参考方案6】：只需添加所有方法
import math
def getMinorIndex(matrixLocal, x, y):
minor = []
for i in range(3):
minorRow = []
if i == x:
continue
for j in range(3):
if j == y:
continue
minorRow.append(matrixLocal[i][j])
minor.append(minorRow)
return minor
def getDeterminant2By2(matrixLocal):
determinant = matrixLocal[0][0] * matrixLocal[1][1] - matrixLocal[0][1] * matrixLocal[1][0]
return determinant
def getDeterminant(matrixLocal):
determinant = 0
for x in range(3):
t = getDeterminant2By2(getMinorIndex(matrixLocal, 0, x))
e = matrixLocal[0][x]
determinant += (t * e * math.pow(-1, x))
return determinant
def getCofactorMatrix(matrixLocal):
cofactorMatrix = []
for i in range(3):
row = []
for j in range(3):
e = matrixLocal[i][j]
t = getDeterminant2By2(getMinorIndex(matrixLocal, i, j))
row.append(t * math.pow(-1, i + j))
cofactorMatrix.append(row)
return cofactorMatrix
def transpose(matrixLocal):
transposeMatrix = []
for i in range(3):
row = []
for j in range(3):
e = matrixLocal[j][i]
row.append(e)
transposeMatrix.append(row)
return transposeMatrix
def divideMatrix(matrixLocal, divisor):
ansMatrix = []
for i in range(3):
row = []
for j in range(3):
e = matrixLocal[i][j]/divisor
row.append(e)
ansMatrix.append(row)
return ansMatrix
cofactor = getCofactorMatrix(matrix)
adjoint = transpose(cofactor)
det = getDeterminant(matrix)
inverse = divideMatrix(adjoint, det)
inverse
【讨论】：
【参考方案7】：我发现高斯乔丹消除算法在尝试此操作时有很大帮助。如果您要使用给定的矩阵（任何大小，即 5x5），其核心公式为 49 页长。最好用这个。要将矩阵求逆，请将其放置为二维数组，然后运行 Inverse 函数
# Python test Guassion Jordan Elimination
# Inputs are 2D array not matrix
Test_Array = [[3,3,2,1,1],[2,1,3,2,3],[1,3,3,2,2],[2,3,3,1,1],
[3,1,2,1,2]]
# Creating storage & initalizing for augmented matrix
# this is the same as the np.zeros((n,2*n)) function
def nx2n(n_Rows, n_Columns):
Zeros = []
for i in range(n_Rows):
Zeros.append([])
for j in range(n_Columns*2):
Zeros[i].append(0)
return Zeros
# Applying matrix coefficients
def update(inputs, n_Rows, n_Columns, Zero):
for i in range(n_Rows):
for j in range(n_Columns):
Zero[i][j] = inputs[i][j]
return Zero
# Augmenting Identity Matrix of Order n
def identity(n_Rows, n_Columns, Matrix):
for i in range(n_Rows):
for j in range(n_Columns):
if i == j:
Matrix[i][j+n_Columns] = 1
return Matrix
# Applying & implementing the GJE algorithm
def Gussain_Jordan_Elimination(n_Rows, n_Columns, Matrix):
for i in range(n_Rows):
if Matrix[i][i] == 0:
print('error cannot divide by "0"')
for j in range(n_Columns):
if i != j:
ratio = Matrix[j][i]/Matrix[i][i]
for k in range(2*n_Columns):
Matrix[j][k] = Matrix[j][k] - ratio * Matrix[i][k]
return Matrix
# Row Operation to make Principal Diagonal Element to '1'
def row_op(n_Rows, n_Columns, Matrix):
for i in range(n_Rows):
divide = Matrix[i][i]
for j in range(2*n_Columns):
Matrix[i][j] = Matrix[i][j]/divide
return Matrix
# Display Inversed Matix
def Inverse(Matrix):
returnable = []
number_Rows = int(len(Matrix))
number_Columns = int(len(Matrix[0]))
Inversed_Matrix = (row_op(number_Rows, number_Columns,
Gussain_Jordan_Elimination(number_Rows, number_Columns,
identity(number_Rows, number_Columns,
update(Matrix, number_Rows, number_Columns,
nx2n(number_Rows, number_Columns))))))
for i in range(number_Rows):
returnable.append([])
for j in range(number_Columns, 2*number_Columns):
returnable[i].append(Inversed_Matrix[i][j])
return returnable
print(Inverse(Test_Array))
【讨论】：
【参考方案8】：我使用http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html 中的公式编写了对 4x4 矩阵求逆的函数：
import numpy as np
def myInverse(A):
detA = np.linalg.det(A)
b00 = A[1,1]*A[2,2]*A[3,3] + A[1,2]*A[2,3]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] - A[1,3]*A[2,2]*A[3,1]
b01 = A[0,1]*A[2,3]*A[3,2] + A[0,2]*A[2,1]*A[3,3] + A[0,3]*A[2,2]*A[3,1] - A[0,1]*A[2,2]*A[3,3] - A[0,2]*A[2,3]*A[3,1] - A[0,3]*A[2,1]*A[3,2]
b02 = A[0,1]*A[1,2]*A[3,3] + A[0,2]*A[1,3]*A[3,1] + A[0,3]*A[1,1]*A[3,2] - A[0,1]*A[1,3]*A[3,2] - A[0,2]*A[1,1]*A[3,3] - A[0,3]*A[1,2]*A[3,1]
b03 = A[0,1]*A[1,3]*A[2,2] + A[0,2]*A[1,1]*A[2,3] + A[0,3]*A[1,2]*A[2,1] - A[0,1]*A[1,2]*A[2,3] - A[0,2]*A[1,3]*A[2,1] - A[0,3]*A[1,1]*A[2,2]
b10 = A[1,0]*A[2,3]*A[3,2] + A[1,2]*A[2,0]*A[3,3] + A[1,3]*A[2,2]*A[3,0] - A[1,0]*A[2,2]*A[3,3] - A[1,2]*A[2,3]*A[3,0] - A[1,3]*A[2,0]*A[3,2]
b11 = A[0,0]*A[2,2]*A[3,3] + A[0,2]*A[2,3]*A[3,0] + A[0,3]*A[2,0]*A[3,2] - A[0,0]*A[2,3]*A[3,2] - A[0,2]*A[2,0]*A[3,3] - A[0,3]*A[2,2]*A[3,0]
b12 = A[0,0]*A[1,3]*A[3,2] + A[0,2]*A[1,0]*A[3,3] + A[0,3]*A[1,2]*A[3,0] - A[0,0]*A[1,2]*A[3,3] - A[0,2]*A[1,3]*A[3,0] - A[0,3]*A[1,0]*A[3,2]
b13 = A[0,0]*A[1,2]*A[2,3] + A[0,2]*A[1,3]*A[2,0] + A[0,3]*A[1,0]*A[2,2] - A[0,0]*A[1,3]*A[2,2] - A[0,2]*A[1,0]*A[2,3] - A[0,3]*A[1,2]*A[2,0]
b20 = A[1,0]*A[2,1]*A[3,3] + A[1,1]*A[2,3]*A[3,0] + A[1,3]*A[2,0]*A[3,1] - A[1,0]*A[2,3]*A[3,1] - A[1,1]*A[2,0]*A[3,3] - A[1,3]*A[2,1]*A[3,0]
b21 = A[0,0]*A[2,3]*A[3,1] + A[0,1]*A[2,0]*A[3,3] + A[0,3]*A[2,1]*A[3,0] - A[0,0]*A[2,1]*A[3,3] - A[0,1]*A[2,3]*A[3,0] - A[0,3]*A[2,0]*A[3,1]
b22 = A[0,0]*A[1,1]*A[3,3] + A[0,1]*A[1,3]*A[3,0] + A[0,3]*A[1,0]*A[3,1] - A[0,0]*A[1,3]*A[3,1] - A[0,1]*A[1,0]*A[3,3] - A[0,3]*A[1,1]*A[3,0]
b23 = A[0,0]*A[1,3]*A[2,1] + A[0,1]*A[1,0]*A[2,3] + A[0,3]*A[1,1]*A[2,0] - A[0,0]*A[1,1]*A[2,3] - A[0,1]*A[1,3]*A[2,0] - A[0,3]*A[1,0]*A[2,1]
b30 = A[1,0]*A[2,2]*A[3,1] + A[1,1]*A[2,0]*A[3,2] + A[1,2]*A[2,1]*A[3,0] - A[1,0]*A[2,1]*A[3,2] - A[1,1]*A[2,2]*A[3,0] - A[1,2]*A[2,0]*A[3,1]
b31 = A[0,0]*A[2,1]*A[3,2] + A[0,1]*A[2,2]*A[3,0] + A[0,2]*A[2,0]*A[3,1] - A[0,0]*A[2,2]*A[3,1] - A[0,1]*A[2,0]*A[3,2] - A[0,2]*A[2,1]*A[3,0]
b32 = A[0,0]*A[1,2]*A[3,1] + A[0,1]*A[1,0]*A[3,2] + A[0,2]*A[1,1]*A[3,0] - A[0,0]*A[1,1]*A[3,2] - A[0,1]*A[1,2]*A[3,0] - A[0,2]*A[1,0]*A[3,1]
b33 = A[0,0]*A[1,1]*A[2,2] + A[0,1]*A[1,2]*A[2,0] + A[0,2]*A[1,0]*A[2,1] - A[0,0]*A[1,2]*A[2,1] - A[0,1]*A[1,0]*A[2,2] - A[0,2]*A[1,1]*A[2,0]
Ainv = np.array([[b00, b01, b02, b03], [b10, b11, b12, b13], [b20, b21, b22, b23], [b30, b31, b32, b33]]) / detA
return Ainv
【讨论】：
您不想使用np.linalg.inv 但np.linalg.det 可以吗？这是一个非常尴尬的要求......
当然，还需要为行列式计算编写另一个“蛮力”实现。或者只是在 Numba 函数之外计算 det 并将其作为参数传递
@sebastian np.linalg.inv 不准确以上是关于没有 Numpy 的矩阵求逆的主要内容，如果未能解决你的问题，请参考以下文章