Speeding up matrix operations

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.

[Cameron]

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]

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]

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

[bil]

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

Answer 2

[Ashu]

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!

[bil]

Hi Ashu,

Thanks for the response. i do agree this would help but I was wondering if there was some cleaner way that could somehow do away with the for loops entirely?

[Ashu]

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]

Thank you, this response was very helpful!

Answer 3

[Bruno Luong]

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