(See also Computational Physics (Landau, Páez, Bordeianu), Chapter 3)
These slides include material from Computational Physics. eTextBook Python 3rd Edition. Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License.
"algorithmic errors"
Approximate $\sin(x)$ with its truncated series expansion: \begin{align} \sin x &= \sum_{n=1}^{+\infty} \frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!}\\ &\approx \sum_{n=1}^{N} \frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!} + \mathcal{E}(x, N) \end{align}
Assume you can only store four decimals:
\begin{align} \text{storage}:&\quad \frac{1}{3} = 0.3333_c \quad\text{and}\quad \frac{2}{3} = 0.6667_c\\ \text{exact}:&\quad 2\times\frac{1}{3} - \frac{2}{3} = 0\\ \text{computer}:&\quad 2 \times 0.3333 - 0.6667 = -0.0001 \neq 0 \end{align}... now imagine adding "$2\times\frac{1}{3} - \frac{2}{3}$" in a loop 100,000 times.
Model the computer representation $x_c$ of a number $x$ as
$$ x_c \simeq x(1+\epsilon_x) $$with the relative error $|\epsilon_x| \approx \epsilon_m$ (similar to machine precision).
Note: The absolute error is $\Delta x = x_c - x$ and is related to the relative error by $\epsilon_x = \Delta x/x$.
What happens when we subtract two numbers $b$ and $c$:
$$a = b - c$$No guarantee that the errors cancel, and the relative error on $a$
$$ \epsilon_a = \frac{a_c}{a} - 1 = \frac{b}{a}\epsilon_b - \frac{c}{a} \epsilon_c $$can be huge for small $a$!
is bad!
i.e. the large number $b/a$ magnifies the error.
If you subtract two large numbers and end up with a small one, then the small one is less significant than any of the large ones.
Repeated calculations of quantities with errors beget new errors: In general, analyze with the rules of error propagation: function $f(x_1, x_2, \dots, x_N)$ with absolute errors on the $x_i$ of $\Delta x_i$ (i.e., $x_i \pm \Delta x_i$): $$ \Delta f(x_1, x_2, \dots; \Delta x_1, \Delta x_2, \dots) = \sqrt{\sum_{i=1}^N \left(\Delta x_i \frac{\partial f}{\partial x_i}\right)^2} $$
Note: relative error $$\epsilon_i = \frac{\Delta x_i}{x_i}$$
Example: division $a = b/c$ (... with short cut) \begin{align} a_c &= \frac{b_c}{c_c} = \frac{b(1+\epsilon_b)}{c(1+\epsilon_b)} \\ \frac{a_c}{a} &= \frac{1+\epsilon_b}{1+\epsilon_c} = \frac{(1+\epsilon_b)(1-\epsilon_c)}{1-\epsilon_c^2} \approx (1+\epsilon_b)(1-\epsilon_c)\\ &\approx 1 + |\epsilon_b| + |\epsilon_c|\\ \epsilon_a = \frac{a_c}{a} - 1 &\approx |\epsilon_b| + |\epsilon_c| \end{align}
(neglected terms of order $\mathcal{O}(\epsilon^2)$); and same for multiplication.
Errors accumulate with every operation.
View error in each calculation as a step in a random walk. The total "distance" (i.e. total error) $R$ over $N$ steps of length $r$ (the individual, "random" errors), is on average
Total relative error $\epsilon_{\text{ro}}$ after $N$ calculations with error of the order of the machine precision $\epsilon_m$ is
$$ \epsilon_{\text{ro}} \approx \sqrt{N} \epsilon_m $$(Only a model, depending on algorithm may be less or even $N!$...)
What you need to know to evaluate an algorithm:
The total error contains approximation and round off errors:
\begin{gather} \epsilon_\text{tot} = \epsilon_\text{app} + \epsilon_\text{ro} \end{gather}Model for the approximation error for an algorithm that takes $N$ steps (operations) to find a "good" answer:
$$ \epsilon_\text{app} \simeq \frac{\alpha}{N^\beta} $$and round off error as
$$ \epsilon_{\text{ro}} \approx \sqrt{N} \epsilon_m $$Model for total error: $$ \epsilon_\text{tot} = \frac{\alpha}{N^\beta} + \sqrt{N} \epsilon_m $$
Analyze $\log_{10} $ of the relative error (direct readout of number of significant decimals).
Image from Computational Physics. eTextBook Python 3rd Edition. Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License.
Total error is a minimum for
\begin{gather} \frac{d\epsilon_\text{tot}}{dN} = -\frac{2}{N^{3}} + \frac{1}{2}\frac{\epsilon_m}{\sqrt{N}} = 0, \quad\text{thus} \quad N^{5/2} = 4 \epsilon_m^{-1}\\ N = \left(\frac{4}{\epsilon_m}\right)^{2/5} \end{gather}What is the best $N$ for single precision $\epsilon_m \approx 10^{-7}$?
import math
def N_opt(eps_m):
return round(math.pow(4./eps_m, 2./5.))
def eps_app(N):
return 1./(N*N)
def eps_ro(N, eps_m):
return math.sqrt(N)*eps_m
epsilon_m = 1e-7 # single precision
N = N_opt(epsilon_m)
err_app = eps_app(N)
err_ro = eps_ro(N, epsilon_m)
print("best N = {0} (for eps_m={1})".format(N, epsilon_m))
print("eps_tot = {0:.3g}".format(err_app + err_ro))
print("eps_app = {0:.3g}, eps_ro = {1:.3g}".format(err_app, err_ro))
best N = 1099 (for eps_m=1e-07) eps_tot = 4.14e-06 eps_app = 8.28e-07, eps_ro = 3.32e-06
Single precision $\epsilon_m \approx 10^{-7}$:
$$ N \approx 1099\\ \epsilon_\text{tot} \approx 4 \times 10^{-6} \\ \epsilon_\text{app} = 8.28 \times 10^{-7} \\ \epsilon_\text{ro} = 3.32 \times 10^{-6} $$Here, most of the error is round-off error! What can you do?
Better algorithms are always a good idea :-)
Remember: trade-off between approximation error and rounding error.