numpy中的三张量产品

Question

我正在将一些代码从C转换为numpy,这个特定的代码段显示出一些显着的减速.似乎只是用numpy数组中的和来替换内部循环是不够的.有没有人建议恢复一些原始性能？

我宁愿避免使用scipy.weave.我想知道我是否能以某种方式使用卷积函数,但文档有点薄,并且似乎有几个选项(scipy.signal,ndimage,numpy.convolve等).是否可以使用爱因斯坦求和函数(einsum)或tensordot？对于tensordot我认为我需要分别做X,Y,Z？

我已经将这些值设置为随机,因此可以剪切和粘贴,但我应该注意MX,MY和MZ的某些结构M(i,j)是|i-j|和M(i,j)=M(j,i)

一般来说,回应一些建议的编辑 n[0] != n[1] != n[2]

n = np.array((50,50,50))
np.random.seed(1)
A = np.random.random(n)
C = np.zeros_like(A)
MX = np.random.random((n[0],n[0]))
MY = np.random.random((n[1],n[1]))
MZ = np.random.random((n[2],n[2]))

for i in xrange(n[0]):
    for j in xrange(n[1]):
        for k in xrange(n[2]):
            C[i,j,k]=(MX[i,:]*A[:,j,k]).sum() + (MY[j,:]*A[i,:,k]).sum() +(MZ[k,:]*A[i,j,:]).sum()

Answer 1

一种选择是使用np.einsum:

D = (  np.einsum('il,ljk->ijk', MX, A)
     + np.einsum('jl,ilk->ijk', MY, A)
     + np.einsum('kl,ijl->ijk', MZ, A))

结果是一样的:

In [14]: np.allclose(C, D)
Out[14]: True

这是一个基准测试,显示它快约170倍:

In [16]: %%timeit
   ....: for i in xrange(n[0]):
    for j in xrange(n[1]):
        for k in xrange(n[2]):
            C[i,j,k]=((MX[i,:]*A[:,j,k]).sum()
                      + (MY[j,:]*A[i,:,k]).sum()
                      + (MZ[k,:]*A[i,j,:]).sum())
1 loops, best of 3: 2.95 s per loop


In [21]: %timeit D = (np.einsum('il,ljk->ijk',MX,A) + np.einsum('jl,ilk->ijk',MY,A) + np.einsum('kl,ijl->ijk',MZ,A))
100 loops, best of 3: 17.4 ms per loop