F0D-07-DecompositionMethods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FLIP (00): Data Science \n",
    "**(Module 03: Linear Algebra)**\n",
    "\n",
    "---\n",
    "- Materials in this module include resources collected from various open-source online repositories.\n",
    "- You are free to use,but NOT allowed to change and distribute this package.\n",
    "\n",
    "Prepared by and for \n",
    "**Student Members** |\n",
    "2006-2023 [TULIP Lab](http://www.tulip.org.au), Australia\n",
    "\n",
    "---\n",
    "\n",
    "## Session 15 Decomposition Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decomposition Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Gaussian elimination we performed row (or column) operations to reduce the linear system\n",
    "\n",
    "\\begin{equation}\n",
    "  A {\\bf x} = {\\bf b}\n",
    "\\end{equation}\n",
    "\n",
    "to a form that was easy to solve, essentially by transforming to an equivalent system where the matrix is triangular."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In *decomposition* methods we instead write the matrix as a product of other matrices in such a way that we have a set of easy-to-solve problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, suppose we can write\n",
    "\n",
    "\\begin{equation}\n",
    "  A = L U\n",
    "\\end{equation}\n",
    "\n",
    "where $L$ is a lower triangular matrix and $U$ upper triangular. In this case, the original linear system is equivalent to the two systems\n",
    "\n",
    "\\begin{equation}\n",
    "  A {\\bf x} = {\\bf b} \\qquad \\Leftrightarrow \\qquad \\left\\{ \\begin{aligned} L {\\bf y} & = {\\bf b}, \\\\ U {\\bf x} & = {\\bf y}. \\end{aligned} \\right. \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As both $L$ and $U$ are triangular, the two systems they define can easily be solved using forwards and backwards substitution respectively. Therefore, *if* we can decompose the matrix $A$ as the product of triangular matrices, the linear system is easily solved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the system\n",
    "\n",
    "\\begin{equation}\n",
    "  \\begin{pmatrix} 2 & 1 & -1 \\\\ 4 & 1 & 0 \\\\ -2 & -3 & 8 \\end{pmatrix} {\\bf x} = \\begin{pmatrix} 0 \\\\ 6 \\\\ 12 \\end{pmatrix}.\n",
    "\\end{equation}\n",
    "\n",
    "We will see later that the matrix can be decomposed as\n",
    "\n",
    "\\begin{equation}\n",
    "  \\begin{pmatrix} 2 & 1 & -1 \\\\ 4 & 1 & 0 \\\\ -2 & -3 & 8 \\end{pmatrix} = \\begin{pmatrix} 1 & 0 & 0 \\\\ 2 & 1 & 0 \\\\ -1 & 2 & 1 \\end{pmatrix} \\begin{pmatrix} 2 & 1 & -1 \\\\ 0 & -1 & 2 \\\\ 0 & 0 & 3 \\end{pmatrix}.\n",
    "\\end{equation}\n",
    "\n",
    "We can quickly check that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2.0, 1.0, -1.0] ,[4.0,  1.0, 0.0], [-2.0, -3.0, 8.0]])\n",
    "L = np.array([[1.0, 0.0,  0.0] ,[2.0,  1.0, 0.0], [-1.0,  2.0, 1.0]])\n",
    "U = np.array([[2.0, 1.0, -1.0] ,[0.0, -1.0, 2.0], [ 0.0,  0.0, 3.0]])\n",
    "np.dot(L, U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then solve $L {\\bf y} = {\\bf b}$ by forward substitution:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\begin{pmatrix}\n",
    "    1 & 0 & 0 \\\\\n",
    "    2 & 1 & 0 \\\\\n",
    "   -1 & 2 & 1\n",
    "  \\end{pmatrix}\n",
    "  \\begin{pmatrix}\n",
    "    y_1 \\\\ y_2 \\\\ y_3\n",
    "  \\end{pmatrix} =\n",
    "  \\begin{pmatrix}\n",
    "    0 \\\\ 6 \\\\ -12\n",
    "  \\end{pmatrix}\n",
    "  \\Rightarrow \\left\\{\n",
    "  \\begin{aligned}\n",
    "    y_1 & = 0, \\\\\n",
    "    y_2 & = 6, \\\\\n",
    "    y_3 & = -24.\n",
    "  \\end{aligned}\n",
    "  \\right.\n",
    "\\end{equation}\n",
    "\n",
    "Finally solve $U {\\bf x} = {\\bf y}$ by backward substitution:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\begin{pmatrix}\n",
    "    2 & 1 & -1 \\\\\n",
    "    0 & -1 & 2 \\\\\n",
    "    0 & 0 & 3\n",
    "  \\end{pmatrix}\n",
    "  \\begin{pmatrix}\n",
    "    x_1 \\\\ x_2 \\\\ x_3\n",
    "  \\end{pmatrix} =\n",
    "  \\begin{pmatrix}\n",
    "    0 \\\\ 6 \\\\ -24\n",
    "  \\end{pmatrix}\n",
    "  \\Rightarrow \\left\\{\n",
    "  \\begin{aligned}\n",
    "    x_1 & = 7, \\\\\n",
    "    x_2 & = -22, \\\\\n",
    "    x_3 & = -8.\n",
    "  \\end{aligned}\n",
    "  \\right.\n",
    "\\end{equation}\n",
    "\n",
    "Again we can cross-check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([0.0, 6.0, -12.0])\n",
    "la.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uniqueness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "  The factorisation is not unique: $L$ and $U$ together have $n^2+n$\n",
    "  free coefficients; $A$ only has $n^2$.  Can freely choose $n$\n",
    "  coefficients.   Two obvious choices are\n",
    "  \n",
    "1. Doolittle's factorisation: The diagonal entries of $L$ are 1.\n",
    "2. Crout's factorisation: The diagonal entries of $U$ are 1.\n",
    "\n",
    "  Then use the explicit matrix multiplication formula, and work from\n",
    "  the top left corner. Each coefficient obeys\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{i j} = \\sum_{s=1}^n \\ell_{i s} u_{s j} = \\sum_{s=1}^{\\min(i,j)}\n",
    "     \\ell_{i s} u_{s j}\n",
    "  \\end{equation}\n",
    "  \n",
    "  where the last equality holds because $L$ and $U$ are triangular:\n",
    "  $\\ell_{i s} = 0$ for $s > i$ and $u_{s j} = 0$ for $s > j$.\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hand factorisation example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "  The previous example had the matrix\n",
    "  \n",
    "  \\begin{equation}\n",
    "    A =\n",
    "    \\begin{pmatrix}\n",
    "      2 & 1 & -1 \\\\\n",
    "      4 & 1 & 0 \\\\\n",
    "      -2 & -3 & 8\n",
    "    \\end{pmatrix}\n",
    "  \\end{equation}\n",
    "  \n",
    "  which we decompose using\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{i j} = \\sum_{s=1}^{\\min(i, j)} \\ell_{i s} u_{s j}.\n",
    "  \\end{equation}\n",
    "\n",
    "  First coefficient, $i = j = 1$, obeys\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{1 1}  = 2 = \\ell_{1 1} u_{1 1}.\n",
    "  \\end{equation}\n",
    "  \n",
    "  Use our freedom to set some coefficients to choose $\\ell_{1 1} = 1$, for example, which gives $u_{1 1} = a_{1 1} = 2$.\n",
    "  \n",
    "  \n",
    "  Now consider the first row of $U$ ($i = 1$, $j$ free) and the first column of $L$ ($i$ free, $j = 1$). The first row of $U$ obeys\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{1 j} = \\ell_{1 1} u_{1 j} = u_{1 j}\n",
    "  \\end{equation}\n",
    "  \n",
    "  and the first column of $L$ obeys\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{i 1} = \\ell_{i 1} u_{1 1} = 2 \\ell_{i 1}.\n",
    "  \\end{equation}\n",
    "  \n",
    "  Therefore we know that\n",
    "  \n",
    "  \\begin{equation}\n",
    "    L =\n",
    "    \\begin{pmatrix}\n",
    "      1 & 0 & 0 \\\\\n",
    "      \\tfrac{4}{2} & ?? & 0 \\\\\n",
    "      \\tfrac{-2}{2} & ?? & ??\n",
    "    \\end{pmatrix}, \\quad\n",
    "    U =\n",
    "    \\begin{pmatrix}\n",
    "      2 & 1 & -1 \\\\\n",
    "      0 & ?? & ?? \\\\\n",
    "      0 & 0 & ??\n",
    "    \\end{pmatrix}.\n",
    "  \\end{equation}\n",
    "\n",
    "  Go to second row and column.\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{2 2} = \\ell_{2 1} u_{1 2} + \\ell_{2 2} u_{2 2}.\n",
    "  \\end{equation}\n",
    "  \n",
    "  Already computed entries $\\ell_{2 1} = 2$, $u_{1 2} = 1$. Again set $\\ell_{2 2} = 1$, giving\n",
    "  \n",
    "  \\begin{equation}\n",
    "    a_{2 2} = 1 = 2 + u_{2 2} \\Rightarrow u_{2 2} = -1.\n",
    "  \\end{equation}\n",
    "\n",
    "  As in the previous step, consider the second row of $U$ and\n",
    "  second column of $L$, finding\n",
    "  \n",
    "  \\begin{align}\n",
    "    u_{2 j}    & = a_{2 j} - \\ell_{2 1} u_{1 j}, \\\\\n",
    "    \\ell_{i 2} & = \\left(a_{i 2} - \\ell_{i 1} u_{1 2} \\right) / u_{2 2},\n",
    "  \\end{align}\n",
    "  \n",
    "  which implies that\n",
    "  \n",
    "  \\begin{align}\n",
    "    u_{2 3}    & = 0 - 2 \\times (-1) = 2, \\\\\n",
    "    \\ell_{3 2} & = \\left(-3 - (-1) \\times 1\\right) / (-1) = 2.\n",
    "  \\end{align}\n",
    "  \n",
    "  We now have the first two rows and columns:\n",
    "  \n",
    "  \\begin{equation}\n",
    "    L =\n",
    "    \\begin{pmatrix}\n",
    "      1 & 0 & 0 \\\\\n",
    "      2 & 1 & 0 \\\\\n",
    "      -1 & 2 & ??\n",
    "    \\end{pmatrix}, \\quad\n",
    "    U =\n",
    "    \\begin{pmatrix}\n",
    "      2 & 1 & -1 \\\\\n",
    "      0 & -1 & 2 \\\\\n",
    "      0 & 0 & ??\n",
    "    \\end{pmatrix}.\n",
    "  \\end{equation}\n",
    "  \n",
    "  Continuing as above, use last free choice to set $\\ell_{3 3} =\n",
    "  1$. Finally, as above, compute that $u_{3 3} = 3$. So we finally\n",
    "  have\n",
    "  \n",
    "  \\begin{equation}\n",
    "    L =\n",
    "    \\begin{pmatrix}\n",
    "      1 & 0 & 0 \\\\\n",
    "      2 & 1 & 0 \\\\\n",
    "      -1 & 2 & 1\n",
    "    \\end{pmatrix}, \\quad\n",
    "    U =\n",
    "    \\begin{pmatrix}\n",
    "      2 & 1 & -1 \\\\\n",
    "      0 & -1 & 2 \\\\\n",
    "      0 & 0 & 3\n",
    "    \\end{pmatrix}.\n",
    "  \\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $LU$ Decomposition Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Write down the explicit formula for the matrix coefficients.\n",
    "\n",
    "  1. Work from the first row / column to the last.\n",
    "  2. Look at the diagonal entry of $A$: use the freedom to choose the value of either $L$ or $U$'s diagonal entry.\n",
    "  3. The appropriate row of $U$ and column of $L$ follows from the explicit formula as all other entries are known.\n",
    "  \n",
    "* If either $u_{k k}$ or $\\ell_{k k}$ are zero the simple algorithm fails; however, an $LU$ decomposition may still exist.\n",
    "* $U$ is the matrix that would be found using Gaussian elimination. $LU$ decomposition has the advantage that it holds for any ${\\bf b}$ in $A {\\bf x} = {\\bf b}$, whereas with Gaussian elimination the whole algorithm needs repeating.\n",
    "* As with Gaussian elimination, pivoting for accuracy is necessary for general matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example LU decomposition algorithm, taken from the solution to the worksheets. Note that, as it doesn't use pivoting, it may well fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LU_decomposition(A):\n",
    "    \"\"\"Perform LU decomposition using the Doolittle factorisation.\"\"\"\n",
    "    \n",
    "    L = np.zeros_like(A)\n",
    "    U = np.zeros_like(A)\n",
    "    N = np.size(A, 0)\n",
    "    \n",
    "    for k in range(N):\n",
    "        L[k, k] = 1\n",
    "        U[k, k] = (A[k, k] - np.dot(L[k, :k], U[:k, k])) / L[k, k]\n",
    "        for j in range(k+1, N):\n",
    "            U[k, j] = (A[k, j] - np.dot(L[k, :k], U[:k, j])) / L[k, k]\n",
    "        for i in range(k+1, N):\n",
    "            L[i, k] = (A[i, k] - np.dot(L[i, :k], U[:k, k])) / U[k, k]\n",
    "    \n",
    "    return L, U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check this code against the example above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L2, U2 = LU_decomposition(A)\n",
    "print(\"Compare the code and hand answers for L:\\n{},\\n{}\\nand for U:\\n{},\\n{}\".format(L2, L, U2, U))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditions for factorisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "  We *assumed* an $LU$ decomposition is always possible. However,\n",
    "  we then noted that the algorithm works only if both $u_{k k}$ and\n",
    "  $\\ell_{k k}$ are non-zero, which is not known in advance. \n",
    "  \n",
    "\n",
    "  Two key concepts used in determining if a factorisation exists.\n",
    "  \n",
    "1.  The *principal minor* of order $k$ of a matrix $A$ is the (sub-) matrix\n",
    "    \\begin{equation*}\n",
    "      \\begin{pmatrix}\n",
    "        a_{1 1} & \\cdots & a_{1 k} \\\\\n",
    "        \\vdots & \\ddots & \\vdots \\\\\n",
    "        a_{k 1} & \\cdots & a_{k k}\n",
    "      \\end{pmatrix} .\n",
    "    \\end{equation*}\n",
    "    \n",
    "2. A matrix is *strictly diagonally dominant* if\n",
    "    \\begin{equation*}\n",
    "      | a_{i i} | > \\sum_{\\substack{j = 1 \\\\ j \\ne i}}^n | a_{i j} |,\n",
    "        \\quad (1 \\le i \\le n).\n",
    "    \\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Theorems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "  Summarize the two key theorems; see references given in the notes:\n",
    "\n",
    "  **Theorem 1.1**: If all $n - 1$ leading principal minors of the $n\n",
    "  \\times n$ matrix $A$ are nonsingular, then an $LU$ decomposition of\n",
    "  $A$ exists. \n",
    "            \n",
    "\n",
    "  The relation between $LU$ decomposition and Gaussian elimination is\n",
    "  exploited to prove this theorem. It also shows that, if the matrix\n",
    "  is nonsingular, its rows can be permuted so that an $LU$\n",
    "  decomposition of the permuted matrix can be found.\n",
    "  \n",
    "  **Theorem 1.2}**: Every strictly diagonally dominant matrix is\n",
    "  nonsingular and has an $LU$ decomposition. \n",
    "  \n",
    "\n",
    "  Diagonal dominance is a very simple condition to check. Most\n",
    "  matrices involved in numerical methods for solving PDEs are\n",
    "  diagonally dominant, leading to the practical utility of this\n",
    "  theorem.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative decomposition methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cholesky factorisation - example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to decompose the symmetric, positive definite matrix\n",
    "\n",
    "  \\begin{equation}\n",
    "    A =\n",
    "    \\begin{pmatrix}\n",
    "      2 & 1 \\\\\n",
    "      1 & 3\n",
    "    \\end{pmatrix}\n",
    "  \\end{equation}\n",
    "  \n",
    "  as $A = L L^T$ where $L$ is lower triangular with positive diagonal\n",
    "  elements. Again we write the coefficients out explicitly and\n",
    "  look at the diagonal elements first:\n",
    "  \n",
    "  \\begin{align}\n",
    "                && A & = L L^T \\\\\n",
    "    \\Rightarrow &&\n",
    "    \\begin{pmatrix}\n",
    "      2 & 1 \\\\\n",
    "      1 & 3\n",
    "    \\end{pmatrix}\n",
    "    & =\n",
    "    \\begin{pmatrix}\n",
    "      \\ell_{1 1}^2 & \\ell_{1 1} \\ell_{2 1} \\\\\n",
    "      \\ell_{1 1} \\ell_{2 1} & \\ell_{2 1}^2 +  \\ell_{2 2}^2\n",
    "    \\end{pmatrix} \\\\ \n",
    "    \\Rightarrow && &\\left\\{ \n",
    "      \\begin{aligned}\n",
    "        \\ell_{1 1} & = \\sqrt{2} \\\\\n",
    "        \\ell_{2 1} & = 1 / \\sqrt{2} \\\\\n",
    "        \\ell_{2 2} & = \\sqrt{3 - \\ell_{2 1}^2} \\\\\n",
    "                   & = \\sqrt{5/2}\n",
    "      \\end{aligned}\\right. .\n",
    "  \\end{align}\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}