Composition of two matrices with the same number of columns, one of them is integer-valued
1 view (last 30 days)
Show older comments
I have two matrices A, B with the same number of columns (d). The entries of B are positive integers. I want to create a matrix C defined by C(i,j)=A(B(i,j),j) and find the products of the rows of C.
I tried the for loop below but it is quite slow. Is there any MATLAB operation to compute this kind of composition faster?
B = zeros(n,d);
m = max(B(:));
A = zeros(m,d);
C = zeros(n,d);
for ii = 1:n
for jj = 1:d
C(ii,jj) = A(B(ii,jj),jj);
end
end
prod_C = prod(C,2);
2 Comments
Accepted Answer
Jan
on 9 Dec 2021
Edited: Jan
on 11 Dec 2021
n = 1e4;
m = 1e4;
d = 1e3;
B = randi([1, m], n, d);
A = randn(m, d) * 2;
% Version 1: The original code
tic
C = zeros(n, d);
for ii = 1:n
for jj = 1:d
C(ii,jj) = A(B(ii,jj),jj);
end
end
prod_C = prod(C,2);
toc
% Version 2: Change the order of loops
tic
C = zeros(n, d);
for jj = 1:d
for ii = 1:n
C(ii,jj) = A(B(ii,jj),jj);
end
end
prod_C2 = prod(C,2);
toc
% Version 3: Omit the inner loop
tic
C = zeros(n, d);
for jj = 1:d
C(:, jj) = A(B(:, jj), jj);
end
prod_C3 = prod(C, 2);
toc
% Version 4: Multiply on the fly (no large intermediate array!)
tic
prod_C4 = ones(n, 1);
for jj = 1:d
prod_C4 = prod_C4 .* A(B(:, jj), jj);
end
toc
10 times faster finally.
A nit-picking improvement of version: 4
% Version 4b: Should not change the speed measureably
tic
prod_C4 = A(B(:, 1), 1); % Save one multiplication
for jj = 2:d % Start at jj=2
prod_C4 = prod_C4 .* A(B(:, jj), jj);
end
toc
% [EDITED] 4c: Mulitply on the fly in a loop:
tic
prod_C4 = ones(n, 1);
for jj = 1:d
for ii = 1:n
prod_C4(ii) = prod_C4(ii) * A(B(ii,jj),jj);
end
end
toc
% Version 5 - linear indexing, no loop, but large intermediate arrays:
tic;
prod_C5 = prod(A(B + (0:d-1)*m), 2);
toc
% Are all results equal?
isequal(prod_C, prod_C2, prod_C3, prod_C4, prod_C5)
Choose version 4.
[EDITED] In my local Matlab 2018b, version 4c is 30% faster than version 4.
2 Comments
More Answers (0)
See Also
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!