Objectives
-
Find the \(LU\)-factorization of a matrixDetermine when an \(LU\)-factorization does not exist.Use a \(LU\) factorization to solve a linear system.Use a \(LU\) factorization to find an inverse matrix.
Section 5.2 LU Factorization
Let \(A\) be a square matrix. We seek to write \(A=LU\text{,}\) where \(L\) is a lower triangular matrix and \(U\) is an upper triangular matrix. The reason for doing such a factorization is often to solve \(A\vec{x}=\vec{b}\) for multiple right hand sides.
Recall that an upper triangle matrix is a matrix (not-necessarily square) such that all elements below the diagonal are 0. A lower triangular matrix is a matrix such that all elements above the diagonal are 0. More precisely,
The following example shows that there is a lower-triangular matrix \(L\) and an upper triangular matrix \(U\) whose product is the original matrix \(A\text{.}\) Later we will show where this arises.
Example 5.2.1.
Let
\begin{align*}
L \amp= \begin{bmatrix}
1 \amp 0 \amp 0 \\ 2 \amp 1 \amp 0 \\ -2 \amp -3 \amp 1
\end{bmatrix}
\amp U \amp = \begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
0 \amp 0 \amp 1
\end{bmatrix} \amp
A \amp= \begin{bmatrix}
3 \amp -2 \amp 1 \\
6 \amp 0 \amp 4 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix}
\end{align*}
Show that \(A=LU\)
Solution.
Using standard matrix multiplication:
\begin{align*}
A \amp = \begin{bmatrix}
1(3)+0(0)+0(0) \amp 1(-2) +0(4)+0(0) \amp 1(1)+0(2)+0(1) \\
2(3)+1(0)+0(0) \amp 2(-2)+1(4)+0(0) \amp 2(1)+1(2)+0(1) \\
-2(3)-3(0)+1(0) \amp -2(-2)-3(4)+1(0)\amp-2(1)-3(2)+1(1)
\end{bmatrix} \\
\amp \begin{bmatrix}
3 \amp -2 \amp 1 \\ 6 \amp 0 \amp 4 \\ -6 \amp -8 \amp -7
\end{bmatrix}
\end{align*}
This is nice, but the goal of LU factorization is to take any matrix \(A\) and find matrices \(L\) and \(U\text{.}\) We first, show how to do this for a 3 by 3 matrix and then extend it.
Subsection 5.2.1 Do All Matrices have an \(LU\)-factorization?
We first ask the question whether or not all matrices have a factorization and we will see that not all do not. We can try to factor
\begin{equation*}
\begin{bmatrix} 0 \amp 1 \\ 1 \amp 0 \end{bmatrix}
\end{equation*}
by writing down two matrices as follows.
\begin{equation*}
\begin{bmatrix} 0 \amp 1 \\ 1 \amp 0 \end{bmatrix} =
\begin{bmatrix} a \amp 0 \\ b \amp c \end{bmatrix}
\begin{bmatrix} x \amp y \\ 0 \amp z \end{bmatrix}
\end{equation*}
and then multiplying the two matrices on the right to get
\begin{equation*}
\begin{bmatrix} 0 \amp 1 \\ 1 \amp 0 \end{bmatrix} =
\begin{bmatrix} ax \amp ay \\ bx \amp by + cz \end{bmatrix}
\end{equation*}
Notice that the upper right and lower left corners are \(ay=1\) and \(bx=1\text{.}\) However, the upper left leads to \(ax=0\text{,}\) so either \(a=0\) or \(x=0\text{,}\) which results in a contradiction.
Theorem 5.2.2.
If a matrix \(A\) can be row-reduced without any row swaps, then there exists matrices \(L\) and \(U\) such that \(A=LU\text{.}\)
Subsection 5.2.2 LU fatorization and elementary matrices
We will see that knowing the elementary matrices related to reducing a matrix will lead naturally to its \(LU\) factorization. Consider the following example:
\begin{equation*}
A = \begin{bmatrix}
3 \amp -2 \amp 1 \\
6 \amp 0 \amp 4 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix}
\end{equation*}
We will row-reduce \(A\) using elementary row operations. Since the first step would be \(-2R_1+R_2\text{,}\) then
\begin{equation*}
E_{-2R_1+R_2} A =
\begin{bmatrix}
1 \amp 0 \amp 0 \\
-2 \amp 1 \amp 0 \\
0 \amp 0 \amp 1
\end{bmatrix}
\begin{bmatrix}
3 \amp -2 \amp 1 \\
6 \amp 0 \amp 4 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix}
= \begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix}
\end{equation*}
Next, we elimnate the \(-6\) in the lower left corner with
\begin{equation*}
E_{2R_1+R_3}E_{-2R_1+R_2} A =
\begin{bmatrix}
1 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \\
2 \amp 0 \amp 1
\end{bmatrix}
\begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix}
= \begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
0 \amp -12 \amp -5
\end{bmatrix}
\end{equation*}
And the next step is to perform \(3R_2+R_3\text{,}\) so
\begin{equation*}
E_{3R_2+R_3} E_{2R_1+R_3}E_{-2R_1+R_2} A =
\begin{bmatrix}
1 \amp 0 \amp 0 \\
0 \amp 1 \amp 0 \\
0 \amp 3 \amp 1
\end{bmatrix}
\begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
0 \amp -12 \amp -5
\end{bmatrix}
= \begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
0 \amp 0 \amp -5
\end{bmatrix}
\end{equation*}
and since this is an upper triangular matrix, then note that we have the situation that \(B A = U\) and mutliplying through by \(B^{-1}\) would result in \(A= B^{-1} U\text{.}\) Since \(B\) is the product of lower triangular elementary matrices, then the inverse is the product of these matrices, itself a lower triangular matrix.
More specifically, since
\begin{equation*}
B = E_{3R_2+R_3} E_{2R_1+R_3}E_{-2R_1+R_2}
\end{equation*}
then
\begin{equation*}
B^{-1} = E_{2R_1 + R_2} E_{-2R_1+R_3} E_{-3R_2 + R_3}
\end{equation*}
where the property \((E_1 E_2)^{-1} = E_2^{-1} E_1^{-1} \) from Theorem 2.3.16 is used. This results in
\begin{align*}
B^{-1} \amp = \begin{bmatrix} 1 \amp 0 \amp 0 \\ 2 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1 \end{bmatrix}
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ -2 \amp 0 \amp 1 \end{bmatrix}
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 0 \amp -3 \amp 1 \end{bmatrix} \\
\amp = \begin{bmatrix}
1 \amp 0 \amp 0 \\ 2 \amp 1 \amp 0 \\ -2 \amp -3 \amp 1
\end{bmatrix}
\end{align*}
which is the matrix \(L\text{.}\) Therefore, we can write
\begin{equation*}
\begin{bmatrix}
3 \amp -2 \amp 1 \\
6 \amp 0 \amp 4 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix} =
\begin{bmatrix}
1 \amp 0 \amp 0 \\ 2 \amp 1 \amp 0 \\ -2 \amp -3 \amp 1
\end{bmatrix}
\begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 0 \\
0 \amp 0 \amp -5
\end{bmatrix}
\end{equation*}
Remark 5.2.3. Finding the LU factorization.
Consider a matrix \(A\) which we desire to write as \(A=LU\) for lower triangular matrix \(L\) and upper triangular matrix \(U\text{.}\) If \(A\) can be row reduced to an upper matrix \(U\text{,}\) with the series of elementary matrices \(E_1, E_2, \ldots E_k\) then
\begin{equation*}
L = E_k^{-1}, E_{k-1}^{-1} \cdots E_2^{-1} E_1^{-1}
\end{equation*}
Also the example given above was a square matrix, this does not have to be the case. Consider the next example.
Example 5.2.4.
Find the \(LU\)-decomposition of
\begin{equation*}
A = \begin{bmatrix} 1 \amp 3 \amp 0 \amp -6 \\
2 \amp 4 \amp -1 \amp 0 \\
0 \amp 1 \amp 0 \amp 7
\end{bmatrix}
\end{equation*}
Solution.
We will perform row operations on \(A\) to get to row-reduced form
\begin{align*}
-2R_1+R_2 \to R_2 \amp \qquad
\begin{bmatrix} 1 \amp 3 \amp 0 \amp -6 \\
0 \amp -2 \amp -1 \amp 12 \\
0 \amp 1 \amp 0 \amp 7
\end{bmatrix}\\
-\frac{1}{2} R_2 \to R_2 \amp \qquad
\begin{bmatrix} 1 \amp 3 \amp 0 \amp -6 \\
0 \amp 1 \amp 1/2 \amp -6 \\
0 \amp 1 \amp 0 \amp 7
\end{bmatrix}\\
-R_2 + R_3 \to R_3 \amp \qquad
\begin{bmatrix} 1 \amp 3 \amp 0 \amp -6 \\
0 \amp 1 \amp 1/2 \amp -6 \\
0 \amp 0 \amp -1/2 \amp 13
\end{bmatrix}
\end{align*}
And this shows that the product of the elementary matrices \(E_{-2R_1+R_2}, E_{(-1/2)R_2}, E_{-R_2+R_3}\) would reduce the matrix \(A\) to row-reduced form. Thus the lower triangular matrix is
\begin{align*}
L \amp = E^{-1}_{-R_2+R_3} E^{-1}_{(-1/2)R_2} E^{-1}_{-2R_1+R2} \\
\amp = E_{R_2+R_3} E_{-2 R_2} E_{2R_1+R_2} \\
\amp =
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 0 \amp 1 \amp 1 \end{bmatrix}
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp -2 \amp 0 \\ 0 \amp 0 \amp 1 \end{bmatrix}
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 2 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1 \end{bmatrix}\\
\amp =
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 2 \amp -2 \amp 0 \\ 0 \amp 1 \amp 1 \end{bmatrix}
\end{align*}
Note that since elementary matrices are square, and \(L\) is the product of the elementary matrices that \(L\) is also square. For a general \(m \times n\) matrix, \(L\) is a \(m \times m\) matrix and \(U\) is a \(m \times n\) matrix.
Subsection 5.2.3 Generalization of LU Decomposition
As we saw above in Subsection 5.2.1 not all matrices have an LU decomposition. However, if we generalize this a bit, then we can.
Additionally, an square matrix \(A\) can be written as \(PA=LU\text{,}\) where \(L\) is a lower-diagonal, \(U\) is upper-diagonal and \(P\) is a permutation matrix.
Example 5.2.5.
Find the factorization \(PA=LU\) for
\begin{equation*}
A = \begin{bmatrix} 0 \amp 2 \amp 1 \\ 1 \amp 4 \amp 0 \\ 2 \amp 0 \amp -3 \end{bmatrix}
\end{equation*}
Solution.
If we perform row operations on \(A\) as follows
\begin{align}
R_1 \leftrightarrow R_2 \amp \qquad
\begin{bmatrix} 1 \amp 4 \amp 0 \\ 0 \amp 2 \amp 1 \\ 2 \amp 0 \amp -3 \end{bmatrix} \notag\\
-2 R_1 + R_3 \to R_3 \amp \qquad
\begin{bmatrix} 1 \amp 4 \amp 0 \\ 0 \amp 2 \amp 1 \\ 0 \amp -8 \amp -3 \end{bmatrix} \notag\\
4R_2 + R_3 \to R_3 \amp \qquad
\begin{bmatrix} 1 \amp 4 \amp 0 \\ 0 \amp 2 \amp 1 \\ 0 \amp 0 \amp 1 \end{bmatrix} \tag{5.2.1}
\end{align}
And since the first step was a permutation matrix, then
\begin{equation}
P = E_{1 \leftrightarrow 2}
= \begin{bmatrix} 0 \amp 1 \amp 0 \\ 1 \amp 0 \amp 0 \\ 0 \amp 0 \amp 1\end{bmatrix}\tag{5.2.2}
\end{equation}
\begin{align*}
L \amp = E^{-1}_{4R_2+R_3} E^{-1}_{-2R_1+R_3} = E_{-4R_2+R_3}E_{2R_1+R_3} \\
\amp = \begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 0 \amp -4 \amp 1 \end{bmatrix}
\begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 2 \amp 0 \amp 1 \end{bmatrix}\\
\amp = \begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 2 \amp -4 \amp 1 \end{bmatrix}
\end{align*}
The matrix \(L\) is given above, then \(U\) is the matrix in (5.2.1) and \(P\) is given in (5.2.2). The matrix factorization thus is
\begin{equation*}
PA=LU
\end{equation*}
Subsection 5.2.4 Solving linear systems using LU Decomposition
The main application of LU decomposition is that of solving linear systems. Consider the matrix equation \(A\vec{x} = \vec{b}\) and assume that \(A=LU\) (that is the needed permutation matrix is \(P=I\)).
First, this works, because if \(\vec{y} = U\vec{x}\text{,}\) then
\begin{align*}
L\vec{y} \amp = \vec{b} \\
L U\vec{x} \amp = \vec{b} \\
A \vec{x} \amp = \vec{b}
\end{align*}
Example 5.2.6.
\begin{align*}
A \amp = \begin{bmatrix}
3 \amp -2 \amp 1 \\
6 \amp 0 \amp 4 \\
-6 \amp -8 \amp -7 \\
\end{bmatrix} \amp \vec{b} \amp = \begin{bmatrix}
14 \\22 \\ -9
\end{bmatrix}
\end{align*}
Solution.
\begin{align*}
L \amp = \begin{bmatrix}
1 \amp 0 \amp 0 \\
2 \amp 1 \amp 0 \\
-2 \amp -3 \amp 1
\end{bmatrix}
\amp U \amp = \begin{bmatrix}
3 \amp -2 \amp 1 \\
0 \amp 4 \amp 2 \\
0 \amp 0 \amp 1
\end{bmatrix}
\end{align*}
To use this to solve \(A\vec{x}=\vec{b}\text{,}\) first we solve \(L \vec{y} = \vec{b}\) via forward substitution.
\begin{align*}
y_1 \amp = 14 \\
2y_1 + y_2 \amp = 22 \\
\text{so}\qquad\qquad\amp\\
y_2 \amp = 22-2y_1 = 22-28 =-6 \\
\text{and the last equation is}\qquad\qquad\amp\\
-2 y_1 -3y_2 + y_3 \amp = -9 \\
\text{or}\qquad\qquad\amp\\
y_3 \amp = -9 + 2(14)+3(-6) = 1
\end{align*}
So the solution to \(L\vec{y}=\vec{b}\) is \(\vec{y}=[14,-6,1]^{\intercal}\text{.}\) Next, we solve \(U\vec{x}=\vec{y}\) by back substitution,
\begin{align*}
x_3 \amp = 1 \\
4 x_2 + 2 x_3 \amp = -6 \\
\amp \qquad \qquad \text{or}\\
x_2 \amp = (-6 -2 (1))/4 = -2 \\
\amp \qquad \qquad \text{and the first equation is}\\
3 x_1 -2x_2+x_3 \amp = 14 \\
\amp \qquad \qquad \text{or}\\
3x_1 \amp = 14+2(-2)-x_3 = 9\\
x_1 \amp = 3
\end{align*}
So the solution is
\begin{equation*}
\vec{x} = \begin{bmatrix}
1 \\ -2 \\ 3
\end{bmatrix}
\end{equation*}
You can see from the example above that once the matrix \(A\) factored that it is relatively simple to find the solutions to \(L\vec{y}=\vec{b}\) and \(A\vec{x}=\vec{y}\) and there are relatively few computations to perform.
In fact, this is true in general. In that if one actually finds the factorization of \(A\) and then solves this in the manner show that about 1/2 of the operations are done then solving \(A\vec{x}=\vec{b}\) directly, say by reducing the matrix to reduced row echelon form.
Subsection 5.2.5 Inverting a Matrix
The same idea can be use to find inverse of \(A\) by repeated solving \(LU\vec{x} =
\vec{b}\) by repeating this for \(\vec{b}\) the columns of the identity matrix. The example below shows this without all of the details of finding the factorization:
Example 5.2.7.
Find the inverse of the matrix in Example 2.3.12,
\begin{equation*}
A = \begin{bmatrix}
3 \amp 3 \amp -1 \\ 2 \amp 1 \amp -3 \\ 0 \amp 2 \amp 5
\end{bmatrix}
\end{equation*}
using the LU-decomposition.
Solution.
Following the steps above, the LU decomposition is
\begin{align*}
L \amp = \begin{bmatrix}
1 \amp 0 \amp 0 \\ 2/3 \amp 1 \amp 0 \\
0 \amp -2 \amp 1
\end{bmatrix} \amp U \amp = \begin{bmatrix}
3 \amp 3 \amp -1 \\
0 \amp -1 \amp -7/3 \\
0 \amp 0 \amp 1/3
\end{bmatrix}
\end{align*}
and we will solve \(A\vec{x}=LU\vec{x}=\vec{e}_j\text{,}\) where \(\vec{e}_j\) is the \(j\)th column of the 3 by 3 identity matrix.
\begin{equation*}
\vec{y_1} = \begin{bmatrix}
1 \\ -2/3 \\ -4/3
\end{bmatrix}
\end{equation*}
\begin{equation*}
\vec{x}_1 = \begin{bmatrix}
-11 \\ 10 \\ -4
\end{bmatrix}
\end{equation*}
Repeating this for \(\vec{b}=\vec{e}_2\text{,}\) first solving \(L\vec{y}_2 =
(0,1,0)^{\intercal}\) or
\begin{equation*}
\vec{y}_2 = \begin{bmatrix}
0 \\ 1 \\ 2
\end{bmatrix}
\end{equation*}
\begin{equation*}
\vec{x}_2 = \begin{bmatrix}
17 \\ -15 \\ 6
\end{bmatrix}
\end{equation*}
and lastly, solve \(L\vec{y}_3 = (0,0,1)^{\intercal}\) or
\begin{equation*}
\vec{y}_3 = \begin{bmatrix}
0 \\ 0 \\ 1
\end{bmatrix}
\end{equation*}
and solving \(U\vec{x}_3 = \vec{y}_3\) results in
\begin{equation*}
\vec{x}_3 = \begin{bmatrix}
8 \\ -7 \\ 3
\end{bmatrix}
\end{equation*}
The inverse matrix is then the matrix with \(\vec{x}_j\) as the columns or
\begin{equation*}
A^{-1} = \begin{bmatrix}
-11 \amp 17 \amp 8 \\
10 \amp -15 \amp -7 \\
-4 \amp 6 \amp 3
\end{bmatrix}
\end{equation*}
