Filter löschen
Filter löschen

How to eliminate for loop and perform the computation parallely?

4 Ansichten (letzte 30 Tage)
Hi, I am looking a way to reduce the time for computation of the problem of this size. This code is a larger part of a big code where this code is ran several 100 times.
The data will be as large as shown in the code. I am looking a way to eliminate the loop and perform complete problem in a sigle go, to make this step fast. Could some one propose a solution without disturbing the fomulation through matrix multiplication istead of loop?
clear all;
close all;
fs = 35000;
T = 100; % s
Nt = fs*T;
M = 60;
a1 = randn(M,Nt);
a2 = randn(M,Nt);
lambda = randn(1,Nt);
B = randn(M,Nt);
result = zeros(2,Nt);
tic
est = zeros(M,Nt);
for i= 1:Nt
A = [a1(:,1) a2(:,1)];
result(:,i) = (A'*A + lambda(1,i)*[0 1;0 1])\(A'*B(:,i));
est(:,i) = a1(:,i)*result(1,i)+a2(:,i)*result(2,i);
end
toc

Akzeptierte Antwort

chicken vector
chicken vector am 19 Apr. 2023
Bearbeitet: chicken vector am 19 Apr. 2023
%% Data:
fs = 35000;
T = 100;
Nt = fs*T;
M = 60;
a1 = randn(M,Nt);
a2 = randn(M,Nt);
lambda = randn(1,Nt);
B = randn(M,Nt);
result = zeros(2,Nt);
est = zeros(M,Nt);
A = [a1(:,1) a2(:,1)];
%% Vectorised method:
tic
num = (repmat(A'*A, 1, 1, Nt) + repmat(reshape([zeros(size(lambda)); lambda], 1, 2, []), 2, 1)); % \ (A' * B)
dem = reshape((A' * B), 2, 1, []);
resultVectorised = reshape(pagemldivide(num,dem), 2, []);
estVectorised = a1 .* result2(1,:) + a2 .* result2(2,:);
timeVectorised = toc;
%% Loop method:
tic
for i = 1 : Nt
result(:,i) = (A'*A + lambda(1,i) * [0 1 ; 0 1]) \ (A' * B(:,i));
est(:,i) = a1(:,i) * result(1,i) + a2(:,i) * result(2, i);
end
timeLoop = toc;
%% Output:
fprintf('Vectorisation: %.3f s\n', timeVectorised)
fprintf('Loop: %.3f s\n', timeLoop)
fprintf('Improvement: %.3f\n', timeLoop / timeVectorised)
fprintf('Absolute error: %e\n', max(max(abs(est - estVectorised))))
Vectorisation: 1.080 s
Loop: 20.117 s
Improvement: 18.632
Absolute error: 1.776357e-15
I think the vectorised method can be improved since, personally, I found hard to implement it without precisely know what I am doing.
Moreover, I think the small numerical error occurs due to pagemldivide which might use a different method from \, but I am not sure.
  4 Kommentare
Kalasagarreddi Kottakota
Kalasagarreddi Kottakota am 21 Apr. 2023
Bearbeitet: Kalasagarreddi Kottakota am 21 Apr. 2023
Hi @chicken vector, thanks and I have one final question. In the code if my AA is defined like the below, how can I write the dem? I modified the loop appropriatly.
clear all:
%% Data:
fs = 2;
T = 3;
Nt = fs*T;
M = 5;
a1 = randn(M,Nt);
a2 = randn(M,Nt);
lambda = randn(1,Nt);
B = randn(M,Nt);
result = zeros(2,Nt);
est = zeros(M,Nt);
%% Vectorised method:
tic
LAM = zeros(2, 2, length(lambda));
LAM(2,2,:) = lambda;
%%% modification %%%%%
AA = zeros(M, 2, length(lambda));
for i=1:Nt
AA(:,:,i) = [a1(:,i) a2(:,i)];
end
num = pagemtimes(pagectranspose(AA),AA) + LAM; % \ (AA' * B)
dem = reshape(pagemtimes(pagectranspose(AA),B), 2, 1, []);%%%%-------------?
resultVectorised = reshape(pagemldivide(num,dem), 2, []);
estVectorised = a1 .* resultVectorised(1,:) + a2 .* resultVectorised(2,:);
timeVectorised = toc;
%% Loop method:
tic
for i = 1 : Nt
A = [a1(:,i) a2(:,i)]; %%% modification
result(:,i) = (A'*A + lambda(1,i) * [0 0 ; 0 1]) \ (A' * B(:,i));
est(:,i) = a1(:,i) * result(1,i) + a2(:,i) * result(2, i);
end
timeLoop = toc;
%% Output:
fprintf('Vectorisation: %.3f s\n', timeVectorised)
fprintf('Loop: %.3f s\n', timeLoop)
fprintf('Improvement: %.3f\n', timeLoop / timeVectorised)
fprintf('Absolute error: %e\n', max(max(abs(est - estVectorised))))
chicken vector
chicken vector am 21 Apr. 2023
Bearbeitet: chicken vector am 21 Apr. 2023
Sorry but I am not understand what you're asking.
The loop is the same as before, is the vectorised method that is changed.
Can you articulate?

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Creating and Concatenating Matrices finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by