python中使用numpy的多线性映射

Question

在 python 中，我有一个大小为(n,n,n) 的三维数组T 和一个大小为(n,k) 的二维数组。

在线性代数中，由T 定义并应用于W 的多线性映射在代码中为：

X3 = np.zeros((k,k,k))
    for i in xrange(k):
       for j in xrange(k):
            for t in xrange(k):
               for l in xrange(n):
                 for m in xrange(n)
                    for h in xrange(n):
                        X3[i, j, t] += M3[l, m, h] * W[l, i] * W[m, j] * W[h, t]

见

https://en.wikipedia.org/wiki/Multilinear_map

供参考。

这很慢。我想知道 numpy 中是否存在任何替代或任何预构建功能可以加快操作速度。

Answer 1

这是一种使用一系列点积的方法 -

# Get partial products and thus reach to final output
p1 = np.tensordot(W,M3,axes=(0,0))
p2 = np.tensordot(p1,W,axes=(1,0))
X3out = np.tensordot(p2,W,axes=(1,0))

Answer 2

爱因斯坦求和约定？使用np.einsum！

X3 = np.einsum('lmh,li,mj,ht->ijt', M3, W, W, W)

编辑：

出于好奇，我刚刚运行了一些基准测试，将 np.einsum 与 Divakar 的方法进行了比较。与np.einsum的区别很大！

import numpy as np

def approach1(a, b):
    x = np.einsum('lmh,li,mj,ht->ijt', a, b, b, b)

def approach2(a, b):
    p1 = np.tensordot(b, a, axes=(0,0))
    p2 = np.tensordot(p1, b, axes=(1,0))
    x = np.tensordot(p2, b, axes=(1,0))

n = 100
k = 10
a = np.random.random((n, n, n))
b = np.random.random((n, k))
%timeit approach1(a, b)      # => 1 loop, best of 3: 26 s per loop
%timeit approach2(a, b)      # => 100 loops, best of 3: 4.23 ms per loop

this question 对此进行了一些讨论。这一切似乎都归结为np.einsum 试图实现的普遍性——以能够将计算卸载到低级线性代数包为代价。