{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 08 Errors\n",
"\n",
"(See also *Computational Physics* (Landau, Páez, Bordeianu), Chapter 3)\n",
"\n",
"\n",
"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. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Stupidity or Incompetence\n",
"(e.g., [PEBCAK](https://en.wiktionary.org/wiki/PEBCAK))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Random errors \n",
"- cosmic rays\n",
"- random bit flips"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Approximation errors\n",
"\n",
"\"**algorithmic errors**\"\n",
" - simplifying and adapting mathematics to the computer\n",
" - should decrease as $N$ increases\n",
"\n",
"#### Example:\n",
"Approximate $\\sin(x)$ with its truncated series expansion:\n",
"\\begin{align}\n",
"\\sin x &= \\sum_{n=1}^{+\\infty} \\frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!}\\\\\n",
" &\\approx \\sum_{n=1}^{N} \\frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!} + \\mathcal{E}(x, N)\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Round-off errors\n",
"- finite precision for storing floating-point numbers (32 bit, 64 bit)\n",
"- not known exactly (treat as uncertainty)\n",
"- can *accumulate* and lead to *garbage*\n",
"\n",
"#### Example: \n",
"Assume you can only store four decimals:\n",
"\n",
"\\begin{align}\n",
" \\text{storage}:&\\quad \\frac{1}{3} = 0.3333_c \\quad\\text{and}\\quad \\frac{2}{3} = 0.6667_c\\\\\n",
" \\text{exact}:&\\quad 2\\times\\frac{1}{3} - \\frac{2}{3} = 0\\\\\n",
" \\text{computer}:&\\quad 2 \\times 0.3333 - 0.6667 = -0.0001 \\neq 0\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"... now imagine adding \"$2\\times\\frac{1}{3} - \\frac{2}{3}$\" in a loop 100,000 times."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## The problems with *subtractive cancelation* \n",
"Model the computer representation $x_c$ of a number $x$ as\n",
"\n",
"$$\n",
"x_c \\simeq x(1+\\epsilon_x)\n",
"$$\n",
"\n",
"with the *relative* error $|\\epsilon_x| \\approx \\epsilon_m$ (similar to machine precision).\n",
"\n",
"Note: The *absolute* error is $\\Delta x = x_c - x$ and is related to the relative error by $\\epsilon_x = \\Delta x/x$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"What happens when we subtract two numbers $b$ and $c$: \n",
"\n",
"$$a = b - c$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"\\begin{gather}\n",
"a_c = b_c - c_c = b(1+\\epsilon_b) - c(1+\\epsilon_c)\\\\\n",
"\\frac{a_c}{a} = 1 + \\frac{b}{a}\\epsilon_b - \\frac{c}{a} \\epsilon_c\n",
"\\end{gather}\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"No guarantee that the errors cancel, and the relative error on $a$\n",
"\n",
"$$\n",
"\\epsilon_a = \\frac{a_c}{a} - 1 = \\frac{b}{a}\\epsilon_b - \\frac{c}{a} \\epsilon_c\n",
"$$ \n",
"\n",
"can be huge for small $a$!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Subtracting two nearly equal numbers\n",
"\n",
"$$b \\approx c$$ is bad!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"\\begin{align}\n",
"\\frac{a_c}{a} &= 1 + \\frac{b}{a}(\\epsilon_b - \\epsilon_c) \\\\\n",
"\\left| \\frac{a_c}{a} \\right| &\\leq 1 + \\left| \\frac{b}{a} \\right| (|\\epsilon_b| + |\\epsilon_a|)\n",
"\\end{align}\n",
"\n",
"i.e. the large number $b/a$ magnifies the error. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# Beware of subtractions!\n",
"\n",
"**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.**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Round-off errors\n",
"\n",
"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$):\n",
"$$\n",
"\\Delta f(x_1, x_2, \\dots; \\Delta x_1, \\Delta x_2, \\dots) =\n",
" \\sqrt{\\sum_{i=1}^N \\left(\\Delta x_i \\frac{\\partial f}{\\partial x_i}\\right)^2}\n",
"$$\n",
"\n",
"Note: relative error $$\\epsilon_i = \\frac{\\Delta x_i}{x_i}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Example: division $a = b/c$ (... with short cut)\n",
"\\begin{align}\n",
"a_c &= \\frac{b_c}{c_c} = \\frac{b(1+\\epsilon_b)}{c(1+\\epsilon_b)} \\\\\n",
"\\frac{a_c}{a} &= \\frac{1+\\epsilon_b}{1+\\epsilon_c} \n",
" = \\frac{(1+\\epsilon_b)(1-\\epsilon_c)}{1-\\epsilon_c^2} \\approx (1+\\epsilon_b)(1-\\epsilon_c)\\\\\n",
" &\\approx 1 + |\\epsilon_b| + |\\epsilon_c|\\\\\n",
"\\epsilon_a = \\frac{a_c}{a} - 1 &\\approx |\\epsilon_b| + |\\epsilon_c|\n",
"\\end{align}\n",
"\n",
"(neglected terms of order $\\mathcal{O}(\\epsilon^2)$); and same for multiplication.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**Errors accumulate with every operation.**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Model for round-off error accumulation\n",
"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"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"$$ R \\approx \\sqrt{N} r $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Total relative error $\\epsilon_{\\text{ro}}$ after $N$ calculations with error of the order of the machine precision $\\epsilon_m$ is\n",
"\n",
"$$ \\epsilon_{\\text{ro}} \\approx \\sqrt{N} \\epsilon_m $$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"(Only a model, depending on algorithm may be less or even $N!$...)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Total error of an algorithm\n",
"What you need to know to evaluate an algorithm:\n",
"1. Does it converge? (What $N$ do I need?)\n",
"2. How precise are the converged results (What is the error $\\epsilon_\\text{tot}$?)\n",
"3. What is its run time? (How fast is it for a given problem size?)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The total error contains *approximation* and *round off* errors:\n",
"\n",
"\\begin{gather}\n",
"\\epsilon_\\text{tot} = \\epsilon_\\text{app} + \\epsilon_\\text{ro}\n",
"\\end{gather}\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Model for the approximation error for an algorithm that takes $N$ steps (operations) to find a \"good\" answer:\n",
"\n",
"$$\n",
"\\epsilon_\\text{app} \\simeq \\frac{\\alpha}{N^\\beta}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"and round off error as\n",
"\n",
"$$\n",
"\\epsilon_{\\text{ro}} \\approx \\sqrt{N} \\epsilon_m\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Model for total error:\n",
"$$\n",
"\\epsilon_\\text{tot} = \\frac{\\alpha}{N^\\beta} + \\sqrt{N} \\epsilon_m\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Analyze $\\log_{10} $ of the relative error (direct readout of number of significant decimals).\n",
"\n",
"\n",
"\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Example analysis\n",
"\\begin{gather}\n",
"\\epsilon_\\text{app} = \\frac{1}{N^2}, \\quad \\epsilon_\\text{ro} = \\sqrt{N}\\epsilon_m\\\\\n",
"\\epsilon_\\text{tot} = \\frac{1}{N^2} + \\sqrt{N}\\epsilon_m\n",
"\\end{gather}"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Total error is a *minimum* for\n",
"\n",
"\\begin{gather}\n",
"\\frac{d\\epsilon_\\text{tot}}{dN} = -\\frac{2}{N^{3}} + \\frac{1}{2}\\frac{\\epsilon_m}{\\sqrt{N}} = 0, \\quad\\text{thus} \\quad\n",
"N^{5/2} = 4 \\epsilon_m^{-1}\\\\\n",
"N = \\left(\\frac{4}{\\epsilon_m}\\right)^{2/5}\n",
"\\end{gather}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"What is the best $N$ for single precision $\\epsilon_m \\approx 10^{-7}$?"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"import math\n",
"def N_opt(eps_m):\n",
" return round(math.pow(4./eps_m, 2./5.))\n",
"def eps_app(N):\n",
" return 1./(N*N)\n",
"def eps_ro(N, eps_m):\n",
" return math.sqrt(N)*eps_m"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best N = 1099 (for eps_m=1e-07)\n",
"eps_tot = 4.14e-06\n",
"eps_app = 8.28e-07, eps_ro = 3.32e-06\n"
]
}
],
"source": [
"epsilon_m = 1e-7 # single precision\n",
"\n",
"N = N_opt(epsilon_m)\n",
"err_app = eps_app(N)\n",
"err_ro = eps_ro(N, epsilon_m)\n",
"print(\"best N = {0} (for eps_m={1})\".format(N, epsilon_m))\n",
"print(\"eps_tot = {0:.3g}\".format(err_app + err_ro))\n",
"print(\"eps_app = {0:.3g}, eps_ro = {1:.3g}\".format(err_app, err_ro))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Single precision $\\epsilon_m \\approx 10^{-7}$:\n",
"\n",
"$$\n",
"N \\approx 1099\\\\\n",
"\\epsilon_\\text{tot} \\approx 4 \\times 10^{-6} \\\\\n",
"\\epsilon_\\text{app} = 8.28 \\times 10^{-7} \\\\\n",
"\\epsilon_\\text{ro} = 3.32 \\times 10^{-6}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Here, most of the error is round-off error! What can you do?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"* use double precision (delay round-off error)\n",
"* use a better algorithm, e.g. $\\epsilon_\\text{app}\\simeq \\frac{2}{N^4}$ (uses fewer steps)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**Better algorithms are always a good idea :-)**\n",
"\n",
"Remember: trade-off between **approximation error** and **rounding error**."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}