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?

Improve this question
$\endgroup$
12
  • $\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, 2016 at 19:40
  • 3
    $\begingroup$ I can make it slower by changing Xmn[m] to Xmn[m,n]. $\endgroup$
    – wxffles
    Jan 18, 2016 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, 2016 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, 2016 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, 2016 at 20:57

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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