How to perform a sum over matrices in several dimensions?

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?

[Matt J]

Is there any sparsity in these arrays at all?

[Henrik]

Unfortunately no, as they stand here, all the elements of the matrices are non-zero numbers

[Patrik Ek]

Are the matrices calculated in the same loop? I mean, the most time consuming step is most likely to calculate the matrix M5, is it possible to perform the matrix calculations in another way? I am not really sure of how you want to calculate the matrix. Are there any way to clarify that?

[Walter Roberson]

Do some factoring. Your M1 and M2 indices are independent of qx and qy, so you can calculate the M1 * M2 part independently.

[Henrik]

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.

Answer 1

[Matt J]

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

[Henrik]

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!