03_Numerical_Linear_Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](figs/title.png)\n",
    "\n",
    "<h1><center>Module 06: Numerical Linear Algebra</center></h1>\n",
    "\n",
    "In this chapter, we will discuss how to solve linear systems of equations numerically, and also introduce some basic parallel computing algorithms using OpenMP. We start with simple matrix multiplication, and then introduce performing Gauss elimination numerically, which is similar to how human solves 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": [
    "## Matrix Multiplication and its Parallel Computing Algorithm\n",
    "\n",
    "Matrix multiplication serves as the fundamental of linear algebra. The numerical method to calculate a matrix multiplication is fairly straightforward:\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_{m1} & a_{m2} & ... & a_{mn}\n",
    "\\end{matrix} \\right], \\ \\ \n",
    "B = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "b_{11} & b_{12} & ... & b_{1p} \\\\\n",
    "b_{21} & b_{22} & ... & b_{2p} \\\\\n",
    ". & . & . & . \\\\\n",
    "b_{n1} & b_{n2} & ... & b_{np}\n",
    "\\end{matrix} \\right], \\ \\ \n",
    "C = AB = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "c_{11} & c_{12} & ... & c_{1p} \\\\\n",
    "c_{21} & c_{22} & ... & c_{2p} \\\\\n",
    ". & . & . & . \\\\\n",
    "c_{m1} & c_{m2} & ... & c_{mp}\n",
    "\\end{matrix} \\right], \\ \\ \n",
    "$$\n",
    "where\n",
    "$$\n",
    "c_{ij} = \\sum_{k=1}^{n}a_{ik}b_{kj}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix multiplication has a computational complexity of $O(n^3)$. Because each $c_{ij}$ can be calculated independently, matrix multiplication can be highly parallelized in the OpenMP environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Direct implement: 5.409850597381592 s\n",
      "Vectorization implement: 0.23753571510314941 s\n",
      "Package implement: 0.010262727737426758 s\n"
     ]
    }
   ],
   "source": [
    "# Python realization of matrix multiplication\n",
    "from time import time\n",
    "N = 200;M = 200\n",
    "A = np.random.rand(M,N)\n",
    "B = np.random.rand(N,M)\n",
    "C = np.random.rand(M,M)\n",
    "\n",
    "# direct implement\n",
    "tic = time()\n",
    "for i in range(M):\n",
    "    for j in range(M):\n",
    "        for k in range(N):\n",
    "            C[i,j] += A[i,k]*B[k,j]\n",
    "toc = time()\n",
    "print('Direct implement: {} s'.format(toc-tic))\n",
    "\n",
    "# vectorization\n",
    "tic = time()\n",
    "for i in range(M):\n",
    "    for j in range(M):\n",
    "        C[i,j] = np.sum(A[i,:]*B[:,j])\n",
    "toc = time()\n",
    "print('Vectorization implement: {} s'.format(toc-tic))\n",
    "\n",
    "# for package usage\n",
    "tic = time()\n",
    "C = A@B\n",
    "toc = time()\n",
    "print('Package implement: {} s'.format(toc-tic))\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# C and OpenMP realization of matrix multiplication\n",
    "void matrix_mul_parallel(int m, int n, int p, float A[m][n], float B[n][p], float C[m][p], int thread_num){\n",
    "    omp_set_num_threads(thread_num);\n",
    "    #pragma omp parallel for \n",
    "    for (int i=0;i<m;i++)\n",
    "        for (int j=0;j<p;j++){\n",
    "            float sum = 0;\n",
    "            for (int k=0;k<n;k++)\n",
    "                sum += A[i][k]*B[k][j];\n",
    "            C[i][j] = sum;\n",
    "        }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Strassen Algorithm\n",
    "\n",
    "Strassen algorithm replaces a large matrix multiplication operation with a series of smaller sub-matrices multiplications and summations. It is slightly better than the direct multiplication, but the idea promoted more research on matrix multiplication which led to faster methods. Check [wikipedia](https://en.wikipedia.org/wiki/Strassen_algorithm) for more details. \n",
    "\n",
    "$$ C=AB $$\n",
    "$$\n",
    "A = \\left[\n",
    "\\begin{matrix}\n",
    "A_{1,1} & A_{1,2} \\\\\n",
    "A_{2,1} & A_{2,2}\n",
    "\\end{matrix}\n",
    "\\right],\n",
    "B = \\left[\n",
    "\\begin{matrix}\n",
    "B_{1,1} & B_{1,2} \\\\\n",
    "B_{2,1} & B_{2,2}\n",
    "\\end{matrix}\n",
    "\\right],\n",
    "C = \\left[\n",
    "\\begin{matrix}\n",
    "C_{1,1} & C_{1,2} \\\\\n",
    "C_{2,1} & C_{2,2}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "The naive algorithm would be:\n",
    "$$\n",
    "C_{1,1} = A_{1,1}B_{1,1} + A_{1,2}B_{2,1} \\\\\n",
    "C_{1,2} = A_{1,1}B_{1,2} + A_{1,2}B_{2,2} \\\\\n",
    "C_{2,1} = A_{2,1}B_{1,1} + A_{2,2}B_{2,1} \\\\\n",
    "C_{2,2} = A_{2,1}B_{1,2} + A_{2,2}B_{2,2} \n",
    "$$\n",
    ", which still need 8 multiplications to calculate the $C_{i,j}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Strassen algorithm defines instead new matrices:\n",
    "$$\n",
    "M_1 := (A_{1,1}+A_{2,2})(B_{1,1} + B_{2,2}) \\\\ \n",
    "M_2 := (A_{2,1}+A_{2,2})B_{1,1} \\\\\n",
    "M_3 := A_{1,1}(B_{1,2} - B_{2,2}) \\\\\n",
    "M_4 := A_{2,2}(B_{2,1} - B_{1,1}) \\\\\n",
    "M_5 := (A_{1,1}+A_{1,2})B_{2,2} \\\\\n",
    "M_6 := (A_{2,1}-A_{1,1})(B_{1,1} + B_{1,2}) \\\\ \n",
    "M_7 := (A_{1,2}-A_{2,2})(B_{2,1} + B_{2,2})\n",
    "$$\n",
    "in this case, only 7 multiplications are used instead of 8. We may now express the $C$ in terms of $M$:\n",
    "$$\n",
    "C_{1,1} = M_1+M_4-M_5+M_7 \\\\\n",
    "C_{1,2} = M_3+M_5 \\\\\n",
    "C_{2,1} = M_2+M_4 \\\\\n",
    "C_{2,2} = M_1-M_2+M_3+M_6 \n",
    "$$\n",
    "\n",
    "The complexity of Strassen algorithm is $\\approx O(N^{2.8})$. The fastest known algorithm is called [Coppersmith-Winograd algorithm](https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm), which has a complexity of $O(n^{2.375477})$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix Inversion\n",
    "\n",
    "In this section, we focuses on solving a linear system of equations using matrix inversion and numerical methods. \n",
    "\n",
    "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 **coefficient** $a_{jk}$ and $b_j$ are given. \n",
    "\n",
    "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": [
    "### Gauss Elimination\n",
    "\n",
    "[Gauss elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) and [back substitution](https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution) are common methods to solve linear equation set using matrix method. In this section we discuss the numerical realization of these methods. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform Gauss elimination, we first need to form an **augmented matrix** of the system by combining matrix $A$ and vector $b$:\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",
    "\n",
    "Similar to the linear equation sets, there are **three elementary row operations** we can perform on the augmented matrix $\\tilde A$. Applying these rules 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 rules for linear equation system: \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": [
    "The goal of Gauss Elimination is to reduce the coefficient matrix $A$ from a full matrix to a **triangular** matrix, through a combination of these three elementary row operations. Then the triangular system can be easily solved by **back substitution** operation. \n",
    "\n",
    "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",
    "\n",
    "We can solve this entire system from the last equation, and step by step move backwards using **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$ in the Augumented matrix $\\tilde{A}$ into a triangular form, we eliminate coefficient of $x_1$ from row 2 to row $n$ of $\\tilde A$ first by subtracting the suitable multiples of row 1 from the other rows. In this step, the first row is called the **pivot equation** and the coefficient $a_{11}$ is called the **pivot**. We changes all the rows of $\\tilde A$, except for the pivot row. After this step, all the elements of the first column of $\\tilde A$ become zeros except for the first row. \n",
    "\n",
    "Then we repeat this process and eliminate the coefficient of $x_2$ from row 3 to row $n$ by subtracting suitable multiples of row 2 from all the rows below. At this point, all elements of the first and second columns below the matrix diagnal are zeros. \n",
    "\n",
    "We keep repeating this process for $n$ times, and eventually reduce the matrix $A$ to a triangular matrix.\n",
    "\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. When the problem is solved numerically, it **should be** as large as possible in absolute value to avoid the numeric instability. \n",
    "\n",
    "As a result, 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": 1,
   "metadata": {},
   "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.986541509628296"
      ]
     },
     "execution_count": 4,
     "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": 5,
   "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": 6,
   "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 include\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 $(n-k)(n-k+1)$ subtractions. This operation is performed $n-1$ steps, so $k$ goes from 1 to $n-1$. 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 represent it as\n",
    "$$\n",
    "f(n) = O(n^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the back substituion, each step we make $n-i$ multiplications and as many subtractions, as well as 1 division. Hence the number of operations 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 on average, 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 is through **LU-Factorization**. \n",
    "\n",
    "An LU-factorization of a given square matrix $A$ is of the form\n",
    "$$ A=LU $$\n",
    "where $L$ is a **lower triangular** matrix and $U$ is a **upper triangular** matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be proved that for any nonsingular full matrix, the rows can be reordered so that the resulting matrix A has an LU-factorization. It turns out that $L$ is 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; \\;\\; j\\geq 2\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to Gauss elimination, the advantage of LU factorization is that $L$ and $U$ can be computed directly, without solving simultaneous equations. The LU factorization is independent from the vector $b$. It has particular advantage if the same matrix needs to be inverted multiple times with different $b$.\n",
    "\n",
    "Once $L$ and $U$ are calculated, 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 donot 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 to 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** with $A$ and an identity matrix $I$:\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": 7,
   "metadata": {},
   "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": 8,
   "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 other methods we discussed before are **direct methods** for solving linear systems of equations. These methods provide solutions after an certain amount of computation that can be specified in advance. In contrast, there are a large group of **indirect** or **iterative methods**, where 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 necessary for achieving a required accuracy. \n",
    "Iterative methods in general have less computational cost comparing to direct methods, especially for very large or very sparse matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gauss-Seidel Iteration Method\n",
    "Gauss-Seidel method is an iterative method of great importance for solving linear systems. In this section, we first derive the general formulas for the method.\n",
    "\n",
    "We can rearrange and scale the equations such that the diagonal coefficients $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 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": [
    "#### Jacobi Method for Parallel Computing\n",
    "\n",
    "It can be observed that Gaussian-Seidel iteration method is very difficult to parallelize, because of the input dependency at each step. Jacobi method is a slightly altered version of Gaussian-Seidel method, which is more parallel friendly. \n",
    "\n",
    "For Jacobi method, each iteration can be presented as:\n",
    "$$\n",
    "x_i^{(k+1)} = \\frac{1}{a_{ii}} \\left(\n",
    "b_i - \\sum_{j=1}^{i-1} a_{ij}x_j^{(k)} - \\sum_{j=i+1}^n a_{ij}x_j^{(k)}\n",
    "\\right),\\ \\ \\ \\ i=1,2,...,n\n",
    "$$\n",
    "\n",
    "For Jacobi method, the new iteration $x^{(k+1)}$ only depends on the value of current iteration $x^{(k)}$, thus each element of $x^{(k+1)}$ can be calculated at different nodes to accelerate the computation.\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](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) of $C$ have the absolute value less than 1. However, calculating eigenvalues of $C$ is a more computational expensive operation than the iteration itself. As a result, 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](https://en.wikipedia.org/wiki/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 $$\n",
    "\n",
    "We will discuss matrix and vector norm in the following section in more details. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.97 ms ± 187 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "65.9 ms ± 1.83 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Maximum Error of iteration solution:  1.9403455026323346e-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) # increase the diagnal elements to stablize the iteration\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 (vector $b$) can cause \"large\" changes in the solution (unknown vector $x$). 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 and relative. 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": [
    "We demonstrate this concept using 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.\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$. It can be observed that the error in the solution $x$ is significantly magnified by the error in the data $b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this section is to show that the ill-conditioning of a linear system and of its coefficient matrix $A$ can be measured by a quantitative 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. We first discuss the norm of vector and matrix, then discuss the defination of 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 vector length. 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 for computation purposes are 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$ a vector with $n$ components, then $Ax$ is a vector with $n$ components. We can 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 max of column sum norm of $A$.\n",
    "+ for the $l_\\infty$-norm one gets the max of 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 **condition number** $\\kappa(A)$ of a (nonsingular) square matrix $A$, which is the key concept of ill-conditioning discussion. The condition number is 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": [
    "Let's assume that 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 boundary between \"well-conditioned\" and \"ill-conditioned\" matrices, but generally the solution will get less accurate 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",
    "Some of this teaching material is from:\n",
    "\n",
    "Kreyszig, E., 2018. Advanced Engineering Mathematics, 10-th edition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "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": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}