15

通常我会在一个循环中反转一个 3x3 矩阵的数组，for如下例所示。不幸的是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])
4

4 回答 4

16

事实证明，您在 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”矩阵被覆盖并最终包含结果。这段代码现在给出了一致的结果。

于 2012-08-17T02:34:07.770 回答
8

自从提出并回答了这个问题以来，一些事情发生了变化，现在numpy.linalg.inv支持多维数组，将它们作为矩阵堆栈处理，矩阵索引在最后（换句话说，数组 shape (...,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)

请注意后一个函数中的两个转置：必须将 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)
于 2018-03-30T23:44:05.057 回答
3

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

于 2012-08-17T01:10:49.573 回答
1

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

在您必须计算大量逆数、特征值、小型 3x3 矩阵或类似情况的点积的问题上，我使用的 numpy-MKL 通常会以相当大的优势胜过。

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

请记住，Numpy 默认使用 C 排序数组（最后一个维度变化最快）。对于此示例，我从Matrix inversion (3,3) python - hard coded 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
于 2018-04-02T15:11:40.840 回答