05-basic-numerical-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Basic Numerical Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.1.1.** The standard Python library for linear algebra is [`numpy`](http://www.numpy.org/). In this section, we cover just enough of the numpy API to implement a few algorithms in the following sections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.1.2.** The basic object in `numpy` is `ndarray`, an $n$-dimensional generalization of Python's list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4]\n",
      "shape of a: (4,)\n",
      "\n",
      "[[1 2 3 4]]\n",
      "shape of b: (1, 4)\n",
      "\n",
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n",
      "shape of c: (4, 1)\n",
      "\n",
      "[[1 2]\n",
      " [3 4]]\n",
      "shape of d: (2, 2)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "a = np.array([1, 2, 3, 4])\n",
    "b = np.array([[1, 2, 3, 4]])\n",
    "c = np.array([[1], [2], [3], [4]])\n",
    "d = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "print(a)\n",
    "print('shape of a: {}'.format(a.shape))\n",
    "print()\n",
    "\n",
    "print(b)\n",
    "print('shape of b: {}'.format(b.shape))\n",
    "print()\n",
    "\n",
    "print(c)\n",
    "print('shape of c: {}'.format(c.shape))\n",
    "print()\n",
    "\n",
    "print(d)\n",
    "print('shape of d: {}'.format(d.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **shape** of an `ndarray` gives us the dimensions. `b` is a 1-by-4 matrix, or a row vector. `c` is a 2-by-2 vector, or a column vector. `d` is a 2-by-2 matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.1.3.** Entering a column vector requires many brackets, which can be tedious. We can instead use `transpose`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3 4]]\n",
      "shape of b: (1, 4)\n",
      "\n",
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n",
      "shape of b.transpose(): (4, 1)\n"
     ]
    }
   ],
   "source": [
    "print(b)\n",
    "print('shape of b: {}'.format(b.shape))\n",
    "print()\n",
    "\n",
    "print(b.transpose())\n",
    "print('shape of b.transpose(): {}'.format(b.transpose().shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, a matrix can be entered by first typing out a one-dimensional array and then putting the array through `reshape`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3 4]]\n",
      "shape of b: (1, 4)\n",
      "\n",
      "[[1 2]\n",
      " [3 4]]\n",
      "shape of b.reshape((2,2)): [[1 2]\n",
      " [3 4]]\n",
      "\n",
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n",
      "shape of b.reshape((4,1)): [[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n"
     ]
    }
   ],
   "source": [
    "print(b)\n",
    "print('shape of b: {}'.format(b.shape))\n",
    "print()\n",
    "\n",
    "print(b.reshape((2,2)))\n",
    "print('shape of b.reshape((2,2)): {}'.format(b.transpose().reshape((2,2))))\n",
    "print()\n",
    "\n",
    "print(b.reshape((4,1)))\n",
    "print('shape of b.reshape((4,1)): {}'.format(b.transpose().reshape((4,1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.1.4.** There are a number of pre-built functions for special kinds of vectors and matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4]\n",
      "\n",
      "[2 3 4 5 6 7]\n",
      "\n",
      "[ 2  5  8 11 14]\n",
      "\n",
      "[[1.]]\n",
      "\n",
      "[[1. 0.]\n",
      " [0. 1.]]\n",
      "\n",
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n",
      "\n",
      "[0.]\n",
      "\n",
      "[0. 0.]\n",
      "\n",
      "[0. 0. 0.]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(np.arange(5))\n",
    "print()\n",
    "\n",
    "print(np.arange(2, 8))\n",
    "print()\n",
    "\n",
    "print(np.arange(2, 15, 3))\n",
    "print()\n",
    "\n",
    "print(np.eye(1))\n",
    "print()\n",
    "\n",
    "print(np.eye(2))\n",
    "print()\n",
    "\n",
    "print(np.eye(3))\n",
    "print()\n",
    "\n",
    "print(np.zeros(1))\n",
    "print()\n",
    "\n",
    "print(np.zeros(2))\n",
    "print()\n",
    "\n",
    "print(np.zeros(3))\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.1.5.** Matrix multiplications are performed via `dot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1  1]\n",
      " [-2  2]]\n",
      "\n",
      "[[-1  1]\n",
      " [-2  2]]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([[2, 3], [5, 7]])\n",
    "y = np.array([[1, -1], [-1, 1]])\n",
    "\n",
    "print(np.dot(x,y))\n",
    "print()\n",
    "print(x.dot(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are running Python 3.5 or higher, the binary operator `@` may be used to denote matrix multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1  1]\n",
      " [-2  2]]\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "\n",
    "version_major, version_minor = sys.version_info[0:2]\n",
    "\n",
    "if version_major >= 3 and version_minor >= 5:\n",
    "    print(x @ y)\n",
    "else:\n",
    "    print('unsupported operation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.1.6.** `ndarray` supports coordinatewise operations. Observe, in particular, that `*` does not result in matrix multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4 6 8]\n",
      "\n",
      "[[3 2]\n",
      " [1 0]]\n",
      "\n",
      "[[ 1  4  3]\n",
      " [12 -5 -6]]\n",
      "\n",
      "[[0.5]\n",
      " [1.5]]\n"
     ]
    }
   ],
   "source": [
    "print(np.array([1, 2, 3]) + np.array([3,4,5]))\n",
    "print()\n",
    "\n",
    "print(np.array([[4,3],[2,1]]) - np.array([[1,1],[1,1]]))\n",
    "print()\n",
    "\n",
    "print(np.array([[1, 2, 3], [4, 5, 6]]) * np.array([[1, 2, 1], [3, -1, -1]]))\n",
    "print()\n",
    "\n",
    "print(np.array([[1], [3]]) / np.array([[2], [2]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operations with mismatching dimensions sometimes result in legitimate operation, via [broadcasting](https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html). For example, scalar multiplication works fine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6, 12,  9],\n",
       "       [ 3,  6, 15],\n",
       "       [-3, -3, -3]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3 * np.array([[2,4,3], [1,2,5], [-1, -1, -1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Row-wise operations and column-wise operations are also possible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 2]\n",
      " [3 4 5]\n",
      " [6 7 8]]\n",
      "\n",
      "[[ 5  6  7]\n",
      " [ 2  3  4]\n",
      " [ 9 10 11]]\n",
      "\n",
      "[[ 5  0  5]\n",
      " [ 8  3  8]\n",
      " [11  6 11]]\n"
     ]
    }
   ],
   "source": [
    "x = np.array([5, -1, 3])\n",
    "y = np.arange(9).reshape((3, 3))\n",
    "\n",
    "print(y)\n",
    "print()\n",
    "\n",
    "print(x.reshape((3, 1)) + y)\n",
    "print()\n",
    "\n",
    "print(x.reshape((1, 3)) + y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. Floating-Point Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.2.1.** The real and complex number systems that we use in mathematics unfortunately cannot be stored in a computer. On a computer, we use a finite approximation known as the floating-point number system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **floating-point number system** in base $\\beta$ with precision $p$, minimum exponent $e_\\min$, and maximum exponent $e_\\max$ is the discrete subset $\\mathscr{F} = \\mathscr{F}(\\beta,p,e_\\min,e_\\max)$ consisting of 0, $\\infty$, $-\\infty$, `NaN`, and numbers of the form\n",
    "\n",
    "$$\\pm(d_0 + d_1 / \\beta + \\cdots + d_{p-1} / \\beta^{p-1}) \\beta^e,$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "- $d_0,\\ldots,d_{p-1},e$ are integers,\n",
    "- $1 \\leq d_0 < \\beta$,\n",
    "- $0 \\leq d_i < \\beta$ for all $1 \\leq i \\leq p-1$, and\n",
    "- $e_\\min \\leq e \\leq e_\\max$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a floating-point number of the above form, we call the number $d_0 + d_1/\\beta + \\cdots + d_{p-1}/\\beta^{p-1}$ the **significand**, or the **mantissa**. The number $e$ is called the **exponent**. The sign $\\pm$ is called, naturally, the **sign**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The condition $d_0 \\geq 1$ makes the floating-point representation of a number unique. For this reason, a floating-point number with $d_0 \\geq 1$ is called **normalized**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.2.2.** How do we compute the floating-point representation $\\operatorname{fl}(x)$ of a real number $x$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It makes sense to define $\\operatorname{fl}(x)$ as the floating-point number *nearest* to $x$. To approximate $\\operatorname{fl}(x)$, we observe that $x$ admits a unique infinite-sum representation\n",
    "\n",
    "$$x = \\pm (d_0 + d_1/\\beta + \\cdots  + d_{p-1}/\\beta_{p-1} + d_p /\\beta_p + \\cdots) \\beta^e,$$\n",
    "\n",
    "provided that $d_0 \\geq 1$ and $d_i \\neq \\beta - 1$ for at least one $i \\geq 1$. The infinite-sum representation of $x$ shows that $\\operatorname{fl}(x)$ must lie between\n",
    "\n",
    "$$\\pm (d_0 + d_1/\\beta + \\cdots  + d_{p-1}/\\beta_{p-1}) \\beta^e - \\beta_{p-1}\\beta^e$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\pm (d_0 + d_1/\\beta + \\cdots  + d_{p-1}/\\beta_{p-1}) \\beta^e + \\beta_{p-1}\\beta^e$$\n",
    "\n",
    "It follows that\n",
    "\n",
    "$$\\vert\\operatorname{fl}(x) - x\\vert \\leq 2\\beta^{-(p-1)}\\beta^e.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain an error bound independent of the choice of $x$, we normalize the difference:\n",
    "\n",
    "$$\\delta_x = \\frac{\\vert \\operatorname{fl}(x) - x\\vert}{\\vert x \\vert}.$$\n",
    "\n",
    "Since $d_0$ in the infinite-series representation of $x$ must at least be 1, it follows that $\\vert x \\vert \\geq \\beta^e$. We now see that\n",
    "\n",
    "$$\\delta_x \\leq \\frac{2 \\beta^{-(p-1)\\beta^e}}{\\beta^e} = 2 \\beta^{-(p-1)}$$\n",
    "\n",
    "regardless of the choice of $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there exists at least one bound independent of the choice of $x$, we can find the smallest number $\\epsilon_{\\beta, p}$ such that\n",
    "\n",
    "$$\\delta_x \\leq \\epsilon_{\\beta, p}$$\n",
    "\n",
    "for every real number $x$. The bound $\\epsilon_{\\beta, p}$ is called the **machine epsilon** with respect to base $\\beta$ and precision $p$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of the optimal bound, many references settle for the bound\n",
    "\n",
    "$$\\delta_x \\leq \\frac{1}{2} \\beta^{-(p-1)},$$\n",
    "\n",
    "which is still four times better than the crude bound we have obtained above. It is common to see this quantity referred to see as the machine epsilon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.2.3.** We typically take the base $\\beta$ to be 2. For notational convenience, let us write\n",
    "\n",
    "$$(a_k \\ldots a_1.b_1 \\ldots b_m)_2$$\n",
    "\n",
    "to denote the sum\n",
    "\n",
    "$$\\sum_{i=1}^k a_i 2^i + \\sum_{j=1}^k b_i 2^{-j}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A typical binary representation of a floating-point number is\n",
    "\n",
    "$$\\underbrace{s}_{\\text{sign}} \\mid \\underbrace{e_1 \\cdots e_{q}}_{\\text{exponent}} \\mid \\underbrace{d_1 \\cdots d_{p-1}}_{\\text{significand}},$$\n",
    "\n",
    "which codifies\n",
    "\n",
    "$$\\pm(1.d_1 \\ldots d_{p-1})_2 \\times 2^{(e_1,\\ldots,e_q)_2 - (\\beta^{e_\\max} - 1)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is $d_0$ not included? Since $\\beta = 2$, normalization implies that $d_0$ always equals 1. There is no need to record it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of bits available for the exponent determines $e_\\min$ and $e_\\max$. Typically, [two's complement](https://en.wikipedia.org/wiki/Two%27s_complement) representation is used. If $e_1 = \\cdots = e_{p-1} = 1$, then the above binary representation codifies $\\infty$ if all significand bits are zero, and `NaN` otherwise. This practice exists to ensure consistent computational results across machines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.2.4.** The **IEEE Standard for Floating-Point Arithmetic** (IEEE 754) stipulates how a floating-point number system should be implemented in order to guarantee consistent computational results across machines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **single format** uses 1 bit for the sign, 8 bits for the exponent, and 23 bits for the significand, totalling 32 bits. The **double format** uses 1 bit for sign, 11 bits for the exponent, and 52 bits for the significand, totalling 64 bits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The floating-point representation of real numbers follows *rounding modes* specified in the standard. The most popular rounding mode is, of course, rounding to the nearest floating-point number (§**5.2.3**). The machine epsilon for the double format is approximately $2^{-52} \\approx 10^{-16}$. We can check this with a simple binary search:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.220446049250313e-16\n"
     ]
    }
   ],
   "source": [
    "e = 1.0\n",
    "while (1.0 + 0.5 * e) != 1.0:\n",
    "    e = 0.5 * e\n",
    "print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another important aspect of IEEE 754 is the precision requirements on floating-point arithmetic operations. Specifically, the standard stipulates that\n",
    "\n",
    "$$\\operatorname{fl}(x \\star y) = \\operatorname{fl}(x) \\star \\operatorname{f}(y)$$\n",
    "\n",
    "whenever $x$ and $y$ are real numbers and $\\star$ is addition, subtraction, multiplication, or division. From this stipulation follows the so-called **fundamental axiom of floting-point arithmetic**: *every operation of floating-point arithmetic is exact up to a relative error size of at most the machine epsilon.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.2.5.** The limitations of the floating-point arithmetic implies that we cannot expect a computer program to produce an exact answer to every problem. To quantify this, we let $X$ and $Y$ be vector spaces and define **mathematical problem** to be a function $f:X \\to Y$ that maps each *data point* $X$ to a *solution* in $Y$. The floating-point computed analogue of $f$ is denoted by $\\tilde{f}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the normalized error\n",
    "\n",
    "$$\\delta_x = \\frac{\\vert \\operatorname{fl}(x) - x\\vert}{\\vert x \\vert}$$\n",
    "\n",
    "for approximating a real number with a floating-point number is bounded by the machine epsilon $\\epsilon$ (§**5.2.2**). Since $\\tilde{f}$ can be thought of as a finite process of such approximations, we might hope for \n",
    "\n",
    "$$\\frac{\\|\\tilde{f}(x) - f(x)\\|}{\\|f(x)\\|} = O(\\epsilon)$$\n",
    "\n",
    "for all $x \\in X$, where $\\|\\cdot\\|$ is an appropriate norm defined on $Y$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is nearly always too strict, so we settle for\n",
    "\n",
    "$$\\frac{\\|\\tilde{f}(x) - f(\\tilde{x})\\|}{\\|f(\\tilde{x})\\|} = O(\\epsilon)$$\n",
    "\n",
    "for at least one $\\tilde{x}$, depending on the choice of $x$, such that\n",
    "\n",
    "$$\\frac{\\|\\tilde{x} - x\\|}{\\|x\\|} = O(\\epsilon).$$\n",
    "\n",
    "An algorithm that satisfies this condition is called **stable**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we allow for rounding errors in the input data, it is possible that, for all $x \\in X$, we have the equality\n",
    "\n",
    "$$\\|\\tilde{f}(x) = f(\\tilde{x})\\|$$\n",
    "\n",
    "for at least one $\\tilde{x}$ such that\n",
    "\n",
    "$$\\frac{\\|\\tilde{x} - x\\|}{\\|x\\|} = O(\\epsilon).$$\n",
    "\n",
    "Such an algorithm is called **backward stable**."
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}