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 imperitive convolution algorithm I have to a recursive & lazy version that is actually fast enough to use, and 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 both 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 imperitive language in a Lisp. But perhaps I'm doing it wrong!

EDIT: Just to show how easy it is to express in an imperitive 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 imperitive 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 imperitive 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 solution 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

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 2

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.

Answer 3

(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 4

(defn convolve [xs is]
  (if (> (count xs) (count is))
    (convolve is xs)
    (let [os (dec (+ (count xs) (count is)))
          lxs (count xs)
          lis (count is)]
      (for [i (range os)]
        (let [[start-x end-x] [(- lxs (min lxs (- os i))) (min i (dec lxs))]
              [start-i end-i] [(- lis (min lis (- os i))) (min i (dec lis))]]
          (reduce + (map *
                         (rseq (subvec xs start-x (inc end-x)))
                         (subvec is start-i (inc end-i)))))))))

It is possible to express a lazy, functional solution in concise terms. Alas, the performance for > 2k is impractical. I'm interested to see if there are ways to speed it up without sacrificing readability.

Edit: After reading drcabana's informative post on the topic (http://erl.nfshost.com/2010/07/17/discrete-convolution-of-finite-vectors/), I've updated my code to accept different sized vectors. His implementation is better performing: for xs size 3, is size 1000000, ~2s vs ~3s

Edit 2: Taking drcabana's ideas of simply reversing xs and padding is, I arrived at:

(defn convolve [xs is]
  (if (> (count xs) (count is))
    (convolve is xs)
    (let [is (concat (repeat (dec (count xs)) 0) is)]
      (for [s (take-while not-empty (iterate rest is))]
         (reduce + (map * (rseq xs) s))))))

This is probably as concise as it's going to get, but it is still slower overall, likely due to take-while. Kudos to the blog author for a well considered approach. The only advantage here is that the above is truly lazy in that if I ask (nth res 10000), it will only need the first 10k calculations to arrive at a result.