MATLAB 的 bsxfun 是不是最好的？ Python的numpy.einsum？

Question

我需要尽可能高效地实现一个非常大的乘法和求和运算。到目前为止，我发现的最佳方法是 MATLAB 中的bsxfun，我将问题表述为：

L = 10000;
x = rand(4,1,L+1);
A_k = rand(4,4,L);
tic
for k = 2:L
    i = 2:k;
    x(:,1,k+1) = x(:,1,k+1)+sum(sum(bsxfun(@times,A_k(:,:,2:k),x(:,1,k+1-i)),2),3);
end
toc

请注意，L 在实践中会更大。有没有更快的方法？奇怪的是，我需要先将单例维度添加到x，然后在其上添加sum，但否则我无法让它工作。

它仍然比我尝试过的任何其他方法快得多，但对于我们的应用程序来说还不够。听说 Python 函数 numpy.einsum 可能效率更高，但在考虑移植我的代码之前，我想先在这里问一下。

我正在使用 MATLAB R2017b。

Answer 1

我相信你的两个总结都可以删除，但我暂时只删除了更容易的一个。第二维的求和是微不足道的，因为它只影响A_k 数组：

B_k = sum(A_k,2);
for k = 2:L
    i = 2:k;
    x(:,1,k+1) = x(:,1,k+1) + sum(bsxfun(@times,B_k(:,1,2:k),x(:,1,k+1-i)),3);
end

通过这一单一更改，我的笔记本电脑上的运行时间从约 8 秒减少到约 2.5 秒。

也可以通过将 times+sum 转换为矩阵向量积来删除第二个求和。它需要一些单一的摆弄才能获得正确的维度，但是如果您定义一个辅助数组 B_k 并反转第二个维度，您可以使用此辅助数组 C_k 生成剩余总和为 ~x*C_k，给出或给reshape打几个电话。

---

所以仔细观察后，我意识到我最初的评估过于乐观：你在剩余任期内的两个维度都有乘法，所以它不是一个简单的矩阵乘积。无论如何，我们可以将该项重写为矩阵乘积的对角线。这意味着我们正在计算一堆不必要的矩阵元素，但这似乎仍然比bsxfun 方法稍快，我们也可以摆脱你讨厌的单例维度：

L = 10000;
x = rand(4,L+1);
A_k = rand(4,4,L);
B_k = squeeze(sum(A_k,2)).';

tic
for k = 2:L
    ii = 1:k-1;
    x(:,k+1) = x(:,k+1) + diag(x(:,ii)*B_k(k+1-ii,:));
end
toc

这在我的笔记本电脑上运行大约 2.2 秒，比之前获得的大约 2.5 秒要快一些。

Answer 2

由于您使用的是新版本的 Matlab，您可以尝试广播 /implicit expansion 而不是 bsxfun：

x(:,1,k+1) = x(:,1,k+1)+sum(sum(A_k(:,:,2:k).*x(:,1,k-1:-1:1),3),2);

我还更改了求和顺序并删除了i 变量以进一步改进。在我的机器上，使用 Matlab R2017b，L = 10000 的速度提高了大约 25%。