4

私は1年以上前にここで質問者とほぼ同じ状況にいます: kxnxn行列を反転またはドットする高速な方法

したがって、次元 (N、M、M) のインデックス a[n、i、j] を持つテンソルがあり、N の各 n に対して M*M 正方行列部分を反転したいと考えています。

たとえば、私が持っているとします

In [1]:    a = np.arange(12)
           a.shape = (3,2,2)
           a

Out[1]: array([[[ 0,  1],
                  [ 2,  3]],

                  [[ 4,  5],
                  [ 6,  7]],

                  [[ 8,  9],
                  [10, 11]]])

次に、for ループの反転は次のようになります。

In [2]: inv_a = np.zeros([3,2,2])
        for m in xrange(0,3):
            inv_a[m] = np.linalg.inv(a[m])
        inv_a

Out[2]: array([[[-1.5,  0.5],
                  [ 1. ,  0. ]],

                  [[-3.5,  2.5],
                  [ 3. , -2. ]],

                  [[-5.5,  4.5],
                 [ 5. , -4. ]]])

githubのこの問題によると、これは NumPy 2.0 で実装されるようです...

seberg が github issue thread で指摘したように、開発バージョンをインストールする必要があると思いますが、今すぐベクトル化された方法でこれを行う別の方法はありますか?

4

1 に答える 1

4

更新: NumPy 1.8 以降では、関数numpy.linalgは一般化されたユニバーサル関数です。つまり、次のようなことができるようになりました。

import numpy as np
a = np.random.rand(12, 3, 3)
np.linalg.inv(a)

これにより、各 3x3 配列が反転され、結果が 12x3x3 配列として返されます。numpy 1.8 リリース ノートを参照してください。

元の回答:

は比較的小さいのでN、一度にすべての行列に対して手動で LU 分解を計算するのはどうでしょうか。これにより、関連する for ループが比較的短くなります。

通常の NumPy 構文でこれを行う方法は次のとおりです。

import numpy as np
from numpy.random import rand

def pylu3d(A):
    N = A.shape[1]
    for j in xrange(N-1):
        for i in xrange(j+1,N):
            #change to L
            A[:,i,j] /= A[:,j,j]
            #change to U
            A[:,i,j+1:] -= A[:,i,j:j+1] * A[:,j,j+1:]

def pylusolve(A, B):
    N = A.shape[1]
    for j in xrange(N-1):
        for i in xrange(j+1,N):
            B[:,i] -= A[:,i,j] * B[:,j]
    for j in xrange(N-1,-1,-1):
        B[:,j] /= A[:,j,j]
        for i in xrange(j):
            B[:,i] -= A[:,i,j] * B[:,j]

#usage
A = rand(1000000,3,3)
b = rand(3)
b = np.tile(b,(1000000,1))
pylu3d(A)
# A has been replaced with the LU decompositions
pylusolve(A, b)
# b has been replaced to the solutions of
# A[i] x = b[i] for each A[i] and b[i]

私が書いたようにpylu3d、LU 分解を計算するために A をその場で変更します。N各xN行列をその LU 分解で置き換えた後、行列システムの右辺を表すx配列pylusolveを解くために使用できます。その場で変更し、システムを解決するために適切な後方置換を行います。書かれているように、この実装にはピボットが含まれていないため、数値的に安定していませんが、ほとんどの場合、十分に機能するはずです。MNbb

配列がメモリ内でどのように配置されているかにもよりますが、Cython を使用する方がおそらく少し高速です。同じことを行う 2 つの Cython 関数を次に示しますが、M最初に反復処理を行います。ベクトル化されていませんが、比較的高速です。

from numpy cimport ndarray as ar
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def lu3d(ar[double,ndim=3] A):
    cdef int n, i, j, k, N=A.shape[0], h=A.shape[1], w=A.shape[2]
    for n in xrange(N):
        for j in xrange(h-1):
            for i in xrange(j+1,h):
                #change to L
                A[n,i,j] /= A[n,j,j]
                #change to U
                for k in xrange(j+1,w):
                    A[n,i,k] -= A[n,i,j] * A[n,j,k]

@cython.boundscheck(False)
@cython.wraparound(False)
def lusolve(ar[double,ndim=3] A, ar[double,ndim=2] b):
    cdef int n, i, j, N=A.shape[0], h=A.shape[1]
    for n in xrange(N):
        for j in xrange(h-1):
            for i in xrange(j+1,h):
                b[n,i] -= A[n,i,j] * b[n,j]
        for j in xrange(h-1,-1,-1):
            b[n,j] /= A[n,j,j]
            for i in xrange(j):
                b[n,i] -= A[n,i,j] * b[n,j]

Numba を使用することもできますが、この場合は Cython ほど高速に実行できませんでした。

于 2013-07-29T16:01:44.757 に答える