Response

28 views
Skip to first unread message

Isaac Hodes

unread,
Jul 17, 2010, 4:43:33 PM7/17/10
to Clojure
Thanks everyone for the input! I've also got some nice replies on
Stack Overflow – there's some sort of a discussion there as well.

I thought I'd add my voice to the crowd. First off, purely imperative
C is fairly elegant itself. I've made it a simple function (along with
the Python version on the SO post), though I could easily use malloc
or some other way to make ys[] dynamic and not have it required as a
function argument, and I could have zeroed it in the function as well,
but I left that off for simplicity's sake.

http://pastie.org/1048769 for formatted version of:

double convolve(double *xs, double *is, double *ys){
int i,j;
for(i=0; i<len(xs); i++){
for(j=0; j<len(is); j++){
ys[i+j] = ys[i+j] + (xs[i]*is[j]);
}
}
return ys;
}

As for another Clojure version, this one is short and dare I say a
little bit elegant, but uses a couple utility functions I have around
in my utils.clj file for thie project anyway, and it is impractical
and slow for even 1000 points, though as you can see it's a naive
implementation. http://paste.lisp.org/display/112570#1

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

(defn vec-add
"Adds vectors
together.

e.g. (vec-add [1 2 3] [6 4 4]) => [7 6
7]
(vec-add [1 2] [5 5 5 5]) => [6 7 5
5]
"
([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))))

(defn pad-l
"Pads 'xs on the left with 'e, to length 'len, if it's less than
'len."
([xs len] (pad-l xs len 0))
([xs len e]
(if (<= len (count xs))
xs
(recur (cons e xs) len e))))

These and other answers that are more lispy don't quite cut it, but as
someone mentioned above me, I could use pmap and parallelize the
procedure some.

I think my question has been answered: there are ways to do it in
Clojure, but to get efficiency you have to drop down into the almost-
Java level.

It's just a shame, it seems to me, that there is such a nice way to
represent the procedure in Python or even C, yet Clojure (or any Lisp
really) struggles to idiomatically answer this question of
convolution.

Interesting, I think, is the fact that Haskell can more idiomatically
handle convolution; though I can't quite read the examples that are
around the web (I need to get on that…). I wonder if there's something
that could be added to Clojure to do what Haskell does, if it's really
doing anything differently.

Philosophically, it's interesting that a problem like convolution, in
which a signal is being *mutated*, is difficult to express in an
functional manner without some struggle or awkwardness. It is a rare
case where something is actually being changed, mutated if you will,
but not over time so much as *at* the time you look at it. I'm not
sure I'm expressing myself well here, but there seems to be a
distinction between the problem of convolution and that other other,
more conventional and functional problems.

Either way, thank you all again.



David Nolen

unread,
Jul 17, 2010, 7:04:28 PM7/17/10
to clo...@googlegroups.com
(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))

I don't think this is so bad and it can do a million points in ~400ms on my machine. 100,000 points in ~25ms. I don't consider this dropping to the Java level at all.

Michael Gardner

unread,
Jul 17, 2010, 7:11:19 PM7/17/10
to clo...@googlegroups.com
On Jul 17, 2010, at 3:43 PM, Isaac Hodes wrote:

> It's just a shame, it seems to me, that there is such a nice way to
> represent the procedure in Python or even C, yet Clojure (or any Lisp
> really) struggles to idiomatically answer this question of
> convolution.

No, it's pretty easy to do convolution idiomatically, as demonstrated by some of the replies you got. What makes it seem difficult is that the literature is all written with imperative languages in mind.

What Clojure does "struggle" with is performance, in that the obvious way of writing something often results in poor performance, especially for heavy numeric code. This applies to most (all?) dynamic languages, though; you'd be better off with C if high-performance numerics are your top priority.

Carson

unread,
Jul 17, 2010, 7:17:25 PM7/17/10
to Clojure
Hi David,
Would appreciate if you could try out my attempt at this [1] on your
machine. My machine ain't as fast as yours apparently...

[1]
https://groups.google.com/group/clojure/tree/browse_frm/thread/ee4169bc292ab572/6d03461efde166ad?rnum=11&_done=%2Fgroup%2Fclojure%2Fbrowse_frm%2Fthread%2Fee4169bc292ab572%2Ff769ca6eb6fb8622%3Ftvc%3D1%26#doc_6d03461efde166ad

Thanks,
Carson

David Nolen

unread,
Jul 17, 2010, 7:58:06 PM7/17/10
to clo...@googlegroups.com
I tried this, it runs in about the same amount of time as it did for you.

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

j-g-faustus

unread,
Jul 17, 2010, 8:33:57 PM7/17/10
to Clojure
On 17 Jul, 22:43, Isaac Hodes <iho...@mac.com> wrote:
> It's just a shame, it seems to me, that there is such a nice way to
> represent the procedure in Python or even C, yet Clojure (or any Lisp
> really) struggles to idiomatically answer this question of
> convolution.

I'll have to disagree with the conclusion. I agree that there are many
algorithms that can be elegantly expressed as iteration over mutable
arrays, but in any Lisp with arrays (including at least Clojure and
Common Lisp) you can use those algorithms in Lisp just as much as you
can in Python or C. I can't really see that it is less idiomatic to
use arrays in Lisp than in C when it suits the problem.

The only difference I can tell between David Nolen's version and the
C, Python or DSP Guide algorithm is that Clojure has a more verbose
syntax for accessing arrays - the logic, flow and operations are
exactly the same. (And you don't need to initialize the result array
to 0.)

Philosophically, I think it's mostly a matter of the set of 'cheap' vs
'expensive' operations being different in imperative and functional
code, so for a well-performing functional version you normally need a
different algorithm, which can be a non-trivial task to discover. But
you did get a couple of nice looking and decently performing
functional versions among the responses.

Practically, I go for imperative algorithms most of the time since I
can copy them pretty much verbatim from the literature and don't need
to spend hours on each of them trying to find an efficient functional
formulation. Plus of course that the imperative versions are faster.
I'm writing a fair amount of this kind of code, and I'm using Clojure
since the number crunching can be an order of magnitude faster than
Python or Ruby (for my purposes Python is simply too slow, while
imperative and type-hinted Clojure is fast enough), and Clojure is
much nicer to use than Java or C.

But certainly, if you need numeric performance over everything C or
Fortran would be the language of choice.


jf

Frederick Polgardy

unread,
Jul 17, 2010, 8:44:09 PM7/17/10
to clo...@googlegroups.com
This example is beside the point of the original question. It uses mutable arrays. It's very much dropping to the Java level. Am I missing something?

-Fred

--
Science answers questions; philosophy questions answers.

Carson

unread,
Jul 17, 2010, 9:07:17 PM7/17/10
to Clojure
Have you taken a look at my attempt at a solution (on the other/
original thread)?

https://groups.google.com/group/clojure/tree/browse_frm/thread/ee4169bc292ab572/6d03461efde166ad?rnum=11&_done=%2Fgroup%2Fclojure%2Fbrowse_frm%2Fthread%2Fee4169bc292ab572%3F#doc_6d03461efde166ad

I don't know if it can be called idiomatic, but it looks functional to
me, doesn't use any type hints or mutable arrays, and it seems pretty
fast. Convolves two 1000000 size vectors in roughly ~600 ms I
think...

Carson

Frederick Polgardy

unread,
Jul 17, 2010, 9:13:01 PM7/17/10
to clo...@googlegroups.com
I think it really doesn't get any clearer than this in terms of intent. While I was adept at calculus-level math 20 years ago, I've forgotten the little I knew of matrices. This is the first algorithm that has communicated by visual inspection (to me) exactly what a convolution is.

-Fred

--
Science answers questions; philosophy questions answers.

Carson

unread,
Jul 18, 2010, 12:32:58 AM7/18/10
to Clojure
I guess different people see things differently. I actually didn't
understand the intent of the imperative version, even though it takes
less lines I guess. There's quite a few things called convolution.
My naive implementation was functional and pretty fast. Turns out it
can be even faster with just a little work (still functional):

(defn point-convolve [n ^doubles fs ^doubles gs]
(let [window-size (alength fs)
padding (repeat window-size 0.0)
padded-gs (concat padding gs padding)
window-of-gs (double-array (subvec (vec padded-gs)
(inc n)
(inc (+ n window-size))))
reverse-window-fn (fn [i] (aget window-of-gs
(mod (- -1 i) window-size)))
reverse-window (map reverse-window-fn (range window-size))]
(reduce + (map * fs reverse-window))))

(defn convolve [^doubles fs ^doubles gs]
(map #(point-convolve % fs gs)
(range (+ (alength fs) (alength gs)))))

(let [dxs (double-array (for [i (range 3000000)] (rand 100)))
dys (double-array (for [i (range 3000000)] (rand 100)))]
(time (dotimes [_ 10] (convolve dxs dys))))

"Elapsed time: 0.506432 msecs"

Have a good weekend,
Carson

Carson

unread,
Jul 18, 2010, 12:33:12 AM7/18/10
to Clojure
Thanks David.

On Jul 17, 4:58 pm, David Nolen <dnolen.li...@gmail.com> wrote:
> I tried this, it runs in about the same amount of time as it did for you.
>
> On Sat, Jul 17, 2010 at 7:17 PM, Carson <c.sci.b...@gmail.com> wrote:
> > Hi David,
> > Would appreciate if you could try out my attempt at this [1] on your
> > machine.  My machine ain't as fast as yours apparently...
>
> > [1]
>
> >https://groups.google.com/group/clojure/tree/browse_frm/thread/ee4169...
>
> > Thanks,
> > Carson
>
> > On Jul 17, 4:04 pm, David Nolen <dnolen.li...@gmail.com> wrote:
> > > (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))
>
> > > I don't think this is so bad and it can do a million points in ~400ms on
> > my
> > > machine. 100,000 points in ~25ms. I don't consider this dropping to the
> > Java
> > > level at all.
>
> > --
> > 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<clojure%2Bunsu...@googlegroups.com>

Per Vognsen

unread,
Jul 18, 2010, 1:02:43 AM7/18/10
to clo...@googlegroups.com
On Sun, Jul 18, 2010 at 11:32 AM, Carson <c.sci...@gmail.com> wrote:
> I guess different people see things differently.  I actually didn't
> understand the intent of the imperative version, even though it takes
> less lines I guess.  There's quite a few things called convolution.
> My naive implementation was functional and pretty fast.  Turns out it
> can be even faster with just a little work (still functional):
>
> (defn point-convolve [n ^doubles fs ^doubles gs]
>  (let [window-size (alength fs)
>        padding (repeat window-size 0.0)
>        padded-gs (concat padding gs padding)
>        window-of-gs (double-array (subvec (vec padded-gs)
>                                           (inc n)
>                                           (inc (+ n window-size))))
>        reverse-window-fn (fn [i] (aget window-of-gs
>                                        (mod  (- -1 i) window-size)))
>        reverse-window (map reverse-window-fn (range window-size))]
>    (reduce + (map * fs reverse-window))))
>
> (defn convolve [^doubles fs ^doubles gs]
>  (map #(point-convolve % fs gs)
>       (range (+ (alength fs) (alength gs)))))
>
> (let [dxs (double-array (for [i (range 3000000)] (rand 100)))
>      dys (double-array (for [i (range 3000000)] (rand 100)))]
>  (time (dotimes [_ 10] (convolve dxs dys))))
>
> "Elapsed time: 0.506432 msecs"

How about some common sense? The convolution algorithm everyone is
discussing runs in O(n m) time. So your measured time would imply that
a single of your computer's core is capable of 10 * 3,000,000^2 = 9 *
10^14 FLOPS, or about one peta-FLOPS. That should have tipped you off
you aren't actually measuring anything useful.

Wrap your (convolve dxs dys) calls with (doall ...). And you will want
to reduce the three million number to something much smaller, or
you'll be waiting for a very long time.

By the way, no-one actually uses this quadratic-time algorithm for
very wide kernels. For that you use FFTs and pointwise multiplication
(look up the convolution theorem), which yields an O(n log n)
algorithm.

-Per

Carson

unread,
Jul 18, 2010, 2:18:59 AM7/18/10
to Clojure
Hi Per,
woh, take it easy. I don't claim to be an expert. Thanks for showing
me that though. It certainly didn't seem right at first, but I had
trouble figuring out the laziness in clojure, me being new to it.

Anyway, have a good weekend!
Carson

Isaac Hodes

unread,
Jul 17, 2010, 10:46:13 PM7/17/10
to Clojure
First: sorry for splitting the conversation. That was my fault: don't
know how it happened. For those following along: the thread this
started with is here: http://groups.google.com/group/clojure/browse_thread/thread/ee4169bc292ab572

Second, I don't think I expressed myself well in the original post of
this thread.

The performance isn't a huge deal to me: I'm sorry if I left the
impression that this was a pressing issue. I will eventually be
implementing the FFT in Clojure (code to be posted if it's decent: I
think I might have a better time expressing that functionally than
this) and that speeds up convolution by orders of magnitude anyway.
I'll also, ideally, be parallelizing it. This was more of an exercise
for myself, and than it turned out to be so interesting that I thought
I'd pose it to the community, as my own solutions seemed lacking or
rather imperative.

I am impressed with the beauty of some of the functional solutions;
the issue is that the ones that work well are decidedly less
functional than would be ideal. Ideal only in that, well, it'd be nice
to have a functional solution that is speedy: I think it'd be neat.

The solutions here ranged from functional (slow), to procedural but
still in Clojure (could get rather fast) to basically Java (a little
faster than the latter).

I'm not unhappy with this: I learned a lot from the answers here, and
I got an idea of some limitation of purely functional Clojure (for
now: Stuart Sierra says this might get better soon?) I want to say
thank you to all of you who posted. I'll be posting more soon, about
the FFT, when the time comes: should be fun!

--individual responses--

@David Nolan:
Sent you a Twitter message–I am having trouble running your solution,
but that looks nicely Clojuresque and speedy. I was unaware of
`dotimes`.

@Michael Gardener:
You're right: it's idiomatic, though not functional. I misspoke there:
Clojure is, reluctantly it seems, multi-paradigm. And you're right, if
I just wanted raw, single-threaded, speed, I would go with C. But I'd
like to do this in Clojure!

@Carson:
I did check out your response: it's rather fast on my machine. it's
not really functional, though, as you use the `let` macro as a way of
procedurally executing a lot of functions. This isn't bad at all, but
you're not composing functions.

@j-g-faustus:
I agree with you completely: as I said above, it can be expressed in
the language, and can even be rather elegant. What you don't get is
speed, and not too much functional-style programming going on.

As for what's practical, yes, Clojure might not be it if all I want is
the raw speed. And translating imperative algorithms to functional
isn't always easy: but it's fun! And the other advantages that Clojure
has over C and its ilk make it really nice for the project(s) I'm
working on.

-----

Thank you all again!

Michael Gardner

unread,
Jul 18, 2010, 11:47:37 AM7/18/10
to clo...@googlegroups.com
On Jul 17, 2010, at 9:46 PM, Isaac Hodes wrote:

> I did check out your response: it's rather fast on my machine. it's
> not really functional, though, as you use the `let` macro as a way of
> procedurally executing a lot of functions. This isn't bad at all, but
> you're not composing functions.

'Functional' means simply avoiding mutation and side-effects. How the code is organized is a separate issue.

In case you're interested, I made a small change to the code I posted before that makes it significantly faster. I don't know what you consider 'slow', but it takes 2-3x as long as your imperative Python version on my hardware (less if I replace the outermost map with pmap, depending on the input sizes). I think that's quite acceptable, especially given Clojure's relative youth as a language.

(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."
(map #(vector % (- i %))
(range
(max 0 (inc (- i max-n)))
(min (inc i) max-m))))

(defn convolve-1 [i ms ns]
"Returns the ith value of the convolution of ms with ns."
(reduce +
(map (fn [[m n]] (* (ms m) (ns n)))
(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))))))

Jeffrey Schwab

unread,
Jul 18, 2010, 2:14:37 PM7/18/10
to clo...@googlegroups.com
On 7/17/10 4:43 PM, Isaac Hodes wrote:

Apologies in advance if this is somewhat OT in a Clojure group.

> double convolve(double *xs, double *is, double *ys){
> int i,j;
> for(i=0; i<len(xs); i++){
> for(j=0; j<len(is); j++){

How are you getting the size of the collections from the pointers, and
why are you checking the lengths on every loop pass (rather than once, a
priori)?

> ys[i+j] = ys[i+j] + (xs[i]*is[j]);

The following would be more idiomatic:

ys[i + j] += xs[i] * is[j];

Less repetition, less line noise, and in C++, it may avoid some
unnecessary temporaries if you ever abstract from raw ints and doubles
to something less likely to be elided by the compiler (e.g.
Clojure-style, overflow-detecting types).

> }
> }
> return ys;
> }

Reply all
Reply to author
Forward
0 new messages