Composition of two matrices with the same number of columns, one of them is integer-valued
6 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Jan
2021 年 12 月 9 日
編集済み: Jan
2021 年 12 月 11 日
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 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Loops and Conditional Statements についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!