有没有一种方法可以用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

自从这个问题被提出和回答以来，有些事情已经发生了变化，现在^{}支持多维数组，将它们作为矩阵索引最后的矩阵堆栈来处理（换句话说，形状数组(...,M,N,N)）。好像是introduced in 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)

注意后一个函数中的两个转置：shape (N,N,M)的输入必须转置到shape (M,N,N)才能使np.linalg.inv工作，然后结果必须重新排列到shape (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 2

结果你在numpy.linalg代码中被烧掉了两层。如果你看看numpy.linalg.inv，你会发现它只是一个对numpy.linalg.solve的调用（a，inv（a.shape[0]）。这会在for循环的每次迭代中重新创建标识矩阵。因为所有的数组都是相同的大小，这是浪费时间。通过预先分配身份矩阵跳过此步骤可节省约20%的时间（快速逆）。我的测试表明，预先分配数组或从结果列表中分配数组没有多大区别。

再深入一层，您会发现对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

编辑：我没有仔细查看solve返回的内容。结果是“b”矩阵被覆盖，最后包含结果。这段代码现在给出了一致的结果。

Answer 3

因为循环并不一定比其他选择慢得多，而且在这种情况下，它也不会对您有多大帮助。但这里有一个建议：

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