How to aviod the following for loops

1 visualización (últimos 30 días)
Qiu Xu
Qiu Xu el 6 de Abr. de 2022
Comentada: Qiu Xu el 12 de Abr. de 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!

Respuesta aceptada

Voss
Voss el 6 de Abr. de 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)')

Más respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by