Python 没有 Numpy 的矩阵求逆
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/32114054/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Matrix inversion without Numpy
提问by Alessandro Vianello
I want to invert a matrix without using numpy.linalg.inv.
我想在不使用numpy.linalg.inv 的情况下反转矩阵。
The reason is that I am using Numba to speed up the code, but numpy.linalg.inv is not supported, so I am wondering if I can invert a matrix with 'classic' Python code.
原因是我使用 Numba 来加速代码，但不支持 numpy.linalg.inv，所以我想知道我是否可以使用“经典”Python 代码反转矩阵。
With numpy.linalg.invan example code would look like that:
使用numpy.linalg.inv示例代码如下所示：
import numpy as np
M = np.array([[1,0,0],[0,1,0],[0,0,1]])
Minv = np.linalg.inv(M)
采纳答案by Alessandro Vianello
I used the formula from http://cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.htmlto write the function that does the inversion of a 4x4 matrix:
我使用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
回答by nigel222
For a 4 x 4 matrix it's probably just about OK to use the mathematical formula, which you can find using Googling "formula for 4 by 4 matrix inverse". For example here (I can't vouch for its accuracy):
对于 4 x 4 矩阵，使用数学公式可能就可以了，您可以使用谷歌搜索“4 x 4 矩阵求逆公式”找到该公式。例如这里（我不能保证其准确性）：
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
http://www.cg.info.hiroshima-cu.ac.jp/~miyazaki/knowledge/teche23.html
In general inverting a general matrix is not for the faint-hearted. You have to be aware of all the mathematically difficult cases and know why they won't apply to your usage, and catch them when you are supplied with mathematically pathological inputs (that, or return results of low accuracy or numerical garbage in the knowledge that it won't matter in your usage case provided you don't actually end up dividing by zero or overflowing MAXFLOAT ... which you might catch with an exception handler and present as "Error: matrix is singular or very close thereto").
一般来说，反转一般矩阵不适合胆小的人。您必须了解所有数学上的困难情况，并知道为什么它们不适用于您的用法，并在为您提供数学病态输入时捕获它们（或者在知识中返回低准确度或数字垃圾的结果）如果您实际上最终没有除以零或溢出 MAXFLOAT，这在您的用例中无关紧要......您可能会用异常处理程序捕获并显示为“错误：矩阵是奇异的或非常接近”）。
It's generally better as a programmer to use library code written by numerical mathematics experts, unless you are willing to spend time understanding the physical and mathematical nature of the particular problem that you are addressing and become your own mathematics expert in your own specialist field.
作为程序员，使用由数值数学专家编写的库代码通常会更好，除非您愿意花时间了解您正在解决的特定问题的物理和数学性质，并成为您自己专业领域的数学专家。
回答by stackPusher
Here is a more elegant and scalable solution, imo. It'll work for any nxn matrix and you may find use for the other methods. Note that getMatrixInverse(m) takes in an array of arrays as input. Please feel free to ask any questions.
这是一个更优雅和可扩展的解决方案，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
回答by webelo
As of at least July 16, 2018 Numba has a fast matrix inverse. (You can see how they overload the standard NumPy inverse and other operations here.)
至少截至 2018 年 7 月 16 日，Numba 具有快速矩阵求逆。（您可以在此处查看它们如何重载标准 NumPy 逆运算和其他操作。）
Here are the results of my benchmarking:
以下是我的基准测试结果：
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)
For small matrices it is particularly fast:
对于小矩阵，它特别快：
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
[Out]
[出去]
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)
Notice that the speedup only works for NumPy inverse, not SciPy (as expected).
请注意，加速仅适用于 NumPy 逆，而不适用于 SciPy（如预期）。
Slightly larger matrix:
稍微大一点的矩阵：
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
[Out]
[出去]
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)
So there's still a speedup here but SciPy is catching up.
所以这里仍然有加速，但 SciPy 正在迎头赶上。
回答by hLk
Here is another way, using gaussian elimination instead:
这是另一种方法，使用高斯消元法：
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:
print("MATRIX NOT INVERTIBLE")
return -1
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
