Speeding up matrix operations

Question

bil 2023-3-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1933785-speeding-up-matrix-operations

评论： bil 2023-3-24

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

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

在 MATLAB Online 中打开

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 2023-3-23

编辑：bil 2023-3-23

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.

请先登录，再进行评论。

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

Answer 1

Bruno Luong 2023-3-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1933785-speeding-up-matrix-operations#answer_1199734

在 MATLAB Online 中打开

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

bil 2023-3-24

Great! Looking directly into the integration formula is a good idea, thanks!

请先登录，再进行评论。

Answer 2

Ashu 2023-3-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1933785-speeding-up-matrix-operations#answer_1198925

编辑：Ashu 2023-3-23

在 MATLAB Online 中打开

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 :

https://www.mathworks.com/help/matlab/ref/trapz.html

To know more about parallel for loops, you can refer:

https://www.mathworks.com/help/matlab/ref/parfor.html

Hope it helps!

3 个评论
显示 1更早的评论隐藏 1更早的评论

Ashu 2023-3-23

编辑：Ashu 2023-3-23

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 2023-3-24

Thank you, this response was very helpful!

请先登录，再进行评论。

Answer 3

Bruno Luong 2023-3-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1933785-speeding-up-matrix-operations#answer_1199645

编辑：Bruno Luong 2023-3-23

在 MATLAB Online 中打开

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