accumarray is faster than sparse

Alec Jacobson

June 10, 2019

weblog/

Optimizing matlab code is a puzzle-solving game I've come to really enjoy. The language is such a mess that there's always interesting, unusual ways to squeeze out extra performance. Recently, I discovered that accumarray is faster than sparse for map-reduce paradigms. This applies pretty generally. If you're using sparse to map-reduce into a vector, then it's a straightforward switch. The only different seems to be that accumarray will output a full (aka dense) vector and it needs column-vector inputs:

u_s = sparse(I,1,V);
% is equivalent to:
vec = @(X) reshape(X,[],1);
u_a = accumarray(vec(I),vec(V));

Here's another case that is particularly useful. Suppose you have a calculation c = A * b where A is a sparse matrix, and you're only building it so that you can multiply it once against b. (This happens quite a bit in gradient-based optimization on FEM meshes). Using accumarray you can skip constructing A entirely, and directly compute c:

c_s = sparse(I,1,V) * b;
% is equivalent to:
vec = @(X) reshape(X,[],1);
c_a = accumarray(vec(I),vec(V).*vec(b(J)));

This can be significantly faster. Here's a little timing script:

m = 1000000;
n = 100000;
nnz = 10*m;
I = ceil(m*rand(nnz,4));
J = ceil(m*rand(nnz,4));
V = rand(nnz,4);
A = sparse(I,J,V);
vec = @(X) reshape(X,[],1);
b = rand(m,1);
max(abs( sparse(I,J,V)*b - accumarray(vec(I),vec(V).*vec(b(J))) ))
timeit(@() sparse(I,J,V)*b)
timeit(@() accumarray(vec(I),vec(V).*vec(b(J))))

On my laptop this prints:

ans =
   1.7764e-14
ans =
       13.357
ans =
       1.5403

More than an 8x speedup for a simple syntax change.