Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a recursive function defined like this:

Clear[MyFnc];
MyFnc[X_] := Block[{n, val, Xm, Xmn},
    n = Length[X]; If[n == 1, Return[{2, 1}]];
    Xm[mm_] := X[[1 ;; mm]];
    Xmn[mm_, nn_] := X[[mm + 1 ;; nn]];
    val = Table[Block[{XmT = Total[Xm[m]], XmnT = Total[Xmn[m, n]]},
          (XmT.XmnT)/(XmT.XmT)*{m^2, 2 m}*(MyFnc[Xm[m]][[1]])*(MyFnc[Xmn[m,n]][[2]])],
          {m, 1, n - 1}] // Total;
    Return[val/{n^2, n}]];

Function takes a list of vectors of arbitrary length $n$, e.g. for $n=5$ and 3D vectors;

arg := Table[RandomReal[{-1, 1}, 3], {5}]

and gives a list of two numbers, i.e. MyFnc[arg] gives for example {0.24443, 1.10547}.

Since this function is to be used as an integrand in numerical integration, it would need to be called many times. So evaluation time is important,

(Table[MyFnc[arg], {10^4}] // AbsoluteTiming)[[1]]  

(* ==> 18.5959 *)

Is there a way to significantly speed up this function? Can such recursive functions be compiled efficiently?

share|improve this question
    
It is not possible to compile recursive functions (in general) because the virtual machine does not support recursion (for some reason). However, if the function can be adjusted to become iterative instead of recursive, then it can be compiled. – Oleksandr R. Jan 18 at 19:40
3  
I can make it slower by changing Xmn[m] to Xmn[m,n]. – wxffles Jan 18 at 19:56
    
It might be a good idea to re-write this function so the argument X is global, and the recursive calls are over the row indices of X. – Anton Antonov Jan 18 at 20:30
3  
@OleksandrR. Compile has handled recursion (with some help needed at times for type inferencing) since around version 7 or maybe 8. – Daniel Lichtblau Jan 18 at 20:33
    
@DanielLichtblau then why does With[{fibonacci = Function[Null, If[#1 == 1 || #1 == 2, 1, #0[#1 - 1] + #0[#1 - 2]]]}, Compile[{{i, _Integer, 0}}, fibonacci[i]]] crash the kernel? Apparently Compile is trying to unroll this recursion (the only way the VM can handle it apart from the opcode 43 mechanism, AFAIK) and gets into difficulty. If you are talking about "recursively call out of the VM and instantiate a new VM instance for every call", yes it can handle that, I agree, but I don't think it's useful in most cases because of the huge overhead per call. – Oleksandr R. Jan 18 at 20:57
up vote 7 down vote accepted

You can get it through Compile as below. Note that I have not tested for correctness.

myFncC = Compile[{{X, _Real, 2}}, Block[
    {n, val},
    n = Length[X];
    If[n == 1, Return[{2., 1.}]];
    val = Total[
      Table[
       Block[
        {XmT = Total[X[[1 ;; m]]], XmnT = Total[X[[m + 1 ;; n]]], 
         mFm = myFncC[X[[1 ;; m]]], mFmn = myFncC[X[[m + 1 ;; n]]]},
        (XmT.XmnT)/(XmT.XmT)*{m^2, 2 m}*mFm[[1]]*mFmn[[2]]], {m, 1, 
        n - 1}]];
    val/{n^2, n}]];

arg := Table[RandomReal[{-1, 1}, 3], {5}]
(Table[myFncC[arg], {10^4}] // AbsoluteTiming)[[1]]

(* Out[67]= 1.13199 *)

Adding "CompilationTarget" -> "C" brings it down a hair more, to .75 seconds.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.