有没有办法用 numpy 有效地反转矩阵数组

Question

通常我会在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])

Answer 1

事实证明，你在numpy.linalg代码中被烧毁了两级。 如果你看一下numpy.linalg.inv，你会发现它只是对numpy.linalg.solve（A， inv（A.shape[0]）的调用。 这具有在 for 循环的每次迭代中重新创建单位矩阵的效果。 由于所有数组的大小都相同，因此这是浪费时间。 通过预先分配单位矩阵跳过此步骤可将时间缩短 ~20% 的时间 （fast_inverse）。 我的测试表明，预先分配数组或从结果列表中分配数组并没有太大区别。

再深入一层，你会发现对 lapack 例程的调用，但它包含在几个健全性检查中。 如果你把所有这些去掉，只在你的for循环中调用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

编辑：我没有仔细研究解决中返回的内容。 事实证明，"b"矩阵被覆盖并最终包含结果。 此代码现在提供一致的结果。

Answer 2

自从这个问题被提出和回答以来，一些事情发生了变化，现在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)

Answer 3

Numpy-Blas调用并不总是最快的可能

在必须计算大量逆、特征值、小 3x3 矩阵的点积或类似情况的问题上，我使用的 numpy-MKL 通常可以远远超过

。

这个外部 Blas 例程通常是为较大矩阵的问题而设计的，对于较小的矩阵，您可以编写标准算法或查看例如。英特尔 IPP。

还请记住，Numpy 默认使用 C 排序数组（最后一个维度变化最快）。对于这个例子，我从矩阵反转（3,3）python中获取了代码 - 硬编码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

Answer 4

对于循环确实不一定比替代方案慢很多，而且在这种情况下，它对您没有多大帮助。但这里有一个建议：

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 例程"矢量化"，希望创建一个更干净的解决方案，但我必须重新研究一下。数组 A 中排序的变化可能是最重要的变化，因为它利用了 numpy 数组按列排序的事实，因此数据的线性读出以这种方式略快。