Sparse Matrix - More Efficient Assignment Operation

I have a large sparse matrix for which I am attempting to assign certain segments (simplified example as per below snippet) to some identity matrices. The assignment method below is extremely inefficient, taking up a large amount of memory and too much time (10sec for example below on my server, but my real dataset is much bigger). Looking at the screenshot from taskmanager, we can see clearly a large spike in memory usage while the operation is being carried out.

For somebody used to dealing with large sparse matrices, there should be a more efficient way of carrying out such an assignment operation. Probably a pretty simple solution but I have little experience with sparse-matrix syntax.

A = sparse(20000,60000);
A(10001:end,20001:50000) = ([speye(10000,10000),-speye(10000,10000),speye(10000,10000)]);

Answer 1

[the cyclist]

That operation took less than 3 milliseconds on my machine.

The resulting array is 960,008 bytes.

[Mark Whirdy]

Hi there

Sure enough on my local box (2012b much lower spec), it takes 3ms. Windows Server takes 10sec with 64bit 32G R2010b. My colleagues find the same.

What version of matlab are you on?

[the cyclist]

MATLAB R21013a on a Macbook Pro running Mountain Lion

[Mark Whirdy]

We ran on linux server & got 15sec on a 2009 matlab. Maybe there have been some improvements in sparse handling between versions at some point - any mathworks developers reading this have any insights?

[Sean de Wolski]

Hi Mark,

I can run the timings for you on any and all releases when I return to Natick on Friday.

Answer 2

[Cedric]

The structure of sparse matrices in memory makes this kind of indexing operations slower than building the sparse matrix directly using vectors of row/col IDs and values.

 I = [10001:20000, 10001:20000, 10001:20000] ;
 J = [20001:50000] ;
 V = [ones(1, 1e4), -ones(1, 1e4), ones(1, 1e4)] ;
 S = sparse(I, J, V, 2e4, 6e4) ;
 spy(S) ;    % Check structure.

You can easily make it more flexible/concise. I can't test it now though.

[Mark Whirdy]

Hi Cedric, thanks for your answer. You are right in a general sense, but for my particular use-case the elements of matrix A need to be assigned different values throughout the file so I think that re-allocating memory continuously for an entire S (by using S = sparse(I, J, V, 2e4, 6e4); for each assignment) will end up partitioning the memory etc. Maybe I'm wrong.

It seems like there's maybe just some version issue.

Thanks, All the best, Mark

[Jan]

Creating the 4(!) intermediate arrays wastes time:

[speye(10000,10000), -speye(10000,10000), speye(10000,10000)]

The following copy to a submatrix needs more time also than setting the values directly and pre-allocating the number of expected values. Therefore I assume Cedric's asnwer should be much more efficient.

[Cedric]

I can't answer for your case specifically as I don't know the details. On my side though, I almost always build index/value arrays during the "build phase" of the sparse matrix content, or at least until I need the sparse matrix for computation purpose. Even when I have some loop that defines blocks of a final matrix, I build blocks of indices actually. The reason is the column compressed internal structure of sparse matrices, that makes block-indexing inefficient (to some extent).

[Mark Whirdy]

ok, thanks guys

[James Tursa]

Slight change to Cedric's method to get the size as desired (added a 0 element at the end):

 I = [10001:20000, 10001:20000, 10001:20000 20000] ;
 J = [20001:50000 60000] ;
 V = [ones(1, 1e4), -ones(1, 1e4), ones(1, 1e4) 0] ;

Also, there IS a version specific issue here. The original code OP posted generates an out-of-memory on my Win32 XP machine through R2010b, but is fine for R2011a onwards. Which leads me to believe that there is something in the indexing or concatenating that changed between these versions with regards to sparse matrices and how they use temporary memory in the background. The final result isn't that large.

Regarding OP's comment that he needs to insert different values (other than 1's and -1's I presume), we would need to see some specifics in order to offer any suggestions.

[Cedric]

It would be interesting to test with a smaller size matrix that would fit in memory as a dense matrix. There might be some conversion to dense in particular cases of indexing that could be spotted this way. I've been using quite intensely sparse matrices in MATLAB since ~2006 I guess, and I already had troubles with unsuspected conversions to dense a few times (but not while CAT-ing sparse matrices, as it is not an operation that I am doing frequently).

Answer 3

[Teja Muppirala]

The performance of sparse matrix indexing was enhanced in R2011a.

http://www.mathworks.com/help/matlab/release-notes.html

If you can't upgrade, as somewhat of a workaround, you should be able to get away with something like this:

A = sparse(20000,60000);
A = spreplace(A,10001:20000,20001:50000,[speye(10000,10000),-speye(10000,10000),speye(10000,10000)]);

Where SPREPLACE is the following general purpose function:

function A = spreplace(A,I,J,B)
% Equivalent to 
% >> A(I,J) = B
% But does not support "end" indexing in I and J

I = I(:);
J = J(:);

[iA,jA,sA] = find(A);
[iB,jB,sB] = find(B);

trimA = ~(ismember(iA,I) & ismember(jA,J));

A = sparse([iA(trimA); I(iB)],...
           [jA(trimA); J(jB)],...
           [sA(trimA); sB],...
            size(A,1),size(A,2));

When I try this in R2007a, Your original code takes 15 seconds, The workaround above gets it done in about 5 milliseconds.