In clojure, what's the efficient way to convolve a vector by a kernel?

Question

I came up with this:

(def kernel [0 1 1 2 3 3 0 0 0 0 0 0])
(def data [1 5 7 4 8 3 9 5 6 3 2 1 1 7 4 9 3 2 1 8 6 4])

(defn capped+ [a b c] (let [s (+ a b)] (if (> s c) c s)))

(defn *+ [a b]
  (if (> (count a) (count b))  
    (reduce + (map-indexed (fn _ [i x] (* (a i) (b i))) b))
    (reduce + (map-indexed (fn _ [i x] (* (a i) (b i))) a))))

(defn slide*i [d k] 
 (let [ki (into [] (reverse k)) kl (count k) dl (count d)]
   (map-indexed 
      (fn [idx item] 
        (/ (*+ ki (subvec d idx (capped+ idx kl dl))) 
           (reduce + ki))) 
      d)))

(def s (slide*i data kernel))

It's not the most elegant code, but it works fine. I actually want to use it to smooth some very spiky! data.

Any suggestions to make this more beautiful or more efficient or more accurate (personally I don't care about the tail being inaccurate because in my case I never use it) are welcomed.

Answer 1

You can certainly improve the performance of this operation significantly. The good news is that you don't need to drop into Java for this: Clojure is extremely fast if you optimise it correctly and in most instances can produce the same speed as pure Java.

For maximum performance of numerical code in Clojure you will want to use:

• arrays, because you want mutable storage with very fast writes and lookup. Clojure sequences and vectors are beautiful, but they come with overheads that you probably want to avoid for truly performance-critical code
• double primitives, because they offer much faster maths.
• aset / aget / areduce - these are extremely fast operations designed for arrays and basically give you the same bytecode as pure Java equivalents.
• imperative style - although it's unidiomatic in Clojure, it gets the fastest results (mainly because you can avoid overheads from memory allocations, boxing and function calls). An example would be using dotimes for a fast imperative loop.
• (set! *warn-on-reflection* true) - and eliminate any warnings your code produces, because reflection is a big performance killer.

The following should be along the right lines and will probably get you roughly equivalent performance to Java:

(def kernel (double-array [0 1 1 2 3 3 0 0 0 0 0 0]))
(def data (double-array [1 5 7 4 8 3 9 5 6 3 2 1 1 7 4 9 3 2 1 8 6 4]))

(defn convolve [^doubles kernel-array ^doubles data-array]
  (let [ks (count kernel-array)
        ds (count data-array)
        output (double-array (+ ks ds))
        factor (/ 1.0 (areduce kernel-array i ret 0.0 (+ ret (aget kernel-array i))))]    
    (dotimes [i (int ds)]
      (dotimes [j (int ks)]
        (let [offset (int (+ i j))]
          (aset output offset (+ (aget output offset) (* factor (* (aget data-array i) (aget kernel-array j))))))))
    output))

(seq (convolve kernel data))
=> (0.0 0.1 0.6 1.4 2.4 4.4 5.5 6.1000000000000005 5.600000000000001 6.200000000000001 5.499999999999999 5.9 4.199999999999999 3.3000000000000003 2.5 2.2 3.3 4.4 5.6000000000000005 4.8 4.8999999999999995 3.1 3.5 4.300000000000001 5.0 3.0 1.2000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0)

I've not trimmed the output array or done any bounding so you'll probably need to hack this solution a bit to get exactly the output you want, but hopefully you get the idea.....

Some very rough benchmarking:

(time (dotimes [i 1000] (seq (convolve kernel data))))
=> "Elapsed time: 8.174109 msecs"

i.e. that's about 30ns per kernel / data pair combination - I expect that's pretty much hitting the bounds of cached memory access.