MATLAB: Sparse( ix, jx, sx, rx, cx ) implausibly slow

sparse

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

Best Answer

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