How to improve ( mathematica ) performance when dealing with recursive functions?

Question

Background. I want to print a table of convergents for 31^(1/2). I made the following recursive definition of the table. ( Exchange 31^(1/2) with the golden ratio and the table below would contain the Fibonacci series ).

 cf := ContinuedFraction
 tf := TableForm
 p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
 a[k_] := cf[Sqrt[31], k][[k]]
 p[k_] := a[k]*p[k - 1] + p[k - 2]
 q[k_] := a[k]*q[k - 1] + q[k - 2]
 s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}]] // tf

Timing increases exponentially fast. I had to alt+. (abort) for s[4].

Question: How to improve ( mathematica ) performance when dealing with recursive functions?

Well, I have no idea how mathematica works but dynamic programming. (Search it up) — flight, Sep 01 '11 at 09:46
The appropriate page of Mathematica documentation for this is [Functions That Remember Values They Have Found](http://reference.wolfram.com/mathematica/tutorial/FunctionsThatRememberValuesTheyHaveFound.html). — Simon, Sep 01 '11 at 10:40
Also see the SO questions: [Performance of Fibonacci](http://stackoverflow.com/questions/4130161/performance-of-fibonacci) and the good discussions in [The best way to construct a function with memory](http://stackoverflow.com/questions/5287817/the-best-way-to-construct-a-function-with-memory) — Simon, Sep 01 '11 at 10:42
Related: "[Tail call optimization in Mathematica](http://stackoverflow.com/questions/4481301/tail-call-optimization-in-mathematica)." — Alexey Popkov, Sep 01 '11 at 13:38

Szabolcs · Accepted Answer · 2011-09-01T09:49:56.173

8

From a quick (not thorough, to admit it) look at your code, it looks like both p and q are defined recursively in terms of two previous values. This means that to calculate the nth value of p, ~2^n evaluations are needed (every step doubles the number). So yes, complexity is exponential, regardless of Mathematica or any other language being used.

If you insist on using a recursive formulation of the problem (e.g. for simplicity), then the simplest way to reduce the performance penalty is using memoization, i.e. doing something like

p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]

Don't forget to Clear[p] before any redefinition.

In short, memoization means having the function remember the computation result for each input, so subsequent evaluations are faster. It's very likely faster, but more complicated, to compute two values (p_(n+1) and p_(n)) from two previous values (p_(n) and p_(n-1)), then the complexity will be linear instead of exponential.

_{I hope this helps. I don't have Mathematica here to test right now.}

edited Sep 01 '11 at 09:49

answered Sep 01 '11 at 09:44

Szabolcs

24,728
9
85
174

+1 for Memoization, although I suspect it would be well worth looking for an iterative solution to your problem instead. e.g. for Fibonacci sequence it's much much faster to build up from the bottom than to use recursion (even with memoization). – ForbesLindesay Sep 01 '11 at 09:50
@Tuskan360, certainly, that's faster (also see my edit). The point was that memoization the minimal modification to the existing implementation that makes the performance very usable. The "iterative" solution is equivalent to using pairs of values for calculation (as in my edit). This can be done using recursion, `Nest`, or procedural programming. – Szabolcs Sep 01 '11 at 09:52
@ Szabolcs - That worked great! s[4] down from minutes to milliseconds. – nilo de roock Sep 01 '11 at 10:07

Daniel Lichtblau · Answer 2 · 2011-09-02T18:09:18.403

Here is a small further refinement. Since this is a quadratic irrational you can also compute the a[k] coefficients more directly.

In[499]:= Clear[a, p, q, cf]
cf = ContinuedFraction[Sqrt[31]];
cf2len = Length[cf[[2]]];
a[1] = cf[[1]];
a[k_] := cf[[2, Mod[k - 1, cf2len, 1]]]
p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
q[k_] := q[k] = a[k]*q[k - 1] + q[k - 2]
s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}];]

In[508]:= s[1000]

Out[508]= {0.12, Null}

In[509]:= Clear[a, p, q, cf]
cf := ContinuedFraction
p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
a[k_] := a[k] = cf[Sqrt[31], k][[k]]
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
q[k_] := q[k] = a[k]*q[k - 1] + q[k - 2]
s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}];]

In[516]:= s[1000]

Out[516]= {6.08, Null}

Also you can get a[k] in closed form, though it is not terribly pretty.

In[586]:= Clear[a];
asoln[k_] = 
 FullSimplify[
  a[k] /. First[
    RSolve[Join[
      Table[a[k] == cf[[2, Mod[k - 1, cf2len, 1]]], {k, 
        cf2len}], {a[k] == a[k - 8]}], a[k], k]], Assumptions -> k > 0]

Out[587]= (1/(8*Sqrt[2]))*(4*(Cos[(k*Pi)/4] + Sin[(k*Pi)/4])*
        (-2*Sqrt[2] + (5 + 2*Sqrt[2])*Sin[(k*Pi)/2]) + 
      Sqrt[2]*(25 - 9*Cos[k*Pi] + 26*Sin[(k*Pi)/2] - 9*I*Sin[k*Pi]))

Offhand I do not know whether this can be used to get a direct solution for p[k] and q[k]. RSolve seems unable to do that.

--- edit ---

As others have mentioned, it can be cleaner to just build the list from first to last. Here is the handling of p[k], using memoization as above vs NestList.

Clear[a, p, q, cf]
cf = ContinuedFraction[Sqrt[31]];
cf2len = Length[cf[[2]]];
a[1] = cf[[1]];
a[k_] := cf[[2, Mod[k - 1, cf2len, 1]]]

p[-1] = 0; p[0] = 1;
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
s[n_] := Timing[Table[p[k], {k, n}];]

In[10]:= s[100000]

Out[10]= {1.64, Null}

In[153]:= s2[n_] := Timing[ll = Module[{k = 0},
      NestList[(k++; {#[[2]], a[k]*#[[2]] + #[[1]]}) &, {0, 1}, 
       n]][[All, 2]];]

In[154]:= s2[100000]

Out[154]= {0.78, Null}

In addition to being somewhat faster, this second approach does not keep a large number of definitions around. And you do not really need them in order to generate more elements, because this iteration can be resumed using a pair from the last elements (make sure they start at 0 and 1 modulo 8).

I will mention that one can obtain a closed form for p[k]. I found it convenient to break the solution into 8 (that is, cf2len) pieces and link them via recurrences. The reasoning behind the scenes comes from basic generating function manipulation. I did some slightly special handling of one equation and one initial condition to finesse the fact that a[1] is not part of the repeating sequence.

In[194]:= func = Array[f, cf2len];
args = Through[func[n]];
firsteqns = {f[2][n] == a[2]*f[1][n] + f[cf2len][n - 1], 
   f[1][n] == a[9]*f[cf2len][n - 1] + f[cf2len - 1][n - 1]};
resteqns = 
  Table[f[j][n] == a[j]*f[j - 1][n] + f[j - 2][n], {j, 3, cf2len}];
inits = {f[8][0] == 1, f[1][1] == 5};
eqns = Join[firsteqns, resteqns, inits];

In[200]:= 
soln = FullSimplify[args /. First[RSolve[eqns, args, n]], 
   Assumptions -> n > 0];

In[201]:= FullSimplify[Table[soln, {n, 1, 3}]]

Out[201]= {{5, 6, 11, 39, 206, 657, 863, 1520}, {16063, 17583, 33646, 
  118521, 626251, 1997274, 2623525, 4620799}, {48831515, 53452314, 
  102283829, 360303801, 1903802834, 6071712303, 7975515137, 
  14047227440}}

Quick check:

In[167]:= s2[16]; ll

Out[167]= {1, 5, 6, 11, 39, 206, 657, 863, 1520, 16063, 17583, 33646, \
118521, 626251, 1997274, 2623525, 4620799}

We can now define a function from this.

In[165]:= 
p2[k_Integer] := soln[[Mod[k, cf2len, 1]]] /. n -> Ceiling[k/cf2len]

In[166]:= Simplify[p2[4]]

Out[166]= 39

I do not claim that this is particularly useful, just wanted to see if I could actually get something to work.

--- end edit ---

Daniel Lichtblau

Creating a closed form for a[k] with recurrences / generating functions is an excellent idea. Never seen it thus far in the literature on continued fractions. — nilo de roock, Sep 02 '11 at 20:12
@nilo I'm guessing you meant p[k]. The recurrence for a[k] (disregarding the first element) is simple since it is a repeating sequence: a[k] = a[k-8]. The generating function is not too hard to obtain. In notation I am using it should be Sum[x^k*cf[[2, k]]/(1 - x^cf2len), {k, cf2len}]. But I was not able to find a way to use this to get p or q as generating functions. — Daniel Lichtblau, Sep 02 '11 at 20:42
p[k]. Of course. - Well, generalized the sequence of convergents of the continued fraction of sqrt[n] in relation to solving the pell equation x^2-n*y^2=1 with (x,y) in N. - Interesting topic. Could spent months on it, if I hadn't so many competing activities. ;-) I just found this btw. Not sure if it is relevant. http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P68.pdf - — nilo de roock, Sep 03 '11 at 07:45

How to improve ( mathematica ) performance when dealing with recursive functions?

Question: How to improve ( mathematica ) performance when dealing with recursive functions?

2 Answers2

Linked