{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Gaussian elimination with elimination matrices"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 56,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 57,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[-2.,  2., -1.],\n",
              "       [-3.,  1., -9.],\n",
              "       [-5., -5., -2.]])"
            ]
          },
          "execution_count": 57,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "n = 3\n",
        "\n",
        "np.random.seed(15)\n",
        "A = np.round(5*np.random.randn(n, n))\n",
        "\n",
        "A"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`U` is the copy of `A` that we'll modify:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 58,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "U = A.copy()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now eliminate `U[1,0]`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 59,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1. ,  0. ,  0. ],\n",
              "       [-1.5,  1. ,  0. ],\n",
              "       [ 0. ,  0. ,  1. ]])"
            ]
          },
          "execution_count": 59,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M1 = np.eye(n)\n",
        "M1[1,0] = -U[1,0]/U[0,0]\n",
        "M1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 60,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[-2. ,  2. , -1. ],\n",
              "       [ 0. , -2. , -7.5],\n",
              "       [-5. , -5. , -2. ]])"
            ]
          },
          "execution_count": 60,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U = M1.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now eliminate `U[2,0]`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 61,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "M2 = np.eye(n)\n",
        "M2[2,0] = -U[2,0]/U[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 62,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -2. ,   2. ,  -1. ],\n",
              "       [  0. ,  -2. ,  -7.5],\n",
              "       [  0. , -10. ,   0.5]])"
            ]
          },
          "execution_count": 62,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U = np.dot(M2, U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now eliminate `U[2,1]`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 63,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "M3 = np.eye(n)\n",
        "M3[2,1] = -U[2,1]/U[1,1]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 64,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -2. ,   2. ,  -1. ],\n",
              "       [  0. ,  -2. ,  -7.5],\n",
              "       [  0. ,   0. ,  38. ]])"
            ]
          },
          "execution_count": 64,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U = M3.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "------------------\n",
        "\n",
        "Try inverting one of the `M`s:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 65,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[ 1.   0.   0. ]\n",
            " [ 0.   1.   0. ]\n",
            " [-2.5  0.   1. ]]\n",
            "[[ 1.  -0.  -0. ]\n",
            " [ 0.   1.   0. ]\n",
            " [ 2.5  0.   1. ]]\n"
          ]
        }
      ],
      "source": [
        "print(M2)\n",
        "print(la.inv(M2))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "----------------\n",
        "\n",
        "So we've built `M3*M2*M1*A=U`. Test:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 66,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -2. ,   2. ,  -1. ],\n",
              "       [  0. ,  -2. ,  -7.5],\n",
              "       [  0. ,   0. ,  38. ]])"
            ]
          },
          "execution_count": 66,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U2 = M3.dot(M2.dot(M1.dot(A)))\n",
        "U2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 67,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -2. ,   2. ,  -1. ],\n",
              "       [  0. ,  -2. ,  -7.5],\n",
              "       [  0. ,   0. ,  38. ]])"
            ]
          },
          "execution_count": 67,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now define `L`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 68,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1. ,  0. ,  0. ],\n",
              "       [ 1.5,  1. ,  0. ],\n",
              "       [ 2.5,  5. ,  1. ]])"
            ]
          },
          "execution_count": 68,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "L = la.inv(M1).dot(la.inv(M2).dot(la.inv(M3)))\n",
        "L"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Observations? (Shape? Diagonal values?)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 69,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.]])"
            ]
          },
          "execution_count": 69,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "np.dot(L, U) - A"
      ]
    }
  ],
  "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 0
}