4

Python/Cython/Numpy を使用して多くの 4x4 行列を非常に迅速に乗算する方法を探しています。

私の現在の試みを示すために、計算する必要があるアルゴリズムがあります

A_1 * A_2 * A_3 * ... * A_N

どこでも

A_i != A_j

Python での例:

means = array([0.0, 0.0, 34.28, 0.0, 0.0, 3.4])
stds = array([ 4.839339, 4.839339, 4.092728, 0.141421, 0.141421, 0.141421])

def fn():
    steps = means+stds*numpy.random.normal(size=(60,6))
    A = identity(4)
    for step in steps:
        A = dot(A, transform_step_to_4by4(step))
%timeit fn()

1000 loops, best of 3: 570 us per loop

Cython/Numpy でこのアルゴリズムを実装すると、Eigen/C++ をすべて最適化して使用した同等のコードよりも約 100 倍遅くなります。しかし、私は本当にC++を使いたくありません。

4

3 に答える 3

6

乗算する各行列を生成するために Python 関数呼び出しを行う必要がある場合は、基本的にパフォーマンスが低下します。ただし、関数をベクトル化してtransform_step_to_4by4、形状(n, 4, 4)を含む配列を返すことができる場合は、次を使用して時間を節約できますmatrix_multiply。

import numpy as np
from numpy.core.umath_tests import matrix_multiply

matrices = np.random.rand(64, 4, 4) - 0.5

def mat_loop_reduce(m):
    ret = m[0]
    for x in m[1:]:
        ret = np.dot(ret, x)
    return ret

def mat_reduce(m):
    while len(m) % 2 == 0:
        m = matrix_multiply(m[::2], m[1::2])
    return mat_loop_reduce(m)

In [2]: %timeit mat_reduce(matrices)
1000 loops, best of 3: 287 us per loop

In [3]: %timeit mat_loop_reduce(matrices)
1000 loops, best of 3: 721 us per loop

In [4]: np.allclose(mat_loop_reduce(matrices), mat_reduce(matrices))
Out[4]: True

n の代わりに log(n) Python 呼び出しがあり、2.5 倍のスピードアップに適しています。n = 1024 の場合は 10 倍に近くなります。明らかmatrix_multiplyに ufunc であり、.reduceコードを許可するメソッドがあります。 Python でループを実行しないようにします。私はそれを実行することができませんでしたが、不可解なエラーが発生し続けます:

In [7]: matrix_multiply.reduce(matrices)
------------------------------------------------------------
Traceback (most recent call last):
  File "<ipython console>", line 1, in <module>
RuntimeError: Reduction not defined on ufunc with signature
于 2013-04-14T08:58:54.063 に答える
2

(60,6)配列を に変換する方法がわからないため、速度をメソッドと比較することはできません(4,4)が、これはシーケンスのドットを取得するために機能します。

A = np.arange(16).reshape(4,4)
B = np.arange(4,20).reshape(4,4)
C = np.arange(8,24).reshape(4,4)

arrs = [A, B, C]

P = reduce(np.dot, arrs)

これは以下と同等ですが、より高速である必要があります。

P = np.identity(4, dtype = arrs[0].dtype)
for arr in arrs:
    P = np.dot(P, arr)

タイミングテスト:

In [52]: def byloop(arrs):
   ....:     P = np.identity(4)
   ....:     for arr in arrs:
   ....:         P = np.dot(P, arr)
   ....:     return P
   ....: 

In [53]: def byreduce(arrs):
   ....:     return reduce(np.dot, arrs)
   ....: 

In [54]: byloop(arrs)
Out[54]: 
array([[  5104,   5460,   5816,   6172],
       [ 15728,  16820,  17912,  19004],
       [ 26352,  28180,  30008,  31836],
       [ 36976,  39540,  42104,  44668]])

In [55]: byreduce(arrs)
Out[55]: 
array([[ 5104,  5460,  5816,  6172],
       [15728, 16820, 17912, 19004],
       [26352, 28180, 30008, 31836],
       [36976, 39540, 42104, 44668]])

ここでlen(arrs) = 1000:

In [56]: timeit byloop(arrs)
1000 loops, best of 3: 1.26 ms per loop

In [57]: timeit byreduce(arrs)
1000 loops, best of 3: 656 us per loop
于 2013-04-13T20:46:14.647 に答える
1

このように小さい行列の場合、BLAS をまったく使用しない方が有利な場合があります。libSMMなどの小さな行列乗算ライブラリがあります(他にもあります)。

それらを f2py や cython でラップすることもできますし、 --- cython や fortran/f2py で独自に作成する方が簡単かもしれません。

于 2013-04-22T22:27:27.257 に答える