{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Solving $ A x = b $ via LU factorization and triangular solves" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, you will implement an LU factorization, solve a system with a unit lower triangular matrix, and solve a system with an upper triangular matrix. This notebook culminates in a routine that combines these three steps into a routine that solves $ A x = b $.\n", "\n", " Be sure to make a copy!!!! \n", "\n", "

Preliminaries

\n", "\n", "Here is a list of laff routines that you might want to use in this notebook:\n", "\n", "\n", "laff.invscal( alpha, x) 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", "\n", "

First, let's create a matrix $ A $ and right-hand side $ b $

" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import laff\n", "import flame\n", "\n", "# Applying the LU Factorization to a matrix is the process of\n", "# computing a unit lower triangular matrix, L, and upper \n", "# triangular matrix, U, such that A = L U. To avoid nasty fractions\n", "# in these caclulations, we create our matrix A from two matrices, L and U, \n", "# whose elements we know to be integer valued.\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 )\n", "\n", "# Much later, we are also going to solve A x = y. Here we create that y:\n", "x2 = np.matrix( '1;\\\n", " 2;\\\n", " -2;\\\n", " 1' )\n", "y = A * x2\n", "print( 'y = ' )\n", "print( y )" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Implement the LU factorization routine from 6.3.1

\n", "\n", "Here is the algorithm:\n", "\n", "\"LU\n", " \n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", " LU_unb_var5( A ) \n", "with the Spark webpage for the algorithm given above." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Insert code here" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Test the routine

\n", "\n", " 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. " ] }, { "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", "LU_unb_var5( 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 $. \n", "\n", " 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. \n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "print( 'A after factorization' )\n", "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": [ "

Implement the routine Ltrsv_unb_var1 from 6.3.2

\n", "\n", "(if you have not yet visited Unit 6.3.2, do so now!)\n", "\n", "Here is the algorithm:\n", "\n", "\"Unit\n", " \n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", " Ltrsv_unb_var1( L, b )\n", "with the Spark webpage for the algorithm given above." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Insert code here" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Test Ltrsv_unb_var1

\n", "\n", "Take the output from LU_unb_var5, and use it to solve $ L z = b $, where $ L $ is unit lower triangular with its strictly lower triangular part stored in the strictly lower triangular part of $ A $.\n", "\n", " 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. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "print( A )\n", "print( b )\n", "Ltrsv_unb_var1( A, b )\n", "\n", "print( 'updated b' )\n", "print( b )" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To be able to perform the back substitution before we implement the upper triangular solve routine described in 6.3.3, we do it the hard way here. This allows us to see if the LU factorization followed by the solve with the unit lower triangular system gives the correct intermediate result.\n", "\n", " 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. " ] }, { "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": [ " x - xold should yield a zero vector" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Implement the routine Utrsv_unb_var1 from 6.3.3

\n", "\n", "(if you have not yet visited Unit 6.3.3, do so now!)\n", "\n", "Here is the algorithm:\n", "\n", "\"Upper\n", " \n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", "Utrsv_unb_var1( U, b )\n", "with the Spark webpage for the algorithm given above.\n", "\n", " Hint: Implement $ \\beta_1 := \\beta - u_{12}^T b_2 $ as laff.dots( -u12t, b2, beta1 ) " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Insert code here" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Test Utrsv_unb_var1

\n", "\n", "Take the output from LU_unb_var5 and Ltrsv_unb_var1 and use it to solve $ U x = b $, where $ U $ is upper triangular and stored in the upper triangular part of $ A $.\n", "\n", " 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. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# just to be sure, let's start over. We'll recreate A, x, and b, run all the routines, and\n", "# then compare the updated b to the original vector x.\n", "\n", "A = L * U\n", "b = A * x\n", "\n", "LU_unb_var5( A )\n", "Ltrsv_unb_var1( A, b )\n", "Utrsv_unb_var1( A, b )\n", "\n", "print( 'updated b' )\n", "print( b )\n", "print( 'original x' )\n", "print( x )\n", "print( 'b - x' )\n", "print( b - x )\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In theory, b - x should yield a zero vector..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Implement the routine Solve from 6.3.4

\n", "\n", "(if you have not yet visited Unit 6.3.4, do so now!)\n", "\n", "This time, we do NOT use Spark! What we need to do is write a routine that, when given a matrix $ A $ and right-hand side vector $ b $, solves $ A x = b $, overwriting $ A $ with the LU factorization and overwriting $ b $ with the solution vector $ x $:\n", "\n", "\n", "\n", " 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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the routine\n", " Solve( A, b ) " ] }, { "cell_type": "code", "collapsed": false, "input": [ "def Solve( A, b ):\n", " \n", " # insert appropriate calls to routines you have written here!\n", " \n", " " ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

Test Solve

\n", "\n", " 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. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# just to be sure, let's start over. We'll recreate A, x, and b, run all the routines, and\n", "# then compare the updated b to the original vector x.\n", "\n", "A = L * U\n", "b = A * x\n", "\n", "Solve( A, b )\n", "\n", "print( 'updated b' )\n", "print( b )\n", "print( 'original x' )\n", "print( x )\n", "print( 'b - x' )\n", "print( b - x )\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In theory, b - x should yield a zero vector..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

What if a new right-hand side comes along?

\n", "\n", "What if we are presented with a new right-hand side, call it $ y $, with which we want to solve $ A x = y $, overwriting $ y $ with the solution? (We created such a $ y $ at the top of this notebook.)\n", "Notice that you can take the matrix $A $ that was modified by Solve and use it with Ltrsv_unb_var1 and Utrsv_unb_var1:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# insert appropriate calls here.\n", "\n", "\n", "print( 'y = ' )\n", "print( y )\n", "\n", "print( 'x2 - y =' )\n", "print( x2 - y )" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "x2 - y should evaluate to the zero vector.\n", "\n", "\n", "

Important: you should not refactor $ A $!!!!

" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }