Composition of two matrices with the same number of columns, one of them is integer-valued

Question

Quang Huy Pham 2021-12-9

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1606950-composition-of-two-matrices-with-the-same-number-of-columns-one-of-them-is-integer-valued

编辑： Jan 2021-12-11

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

What are typical sizes of n, m and d? Dozens or billions?

Quang Huy Pham 2021-12-9

These matrices are not big: d is dozens, n is thousands, m is less than 10. But I need to run the code repeatedly millions of times for simulation. The index matrix B is fixed throughout, A changes over time.

请先登录，再进行评论。

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

Answer 1

Jan 2021-12-9

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1606950-composition-of-two-matrices-with-the-same-number-of-columns-one-of-them-is-integer-valued#answer_851615

编辑：Jan 2021-12-11

在 MATLAB Online 中打开

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
Elapsed time is 0.601715 seconds.
% 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
Elapsed time is 0.170638 seconds.
% 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
Elapsed time is 0.114302 seconds.
% 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
Elapsed time is 0.073988 seconds.

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
Elapsed time is 0.072217 seconds.
% [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
Elapsed time is 0.121049 seconds.
% Version 5 - linear indexing, no loop, but large intermediate arrays:
tic;
prod_C5 = prod(A(B + (0:d-1)*m), 2);
toc
Elapsed time is 0.182863 seconds.
% Are all results equal?
isequal(prod_C, prod_C2, prod_C3, prod_C4, prod_C5)
ans = logical
   1