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

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?

[JESUS DAVID ARIZA ROYETH]

please attach the dependency of Ri also

[Tudor]

I updated the question with some additional details.

[Matt J]

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]

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]

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

Answer 1

[Matt J]

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]

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]

If that's the case, my intuitions says that your current approach is probably optimal.