1

問題のセットアップ: サイズ n1、n2、n3=nx、ny、nz のデータの 3D (空間) グリッドがあり、a または b に対して nz=1 である可能性があります。そのグリッド内の各ポイントには、グリッドポイントごとにサイズ NDIM (通常は =4) のデータのベクトル (a) と、グリッドポイントごとにサイズ NDIMxNDIM の別の行列 (b) があります。ab や ba などの (ポイントごとの) ものを、メモリと CPU の両方で最も効率的に計算したいと考えています。

本質的に、私は一般化したいと思いますA loopless 3D matrix multiplication in python。効果があったように思います。しかし、私はそれを理解していません。グーグルとスタックオーバーフローを検索しても何の助けにもなりません。説明してさらに一般化してください！ありがとう！

import numpy as np
# gives a.b per point:
nx=5
ny=8
nz=3
a = np.arange(nx*ny*nz*4).reshape(4, nx,ny,nz)
b = np.arange(nx*ny*1*4*4).reshape(4, 4, nx,ny,1)
ctrue=a*0.0
for ii in np.arange(0,nx):
    for jj in np.arange(0,ny):
        for kk in np.arange(0,nz):
           ctrue[:,ii,jj,kk] = np.tensordot(a[:,ii,jj,kk],b[:,:,ii,jj,0],axes=[0,1])

c2 = (a[:,None,None,None] * b[:,:,None,None,None]).sum(axis=1).reshape(4,nx,ny,nz)
np.sum(ctrue-c2)
# gives 0 as required


# gives b.a per point:
ctrue2=a*0.0
for ii in np.arange(0,nx):
    for jj in np.arange(0,ny):
        for kk in np.arange(0,nz):
                ctrue2[:,ii,jj,kk] = np.tensordot(a[:,ii,jj,kk],b[:,:,ii,jj,0],axes=[0,0])

btrans=np.transpose(b,(1,0,2,3,4))
c22 = (a[:,None,None,None] * btrans[:,:,None,None,None]).sum(axis=1).reshape(4,nx,ny,nz)
np.sum(ctrue2-c22)
# gives 0 as required

# Note that only the single line for c2 and c22 are required -- the rest of the code is for testing/comparison to see if that line works.

# Issues/Questions:
# 1) Please explain why those things work and further generalize!
# 2) After reading about None=np.newaxis, I thought something like this would work:
    c22alt = (a[:,None,:,:,:] * btrans[:,:]).sum(axis=1).reshape(4,nx,ny,nz)
    np.sum(ctrue2-c22alt)
# but it doesn't.
# 3) I don't see how to avoid assignment of a separate btrans.  An np.transpose on b[:,:,None,None,None] doesn't work.

その他の関連リンク: Numpy: Multiplying a matrix with a 3d tensor -- Suggestion How to use numpy with 'None' value in Python?

4

1 に答える 1

1

そもそも、あなたのコードはひどく複雑すぎます。製品a.bとは次のb.aように簡略化できます。

c2 = (a * b).sum(axis=1)
c22 = (a * b.swapaxes(0, 1)).sum(axis=1)

代わりに;np.sum(ctrue - c2)を使用する必要があることに注意してください。np.all(ctrue == c2)2つの方法で同じ合計の結果が得られた場合、前者は間違った結果をもたらす可能性があります。

なぜこれが機能するのですか？単一の要素を考えてみましょう。

a0 = a[:, 0, 0, 0]
b0 = b[:, :, 0, 0, 0]

テンソルドットを取ることは。np.tensordot(a0, b0, axes=(0, 1))と同等(a0 * b0).sum(axis=1)です。これは放送のためです。の(4, )形状はの形状にa0ブロードキャストされ、配列は要素ごとに乗算されます。次に、軸を合計するとテンソルドットが得られます。(4, 4)b01

他の内積の場合、はと同じnp.tensordot(a0, b0, axes=(0, 0))です。はと同じです。はと同じです。転置することにより、のもう一方の軸に対して効果的にブロードキャストします。によって同じ結果を得ることができます。(a0 * b0.T).sum(axis=1)b0.Tb0.transpose()b0.swapaxes(0, 1)b0a0b0(a0[:, None] * b0).sum(axis=0)

NumPyの要素ごとの操作の優れている点は、高軸が対応している、またはブロードキャストできる限り、高軸を完全に無視できることです。したがって、何が機能しa0、b0（ほとんど）機能しaますb。

最後に、アインシュタインの縮約を使用して、これをはるかに明確にすることができます。

c2 = np.einsum('i...,ji...->j...', a, b)
c22 = np.einsum('i...,ij...->j...', a, b)
于 2012-08-14T14:57:53.243 に答える