{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Matrix factorization"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Let $L$ be a lower-triangular matrix ($n \\times n$) and $U$ be an upper-triangular matrix (also $n \\times n$).\n",
"\n",
"We've discussed how to solve $Ux = y$ with backward substitution:\n",
"\n",
"$$ x_i = \\frac{\\displaystyle y_i - \\sum_{j = i+1}^n u_{ij}x_j}{u_{ii}},\\quad i = n, n-1, \\ldots, 2,1.$$\n",
"\n",
"We also discussed how to solve $Ly = b$ with forward substitution:\n",
"\n",
"$$ y_i = \\frac{\\displaystyle b_i - \\sum_{j=1}^{i-1} l_{ij} y_j }{l_{ii}}, \\quad i = 1,2,\\ldots,n-1,n.$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Using an LU factorization\n",
"\n",
"Assume that we can express $A$ as $LU$ and we want to solve $Ax=b$:\n",
"\n",
"$$A x = LU x = b = L(Ux) = b.$$\n",
"\n",
"1. Solve $Ly = b$ for $y = Ux$ by forward substitution.\n",
"2. Solve $Ux = y$ for $x$ by backward substitution.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Operation count\n",
"\n",
"Recall from Lecture 9 that backward substitution requires $O(n^2)$ operations. The same is true of forward substitution."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"To solve $Ax = b$ when a factorization $A = LU$ is known therefore also requires $O(n^2)$ operations. \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"This should be compared with Gaussian elimination with backward substitution (or worse, matrix inversion), which requires $O(n^3)$ operations."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"But what if the $LU$ factorization is not known in advance? How much effort is required to compute it? The answer is $O(n^3)$ because Gaussian elimination is the means by which we $LU$-factorize, as we will see next. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Thus $LU$ factorization and Gaussian elimination are just two sides of the same coin. "
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"#### The LU factorization algorithm\n",
"\n",
"Assume $A$ is an $n\\times n$ matrix that can be put in upper-triangular form without row swaps (we'll deal with those next lecture).\n",
"\n",
"Consider the $3\\times 3$ case first\n",
"\n",
"$$A = \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix}$$\n",
"\n",
"$$ R_2 - \\frac{a_{21}}{a_{11}}R_1 \\to R_2$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"$$ R_2 - \\frac{a_{21}}{a_{11}}R_1 \\to R_1$$\n",
"\n",
"$$A = \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix}$$\n",
"\n",
"Here $a^{(1)}_{22} = a_{22} - a_{12} \\frac{a_{21}}{a_{11}}$ and $a^{(1)}_{23}= a_{23} - a_{13} \\frac{a_{21}}{a_{11}}$. The question that will lead us to an $LU$ factorization is: What type of matrix does this row operation correspond to?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"$$ \\begin{bmatrix} 1 & 0 & 0 \\\\ -a_{21}/a_{11} & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix}\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Recall that the goal is to express $A = LU$, so we find the inverse of the lower-triangular matrix:\n",
"\n",
"$$\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix}\\begin{bmatrix} 1 & 0 & 0 \\\\ -a_{21}/a_{11} & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix} = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix} $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"So, after a single step, we have found\n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix}\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix}$$\n",
"\n",
"This is closer to a lower/upper factorization."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Next we perform the row operation $R_{3} - \\frac{a_{31}}{a_{11}} R_1 \\to R_3$ which corresponds to\n",
"\n",
"$$\\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ - a_{31}/a_{11} & 0 & 1 \\end{bmatrix}\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} = \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix}$$\n",
"\n",
"where $a^{(1)}_{32} = a_{32} - a_{12} \\frac{a_{31}}{a_{11}}$ and $a^{(1)}_{33} = a_{33} - a_{13} \\frac{a_{31}}{a_{11}}$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"The inverse of the lower-triangular factor can be confirmed\n",
"\n",
"$$\\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ a_{31}/a_{11} & 0 & 1 \\end{bmatrix}\\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ - a_{31}/a_{11} & 0 & 1 \\end{bmatrix} = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0& 0 & 1 \\end{bmatrix}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Thus:\n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ a_{31}/a_{11} & 0 & 1 \\end{bmatrix}\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix}$$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Combine the two row operations to yield:\n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix} \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ a_{31}/a_{11} & 0 & 1 \\end{bmatrix}\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix}$$\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"which can be written as: \n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ a_{31}/a_{11} & 0 & 1 \\end{bmatrix} \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix}$$\n",
"\n",
"This brings us very close to an $LU$ factorization"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"For the last step, we need to perform $R_3 - \\frac{a^{(1)}_{32}}{a^{(1)}_{22}} R_2 \\to R_3$\n",
"\n",
"$$ \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & -a^{(1)}_{32}/a^{(1)}_{22} & 1 \\end{bmatrix} \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix} = \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & 0 & a^{(2)}_{33} \\end{bmatrix}$$\n",
"\n",
"where $a^{(2)}_{33} = a^{(1)}_{33} - a^{(1)}_{23} \\frac{a^{(1)}_{32}}{a^{(1)}_{22}}$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Again, you can check that the inverse of the matrix on the left is simple\n",
"\n",
"$$ \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & a^{(1)}_{32}/a^{(1)}_{22} & 1 \\end{bmatrix} \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & -a^{(1)}_{32}/a^{(1)}_{22} & 1 \\end{bmatrix} = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix} $$\n",
"\n",
"so that we get: \n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix} = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & a^{(1)}_{32}/a^{(1)}_{22} & 1 \\end{bmatrix} \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & 0 & a^{(2)}_{33} \\end{bmatrix}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Inserting this into the existing factorization:\n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ a_{31}/a_{11} & 0 & 1 \\end{bmatrix} \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & a^{(1)}_{32} & a^{(1)}_{33} \\end{bmatrix}$$\n",
"\n",
"yields: \n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ a_{31}/a_{11} & 0 & 1 \\end{bmatrix} \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & a^{(1)}_{32}/a^{(1)}_{22} & 1 \\end{bmatrix} \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & 0 & a^{(2)}_{33} \\end{bmatrix}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"The inner factors simplify and we have an LU factorization\n",
"\n",
"$$\\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ a_{21} & a_{22} & a_{23} \\\\ a_{31} & a_{32} & a_{33} \\end{bmatrix} =\\begin{bmatrix} 1 & 0 & 0 \\\\ a_{21}/a_{11} & 1 & 0 \\\\ a_{31}/a_{11} & a^{(1)}_{32}/a^{(1)}_{22} & 1 \\end{bmatrix} \\begin{bmatrix} a_{11} & a_{12} & a_{13} \\\\ 0 & a^{(1)}_{22} & a^{(1)}_{23} \\\\ 0 & 0 & a^{(2)}_{33} \\end{bmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"The LU algorithm (no row interchanges)\n",
"\n",
"INPUT a n x n matrix A\n",
"\n",
"OUTPUT the LU factorization of A, or a message of failure\n",
"\n",
"STEP 1: Set L to be the n x n identity matrix;\n",
"STEP 2: For j = 1, 2, ...,n do STEPS 3-4\n",
" STEP 3: If A(j,j) = 0\n",
" OUTPUT('LU Factorization failed')\n",
" STOP.\n",
" STEP 4: For i = j+1, j+2, ... , n do STEPS 5-6\n",
" STEP 5: Set L(i,j) = A(i,j)/A(j,j);\n",
" STEP 6: Perform row operation Ri - L(i,j)*Rj --> Ri on A;\n",
"STEP 7: OUTPUT([L,A]);\n",
" "
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 1
}