notebooks/Chapter 4 - LU Factorization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import sympy as sy\n",
    "sy.init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> LU Factorization</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LU factorization arises due its computational efficiency, it mainly facilitates solving the system of linear equations, though the flops (number of floating operations) of LU is still higher than Guassian-Jordon. \n",
    "\n",
    "Nonetheless LU factorization is particularly handy if you are computing multiple solutions of $A x =b$. \n",
    "\n",
    "One example is, if you have a set of $\\{b_i,\\ i \\in \\mathbb Z\\}$ to substitute in the system, such that\n",
    "$$\n",
    "Ax=b_1,\\quad Ax=b_2,\\quad Ax=b_3,\\quad Ax=b_4, \\quad...\n",
    "$$\n",
    "we only decompose $A$ once, but Guassian-Jordon algorithm will have to re-do operations with every augmented $[A |b_i]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> LU Algorithm</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We aim to decompose a matrix $A$ into $L$ and $U$, which represent _lower_ and _upper triangular matrix_ respectively. And $L$ must have $1$'s on its principal diagonal.\n",
    "$$\n",
    "A = LU\n",
    "$$\n",
    "For instance,\n",
    "$$\n",
    "A=\\underbrace{\\left[\\begin{array}{cccc}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "* & 1 & 0 & 0 \\\\\n",
    "* & * & 1 & 0 \\\\\n",
    "* & * & * & 1\n",
    "\\end{array}\\right]}_{L}\\underbrace{\\left[\\begin{array}{ccccc}\n",
    "\\blacksquare & * & * & * & * \\\\\n",
    "0 & \\blacksquare & * & * & * \\\\\n",
    "0 & 0 & 0 & \\blacksquare & * \\\\\n",
    "0 & 0 & 0 & 0 & 0\n",
    "\\end{array}\\right]}_{U}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plausibility of decomposition hinges on if $A$ can be converted into a upper triangular matrix after a series of row operations, i.e. \n",
    "\n",
    "$$E_p...E_2E_1  A = U$$\n",
    "\n",
    "Then \n",
    "\n",
    "$$A = (E_p...E_2E_1)^{-1}U = LU$$\n",
    "\n",
    "where \n",
    "\n",
    "$$L = (E_p...E_2E_1)^{-1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will hand calculate an example here:\n",
    "\n",
    "$$\n",
    "A =\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\3 & 4 & 6\\\\0 & 8 & 8\\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\-\\frac{1}{3} & 1 & 0\\\\0 & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\3 & 4 & 6\\\\0 & 8 & 8\\\\ \n",
    "\\end{bmatrix}}_{A}\n",
    "=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 8 & 8\\\\ \n",
    "\\end{bmatrix}}_{E_1A}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\0 &1  & 0\\\\0 & -\\frac{8}{3} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 8 & 8\\\\ \n",
    "\\end{bmatrix}}_{E_1A}\n",
    "=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 0 & -\\frac{8}{3}\\\\ \n",
    "\\end{bmatrix}}_{E_2E_1A=U}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "L^{-1} = E_2E_1 = \n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\0 &1  & 0\\\\0 & -\\frac{8}{3} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\-\\frac{1}{3} & 1 & 0\\\\0 & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\-\\frac{1}{3} & 1 & 0\\\\\\frac{8}{9} & -\\frac{8}{3} & 1\\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can compute $E_1^{-1}$ and $E_2^{-1}$ directly, because\n",
    "$$\n",
    "L = E_1^{-1}E_2^{-1}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct augmented matrices for $E_2$ and $E_1$\n",
    "$$\n",
    "[E_1| I]=\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "-\\frac{1}{3} & 1 & 0 & 0  & 1 & 0\\\\\n",
    "0 & 0 & 1 & 0 & 0 & 1\n",
    "\\end{array}\\right]\n",
    "\\sim\n",
    "\\left[\\begin{array}{ccc|ccc}\n",
    "1 &0 &0 &1 & 0 & 0\\\\\n",
    "0& 1& 0& \\frac{1}{3} & 1 & 0\\\\\n",
    "0 & 0 & 1 &0 & 0 & 1\n",
    "\\end{array}\\right]\n",
    "=[I|E_1^{-1}]\n",
    "$$\n",
    "\n",
    "$$\n",
    "[E_2| I]=\n",
    "\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0  & 1 & 0\\\\\n",
    "0 & -\\frac{8}{3} & 1 & 0 & 0 & 1\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\sim\n",
    "\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0  & 1 & 0\\\\\n",
    "0 & 0 & 1 & 0 &\\frac{8}{3}& 1\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=[I|E_2^{-1}]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we get the inverse of $E_1$ and $E_2$:\n",
    "$$\n",
    "E_1^{-1} = \\left[\\begin{array}{ccc} 1& 0& 0 \\\\ \\frac{1}{3}& 1& 0\\\\ 0& 0 &1  \\end{array}\\right]\\\\\n",
    "E_2^{-1} = \\left[\\begin{array}{ccc} 1& 0& 0 \\\\ 1& 1& 0\\\\ 0& \\frac{8}{3} &1  \\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "L = E_1^{-1}E_2^{-1}= \n",
    "\\left[\\begin{array}{ccc} 1& 0& 0 \\\\ \\frac{1}{3}& 1& 0\\\\ 0& 0 &1  \\end{array}\\right]\n",
    "\\left[\\begin{array}{ccc} 1& 0& 0 \\\\ 1& 1& 0\\\\ 0& \\frac{8}{3} &1 \\end{array}\\right]\n",
    "=\n",
    "\\left[\\begin{array}{ccc} 1& 0& 0 \\\\ \\frac{4}{3}&1& 0\\\\ 0& \\frac{8}{3} &1 \\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final result of LU decomposition:\n",
    "$$\n",
    "A = LU =\n",
    "\\underbrace{\n",
    "\\left[\\begin{array}{ccc} 1& 0& 0 \\\\ \\frac{4}{3}&1& 0\\\\ 0& \\frac{8}{3} &1 \\end{array}\\right]}_{L}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 0 & -\\frac{8}{3}\\\\ \n",
    "\\end{bmatrix}}_{E_2E_1A=U}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a look at SciPy results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9, 3, 6],\n",
       "       [3, 4, 6],\n",
       "       [0, 8, 8]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[9, 3, 6], [3, 4, 6], [0, 8, 8]]); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 1., 0.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        , 0.        ],\n",
       "       [0.        , 1.        , 0.        ],\n",
       "       [0.33333333, 0.375     , 1.        ]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[9., 3., 6.],\n",
       "       [0., 8., 8.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P, L, U = sp.linalg.lu(A)\n",
    "P\n",
    "L\n",
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9., 3., 6.],\n",
       "       [3., 4., 6.],\n",
       "       [0., 8., 8.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P@L@U # this is A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to notice the SciPy ```lu``` function give us more than $L$ and $U$, but also $P$, which is a _permutation matrix_. Theoretically, we don't need row switches to convert $A$ into $U$, but in some situations we must make row switches beforehand, otherwise decomposition will fail to materialize.\n",
    "\n",
    "Thus Scipy uses $PLU$ decomposition to make the procedure always possible\n",
    "$$\n",
    "A  = PLU\n",
    "$$\n",
    "\n",
    "Also $P = P^{-1}$, why? It's easier to analyze in augmented matrix expression, the inverse of row-switched elementary matrices are themselves.\n",
    "\n",
    "$$\n",
    "[P| I]=\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 0 & 1 & 0 & 1 & 0\\\\\n",
    "0 & 1 & 0 & 0 & 0 & 1\n",
    "\\end{array}\\right]\n",
    "\\sim\n",
    "\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0 & 0 & 1\\\\\n",
    "0 & 0 & 1 & 0 & 1 & 0\n",
    "\\end{array}\\right]=[I| P^{-1}]  =[P^{-1}|I]\n",
    "$$\n",
    "\n",
    "With these knowledge, what we are decomposing is actucally\n",
    "\n",
    "$$\n",
    "PA = LU\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\0 & -\\frac{3}{8} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\-\\frac{1}{3} & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 0 & 1\\\\0 & 1 & 0\\\\\n",
    "\\end{bmatrix}}_{E_0}\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\3 & 4 & 6\\\\0 & 8 & 8\\\\\n",
    "\\end{bmatrix}\\\\\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\0 & -\\frac{3}{8} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\-\\frac{1}{3} & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9& 3& 6\\\\0&8 &8\\\\ 3& 4& 6\n",
    "\\end{bmatrix}}_{E_0A}\\\\\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\0 & -\\frac{3}{8} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 &3 &6\\\\0 &8& 8 \\\\0& 3& 4\n",
    "\\end{bmatrix}}_{E_1E_0A}\\\\\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 &3 &6\\\\0 &8& 8 \\\\0& 0& -1\n",
    "\\end{bmatrix}}_{E_2E_1E_0A}=U\n",
    "\\end{align}\n",
    "$$\n",
    "Rearrange that we can see $PL$\n",
    "$$\n",
    "A = \\underbrace{E_0^{-1}}_{P} \\underbrace{(E_1^{-1}E_2^{-1})}_{L}U\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Solving Linear System by Using LU Factorization</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the linear system：\n",
    "$$\n",
    "\\begin{align}\n",
    "3x_1-7x_2 -2x_3+2x_4&=-9\\\\\n",
    "-3x_1+5x_2 +x_3 &=5\\\\\n",
    "6x_1-4x_2 -5x_4&=7\\\\\n",
    "-9x_1+5x_2 -5x_3+12x_4&=11\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "In matrix form:\n",
    "$$\n",
    "\\underbrace{\\left[\\begin{array}{rrrr}\n",
    "3 & -7 & -2 & 2 \\\\\n",
    "-3 & 5 & 1 & 0 \\\\\n",
    "6 & -4 & 0 & -5 \\\\\n",
    "-9 & 5 & -5 & 12\n",
    "\\end{array}\\right]}_{A}\n",
    "\\left[\\begin{array}{r}\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "x_3 \\\\\n",
    "x_4\n",
    "\\end{array}\\right]\n",
    "=\n",
    "\\underbrace{\\left[\\begin{array}{r}\n",
    "-9 \\\\\n",
    "5 \\\\\n",
    "7 \\\\\n",
    "11\n",
    "\\end{array}\\right]}_{b}\n",
    "$$\n",
    "Perform $LU$ decomposition:\n",
    "\n",
    "$$\\underbrace{\\left[\\begin{array}{rrrr}\n",
    "3 & -7 & -2 & 2 \\\\\n",
    "-3 & 5 & 1 & 0 \\\\\n",
    "6 & -4 & 0 & -5 \\\\\n",
    "-9 & 5 & -5 & 12\n",
    "\\end{array}\\right]}_{A}\n",
    "=\\underbrace{\n",
    "\\left[\\begin{array}{rrrr}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "-1 & 1 & 0 & 0 \\\\\n",
    "2 & -5 & 1 & 0 \\\\\n",
    "-3 & 8 & 3 & 1\n",
    "\\end{array}\\right]}_{L}\\underbrace{\\left[\\begin{array}{rrrr}\n",
    "3 & -7 & -2 & 2 \\\\\n",
    "0 & -2 & -1 & 2 \\\\\n",
    "0 & 0 & -1 & 1 \\\\\n",
    "0 & 0 & 0 & -1\n",
    "\\end{array}\\right]}_{U}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace $A$ by $LU$, we get $L(Ux) = b$, now solve this pair of equations\n",
    "\n",
    "$$\n",
    "Ly = b\\\\\n",
    "Ux = y\n",
    "$$\n",
    "\n",
    "Construct augmented matrix $[L|b]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{\n",
    "\\left[\\begin{array}{rrrr}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "-1 & 1 & 0 & 0 \\\\\n",
    "2 & -5 & 1 & 0 \\\\\n",
    "-3 & 8 & 3 & 1\n",
    "\\end{array}\\right]}_{L}\n",
    "\\underbrace{\\left[\\begin{array}{r}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "y_3 \\\\\n",
    "y_4\n",
    "\\end{array}\\right]}_{y}\n",
    "=\n",
    "\\underbrace{\\left[\\begin{array}{r}\n",
    "-9 \\\\\n",
    "5 \\\\\n",
    "7 \\\\\n",
    "11\n",
    "\\end{array}\\right]}_{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{array}{rrrr|r}\n",
    "1 & 0 & 0 & 0 & -9 \\\\\n",
    "-1 & 1 & 0 & 0 & 5\\\\\n",
    "2 & -5 & 1 & 0 & 7\\\\\n",
    "-3 & 8 & 3 & 1 & 11\n",
    "\\end{array}\\right]\\sim\n",
    "\\left[\\begin{array}{rrrr|r}\n",
    "1 & 0 & 0 & 0 & -9 \\\\\n",
    "0 & 1 & 0 & 0 & -4\\\\\n",
    "0 & 0 & 1 & 0 & 5\\\\\n",
    "0 & 0 & 0 & 1 & 1\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{array}{r}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "y_3 \\\\\n",
    "y_4\n",
    "\\end{array}\\right]=\n",
    "\\left[\\begin{array}{r}\n",
    "-9 \\\\\n",
    "-4 \\\\\n",
    "5 \\\\\n",
    "1\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we solve $Ux = y$, to show this in NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3.,  4., -6., -1.])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.array([[3, -7, -2, 2], [0, -2, -1, 2], [0, 0, -1, 1], [0, 0, 0,-1]])\n",
    "y = np.array([-9, -4, 5, 1])\n",
    "np.linalg.solve(U, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the process is correct, this is the solution of the system and we can verify results by invoking ```np.linalg.solve()```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3.,  4., -6., -1.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[3, -7, -2, 2], [-3, 5, 1, 0], [6, -4, 0, -5], [-9, 5, -5, 12]])\n",
    "b = np.array([-9, 5, 7, 11])\n",
    "np.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are the same, $LU$ decomposition works!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Cholesky Factorization</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cholesky decomposition is analogous to taking the square root of a matrix. A common example is the decomposition of a covariance matrix. In this context, the $L$ matrix plays a role similar to that of the correlation matrix.\n",
    "$$\n",
    "\\Sigma = L L^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we have a vector of uncorrelated random variables $\\mathbf{z}$. Each element of $\\mathbf{z}$ is a standard normal random variable with zero mean and unit variance. In mathematical terms, $\\mathbf{z} \\sim \\mathcal{N}(0, I)$, where $I$ is the identity matrix.\n",
    "$$\n",
    "\\mathbf{z}=\\left(\\begin{array}{c}\n",
    "z_1 \\\\\n",
    "z_2 \\\\\n",
    "\\vdots \\\\\n",
    "z_n\n",
    "\\end{array}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can transform non-correlated series $x$ into corelated $z$ \n",
    "$$\n",
    "x = Lz\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand this, we examine the covariance matrix of the transformed vector \\( x \\):\n",
    "$$\n",
    "\\text{Cov}(x) = \\text{Cov}(Lz) = L \\text{Cov}(z) L^T = L I L^T = L L^T = \\Sigma\n",
    "$$\n",
    "This shows that $\\text{Cov}(x)$ has the desired covariance structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see a numerical example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.99424759, 0.80479653, 0.71424917],\n",
       "       [0.80479653, 1.00962669, 0.81794008],\n",
       "       [0.71424917, 0.81794008, 1.02650048]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# use Antithetic Variates to reduce variance\n",
    "n = 1000\n",
    "z = np.random.normal(size=(1000, 3))\n",
    "z2 = -z\n",
    "z = np.concatenate((z, z2), axis=0)\n",
    "\n",
    "cov_matrix = np.array([[1, 0.8, 0.7], [0.8, 1, 0.8], [0.7, 0.8, 1]])\n",
    "L = np.linalg.cholesky(cov_matrix)\n",
    "\n",
    "z = z @ L.T\n",
    "\n",
    "np.cov(z.T)"
   ]
  }
 ],
 "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.12.4"
  },
  "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": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}