Most efficient way to add multiple sparse matrices in a loop in MATLAB

I have a code that repeatedly calculates a sparse matrix in a loop (it performs this calculation 13472 times to be precise). Each of these sparse matrices is unique.

After each execution, it adds the newly calculated sparse matrix to what was originally a sparse zero matrix.

When all 13742 matrices have been added, the code exits the loop and the program terminates.

The code bottleneck occurs in adding the sparse matrices. I have made a dummy version of the code that exhibits the same behavior as my real code. It consists of a MATLAB function and a script given below.

(1) Function that generates the sparse matrix:

function out = test_evaluate_stiffness(n)

ind = randi([1 n*n],300,1);

val = rand(300,1);

[I,J] = ind2sub([n,n],ind);

out = sparse(I,J,val,n,n);

end

(2) Main script (program)

% Calculate the stiffness matrix

n=1000;

K=sparse([],[],[],n,n,n^2);

tic

for i=1:13472

      temp=rand(1)*test_evaluate_stiffness(n);

      K=K+temp;

end

fprintf('Stiffness Calculation Complete\nTime taken = %f s\n',toc)

I'm not very familiar with sparse matrix operations so I may be missing a critical point here that may allow my code to be sped up considerably.

Am I handling the updating of my stiffness matrix in a reasonable way in my code? Is there another way that I should be using sparse that will result in a faster assembly?

Note that I have coupling in my system so the the values by which I augment my global stiffness matrix in each loop execution aren't necessarily on the diagonal.

A profiler report is also provided below:

%

Answer 1

[John D'Errico]

The problems is? You are doing it wrong. Flat out, wrong.

NEVER build a sparse matrix up like this. NEVER. NEVER. If I say it 3 times, it must be true.

The problem is that every time you add new entries to the matrix, it must move those elements around in memory. EVERY time. Plus, there are more elements that are non-zero in the final matrix for every iteration, so it gets a little larger every time. So more memory must be allocated. This is the classic problem of dynamic allocation.

The answer? Do it the right way! There is a far better way to solve this problem.

Don't create a sparse matrix on every iteration. Instead, just store the non-zero elements. I would suggest storing them in a simple flat array, with 3 columns for each iteration. Since you know there will be 13472 sub-arrays, save them as cells of a cell array.

C = cell(13472,1);

Now, on each iteration, stuff one cell with the information needed to create the sub-array, thus

   for i = 1:13472
    .... stuff ...

    C{i} = [I,J,Val];
  end

Now, AFTER the loop is done, concatenate all of those cells into ONE array.

IJV = cell2mat( C );

Finally, ONLY now do you want to call sparse. CALL IT ONE TIME. I'll assume the final size of the array is known to be n by m.

A = sparse(IJV(:,1),IJV(:,2),IJV(:,3),n,m);

In the event that you do not know how many sub-arrays you wish to sum, there are still good ways to solve the problem, accumulating the information in cell array elements. You can find tools to do this on the file exchange, in the form of my growdata tools.

Learn to build sparse matrices efficiently. There are god ways, and there are obscenely bad ways to do almost anything.

[Kobye]

John, you're absolutely correct of course. Thanks!

[John D'Errico]

I've been building sparse matrices for many years now, and I've made all possible mistakes already. Well, I think I have. :)

[Xavier DaNNEELS]

Thank you for this John! Helped me out too!

[Stephen23]

@John D'Errico: how does this add the values if two cells address the same location in the sparse array?

[Zoltán Csáti]

@Stephen Cobeldick: the sparse command has this assumption. See: https://nl.mathworks.com/help/matlab/ref/sparse.html#bul62_1

Answer 2

[James Tursa]

Unless you know in advance the pattern of the non-zero spots and can exploit that, then I think you are stuck. At every iteration of K = K + temp, the entire data area of K in general will need to be copied to a new memory location and the temp data will need to be added/merged in as well.

If you do know the pattern in advance, then possibly something could be done to minimize the data copying.

How big is n in your real problem? Can you do some calculations in full and then convert to sparse later for your downstream needs?

[Kobye]

James, in my real problem n = 7910. Thus my "temp" matrices are 7910x7910. Each has about 8000 nonzero elements. However, while these elements may be off-diagonal, they generally do have bandedness. This is demonstrated in the following image which was generated using spy(K) after running the real program (note that the dummy code won't give a banded K matrix, although I find the same ratio of time is taken by the K+temp operation in that code, about 95% of the total computation effort).

With regards to performing "full" operations, that massively slows down the program and instead of 70 seconds it will take almost ten times as long!

[James Tursa]

Two strategies to consider:

(1) Have K full and keep temp sparse. If MATLAB is smart enough, should be able to add temp to K without turning temp full first. If MATLAB is not smart enough a mex routine could be written to do this in-place. Drawback is that K memory is still accessed in fragmented fashion when adding temp, and K is very large, so there may not be any benefit to this approach.

(2) Create K as a banded sparse matrix in the pattern that you know the result to be, putting realmin values in the banded spots. The spdiags function can help with this. Then, if MATLAB is smart enough (no guarantee of this), the K = K + temp will not require any data copying. If MATLAB is not smart enough for this, a mex routine could easily be written to do the K = K + temp calculation in-place to avoid the data copy. Any leftover realmin values could be cleaned out at the end.

UPDATE: I just did a test and MATLAB is not smart enough to avoid the data copy in (2) above. So your only option for this method would be a mex routine.

[Kobye]

James, I don't know if it is allowed on these forums to link external websites, but I asked this question on another forum and a user there was able to speed up my code by 27x! The solution was to not add the temp matrix in each loop execution, but save I,J,val in an array and append that to an overall I,J,V array for each execution. The sparse(I,J,V) command is then only used once the loop has been exited, and automatically adds repeated entries. This achieved my goal of massively speeding up this piece of code!

[John D'Errico]

See my answer. There are good ways to do this without resorting to mex.