01-Scientific-Computing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scientific Computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 01.01 Scientific Computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given mathematical relationship $y = f(x)$, do one of the following:\n",
    "* evaluate\n",
    "    * given input $x$, then compute output $y$\n",
    "* solve\n",
    "    * find an input $x$ that produces output $y$\n",
    "* optimize\n",
    "    * find an input $x$ that produces an extreme value (max or min) $y$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to analyze a numerical solution:\n",
    "* discrete vs. continuous\n",
    "* linear vs. nonlinear\n",
    "* finite or infinite dimensional\n",
    "* purely algebraic or via derivatives or via integrals\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**General Approach** Replace a difficult problem by an easier one having the same or closely related solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 01.02 Approximations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sources of approximations.\n",
    "* Before computation.\n",
    "    * modelling\n",
    "    * measurements\n",
    "    * inherited from previous computations used as input\n",
    "* During computation.\n",
    "    * truncation / discretization\n",
    "    * rounding\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measuring error.\n",
    "* **absolute error** = approximate value - true value\n",
    "* **relative error** = absolute error / true value\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Types of error.\n",
    "\n",
    "$\n",
    "\\begin{aligned}\n",
    "\\text{total error} &= \\hat{f}(\\hat{x}) - f(x) \\\\\n",
    "&= \\text{computational error} + \\text{propagated data error} \\\\\n",
    "\\\\\n",
    "\\text{computational error} &= \\hat{f}(\\hat{x}) - f(\\hat{x}) \\\\\n",
    "\\text{propagated data error} &= f(\\hat{x}) - f(x) \\\\\n",
    "\\end{aligned}\n",
    "$\n",
    "\n",
    "where\n",
    "* $\\hat{f}(\\hat{x})$ approximate function with approximate input\n",
    "* $f(x)$ true function with true input\n",
    "* $f(\\hat{x})$ true function with approximate input\n",
    "\n",
    "\n",
    "The *computational error* is difference in result obtained using approximate function and true function.\n",
    "\n",
    "The *propagated data error* is difference in result obtained using approximate input and true input.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Types of computational error.\n",
    "\n",
    "* truncation error\n",
    "    * difference between true result and result using exact arithmetic\n",
    "        * source: mathematical approximations such as truncating series or discrete approximations\n",
    "* rounding error\n",
    "    * difference between result of approximate function using exact arithmetic and using limited precision\n",
    "        * source: inexact representation of real numbers and arithmetic\n",
    "\n",
    "In practice, usually **one** of these dominates.\n",
    "   * truncation error dominates in continuous problems\n",
    "   * rounding error dominates in algebraic problems\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 01.03 Forward and Backward Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Forward Error: $\\Delta{y} = \\hat{y} - y$\n",
    "\n",
    "Backward Error: $\\Delta{x} = \\hat{x} - x$\n",
    "\n",
    "where\n",
    "* $\\hat{y}$ is the computed result\n",
    "* $y$ is the true result\n",
    "* $\\hat{x}$ is the approximate input\n",
    "* $x$ is the true input\n",
    "\n",
    "In practice, **backward error** is easier to compute than **forward error**.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 01.04 Conditioning, Stability, and Accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A problem is **well-posed** if solution:\n",
    "* exists\n",
    "* unique\n",
    "* depends continuously on problem data\n",
    "\n",
    "A problem is **ill-conditioned** if relative change in solution can be much larger than input data.\n",
    "* The term **sensitivity** refers to how ill-conditioned a problem is.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Condition Number\n",
    "\n",
    "The condition number relates backward error to the forward error.\n",
    "\n",
    "$$\n",
    "|\\text{relative forward error}| = \\text{cond} \\times |\\text{relative backward error}|\n",
    "$$\n",
    "\n",
    "More formally.\n",
    "\n",
    "$\n",
    "\\begin{aligned}\n",
    "\\text{cond} &= \\frac{|[f(\\hat{x}) - f(x)] / f(x)|}{|(\\hat{x} - x) / x|} \\\\\n",
    "&= \\frac{|\\Delta{y} / y|}{|\\Delta{x} / x|} \\\\\n",
    "&= \\frac{|x f'(x)|}{|f(x)|}\n",
    "\\end{aligned}\n",
    "$\n",
    "\n",
    "where\n",
    "* *cond* is the condition number\n",
    "* $f(\\hat{x})$ is the true function with approximate input\n",
    "* $f(x)$ is the true function with true input\n",
    "* $\\hat{x}$ is the approximate input\n",
    "* $x$ is the true input\n",
    "\n",
    "A problem is **ill-conditioned** when $\\text{cond} >> 1$\n",
    "* The problem is **ill-conditioned** when the relative change in the output is much larger than input.\n",
    "\n",
    "Condition number of inverse of $f$ is reciprocal of the condition number of $f$.\n",
    "\n",
    "$\n",
    "\\text{cond}(f^{-1}) = \\frac{1}{\\text{cond}(f)}\n",
    "$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stability** is analogous to conditioning, but in context of *computational error*.\n",
    "* Computational error refers to effect on the result computed by an algorithm.\n",
    "* In contrast, conditioning refers to the effects of data error on solution to problem.\n",
    "\n",
    "In terms of error analysis, we say an algorithm is stable when the solution produced has relatively small backward error.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Accuracy** is the closeness of the computed solution to true solution and depends on:\n",
    "* conditioning of problem\n",
    "* stability of algorithm\n",
    "\n",
    "\n",
    "| Algorithm Stability | Problem Conditioning | Result | \n",
    "| ------------------- | -------------------- | ------ |\n",
    "| stable | well-conditioned | accurate |\n",
    "| stable | ill-conditioned | not accurate |\n",
    "| unstable | well-conditioned | not accurate |\n",
    "| unstable | ill-conditioned | not accurate |\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 01.05 Floating-Point Numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Representation\n",
    "\n",
    "Components of a floating point number system $\\mathbb{F}$:\n",
    "* $\\beta$ base or radix\n",
    "* $p$ precision\n",
    "* $[L, U]$ lower and upper exponent\n",
    "\n",
    "The floating point number $x \\in \\mathbb{F}$ has the form:\n",
    "\n",
    "$\n",
    "x = \\pm \\left( d_0 + \\frac{d_1}{\\beta} + \\frac{d_2}{\\beta^2} + ... + \\frac{d_{p-1}}{\\beta^{p-1}} \\right) \\beta^{E}\n",
    "$\n",
    "\n",
    "where\n",
    "* $d_i$ is an integer $0 \\leq d_i \\leq \\beta - 1, \\qquad i=0,...,p-1$\n",
    "* $E$ is an integer $L \\leq E \\leq U$\n",
    "\n",
    "The **mantissa**, m, refers to the parenthesized expression:\n",
    "\n",
    "$\n",
    "m = \\displaystyle \\sum_{i=0}^{p - 1} \\left( \\frac{d_i}{\\beta^i} \\right)\n",
    "$\n",
    "\n",
    "The mantissa is **normalized** when $1 \\leq m \\lt \\beta$.  This provides the following benefits:\n",
    "* Makes each number unique.\n",
    "* Eliminates any leading zeros, maximizing available precision. \n",
    "\n",
    "The **exponent** field, $E = U - L + 1$, determines range of representable magnitudes.\n",
    "\n",
    "* Number of normalized floating point numbers, N, in a system:\n",
    "\n",
    "$\n",
    "N = 2 (\\beta - 1) \\beta^{p -1} (U - L + 1) +1\n",
    "$\n",
    "\n",
    "* Underflow, UFL, smallest possible normalized number:\n",
    "\n",
    "$\n",
    "UFL = \\beta^L\n",
    "$\n",
    "\n",
    "* Overflow, OFL, largest possible normalized number:\n",
    "\n",
    "$\n",
    "OFL = \\beta^{U+1} (1 - \\beta^{-p})\n",
    "$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The components of the IEEE standard are listed in the table below:\n",
    "\n",
    "\n",
    "| System | $\\beta$ | $p$ | L | U |\n",
    "| ------ | ------- | ------ | - | - |\n",
    "| IEEE SP | 2 | 24 | -126 | 127 |\n",
    "| IEEE DP | 2 | 53 | -1022 | 1023 |\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the number of numbers for single and double precision:\n",
    "* single precision ~ $10^{9}$ aka giga\n",
    "* double precision ~ $10^{18}$ aka exa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the number-of-numbers for a given representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "61,441\n",
      "4,261,412,865\n",
      "18,428,729,675,200,069,633\n"
     ]
    }
   ],
   "source": [
    "def number_of_numbers(beta, p, L, U):\n",
    "    \"\"\"\n",
    "    Return the number of numbers in the given floating point system.\n",
    "    \"\"\"\n",
    "    return 2 * (beta-1) * beta**(p-1) * (U-L+1) + 1\n",
    "\n",
    "\n",
    "ieee_half_precision = {\n",
    "    'beta': 2,\n",
    "    'p': 11,\n",
    "    'L': -14,\n",
    "    'U': 15\n",
    "}\n",
    "\n",
    "ieee_single_precision = {\n",
    "    'beta': 2,\n",
    "    'p': 24,\n",
    "    'L': -126,\n",
    "    'U': 127,\n",
    "}\n",
    "\n",
    "ieee_double_precision = {\n",
    "    'beta': 2,\n",
    "    'p': 53,\n",
    "    'L': -1022,\n",
    "    'U': 1023,\n",
    "}\n",
    "\n",
    "print('{:,}'.format(number_of_numbers(**ieee_half_precision)))\n",
    "print('{:,}'.format(number_of_numbers(**ieee_single_precision)))\n",
    "print('{:,}'.format(number_of_numbers(**ieee_double_precision)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert a binary string to (single precision) floating point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25.033218383789062\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "def binary_to_float(binary, beta=2, p=24, L=-126, U=127):\n",
    "    \"\"\"\n",
    "    Return the floating point number from a binary representation.\n",
    "    \"\"\"\n",
    "    # Parse sign bit.\n",
    "    sign = int(binary[0], base=beta)\n",
    "    # Parse exponent as unsigned and rescale by subtracting L.\n",
    "    exp_bits = math.ceil(math.log(U-L+1, beta))\n",
    "    exp = int(binary[1:1+exp_bits], base=beta) + (L-1)\n",
    "    # Parse fraction aka mantissa and add implicit bit.\n",
    "    mantissa = '1' + binary[1+exp_bits:1+exp_bits+p]\n",
    "    sig, denom = 0., 1.\n",
    "    for d in mantissa:\n",
    "        if d == '1':\n",
    "            sig += denom\n",
    "        denom /= 2.\n",
    "    return sig * beta**(exp)\n",
    "\n",
    "binary_str = '01000001110010000100010000001000'\n",
    "print(binary_to_float(binary_str))\n",
    "\n",
    "# Compare to library method for converting binary to single precision.\n",
    "import struct\n",
    "expected = struct.unpack('f', struct.pack('I', int(binary_str,2)))[0]\n",
    "assert(expected == binary_to_float(binary_str))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the machine epsilon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.220446049250313e-16 53\n"
     ]
    }
   ],
   "source": [
    "def epsilon():\n",
    "    \"\"\"\n",
    "    Compute the machine epsilon and precision.\n",
    "    \"\"\"\n",
    "    eps, p = 1.0, 0\n",
    "    while (1 + eps) > 1:\n",
    "        eps /= 2.\n",
    "        p += 1\n",
    "    return eps*2., p  # Rescale epsilon to last true value.\n",
    "\n",
    "eps, p = epsilon()\n",
    "print(eps, p)\n",
    "\n",
    "# Compare to library constants.\n",
    "import sys\n",
    "assert(eps == sys.float_info.epsilon)\n",
    "assert(p == sys.float_info.mant_dig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exceptional Values\n",
    "\n",
    "**Inf** divide any finite number by zero eg 1/0\n",
    "\n",
    "**NaN** undefined operation eg 0/0, $0 \\times \\text{Inf}$, $\\text{Inf} / \\text{Inf}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 01.06 Floating-Point Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How is floating-point arithmetic performed?\n",
    "* addition\n",
    "  * shift mantissa until exponents match\n",
    "  * possible loss of digits of smaller number\n",
    "* multiplication\n",
    "  * product of 2 p-digit mantissas\n",
    "  * possible loss of digits if $p_i + p_j$ > machine precision\n",
    "* division\n",
    "  * quotient of 2 p-digit mantissas\n",
    "  * possible loss of digits if $\\frac{p_i}{p_j}$ > machine precision\n",
    "\n",
    "In general, **overflow** is worse than **underflow**.\n",
    "  * Overflow: No good approximations to arbitrarily large magnitudes.\n",
    "  * Underflow: Zero is reasonable approximation to small magnitudes.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Floating point arithmetic is **not** associative.\n",
    "\n",
    "$\n",
    "(1 + \\epsilon) + \\epsilon = 1 \\\\\n",
    "1 + (\\epsilon + \\epsilon) > 1\n",
    "$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cancellation** is result of subtracting numbers of similar magnitudes.\n",
    "  * The most significant aka leading digits of the results are lost.\n",
    "  * Compare to **rounding** where least significant aka trailing digits are lost.\n",
    "\n",
    "Demonstrations.\n",
    "\n",
    "* Example 1.\n",
    "\n",
    "Subtract two numbers which differ by $\\epsilon$, answer is $2 \\epsilon$ in real arithmetic.\n",
    "\n",
    "$\n",
    "(1 + \\epsilon) - (1 - \\epsilon) = 1 - 1 = 0\n",
    "$\n",
    "\n",
    "* Example 2.\n",
    "\n",
    "Summing alternating series when $x < 0$.\n",
    "\n",
    "$\n",
    "e^x = 1 + x + \\frac{x^2}{2!} + ...\n",
    "$\n",
    "\n",
    "\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}