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