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

3 views (last 30 days)
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?

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 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 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...

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);

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 on 6 Dec 2019
If that's the case, my intuitions says that your current approach is probably optimal.