# Improve speed for calculating a recursive sequence

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!

-
Have you tried caching? alpha[i_] := alpha[i] = ArcCos[Cos[alpha[i-1]]*Cos[beta[[i]]] + Sin[alpha[i-1]]*Sin[beta[[i]]]*Cos[gamma[[i]]]] – J. M. May 26 '12 at 12:44
I am such a *** idiot! Thank you so much for this comment! I totally forgot the caching although I used it some time ago (5 years, perhaps too long time ago). Unfortunately, I did not see this trick when I searched for a solution, otherwise I would not have asked this question... Thank you anyway! You saved my weekend! How can I accept your comment as answer? – partial81 May 26 '12 at 12:58
@J.M. hmm, in my amazement that nobody had answered this, I wrote my answer without seeing your comment... Hope you were not writing up one also! – acl May 26 '12 at 13:07
You can accept acl's answer. :) I didn't answer because I can't test at the moment. – J. M. May 26 '12 at 13:07
Ok, if you do not mind! Thanks again! – partial81 May 26 '12 at 13:12

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}];


-
Thanks acl and @J.M. for answering my question so quickly! – partial81 May 26 '12 at 13:14
There was a mistake: you were importing strings (you didn't convert to numbers after StringSplit.) I fixed it. – Szabolcs May 26 '12 at 13:15
@Szabolcs yes, I had just realised it and our edits almost overlapped... thanks – acl May 26 '12 at 13:16
Okay, I'll remove the comments then. – Szabolcs May 26 '12 at 13:20

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).

-
With this particular formulation, one could even simplify the FoldList[] expression as FoldList[ArcCos[{Cos[#1], Sin[#1]}.#2] &, 0, Transpose[{Cos[beta], Sin[beta] Cos[gamma]}]]. One might even consider using the Weierstrass substitution here, since the beta values are all $< \pi/2$... – J. M. May 26 '12 at 18:47
@J.M. That's correct - the dot product may give a tiny additional speed boost but I didn't check that. – Jens May 27 '12 at 0:50
Thanks @Jens for this additional answer. FoldList is something I have never used so far – that will change now (although the additional speed boost is just small). – partial81 May 27 '12 at 11:24
@ Jens, Let’s say: alpha[i_]:=alpha[i]=ArcCos[Cos[alpha[i - 1]]*Cos[beta[[i+1]]]+Sin[alpha[i - 1]]*Sin[beta[[i]]]*Cos[gamma[[i]]]]. How do I have to adapt your solution with FoldList for this case? I really have problems to see a solution because now I would need in my iteration the i+1 element of beta. I would be happy if you could give me a hint again! – partial81 May 31 '12 at 13:10
@partial81 The term you changed is the first one in Transpose[{Cos[beta], Sin[beta] Cos[gamma]}], and the change is simply a shift in the index. The fastest way this can be implemented is to use RotateLeft: Transpose[{RotateLeft[Cos[beta]], Sin[beta] Cos[gamma]}] – Jens May 31 '12 at 15:39

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).

-
Thanks @WReach for answering the other part of my question. Now I understand why it was so slowly at the beginning. And it is nice that you give a link to the useful functional iterations. – partial81 May 27 '12 at 20:53