Summation of dense and sparse vectors

Question

I need to sum up a dense and a sparse vectors in Matlab, and the naive way of doing it, that is:

w = rand(1e7,1);
d = sprand(1e7,1,.001);
tic
for k = 1 : 100
    w = w + d;
end
toc

takes about 3.5 seconds, which is about 20 times slower then the way I'd expect Matlab to implement it behind the hood:

for k = 1 : 100
    ind = find(d);
    w(ind) = w(ind) + d(ind);
end

(of course the timing of this faster version depends on the sparsity).

So, why doesn't Matlab do it 'the fast way'? My experience with Matlab till now suggests that it is quite good at utilizing sparsity.

More important, are there other 'sparse' operations I should suspect as being not efficient?

Answer 1

I don't know the answer for sure, but I will give you my guess about what is happening. I don't know Fortran, but from a C++ perspective, what you are showing makes sense, we just need to deconstruct that statement.

The pseudo-code translation of a = b + c where a,b are full and c is sparse would look something like a.operator= ( b.operator+ (c) ).

In all likelihood, full matrix containers in Matlab should have specialised arithmetic operators to deal with sparse inputs, ie something like full full::operator+ ( const sparse& ). The main thing to notice here is that the result of a mixed full/sparse arithmetic operation has to be full. So we're going to need to create a new full container to store the result, even though there are few updated values. [ Note: the returned full container is a temporary, and so the assignment a.operator= ( ... ) may avoid a full additional copy, eg with full& full::operator= ( full&& ). ]

Unfortunately, there is no way around returning a new full container because there is no arithmetic compound operations in Matlab (ie operator +=). This is why Matlab cannot exploit the fact that in your example a is the same as b (try to time your loop with x = w + d instead, there's no runtime difference), and this is where the overhead comes from IMO. [ Note: even when there is no explicit assignment, eg b+c;, the generic answer variable ans is assigned. ]

Interestingly, there seems to be a notable difference between full full::operator+ ( const sparse& ) and full sparse::operator+ ( const full& ), that is between a = b + c and a = c + b; I can't say more about why that is, but the latter seems to be faster.

In any case, my short answer is "because Matlab doesn't have arithmetic compound operators", which is unfortunate. If you know you are going to do these operations a lot though, it shouldn't be too hard to implement ad hoc optimized versions like the one you proposed. Hope this helps!