06.2.5 Gaussian Elimination - Answer.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Gaussian Elimination"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this notebook, we walk you through the steps that take you from applying Gauss transforms to a matrix to an algorithm that performs these steps.\n",
      "\n",
      "<font color=red> Be sure to make a copy!!!! </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2>Preliminaries</h2>\n",
      "\n",
      "Here is a list of laff routines that you might want to use in this notebook:\n",
      "<ul>\n",
      "<li> <code>laff.invscal( alpha, x )</code> $x := x / \\alpha$ (See note below)\n",
      "<li> <code>laff.axpy( alpha, x, y )</code> $y := \\alpha x + y$\n",
      "<li> <code>laff.ger( alpha, x, y, A )</code> $A := \\alpha x y^T + A$\n",
      "</ul>\n",
      "\n",
      "<code>laff.invscal( alpha, x)</code> was added after the course started. If you're having problems with git and are manually downloading the notebooks, check the wiki page for notebooks for help on getting updated laff routines.\n",
      "\n",
      "<h2> First, let's create a matrix </h2>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import laff\n",
      "import flame\n",
      "\n",
      "# Later this week, you will find out that applying Gaussian Elimination to a matrix\n",
      "# is equivalent to computing a unit lower triangular matrix, L, and upper \n",
      "# triangular matrix, U, such that A = L U.  Here we use that fact to create a \n",
      "# matrix with which to perform Gaussian elimination that doesn't have nasty fractions.\n",
      "\n",
      "L = np.matrix( ' 1, 0, 0. 0;\\\n",
      "                -2, 1, 0, 0;\\\n",
      "                 1,-3, 1, 0;\\\n",
      "                 2, 3,-1, 1' )\n",
      "\n",
      "U = np.matrix( ' 2,-1, 3,-2;\\\n",
      "                 0,-2, 1,-1;\\\n",
      "                 0, 0, 1, 2;\\\n",
      "                 0, 0, 0, 3' )\n",
      "\n",
      "A = L * U\n",
      "\n",
      "print( 'A = ' )\n",
      "print( A )\n",
      "\n",
      "# create a solution vector x\n",
      "\n",
      "x = np.matrix( '-1;\\\n",
      "                 2;\\\n",
      "                 1;\\\n",
      "                -2' )\n",
      "\n",
      "# store the original value of x\n",
      "\n",
      "xold = np.matrix( np.copy( x ) )\n",
      "\n",
      "# create a solution vector b so that A x = b\n",
      "\n",
      "b = A * x\n",
      "print( 'b = ' )\n",
      "print( b )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Solution via Gauss transforms </h2>\n",
      "\n",
      "Let's use a sequence of Gauss transforms to reduce the matrix to an upper triangular matrix.  We then apply the Gauss transforms to the right-hand side $ b $.  Finally, we perform backsubstition.  Well, actually, here we have Python do the work for us."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Step 1: \n",
      "\n",
      "$ \\left( \\begin{array}{c c c c}\n",
      "1 & 0 & 0 & 0 \\\\\n",
      "-\\lambda_{1,0} & 1 & 0 & 0 \\\\\n",
      "-\\lambda_{2,0} & 0 & 1 & 0 \\\\\n",
      "-\\lambda_{3,0} & 0 & 0 & 1 \n",
      "\\end{array} \n",
      "\\right) \n",
      "\\left( \\begin{array}{ c c c c}\n",
      "\\alpha_{0,0} & \\alpha_{0,1} & \\alpha_{0,2} & \\alpha_{0,3} \\\\\n",
      "\\alpha_{1,0} & \\alpha_{1,1} & \\alpha_{1,2} & \\alpha_{1,3} \\\\\n",
      "\\alpha_{2,0} & \\alpha_{2,1} & \\alpha_{2,2} & \\alpha_{2,3} \\\\\n",
      "\\alpha_{3,0} & \\alpha_{3,1} & \\alpha_{3,2} & \\alpha_{3,3} \n",
      "\\end{array} \\right) \n",
      "\\rightarrow \n",
      "\\left( \\begin{array}{ c c c c}\n",
      "\\upsilon_{0,0} & \\upsilon_{0,1} & \\upsilon_{0,2} & \\upsilon_{0,3} \\\\\n",
      "0 & \\alpha_{1,1} & \\alpha_{1,2} & \\alpha_{1,3} \\\\\n",
      "0 & \\alpha_{2,1} & \\alpha_{2,2} & \\alpha_{2,3} \\\\\n",
      "0 & \\alpha_{3,1} & \\alpha_{3,2} & \\alpha_{3,3} \n",
      "\\end{array} \\right) $ \n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L0 = np.matrix( ' 1, 0, 0, 0;\\\n",
      "                  2, 1, 0, 0;\\\n",
      "                 -1, 0, 1, 0;\\\n",
      "                 -2, 0, 0, 1' );\n",
      "A = L0 * A\n",
      "print( A )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Step 2: \n",
      "\n",
      "$ \\left( \\begin{array}{c c c c}\n",
      "1 & 0 & 0 & 0 \\\\\n",
      "0 & 1 & 0 & 0 \\\\\n",
      "0 & -\\lambda_{2,1} & 1 & 0 \\\\\n",
      "0 & -\\lambda_{3,1} & 0 & 1 \n",
      "\\end{array} \n",
      "\\right) \n",
      "\\left( \\begin{array}{ c c c c}\n",
      "\\upsilon_{0,0} & \\upsilon_{0,1} & \\upsilon_{0,2} & \\upsilon_{0,3} \\\\\n",
      "0 & \\alpha_{1,1} & \\alpha_{1,2} & \\alpha_{1,3} \\\\\n",
      "0 & \\alpha_{2,1} & \\alpha_{2,2} & \\alpha_{2,3} \\\\\n",
      "0 & \\alpha_{3,1} & \\alpha_{3,2} & \\alpha_{3,3} \n",
      "\\end{array} \\right)\n",
      "\\rightarrow \n",
      "\\left( \\begin{array}{ c c c c}\n",
      "\\upsilon_{0,0} & \\upsilon_{0,1} & \\upsilon_{0,2} & \\upsilon_{0,3} \\\\\n",
      "0 & \\upsilon_{1,1} & \\upsilon_{1,2} & \\upsilon_{1,3} \\\\\n",
      "0 & 0 & \\alpha_{2,2} & \\alpha_{2,3} \\\\\n",
      "0 & 0 & \\alpha_{3,2} & \\alpha_{3,3} \n",
      "\\end{array} \\right) $ \n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>\n",
      "\n",
      "\n",
      "YOU fill in the \"?\" in L1 given the matrix you just produced from the code block above:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L1 = np.matrix( ' 1, 0, 0, 0;\\\n",
      "                  0, 1, 0, 0;\\\n",
      "                  0, 3, 1, 0;\\\n",
      "                  0,-3, 0, 1' );\n",
      "A = L1 * A\n",
      "print( A )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Step 3: \n",
      "\n",
      "$ \\left( \\begin{array}{c c c c}\n",
      "1 & 0 & 0 & 0 \\\\\n",
      "0 & 1 & 0 & 0 \\\\\n",
      "0 & 0 & 1 & 0 \\\\\n",
      "0 & 0 & -\\lambda_{3,2} & 1 \n",
      "\\end{array} \n",
      "\\right) \n",
      "\\left( \\begin{array}{ c c c c}\n",
      "\\upsilon_{0,0} & \\upsilon_{0,1} & \\upsilon_{0,2} & \\upsilon_{0,3} \\\\\n",
      "0 & \\upsilon_{1,1} & \\upsilon_{1,2} & \\upsilon_{1,3} \\\\\n",
      "0 & 0 & \\alpha_{2,2} & \\alpha_{2,3} \\\\\n",
      "0 & 0 & \\alpha_{3,2} & \\alpha_{3,3} \n",
      "\\end{array} \\right)\n",
      "\\rightarrow \n",
      "\\left( \\begin{array}{ c c c c}\n",
      "\\upsilon_{0,0} & \\upsilon_{0,1} & \\upsilon_{0,2} & \\upsilon_{0,3} \\\\\n",
      "0 & \\upsilon_{1,1} & \\upsilon_{1,2} & \\upsilon_{1,3} \\\\\n",
      "0 & 0 & \\upsilon_{2,2} & \\upsilon_{2,3} \\\\\n",
      "0 & 0 & 0 & \\upsilon_{3,3} \n",
      "\\end{array} \\right) $ \n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>\n",
      "\n",
      "\n",
      "YOU fill in the \"?\" in L2 given the matrix you just produced from the code block above:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L2 = np.matrix( ' 1, 0, 0, 0;\\\n",
      "                  0, 1, 0, 0;\\\n",
      "                  0, 0, 1, 0;\\\n",
      "                  0, 0, 1, 1' );\n",
      "A = L2 * A\n",
      "print( A )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, apply the Gauss transforms to the right-hand side, $b$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b = L0 * b\n",
      "b = L1 * b\n",
      "b = L2 * b\n",
      "\n",
      "print( 'b after forward substitution (application of Gauss transforms):' )\n",
      "print( b )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, perform back substitution\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[ 3 ] = b[ 3 ] / A[ 3,3 ]\n",
      "x[ 2 ] = ( b[ 2 ] - A[ 2,3 ] * x[ 3 ] ) / A[ 2,2 ]\n",
      "x[ 1 ] = ( b[ 1 ] - A[ 1,2 ] * x[ 2 ] - A[ 1,3 ] * x[ 3 ] ) / A[ 1,1 ]\n",
      "x[ 0 ] = ( b[ 0 ] - A[ 0,1 ] * x[ 1 ] - A[ 0,2 ] * x[ 2 ]- A[ 0,3 ] * x[ 3 ] ) / A[ 0,0 ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( 'x = ' )\n",
      "print( x )\n",
      "\n",
      "print( 'x - xold' )\n",
      "print( x - xold )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<code> x - xold </code> should yield a zero vector"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Now, let's implement the Gaussian Elimination routine from 6.2.5 </h2>\n",
      "\n",
      "Here is the algorithm:\n",
      "\n",
      "<img src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/6.2.5_GE.png\" alt=\"Gaussian elimination algorithm\" width=50%>\n",
      "    \n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create the routine\n",
      "<code> GaussianElimination( A ) </code>\n",
      "with the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> for the algorithm given above.\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "Recall from the explanation that a Gauss transform needs to be computed so that\n",
      "<br>\n",
      "$ \n",
      "\\left( \\begin{array}{c | c }\n",
      "1 & 0 \\\\ \\hline\n",
      "-l_{21} & I \n",
      "\\end{array} \\right)\n",
      "\\left( \\begin{array}{c | c }\n",
      "\\alpha_{11} & a_{12}^T \\\\ \\hline\n",
      "a_{21} & A_{22} \n",
      "\\end{array} \\right)\n",
      "= \n",
      "\\left( \\begin{array}{c | c }\n",
      "\\alpha_{11} & a_{12}^T \\\\ \\hline\n",
      "0 & A_{22} - l_{21} a_{12}^T\n",
      "\\end{array} \\right)\n",
      "$.  \n",
      "Notice that it must be true that $ l_{21} := a_{21} / \\alpha_{11} $ in order to get a zero in the bottom left quadrant of the resulting matrix.\n",
      "\n",
      "It follows from partitioned matrix multiplication that we must update $ A_{22} := A_{22} - l_{21} * a_{12}^T $.\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "import laff\n",
      "\n",
      "def GaussianElimination_unb(A):\n",
      "\n",
      "    ATL, ATR, \\\n",
      "    ABL, ABR  = flame.part_2x2(A, \\\n",
      "                               0, 0, 'TL')\n",
      "\n",
      "    while ATL.shape[0] < A.shape[0]:\n",
      "\n",
      "        A00,  a01,     A02,  \\\n",
      "        a10t, alpha11, a12t, \\\n",
      "        A20,  a21,     A22   = flame.repart_2x2_to_3x3(ATL, ATR, \\\n",
      "                                                       ABL, ABR, \\\n",
      "                                                       1, 1, 'BR')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.invscal( alpha11, a21 )        #  a21 := a21 / alpha11\n",
      "        laff.ger( -1.0, a21, a12t, A22 )    #  A22 := A22 - a21 * a12t\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        ATL, ATR, \\\n",
      "        ABL, ABR  = flame.cont_with_3x3_to_2x2(A00,  a01,     A02,  \\\n",
      "                                               a10t, alpha11, a12t, \\\n",
      "                                               A20,  a21,     A22,  \\\n",
      "                                               'TL')\n",
      "\n",
      "    flame.merge_2x2(ATL, ATR, \\\n",
      "                    ABL, ABR, A)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3> Test the routine </h3>\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# recreate matrix A\n",
      "A = L * U\n",
      "\n",
      "# recreate the right-hand side\n",
      "b = A * xold\n",
      "\n",
      "# apply Gaussian elimination to matrix A\n",
      "GaussianElimination_unb( A )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compare the overwritten matrix, $ A $, to the original matrices, $ L $ and $ U $.  The upper triangular part of $ A $ should equal $ U $ and the strictly lower triangular part of $ A $ should equal the strictly lower triangular part of $ L $.  The reason for this will become clear in Unit 6.3.1.\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( A )\n",
      "print( 'Original L' )\n",
      "print( L )\n",
      "print( 'Original U' )\n",
      "print( U )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2> Now, let's implement the forward substitution routine from 6.2.5 </h2>\n",
      "\n",
      "Here is the algorithm:\n",
      "\n",
      "<img src=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/6.2.5_FS.png\" alt=\"Forward substitution algorithm\" width=75%>\n",
      "    \n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Create the routine\n",
      "<code> ForwardSubstitution_unb( A, b ) </code>\n",
      "with the <a href=\"https://studio.edx.org/c4x/UTAustinX/UT.5.01x/asset/index.html\"> Spark webpage</a> for the algorithm."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import flame\n",
      "\n",
      "def ForwardSubstitution_unb(A, b):\n",
      "\n",
      "    ATL, ATR, \\\n",
      "    ABL, ABR  = flame.part_2x2(A, \\\n",
      "                               0, 0, 'TL')\n",
      "\n",
      "    bT, \\\n",
      "    bB  = flame.part_2x1(b, \\\n",
      "                         0, 'TOP')\n",
      "\n",
      "    while ATL.shape[0] < A.shape[0]:\n",
      "\n",
      "        A00,  a01,     A02,  \\\n",
      "        a10t, alpha11, a12t, \\\n",
      "        A20,  a21,     A22   = flame.repart_2x2_to_3x3(ATL, ATR, \\\n",
      "                                                       ABL, ABR, \\\n",
      "                                                       1, 1, 'BR')\n",
      "\n",
      "        b0,    \\\n",
      "        beta1, \\\n",
      "        b2     = flame.repart_2x1_to_3x1(bT, \\\n",
      "                                         bB, \\\n",
      "                                         1, 'BOTTOM')\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        laff.axpy( -beta1, a21, b2 )\n",
      "\n",
      "        #------------------------------------------------------------#\n",
      "\n",
      "        ATL, ATR, \\\n",
      "        ABL, ABR  = flame.cont_with_3x3_to_2x2(A00,  a01,     A02,  \\\n",
      "                                               a10t, alpha11, a12t, \\\n",
      "                                               A20,  a21,     A22,  \\\n",
      "                                               'TL')\n",
      "\n",
      "        bT, \\\n",
      "        bB  = flame.cont_with_3x1_to_2x1(b0,    \\\n",
      "                                         beta1, \\\n",
      "                                         b2,    \\\n",
      "                                         'TOP')\n",
      "\n",
      "    flame.merge_2x1(bT, \\\n",
      "                    bB, b)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Apply the Gauss transforms to the right-hand side\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( A )\n",
      "print( b )\n",
      "ForwardSubstitution_unb( A, b )\n",
      "\n",
      "print( 'updated b' )\n",
      "print( b )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, perform back substitution\n",
      "\n",
      "<font color=red> Important: if you make a mistake, rerun ALL cells above the cell in which you were working before you rerun the one in which you are working. </font>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x[ 3 ] = b[ 3 ] / A[ 3,3 ]\n",
      "x[ 2 ] = ( b[ 2 ] - A[ 2,3 ] * x[ 3 ] ) / A[ 2,2 ]\n",
      "x[ 1 ] = ( b[ 1 ] - A[ 1,2 ] * x[ 2 ] - A[ 1,3 ] * x[ 3 ] ) / A[ 1,1 ]\n",
      "x[ 0 ] = ( b[ 0 ] - A[ 0,1 ] * x[ 1 ] - A[ 0,2 ] * x[ 2 ]- A[ 0,3 ] * x[ 3 ] ) / A[ 0,0 ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print( 'x = ' )\n",
      "print( x )\n",
      "\n",
      "print( 'x - xold' )\n",
      "print( x - xold )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<code> x - xold </code> should yield a zero vector"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}