How to perform a sum over matrices in several dimensions?

1 次查看（过去 30 天）
Henrik
Henrik 2014-1-13
评论： Henrik 2014-1-13
Hello I have a rather large sum that I want to evaluate in MATLAB. Symbolically it looks like this:
huge_sum(qx,qy)=sum_n(sum_m(sum_kx(sum_ky(sum_i(sum_j( M1(n,m,kx,ky)*M2(n,kx,ky,i,j)*M3(m,kx+qx,ky+qy,i,j)*M4(qx,qy,i,j)
) ) ) ) ) ),
where each of the M's is a matrix. The indices run from 1 to
qx,qy,kx,ky: 25
i,j,m,n: 16
My current implementation is rather slow; it's simply a bunch of nested for loops, defining a matrix,M5(m,n,kx,ky,i,j), which I then sum for each value of qx.
The problem is that this code will take around 100 days to run, which is of course too much! Can anyone suggest a smarter way to do this?
  5 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2014-1-13
Do some factoring. Your M1 and M2 indices are independent of qx and qy, so you can calculate the M1 * M2 part independently.
Henrik
Henrik 2014-1-13
Calculating the matrices M1-M4 is rather quick and is done beforehand. I know it's the calculation of M5 which is slow; at the moment I calculate its values one entry at a time. I would like to vectorize it somehow, but I can't see how to do that when the matrices have different sizes..
@Walter thanks, but I actually realized I made a mistake in what I wrote here - M1 depends on qx and qy as well.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2014-1-13
编辑：Matt J 2014-1-13
Warning. Not tested.
op=@(A,B,dim) squeeze(sum(bsxfun(@times,A,B),dim))
M2=reshape(M2,16,1,25,25,256); %M2(n,1,kx,ky,z)
T=op(M1,M2,1); %T(m,kx,ky,z)
M3=M3(:,:,:,:); %M3(m,kx+qx,ky+qy,z)
M4=M4(:,:,:); %M4(qx,qy,z)
sx=size(M3,2);
sy=size(M3,3);
Tr=reshape(T,16,25,25,1,1,256);
M3r=reshape(M3,16,1,1,sx,sy,256);
U=op(Tr,M3r,1); %U(kx,ky,kx+qx,ky+qy,z)
U=permute(U,[1,3,2,4,5]); %U(kx,kx+qx,ky,ky+qy,z)
U=reshape(U,25*sx,25*sy, 256);
K=1:25;
for qx=K
idx=sub2ind([25,sx],K,K+qx);
for qy=K
idy=sub2ind([25,sy],K,K+qy);
u=reshape(U(idx,idy,:)[],256);
v=reshape(M3(qx,qy,:),[],256);
M5(qx,qy)=sum(u,1)*v.';
end
end
  1 个评论
Henrik
Henrik 2014-1-13
I can't say I understand exactly how this code works, but it seems to be doing the right thing. I will try and think of a way to test it. Thanks for the effort!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by