Re: Why is recursion so slow in Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg100568] Re: Why is recursion so slow in Mathematica?
- From: Daniel <janzon at gmail.com>
- Date: Mon, 8 Jun 2009 02:07:03 -0400 (EDT)
- References: <h0d6s8$spq$1@smc.vnet.net> <h0fvn1$ree$1@smc.vnet.net>
On 7 Juni, 11:02, Szabolcs Horv=E1t <szhor... at gmail.com> wrote:
> Daniel 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=
ist
> > (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=
ica
> > *)
>
> So you discovered that appending to arrays is slow while appending to
> linked lists is fast :)
>
> Actually you are not comparing the exact same algorithms. The data
> structure you used in those languages was a linked list. Removing or
> appending elements takes a constant time for linked lists.
> Mathematica's List is a vector-type data structure, so
> removing/appending elements takes a time proportional to the vector's
> length.
>
> One data structure is not better than the other, of course, they're just
> useful for different purposes. For example, a linked list is unsuitabl=
e
> for applications where the elements need to be accessed randomly instead
> of sequentially. Random access is often needed for numerical
> computations/simulations.
>
> For a fair comparison, use a linked list in Mathematica too. This coul=
d
> look like {1,{2,{3,{}}}}
>
> myselect[{}, test_] = {};
> myselect[{head_, tail_}, test_] :=
> If[test[head], {head, myselect[tail, test]}, myselect[tail, test]]
>
> toLinkedList[list_List] := Fold[{#2, #1} &, {}, list]
>
> data = toLinkedList@Table[Random[], {20000}];
>
> Block[{$RecursionLimit = \[Infinity]},
> Timing[myselect[data, # > 0.5 &];]]
>
> This runs in 0.2 sec on my (not very fast) system.
>
> And you can check that the function works correctly:
>
> Block[{$RecursionLimit = \[Infinity]},
> Flatten@myselect[data, # > 0.5 &] ===
> Select[Flatten[data], # > 0.5 &]]
>
> About the thing that you call "naive C-like behavior": you probably mean
> that Mathematica doesn't automatically perform the tail-call
> optimization. That is correct, you need to transform your recursions
> explicitly.
>
> Here's an example:
>
> myselect2[dest_, {head_, tail_}, test_] :=
> If[test[head],
> myselect2[{head, dest}, tail, test],
> myselect2[dest, tail, test]]
>
> myselect2[dest_, {}, test_] := dest
>
> Note that it is $IterationLimit that we need to increase now and this
> version doesn't fill up the evaluation stack.
>
> Also note that this produces the result in the reverse order:
>
> Block[{$IterationLimit = \[Infinity]},
> Reverse@Flatten@myselect2[{}, data, # > 0.5 &] ===
> Select[Flatten[data], # > 0.5 &]]
>
> We could produce the result in the same order, but reverse linking (like
> {{{{},1},2},3}).
>
> I couldn't do any better than this using a singly-linked list. It look=
s
> like OCaml (which I don't know BTW) uses a singly linked list too, and
> can't do the tail-call optimization, so its stack gets filled up.
>
> One more thing: even when using the same data structure and same
> algorithm, it is of course expected that a very high level interpreted
> language like Mathematica is not going to perform nearly as well as a
> low level compiled language like OCaml. In high level languages the
> usual solution to this performance problem is built-in functions: don't
> implement Select in Mathematica itself. Use the built-in one, which is
> written in C and much faster than a Mathematica implementation could ever=
be.
>
> Hope this helps,
> Szabolcs
Oh, so they are arrays! Does that mean Mathematica has to reallocate
memory to fit in new elements?
Just a check, what's the space requirment for your representations of
linked list?
ByteCount[Fold[{#2,#1}&,{},Table[RandomReal[],{1000}]]
48024 bytes for 1000 eight-byte doubles. Same byte count for
RandomInteger, BTW! In C I guess it would be something like 8*1000
bytes for the data plus 4*1000 bytes for the pointers to the next list
elment. The Mathematica overhead seems reasonably ok to me, and it
seems to be linear with the number of elements.
Tail-recursion was a good tip, thanks a lot! :) Just for fun I tried
a tail-recursive version of the standard list implementation of
myselect (i.e. using arrays). The time performance is almost the same,
but the tail-recursive versions doesn't eat memory like a monster.
In summary, with your help I feel confident to continue my journey
into Mathematica land :) I hope one day I will be able to contribute
back to this group and the community.
Regards,
Daniel