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
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.
Comments
Comments are pre-moderated. Please be patient and your comment will appear soon.
There are currently no comments
New Comment