Muller's Recurrence

The recurrence devised by Jean-Michel Muller[1] and described in a slightly modified form by William Kahan in a 2006 monograph[2] provides an interesting example of how the finite precision of floating point arithmetic can lead to wildly incorrect results in some circumstances.

Consider the function $$E(y, z) = 108 - \frac{815 - \frac{1500}{z}}{y}$$ and define the recurrence $x_{n+1} = E(x_n, x_{n-1})$, with initial values $x_0=4, x_1=4.25$. Does this recurrence have a limit (that is, does $x_n \rightarrow L$ for some finite $L$ as $n\rightarrow\infty$), and if so what is it? The mathematical analysis described below demonstrates that the limit does exist and that $L=5$. A numerical approach would be to calculate $x_n$ for increasing $n$ until $|x_n - x_{n-1}| < \epsilon$, for some tiny value of $\epsilon$. This can be implemented with floating point arthimetic, and the recurrence does indeed seem to converge by about $n=80$. However, it converges on the wrong limit: $L=100$.

Python has a library module, fraction for exact rational arithmetic, and this can be used to find the correct limit. The code below demonstrates this.; the output columns are the rational arithmetic ("exact") result, the naive floating point implementation result, and the output of a more numerically-stable floating point algorithm derived below.

import numpy as np
from fractions import Fraction

def E_float(y, z):
"""Naive implementation of Muller's recurrence."""
return 108 - (815 - 1500 / z) / y

def E_float2(y, _):
"""Numerically stable implementation of Muller's recurrence."""
return 8 - 15/y

def E_fraction(y, z):
"""Muller's recurrence implemented with rational arithmetic."""
return 108 - Fraction((815 - Fraction(1500, z)), y)

def do_recurrence(E, x):
for n in range(2, N+1):
x[n] = E(x[n-1], x[n-2])

N = 80
x_fraction = [0]*(N+1)
x_fraction[0], x_fraction[1] = 4, Fraction(17, 4)
do_recurrence(E_fraction, x_fraction)

x_float = [0]*(N+1)
x_float[0], x_float[1] = 4, 4.25
do_recurrence(E_float, x_float)

x_float2 = [0]*(N+1)
x_float2[0], x_float2[1] = 4, 4.25
do_recurrence(E_float2, x_float2)

print(' n {:>22s} {:^22s} {:^22s}'.format('rational arithmetic',
'floating point 1', 'floating point 2'))
print('--', '-'*22, '-'*22, '-'*22)
for n in range(N+1):
print('{:2d} {:22.15f} {:22.15f} {:22.15f}'.format(n, float(x_fraction[n]),
x_float[n], x_float2[n]))


Output (note that for clarity the rational expression has been converted back to the nearest representable floating point number in the first column):

 n    rational arithmetic    floating point 1       floating point 2
-- ---------------------- ---------------------- ----------------------
0      4.000000000000000      4.000000000000000      4.000000000000000
1      4.250000000000000      4.250000000000000      4.250000000000000
2      4.470588235294118      4.470588235294116      4.470588235294118
3      4.644736842105263      4.644736842105218      4.644736842105264
4      4.770538243626063      4.770538243625083      4.770538243626063
5      4.855700712589074      4.855700712568563      4.855700712589074
6      4.910847499082793      4.910847498660630      4.910847499082793
7      4.945537404123916      4.945537395530508      4.945537404123916
8      4.966962581762701      4.966962408040999      4.966962581762701
9      4.980045701355631      4.980042204293014      4.980045701355632
10      4.987979448478392      4.987909232795786      4.987979448478393
11      4.992770288062068      4.991362641314552      4.992770288062069
12      4.995655891506634      4.967455095552268      4.995655891506635
13      4.997391268381344      4.429690498308830      4.997391268381344
14      4.998433943944817     -7.817236578459315      4.998433943944817
15      4.999060071970894    168.939167671064581      4.999060071970893
16      4.999435937146839    102.039963152059272      4.999435937146838
17      4.999661524103767    100.099947516249699      4.999661524103767
18      4.999796900713418    100.004992040972439      4.999796900713418
19      4.999878135477931    100.000249579237305      4.999878135477931
20      4.999926879504600    100.000012478620164      4.999926879504599
21      4.999956127061158    100.000000623921608      4.999956127061157
22      4.999973676005713    100.000000031195796      4.999973676005713
23      4.999984205520272    100.000000001559783      4.999984205520272
24      4.999990523282228    100.000000000077989      4.999990523282227
25      4.999994313958560    100.000000000003894      4.999994313958560
26      4.999996588371256    100.000000000000199      4.999996588371256
27      4.999997953021357    100.000000000000014      4.999997953021357
28      4.999998771812312    100.000000000000000      4.999998771812312
29      4.999999263087206    100.000000000000000      4.999999263087206
30      4.999999557852258    100.000000000000000      4.999999557852258
31      4.999999734711332    100.000000000000000      4.999999734711332
32      4.999999840826790    100.000000000000000      4.999999840826790
33      4.999999904496072    100.000000000000000      4.999999904496072
34      4.999999942697642    100.000000000000000      4.999999942697642
35      4.999999965618585    100.000000000000000      4.999999965618585
36      4.999999979371150    100.000000000000000      4.999999979371150
37      4.999999987622690    100.000000000000000      4.999999987622690
38      4.999999992573614    100.000000000000000      4.999999992573614
39      4.999999995544169    100.000000000000000      4.999999995544169
40      4.999999997326501    100.000000000000000      4.999999997326501
41      4.999999998395901    100.000000000000000      4.999999998395901
42      4.999999999037541    100.000000000000000      4.999999999037541
43      4.999999999422524    100.000000000000000      4.999999999422524
44      4.999999999653514    100.000000000000000      4.999999999653514
45      4.999999999792109    100.000000000000000      4.999999999792109
46      4.999999999875265    100.000000000000000      4.999999999875265
47      4.999999999925159    100.000000000000000      4.999999999925159
48      4.999999999955095    100.000000000000000      4.999999999955095
49      4.999999999973057    100.000000000000000      4.999999999973057
50      4.999999999983834    100.000000000000000      4.999999999983834
51      4.999999999990300    100.000000000000000      4.999999999990301
52      4.999999999994181    100.000000000000000      4.999999999994181
53      4.999999999996509    100.000000000000000      4.999999999996509
54      4.999999999997905    100.000000000000000      4.999999999997906
55      4.999999999998743    100.000000000000000      4.999999999998743
56      4.999999999999246    100.000000000000000      4.999999999999246
57      4.999999999999547    100.000000000000000      4.999999999999547
58      4.999999999999728    100.000000000000000      4.999999999999728
59      4.999999999999837    100.000000000000000      4.999999999999837
60      4.999999999999902    100.000000000000000      4.999999999999902
61      4.999999999999941    100.000000000000000      4.999999999999941
62      4.999999999999964    100.000000000000000      4.999999999999964
63      4.999999999999979    100.000000000000000      4.999999999999979
64      4.999999999999988    100.000000000000000      4.999999999999988
65      4.999999999999992    100.000000000000000      4.999999999999993
66      4.999999999999996    100.000000000000000      4.999999999999996
67      4.999999999999997    100.000000000000000      4.999999999999997
68      4.999999999999998    100.000000000000000      4.999999999999998
69      4.999999999999999    100.000000000000000      4.999999999999999
70      4.999999999999999    100.000000000000000      5.000000000000000
71      5.000000000000000    100.000000000000000      5.000000000000000
72      5.000000000000000    100.000000000000000      5.000000000000000
73      5.000000000000000    100.000000000000000      5.000000000000000
74      5.000000000000000    100.000000000000000      5.000000000000000
75      5.000000000000000    100.000000000000000      5.000000000000000
76      5.000000000000000    100.000000000000000      5.000000000000000
77      5.000000000000000    100.000000000000000      5.000000000000000
78      5.000000000000000    100.000000000000000      5.000000000000000
79      5.000000000000000    100.000000000000000      5.000000000000000
80      5.000000000000000    100.000000000000000      5.000000000000000


Analysis

This section fills in some of the gaps in Kahan's explanation.

Let $x_n = \frac{y_{n+1}}{y_n}$. Then, $$\frac{y_{n+2}}{y_{n+1}} = 108 - \frac{815 - \frac{1500}{y_n}y_{n-1}}{y_{n+1}}y_n.$$ Therefore, $$y_{n+2} = 108 y_{n+1} - 815y_n + 1500y_{n-1}$$ This is a linear recurrence which may be written in closed form using the roots of its characteristic polynomial: $$P(z) = z^3 - 108z^2 + 815z - 1500 = (z-3)(z-5)(z-100).$$ which implies that $$y_n = \alpha 3^n + \beta 5^n + \gamma 100^n$$ and so $$x_n = \frac{\alpha 3^{n+1} + \beta 5^{n+1} + \gamma 100^{n+1}}{\alpha 3^n + \beta 5^n + \gamma 100^n}.$$ For the prescribed $x_0=4$ we have $4 = (3\alpha + 5\beta + 100\gamma)/(\alpha+\beta+\gamma)$. Clearly $\alpha$ and $\beta$ cannot both be zero. Without loss of generality (we have three unknowns but two initial conditions), choose $\alpha=1$. Then we have $$x_0 = 4 = \frac{3 + 5\beta + 100\gamma}{1+\beta+\gamma}, \quad x_1 = \frac{17}{4} = \frac{9+25\beta+10000\gamma}{3+5\beta+100\gamma},$$ from which it is easy to show that $\beta=1$ and $\gamma=0$. Therefore, \begin{align*} x_n &= \frac{3^{n+1}+5^{n+1}}{3^n + 5^n} = \frac{(5-2)\cdot3^n + 5\cdot 5^n}{3^n+5^n}\\ &=5 - 2\frac{3^n}{3^n+5^n} = 5 - \frac{2}{1+\left(\frac{5}{3}\right)^n} \end{align*} Clearly $\lim_{n\to \infty}x_n = 5$. This expression can help cast the recurrence back into a more numerically stable form, since $$x_{n-1} = 5 - \frac{2}{1+\left(\frac{5}{3}\right)^{n-1}} \Rightarrow \left(\frac{5}{3}\right)^n = \frac{5}{3}\frac{x_{n-1}-3}{5-x_{n-1}}.$$ and so $$x_n = 5 - \frac{2}{1 + \frac{5}{3}\frac{x_{n-1}-3}{5-x_{n-1}}} = 8 - \frac{15}{x_{n-1}}.$$

How does the floating point implementation of the original recurrence get it so wrong? To see this, note that in the case of a small rounding error, the closed form expression for $x_n$ becomes $$x_n \approx \frac{3^{n+1}+5^{n+1}+\gamma_n 100^{n+1}}{3^n + 5^n + \gamma_n 100^n}$$ where $\gamma_n$ is a small number on the order of the machine epsilon. The numerator may be rewritten \begin{align*} 3^{n+1}+5^{n+1}+\gamma_n 100^{n+1} &= (100-97)3^n + (100-95)5^n + 100\gamma_n 100^n\\ &= 100(3^n + 5^n + \gamma_n 100^n) - 97\cdot 3^n - 95\cdot 5^n. \end{align*} This leads to an alternative expression for $x_n$: $$x_n = 100 - \frac{97\cdot 3^n + 95\cdot 5^n}{3^n + 5^n + \gamma_n 100^n} = 100 - \frac{95 + 97\left(\frac{3}{5}\right)^n}{20^n\gamma_n + 1 + \left(\frac{3}{5}\right)^n}.$$ If $\gamma_n=0$, the correct limit $\lim_{n\to \infty}x_n = 5$ is obtained. However, if $\gamma_n \ne 0$ the growing term $20^n\gamma_n$ in the denominator leads to the change in the limit to 100.

References

1. Jean-Michel Muller's website.
2. W. Kahan, 2006, "How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?" [pdf].
3. This post was inspired by this math.stackexchange post which I happened upon.
Current rating: 4.3