# How to eliminate for loop and perform the computation parallely?

3 vues (au cours des 30 derniers jours)
Kalasagarreddi Kottakota le 19 Avr 2023
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
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

chicken vector le 19 Avr 2023
Modifié(e) : chicken vector le 19 Avr 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 commentairesAfficher 2 commentaires plus anciensMasquer 2 commentaires plus anciens
Kalasagarreddi Kottakota le 21 Avr 2023
Modifié(e) : Kalasagarreddi Kottakota le 21 Avr 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 le 21 Avr 2023
Modifié(e) : chicken vector le 21 Avr 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?

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Creating and Concatenating Matrices dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by