How to aviod the following for loops

1 回表示 (過去 30 日間)
Qiu Xu
Qiu Xu 2022 年 4 月 6 日
コメント済み: Qiu Xu 2022 年 4 月 12 日
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!

採用された回答

Voss
Voss 2022 年 4 月 6 日
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)')
  1 件のコメント
Qiu Xu
Qiu Xu 2022 年 4 月 12 日
Thanks very much for your helps

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLogical についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by