How to vectorize the following code for matrix insertion?

Question

Jordan Tryggvesson 2021-5-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/827085-how-to-vectorize-the-following-code-for-matrix-insertion

评论： Jordan Tryggvesson 2021-5-13

I have an nxn matrix A and a 3xm array g. In general, n will be large, so A will be pre-allocated,

I want to carry out the following insertion, where I pick index pairs from the first two rows of g and values from the last row of g, in order to form a 2x2 matrix which I insert in A:

for i = 1:m
    A(g(1:2,i),g(1:2,i)) = g(3,i)*[1 2; 3 4];
end

Note that [1 2; 3 4] is some fixed matrix of different size from A. The purpose of the above code is to insert the matrix [1 2; 3 4] , scaled by g(3,i), into the matrix A at position A(g(1:2,i),g(1:2,i)).

Is there any way to vectorize this code? For example, I tried the following (ugly) solution, but it doesn't work:

A(g(1:2,:),g(1:2,:)) = [1*g(3,:) 2*g(3,:);3*g(3,:) 4*g(3,:)];

Any suggestions?

4 个评论
显示 2更早的评论隐藏 2更早的评论

J. Alex Lee 2021-5-11

You might look into "sparse", which has well-defined syntax for making such insertions, not even limited to contiguous sub-matrices. No idea if it will be faster though

Jordan Tryggvesson 2021-5-11

编辑：Jordan Tryggvesson 2021-5-11

Alright, thanks for the suggestion!

请先登录，再进行评论。

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

Answer 1

J. Alex Lee 2021-5-11

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/827085-how-to-vectorize-the-following-code-for-matrix-insertion#answer_697720

编辑：J. Alex Lee 2021-5-11

在 MATLAB Online 中打开

There is a minor point in David's answer, which is that not all index collisions are taken care of in how the g matrix is constructed...

So the "sparse" method I suggested, which is basically the same in spirit as David's answer except it doesn't need sub2ind and assumes we want sparse, will give a different answer because instead of overwriting repeated index (as both above methods do), using sparse with repeated subscripts will add the values that are repeated.

So slightly modifying the original loop version to be additive (I know it changes the rules, but just to keep this answer simpler), you can compare timing. On my computer the sparse method is faster as m increases. I observe that David's vectorized answer does see some speed gain versus looping, but advantage scales differently than the sparse method...I'm speculating now, but I think that's because the bottle neck is indexing into full matrices...but i'm no computer scientist so I'll let better minds than mine chime in on that.

% make up some data using David's answer
N = 1000;
m = 200;
g = randi(N,3,2*m);
% eliminate instances where g(1:2,k) is a repeated index
irep = (g(1,:)-g(2,:))==0;  
g(:,irep) = [];
g = g(:,1:m);
% original, but modified to add instead of overwrite
tic
A = zeros(N);
for k = 1:m
    A(g(1:2,k),g(1:2,k)) = A(g(1:2,k),g(1:2,k)) + g(3,k)*[1,2;3 4];
end
toc
% using sparse
tic
ga = g(1:2,:);
sub1 = [ga(1,:); ga; ga(2,:)];
sub2 = [ga; ga];
A2 = sparse(sub1,sub2,[1 2 3 4]'*gb,N,N);
toc
err = max(abs(A-A2),[],'all')  % should be zero

A better apples-to-apples comparison would be if the matrix g was constructed such that all index collisions are avoided, or if the original application doesn't prohibit that, need some clarification on intended behavior.

2 个评论
显示无隐藏无

J. Alex Lee 2021-5-11

在 MATLAB Online 中打开

Here's a minimal example that I had to do for myself to understand why my method was giving different answers

C = zeros(1,6); C([1,2,3,2]) = [1 2 3 4]
C = 1×6
     1     4     3     0     0     0

vs

full(sparse([1 1 1 1],[1,2,3,2],[1 2 3 4],1,6))
ans = 1×6
     1     6     3     0     0     0

Jordan Tryggvesson 2021-5-13

Thank you very much, and sorry for my late reply!

请先登录，再进行评论。

Answer 2

David Goodmanson 2021-5-11

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/827085-how-to-vectorize-the-following-code-for-matrix-insertion#answer_697650

编辑：David Goodmanson 2021-5-11

在 MATLAB Online 中打开

Hi Jordan,

Here is one method. Whether or not it is faster than the for loop method is another question. The code does exactly the same as your code, inserting elements into A. Consequently any elements of A that are nonzero (from previous actions of the for loop) get overwritten. The elements are not added.

% make up some data
N = 100;
A = zeros(N,N);
m = 88;
g = randi(N,3,2*m);
% eliminate instances where g(1:2,k) is a repeated index
irep = (g(1,:)-g(2,:))==0;  
g(:,irep) = [];
g = g(:,1:m);
% original  
for k = 1:m
    A(g(1:2,k),g(1:2,k)) = g(3,k)*[1 2; 3 4];
end
% vectorized
A1 = zeros(N,N);
ga = g(1:2,:);
gb = g(3,:);
sub1 = [ga(1,:); ga; ga(2,:)];
sub2 = [ga; ga];
ind = sub2ind([N N],sub1,sub2);
A1(ind) = [1 2 3 4]'*gb;
max(abs(A-A1),[],'all')  % should be zero