Re: Why is recursion so slow in Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg100569] Re: Why is recursion so slow in Mathematica?
- From: Daniel <janzon at gmail.com>
- Date: Mon, 8 Jun 2009 02:07:14 -0400 (EDT)
- References: <h0d6s8$spq$1@smc.vnet.net> <h0fvo1$rfd$1@smc.vnet.net>
Leonid,
Thanks a lot for your reply. It seems like Wagner's book (Power
Programming in Mathematica) is out of print, too bad! But your book
seems to be interesting as well. I will definately read it the whole
thing; the section on why you think Mathematica is so good was
interesting. I guess I need to learn about reap/sow as well...
I do understand that Mathematica is what it is and is not supposed to
compete with low-level languages on speed.
All the best,
Daniel
On 7 Juni, 11:02, lsh... at gmail.com wrote:
> Daniel,
>
> Let me tell you straight away that you won't get efficiency of OCaml
> in Mathematica,
> the latter really being a CAS, research and prototyping tool, but not
> an efficiency-oriented production language.
>
> Regarding your question, the short answer is that recursion is slow
> because lists are implemented as arrays in Mathematica. When you
> pattern-match, the whole array <tail> is copied every time. This
> explains both the slowdown - the time becomes effectively quadratic
> (assuming that the complexity of array copying is linear), and the
> (stack) memory use explosion.
>
> So, IMO, indeed Mathematica is unfortunately less suited for the
> straightforward recursive solutions like yours. The closest thing you
> can probably do to stay in the recursive mindset is to use linked
> lists. Here is an implementation that does just that:
>
> Clear[myselectAux];
> myselectAux[{}, predicate_] = {};
> myselectAux[{head_, tail_List}, predicate_] :=
> If[predicate[head], {head, myselectAux[tail, predicate]},
> myselectAux[tail, predicate]];
>
> Clear[toLinkedList];
> toLinkedList[x_List] := Fold[{#2, #1} &, {}, Reverse[x]];
>
> Clear[mySelectLL];
> mySelectLL[x_List, predicate_] :=
> Flatten@myselectAux[toLinkedList[x], predicate];
>
> You can find more info on linked lists in Mathematica in the book of
> David Wagner,
> in the notes by Daniel Lichtblau (see Wolfram library archive, should
> be named
> something like "efficient data structures in Mathematica"), and also
> in my book
> (http://www.mathprogramming-intro.org/book/node525.html).
>
> For example:
>
> In[1] = test = RandomInteger[{1, 10}, 10]
>
> Out[1] = {6, 10, 1, 9, 4, 6, 1, 9, 10, 4}
>
> In[2] = toLinkedList@test
>
> Out[2] = {6, {10, {1, {9, {4, {6, {1, {9, {10, {4, {}}}}}}}}}}}
>
> In[3] = mySelectLL[test, # > 5 &]
>
> Out[3] = {6, 10, 9, 6, 9, 10}
>
> Now, the benchmarks: I use yours, on my 5 years old 1.8 GHz P IV, and
> I get:
>
> In[4] =
> data = Table[Random[], {20000}];
> Block[{$RecursionLimit = Infinity},
> Timing[data2 = mySelectLL[data, # > 0.5 &];]]
>
> Out[4] = {0.2, Null}
>
> Now, this is not the most efficient implementation. David Wagner in
> his excellent book
> discussed in great detail this sort of problems with the mergesort
> algorithm as an archetypal example. I believe that you can squeeze
> another factor of 2-5 by further optimizing the implementation
> (keeping it recursive) with several performance-tuning techniques
> available in Mathematica.
>
> By the way, I tested the above code with larger lists and it seems to
> continue be linear complexity up to 100 000 elements, but the kernel
> (v.6) just crashes after that - this
> behavior I can not explain, since the memory usage seems quite decent.
>
> So, to summarize my opinion: if you prefer to stay in a recursive
> mode, you *can* use recursion in Mathematica rather efficiently by
> switching to linked lists (and this will require some extra work but
> not so much of it). This may be enough to move you from the domain of
> toy problems to that of reasonably complex prototypes. But don't
> expect efficiency of the highly optimized production language with a
> native-code compiler (or even byte code interpreter) from a rewrite
> system such as Mathematica. Also, I'd say that other programming
> styles such as structural operations (APL - like) will bring you
> closer to efficiency of Ocaml etc. since then you can leverage the
> highly-optimized vectorized operations built in Mathematica.
>
> Hope this helps.
>
> Regards,
> Leonid
>
> On Jun 6, 12:46 am, Daniel <jan... at gmail.com> wrote:
>
> > This post is about functional programming in Mathematica versus other
> > functional languages such as OCaml, SML or Haskell. At least a naive
> > use of functional constructs in Mathematica is horrendously slow. Am I
> > doing something wrong? Or isn't Mathematica really suitable for
> > functional programming beyond toy programs? Couldn't the Wolfram team
> > make a more efficient implementation for recursion, as other
> > functional languages has done? (Instead of the naive C-like behavior
> > for recursively defined functions.)
>
> > As grounds for my question/argument, I wrote my own version of select,
> > as below
>
> > myselect[{}, predicate_] = {}
> > myselect[{head_, tail___}, predicate_] := If[predicate[head],
> > Join[{head}, myselect[{tail}, predicate]],
> > myselect[{tail}, predicate]
> > ]
>
> > Then I tried this function on a 20.000 element vector with machine
> > size floats:
>
> > data = Table[Random[], {20000}];
> > $RecursionLimit = 100000;
> > Timing[data2 = myselect[data, # > 0.5 &];]
>
> > The result is {7.05644, Null}, and hundreds of MB of system memory are
> > allocated. On 1.7 GHZ dual core Intel machine with 1 GB of RAM. For
> > 20.000 floats! It's just a megabyte!
>
> > The following OCaml program executes in apparently no-time. It is not
> > compiled and does the same thing as the above Mathematica code. After
> > increasing the list by a factor of ten to 200.000 elements, it still
> > executes in a blink. (But with 2.000.000 elements my OCaml interpreter
> > says Stack overflow.)
>
> > let rec randlist n = if n=0 then [] else Random.float(1.0) :: randl=
is=
> t
> > (n-1);;
>
> > let rec myselect = function
> > [],predicate -> []
> > | x::xs,predicate -> if predicate(x) then x::myselect(xs,predica=
te=
> )
> > else myselect(xs,predicate);;
>
> > let mypred x = x>0.5;;
>
> > let l=randlist(20000);;
> > let l2=myselect(l,mypred);; (* lightning-fast compared to Mathemat=
ic=
> a
> > *)