Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
5  
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
1  
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. Nov 20 '12 at 13:56
    
@J.M. sure: how Sparse is the convolution to start with. –  chris Nov 20 '12 at 14:04
show 2 more comments

1 Answer

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

Mathematica graphics

Lets convolve this 2D field by the above defined sparsified matrix using sparse multiplication:

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

Mathematica graphics

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.