Improve speed for calculating a recursive sequence

Question

I want to calculate following recursive sequence:

$\alpha_{0}=0,\\ \cos(\alpha_{i})=\cos(\alpha_{i-1})\cdot\cos(\beta_{i})+\sin(\alpha_{i-1})\cdot\sin(\beta_{i})\cdot\cos(\gamma_{i}).$

In mathematica 8.0.1 I try to do this with:

endRecursion = 20;
beta = Import["beta.txt", "List"];
beta = beta[[1;;endRecursion]];
gamma = RandomReal[{0,2*Pi}, endRecursion];

alpha[0]:= 0;
alpha[i_]:= 
   ArcCos[Cos[alpha[i-1]]*Cos[beta[[i]]]+ 
   Sin[alpha[i-1]]*Sin[beta[[i]]]*Cos[gamma[[i]]]];

result=Table[alpha[i], {i, 1, endRecursion}];

The values for beta are imported because they are generated by another formula which is not important for my question here. Values for beta can be found here.

The problem is that the calculation of the result is very slowly. It takes about 14s for just 20 iterations! And I need about 300 iterations (times ~10000, because I have several lists of beta values).

I recognized that the calculation is much faster if I skip the Cos[alpha[i-1]] in the first or the Sin[alpha[i-1]] in the 2nd summand. Then it takes just ~ 0.01 s for 20 iterations. But unfortunately this falsifies my result. I cannot understand why it is so much faster if there is just one time alpha[i-1] on the right side of the equation. Am I doing something wrong? Is there any way to increase the speed of this calculation? I would be glad about any help!

Answer 1

Although you can always implement recursion by hand the way you are trying to do, there are also some specialized functions that are designed to make your life a little easier. In particular, there is FoldList.

This approach is also slightly faster than the manual recursion approach (assuming you start from a clean slate):

endRecursion = 20;
beta = Import["beta.txt", "List"];
beta = beta[[1 ;; endRecursion]];
gamma = RandomReal[{0, 2*Pi}, endRecursion];
FoldList[ArcCos[Cos[#1]*#2[[1]] + Sin[#1]*#2[[2]]] &, 0,
 Transpose[{Cos[beta], Sin[beta] Cos[gamma]}]]

After initializing the lists as in your example, I need only a single statement to generate the list.

The folding function is what appears in the original recursive formulation, but it is written as a pure function with two formal arguments #1 and #2.

The first argument, #1, refers to your alpha and its starting value is passed to FoldList as the second argument (0). The second argument in the pure function (#2) refers to the elements of the list Transpose[{Cos[beta], Sin[beta] Cos[gamma]}] in which I collected the pre-existing lists beta and gamma (but in a form that is already pre-processed to avoid having to evaluate cosines and sines individually in the recursion).

Answer 2

You can speed it up by "memoization", i.e., remembering previous values of alpha[i]:

endRecursion = 20;
beta = Import["beta.txt", "List"];
beta = beta[[1 ;; endRecursion]];
gamma = RandomReal[{0, 2*Pi}, endRecursion];
alpha[0] := 0;
alpha[i_] := 
  alpha[i] = 
   ArcCos[Cos[alpha[i - 1]]*Cos[beta[[i]]] + 
     Sin[alpha[i - 1]]*Sin[beta[[i]]]*Cos[gamma[[i]]]];
result = Table[alpha[i], {i, 1, endRecursion}];

See also the docs on this.

Answer 3

The function runs slowly because each iteration calls alpha twice the number of times of the previous iteration. Thus, $n$ iterations invoke alpha $2^n$ times -- exponential run time! Fortunately, the two nested calls to alpha are identical. If we change alpha to call itself only once instead of twice, the algorithm runs in linear time while still retaining simple recursive form:

endRecursion = 20;
beta = Import["beta.txt", "List"];
beta = beta[[1;;endRecursion]];
gamma = RandomReal[{0,2*Pi}, endRecursion];

alpha[0] = 0;
alpha[i_]:=
  Module[{a = alpha[i-1], b = beta[[i]], g = gamma[[i]]}
  , ArcCos[Cos[a]*Cos[b]+Sin[a]*Sin[b]*Cos[g]]
  ]

result = Table[alpha[i], {i, 1, endRecursion}];

Even so, I agree with @Jens that functional iteration is frequently preferred over recursion in order to avoid hitting the $RecursionLimit (or $IterationLimit).