MATLAB的bsxfun最好吗？Python的numpy.einsum？

Question

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

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

由于您使用的是新版本的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,这个速度提高了大约25%L = 10000.

Answer 2

我相信你的两个总结都可以删除,但我暂时只删除了更简单的总结.第二维的总和是微不足道的,因为它只影响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秒.

通过将时间+和转换为矩阵向量乘积,也可以去除第二个求和.它需要一些单独摆弄以获得正确的尺寸,但是如果你定义一个B_k第二个维度反转的辅助数组,你可以x*C_k用这个辅助数组生成剩余的和〜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秒快一些.