sparse( ix, jx, sx, rx, cx ) implausibly slow

My understanding was that the five argument sparse call was designed to be a fast way of creating sparse matrices. But in some Kronecker product code I'm using (full code at the bottom, code indirectly derived from here), the profiler consistently reports it as the bottleneck.

Bizarrely, even doing sortrows( [ix,jx,sx] ) first does not speed up the call to sparse, despite the fact that this is basically the internal storage representation for sparse matrices used by Matlab.

Do you have any suggestions for speeding up the call to sparse?

Thanks in advance,

Tom

-------------------------------------------------------------------------------------------------

Suggested test:

 n=300;A=sprandn(n,n,0.1);B=sprandn(n,n,0.1);
 profile on;X=spkron(A,B);profile off;profile viewer

Required functions:

 function X = spkron( A, B )

    global spkron_use_mex

    [I, J] = size(A);
    [K, L] = size(B);

    [ia,ja,sa] = find( A );
    [ib,jb,sb] = find( B );

    a = double( [ia,ja,sa] );
    b = double( [ib,jb,sb] );

    if isempty( spkron_use_mex )
        [ ix, jx, sx ] = spkron_internal( K,a, L,b );
    else
        [ ix, jx, sx ] = spkron_internal_mex_mex( int32(K),a, int32(L),b );
    end
    X = sparse( ix, jx, sx, I*K, J*L );

 end

 function [ ix, jx, sx ] = spkron_internal( K,a, L,b )

    % derived from alt_kron.m

    ma = max( abs( a(:,3) ) ) * eps;
    mb = max( abs( b(:,3) ) ) * eps;

    a( abs(a(:,3))<mb, : ) = [];
    b( abs(b(:,3))<ma, : ) = [];

    ix = bsxfun(@plus,b(:,1),K*(a(:,1)-1).');
    jx = bsxfun(@plus,b(:,2),L*(a(:,2)-1).');

    sx = bsxfun(@times,b(:,3),a(:,3).');
    sx( abs( sx ) < eps ) = 0;

 end

[Matt J]

In your suggested test, your "sparse" matrices are really rather dense (0.1). If this density is representative of the work you will be doing, and if the Kronecker products are going to be used for matrix-vector multiplications and mldivides, then you will get more speed using KRONPROD, even excluding the time to do the kron operation. Compare:

    n=300;
    A=sprandn(n,n,0.1);
    B=sprandn(n,n,0.1);
    x=rand(n^2,1);

      K=kron(A,B); 
      tic 
       y=K*x;
      toc
      %Elapsed time is 0.179062 seconds.

      tic;
       Ko=KronProd({B,A});
       y=Ko*x;
      toc;
      %Elapsed time is 0.008858 seconds.

[Tom Holden]

Don't worry, it's not representative... I am aware of your (very nice) library though, and it's actually being distributed with a related project to the one I'm currently working on. (See: https://github.com/lanhongken/nlma).

Indeed, I currently have a comment in one bit of my code that reads something like "TODO: convert this to use the KronProd class", in an area of the code where the matrices are denser.

Answer 1

[Chris Turnes]

Using 5 input arguments with the sparse command is indeed an efficient way to construct a sparse matrix. However, you most likely will not be able to speed up the call to sparse, because no matter how the matrix is constructed there is some overhead cost in storing it.

Specifically, MATLAB uses an efficient storage format for sparse matrices. Rather than store the column indices of the non-zero elements, sparse matrices instead store a vector that holds, for each column, the total number of nonzero elements in all preceding columns. The purpose of this is to reduce the overall storage, as in most cases the total number of non-zeros is greater than the number of columns. However, it also means that sparse has to compute the entries of this vector, which introduces overhead.

The suggested test for your code demonstrates quite well why this system is beneficial. When I run the test, I produce the following sparse matrix X:

>> whos X
  Name          Size                    Bytes  Class     Attributes

    X         90000x90000            1175152232  double    sparse

The memory used by X is given by the formula

>> int32(2*8*nnz(X) + (size(X,2)+1)*8)
ans =
  1175152232

The first term is the cost of storing the values of the entries and their row indices, and the second is the cost of storing the column data. Compare this with the cost of storing the both the column and row indices:

>> int32(3*8*nnz(X))
ans =
  1761648336

So MATLAB's format saves 586 MB of memory, but it comes at the cost of a more time-consuming initialization.

[Matt J]

Tom Holden Commented:

Very clear, thanks. I had misunderstood the internal storage format for sparse matrices it seems.