How to speed up sparse matrix index operation in Matlab?

Question

I need to create spare matrices with variable elements. Unfortunately, sparse matrix index operations are very slow.

Is there any way to speed up the process? Maybe there are some tricks that I don't know of?

Here is a minimal working example:

N = 400;
A = eye(7)+  [ zeros(3,3) eye(3)     zeros(3,1) ;
               zeros(3,3) zeros(3,3) zeros(3,1) ;
               zeros(1,3) zeros(1,3) zeros(1,1) ] ;
B = [ zeros(3,3) zeros(3,1) ;
      eye(3)     zeros(3,1) ;
      zeros(1,3) -2     ] ;

Su = sparse(7*N, 4*N);

for i=1:N
    for j=0:i-1
        Su(((i)*7 + 1) : ((i+1) * 7), (j*4 + 1) : ...
            ((j+1) * 4)) = A^(i-j-1) * B;
    end
end

Sx = sparse(N*7, N*7);

for i=0:N
    Sx((1 + i*7 : (i+1)*7), (1 + i*7 : (i+1)*7)) =  A^i;
end

Su and Sx are matrices that are partitioned by (products of) A and B.

Answer 1

MATLABs code analyzer already mentions that this method of assigning sparse matrices is very slow as the nonzero pattern is changed which incurs significant overhead. Forgoing the sparse matrices and just using Su = zeros(.. instead makes the function execute in a reasonable time, even if you convert the matrices to sparse matrices afterwards.

If the memory used is of any concern, you could consider defining the zeros matrix as a different type, e.g. int16, as all values appear to be signed integers smaller than 400.

Otherwise, if you can predict the nonzero pattern in you're sparse matrix you can define it correctly upfront and won't suffer performance problems. You don't need to know the exact values, just whether a value is nonzero or not. Taking the following from MATLAB's Code Analyzer:

Explanation

Code Analyzer detects an indexing pattern for a sparse array that is likely to be slow. An assignment that changes the nonzero pattern of a sparse array can cause this error because such assignments result in considerable overhead.

Suggested Action

If possible, build sparse arrays using sparse as follows, and do not use indexed assignments (such as C(4) = B) to build them:

1. Create separate index and value arrays.
2. Call sparse to assemble the index and value arrays.

If you must use indexed assignments to build sparse arrays, you can optimize performance by first preallocating the sparse array with spalloc.

If the code changes only array elements that are already nonzero, then the overhead is reasonable. Suppress this message as described in Adjust Code Analyzer Message Indicators and Messages. For more information, see “Constructing Sparse Matrices”.

So, if you can create an array of indices that will contain nonzero values, you can call sparse(i,j,s[,m,n,nzmax]) to predefine the nonzero pattern. If you only know the number of nonzero values, but not their pattern, using spalloc should also improve your performance.