python - 有没有办法用 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

编辑: 我没有仔细查看 solve 返回的内容。事实证明，'b' 矩阵被覆盖并包含最后的结果。此代码现在提供一致的结果。