A functional, efficient, convolution function. Is imperitive-style the answer?

[Isaac Hodes]

I posted this on StackOverflow yesterday, but to no avail: I don't
think many people looked at it, or least I didn't get much feedback.

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.

It is turning out to be very difficult. I have four or five versions
in my buffer right now, and none of them are acceptable. I've even
tried a number of versions using into-array etc, e.g. mutables.

Instead of posting a lot of links to pastie, I'm just copying the SO
link here, where my code and the algorithms can be found.

Right now, it seems as though this is something Clojure cannot do.
David Nolen mentioned that the appropriate Java-interop functionality
may come in Clojure 1.3, but aside from eschewing the functional style
by using arrays, it seems as though there should be a better answer.
(As an aside, where can I find out about 1.3?)

I have a Python implementation to show how nicely it can be expressed
in an imperitive language (though I suppose I dabbled in the
functional style initializing my return array).

Is it possible to do this, and do it well?

Thank you so much, in advance.

[Frederick Polgardy]

Where's the link? :)

-Fred

--
Science answers questions; philosophy questions answers.

> --
> You received this message because you are subscribed to the Google
> Groups "Clojure" group.
> To post to this group, send email to clo...@googlegroups.com
> Note that posts from new members are moderated - please be patient with your first post.
> To unsubscribe from this group, send email to
> clojure+u...@googlegroups.com
> For more options, visit this group at
> http://groups.google.com/group/clojure?hl=en

[Michael Gardner]

I'm assuming the StackOverflow link you refer to is http://stackoverflow.com/questions/3259825/trouble-with-lazy-convolution-fn-in-clojure.

I would think about the problem this way: to compute the value at index i in the output list, you multiply together each pair of values in the input lists whose indices add up to i, then add them all together. So I would first write a small helper function to return a list of such index pairs for a given output index, then use it with map to compute each output value in one go.

(defn convolve-indices [i max-m max-n]
"Lists all index pairs adding up to i, where the first index is less than max-m and the second is less than max-n."
(filter #(< (first %) max-m)
(filter #(< (second %) max-n)
(map #(vector % (- i %))
(range (inc i))))))

(defn convolve-1 [i ms ns]
"Returns the ith value of the convolution of ms with ns."
(reduce +
(map #(* (ms (first %)) (ns (second %)))
(convolve-indices i (count ms) (count ns)))))

(defn convolve [ms ns]
"Convolves ms with ns."
(map #(convolve-1 % ms ns)
(range (dec (+ (count ms) (count ns))))))

I would expect this to be much slower than the imperative approach, due to all the temporary objects and the lack of type hints. There should be a much less wasteful way to write convolve-indices, for instance. On the plus side, you may be able to use pmap to parallelize the work with little effort.

[James Reeves]

On 16 July 2010 20:57, Isaac Hodes <iho...@mac.com> wrote:
> 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.

Perhaps something like this?

(defn convolve [ns ms]
(loop [nums (for [[i n] (indexed ns)
[j m] (indexed ms)] [(+ i j) (+ n m)])
y (transient (vec (repeat (+ (count ns) (count ms)) 0)))]
(if-let [[i s] (first nums)]
(recur (next nums) (assoc! y i (+ (y i) s)))
(persistent! y))))

- James

[Isaac Hodes]

On Jul 16, 11:26 pm, Frederick Polgardy <f...@polgardy.com> wrote:
> Where's the link? :)
>
> -Fred

Whew. Here it is; my bad!

http://stackoverflow.com/questions/3259825/trouble-with-lazy-convolution-fn-in-clojure

[Frederick Polgardy]

This statement is interesting. Part of what makes the Python implementation so expressive is precisely the use of functional constructs like list comprehensions, ranges, and all the awesome stuff in itertools. A purely imperative implementation would probably be as clunky as the purely functional ones, but for different reasons.

-Fred

--
Science answers questions; philosophy questions answers.

On Jul 16, 2010, at 2:57 PM, Isaac Hodes wrote:

[David Cabana]

I tried to run Jame's code, but the compiler (1.2 beta) squawked at me:
Unable to resolve symbol: indexed in this context
[Thrown class java.lang.Exception]

What do I need to pull in to pick up the "indexed" function?

[Wilson MacGyver]

It's in contrib.seq-utils

http://richhickey.github.com/clojure-contrib/seq-utils-api.html#clojure.contrib.seq-utils/indexed

[James Reeves]

On 17 July 2010 18:41, David Cabana <drca...@gmail.com> wrote:
> I tried to run Jame's code, but the compiler (1.2 beta) squawked at me:
> Unable to resolve symbol: indexed in this context
>  [Thrown class java.lang.Exception]
>
> What do I need to pull in to pick up the "indexed" function?

Sorry, I should have mentioned you need Clojure 1.2.0 beta,
Clojure-Contrib 1.2.0 beta, and you need to include "indexed" from the
clojure.contrib.seq namespace:

(use '[clojure.contrib.seq :only (indexed)])

For version 1.1.0, it needs a little adjusting. The following code
should work under 1.1.0 (but I haven't tested this):

(use '[clojure.contrib.seq-utils :only (indexed)])

(defn convolve [ns ms]
(loop [nums (for [[i n] (indexed ns)
[j m] (indexed ms)] [(+ i j) (+ n m)])

y (vec (repeat (+ (count ns) (count ms)) 0))]

(if-let [[i s] (first nums)]

(recur (next nums) (assoc y i (+ (y i) s)))
y)))

This code will be slower, however, as it isn't using transients.
Transients are a way of speeding up manipulating data structures by
assuming you don't need to keep around the old copies.

- James

[David Cabana]

I thought this problem was interesting enough to merit better
treatment than I can give it here, hence a blog post was in order.
Brief version: I think I have a lazy, functional, and idiomatic
implementation with decent performance. Check it out here:

http://erl.nfshost.com/2010/07/17/discrete-convolution-of-finite-vectors/

[Carson]

Hi, Try this:

(defn stationary-convolve [n fs gs]
(let [window-size (count fs)
padding (repeat window-size 0)
padded-gs (concat padding gs padding)
new-center (+ n window-size)
window-of-gs (subvec (vec padded-gs)
(- (inc new-center) window-size)
(inc new-center))
reverse-window-fn (fn [i] (nth window-of-gs
(mod (- -1 i) window-size)))
reverse-window (map reverse-window-fn (range 0 window-size))]
(reduce + (map * fs reverse-window))))

(defn convolve [fs gs]
(map #(stationary-convolve % fs gs)
(range 0 (+ (count fs) (count gs)))))

(let [xs (for [i (range 1000000)] (rand 100))
ys (for [i (range 1000000)] (rand 100))]
(time (dotimes [_ 3] (convolve xs ys))))

"Elapsed time: 2540.967482 msecs"

Seems to run pretty fast and seems to be the kind of convolution you
want.

Thanks for the interesting problem!

Carson

[Carson]

I should add, "oops, ignore what I wrote". :)

see:

https://groups.google.com/group/clojure/tree/browse_frm/thread/33366bccc6df7756/415072576d83b757?rnum=11&_done=%2Fgroup%2Fclojure%2Fbrowse_frm%2Fthread%2F33366bccc6df7756%3F#doc_0d0a3759a0a10328

Carson

[Mark Engelberg]

I have found Clojure to be consistently faster than Python, so I
thought it would be instructive to compare the Python code to the
closest Clojure equivalent.

Here's the originally posted Python code:
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

Now, if you want to write this functionally in Clojure, you've got two
main obstacles in your way. First, you can't "bang a vector in place"
-- you want to use Clojure's built-in persistent vectors. For the
best speed, you can call transient at the beginning of the operations
and persistent! at the end, but you still need to pass the vector
around as an accumulator. This leads to the second obstacle --
because you're not "banging a vector in place" you can't really use
for loops or doseq loops to conveniently iterate through the indices.
You need to unroll these for loops into loop/recur constructs. It's a
bit of a nuisance to do all this, but it's actually a fairly
straightforward translation process. Here's what you get:

(defn convolve [ns ms]
(let [ns (vec ns), ms (vec ms),
count-ns (count ns), count-ms (count ms)]
(loop [n 0,
y (transient (vec (repeat (dec (+ count-ns count-ms)) 0)))]
(if (= n count-ns) (persistent! y)
(recur (inc n)
(loop [m 0, y y]
(if (= m count-ms) y
(let [n+m (+ n m)]
(recur (inc m) (assoc! y n+m
(+ (y n+m) (* (ns n) (ms m)))))))))))))

By my measurements, this Clojure code runs (under java -server) twice
as fast as the comparable Python code. Considering that Clojure's
math is inherently slower than Python's, and Clojure is using more
expensive persistent data structures, it's impressive that Clojure is
faster. Clojure also offers the opportunity to go faster (with Java
arrays), and even to deliver the results lazily if that's useful to
you (but unsurprisingly, the laziness comes at a bit of a cost).
Others have described those solutions, so I won't rehash those here.

But there are two strikes against Clojure here:
1. This Clojure code is definitely harder to read than the Python
equivalent. The transformation obfuscates what's going on.
2. Running the Python code under the psyco JIT compiler speeds the
Python code up by a factor of 10, making it about 5 times faster than
this Clojure transformation. So Python still gets the final word in
speed here (comparing the exact same algorithm on their built-in data
structures).

[gary ng]

J(a language of the APL family) excels at this problem domain. It is functional and very efficient.

 

conv =: +//.@*/

 

It looks like line noise but that is the whole implementation

 

(1 1 ) conv (1 1)

=> 1 2 1

 

http://www.jsoftware.com/