Numerical stability
Contents
Numerical stability¶
In this notebook we look at some issues with the stability and errors that can arise when doing arithmetic on a computer. This will be a very quick pass over this rich topic, we encourage you to read the classic What Every Computer Scientist Should Know About Floating-Point Arithmetic to learn much more about this important topic.
In the last lecture, we discussed how binary floating point arithmetic can be subtle. Consider the follwing three numbers:
a, b, c = 1e-16, 1, -1e-16
Before running any code below, think about what the output of
should be.
Should it be any different whether you compute it as
or
Let’s take a look:
print(f"(a + b) + c = {(a + b) + c}")
print(f"a + (b + c) = {a + (b + c)}")
(a + b) + c = 0.9999999999999999
a + (b + c) = 1.0
Lesson: the finite precision of floating point arithmetic leads to the fact that even simple properties of elementary arithmetic, like the associativity of addition, don’t necessarily hold in all cases.
Now, the above issue is one related to loss of accuracy due to the finite precision and range of binary floating point. These issues run further, however. Consider computing
which obviously we expect to be equal to 0. Is it?
0.1 + 0.1 + 0.1 - 0.3
5.551115123125783e-17
0.1 + 0.2
0.30000000000000004
This problem isn’t one of loss of precision, but one of representability: due to the way binary floating point arithmetic represents numbers, not all decimal fractions can be correctly represented, and thus arithmetic operations incur small errors that can lead to unexpected results.
Decimal Arithmetic in Python¶
A quick aside before we continue. This kind of tool isn’t very commonly used in scientific computing, but the Python Standard Library includes the decimal module, an implementation of Decimal arithmetic (see the full IBM IBM’s General Decimal Arithmetic Specification for more details). Using Decimals, the above issue doesn’t occur:
from decimal import Decimal
a1, b1, c1 = Decimal(1e-16), Decimal(1), Decimal(-1e-16)
print(f"(a + b) + c = {(a1 + b1) + c1}")
print(f"a + (b + c) = {a1 + (b1 + c1)}")
(a + b) + c = 1.000000000000000000000000000
a + (b + c) = 1.000000000000000000000000000
And for our other example:
x, y = Decimal(0.1), Decimal(0.3)
x + x + x - y
Decimal('2.775557561565156540423631668E-17')
mmmhh… What happened here?
x, y = Decimal("0.1"), Decimal("0.3")
x + x + x - y
Decimal('0.0')
Back to numerics¶
Note: For numerical computing, if you need arbitrary precision arithmetic in floating point with support for advanced mathematical operations including special functions, the mpmath library is an excellent resource to know about.
This behavior can have serious implications in a variety of numerical work scenarios.
Consider the seemingly trivial problem of evaluating with a computer the expression
where \(r\) and \(x\) are real numbers with \(r \in [0,4]\) and \(x \in (0,1)\). This expression can also be written in an algebraically equivalent form:
Like above, when using binary floating point these two forms don’t necessarily produce the same answer. First a look at a few simple tests:
def f1(x): return r*x*(1-x)
def f2(x): return r*x - r*x**2
r = 1.9
x = 0.8
print('f1:', f1(x))
print('f2:', f2(x))
f1: 0.30399999999999994
f2: 0.3039999999999998
r = 3.9
x = 0.8
print('f1:', f1(x))
print('f2:', f2(x))
f1: 0.6239999999999999
f2: 0.6239999999999997
The difference is small but not zero:
print('difference:', (f1(x)-f2(x)))
difference: 2.220446049250313e-16
More importantly, this difference begins to accumulate as we perform the same operations over and over. Let’s illustrate this behavior by using the formulas above iteratively, that is, by feeding the result of the evaluation back into the same formula:
We can experiment with different values of \(r\) and different starting points \(x_0\) to observe the different results. We will simply build a python list that contains the results of three different (algebraically equivalent) forms of evaluating the above expression.
Exercise
Write code that computes the iteration of $f(x)$ using three different ways of writing the expression. Store your results and plot them using the `plt.plot()` function (the solution follows).
For completeness, we define three algebraically equivalent formulations:
def f1(x, r=r): return r*x*(1-x)
def f2(x, r=r): return r*x - r*x**2
def f3(x, r=r): return r*(x-x**2)
In order to see the difference between the initial behavior and the later evolution, let’s declare two variables to control our plotting:
num_points = 100 # total number of points to compute
drop_points = 0 # don't display the first drop_points
Using these, we can get the following figure:
%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
x0 = 0.55 # Any starting value
r = 3.9 # Change this to see changes in behavior
fp = (r-1.0)/r
x1 = x2 = x3 = x0
data = []
data.append([x1,x2,x3])
for i in range(num_points):
x1 = f1(x1)
x2 = f2(x2)
x3 = f3(x3)
data.append([x1,x2,x3])
# Display the results
plt.figure()
plt.title('r=%1.1f' % r)
plt.axhline(fp, color='black')
plt.plot(data[drop_points:], '-o', markersize=4)
plt.show()
Exercise
Now, experiment with different values of $r$ as well as different starting points $x_0$. What do you see? What happens when $r$ is small (close to 0)? Experiment with these values of $r$: \[1.9, 2.9, 3.1, 3.5, 3.9\] and think about the behavior of the system as you change $r$.
from ipywidgets import interact
fig, ax = plt.subplots()
@interact(r=(0.01, 4, 0.05))
def _(r=1):
x0 = 0.55 # Any starting value
fp = (r-1.0)/r
x1 = x2 = x3 = x0
data = []
data.append([x1,x2,x3])
for i in range(num_points):
x1 = f1(x1, r)
x2 = f2(x2, r)
x3 = f3(x3, r)
data.append([x1,x2,x3])
# Display the results
ax.clear()
ax.set_title('r=%1.1f' % r)
ax.axhline(fp, color='black')
ax.plot(data[drop_points:], '-o', markersize=4)
fig.canvas.draw_idle()
Once we’ve understood the basic pattern, let’s try to think of the entire evolution of the system as a function of \(r\). First, observe that a sequence generated by an iterative process of the form
will stop producing new values if there is a certain \(x^*\) such that
This special \(x^*\) is called a fixed point of the iterative process. It is easy to show that for our \(f(x)\), the fixed point is
(in fact, that’s the value plotted as a thin black line in the earlier script).
Exercise
Study whether the iteration converges to the fixed point or not by letting it run for each value of r for a few hundred points and discarding those, and then plotting the rest. Make a diagram with these plots as a function of r.
The following code is a simple solution:
def f(x, r):
return r*x*(1-x)
num_points = 500
drop_points = 50
rmin, rmax = 3.4, 4
xmin, xmax = 0, 1
x0 = 0.65
fig, ax = plt.subplots()
ax.set_xlim(rmin, rmax)
ax.set_ylim(xmin, xmax)
for r in np.linspace(rmin, rmax, 2000):
x = np.empty(num_points)
x[0] = x0
for n in range(1, num_points):
x[n] = f(x[n-1], r)
x = x[drop_points:]
rplot = r*np.ones_like(x)
ax.plot(rplot, x, 'b,')
plt.show()
Exercise
Can you relate the features of this figure to the behavior you saw in your earlier plots? Zoom in the region past $r=3$, what finer features do you see? Where is the fixed point we discussed earlier?
# Run this cell any time outputs have been removed to load CSS styles in use
from IPython.display import display, HTML
display(HTML(filename="../style.css"))