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
采纳的回答
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);
10 个评论
Kamran
2021-10-18
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
Kamran
2021-10-18
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
2021-10-18
The images don't tell us quantitatively what the differences are. What is the relative error:
norm(pow0(:)-pow(:))/norm(pow0(:))
Kamran
2021-10-18
Matt J
2021-10-18
Try
F=1./sum( pagemtimes(conj(steer),R).*steer, 2);
Kamran
2021-10-19
Kamran
2021-10-19
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
2021-10-20
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Translated by 
