How to eliminate for loop and perform the computation parallely?

Question

Kalasagarreddi Kottakota 2023-4-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1949798-how-to-eliminate-for-loop-and-perform-the-computation-parallely

编辑： chicken vector 2023-4-21

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

chicken vector 2023-4-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1949798-how-to-eliminate-for-loop-and-perform-the-computation-parallely#answer_1218963

编辑：chicken vector 2023-4-19

在 MATLAB Online 中打开

%% 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 个评论
显示 2更早的评论隐藏 2更早的评论

chicken vector 2023-4-21

编辑：chicken vector 2023-4-21

在 MATLAB Online 中打开

Consider the line:

num = (repmat(A'*A, 1, 1, Nt) + repmat(reshape([zeros(size(lambda)); lambda], 1, 2, []), 2, 1));

[zeros(size(lambda)); lambda]

This defines an array having all zeros in the first row and the values of lambda in the second.

reshape([zeros(size(lambda)); lambda], 1, 2, [])

Then, we reshape this matrix into a 3d array such that we have a collection of N 1x2 arrays, where N is the number of elements inside lambda.

the j array among this collection will be in te form [0 lambda(j)]

repmat(reshape([zeros(size(lambda)); lambda], 1, 2, []), 2, 1)

Finally we repeat these array, doubling the rows and our collection of arrays becomes:

[0 lambda(j) ; 0 lambda(j)]

Doing so we replicated the behaviour of the matrix [0 1 ; 0 1] that you previously define.

As you can see, the new matrix [0 0 ; 0 1] does not share any simmetry as the previous one, so one more passage might be required.

% Define a collection N of 2x2 arrays, filled with zeros:
LAM = zeros(2, 2, length(lambda));
% Place values of lambda in [0 0 ; 0 1], which is position (2,2,:):
LAM(2,2,:) = lambda
% Update variable num:
num = repmat(A'*A, 1, 1, Nt) + LAM;

You can use this approach also with the previous definition by doing:

% Place values of lambda in [0 1 ; 0 1], which is position (:,2,:):
LAM(:,2,:) = repmat(lambda,2,1)

Kalasagarreddi Kottakota 2023-4-21

编辑：Kalasagarreddi Kottakota 2023-4-21

在 MATLAB Online 中打开

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))))