{"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}