matlab：矢量化 4D 矩阵和

Question

我需要在 MATLAB 中执行以下计算：

其中 w 和 v 是具有 N 个元素的向量，A 是一个四维矩阵（N^4 个元素）。这可以通过以下迂腐的代码来实现：

N=10;
A=rand(N,N,N,N);
v=rand(N,1);
w=zeros(N,1);

for pp=1:N
  for ll=1:N
    for mm=1:N
      for nn=1:N
        w(pp)=w(pp)+A(pp,ll,mm,nn)*v(ll)*v(mm)*conj(v(nn));
      end
    end
  end
end

这是非常缓慢的。有没有办法在 MATLAB 中对这种总和进行矢量化？

Answer 1

方法#1

很少有reshape's 和matrix multiplication -

A1 = reshape(A,N^3,N)*conj(v)
A2 = reshape(A1,N^2,N)*v
w = reshape(A2,N,N)*v

---

方法 #2

有一个 bsxfun 、 reshape 和 matrix-multiplication -

A1 = reshape(A,N^3,N)*conj(v)
vm = bsxfun(@times,v,v.')
w = reshape(A1,N,N^2)*vm(:)

---

基准测试

本部分比较了本文中列出的两种方法的运行时，Shai's post 中的第一个测试方法和问题中列出的原始方法。

基准代码

N=100;
A=rand(N,N,N,N);
v=rand(N,1);

disp('----------------------------------- With Original Approach')
tic
%// .... Code from the original post   ...//
toc

disp('----------------------------------- With Shai Approach #1')
tic
s4 = sum( bsxfun( @times, A, permute( conj(v), [4 3 2 1] ) ), 4 ); 
s3 = sum( bsxfun( @times, s4, permute( v, [3 2 1] ) ), 3 );
w2 = s3*v; 
toc

disp('----------------------------------- With Divakar Approach #1')
tic
A1 = reshape(A,N^3,N)*conj(v);
A2 = reshape(A1,N^2,N)*v;
w3 = reshape(A2,N,N)*v;
toc

disp('----------------------------------- With Divakar Approach #2')
tic
A1 = reshape(A,N^3,N)*conj(v);
vm = bsxfun(@times,v,v.');
w4 = reshape(A1,N,N^2)*vm(:);
toc

运行时结果

----------------------------------- With Original Approach
Elapsed time is 4.604767 seconds.
----------------------------------- With Shai Approach #1
Elapsed time is 0.334667 seconds.
----------------------------------- With Divakar Approach #1
Elapsed time is 0.071905 seconds.
----------------------------------- With Divakar Approach #2
Elapsed time is 0.058877 seconds.

结论

这篇文章中的第二种方法似乎比原来的方法提供了大约 80x 加速。

Answer 2

您可以尝试使用bsxfun。

假设 v 是一个 N-by-1 列向量（否则，应该稍微修改排列）。

% sum over n (4th dim)
s4 = sum( bsxfun( @times, A, permute( conj(v), [4 3 2 1] ) ), 4 ); 

现在中间结果只有N-by-N-by-N。

% sum over m (3rd dim)
s3 = sum( bsxfun( @times, s4, permute( v, [3 2 1] ) ), 3 )

继续最后的总和

% sum over l (2nd dim)
w = s3*v; 

---

想一想，您是否考虑过在其 multidim 版本中使用dot？我没有测试它，但它应该可以工作（可能是一些小的更正）。

s4 = dot( A, permute( conj(v), [4 3 2 1] ), 4 );
s3 = dot( s4, permute( v, [3 2 1] ), 3 );
w = s3*v;