2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Oleksandr R. Jan 18 '16 at 19:40
  • 3
    $\begingroup$ I can make it slower by changing Xmn[m] to Xmn[m,n]. $\endgroup$ – wxffles Jan 18 '16 at 19:56
  • $\begingroup$ 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. $\endgroup$ – Anton Antonov Jan 18 '16 at 20:30
  • 3
    $\begingroup$ @OleksandrR. Compile has handled recursion (with some help needed at times for type inferencing) since around version 7 or maybe 8. $\endgroup$ – Daniel Lichtblau Jan 18 '16 at 20:33
  • $\begingroup$ @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. $\endgroup$ – Oleksandr R. Jan 18 '16 at 20:57
7
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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