Trouble with lazy convolution fn in Clojure

Question

I am writing some signal processing software, and I am starting off by writing out a discrete convolution function.

This works fine for the first ten thousand or so list of values, but as they get larger (say, 100k), I begin to get StackOverflow errors, of course.

Unfortunately, I am having a lot of trouble converting the imperative convolution algorithm I have to a recursive & lazy version that is actually fast enough to use (having at least a modicum of elegance would be nice as well).

I am also not 100% sure I have this function completely right, yet – please let me know if I'm missing something/doing something wrong. I think it's correct.

(defn convolve
  "
    Convolves xs with is.

    This is a discrete convolution.

    'xs  :: list of numbers
    'is  :: list of numbers
  "
  [xs is]
  (loop [xs xs finalacc () acc ()]
    (if (empty? xs)
      (concat finalacc acc)
      (recur (rest xs)
             (if (empty? acc)
               ()
               (concat finalacc [(first acc)]))
             (if (empty? acc)
               (map #(* (first xs) %) is)
               (vec-add
                (map #(* (first xs) %) is)
                (rest acc)))))))

I'd be much obliged for any sort of help: I'm still getting my bearings in Clojure, and making this elegant and lazy and/or recursive would be wonderful.

I'm a little surprised how difficult it is to express an algorithm which is quite easy to express in an imperative language in a Lisp. But perhaps I'm doing it wrong!

EDIT: Just to show how easy it is to express in an imperative language, and to give people the algorithm that works nicely and is easy to read, here is the Python version. Aside from being shorter, more concise and far easier to reason about, it executes orders of magnitude faster than the Clojure code: even my imperative Clojure code using Java arrays.

from itertools import repeat

def convolve(ns, ms):
    y = [i for i in repeat(0, len(ns)+len(ms)-1)]
    for n in range(len(ns)):
        for m in range(len(ms)):
            y[n+m] = y[n+m] + ns[n]*ms[m]
    return y

Here, on the other hand, is the imperative Clojure code. It also drops the last, non fully-immersed, values from the convolution. So aside from being slow and ugly, it's not 100% functional. Nor functional.

(defn imp-convolve-1
  [xs is]
  (let [ys (into-array Double (repeat (dec (+ (count xs) (count is))) 0.0))
        xs (vec xs)
        is (vec is)]
     (map #(first %)
          (for [i (range (count xs))]
            (for [j (range (count is))]
              (aset ys (+ i j)
                    (+ (* (nth xs i) (nth is j))
                       (nth ys (+ i j)))))))))

This is so disheartening. Please, someone show me I've just missed something obvious.

EDIT 3:

Here's another version I thought up yesterday, showing how I'd like to be able express it (though other solutions are quite elegant; I'm just putting another one out there!)

(defn convolve-2
  [xs is]
  (reduce #(vec-add %1 (pad-l %2 (inc (count %1))))
         (for [x xs]
           (for [i is]
             (* x i)))))

It uses this utility function vec-add:

(defn vec-add
  ([xs] (vec-add xs []))
  ([xs ys]
     (let [lxs (count xs)
           lys (count ys)
           xs (pad-r xs lys)
           ys (pad-r ys lxs)]
       (vec (map #(+ %1 %2) xs ys))))
  ([xs ys & more]
     (vec (reduce vec-add (vec-add xs ys) more))))
     (vec (reduce vec-add (vec-add xs ys) more))))

Answer 1

(defn ^{:static true} convolve ^doubles [^doubles xs ^doubles is]
  (let [xlen (count xs)
        ilen (count is)
        ys   (double-array (dec (+ xlen ilen)))]
    (dotimes [p xlen]
      (dotimes [q ilen]
        (let [n (+ p q), x (aget xs p), i (aget is q), y (aget ys n)]
          (aset ys n (+ (* x i) y)))))
    ys))

Riffing on j-g-faustus's version if I was doing this in the Clojure equiv branch. Works for me. ~400ms for 1,000,000 points, ~25ms for 100,000 on a i7 Mackbook Pro.

Answer 2

The likely cause of the stack overflow errors is that the lazy thunks are getting too deep. (concat and map are lazy). Try wrapping those calls in doall to force evaluation of their return values.

As for a more functional solution, try something like this:

(defn circular-convolve
  "Perform a circular convolution of vectors f and g"
  [f g]
  (letfn [(point-mul [m n]
        (* (f m) (g (mod (- n m) (count g)))))
      (value-at [n]
        (reduce + (map #(point-mul % n) (range (count g)))))]
    (map value-at (range (count g)))))

Use can use reduce to perform summation easily, and since map produces a lazy sequence, this function is also lazy.

Answer 3

Can't help with a high-performance functional version, but you can get a 100-fold speedup for the imperative version by foregoing laziness and adding type hints:

(defn imp-convolve-2 [xs is]
  (let [^doubles xs (into-array Double/TYPE xs)
        ^doubles is (into-array Double/TYPE is)
        ys (double-array (dec (+ (count xs) (count is)))) ]
    (dotimes [i (alength xs)]
      (dotimes [j (alength is)]
        (aset ys (+ i j)
          (+ (* (aget xs i) (aget is j))
            (aget ys (+ i j))))))
    ys))

With xs size 100k and is size 2, your imp-convolve-1 takes ~6,000ms on my machine when wrapped in a doall, while this one takes ~35ms.

Update

Here is a lazy functional version:

(defn convolve 
  ([xs is] (convolve xs is []))
  ([xs is parts]
    (cond
      (and (empty? xs) (empty? parts)) nil
      (empty? xs) (cons
                    (reduce + (map first parts))
                    (convolve xs is
                      (remove empty? (map rest parts))))
      :else (cons
              (+ (* (first xs) (first is))
                (reduce + (map first parts)))
              (lazy-seq
                (convolve (rest xs) is
                  (cons 
                    (map (partial * (first xs)) (rest is))
                    (remove empty? (map rest parts)))))))))

On sizes 100k and 2, it clocks in at ~600ms (varying 450-750ms) vs ~6,000ms for imp-convolve-1 and ~35ms for imp-convolve-2.

So it's functional, lazy and has tolerable performance. Still, it's twice as much code as the imperative version and took me 1-2 additional hours to find, so I'm not sure that I see the point.

I'm all for pure functions when they make the code shorter or simpler, or have some other benefit over an imperative version. When they don't, I have no objection to switch to imperative mode.

Which is one of the reasons I think Clojure is great, since you can use either approach as you see fit.

Update 2:

I'll amend my "what's the point of doing this functionally" by saying that I like this functional implementation (the second one, further down the page) by David Cabana.

It's brief, readable and times to ~140ms with the same input sizes as above (100,000 and 2), making it by far the best-performing functional implementation of those I tried.

Considering that it is functional (but not lazy), uses no type hints and works for all numeric types (not just doubles), that's quite impressive.