通常情况下，我会像下面的示例中那样在for循环中反转一个3x3矩阵数组。不幸的是，for循环速度较慢。有更快、更有效的方法来完成这个任务吗？
import numpy as np
A = np.random.rand(3,3,100)
Ainv = np.zeros_like(A)
for i in range(100):
Ainv[:,:,i] = np.linalg.inv(A[:,:,i])
有没有一种有效的方法使用numpy反转一个矩阵数组？
16
- katrasnikj
9
4个回答
17
发现问题出在 numpy.linalg的代码中。如果你查看numpy.linalg.inv，你会发现它只是一个对 numpy.linalg.solve(A, inv(A.shape[0]) 的调用。这将在每次循环迭代中重新创建恒等矩阵。由于所有的数组大小都相同，这是一种浪费时间的做法。通过预先分配恒等矩阵来跳过此步骤可以节省约20%的时间（快速求逆矩阵）。我的测试表明，预分配数组或从结果列表中分配数组并没有太大的区别。
深入一级后，你会找到对 lapack例程的调用，但被包裹在数个检查的包装器内。如果你剥离这些检查，在循环中直接调用lapack（因为你已经知道矩阵的维度，而且也许知道它是实数而不是复数），事情将运行得更快。（注意我已经扩大了数组的大小）
import numpy as np
A = np.random.rand(1000,3,3)
def slow_inverse(A):
Ainv = np.zeros_like(A)
for i in range(A.shape[0]):
Ainv[i] = np.linalg.inv(A[i])
return Ainv
def fast_inverse(A):
identity = np.identity(A.shape[2], dtype=A.dtype)
Ainv = np.zeros_like(A)
for i in range(A.shape[0]):
Ainv[i] = np.linalg.solve(A[i], identity)
return Ainv
def fast_inverse2(A):
identity = np.identity(A.shape[2], dtype=A.dtype)
return array([np.linalg.solve(x, identity) for x in A])
from numpy.linalg import lapack_lite
lapack_routine = lapack_lite.dgesv
# Looking one step deeper, we see that solve performs many sanity checks.
# Stripping these, we have:
def faster_inverse(A):
b = np.identity(A.shape[2], dtype=A.dtype)
n_eq = A.shape[1]
n_rhs = A.shape[2]
pivots = zeros(n_eq, np.intc)
identity = np.eye(n_eq)
def lapack_inverse(a):
b = np.copy(identity)
pivots = zeros(n_eq, np.intc)
results = lapack_lite.dgesv(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
if results['info'] > 0:
raise LinAlgError('Singular matrix')
return b
return array([lapack_inverse(a) for a in A])
%timeit -n 20 aI11 = slow_inverse(A)
%timeit -n 20 aI12 = fast_inverse(A)
%timeit -n 20 aI13 = fast_inverse2(A)
%timeit -n 20 aI14 = faster_inverse(A)
结果令人印象深刻：
20 loops, best of 3: 45.1 ms per loop
20 loops, best of 3: 38.1 ms per loop
20 loops, best of 3: 38.9 ms per loop
20 loops, best of 3: 13.8 ms per loop
编辑：我没有仔细看solve返回的内容。结果'b'矩阵被覆盖并且最终包含了结果。这段代码现在会给出一致的结果。
- Carl F.
5
NumPy数组必须是连续的并且按特定顺序（'C'或'F'）吗？ - Juh_
非常好。您能用“eig”做同样的事情吗？ - Juh_
@Juh_: 我认为数组的顺序很重要。使用默认值。最近我实现了一个非常类似的东西，我想要最小特征值。而不是在每个数组上计算它，我查找了2x2的解析解并编码了它。这将使事情加快数百到数千倍。 - Carl F.
如果A是正定的，那么不是更快地将A分解（Cholesky）为LL*，然后使用
solve_triangular或直接使用linalg.lapack.clapack.dtrtri来解决L,I，然后再进行一次矩阵乘法即可得到答案吗？ - dashesyscipy中的linalg.inv似乎使用getri。 - dashesy
11
自此问题被提出并回答以来，一些事情已经发生了变化，现在
numpy.linalg.inv 支持多维数组，将它们处理为矩阵的堆栈，其中矩阵索引是最后一个（换句话说，数组的形状为 (...,M,N,N)）。这似乎已经在 numpy 1.8.0 中引入了。毫不奇怪，这是性能方面迄今为止最佳的选择：import numpy as np
A = np.random.rand(3,3,1000)
def slow_inverse(A):
"""Looping solution for comparison"""
Ainv = np.zeros_like(A)
for i in range(A.shape[-1]):
Ainv[...,i] = np.linalg.inv(A[...,i])
return Ainv
def direct_inverse(A):
"""Compute the inverse of matrices in an array of shape (N,N,M)"""
return np.linalg.inv(A.transpose(2,0,1)).transpose(1,2,0)
注意后面函数中的两次转置操作：形状为(N,N,M)的输入需被转置为形状为(M,N,N)的格式才能使用np.linalg.inv方法，然后将结果重新排列回(M,N,N)的形状。
在Python 3.6和numpy 1.14.0上使用IPython进行检查和计时：
In [5]: np.allclose(slow_inverse(A),direct_inverse(A))
Out[5]: True
In [6]: %timeit slow_inverse(A)
19 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [7]: %timeit direct_inverse(A)
1.3 ms ± 6.39 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
- Andras Deak -- Слава Україні
1
1谢谢。这是使用现代版本的numpy的正确答案。 - Nick Alger
4
Numpy-Blas调用并不总是最快的选择
在需要计算大量逆矩阵、特征值、小型3x3矩阵点积等问题时，我使用的numpy-MKL通常会被相当大幅度地超越。
这些外部Blas例程通常用于较大矩阵的问题，对于较小的矩阵，您可以编写标准算法或查看例如Intel IPP。
请注意，Numpy默认使用C顺序数组（最后一个维度变化最快）。 对于此示例，我从Matrix inversion (3,3) python - hard coded vs numpy.linalg.inv中取出了代码并进行了修改。
import numpy as np
import numba as nb
import time
@nb.njit(fastmath=True)
def inversion(m):
minv=np.empty(m.shape,dtype=m.dtype)
for i in range(m.shape[0]):
determinant_inv = 1./(m[i,0]*m[i,4]*m[i,8] + m[i,3]*m[i,7]*m[i,2] + m[i,6]*m[i,1]*m[i,5] - m[i,0]*m[i,5]*m[i,7] - m[i,2]*m[i,4]*m[i,6] - m[i,1]*m[i,3]*m[i,8])
minv[i,0]=(m[i,4]*m[i,8]-m[i,5]*m[i,7])*determinant_inv
minv[i,1]=(m[i,2]*m[i,7]-m[i,1]*m[i,8])*determinant_inv
minv[i,2]=(m[i,1]*m[i,5]-m[i,2]*m[i,4])*determinant_inv
minv[i,3]=(m[i,5]*m[i,6]-m[i,3]*m[i,8])*determinant_inv
minv[i,4]=(m[i,0]*m[i,8]-m[i,2]*m[i,6])*determinant_inv
minv[i,5]=(m[i,2]*m[i,3]-m[i,0]*m[i,5])*determinant_inv
minv[i,6]=(m[i,3]*m[i,7]-m[i,4]*m[i,6])*determinant_inv
minv[i,7]=(m[i,1]*m[i,6]-m[i,0]*m[i,7])*determinant_inv
minv[i,8]=(m[i,0]*m[i,4]-m[i,1]*m[i,3])*determinant_inv
return minv
#I was to lazy to modify the code from the link above more thoroughly
def inversion_3x3(m):
m_TMP=m.reshape(m.shape[0],9)
minv=inversion(m_TMP)
return minv.reshape(minv.shape[0],3,3)
#Testing
A = np.random.rand(1000000,3,3)
#Warmup to not measure compilation overhead on the first call
#You may also use @nb.njit(fastmath=True,cache=True) but this has also about 0.2s
#overhead on fist call
Ainv = inversion_3x3(A)
t1=time.time()
Ainv = inversion_3x3(A)
print(time.time()-t1)
t1=time.time()
Ainv2 = np.linalg.inv(A)
print(time.time()-t1)
print(np.allclose(Ainv2,Ainv))
性能
np.linalg.inv: 0.36 s
inversion_3x3: 0.031 s
- max9111
1
1对于未来的读者，我认为值得一提的是，在大多数机器上，使用标志
parallel=True和prange而不是range可以显著加快处理如此大的输入矩阵的速度。 - Jérôme Richard3
循环并不一定比其他方法慢，在这种情况下，使用循环也不会帮助你太多。但我建议如下：
import numpy as np
A = np.random.rand(100,3,3) #this is to makes it
#possible to index
#the matrices as A[i]
Ainv = np.array(map(np.linalg.inv, A))
对比你的解决方案和这个解决方案的时间会有轻微但明显的差异：
# The for loop:
100 loops, best of 3: 6.38 ms per loop
# The map:
100 loops, best of 3: 5.81 ms per loop
我尝试使用numpy的“vectorize”例程，希望创建一个更干净的解决方案，但我需要再仔细研究一下。在数组A中改变顺序可能是最重要的改变，因为它利用了numpy数组按列排序的事实，因此以这种方式进行线性读取数据略微更快。
- Dikkemik
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
for循环真的很慢吗？ - Robert Harveyinv函数求解3x3矩阵的逆约需要51.8微秒。执行for i in range(100): pass只需2.89微秒，因此每次求逆所需的循环开销可以完全忽略不计。切片操作所需的时间约为1.2微秒。我认为for循环速度在这里不是一个因素，只有timeit数据才能说服我相信其他可能性。 - DSMtimeit的解释，以及当你有嵌套循环时只应该关注优化最内层循环等等）。 - mgilson