有没有办法用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循环中调用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 中引入的。毫不奇怪，就性能而言，这是迄今为止最好的选择：

\n\n

import numpy as np\n\nA = np.random.rand(3,3,1000)\n\ndef slow_inverse(A):\n    """Looping solution for comparison"""\n    Ainv = np.zeros_like(A)\n\n    for i in range(A.shape[-1]):\n        Ainv[...,i] = np.linalg.inv(A[...,i])\n    return Ainv\n\ndef direct_inverse(A):\n    """Compute the inverse of matrices in an array of shape (N,N,M)"""\n    return np.linalg.inv(A.transpose(2,0,1)).transpose(1,2,0)\n

\n\n

请注意后一个函数中的两次转置： shape 的输入(N,N,M)必须转置为 shape(M,N,N)才能np.linalg.inv工作，然后结果必须重新排列回 shape (M,N,N)。

\n\n

使用 IPython 在 python 3.6 和 numpy 1.14.0 上进行检查和计时结果：

\n\n

In [5]: np.allclose(slow_inverse(A),direct_inverse(A))\nOut[5]: True\n\nIn [6]: %timeit slow_inverse(A)\n19 ms \xc2\xb1 138 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n\nIn [7]: %timeit direct_inverse(A)\n1.3 ms \xc2\xb1 6.39 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 1000 loops each)\n

\n