How to aviod the following for loops

1 view (last 30 days)
Qiu Xu
Qiu Xu on 6 Apr 2022
Commented: Qiu Xu on 12 Apr 2022
Hi, I want to perform the following calculations on Matlab
N=1000;
M=200;
K=256;
U=rand(N,M)+i*rand(N,M); % where i is complex
V=rand(N,M)+i*rand(N,M); % where i is complex
B=rand(N,K);
C=zeros(1,K);
D=zeros(1,K);
const_c=2;
const_d=0.5;
for m=1:M
ux=zeros(1,K);
vx=zeros(1,K);
for n=1:N
ux=ux+U(n,m)*B(n,:);
vx=vx+V(n,m)*B(n,:);
end
C=C+const_c*abs(vx).^2;
D=D+const_d*ux.*conj(vx);
end
There are two loops, for large M and N, the loops is very slowly. So How can I aviod the loops without using parallel computing? I think it may be realized by using 'reshape' function, but I do not know how to use it for this code.
Thanks very much!

Accepted Answer

Voss
Voss on 6 Apr 2022
N=1000;
M=200;
K=256;
U=rand(N,M)+1i*rand(N,M);
V=rand(N,M)+1i*rand(N,M);
B=rand(N,K);
const_c=2;
const_d=0.5;
% no loops
tic
ux = U.'*B;
vx = V.'*B;
C_new = const_c*sum(abs(vx).^2,1);
D_new = const_d*sum(ux.*conj(vx),1);
toc
Elapsed time is 0.013907 seconds.
% with loops
tic
C=zeros(1,K);
D=zeros(1,K);
for m=1:M
ux=zeros(1,K);
vx=zeros(1,K);
for n=1:N
ux=ux+U(n,m)*B(n,:);
vx=vx+V(n,m)*B(n,:);
end
C=C+const_c*abs(vx).^2;
D=D+const_d*ux.*conj(vx);
end
toc
Elapsed time is 1.765365 seconds.
% results are not exactly the same, but the differences are ~1e-8 out of ~1e7:
[max(abs(C-C_new)) max(abs(D-D_new))]
ans = 1×2
1.0e-07 * 0.5960 0.1677
% which is about as accurate as you can get for floating-point numbers around 1e7:
[max(eps(C)) max(eps(D))]
ans = 1×2
1.0e-08 * 0.7451 0.1863
% plots to inspect the differences:
subplot(3,1,1)
plot(C,'--b','LineWidth',2)
hold on
plot(C_new,'-c')
title('C,C\_new')
subplot(3,1,2)
plot(real(D),'--r','LineWidth',2)
hold on
plot(real(D_new),'-m')
title('real(D,D\_new)')
subplot(3,1,3)
plot(imag(D),'--r','LineWidth',2)
hold on
plot(imag(D_new),'-m')
title('imag(D,D\_new)')

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!