有没有办法用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]）的调用。这会在每次迭代中重新创建单位矩阵因为所有的数组都是相同的大小，这是浪费时间。通过预先分配单位矩阵来削减这个步骤可以减少约20％的时间（fast_inverse）。我的测试表明，预先分配数组或者分配数组它从一个结果列表并没有太大的区别

看更深的一层，你会发现拉包例程的调用，但它包裹在几个健全性检查，如果你去掉所有这些，只是调用lapack在你的for循环（因为你已经知道你的矩阵的尺寸，也许知道它是真实的，并不复杂），事情运行得更快（请注意，我已经使我的阵列更大）：

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

For循环确实不一定比替代方法慢得多，在这种情况下，它也不会帮你太多。但这里是一个建议：

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数组按列顺序排列的事实，因此数据的线性读出以这种方式稍微快一点。