In clojure, what is an efficient way to convolve a vector with 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+ [abc] (let [s (+ ab)] (if (> sc) cs))) (defn *+ [ab] (if (> (count a) (count b)) (reduce + (map-indexed (fn _ [ix] (* (ai) (bi))) b)) (reduce + (map-indexed (fn _ [ix] (* (ai) (bi))) a)))) (defn slide*i [dk] (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)) 

This is not the most elegant code, but it works great. I actually want to use it to smooth out the very sharp! data.

Any suggestions to make it more beautiful or more effective or more accurate (it’s not important for me personally that the tail is inaccurate, because in my case I never use it) are welcome.

Answer 1

You can significantly improve the performance of this operation. The good news is that you don’t have to rush into Java for this: Clojure is very fast if you optimize it correctly and in most cases you can create the same speed as pure Java.

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

• arrays because you want to change the repository with very fast write and search. Clojure sequences and vectors are beautiful, but they come with overheads that you probably want to avoid for really critical code
• double because they offer much faster mathematical models.
• aset / aget / areduce are extremely fast operations designed for arrays and basically give you the same bytecode as pure Java equivalents.
• imperative style - although it is uniform in Clojure, it gets the fastest results (mainly because you can avoid the overhead of memory allocation, boxing and function calls). An example is dotimes for a fast imperative cycle.
• (set! * warn-on-reflection * true) - and eliminate any warnings that your code generates, because reflection is a big killer of performance.

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

 (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 (+ ij))] (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 did not trim the output array or make any restrictions, so you probably have to crack this solution a bit to get exactly the result 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 about 30 ns per core / data combination - I expect that I will be able to largely limit access to cached memory.

Answer 2

  ;  unused, but left for demonstration
 (defn capped + [abc]
   (min (+ ab) c))

 (defn * + [ab]
   (reduce + (map * ab)))

 (defn slide * i [dk]
   (let [ki (reverse k)
         kl (count k)
         ks (reduce + k)]
     (map # (/ (* +% 1 ki) ks) (partition kl 1 d))))

With partition results:

  (59/10 21/5 33/10 5/2 11/5 33/10 22/5 28/5 24/5 49/10 31/10) 

But with partition-all you get exactly what your solution led to:

  (59/10 21/5 33/10 5/2 11/5 33/10 22/5 28/5 24/5 49/10 31/10 7/2 43/10 5 3 6/5 0 0 0 0 0 0) 

Answer 3

An effective way to do this is to create a java class that does the convolution and call it from clojure, passing it a java array, if possible. Clojure implementation should work on java arrays if efficiency is a problem.