Speeding up matrix operations
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
Cameron
2023-3-23
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.
采纳的回答
Bruno Luong
2023-3-23
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
AA = A; AA(:,[1 end]) = AA(:,[1 end]) / 2;
E = AA*B
更多回答(2 个)
Ashu
2023-3-23
编辑:Ashu
2023-3-23
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
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 个评论
Bruno Luong
2023-3-23
编辑:Bruno Luong
2023-3-23
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
% is equivalent to
D = reshape(trapz(A.*reshape(B,[1 size(B)]),2),[N M])
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!