4D matrix multiplication

I do the following in 4 loops and it takes ages to complete. Is there a way this code could be made more efficeint, without using parallel processing toolbox?

'steer' is a 136x101x101x16 matrix

'R' is a 136x16x16 matrix

'pow' and 'F' are 101x101 matrices.

pow = zeros(grdpts_y, grdpts_x); %grdpts_y, grdpts_x = 101
for l=1:nf  %nf = 136
    F = zeros(grdpts_y,grdpts_x); 
        for i=1:grdpts_x
            for j=1:grdpts_y
              F(i,j) = F(i,j) + 1./(squeeze(steer(l,i,j,:))'*squeeze(R(l,:,:))*squeeze(steer(l,i,j,:)));  
            end
        end
    F = F.*conj(F);
    pow = pow + F;    
end

Thanks in advance,

Kamran

Answer 1

[Matt J]

steer=reshape( permute(steer,[2,3,4,1]),101^2,[],136  );
R=permute(R,[2,3,1]);

F=1./sum( pagemtimes(conj(steer),R).*steer, 2);

F=reshape( abs(F).^2 ,101,101,[]);

pow=sum(F,3);

[Kamran]

Thanks for the answer. It is very fast, but I get different results from the loop implementation.

[Matt J]

Not me. See the comparison below:

[grdpts_y, grdpts_x] = deal(101);nf = 26;
steer=rand(26,101,101,16);
R=rand(26,16,16);

%Original version
pow = zeros(grdpts_y, grdpts_x); %grdpts_y, grdpts_x = 101
for l=1:nf  %nf = 136
    F = zeros(grdpts_y,grdpts_x); 
    for i=1:grdpts_x
        for j=1:grdpts_y
            F(i,j) = F(i,j) + 1./(squeeze(steer(l,i,j,:))'*squeeze(R(l,:,:))*squeeze(steer(l,i,j,:)));  
        end
    end
    F = F.*conj(F);
    pow = pow + F;    
end

pow0=pow;

%Proposed Version
steer=reshape( permute(steer,[2,3,4,1]),101^2,[],26  );
R=permute(R,[2,3,1]);

F=1./sum( pagemtimes(steer,R).*steer, 2);

F=reshape( abs(F).^2 ,101,101,[]);

pow=sum(F,3);

%Difference

difference=pow(:)-pow0(:);
[max(difference),min(difference)] %min and max differences

ans = 1×2

1.0e+-16 * 0.4163 -0.2776

[Kamran]

Strange. I run your code and got the exact same results as you, but when I run it on my data the results from the two implemntations are different.

It shouldn't matter but 'steer' and 'R' are complex valued! Probably, I am doing something wrong!

[Matt J]

The images don't tell us quantitatively what the differences are. What is the relative error:

norm(pow0(:)-pow(:))/norm(pow0(:))

[Kamran]

strangely enough

norm(pow0(:)-pow(:))/norm(pow0(:)) = 1

[Matt J]

Try

F=1./sum( pagemtimes(conj(steer),R).*steer, 2);

[Kamran]

That did it. Thank you very much.

[Kamran]

Can I ask one more question? I need to do the loop differntly for another method and I haven't quite understood how your solution works. I need to do the same for the following:

for l=1:nf
    F= zeros(grdpts_y, grdpts_x);
    for i=1:grdpts_x
        for j=1:grdpts_y
            [Vec, Val] = eig(squeeze(R(l,:,:)));
            [Val Seq] = sort(max(Val));
            Vec_s = Vec(:,Seq(nstat ,nstat));
            Vec_n = Vec(:,Seq(1:nstat-1));
            F(i, j) = F(i,j) + 1./(squeeze(steer(l,i,j,:))'*Vec_n*Vec_n'*squeeze(steer(l,i,j,:)));
        end
    end
    F = F.*conj(F);
    pow = pow + F;
end

[Matt J]

In your new version, F will always be real, non-negative, so I don't know why you would still be computing conj(F).

steer=reshape( permute(steer,[2,3,4,1]),101^2,[],136  );

Vec_n=cell(1,nf);
for l=1:nf
            [Vec, Val] = eig(squeeze(R(l,:,:)));
            [Val Seq] = sort(max(Val));
            Vec_s = Vec(:,Seq(nstat ,nstat));
            Vec_n{l}= Vec(:,Seq(1:nstat-1));
end
Vec_n=cat(3,Vec_n{:});


F=1./sum( abs(pagemtimes(conj(steer),Vec_n)).^2, 2);

F=reshape( abs(F).^2 ,101,101,[]);

pow=sum(F,3);

[Kamran]

Thank you very much. You are of course right. Thanks again for the prompt help.