MATLAB的bsxfun是最好的吗？Python's 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在实践中会更大。有没有更快的方法？奇怪的是，我需要先将singleton维度添加到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，您可能会尝试广播/隐式扩展，而不是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%。