Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
5  
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
add comment

3 Answers

up vote 9 down vote accepted

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.

share|improve this answer
    
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
add comment

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

share|improve this answer
1  
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
1  
@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
show 1 more comment

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

share|improve this answer
    
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
add comment

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.