Sparse convolution of sparse arrays

Question

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.)

Answer 1

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[10000]; // 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