Isaac Hodes posted a request to both StackOverflow and the Clojure mailing list, asking for some help implementing a discrete convolution function in Clojure.

I am trying to create a lazy/functional/efficient/Clojuresque function to carry out convolution on two lists/vectors of (ideally BigDecimals but that may be inefficient) doubles.

Right now, it seems as though this is something Clojure cannot do.

Let’s cut right to the code. Here’s my implementation:

So, how does this implementation work? And just what does the convolve function compute? Let’s consider the latter question first.

Recall how you multiply polynomials. Suppose we want to multiply 1+x by 2+3x+5x^2. Using the standard rules of algebra, we find the answer by gathering like terms. The order-0 term is 1*2=2. The order-1 term is 1*3x + x*2=5x, the order-2 term is x*3x + 5x^2=8x^2, and the order-3 term is x*5x^2=5x^3. It is notationally convenient to drop the variables and represent the polynomials as vectors of coefficients in standard form. In the case above, we multiply [1 1] by [2 3 5] and get [2 5 8 5]. The convolve function performs exactly the computation required to multiply the polynomials.

So, how does it work? My hope is that some (crude) pictures will tell the story; you be the judge. Let’s consider a particular example, (convolve [1 2 3] [4 5 6 7 ])
whose value is (4 13 28 34 32 21). Here’s the picture:

3---2---1
--------4---5---6---7    ;; 1*4=4

----3---2---1
--------4---5---6---7    ;; 2*4+1*5=13

--------3---2---1
--------4---5---6---7    ;; 3*4+2*5+1*6=28

------------3---2---1    
--------4---5---6---7    ;; 3*5+2*6+1*7=34

----------------3---2---1
--------4---5---6---7    ;; 3*6+2*7=32

--------------------3---2---1
--------4---5---6---7    ;; 3*7=21

As requested, the code is lazy:

user> (class (convolve [1 2 3] [4 5 6 7 ]))
clojure.lang.LazySeq

But does it perform? Let’s take a look.

user> (def foo2 (into [] (range 100)))
#'user/foo2

user> (def foo3 (into [] (range 1000)))
#'user/foo3

user> (def foo4 (into [] (range 10000)))
#'user/foo4

user> (time (do (doall (convolve [1 2 3] foo2 )) nil))
"Elapsed time: 2.827 msecs"
nil

user> (time (do (doall (convolve [1 2 3] foo3 )) nil))
"Elapsed time: 76.275023 msecs"
nil

user> (time (do (doall (convolve [1 2 3] foo4 )) nil))
"Elapsed time: 4101.002754 msecs"

Ouch. The time appears to be quadratic in the length of the second argument. We can do better than that with just a minor tweak. The problem is the performance of the drop function. Let’s get rid of it. Here’s the new version of the code.

And here are the new performance numbers, which look much better. A million points in under three seconds, with no dropping into java-isms. Sweet.

user> (def foo5 (into [] (range 100000)))
#'user/foo5

user> (def foo6 (into [] (range 1000000)))
#'user/foo6

user> (time (do (doall (convolve [1 2 3] foo3 )) nil))
"Elapsed time: 6.909401 msecs"
nil

user> (time (do (doall (convolve [1 2 3] foo4 )) nil))
"Elapsed time: 62.894431 msecs"
nil

user> (time (do (doall (convolve [1 2 3] foo5 )) nil))
"Elapsed time: 255.064899 msecs"
nil

user> (time (do (doall (convolve [1 2 3] foo6) nil)))
"Elapsed time: 2784.131159 msecs"
nil