{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# LU with Partial Pivoting" ] }, { "cell_type": "code", "execution_count": 131, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import numpy.linalg as la\n", "\n", "np.set_printoptions(precision=3, suppress=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Set-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's grab a (admittedly well-chosen) sample matrix `A`:" ] }, { "cell_type": "code", "execution_count": 132, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0., 4., 19., -7.],\n", " [ -1., -2., -10., -0.],\n", " [ 1., 17., 1., -4.],\n", " [ -5., -8., -6., -2.]])" ] }, "execution_count": 132, "metadata": {}, "output_type": "execute_result" } ], "source": [ "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" ] }, { "cell_type": "markdown", "metadata": {}, "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:" ] }, { "cell_type": "code", "execution_count": 133, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What do these matrices look like?" ] }, { "cell_type": "code", "execution_count": 134, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0., 1., 0., 0.],\n", " [ 1., 0., 0., 0.],\n", " [ 0., 0., 1., 0.],\n", " [ 0., 0., 0., 1.]])" ] }, "execution_count": 134, "metadata": {}, "output_type": "execute_result" } ], "source": [ "row_swap_mat(0,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do they work?" ] }, { "cell_type": "code", "execution_count": 135, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ -1., -2., -10., 0.],\n", " [ 0., 4., 19., -7.],\n", " [ 1., 17., 1., -4.],\n", " [ -5., -8., -6., -2.]])" ] }, "execution_count": 135, "metadata": {}, "output_type": "execute_result" } ], "source": [ "row_swap_mat(0,1).dot(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "--------------" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part I" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`U` is the copy of `A` that we'll modify:" ] }, { "cell_type": "code", "execution_count": 136, "metadata": { "collapsed": false }, "outputs": [], "source": [ "U = A.copy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## First column" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create P1 to swap up the right row: " ] }, { "cell_type": "code", "execution_count": 137, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ -5., -8., -6., -2.],\n", " [ -1., -2., -10., 0.],\n", " [ 1., 17., 1., -4.],\n", " [ 0., 4., 19., -7.]])" ] }, "execution_count": 137, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P1 = row_swap_mat(0, 3)\n", "U = P1.dot(U)\n", "U" ] }, { "cell_type": "code", "execution_count": 138, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [-0.2, 1. , 0. , 0. ],\n", " [ 0.2, 0. , 1. , 0. ],\n", " [ 0. , 0. , 0. , 1. ]])" ] }, "execution_count": 138, "metadata": {}, "output_type": "execute_result" } ], "source": [ "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" ] }, { "cell_type": "code", "execution_count": 139, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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. ]])" ] }, "execution_count": 139, "metadata": {}, "output_type": "execute_result" } ], "source": [ "U = M1.dot(U)\n", "U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Second column\n", "\n", "Create `P2` to swap up the right row:" ] }, { "cell_type": "code", "execution_count": 140, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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. ]])" ] }, "execution_count": 140, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P2 = row_swap_mat(2,1)\n", "U = P2.dot(U)\n", "U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make the second-column elimination matrix `M2`:" ] }, { "cell_type": "code", "execution_count": 141, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [ 0. , 1. , 0. , 0. ],\n", " [ 0. , 0.026, 1. , 0. ],\n", " [ 0. , -0.26 , 0. , 1. ]])" ] }, "execution_count": 141, "metadata": {}, "output_type": "execute_result" } ], "source": [ "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" ] }, { "cell_type": "code", "execution_count": 142, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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]])" ] }, "execution_count": 142, "metadata": {}, "output_type": "execute_result" } ], "source": [ "U = M2.dot(U)\n", "U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Third column\n", "\n", "Create `P3` to swap up the right entry:" ] }, { "cell_type": "code", "execution_count": 143, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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]])" ] }, "execution_count": 143, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P3 = row_swap_mat(3, 2)\n", "U = P3.dot(U)\n", "U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make the third-column elimination matrix `M3`:" ] }, { "cell_type": "code", "execution_count": 144, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0. , 0. , 0. ],\n", " [ 0. , 1. , 0. , 0. ],\n", " [ 0. , 0. , 1. , 0. ],\n", " [ 0. , 0. , 0.462, 1. ]])" ] }, "execution_count": 144, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M3 = np.eye(n)\n", "M3[3,2] = -U[3,2]/U[2,2]\n", "M3" ] }, { "cell_type": "code", "execution_count": 145, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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]])" ] }, "execution_count": 145, "metadata": {}, "output_type": "execute_result" } ], "source": [ "U = M3.dot(U)\n", "U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wrap-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we've built $M3P_3M_2P_2M_1P_1A=U$." ] }, { "cell_type": "code", "execution_count": 150, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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]])" ] }, "execution_count": 150, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1).dot(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---------------------\n", "That left factor is anything but lower triangular:" ] }, { "cell_type": "code", "execution_count": 151, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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]])" ] }, "execution_count": 151, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part II" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now try the reordering trick:" ] }, { "cell_type": "code", "execution_count": 160, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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))" ] }, { "cell_type": "code", "execution_count": 155, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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. ]])" ] }, "execution_count": 155, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L3.dot(L2).dot(L1).dot(P3).dot(P2).dot(P1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "--------------\n", "We were promised that all of the `L`*n* are still lower-triangular:" ] }, { "cell_type": "code", "execution_count": 168, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "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" ] } ], "source": [ "print(L1)\n", "print(L2)\n", "print(L3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So their product is, too:" ] }, { "cell_type": "code", "execution_count": 172, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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. ]])" ] }, "execution_count": 172, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Ltemp = L3.dot(L2).dot(L1)\n", "Ltemp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "`P` is still a permutation matrix (but a more complicated one):" ] }, { "cell_type": "code", "execution_count": 174, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0., 0., 0., 1.],\n", " [ 0., 0., 1., 0.],\n", " [ 1., 0., 0., 0.],\n", " [ 0., 1., 0., 0.]])" ] }, "execution_count": 174, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P = P3.dot(P2).dot(P1)\n", "P" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----------------\n", "So to sum up, we've made:" ] }, { "cell_type": "code", "execution_count": 175, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0., 0., 0., 0.],\n", " [ 0., 0., 0., 0.],\n", " [ 0., 0., 0., 0.],\n", " [ 0., 0., 0., 0.]])" ] }, "execution_count": 175, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Ltemp.dot(P).dot(A) - U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "--------------\n", "Multiply from the left by `Ltemp`${}^{-1}$, which is *also* lower triangular:" ] }, { "cell_type": "code", "execution_count": 179, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "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. ]])" ] }, "execution_count": 179, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L = la.inv(Ltemp)\n", "L" ] }, { "cell_type": "code", "execution_count": 180, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0., 0., 0., 0.],\n", " [ 0., 0., 0., 0.],\n", " [ 0., 0., 0., 0.],\n", " [ 0., -0., 0., 0.]])" ] }, "execution_count": 180, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P.dot(A) - L.dot(U)" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 0 }