Sparse convolution of sparse arrays

The documentation for ListConvolve mentions that "ListConvolve works with sparse arrays", which is true. The result, however is never sparse, eg:

a = SparseArray[50 -> x, 1000];
b = SparseArray[40 -> 2, 50];
ListConvolve[b, a]

Is there a way to do "sparse convolution"?

(The background is multiplication of high-degree many variables polynomials.)

• ListConvolve uses Fourier transforms, and it's very difficult to achieve significantly better performance with the FFT for general sparse inputs and outputs than one can get by doing the full transform. So, while the answer to your question is certainly "yes", producing a solution is certainly more than a matter of a simple programming problem. – Oleksandr R. Nov 20 '12 at 13:44
• If these polynomials are very sparse you can perhaps do it as a sparse matrix-vector multiply, which of course will give you a sparse result. But the performance achievable in this fashion might not be very impressive. – Oleksandr R. Nov 20 '12 at 13:46
• isn't sparse multiplication with a circular matrix a convolution? – chris Nov 20 '12 at 13:51
• @chris, yup, but whether that route is more efficient than FFT depends on quite a lot of things... – J. M. will be back soon Nov 20 '12 at 13:56
• @J.M. sure: how Sparse is the convolution to start with. – chris Nov 20 '12 at 14:04

Sparse multiplication with a circular matrix corresponds to a convolution; on a trivial example let us compare:

matrix = SparseArray[{Band[{1, 1}] -> 2, Band[{1, 2}] -> 1, Band[{2, 1}] -> 1}, {15, 15}];

vec = SparseArray[5 -> x, 15]; matrix.vec // Normal

(* ==> {0, 0, 0, x, 2 x, x, 0, 0, 0, 0, 0, 0, 0, 0, 0} *)

versus

a = SparseArray[5 -> x, 15] // Normal; b = {1, 2, 1};
ListConvolve[b, a]

(* ==> {0, 0, x, 2 x, x, 0, 0, 0, 0, 0, 0, 0, 0} *)

So it seems ListConvolve treats edges differently than circular matrix multiplication. Must be described in the documentation of ListConvolve.

In terms of performance, the issue is how diagonal your matrix is or equivalently how long a or b is.

EDIT

Let's do some timing to estimate the performance of convolution as a function of how sparse the matrix is. Consider a band matrix with a custom width of 2p

Clear[matrix];
matrix[p_, n_] := SparseArray[Join[Table[Band[{1, j}] -> 1/p, {j, 1, p}],
Table[Band[{j, 1}] -> 1/p, {j, 2, p}]], {n, n}];

Clear[vec]; vec[n_] := SparseArray[5 -> x, n];

Let us see how the time to carry out the multiplication increases with the width

Table[{p,matrix[p, 10000].vec; // Timing // First},{p, 2, 10}]

(*
==> ({
{2, 0.049756},
{3, 0.094646},
{4, 0.141837},
{5, 0.197723},
{6, 0.244743},
{7, 0.287893},
{8, 0.344092},
{9, 0.403023},
{10, 0.472343}
})
*)

The scaling with the width p seems to be $p^{4/3}$ (this includes the time to build the convolution matrix)

Fit[Log10[%], {1, x}, x]

(* ==> 1.31 x -1.65 *)

Of course scaling could be different if say the Arrays were made of finite precision numbers.

EDIT 2

Note that the method works in arbitrary dimensions; let us define a 2D convolution matrix using splines:

design = Table[BSplineBasis[3, (x - xi + 1/2)] BSplineBasis[
3, (y - yi + 1/2)], {xi, 1, 5, 1/4}, {yi, 1, 5,1/4}, {x, 1, 5, 1/4}, {y, 1, 5, 1/4}];

matrix = Partition[design// Flatten, 17^2] // SparseArray;

(here I am just lazy; one would of course need to fill the matrix directly as a sparse matrix).

Let us now define a 2D field

field = Table[If[x == 3 && y == 3, 1, 0], {x, 1, 5, 1/4}, {y, 1, 5, 1/4}];
field // MatrixPlot Lets convolve this 2D field by the above defined sparsified matrix using sparse multiplication:

Partition[matrix.SparseArray[Flatten[field]], 17] // MatrixPlot 