15

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])
4

4 に答える 4

16

numpy.linalgコードで2レベル下がって焼かれていることがわかりました。numpy.linalg.invを見ると、これはnumpy.linalg.solve（A、inv（A.shape [0]）の呼び出しにすぎないことがわかります。これには、各反復で単位行列を再作成する効果があります。 forループ。すべての配列が同じサイズであるため、時間の無駄です。単位行列を事前に割り当てることでこの手順をスキップすると、時間が最大20％短縮されます（fast_inverse）。私のテストでは、配列の事前割り当てまたは割り当てが示唆されています。結果のリストからそれはあまり違いはありません。

1レベル深く見ると、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多次元配列をサポートするようになりました。これらは、行列のインデックスが最後にある行列のスタック（つまり、形状の配列(...,M,N,N)）として処理されます。これはnumpy1.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)

後者の関数の2つの転置に注意してください。動作するには、形状の入力を形状に(N,N,M)転置する必要があり、その後、結果を形状に並べ替える必要があります。(M,N,N)np.linalg.inv(M,N,N)

python3.6およびnumpy1.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ルーチン「vectorize」を使用しようとしましたが、それをもう一度調べる必要があります。配列Aの順序の変更は、おそらく最も重要な変更です。これは、numpy配列が列方向に順序付けられているという事実を利用しているため、データの線形読み取りがこの方法でわずかに高速になるためです。

于 2012-08-17T01:10:49.573 に答える
1

Numpy-Blasコールが常に最速の可能性であるとは限りません

多くの逆行列、固有値、小さな3x3行列の内積などを計算する必要がある問題では、私が使用するnumpy-MKLは、かなりの差でパフォーマンスが向上することがよくあります。

この外部Blasルーチンは通常、大きな行列の問題のために作成されます。小さな行列の場合は、標準のアルゴリズムを書き出すか、たとえばを見てください。インテルIPP。

NumpyはデフォルトでC順序の配列を使用することにも注意してください（最後の次元が最も速く変化します）。この例では、Matrix inversion（3,3）pythonからコードを取得しました-ハードコードされた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 に答える