MATLAB Answers

Most efficient way to enter values into pre-allocated sparse matrix?

3 views (last 30 days)
Tudor
Tudor on 6 Dec 2019
Edited: Matt J on 6 Dec 2019
I have a problem where I have a sparse matrix of a specific size, 4000 x 4000. The problem also has a tri-diagonal structure such that I know which elements of the sparse matrix will be non-zero. It is blocks of 4 x 4 along the main diagonal, and the first of-diagonals.
My issue is that I have to compute the 4 x 4 block iteratively, and store them in the matrix as I compute them.
A = spalloc(4000,4000,4*4*1000*3); % Sparse matrix allocation
for i = 1:1000
idx = (1:8)+4(i-1); % In each iteration, the indices that change are known
A(idx,idx) = A(idx,idx) + Ri; % The matrix Ri depends on i
end
Exactly how the matrix that is added, Ri, depend on the interation is a bit complicated to explain. It depends on an outside process, but it has the following structure
% Either
Ri = [G B; B.' C];
% or
Ri = [D E; E.' F];
Which one it is in each iteration depends on a random process that is (and has to be) randomly sampled in each iteration.
From what I have read about sparse matrices in Matlab, it appears to be important how sparse matrices are constructed, where a bad way may take much longer than a better way. So, does anyone known the best way to do this?

  5 Comments

Show 2 older comments
Matt J
Matt J on 6 Dec 2019
Ri = [A B; B.' C];
The matrix A in this expression is presumably not the same as the A representing the final sparse matrix that you are building.
Tudor
Tudor on 6 Dec 2019
You are right, it is not. Thanks for pointing that out, I should have been more careful when I wrote that. I edited my question.
Matt J
Matt J on 6 Dec 2019
Tudor's comment moved here:
I did some testing, and found that
A(idx,idx) = full(A(idx,idx)) + Ri;
is actually a bit faster than
A(idx,idx) = A(idx,idx) + Ri;
However, I don't understand why that is...

Sign in to comment.

Accepted Answer

Matt J
Matt J on 6 Dec 2019
Edited: Matt J on 6 Dec 2019
I would just store all the data from the loop calculations in cells. Then use the data to build the sparse matrix after the loop.
[Icell,Jcell,Rcell]=deal(cell(1000,1));
for i = 1:1000
[Icell{i},Jcell{i}]=ndgrid((1:8)+4(i-1));
Icell{i}=Icell{i}(:);
Jcell{i}=Jcell{i}(:);
Rcell{i}=Ri(:); % The matrix Ri depends on i
end
I=cell2mat(Icell);
J=cell2mat(Jcell);
R=cell2mat(Rcell);
A=sparse(I,J,R,4000,4000);

  2 Comments

Tudor
Tudor on 6 Dec 2019
Thanks for the suggestion! A challenge is that in each iteration, I need the part of A that has been constructed so far to solve a least squares problem of the type
A\b % b is a sparse vector
Do you know if your suggestion is efficient if one has to construct the sparse matrix in each iteration?
Matt J
Matt J on 6 Dec 2019
If that's the case, my intuitions says that your current approach is probably optimal.

Sign in to comment.

More Answers (0)

Sign in to answer this question.


Translated by