{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Methods\n",
"\n",
"# Lecture 6: Numerical Linear Algebra II"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.linalg as sl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Learning objectives:\n",
"\n",
"* More on direct methods: LU decomposition\n",
"* Doolittle's algorithm\n",
"* Properties of lower-triangular matrices\n",
"* Partial pivoting"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"Last week we developed and implemented the Gaussian elimination method to solve the linear matrix system ($A\\pmb{x}=\\pmb{b}$). \n",
"\n",
"This week we will consider a closely related solution method: *LU decomposition* or *LU factorisation*.\n",
"\n",
"Both are examples of *direct* solution methods - in the next lecture we will consider the alternate approach to solve linear systems, namely *iterative* or *indirect* methods."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## LU decomposition - motivation\n",
"\n",
"Last week we implemented Gaussian elimination to solve the matrix system with *one* RHS vector $\\pmb{b}$. \n",
"\n",
"We often have to deal with problems where we have multiple RHS vectors, all with the same matrix $A$.\n",
"\n",
"We could call the same code (e.g. our 'upper_triangle' and 'back_substitution' codes from last lecture) multiple times to solve all of these corresponding linear systems, but note that as the elimination algorithm is actually performing operations based upon (the same) $A$ each time, we would actually be wasting time repeating exactly the same operations - this is therefore clearly not an efficient solution to this problem.\n",
"\n",
"We could easily generalise our Gaussian elimination/back substitution algorithms to include multiple RHS column vectors in the augmented system, perform the same sequence of row operations (but now only once) to transform the matrix to upper-triangular form, and then perform back substitution on each of the transformed RHS vectors from the augmented system - cf. the use of Gaussian elimination to compute the inverse to a matrix by placing the identity on the right of the augmented system.\n",
"\n",
"However, it is often the case that each RHS vector depends on the solutions to the matrix systems obtained from some or all of the earlier RHS vectors, and so this generalisation would not work in this case. For example, our discrete system could be of the form\n",
"\n",
"$$A\\pmb{x}^{n+1} = \\pmb{b}(\\pmb{x}^n, \\ldots)$$\n",
"\n",
"where $\\pmb{x}^{n+1}$ is the numerical solution of a (ordinary or partial) differential equation at time level $n+1$, the RHS is a function of the the solution at the previous time level $n$ (and possibly other things such as focing terms, represented by the $\\ldots$ above. $A$ is then a discretisation of the differential equation written in the form that it maps the old solution to the new one. Time stepping this problem invovles solving the same linear system multiple times with different RHSs. \n",
"\n",
"[NB. this is true if our model is *linear*, in the nonlinear case $A$ would be updated evey time step]."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## LU decomposition - theory\n",
"\n",
"To deal with this situation efficiently we *decompose* or *factorise* the matrix $A$ in such a way that it is cheap to compute a new solution vector $\\pmb{x}$ for any given RHS vector $\\pmb{b}$. This decomposition involves a lower- and an upper-triangular matrix, hence the name LU decomposition. These matrices essentially *encode* the steps conducted in Gaussian elimination, so we don't have to explicilty conduct all of the operations again and again.\n",
"\n",
"Mathematically, let's assume that we have already found/constructed a lower-triangular matrix ($L$ - where all entries above the diagonal are zero) and an upper-triangular matrix ($U$ - where all entries below the diagonal are zero) such that we can write\n",
"\n",
"$$ A = LU $$\n",
"\n",
"In this case the matrix system we need to solve for $\\pmb{x}$ becomes\n",
"\n",
"$$ A\\pmb{x} = \\pmb{b} \\iff (LU)\\pmb{x} = L(U\\pmb{x}) = \\pmb{b} $$\n",
"\n",
"Notice that the matrix-vector product $U\\pmb{x}$ is itself a vector, let's call it $\\pmb{c}$ for the time-being (i.e. \n",
"$\\pmb{c}=U\\pmb{x}$).\n",
"\n",
"The above system then reads \n",
"\n",
"$$ L\\pmb{c} = \\pmb{b} $$\n",
"\n",
"where $L$ is a matrix and $\\pmb{c}$ is an unknown. \n",
"\n",
"As $L$ is in lower-triangular form we can use forward substitution (generalising the back subsitution algorithm/code we developed last week) to very easily find $\\pmb{c}$ in relatively few operations (we don't need to call the entire Gaussian elimination algorithm).\n",
"\n",
"Once we know $\\pmb{c}$ we then solve the second linear system \n",
"\n",
"$$ U\\pmb{x} = \\pmb{c} $$\n",
"\n",
"where now we can use the fact that $U$ is upper-triangular to use our back substitution algorithm again very efficiently to give the solution $\\pmb{x}$ we require.\n",
"\n",
"So for a given $\\pmb{b}$ we can find the corresponding $\\pmb{x}$ very efficiently, we can therefore do this repeatedly as each new $\\pmb{b}$ is given to us.\n",
"\n",
"Our challenge is therefore to find the matrices $L$ and $U$ allowing us to perform the decomposition $A=LU$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## LU decomposition - algorithm\n",
"\n",
"Recall the comment above on the $L$ and $U$ matrices encoding the steps taken in Gaussian elimination. Let's see how this works through the development of the so-called *Doolittle algorithm*.\n",
"\n",
"Let's consider an example matrix:\n",
"\n",
"$$\n",
" A=\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
"{\\color{black}5} & {\\color{black}14} & {\\color{black}7} & {\\color{black}10}\\\\\n",
"{\\color{black}20} & {\\color{black}77} & {\\color{black}41} & {\\color{black}48}\\\\\n",
"{\\color{black}25} & {\\color{black}91} & {\\color{black}55} & {\\color{black}67}\\\\\n",
" \\end{bmatrix}\n",
"$$\n",
"\n",
"the first step of Gaussian elimination is to set the\n",
"sub-diagonal elements in the first column to zero by subtracting multiples of\n",
"the first row from each of the subsequent rows. \n",
"\n",
"For this example, using the symbolic notiation from last week\n",
"this requires the row operations\n",
"\n",
"\\begin{align}\n",
"Eq. (2) &\\leftarrow Eq. (2) - 1\\times Eq. (1)\\\\\n",
"Eq. (3) &\\leftarrow Eq. (3) - 4\\times Eq. (1)\\\\\n",
"Eq. (4) &\\leftarrow Eq. (4) - 5\\times Eq. (1)\\\\\n",
"\\end{align}\n",
"\n",
"or mathematically, and for each element of the matrix (remembering that we are adding rows together - while one of\n",
"the entries of a row will end up being zero, this also has the consequence of updating the rest of the values in that row, hence the iteration over $j$ below):\n",
"\n",
"\\begin{align}\n",
"A_{2j} &\\leftarrow A_{2j} - \\frac{A_{21}}{A_{11}} A_{1j} = A_{2j} - \\frac{5}{5} \\times A_{1j}, \\quad j=1,2,3,4\\\\\n",
"A_{3j} &\\leftarrow A_{3j} - \\frac{A_{31}}{A_{11}} A_{1j} = A_{3j} - \\frac{20}{5} \\times A_{1j}, \\quad j=1,2,3,4\\\\\n",
"A_{4j} &\\leftarrow A_{4j} - \\frac{A_{41}}{A_{11}} A_{1j} = A_{4j} - \\frac{25}{5} \\times A_{1j}, \\quad j=1,2,3,4\\\\\n",
"\\end{align}\n",
"\n",
"Notice that we can also write these exact operations on elements in terms of multiplication by a carefully chosen lower-triangular matrix where the non-zero's below the diagonal are restricted to a single column, e.g. for the example above\n",
"\n",
"$$\n",
" \\begin{bmatrix}\n",
" {\\color{black}1} & {\\color{black}0} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{Orange}{-1}} & {\\color{black}1} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{Orange}{-4}} & {\\color{black}0} & {\\color{black}1} & {\\color{black}0}\\\\\n",
" {\\color{Orange}{-5}} & {\\color{black}0} & {\\color{black}0} & {\\color{black}1}\\\\ \n",
" \\end{bmatrix}\\qquad\\times\\qquad\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{black}5} & {\\color{black}14} & {\\color{black}7} & {\\color{black}10}\\\\\n",
" {\\color{black}20} & {\\color{black}77} & {\\color{black}41} & {\\color{black}48}\\\\\n",
" {\\color{black}25} & {\\color{black}91} & {\\color{black}55} & {\\color{black}67}\\\\ \n",
" \\end{bmatrix}\\qquad=\\qquad\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{blue}{0}} & {\\color{blue}{7}} & {\\color{blue}{2}} & {\\color{blue}{1}}\\\\\n",
" {\\color{blue}{0}} & {\\color{blue}{49}} & {\\color{blue}{21}} & {\\color{blue}{12}}\\\\\n",
" {\\color{blue}{0}} & {\\color{blue}{56}} & {\\color{blue}{30}} & {\\color{blue}{22}}\\\\ \n",
" \\end{bmatrix}\n",
"$$\n",
"\n",
"The lower-triangular matrix (let's call this one $L_0$) is thus encoding the first step in Gaussian elimination.\n",
"\n",
"The next step involves taking the second row of the updated matrix as the new pivot (we will ignore partial pivoting for simplicity), and subtracting multiples of this row from those below in order to set the zeros below the diagonal in the second column to zero. \n",
"\n",
"Note that thios step can be achieved here with the multiplication by the following lower-triangular matrix (call this one $L_1$)\n",
"\n",
"\\begin{equation*}\n",
" \\begin{bmatrix}\n",
" {\\color{black}1} & {\\color{black}0} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}0} & {\\color{black}1} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}0} & {\\color{Orange}{-7}} & {\\color{black}1} & {\\color{black}0}\\\\\n",
" {\\color{black}0} & {\\color{Orange}{-8}} & {\\color{black}0} & {\\color{black}1}\\\\\n",
" \\end{bmatrix}\\qquad\\times\\qquad\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{black}0} & {\\color{black}7} & {\\color{black}2} & {\\color{black}1}\\\\\n",
" {\\color{black}0} & {\\color{black}49} & {\\color{black}21} & {\\color{black}12}\\\\\n",
" {\\color{black}0} & {\\color{black}56} & {\\color{black}30} & {\\color{black}22}\\\\\n",
" \\end{bmatrix}\\qquad=\\qquad\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{black}0} & {\\color{black}7} & {\\color{black}2} & {\\color{black}1}\\\\\n",
" {\\color{black}0} & {\\color{blue}{0}} & {\\color{blue}{7}} & {\\color{blue}{5}}\\\\\n",
" {\\color{black}0} & {\\color{blue}{0}} & {\\color{blue}{14}} & {\\color{blue}{14}}\\\\\n",
" \\end{bmatrix}\n",
"\\end{equation*}\n",
"\n",
"\n",
"and finally for this example to get rid the of the leading 14 on the final row:\n",
"\n",
"\\begin{equation*}\n",
" \\begin{bmatrix}\n",
" {\\color{black}1} & {\\color{black}0} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}0} & {\\color{black}1} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}0} & {\\color{black}0} & {\\color{black}1} & {\\color{black}0}\\\\\n",
" {\\color{black}0} & {\\color{black}0} & {\\color{Orange}{-2}} & {\\color{black}{1}}\\\\\n",
" \\end{bmatrix}\\qquad\\times\\qquad\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{black}0} & {\\color{black}7} & {\\color{black}2} & {\\color{black}1}\\\\\n",
" {\\color{black}0} & {\\color{black}0} & {\\color{black}7} & {\\color{black}5}\\\\\n",
" {\\color{black}0} & {\\color{black}0} & {\\color{black}14} & {\\color{black}14}\\\\\n",
" \\end{bmatrix}\\qquad=\\qquad\\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{black}0} & {\\color{black}7} & {\\color{black}2} & {\\color{black}1}\\\\\n",
" {\\color{black}0} & {\\color{black}0} & {\\color{black}7} & {\\color{black}5}\\\\\n",
" {\\color{black}0} & {\\color{black}0} & {\\color{blue}{0}} & {\\color{blue}{4}}\\\\\n",
" \\end{bmatrix}\n",
"\\end{equation*}\n",
"\n",
"where this lower triangular matrix we will call $L_2$, and the RHS matrix is now in upper-triangular form as we expect from Gaussian elimination (call this $U$).\n",
"\n",
"In summary, the above operations (starting from $A$, first multipling pn the left by $L_0$, the multiplying the result on the left by $L_1$ and so on to arrive at an upper trianglualr matrix $U$) can be written as \n",
"\n",
"$$ L_2(L_1(L_0A)) = U $$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that these special lower-triangular matrices $L_0, \\ldots$ are all examples of what is known as *atomic* lower-triangular matrices: a special form of unitriangular matrix - the diagonals are unity, where the off-diagonal entries are all zero apart from in a single column. \n",
"\n",
"Notice that for these special, simple matrices, their inverse is simply the original with the sign of those off-diagnonals changed:\n",
"\n",
"$$\n",
"\\left[\n",
" \\begin{array}{rrrrrrrrr}\n",
" 1 & 0 & \\cdots & & & & & & 0 \\\\\n",
" 0 & 1 & 0 & \\cdots & & & & & 0 \\\\\n",
" 0 & \\ddots & \\ddots & \\ddots & & & & & \\vdots \\\\\n",
" \\vdots & \\ddots & \\ddots & \\ddots & & & & & \\\\\n",
" & & & 0 & 1 & 0 & & & & \\\\ \n",
" & & & 0 & l_{i+1,i} & 1 & \\ddots & & & \\\\ \n",
" & & & 0 & l_{i+2,i} & 0 & \\ddots & & & \\\\ \n",
" & & & \\vdots & \\vdots & \\vdots & \\ddots & & 0 & \\\\ \n",
" 0 & \\cdots & & 0 & l_{n,i} & 0 & \\cdots & 0 & 1 & \\\\ \n",
"\\end{array}\n",
"\\right]^{-1}\n",
"=\n",
"\\left[\n",
" \\begin{array}{rrrrrrrrr}\n",
" 1 & 0 & \\cdots & & & & & & 0 \\\\\n",
" 0 & 1 & 0 & \\cdots & & & & & 0 \\\\\n",
" 0 & \\ddots & \\ddots & \\ddots & & & & & \\vdots \\\\\n",
" \\vdots & \\ddots & \\ddots & \\ddots & & & & & \\\\\n",
" & & & 0 & 1 & 0 & & & & \\\\ \n",
" & & & 0 & -l_{i+1,i} & 1 & \\ddots & & & \\\\ \n",
" & & & 0 & -l_{i+2,i} & 0 & \\ddots & & & \\\\ \n",
" & & & \\vdots & \\vdots & \\vdots & \\ddots & & 0 & \\\\ \n",
" 0 & \\cdots & & 0 & -l_{n,i} & 0 & \\cdots & 0 & 1 & \\\\ \n",
"\\end{array}\n",
"\\right]\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6.1: lower-triangular matrices\n",
"\n",
"Convince yourselves of the following facts:\n",
"\n",
"* The above statement on what the inverse of the special $L$ matrices is - check on a small example.\n",
"\n",
"\n",
"* The multiplication of arbitrary lower-triangular square matrices is also lower-triangular.\n",
"\n",
"\n",
"* $L_2(L_1(L_0A)) = U \\implies A = L_0^{-1}(L_1^{-1}(L_2^{-1}U))$\n",
"\n",
"\n",
"* and hence that $A=LU$ where $U$ is the upper-triangular matrix found at the end of Guassian elimination, and where $L$ is the \n",
"following matrix\n",
"$$ L = L_0^{-1}L_1^{-1}L_2^{-1} $$\n",
"\n",
"\n",
"* Finally, compute this product of these lower-triangular matrices to show that \n",
"$$L = \n",
" \\begin{bmatrix}\n",
" {\\color{black}1} & {\\color{black}0} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}{1}} & {\\color{black}1} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}{4}} & {\\color{black}7} & {\\color{black}1} & {\\color{black}0}\\\\\n",
" {\\color{black}{5}} & {\\color{black}8} & {\\color{black}2} & {\\color{black}1}\\\\ \n",
" \\end{bmatrix}\n",
"$$\n",
"i.e. that the multiplication of these individual atomic matrices (importantly in this order) simply merges the entries from the non-zero columns of each atomic matrix, and hence is both lower-triangular, as well as trivial to compute."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## LU decomposition - implementation\n",
"\n",
"So we can build an LU code easily from our Gaussian elimination code from last lecture which already works out and performs these tasks. \n",
"\n",
"The final $U$ matrix we need here is as was already constructed through Gaussian elimination, the entries of $L$ we need are simply the ${A_{ik}}/{A_{kk}}$ multipliers we computed as part of the elimination, but threw away previously.\n",
"\n",
"For a given pivot row $k$, for each of these multipliers (for every row below the pivot), as we compute them we know that we are going to transform the augmented matrix in order to achieve a new zero below the diagonal - we can store each multiplier in this position before moving on to the following row, we implicitly know that the diagonals of $L$ will be unity and so don't need to store these (and noting that we don't actually have a space for them anyway!). We then move on to the following pivot row, replacing the zeros in the corresponding column we are zero'ing, but again using the now spare space to store the multipliers.\n",
"\n",
"For example, for the case above \n",
"\n",
"$$ A = \n",
" \\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{black}5} & {\\color{black}14} & {\\color{black}7} & {\\color{black}10}\\\\\n",
" {\\color{black}20} & {\\color{black}77} & {\\color{black}41} & {\\color{black}48}\\\\\n",
" {\\color{black}25} & {\\color{black}91} & {\\color{black}55} & {\\color{black}67}\\\\ \n",
" \\end{bmatrix}\\quad\\rightarrow\\quad\n",
" \\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{blue}{1}} & {\\color{black}{7}} & {\\color{black}{2}} & {\\color{black}{1}}\\\\\n",
" {\\color{blue}{4}} & {\\color{black}{49}} & {\\color{black}{21}} & {\\color{black}{12}}\\\\\n",
" {\\color{blue}{5}} & {\\color{black}{56}} & {\\color{black}{30}} & {\\color{black}{22}}\\\\ \n",
" \\end{bmatrix}\\quad\\rightarrow\\quad\n",
" \\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{blue}1} & {\\color{black}7} & {\\color{black}2} & {\\color{black}1}\\\\\n",
" {\\color{blue}4} & {\\color{blue}{7}} & {\\color{black}{7}} & {\\color{black}{5}}\\\\\n",
" {\\color{blue}5} & {\\color{blue}{8}} & {\\color{black}{14}} & {\\color{black}{14}}\\\\\n",
" \\end{bmatrix}\\quad\\rightarrow\\quad\n",
" \\begin{bmatrix}\n",
" {\\color{black}5} & {\\color{black}7} & {\\color{black}5} & {\\color{black}9}\\\\\n",
" {\\color{blue}1} & {\\color{black}7} & {\\color{black}2} & {\\color{black}1}\\\\\n",
" {\\color{blue}4} & {\\color{blue}7} & {\\color{black}7} & {\\color{black}5}\\\\\n",
" {\\color{blue}5} & {\\color{blue}8} & {\\color{blue}{2}} & {\\color{black}{4}}\\\\\n",
" \\end{bmatrix}\n",
" = [\\color{blue}L\\backslash U]\n",
"$$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6.2: LU decomposition\n",
"\n",
"Starting from your Gaussian elimination code produce a new code to compute the LU decomposition of a matrix. First, store $L$ and $U$ in two different matrices; you could then try a version where you store the entries of both $L$ and $U$ in $A$ as described above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Partial pivoting\n",
"\n",
"At the end of last week we commented that a problem could occur in a numeical implementation of our algorithm if we had a situation where the $A_{kk}$ we divide through by in the Gaussian elimination and/or back substitution algorithms might be zero (or close to zero).\n",
"\n",
"Using Gaussian elimination as an example, let's again consider the algorithm mid-way working on an arbitrary matrix system, i.e. assume that the first $k$ rows have already been transformed into upper-triangular form, while the equations/rows below are not yet in this form:\n",
"\n",
"$$\n",
"\\left[\n",
" \\begin{array}{rrrrrrr|r}\n",
" A_{11} & A_{12} & A_{13} & \\cdots & A_{1k} & \\cdots & A_{1n} & b_1 \\\\\n",
" 0 & A_{22} & A_{23} & \\cdots & A_{2k} & \\cdots & A_{2n} & b_2 \\\\\n",
" 0 & 0 & A_{33} & \\cdots & A_{3k} & \\cdots & A_{3n} & b_3 \\\\\n",
" \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
"\\hdashline \n",
" 0 & 0 & 0 & \\cdots & A_{kk} & \\cdots & A_{kn} & b_k \\\\ \n",
" \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
" 0 & 0 & 0 & \\cdots & A_{nk} & \\cdots & A_{nn} & b_n \\\\\n",
"\\end{array}\n",
"\\right]\n",
"$$\n",
"\n",
"Note we have drawn the horizontal dashed line one row higher, as we are not going to blindly asssume that it is wise to take the current row $k$ as the pivot row, and $A_{kk}$ as the so-called pivot element.\n",
"\n",
"*Partial pivoting* selects the best pivot (row or element) as the one where the $A_{ik}$ (for $i\\ge k$) value is largest (relative to the other values of components in its own row $i$), and then swaps this row with the current $k$ row.\n",
"\n",
"To generalise our codes above we would simply need to search for this row, and perform the row swap operation.\n",
"\n",
"Python's `scipy.linalg` library has its own implementation of the LU decomposition, that uses partial pivoting."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 5. 7. 5. 9.]\n",
" [ 5. 14. 7. 10.]\n",
" [20. 77. 41. 48.]\n",
" [25. 91. 55. 67.]]\n"
]
}
],
"source": [
"A=np.array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91. ,55., 67.]])\n",
"print(A)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[0. 1. 0. 0.]\n",
" [0. 0. 0. 1.]\n",
" [0. 0. 1. 0.]\n",
" [1. 0. 0. 0.]]\n"
]
}
],
"source": [
"P,L,U=sl.lu(A)\n",
"\n",
"# P here is a 'permutation matrix' that performs swaps based upon partial pivoting\n",
"print(P) "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 1. 0. 0. 0. ]\n",
" [ 0.2 1. 0. 0. ]\n",
" [ 0.8 -0.375 1. 0. ]\n",
" [ 0.2 0.375 0.33333333 1. ]]\n"
]
}
],
"source": [
"# the lower-triangular matrix\n",
"print(L) "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 25. 91. 55. 67. ]\n",
" [ 0. -11.2 -6. -4.4 ]\n",
" [ 0. 0. -5.25 -7.25 ]\n",
" [ 0. 0. 0. 0.66666667]]\n"
]
}
],
"source": [
"# the upper-triangular matrix\n",
"print(U) "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 5. 7. 5. 9.]\n",
" [ 5. 14. 7. 10.]\n",
" [20. 77. 41. 48.]\n",
" [25. 91. 55. 67.]]\n"
]
}
],
"source": [
"# double check that P*L*U does indeed equal A\n",
"print(P @ L @ U)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Looking at the form of $P$ above, we can re-order the rows in advance and consider the LU decomposition of the matrix where $P=I$, as below. As we haven't bothered implementing pivoting ourselves, check that your LU implementation recreates the $A$, $L$ and $U$ below."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[25. 91. 55. 67.]\n",
" [ 5. 7. 5. 9.]\n",
" [20. 77. 41. 48.]\n",
" [ 5. 14. 7. 10.]]\n"
]
}
],
"source": [
"A=np.array([[25. ,91. ,55. ,67.],\n",
" [ 5., 7., 5., 9.], \n",
" [20., 77., 41., 48.],\n",
" [ 5., 14., 7., 10.]])\n",
"print(A)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1. 0. 0. 0.]\n",
" [0. 1. 0. 0.]\n",
" [0. 0. 1. 0.]\n",
" [0. 0. 0. 1.]]\n"
]
}
],
"source": [
"P,L,U=sl.lu(A)\n",
"# P now should be the identity as pivoting no longer actually actions any row swaps with this A\n",
"print(P)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 1. 0. 0. 0. ]\n",
" [ 0.2 1. 0. 0. ]\n",
" [ 0.8 -0.375 1. 0. ]\n",
" [ 0.2 0.375 0.33333333 1. ]]\n"
]
}
],
"source": [
"print(L)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 25. 91. 55. 67. ]\n",
" [ 0. -11.2 -6. -4.4 ]\n",
" [ 0. 0. -5.25 -7.25 ]\n",
" [ 0. 0. 0. 0.66666667]]\n"
]
}
],
"source": [
"print(U)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[25. 91. 55. 67.]\n",
" [ 5. 7. 5. 9.]\n",
" [20. 77. 41. 48.]\n",
" [ 5. 14. 7. 10.]]\n"
]
}
],
"source": [
"print(P @ L @ U)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6.3: Partial pivoting\n",
"\n",
"Implement partial pivoting."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": false,
"sideBar": false,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "427.513px"
},
"toc_section_display": false,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}