3

多くの 3x1 ベクトル ペアの外積をできるだけ速く計算しようとしています。これ

n = 10000
a = np.random.rand(n, 3)
b = np.random.rand(n, 3)
numpy.cross(a, b)

は正しい答えを与えますが、同様の質問に対するこの答えに動機付けられて、私はそれeinsumが私をどこかに連れて行くと思いました. 両方あることが分かった

eijk = np.zeros((3, 3, 3))
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1

np.einsum('ijk,aj,ak->ai', eijk, a, b)
np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b)

外積を計算しますが、そのパフォーマンスは期待外れです: どちらの方法も よりもはるかに悪いパフォーマンスを発揮しますnp.cross:

%timeit np.cross(a, b)
1000 loops, best of 3: 628 µs per loop

%timeit np.einsum('ijk,aj,ak->ai', eijk, a, b)
100 loops, best of 3: 9.02 ms per loop

%timeit np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b)
100 loops, best of 3: 10.6 ms per loop

sを改善する方法についてのアイデアはありますeinsumか？

4

2 に答える 2

4

の乗算回数einsum()は よりも多くcross()、最新の NumPy バージョンでは、cross()多くの一時配列を作成しません。したがってeinsum()、 より速くなることはできませんcross()。

クロスの古いコードは次のとおりです。

x = a[1]*b[2] - a[2]*b[1]
y = a[2]*b[0] - a[0]*b[2]
z = a[0]*b[1] - a[1]*b[0]

クロスの新しいコードは次のとおりです。

multiply(a1, b2, out=cp0)
tmp = array(a2 * b1)
cp0 -= tmp
multiply(a2, b0, out=cp1)
multiply(a0, b2, out=tmp)
cp1 -= tmp
multiply(a0, b1, out=cp2)
multiply(a1, b0, out=tmp)
cp2 -= tmp

高速化するには、cython または numba が必要です。

于 2016-09-23T14:48:00.207 に答える
4

np.tensordot次のように、最初のレベルで次元の 1 つを失い、次に を使用np.einsumして他の次元を失うために使用して行列乗算を行うことができます-

np.einsum('aik,ak->ai',np.tensordot(a,eijk,axes=([1],[1])),b)

aまたは、とbusingの間でブロードキャストされた要素ごとの乗算を実行し、 を使用np.einsumして一度に 2 つの次元を失うnp.tensordotこともできます。

np.tensordot(np.einsum('aj,ak->ajk', a, b),eijk,axes=([1,2],[1,2]))

のような新しい軸を導入することで、要素ごとの乗算を実行することもできましたが、a[...,None]*b[:,None]速度が低下するようです。

ただし、これらは提案されたnp.einsum唯一のベースのアプローチよりも優れた改善を示していますが、打ち負かすことはできませんnp.cross。

実行時テスト -

In [26]: # Setup input arrays
    ...: n = 10000
    ...: a = np.random.rand(n, 3)
    ...: b = np.random.rand(n, 3)
    ...: 

In [27]: # Time already posted approaches
    ...: %timeit np.cross(a, b)
    ...: %timeit np.einsum('ijk,aj,ak->ai', eijk, a, b)
    ...: %timeit np.einsum('iak,ak->ai', np.einsum('ijk,aj->iak', eijk, a), b)
    ...: 
1000 loops, best of 3: 298 µs per loop
100 loops, best of 3: 5.29 ms per loop
100 loops, best of 3: 9 ms per loop

In [28]: %timeit np.einsum('aik,ak->ai',np.tensordot(a,eijk,axes=([1],[1])),b)
1000 loops, best of 3: 838 µs per loop

In [30]: %timeit np.tensordot(np.einsum('aj,ak->ajk',a,b),eijk,axes=([1,2],[1,2]))
1000 loops, best of 3: 882 µs per loop
于 2016-09-23T14:23:50.370 に答える