Effacer les filtres
Effacer les filtres

Speeding up matrix operations

4 vues (au cours des 30 derniers jours)
bil
bil le 23 Mar 2023
Commenté : bil le 24 Mar 2023
Hey all,
I would like advice on speeding up the following operation:
Suppose we have a filled NxN matrix A, a filled NxM matrix B with no constraints on the values of each element in matrices A and B. Now we also have an empty NxM matrix C whose elements are defined as follows:
for i = 1:N
for j = 1:M
C(i,j) = trapz(A(i,:).*B(:,j).')
end
end
Clearly for very large N and M, this becomes a very slow process and it seems possible to vectorize or do away with the loop but I am unsure how.
Thanks.
  2 commentaires
Cameron
Cameron le 23 Mar 2023
Did you mean
C(i,j) = trapz(A(i,:),B(:,j))
or
C(i,j) = sum(trapz(A(i,:).*B(:,j)))
or something else? Because when I run a sample bit of code like this
x = ([1:10])';
A = x*(1:10); %10 x 10 array
B = x./(1:10); %10 x 10 array
N = 1:size(A,1);
M = 1:size(B,1);
for i = N
for j = M
C(i,j) = trapz(A(i,:).*B(:,j))
end
end
the value for C will be a 1x10 array which cannot fit into your original C(i,j) index.
bil
bil le 23 Mar 2023
Modifié(e) : bil le 23 Mar 2023
Ah, good catch sorry about that. I wanted trapz of the elementwise product between the i'th row of A and the j'th column of B, so it should really be
C(i,j) = trapz(A(i,:).*B(:,j).')
I have fixed it in my original post.

Connectez-vous pour commenter.

Réponse acceptée

Bruno Luong
Bruno Luong le 23 Mar 2023
Instead of calling trapz, use matrix multiplication, and this probably beats anything out there in term of speed and memory
A = rand(4,10);
B = rand(10,5);
N = size(A,1);
M = size(B,2);
C = zeros(N,M);
for i = 1:N
for j = 1:M
C(i,j) = trapz(A(i,:).*B(:,j).');
end
end
C
C = 4×5
2.3135 3.1579 2.4456 2.5958 1.1064 2.1002 3.4828 2.3818 2.5174 1.2254 1.8511 1.9244 1.6724 1.9377 0.8731 2.9940 3.7654 2.9410 3.9024 1.2115
AA = A; AA(:,[1 end]) = AA(:,[1 end]) / 2;
E = AA*B
E = 4×5
2.3135 3.1579 2.4456 2.5958 1.1064 2.1002 3.4828 2.3818 2.5174 1.2254 1.8511 1.9244 1.6724 1.9377 0.8731 2.9940 3.7654 2.9410 3.9024 1.2115
  1 commentaire
bil
bil le 24 Mar 2023
Great! Looking directly into the integration formula is a good idea, thanks!

Connectez-vous pour commenter.

Plus de réponses (2)

Ashu
Ashu le 23 Mar 2023
Modifié(e) : Ashu le 23 Mar 2023
Hi Bil,
I understand that you want to speedup you code. You can try the following approaches for the same.
arraySize = 1000;
x = ([1:arraySize])';
A = x*(1:arraySize);
B = x./(1:arraySize);
N = 1:size(A,1);
M = 1:size(B,1);
C = zeros(arraySize);
D = zeros(arraySize);
E = zeros(arraySize);
% parallelized loop
tic
parfor i = N
for j = M
C(i,j) = trapz(A(i,:).*B(:,j).');
end
end
Starting parallel pool (parpool) using the 'Processes' profile ...
T1 = toc;
% vectorized
tic
for i = N
D(i,:) = trapz((A(i,:).').*B);
end
T2 = toc;
% simple for loops
tic
for i = N
for j = M
E(i,j) = trapz(A(i,:).*B(:,j).');
end
end
T3 = toc;
Here you can compare the Elapsed Time and see that T1<T2<T3.
To vectorise the operation, you need to understand how 'trapz' works.
If Y is a matrix, then 'trapz(Y)' integrates over each column and returns a row vector of integration values.
That is why while vectorizing the inner for loop, you should transpose A(i,:) and not B.
Now to vectorize the outer for loop you can try to understand the order of operations and move ahead.
To know more about 'trapz', you can refer :
To know more about parallel for loops, you can refer:
Hope it helps!
  3 commentaires
Ashu
Ashu le 23 Mar 2023
Modifié(e) : Ashu le 23 Mar 2023
I have added some more information about vectorization in the answer. I have vectorized the inner for loop, you can do some more analysis of how 'trapz' works with matrices and vectorize the outer for loop.
bil
bil le 24 Mar 2023
Thank you, this response was very helpful!

Connectez-vous pour commenter.


Bruno Luong
Bruno Luong le 23 Mar 2023
Modifié(e) : Bruno Luong le 23 Mar 2023
All loops are removed, but the memory requirement might be an issue.
As I don't know what mean "very large N, M" I can't make adapt the code and make any compromise.
A = rand(4,10);
B = rand(10,5);
N = size(A,1);
M = size(B,2);
C = zeros(N,M);
for i = 1:N
for j = 1:M
C(i,j) = trapz(A(i,:).*B(:,j).');
end
end
C
C = 4×5
2.0640 2.3220 1.5734 1.5456 1.0978 1.8355 2.1440 2.2892 1.8441 1.4358 2.2730 2.5273 2.2265 1.9894 1.8680 3.5998 4.1238 3.5027 3.1071 2.6166
% is equivalent to
D = reshape(trapz(A.*reshape(B,[1 size(B)]),2),[N M])
D = 4×5
2.0640 2.3220 1.5734 1.5456 1.0978 1.8355 2.1440 2.2892 1.8441 1.4358 2.2730 2.5273 2.2265 1.9894 1.8680 3.5998 4.1238 3.5027 3.1071 2.6166

Catégories

En savoir plus sur Numerical Integration and Differentiation 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