How to eliminate for loop and perform the computation parallely?

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

Answer 1

[chicken vector]

%% 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.

[Kalasagarreddi Kottakota]

Hi @chicken vector, If I modify the bolded, [0 1; 0 1] to [0 0;0 1], How can I modify the code you provided?

result(:,i) = (A'*A + lambda(1,i) * [0 1 ; 0 1]) \ (A' * B(:,i));

[chicken vector]

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]

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]

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?