{"worksheets": [{"metadata": {}, "cells": [{"cell_type": "markdown", "source": ["# LU with Partial Pivoting"], "metadata": {}}, {"language": "python", "outputs": [], "collapsed": false, "input": ["import numpy as np\n", "import numpy.linalg as la\n", "\n", "np.set_printoptions(precision=3, suppress=True)"], "metadata": {}, "cell_type": "code", "prompt_number": 131}, {"cell_type": "markdown", "source": ["# Set-up"], "metadata": {}}, {"cell_type": "markdown", "source": ["Let's grab a (admittedly well-chosen) sample matrix `A`:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[  0.,   4.,  19.,  -7.],\n", "       [ -1.,  -2., -10.,  -0.],\n", "       [  1.,  17.,   1.,  -4.],\n", "       [ -5.,  -8.,  -6.,  -2.]])"], "prompt_number": 132, "metadata": {}}], "collapsed": false, "input": ["n = 4\n", "\n", "np.random.seed(235)\n", "A = np.round(5*np.random.randn(n, n))\n", "A[0,0] = 0\n", "A[2,1] = 17\n", "A[0,2] = 19\n", "A"], "metadata": {}, "cell_type": "code", "prompt_number": 132}, {"cell_type": "markdown", "source": ["----------------\n", "## Permutation matrices\n", "\n", "Now define a function `row_swap_mat(i, j)` that returns a permutation matrix that swaps row i and j:"], "metadata": {}}, {"language": "python", "outputs": [], "collapsed": false, "input": ["def row_swap_mat(i, j):\n", "    P = np.eye(n)\n", "    P[i] = 0\n", "    P[j] = 0\n", "    P[i, j] = 1\n", "    P[j, i] = 1\n", "    return P"], "metadata": {}, "cell_type": "code", "prompt_number": 133}, {"cell_type": "markdown", "source": ["What do these matrices look like?"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 0.,  1.,  0.,  0.],\n", "       [ 1.,  0.,  0.,  0.],\n", "       [ 0.,  0.,  1.,  0.],\n", "       [ 0.,  0.,  0.,  1.]])"], "prompt_number": 134, "metadata": {}}], "collapsed": false, "input": ["row_swap_mat(0,1)"], "metadata": {}, "cell_type": "code", "prompt_number": 134}, {"cell_type": "markdown", "source": ["Do they work?"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -1.,  -2., -10.,   0.],\n", "       [  0.,   4.,  19.,  -7.],\n", "       [  1.,  17.,   1.,  -4.],\n", "       [ -5.,  -8.,  -6.,  -2.]])"], "prompt_number": 135, "metadata": {}}], "collapsed": false, "input": ["row_swap_mat(0,1).dot(A)"], "metadata": {}, "cell_type": "code", "prompt_number": 135}, {"cell_type": "markdown", "source": ["--------------"], "metadata": {}}, {"cell_type": "markdown", "source": ["# Part I"], "metadata": {}}, {"cell_type": "markdown", "source": ["`U` is the copy of `A` that we'll modify:"], "metadata": {}}, {"language": "python", "outputs": [], "collapsed": false, "input": ["U = A.copy()"], "metadata": {}, "cell_type": "code", "prompt_number": 136}, {"cell_type": "markdown", "source": ["## First column"], "metadata": {}}, {"cell_type": "markdown", "source": ["Create P1 to swap up the right row: "], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5.,  -8.,  -6.,  -2.],\n", "       [ -1.,  -2., -10.,   0.],\n", "       [  1.,  17.,   1.,  -4.],\n", "       [  0.,   4.,  19.,  -7.]])"], "prompt_number": 137, "metadata": {}}], "collapsed": false, "input": ["P1 = row_swap_mat(0, 3)\n", "U = P1.dot(U)\n", "U"], "metadata": {}, "cell_type": "code", "prompt_number": 137}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 1. ,  0. ,  0. ,  0. ],\n", "       [-0.2,  1. ,  0. ,  0. ],\n", "       [ 0.2,  0. ,  1. ,  0. ],\n", "       [ 0. ,  0. ,  0. ,  1. ]])"], "prompt_number": 138, "metadata": {}}], "collapsed": false, "input": ["M1 = np.eye(n)\n", "M1[1,0] = -U[1,0]/U[0,0]\n", "M1[2,0] = -U[2,0]/U[0,0]\n", "M1"], "metadata": {}, "cell_type": "code", "prompt_number": 138}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5. ,  -8. ,  -6. ,  -2. ],\n", "       [  0. ,  -0.4,  -8.8,   0.4],\n", "       [  0. ,  15.4,  -0.2,  -4.4],\n", "       [  0. ,   4. ,  19. ,  -7. ]])"], "prompt_number": 139, "metadata": {}}], "collapsed": false, "input": ["U = M1.dot(U)\n", "U"], "metadata": {}, "cell_type": "code", "prompt_number": 139}, {"cell_type": "markdown", "source": ["## Second column\n", "\n", "Create `P2` to swap up the right row:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5. ,  -8. ,  -6. ,  -2. ],\n", "       [  0. ,  15.4,  -0.2,  -4.4],\n", "       [  0. ,  -0.4,  -8.8,   0.4],\n", "       [  0. ,   4. ,  19. ,  -7. ]])"], "prompt_number": 140, "metadata": {}}], "collapsed": false, "input": ["P2 = row_swap_mat(2,1)\n", "U = P2.dot(U)\n", "U"], "metadata": {}, "cell_type": "code", "prompt_number": 140}, {"cell_type": "markdown", "source": ["Make the second-column elimination matrix `M2`:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n", "       [ 0.   ,  1.   ,  0.   ,  0.   ],\n", "       [ 0.   ,  0.026,  1.   ,  0.   ],\n", "       [ 0.   , -0.26 ,  0.   ,  1.   ]])"], "prompt_number": 141, "metadata": {}}], "collapsed": false, "input": ["M2 = np.eye(n)\n", "M2[2,1] = -U[2,1]/U[1,1]\n", "M2[3,1] = -U[3,1]/U[1,1]\n", "M2"], "metadata": {}, "cell_type": "code", "prompt_number": 141}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n", "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n", "       [  0.   ,   0.   ,  -8.805,   0.286],\n", "       [  0.   ,   0.   ,  19.052,  -5.857]])"], "prompt_number": 142, "metadata": {}}], "collapsed": false, "input": ["U = M2.dot(U)\n", "U"], "metadata": {}, "cell_type": "code", "prompt_number": 142}, {"cell_type": "markdown", "source": ["## Third column\n", "\n", "Create `P3` to swap up the right entry:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n", "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n", "       [  0.   ,   0.   ,  19.052,  -5.857],\n", "       [  0.   ,   0.   ,  -8.805,   0.286]])"], "prompt_number": 143, "metadata": {}}], "collapsed": false, "input": ["P3 = row_swap_mat(3, 2)\n", "U = P3.dot(U)\n", "U"], "metadata": {}, "cell_type": "code", "prompt_number": 143}, {"cell_type": "markdown", "source": ["Make the third-column elimination matrix `M3`:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n", "       [ 0.   ,  1.   ,  0.   ,  0.   ],\n", "       [ 0.   ,  0.   ,  1.   ,  0.   ],\n", "       [ 0.   ,  0.   ,  0.462,  1.   ]])"], "prompt_number": 144, "metadata": {}}], "collapsed": false, "input": ["M3 = np.eye(n)\n", "M3[3,2] = -U[3,2]/U[2,2]\n", "M3"], "metadata": {}, "cell_type": "code", "prompt_number": 144}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n", "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n", "       [  0.   ,   0.   ,  19.052,  -5.857],\n", "       [  0.   ,   0.   ,   0.   ,  -2.421]])"], "prompt_number": 145, "metadata": {}}], "collapsed": false, "input": ["U = M3.dot(U)\n", "U"], "metadata": {}, "cell_type": "code", "prompt_number": 145}, {"cell_type": "markdown", "source": ["## Wrap-up"], "metadata": {}}, {"cell_type": "markdown", "source": ["So we've built $M3P_3M_2P_2M_1P_1A=U$."], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n", "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n", "       [  0.   ,   0.   ,  19.052,  -5.857],\n", "       [  0.   ,   0.   ,   0.   ,  -2.421]])"], "prompt_number": 150, "metadata": {}}], "collapsed": false, "input": ["M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1).dot(A)"], "metadata": {}, "cell_type": "code", "prompt_number": 150}, {"cell_type": "markdown", "source": ["---------------------\n", "That left factor is anything but lower triangular:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 0.   ,  0.   ,  0.   ,  1.   ],\n", "       [ 0.   ,  0.   ,  1.   ,  0.2  ],\n", "       [ 1.   ,  0.   , -0.26 , -0.052],\n", "       [ 0.462,  1.   , -0.094, -0.219]])"], "prompt_number": 151, "metadata": {}}], "collapsed": false, "input": ["M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1)"], "metadata": {}, "cell_type": "code", "prompt_number": 151}, {"cell_type": "markdown", "source": ["# Part II"], "metadata": {}}, {"cell_type": "markdown", "source": ["Now try the reordering trick:"], "metadata": {}}, {"language": "python", "outputs": [], "collapsed": false, "input": ["L3 = M3\n", "L2 = P3.dot(M2).dot(la.inv(P3))\n", "L1 = P3.dot(P2).dot(M1).dot(la.inv(P2)).dot(la.inv(P3))"], "metadata": {}, "cell_type": "code", "prompt_number": 160}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 0.   ,  0.   ,  0.   ,  1.   ],\n", "       [ 0.   ,  0.026,  1.   ,  0.   ],\n", "       [ 1.   , -0.266, -0.26 ,  0.   ],\n", "       [ 0.462,  0.878, -0.094,  0.   ]])"], "prompt_number": 155, "metadata": {}}], "collapsed": false, "input": ["L3.dot(L2).dot(L1).dot(P3).dot(P2).dot(P1)"], "metadata": {}, "cell_type": "code", "prompt_number": 155}, {"cell_type": "markdown", "source": ["--------------\n", "We were promised that all of the `L`*n* are still lower-triangular:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "stream", "text": ["[[ 1.   0.   0.   0. ]\n", " [ 0.2  1.   0.   0. ]\n", " [ 0.   0.   1.   0. ]\n", " [-0.2  0.   0.   1. ]]\n", "[[ 1.     0.     0.     0.   ]\n", " [ 0.     1.     0.     0.   ]\n", " [ 0.    -0.26   1.     0.   ]\n", " [ 0.     0.026  0.     1.   ]]\n", "[[ 1.     0.     0.     0.   ]\n", " [ 0.     1.     0.     0.   ]\n", " [ 0.     0.     1.     0.   ]\n", " [ 0.     0.     0.462  1.   ]]\n"], "stream": "stdout"}], "collapsed": false, "input": ["print(L1)\n", "print(L2)\n", "print(L3)"], "metadata": {}, "cell_type": "code", "prompt_number": 168}, {"cell_type": "markdown", "source": ["So their product is, too:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n", "       [ 0.2  ,  1.   ,  0.   ,  0.   ],\n", "       [-0.052, -0.26 ,  1.   ,  0.   ],\n", "       [-0.219, -0.094,  0.462,  1.   ]])"], "prompt_number": 172, "metadata": {}}], "collapsed": false, "input": ["Ltemp = L3.dot(L2).dot(L1)\n", "Ltemp"], "metadata": {}, "cell_type": "code", "prompt_number": 172}, {"cell_type": "markdown", "source": ["----\n", "`P` is still a permutation matrix (but a more complicated one):"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 0.,  0.,  0.,  1.],\n", "       [ 0.,  0.,  1.,  0.],\n", "       [ 1.,  0.,  0.,  0.],\n", "       [ 0.,  1.,  0.,  0.]])"], "prompt_number": 174, "metadata": {}}], "collapsed": false, "input": ["P = P3.dot(P2).dot(P1)\n", "P"], "metadata": {}, "cell_type": "code", "prompt_number": 174}, {"cell_type": "markdown", "source": ["-----------------\n", "So to sum up, we've made:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 0.,  0.,  0.,  0.],\n", "       [ 0.,  0.,  0.,  0.],\n", "       [ 0.,  0.,  0.,  0.],\n", "       [ 0.,  0.,  0.,  0.]])"], "prompt_number": 175, "metadata": {}}], "collapsed": false, "input": ["Ltemp.dot(P).dot(A) - U"], "metadata": {}, "cell_type": "code", "prompt_number": 175}, {"cell_type": "markdown", "source": ["--------------\n", "Multiply from the left by `Ltemp`${}^{-1}$, which is *also* lower triangular:"], "metadata": {}}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n", "       [-0.2  ,  1.   ,  0.   ,  0.   ],\n", "       [ 0.   ,  0.26 ,  1.   ,  0.   ],\n", "       [ 0.2  , -0.026, -0.462,  1.   ]])"], "prompt_number": 179, "metadata": {}}], "collapsed": false, "input": ["L = la.inv(Ltemp)\n", "L"], "metadata": {}, "cell_type": "code", "prompt_number": 179}, {"language": "python", "outputs": [{"output_type": "pyout", "text": ["array([[ 0.,  0.,  0.,  0.],\n", "       [ 0.,  0.,  0.,  0.],\n", "       [ 0.,  0.,  0.,  0.],\n", "       [ 0., -0.,  0.,  0.]])"], "prompt_number": 180, "metadata": {}}], "collapsed": false, "input": ["P.dot(A) - L.dot(U)"], "metadata": {}, "cell_type": "code", "prompt_number": 180}]}], "nbformat": 3, "metadata": {"name": "", "signature": "sha256:f365635128cf850308f0b940b5eff4f24875499c9a8d133fb1671c21e50c6f85"}, "nbformat_minor": 0}