04_Numerical_Linear_Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>\n",
       "code_show=true; \n",
       "function code_toggle() {\n",
       " if (code_show){\n",
       " $('div.input').hide();\n",
       " } else {\n",
       " $('div.input').show();\n",
       " }\n",
       " code_show = !code_show\n",
       "} \n",
       "$( document ).ready(code_toggle);\n",
       "</script>\n",
       "The raw code for this Jupyter notebook is by default hidden for easier reading.\n",
       "To toggle on/off the raw code, click <a href=\"javascript:code_toggle()\">here</a>."
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "\n",
    "HTML('''<script>\n",
    "code_show=true; \n",
    "function code_toggle() {\n",
    " if (code_show){\n",
    " $('div.input').hide();\n",
    " } else {\n",
    " $('div.input').show();\n",
    " }\n",
    " code_show = !code_show\n",
    "} \n",
    "$( document ).ready(code_toggle);\n",
    "</script>\n",
    "The raw code for this Jupyter notebook is by default hidden for easier reading.\n",
    "To toggle on/off the raw code, click <a href=\"javascript:code_toggle()\">here</a>.''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](figs/COMPGEOP_TITLE.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numeric Linear Algebra\n",
    "This chapter discuss how to solve linear systems of equations numerically. We start with Gauss elimination, which are similar to how human being solving linear equations on the paper. We will also discuss several variants of this method (Doolittle, Cholesky, Gauss-Jordan). Then we will introduce numeric iteration methods to address the problem. These methods are less intuitive, but more computational effective and much more frequently used nowadays. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss Elimination\n",
    "Gauss elimination and back substitution are common methods to solve linear equation set using matrix method. A **linear system of $n$ equations** in $n$ unknowns $x_1,...,x_n$ is a set of equations of the form\n",
    "$$\n",
    "a_{11}x_1 + ... a_{1n}x_n = b_1 \\\\\n",
    "a_{21}x_1 + ... a_{2n}x_n = b_2 \\\\\n",
    "... .........\\\\\n",
    "a_{n1}x_1 + ... a_{nn}x_n = b_n\n",
    "$$\n",
    "where the **coefficient** $a_{jk}$ and the $b_j$ are given numbers. This equation set can be written in matrix form as\n",
    "$$ Ax=b $$\n",
    "where the **coefficient matrix** $A=[a_{jk}]$ is the $n \\times n$ matrix, $x$ and $b$ are $n\\times 1$ vectors.\n",
    "$$\n",
    "A = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a{12} & ... & a_{1n} \\\\\n",
    "a_{21} & a{22} & ... & a_{2n} \\\\\n",
    ". & . & . & . \\\\\n",
    "a_{n1} & a{n2} & ... & a_{nn}\n",
    "\\end{matrix} \\right], \\ \\ x=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n\n",
    "\\end{matrix}\n",
    "\\right], \\ \\ b=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "b_1 \\\\ b_2 \\\\ \\vdots \\\\ b_n\n",
    "\\end{matrix}\n",
    "\\right].\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then form an **argmented matrix** of the system:\n",
    "$$\n",
    "\\tilde{A} = [A\\ b]=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & ... & a_{1n} & b_1 \\\\\n",
    "a_{21} & ... & a_{2n} & b_2 \\\\\n",
    ".&.&.&.\\\\\n",
    "a_{n1} & ... & a_{nn} & b_n \\\\\n",
    "\\end{matrix}\n",
    "\\right].\n",
    "$$\n",
    "There are **three elementary row operations** we can perform on $\\tilde A$ that will not change the solution of the linear equation system:\n",
    "1. Interchange of two rows\n",
    "2. Addition of a constant multiple of one row to another row\n",
    "3. Multiplication of a row by a nonzero constant $c$\n",
    "\n",
    "These three operations correspond to the following \n",
    "1. Interchange of two equations\n",
    "2. Addition of a constant multiple of one equation to another euqation\n",
    "3. Multiplication of an equation by a nonzero constant $c$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gauss Elimination try to reduce the coefficient matrix $A$ from a full matrix to a **triangular** matrix through a combination of these three elementary row operations, which can then be easily solved by **back substitution**. For instance, a triangular system is\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat a_{11} x_1 + \\hat a_{12} x_2 + \\hat a_{13} x_3 + ...+ \\hat a_{1n} x_n &= \\hat b_1 \\\\\n",
    "\\hat a_{22} x_2 + \\hat a_{23} x_3 + ...+ \\hat a_{2n} x_n &=\\hat b_2 \\\\\n",
    "\\hat a_{33} x_3 + ... + \\hat a_{3n} x_n &= \\hat b_3 \\\\\n",
    "...... \\\\\n",
    "\\hat a_{n-1 n-1} x_{n-1} + \\hat a_{n-1n} x_n &=\\hat b_{n-1} \\\\\n",
    "\\hat a_{nn} x_n &= \\hat b_n\n",
    "\\end{align}\n",
    "$$\n",
    "Then from the last equation, we can solve the entire system step by step by **back substitution**:\n",
    "$$\n",
    "\\begin{align}\n",
    "x_n &= \\frac{\\hat b_n}{\\hat a_{nn}} \\\\\n",
    "x_{n-1} &= \\frac{1}{\\hat a_{n-1n-1}}(\\hat b_{n-1} - \\hat a_{n-1n}x_n) \\\\\n",
    "& ...... \\\\\n",
    " x_1  &=\\frac{1}{\\hat a_{11}} (\\hat b_1 - \\hat a_{12} x_2 - \\hat a_{13} x_3 - ...- \\hat a_{1n} x_n)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To transform coefficient matrix $A$ into a triangular form, we first eliminate coefficient of $x_1$ from row 2 to $n$ of $\\tilde A$ by subtracting suitable multiples of row 1 from the other rows. In this step, the first row is called the **pivot equation** and $a_{11}$ is called the **pivot**. We changes all the rows of $\\tilde A$ except for the pivot row. We then eliminate coefficient of $x_2$ from row 3 to $n$ by subtracting suitable multiples of row 2 from all the rows below, and so on. \n",
    "$$\n",
    "\\tilde{A_0} = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} & ... & a_{1n} & b_1 \\\\\n",
    "a_{21} & a_{22} &... & a_{2n} & b_2 \\\\\n",
    ".&.&.&.&.\\\\\n",
    "a_{n1} & a_{n2} &... & a_{nn} & b_n \n",
    "\\end{matrix}\n",
    "\\right].\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\tilde{A_1} = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} & ... & a_{1n} & b_1 \\\\\n",
    "a_{21} - a_{11}\\frac{a_{21}}{a_{11}} & a_{22}- a_{12}\\frac{a_{21}}{a_{11}} & ... \n",
    "& a_{2n} - a_{1n}\\frac{a_{21}}{a_{11}} & b_2 - b_1\\frac{a_{21}}{a_{11}} \\\\\n",
    ".&.&.&.&.\\\\\n",
    "a_{n1} - a_{11}\\frac{a_{n1}}{a_{11}} & a_{n2}- a_{12}\\frac{a_{n1}}{a_{11}} & ... \n",
    "& a_{nn} - a_{1n}\\frac{a_{n1}}{a_{11}} & b_n - b_1\\frac{a_{n1}}{a_{11}} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} & ... & a_{1n} & b_1 \\\\\n",
    "0 & a_{22}^{(1)} &... & a_{2n}^{(1)} & b_2^{(1)} \\\\\n",
    ".&.&.&.&.\\\\\n",
    "0 & a_{n2}^{(1)} &... & a_{nn}^{(1)} & b_n^{(1)} \n",
    "\\end{matrix}\n",
    "\\right]. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\vdots \n",
    "$$\n",
    "$$\n",
    "\\vdots \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde A_n = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} &a_{13} & ... & a_{1n} & b_1 \\\\\n",
    "0 & a_{22}^{(1)} &a_{23}^{(1)} &... & a_{2n}^{(1)} & b_2^{(1)} \\\\\n",
    "0 & 0 & a_{33}^{(2)} &... & a_{3n}^{(2)} & b_2^{(2)} \\\\\n",
    ".&.&.&.&.&.\\\\\n",
    "0 & 0 & 0 &... & a_{nn}^{(n-1)} & b_n^{(n-1)} \n",
    "\\end{matrix}\n",
    "\\right]. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One practical concern about Gauss Elimination is that the pivot $a_{kk}$ **must be** different from zero and **should be** large in absolute value to avoid numeric instability in the elimination. For this reason in each step we choose our pivot row the one with the absolutely largest $a_{jk}$ (with $j>k$) and interchange it with current row $k$. This method is called **partial pivoting** and is quite popular (e.g., Maple)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Gauss_elimination(A,b,print_process=False):\n",
    "    ''' Function that utilizes Gauss Elimination and back substitution to solve\n",
    "    a linear system Ax = b\n",
    "    usage: x = Gauss_elimination(A,b,print_process=False):\n",
    "    input: \n",
    "        A: nxn coefficient matrix, full rank\n",
    "        b: nx1 vector\n",
    "        print_process: boolean, whether to print the step-by-step process \n",
    "                        for Gauss Elimination\n",
    "    output:\n",
    "        solution x\n",
    "    written by Ge Jin, gjin@mines.edu, 06/2019\n",
    "    '''\n",
    "    n = A.shape[0]\n",
    "    # generate argmented matrix\n",
    "    Ab = np.hstack((A,b.reshape(-1,1)))\n",
    "    # Gauss Elimination\n",
    "    for k in range(n-1):\n",
    "        # find the row with max ajk\n",
    "        maxi = np.argmax(np.abs(Ab[k:,k]))\n",
    "        maxi += k\n",
    "        # swap the rows\n",
    "        Ab[[k,maxi]] = Ab[[maxi,k]]\n",
    "        # eliminate ajk from row k+1 to n\n",
    "        for j in range(k+1,n):\n",
    "            Ab[j,:] = Ab[j,:] - Ab[k,:]*Ab[j,k]/Ab[k,k]\n",
    "        if print_process:\n",
    "            print('Ab_{}='.format(k))\n",
    "            print(Ab)\n",
    "\n",
    "    # back substitution\n",
    "    x = np.zeros(n)\n",
    "    x[n-1] = Ab[n-1,n]/Ab[n-1,n-1]\n",
    "    for k in range(len(x)-2,-1,-1):\n",
    "        x[k] = 1/Ab[k,k]*(Ab[k,n]-np.sum(Ab[k,:-1]*x))\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.633556604385376"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from time import time\n",
    "tic = time()\n",
    "N = 1000\n",
    "x = Gauss_elimination(np.random.rand(N,N)+1,np.random.rand(N),print_process=False)\n",
    "toc =time()\n",
    "toc-tic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ab_0=\n",
      "[[ 6.  2.  8. 26.]\n",
      " [ 0.  4. -2. -5.]\n",
      " [ 0.  8.  2. -7.]]\n",
      "Ab_1=\n",
      "[[ 6.   2.   8.  26. ]\n",
      " [ 0.   8.   2.  -7. ]\n",
      " [ 0.   0.  -3.  -1.5]]\n",
      "x= [ 4.  -1.   0.5]\n"
     ]
    }
   ],
   "source": [
    "# an example\n",
    "A = np.array([[6.,2,8],[3,5,2],[0,8,2]])\n",
    "b = np.array([26.,8.,-7.])\n",
    "x = Gauss_elimination(A,b,print_process=True)\n",
    "print('x=',x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ab_0=\n",
      "[[ 3.  6.  3. 12.]\n",
      " [ 0. -1.  0. -1.]\n",
      " [ 0. -3. -1. -4.]]\n",
      "Ab_1=\n",
      "[[ 3.          6.          3.         12.        ]\n",
      " [ 0.         -3.         -1.         -4.        ]\n",
      " [ 0.          0.          0.33333333  0.33333333]]\n",
      "x= [1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "# an example\n",
    "A = np.array([[3.,6,3],[1,1,1],[2,1,1]])\n",
    "b = np.array([12.,3.,4.])\n",
    "x = Gauss_elimination(A,b,print_process=True)\n",
    "print('x=',x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computation cost of Gauss Elimination\n",
    "Generally, the important factors in judging the quality of a numeric method are\n",
    "- Amount of storage\n",
    "- Amount of time (number of operations)\n",
    "- Effect of roundoff error (numerical stability)\n",
    "\n",
    "For Gauss elimination, the operation count for a full matrix is as follows. In Step $k$ we eliminate $x_k$ from $n-k$ equations. This needs $n-k$ divisions in computing the $a_{jk}/a_{kk}$, and $(n-k)(n-k+1)$ multiplications and as many subtractions. Since we do $n-1$ steps, $k$ goes from 1 to $n-1$ and thus the total number of operations in this forward elimination is\n",
    "$$\n",
    "f(n) = \\sum_{k=1}^{n-1}(n-k) + 2\\sum_{k=1}^{n-1}(n-k)(n-k+1)\n",
    "= \\frac{1}{2}(n-1)n + \\frac{2}{3}(n^2-1)n \\approx \\frac{2}{3} n^3\n",
    "$$\n",
    "where $2n^3/3$ is obtained by dropping lower powers of $n$ (when $n$ is large, the only term that matters is the one with highest power of $n$). Because $f(n)$ grows about proportional to $n^3$, we say that $f(n)$ is order $n^3$ and write\n",
    "$$\n",
    "f(n) = O(n^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the back substituion of $x_i$ we make $n-i$ multiplications and as many subtractions, as well as 1 division. Hence the number of operations in the back substitution is\n",
    "$$\n",
    "b(n) = 2\\sum_{i=1}^{n}(n-i) + n = n(n+1)+n = N^2+2n = O(n^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For instance, if an operation takes $10^{-9}$ second, then the times needed are:\n",
    "\n",
    "| Algorithm | $n=1000$ | $n=10000$ |\n",
    "|---------|-------------------|----------------|\n",
    "|Elimination | 0.7 sec | 660 sec |\n",
    "|Back substitution | 0.001 sec | 0.1 sec |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LU-Factorization\n",
    "Another way to solve linear systems $ Ax=b$ numerically if through **LU-Factorization**. An LU-factorization of a given square matrix $A$ is of the form\n",
    "$$ A=LU $$\n",
    "where $L$ is **lower triangular** matrix and $U$ is **upper triangular** matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be proved that for any nonsingular matrix, the rows can be reordered so that the resulting matrix A has an LU-factorization in which $L$ turns out to be the matrix of the multipliers $m_{jk}$ of the Gauss elimination, with main diagonals equal 1, and $U$ is the matrix of the triangular system at the end of the Gauss elimination. \n",
    "$$\n",
    "A = LU = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "1 & 0 & 0 & ... & 0\\\\ \n",
    "m_{21} & 1 & 0  &...  &0 \\\\ \n",
    "m_{31} &m_{32}  &1  &...  &0 \\\\ \n",
    "\\vdots &\\vdots  &\\vdots  &\\vdots  &\\vdots \\\\ \n",
    "m_{n1} &m_{n2}  &m_{n3}  &...  &1 \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "u_{11} &u_{12}  &u_{13}  &...  &u_{1n} \\\\ \n",
    "0 &u_{22}  &u_{23}  &...  &u_{2n} \\\\ \n",
    "0 &0  &u_{33}  &...  &u_{3n} \\\\ \n",
    "\\vdots &\\vdots  &\\vdots  &\\vdots  &\\vdots \\\\ \n",
    "0 &0  &0  &...  &u_{nn} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And $L$ and $U$ can be calculated using **Doolittle method**:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "u_{1i} &= a_{1i}  &i &= 1,...,n \\\\\n",
    "m_{j1} &= \\frac{a_{j1}}{u_{11}}  &j &=2,...,n \\\\\n",
    "u_{jk} &= a_{jk} - \\sum_{s=1}^{j-1} m_{js}u_{sk}   &k &=j,...,n;\\;\\; j \\geq 2 \\\\\n",
    "m_{lj} &= \\frac{1}{u_{jj}} \n",
    "\\left(\n",
    "a_{lj} - \\sum_{s=1}^{j-1} m_{ls}u_{sj}\n",
    "\\right)   &l &=j+1,...,n; \\;\\; k\\geq 2\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The crucial idea is that $L$ and $U$ can be computed directly, without solving simultaneous equations (thus, without suing the Gauss elimiantion). As a count shows, this needs about $n^3/3$ operations, about half as many as the Gauss elimination. Once we have $L$ and $U$, we can use it for solving $Ax=b$ in two steps, involving only about $n^2$ operations.\n",
    "\n",
    "Noting that $Ax = LUx = b$ may be written as\n",
    "$$\n",
    "Ly = b \\; \\; \\; where \\; \\; \\; Ux=y\n",
    "$$\n",
    "We can use forward substitution to solve $y$, and then use backward substitution to solve $x$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A similar method, **Crout's method**, is obtained if $U$ (instead of $L$) is required to have main diagonals equal 1. In either case the factorization is unique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "1. Solve the LU factorization of a general $3 \\times 3$ matrix \n",
    "$$\n",
    "A = \\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} & a_{13} \\\\\n",
    "a_{21} & a_{22} & a_{23} \\\\\n",
    "a_{31} & a_{32} & a_{33} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cholesky's Method\n",
    "If matrix $A$ is symmetric and positive definite ($A=A^T,x^TAx > 0$ for all $x\\neq0$), we can even choose $U=L^T$ to further reduce the computational cost. In this case, we cannot impose conditions on the main diagonal entries. \n",
    "$$\n",
    "A = LL^T=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "l_{11} & 0 & 0 & ... & 0 \\\\\n",
    "l_{21} & l_{22} & 0 & ... & 0 \\\\\n",
    "l_{31} & l_{32} & l_{33} & ... & 0 \\\\\n",
    "\\vdots &\\vdots &\\vdots &\\vdots &\\vdots \\\\\n",
    "l_{n1} & l_{n2} & l_{n3} & ... & l_{nn} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "l_{11} & l_{21} & l_{31} & ... & l_{n1} \\\\\n",
    "0 & l_{22} & l_{32} & ... & l_{n2} \\\\\n",
    "0 & 0 & l_{33} & ... & l_{n3} \\\\\n",
    "\\vdots &\\vdots &\\vdots &\\vdots &\\vdots \\\\\n",
    "0 & 0 & 0 & ... & l_{nn} \n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The formulas to calculate $l_{jk}$ are:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "l_{11} &= \\sqrt{a_{11}}  &  & \\\\\n",
    "l_{j1} &= \\frac{a_{j1}}{l_{11}}  &j &=2,...,n \\\\\n",
    "l_{jj} &= \\sqrt{a_{jj} - \\sum_{s=1}^{j-1} l_{js}^2}   &j &=2,...,n \\\\\n",
    "l_{pk} &= \\frac{1}{l_{jj}} \n",
    "\\left(\n",
    "a_{pj} - \\sum_{s=1}^{j-1} l_{js}l_{ps}\n",
    "\\right)   &p &=j+1,...,n; \\;\\; j\\geq 2\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss-Jordan Elimination: Matrix Inversion\n",
    "Sometimes it is useful to get the inversion of the matrix $A$, and **Gauss-Jordan elimination** is one of the methods calculate it. The inversion of a matrix is defined as\n",
    "$$\n",
    "AA^{-1} = A^{-1}A = I\n",
    "$$\n",
    "To get $A^{-1}$, similar to Gauss elimination, we can form the **argmented matrix**:\n",
    "$$\n",
    "\\tilde{A} = [A\\;I]=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & ... & a_{1n} & 1 & 0 & 0 & ... & 0 \\\\\n",
    "a_{21} & ... & a_{2n} & 0 & 1 & 0 & ... & 0 \\\\\n",
    ".&.&.&.&.&.&.&.\\\\\n",
    "a_{n1} & ... & a_{nn} & 0 & 0 & 0 & ... & 1 \\\\\n",
    "\\end{matrix}\n",
    "\\right].\n",
    "$$\n",
    "We then apply the Gauss elimination to $\\tilde A$, which gives a matrix of the form $[U \\; H]$ with $U$ being an upper triangular matrix. The Gauss-Jordan method reduces U by further elementary row operations to diagonal form $[I\\; K]$, where $K$ is the inversion of $A$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gauss_jordan_matrix_inv(A,print_process=False):\n",
    "    ''' Function that utilizes Gauss-Jordan Elimination \n",
    "    to calculate the inversion of the input matrix\n",
    "    usage: x = gauss_jordan_matrix_inv(A,print_process=False):\n",
    "    input: \n",
    "        A: nxn coefficient matrix, full rank\n",
    "        print_process: boolean, whether to print the step-by-step process \n",
    "    output:\n",
    "        inversion of A\n",
    "    written by Ge Jin, gjin@mines.edu, 06/2019\n",
    "    '''\n",
    "    n = A.shape[0]\n",
    "    step = 0\n",
    "    # generate argmented matrix\n",
    "    Ab = np.hstack((A,np.identity(n)))\n",
    "    # Gauss Elimination\n",
    "    for k in range(n-1):\n",
    "        # find the row with max ajk\n",
    "        maxi = np.argmax(np.abs(Ab[k:,k]))\n",
    "        maxi += k\n",
    "        # swap the rows\n",
    "        Ab[[k,maxi]] = Ab[[maxi,k]]\n",
    "        # eliminate ajk from row k+1 to n\n",
    "        for j in range(k+1,n):\n",
    "            Ab[j,:] = Ab[j,:] - Ab[k,:]*Ab[j,k]/Ab[k,k]\n",
    "        if print_process:\n",
    "            print('Ab_{}='.format(step))\n",
    "            print(Ab)\n",
    "            step += 1\n",
    "    # change diagonal elements to 1 \n",
    "    for k in range(n):\n",
    "        Ab[k] /= Ab[k,k]\n",
    "        if print_process:\n",
    "            print('Ab_{}='.format(step))\n",
    "            print(Ab)\n",
    "            step += 1\n",
    "    # Gauss-Jordan elimination\n",
    "    for k in range(n-1,0,-1):\n",
    "        for i in range(0,k):\n",
    "            Ab[i] -= Ab[k]*Ab[i,k]\n",
    "            if print_process:\n",
    "                print('Ab_{}='.format(step))\n",
    "                print(Ab)\n",
    "                step += 1\n",
    "    invA = Ab[:,n:]\n",
    "    return invA\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ab_0=\n",
      "[[ 3.         -1.          1.          0.          1.          0.        ]\n",
      " [ 0.          0.66666667  2.33333333  1.          0.33333333  0.        ]\n",
      " [ 0.          2.66666667  4.33333333  0.          0.33333333  1.        ]]\n",
      "Ab_1=\n",
      "[[ 3.         -1.          1.          0.          1.          0.        ]\n",
      " [ 0.          2.66666667  4.33333333  0.          0.33333333  1.        ]\n",
      " [ 0.          0.          1.25        1.          0.25       -0.25      ]]\n",
      "Ab_2=\n",
      "[[ 1.         -0.33333333  0.33333333  0.          0.33333333  0.        ]\n",
      " [ 0.          2.66666667  4.33333333  0.          0.33333333  1.        ]\n",
      " [ 0.          0.          1.25        1.          0.25       -0.25      ]]\n",
      "Ab_3=\n",
      "[[ 1.         -0.33333333  0.33333333  0.          0.33333333  0.        ]\n",
      " [ 0.          1.          1.625       0.          0.125       0.375     ]\n",
      " [ 0.          0.          1.25        1.          0.25       -0.25      ]]\n",
      "Ab_4=\n",
      "[[ 1.         -0.33333333  0.33333333  0.          0.33333333  0.        ]\n",
      " [ 0.          1.          1.625       0.          0.125       0.375     ]\n",
      " [ 0.          0.          1.          0.8         0.2        -0.2       ]]\n",
      "Ab_5=\n",
      "[[ 1.         -0.33333333  0.         -0.26666667  0.26666667  0.06666667]\n",
      " [ 0.          1.          1.625       0.          0.125       0.375     ]\n",
      " [ 0.          0.          1.          0.8         0.2        -0.2       ]]\n",
      "Ab_6=\n",
      "[[ 1.         -0.33333333  0.         -0.26666667  0.26666667  0.06666667]\n",
      " [ 0.          1.          0.         -1.3        -0.2         0.7       ]\n",
      " [ 0.          0.          1.          0.8         0.2        -0.2       ]]\n",
      "Ab_7=\n",
      "[[ 1.   0.   0.  -0.7  0.2  0.3]\n",
      " [ 0.   1.   0.  -1.3 -0.2  0.7]\n",
      " [ 0.   0.   1.   0.8  0.2 -0.2]]\n",
      "Inverted Matrix:\n",
      "[[-0.7  0.2  0.3]\n",
      " [-1.3 -0.2  0.7]\n",
      " [ 0.8  0.2 -0.2]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[-1,1,2],[3,-1,1],[-1,3,4]])\n",
    "invA = gauss_jordan_matrix_inv(A,print_process=True)\n",
    "print('Inverted Matrix:')\n",
    "print(invA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Systems: Solution by Iteration\n",
    "The Gauss elimination and its variants we just discussed belong to the **direct methods** for solving linear systems of equations. These methods give solutions after an amount of computation that can specified in advance. In contrast, in an **indirect** or **iterative method**, we start from an approximation to the true solution and, if successful, obtain better and better approximations from a computational cycle repeated as often as may be necessary for achieving a required accuracy. \n",
    "Iterative methods in general have less computational cost comparing to direct methods, especially for very large matrix or very sparse ones. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss-Seidel Iteration Method\n",
    "This is an iterative method of great practical importance for solving linear system $Ax=b$. To obtain the algorithm, let us derive the general formulas for this iteration. \n",
    "\n",
    "We can rearrange the equations so that $a_{jj}>0$ for $j=1,...,n$. We then decompose $A$ into a lower triangular component $L_*$, and a strictly upper triangular component $U$\n",
    "$$\n",
    "A = L_*+U \\;\\; where \\;\\;\n",
    "L_* = \\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & 0 & ... & 0 \\\\\n",
    "a_{21} & a_{22} & ...& 0 \\\\\n",
    "\\vdots &\\vdots &\\ddots &\\vdots \\\\\n",
    "a_{n1} & a_{n2} & ...& a_{nn} \\\\\n",
    "\\end{matrix}\n",
    "\\right], \\;\\; \n",
    "U = \\left[\n",
    "\\begin{matrix}\n",
    "0 & a_{12} & ... & a_{1n} \\\\\n",
    "0 & 0 & ... & a_{2n} \\\\\n",
    "\\vdots &\\vdots &\\ddots &\\vdots \\\\\n",
    "0 & 0 & ... & 0 \n",
    "\\end{matrix}\n",
    "\\right], \\;\\; \n",
    "$$\n",
    "\n",
    "By substituting it into $Ax=b$, we have\n",
    "$$\n",
    "Ax = (L_*+U)x = b.\n",
    "$$\n",
    "$$ L_*x = b-Ux $$\n",
    "The Gauss-Seidel method now solves the left hand side of this expression for $x$, using previous value for $x$ on the right hand side. Analytically, this may be written as:\n",
    "$$ L_* x^{(k+1)} = b-Ux^{(k)}. $$\n",
    "By taking advantage of the triangular form of $L_*$, the elements of $x^{(k+1)}$ can be computed sequentially using [forward substitution](https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution):\n",
    "$$\n",
    "x_i^{k+1} = \\frac{1}{a_{ii}} \\left(\n",
    "b_i - \\sum_{j=1}^{i-1} a_{ij}x_j^{(k+1)} - \\sum_{j=i+1}^n a_{ij}x_j^{(k)}\n",
    "\\right),\\ \\ \\ \\ i=1,2,...,n\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Convergence**\n",
    "\n",
    "One of the common problem for all iteration-based methods is that the iteration may not always converge to the true answer. For Gauss-Seidel method, to get the convergence condition, we have\n",
    "$$\n",
    "x^{(k+1)} = -L_*^{-1}U x^{(k)} + L_*^{-1}b = Cx^{(k)}+ L_*^{-1}b\n",
    "$$\n",
    "The Gauss-Seidel iteration converges for every $x^{(0)}$ if and only if all the eigenvalues of $C$ have absolute value less than 1. However, calculating eigenvalues of $C$ is a more computational expensive operation than the iteration itself. We usually estimate the convergence using **Sufficient Convergence Condition**:\n",
    "$$\n",
    "\\Vert{C}\\Vert< 1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here $\\Vert C \\Vert $ is the **matrix norm**. There are several frequently used ones.\n",
    "\n",
    "**Frobenius norm**: root square of all the element summations.\n",
    "$$ \\Vert C \\Vert = \\sqrt{\\sum_{j=1}^n\\sum_{k=1}^n c_{jk}^2}$$\n",
    "**Column \"sum\" norm**: the greatest of the sums of the elements in a column of $C$.\n",
    "$$ \\Vert C \\Vert = \\max_k \\sum_{j=1}^n \\vert c_{jk} \\vert $$\n",
    "**Row \"sum\" norm**: the greatest of the sums of the elements in a row of $C$.\n",
    "$$ \\Vert C \\Vert = \\max_j \\sum_{k=1}^n \\vert c_{jk} \\vert $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Gauss_Seidel_iteration(A,b,x0,max_iter=100,err_tol = 1e-3,print_process=False):\n",
    "    ''' Function that utilizes Gauss Seidel Iteration to solve\n",
    "    a linear system Ax = b\n",
    "    usage: x = Gauss_Seidel_iteration(A,b,x0,max_iter=100,err_tol = 1e-3,print_process=False)\n",
    "    input: \n",
    "        A: nxn coefficient matrix, full rank\n",
    "        b: nx1 vector\n",
    "        x0: initial gauss of solution\n",
    "        max_iter: maximum number of iterations\n",
    "        err_tol: error tolorate of the solution\n",
    "        print_process: boolean, whether to print the step-by-step process \n",
    "    output:\n",
    "        solution x\n",
    "    written by Ge Jin, gjin@mines.edu, 06/2019\n",
    "    '''\n",
    "    n = A.shape[0]\n",
    "    x = x0.copy()\n",
    "    iter_success = False\n",
    "    for niter in range(max_iter):\n",
    "        if print_process:\n",
    "            print(x)\n",
    "        old_x = x.copy()\n",
    "        for i in range(n):\n",
    "            x[i] = 1/A[i,i]*(b[i]-np.sum(A[i,:]*x)+A[i,i]*x[i])\n",
    "        max_diff = np.max(np.abs(x-old_x))\n",
    "        if max_diff < err_tol:\n",
    "            iter_success = True\n",
    "            break\n",
    "    if not iter_success:\n",
    "        print('Max iteration number reached! Solution may not be accurate!')\n",
    "        \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[100. 100. 100. 100.]\n",
      "[100.   100.    75.    68.75]\n",
      "[93.75   90.625  65.625  64.0625]\n",
      "[89.0625   88.28125  63.28125  62.890625]\n",
      "[87.890625   87.6953125  62.6953125  62.59765625]\n",
      "[87.59765625 87.54882812 62.54882812 62.52441406]\n",
      "[87.52441406 87.51220703 62.51220703 62.50610352]\n",
      "[87.50610352 87.50305176 62.50305176 62.50152588]\n",
      "[87.50152588 87.50076294 62.50076294 62.50038147]\n",
      "[87.50038147 87.50019073 62.50019073 62.50009537]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([87.50009537, 87.50004768, 62.50004768, 62.50002384])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# demo of Gauss-Seidel Iteration\n",
    "A = np.array([[1,-0.25,-0.25,0],[-0.25,1,0,-0.25],[-0.25,0,1,-0.25],[0,-0.25,-0.25,1]])\n",
    "b = np.array([50,50,25,25])\n",
    "x0 = np.ones(4)*100\n",
    "Gauss_Seidel_iteration(A,b,x0,print_process=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.64 ms ± 2.2 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "31.5 ms ± 10.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Maximum Error of iteration solution:  1.9868014191570804e-05\n"
     ]
    }
   ],
   "source": [
    "# A comparison between Gauss Elimination and Gauss-Seidel Iteration\n",
    "N = 100\n",
    "A = np.random.rand(N,N) + 100*np.identity(N)\n",
    "b = np.random.rand(N)*100\n",
    "x0 = np.random.rand(N)\n",
    "\n",
    "%timeit Gauss_Seidel_iteration(A,b,x0)\n",
    "%timeit Gauss_elimination(A,b)\n",
    "\n",
    "x_gs = Gauss_Seidel_iteration(A,b,x0)\n",
    "x_ge = Gauss_elimination(A,b)\n",
    "print('Maximum Error of iteration solution: ',np.max(np.abs(x_gs - x_ge)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Norms and Condition Number\n",
    "A computational problem is called **ill-conditioned** (or ill-posed) if \"small\" changes in the data (the input) cause \"large\" changes in the solution (the output). In contrast, a problem is called **well-conditioned** (or well-posed) if \"small\" changes in the data cause only \"small\" changes in the solution. \n",
    "Keep in mind these concepts are qualitative. The definition of \"small\" and \"large\" depends on the accuracy of the data and the error tolerance of the solution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a simple linear system as an example. Say we would like to calculate the crossing point location of two straight lines. If the two lines are nearly parallel to each other, this problem is ill-conditioned. Because a small change in the intercept can change the location of the crossing point dramatically, as shown in the figure. \n",
    "\n",
    "<img src=\"files/figs/ill_condition.png\" width=\"400\">\n",
    "\n",
    "Here is an example of ill-conditioned linear system:\n",
    "$$\n",
    "0.9999 x_1 + 1.0001 x_2 = 1 \\\\\n",
    "x_1 - x_2 = 1\n",
    "$$\n",
    "By inducing a small error $\\epsilon$ to the intercept on the second equation, we have\n",
    "$$\n",
    "0.9999 x_1 + 1.0001 x_2 = 1 \\\\\n",
    "x_1 - x_2 = 1 + \\epsilon.\n",
    "$$\n",
    "This system has the solution $x_1 = 0.5 + 5000.5 \\epsilon$, $x_2 = -0.5 + 4999.5 \\epsilon$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to show that ill-conditioning of a linear system and of its coefficient matrix $A$ can be measured by a number, the **condition number** $\\kappa(A)$. Other measures for ill-conditioning have also been proposed, but $\\kappa(A)$ is probably the most widely used one. $\\kappa(A)$ is defined in terms of norm. To reach our goal, we discuss in three steps:\n",
    "1. **Vector norms**\n",
    "2. **Matrix norms**\n",
    "3. **Condition number**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vector Norms\n",
    "A **vector norm** for column vectors $x=[x_j]$ with $n$ components is a generalized length of the vector. It is denoted by $\\Vert x \\Vert$ and is defined by four properties,\n",
    "1. $\\Vert x \\Vert$ is a nonnegative real number.\n",
    "2. $\\Vert x \\Vert = 0$ if and only if $x=0$.\n",
    "3. $\\Vert kx \\Vert = |k| \\Vert x \\Vert$ for all $k$. \n",
    "4. $\\Vert x+y \\Vert \\leq \\Vert x \\Vert + \\Vert y \\Vert$\n",
    "\n",
    "The most important norms in connection with computations is the **p-norm** defined by\n",
    "$$\n",
    "\\Vert x \\Vert_p = \\left( \\sum_{i=1}^{n} |x_i|^p \\right)^{1/p}\n",
    "$$\n",
    "where $p$ is a fixed number. The most commonly used ones are \n",
    "$$\n",
    "\\begin{align}\n",
    "\\Vert x \\Vert_1 &= |x_1| + |x_2| + ... + |x_n|  &(l_1\\text{-norm})  \\\\\n",
    "\\Vert x \\Vert_2 &= \\sqrt{|x_1|^2 + |x_2|^2 + ... + |x_n|^2}  &(\\text{Euclidean or } l_2\\text{-norm}) \\\\ \n",
    "\\Vert x \\Vert_\\infty &= max_j |x_j|  &(l_\\infty \\text{-norm})  \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**\n",
    "What is the norms of vector $x^T = [\\begin{matrix}2&-3&0&1&-4\\end{matrix}]$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Norm\n",
    "If $A$ is an $n \\times n$ matrix and $x$ any vector with $n$ components, then $Ax$ is a vector with $n$ components. We now take a vector norm and consider $\\Vert x \\Vert$ and $ \\Vert Ax \\Vert$. One can prove that there is a number $c$ (depending on $A$) such that\n",
    "$$ \\Vert Ax \\Vert \\leq c \\Vert x \\Vert, \\ \\ \\ \\  \\text{for all } x $$\n",
    "The **matrix norm of $A$** is defined as the smallest possible $c$ valid for all $x$ ($\\neq 0$). Thus\n",
    "$$ \\Vert A \\Vert = \\min(c)= \\max \\frac{\\Vert Ax \\Vert}{\\Vert x \\Vert} $$\n",
    "or\n",
    "$$\n",
    "\\Vert A \\Vert = \\max_{\\Vert x \\Vert=1} \\Vert Ax \\Vert\n",
    "$$\n",
    "Note carefully that $\\Vert A \\Vert$ depends on the vector norm that we selected. In particular, one can show that\n",
    "+ for the $l_1$-norm one gets the column sum norm of $A$.\n",
    "+ for the $l_\\infty$-norm one gets the row sum norm of $A$.\n",
    "\n",
    "More importantly, by the norm definition, we have\n",
    "$$ \\Vert Ax \\Vert \\leq \\Vert A \\Vert \\Vert x \\Vert $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Condition Number of a Matrix\n",
    "We are now ready to introduce the key concept in our discussion of ill-conditioning, the **condition number** $\\kappa(A)$ of a (nonsingular) square matrix $A$, defined by\n",
    "$$ \\kappa(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert  $$\n",
    "For a linear system of equations $Ax=b$, if the **condition number of $A$ is small, then the problem is well-conditioned.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To prove this, let's assume the observation $b$ has a small error $\\epsilon$, which causes an inaccurate solution $\\tilde x$. \n",
    "$$ A \\tilde x = b + \\epsilon$$\n",
    "By substituting $Ax=b$, we have \n",
    "$$\\epsilon = A(x-\\tilde x) $$\n",
    "$$x-\\tilde x = A^{-1}\\epsilon$$\n",
    "Taking the norm on both side yields\n",
    "$$ \\Vert x-\\tilde x \\Vert \\leq \\Vert A^{-1} \\epsilon \\Vert \\leq \\Vert A^{-1} \\Vert \\Vert \\epsilon \\Vert $$\n",
    "Division by $\\Vert x \\Vert $ on both side finally gives\n",
    "$$\n",
    "\\frac{\\Vert x- \\tilde x \\Vert}{\\Vert x \\Vert} \\leq\n",
    "\\frac{1}{\\Vert x \\Vert} \\Vert A^{-1} \\Vert\\Vert \\epsilon \\Vert \n",
    "$$\n",
    "Because $b=Ax$, $\\Vert b \\Vert = \\Vert Ax \\Vert \\leq \\Vert A \\Vert \\Vert x \\Vert$, we have\n",
    "$$\n",
    "\\frac{1}{\\Vert x \\Vert} \\leq \\frac{\\Vert A \\Vert}{\\Vert b \\Vert}\n",
    "$$\n",
    "Thus,\n",
    "$$\n",
    "\\frac{\\Vert x- \\tilde x \\Vert}{\\Vert x \\Vert} \\leq\n",
    "\\frac{\\Vert A \\Vert}{\\Vert b \\Vert} \\Vert A^{-1} \\Vert\\Vert \\epsilon \\Vert \n",
    "= \\kappa(A) \\frac{\\Vert \\epsilon \\Vert}{\\Vert b \\Vert}\n",
    "$$\n",
    "Hence if $\\kappa(A)$ is small, a small error $\\Vert \\epsilon \\Vert$ implies a small relative error $\n",
    "{\\Vert x- \\tilde x \\Vert}/{\\Vert x \\Vert} $, so that the system is well-conditioned. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "1. Calculate the condition number of matrix \n",
    "$$\n",
    "A = \\left[\n",
    "\\begin{matrix}\n",
    "5 & 1 & 1 \\\\\n",
    "1 & 4 & 2 \\\\\n",
    "1 & 2 & 4 \n",
    "\\end{matrix}\n",
    "\\right], \n",
    "\\ \\ \\ \\ \\ \\ \\ \n",
    "\\text{which has the inverse}\n",
    "\\ \\ \\ \\ \\ \\ \\ \n",
    "A^{-1} = \\frac{1}{56}\\left[\n",
    "\\begin{matrix}\n",
    "12 & -2 & -2 \\\\\n",
    "-2 & 19 & -9 \\\\\n",
    "-2 & -9 & 19 \n",
    "\\end{matrix}\n",
    "\\right]. \n",
    "$$\n",
    "\n",
    "2. Calculate the condition number of matrix\n",
    "$$\n",
    "A = \\left[\n",
    "\\begin{matrix}\n",
    "0.9999 & -1.0001  \\\\\n",
    "1.0000 & -1.0000  \n",
    "\\end{matrix}\n",
    "\\right], \n",
    "\\ \\ \\ \\ \\ \\ \\ \n",
    "\\text{which has the inverse}\n",
    "\\ \\ \\ \\ \\ \\ \\ \n",
    "A^{-1} = \\left[\n",
    "\\begin{matrix}\n",
    "-5000 & 5000.5 \\\\\n",
    "-5000 & 4999.5 \n",
    "\\end{matrix}\n",
    "\\right]. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Further comments on Condition Numbers**\n",
    "\n",
    "1. There is no sharp dividing line between \"well-conditioned\" and \"ill-conditioned\", but generally the situation will get worse as we go from systems with small $\\kappa(A)$ to systems with larger $\\kappa(A)$. \n",
    "2. Ill-conditioned matrix takes longer to converge using Gauss-Seidel iteration. \n",
    "3. With modern computer science development, numerical algorithms can handle ill-conditioned matrix much better than they used to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Acknowledgement\n",
    "\n",
    "Most of this teaching material is based on:\n",
    "\n",
    "Kreyszig, E., 2018. Advanced Engineering Mathematics, 10-th edition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "242.023px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}