Tools used for Sparse Matrix Efficiency

I have read documentation and the MATLAB community answers regarding this and I found a lot of helpful information about the power of sparse matrices, I just need help applying it to my sceanrio as this is the first time I am working with them.

I currently have a 22,400 x 22,400 sparse identity matrix that I created with speye(n) where n= 24,000. In a loop I have multiple element assignments that I do which iterates n number of times. The total run time for this code is 13.811054 seconds which is very poor considering this being a placeholder, and I plan on running 10e6 x 10e6 size matrices and I appreciate pointers to have more efficient code.

I ran multiple trials with on how the size affects my speed, along with time it took I calculated how much elements were allocated to the matrix (non-zero elements and I tabulated 3 readings below.

Size Time elapsed Allocation

576 x 576 0.128121 seconds 1.042%

2,800 x 2,800 0.657801 seconds 0.2143%

22,400 x 22,400 13.811054 seconds 0.0268%

The skeleton of my code is shown below;

A = -6.*speye(n); %intializaiton of a nxn sparse matrix

for ix = 1:X
    for iy = 1:Y 
        for iz = 1:Z 
   
            %multiple calculations and conditional statements
            
            A(icent,ixp)=1;
            A(icent,iyp)=1;
            A(icent,izp)=1;
            A(icent,ixm)=1;
            A(icent,iym)=1;
            A(icent,izm)=1;
            %allocations to the sparse matrix (poorly done)
                        
        end
    end
end

l wish to know what is the best to improve the efficency of this, I understand how clunky and inefficent my element allocation is and I want to improve it, any pointers are appreciated.

While I was reading answers on MATLAB, I came across John D'Errico's Answer to a similar question which dealt with using spalloc, which was well-worded, however I can't seem to understand how to apply spalloc() to my sceanrio where I already require a declared sparse identity matrix.

I hope to have an efficient code to run large matrix sizes, appreciate the help in advance.

**Edit - the 1s can be allocated only inside the nested loop as there are prior calculations that need to be done, I realize John D'Errico's Answer mentions in the beginning - In general, DON'T DO IT! That is, don't build a sparse matrix one element at a time. Even if you use spalloc, the matrix will be inefficient to build.

Edit - I am working with a symmetric banded matrix, the spine of it will be the speye() identity matrix

[Walter Roberson]

What is the relationship between ix, iy, iz, and icent, ixp, iyp, izp, ixm, iym, izm ? Are the icent and other variables constants inside the loops so you are just assigning to the same place over and over again?

[Walter Roberson]

If you are creating a symmetric banded matrix, then we would expect to see something more like

            A(icent,ixp)=1;
            A(icent,iyp)=1;
            A(icent,izp)=1;
            
            A(ixp,icent)=1;
            A(iyp,icent)=1;
            A(izp,icent)=1;

[Hans123]

^Thanks as always, Walter Roberson

Answer 1

[Walter Roberson]

used_count = ceil(n + X*Y*Z*6);
A = spalloc(n, n, used_count);
for i = 1 : n; A(i,i) = -6; end
for ix = 1:X
    for iy = 1:Y 
        for iz = 1:Z 
            A(icent, [ixp, iyp, izp, ixm, iym, izm]) = 1;
        end
    end
end

Key points here:

• spalloc() specifying non-zero count as the maximum number of non-zero elements that will ever be stored in the array
• If you assign a complete sparse matrix to another variable, the non-zero count is preserved, as that is the standard copy-on-write situation. However, if you do arithmetic with a sparse matrix, the resulting sparse expression will have all of its extra non-zero count removed. Thus, although it is tempting to do the A=spalloc followed by A=A-6*speye(n) unfortunately the sparse matrix calculated on the right between A and something else will trim out the slack and A would no longer have room to grow. This is why I specifically loop assigning into the diagonal
• merge operations when possible, the sparse overhead is not small

A lot of the time, it is much more efficient to create vectors of row and column indices and vectors of associated values, and sparse() the entire array into place at one time.

Answer 2

[John D'Errico]

Whoever gave you that advice, about never building a sparse matrix one element at a time must have been very knowledgable, as it was totally on-target. Oh. Right. That was me. :)

If you seriously don't know where the elements will end up, then build them in a loop, storing them in an array, thus row, column, value. If you know how many final elements there will be in that array, of even an upper limit on the count, then preallocate the array to contain them, as perhaps an Nmax by 3 array. Of course you do not want to grow the size of that array, as that itself is a CPU intensive thing to do.

In that case, you would best use one of the tools I have posted on the FEX, to save elements in a grown array, but to do so efficiently.

https://www.mathworks.com/matlabcentral/fileexchange/8334-incremental-growth-of-an-array-revisited

Once all of the non-zero elements have been created, create the array using one call to sparse.

[Hans123]

Thank you, John. Appreciate your insight into this problem

Answer 3

[Matt J]

Do not update the matrix A iteratively inside a for-loop. Just save the coordinate pairs and then use them after the loop to build A in its entirety:

[I,J]=deal( nan(6,X*Y*Z));

count=0;
for ix = 1:X
    for iy = 1:Y 
        for iz = 1:Z 
   
            %multiple calculations and conditional statements
  
             count=count+1;
            
             I(:,count)=icent;
             J(:,count)=[ixp,iyp,izp, ixm, iym, izm];
                        
        end
    end
end

I(:, any(isnan(I),1)  )=[];
J(:, any(isnan(J),1)  )=[];

A = -6.*speye(n) + sparse(I(:),J(:),1,n,n);