Skip to main content

Print preview

Printing options

Print headers
Print footers

Section 2.5 \(LU\) Factorization

Handout 2.5 \(LU\) Factorization

 
We have learned several techniques for solving the equation \(A\vec x = \vec b\text{.}\) We could solve \(A\vec x=\vec b\) using \(\vec x=A^{-1}\vec b\text{,}\) but this requires computing the inverse of an \(n\times n\) matrix, which becomes difficult for large \(n\text{.}\) Gaussian elimination is better, but still inefficient for large systems.
More efficient and numerically stable methods rely on matrix factorizations. A matrix factorization (or matrix decomposition) is a factorization of a matrix into a product of matrices. Such factorizations help solve \(A\vec x=\vec b\) and reveal matrix structure. Several matrix factorizations appear throughout this course. In this section, we factor a matrix into lower and upper triangular matrices.
A rectangular matrix \(A\) is upper triangular if \(a_{i,j}=0\) whenever \(i>j\text{.}\)
\begin{equation*} \begin{bmatrix} 1 \amp 5 \amp 0\\ 0 \amp 2 \amp 4 \end{bmatrix}, \quad \begin{bmatrix} 1 \amp 0 \amp 0 \amp 1\\ 0 \amp 2 \amp 1 \amp 0\\ 0 \amp 0 \amp 1 \amp 0\\ 0 \amp 0 \amp 0 \amp 1 \end{bmatrix}, \quad \begin{bmatrix} 2\\0\\0\\0 \end{bmatrix} \end{equation*}
A rectangular matrix \(A\) is lower triangular if \(a_{i,j}=0\) whenever \(i<j\text{.}\)
\begin{equation*} \begin{bmatrix} 1 \amp 0 \amp 0\\ 3 \amp 2 \amp 0 \end{bmatrix}, \quad \begin{bmatrix} 3 \amp 0 \amp 0 \amp 0\\ 1 \amp 1 \amp 0 \amp 0\\ 0 \amp 0 \amp 1 \amp 0\\ 0 \amp 2 \amp 0 \amp 1 \end{bmatrix}, \quad \begin{bmatrix} 1\\2\\1\\2 \end{bmatrix} \end{equation*}
Can you have a matrix that is both upper and lower triangular?
Solution.
There are many matrices which are both, such as
\begin{equation*} \begin{bmatrix} 0 \amp 0 \\ 0 \amp 0 \end{bmatrix} \quad \text{and} \quad \begin{bmatrix} 1 \amp 0 \amp 0 \\ 0 \amp 2 \amp 0 \end{bmatrix} \end{equation*}
 

Example 2.31.

If \(A\) is a \(3 \times 2\) matrix, then its LU factorization has the form:
\begin{equation*} A=LU= \begin{bmatrix} 1 \amp 0 \amp 0\\ \ast \amp 1 \amp 0\\ \ast \amp \ast \amp 1 \end{bmatrix} \begin{bmatrix} \ast \amp \ast\\ 0 \amp \ast\\ 0 \amp 0 \end{bmatrix} \end{equation*}
 
How do we find \(L \) and \(U\text{?}\) Suppose \(A\) can be row reduced to echelon form \(U\) without interchanging rows. Then
\begin{equation*} E_p\cdots E_1A=U, \end{equation*}
where the \(E_j\) are elementary row-operation matrices. These matrices are lower triangular and invertible. For example,
\begin{equation*} \begin{bmatrix} 1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0\\ 2 \amp 0 \amp 1 \end{bmatrix}^{-1} = \begin{bmatrix} 1 \amp 0 \amp 0\\ 0 \amp 1 \amp 0\\ -2 \amp 0 \amp 1 \end{bmatrix}. \end{equation*}
Therefore:
\begin{equation*} A = E_1^{-1}\cdots E_p^{-1}U = LU. \end{equation*}
 
To compute an LU factorization:
  1. Reduce \(A\) to an echelon form \(U\) using only row replacement operations. So no row swapping or scaling.
  2. Place entries in \(L\) so that the same sequence of row operations reduce \(L\) to \(I\text{.}\)
Note, in MATH 1554, the only allowed row replacement operation can be written in the form \(R_i - kR_j \rightarrow R_i\text{.}\) The number \(k\) is what should go into the \(i^\text{th}\) row and \(j^\text{th}\) column of \(L\text{,}\) so \(\ell_{i,j} = k\text{.}\)
Caution 1: Note that \(\ell_{i,j}\) is the opposite of the scalar used in the row replacement. This is because the formula \(R_i - kR_j \rightarrow R_i\) has a negative sign built in.
Caution 2: Note that you are only guaranteed to find the correct \(\ell_{i,j}\) if you are working from left to right and from top to bottom. So, start with the leftmost pivot column and use its pivot row to elimnate all nonzero entries below it. Then move to the next pivot column (from the left) and repeat.
 

Example 2.32.

Compute the LU factorization of:
\begin{equation*} A= \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5\\ -16 \amp 12 \amp 2 \amp -17\\ 8 \amp -6 \amp -12 \amp 22 \end{bmatrix}. \end{equation*}
Solution.
First, let’s determine the sizes of \(L\) and \(U\text{.}\) Since \(U\) is an echelon form of \(A\text{,}\) we know that \(U\) has the same dimension as \(A\text{.}\) Thus, \(U\) is a \(3 \times 4\) matrix. Since \(L\) is a square matrix that multiplies \(U\) on the left, \(L\) must be a \(3 \times 3\) matrix.
Now, we determine \(U\) by row reducing \(A\) to echelon form only using row replacements of the form \(R_i - kR_j \rightarrow R_i\text{.}\) Each time we make one of these row replacements, we are finding an entry of \(L\text{,}\) specifically we are finding that \(\ell_{i,j} = k\text{.}\)
I tend to write down the form of \(L\) with blanks in the lower portion and fill them in after each row replacement. So, I start the problem with the following written towards the bottom of the page:
\begin{equation*} L = \begin{bmatrix} 1 \amp 0 \amp 0 \\ \amp 1 \amp 0 \\ \amp \amp 1 \end{bmatrix} \end{equation*}
 
Note the following row operations:
  • \(R_2 + 4R_1 \rightarrow R_2\text{:}\)
    \begin{equation*} A = \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ -16 \amp 12 \amp 2 \amp -17 \\ 8 \amp -6 \amp -12 \amp 22 \end{bmatrix} \sim \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ 0 \amp 0 \amp -2 \amp 3 \\ 8 \amp -6 \amp -12 \amp 22 \end{bmatrix} \end{equation*}
    Note that the row replacement could be written as \(R_2 - (-4)R_1 \rightarrow R_2\text{,}\) so we have that \(\ell_{2,1} = -4\text{.}\)
     
  • \(R_3 - 2R_1 \rightarrow R_3\text{:}\)
    \begin{equation*} \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ 0 \amp 0 \amp -2 \amp 3 \\ 8 \amp -6 \amp -12 \amp 22 \end{bmatrix}\sim \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ 0 \amp 0 \amp -2 \amp 3 \\ 0 \amp 0 \amp -10 \amp 12 \end{bmatrix} \end{equation*}
    Since we used the row operation \(R_3 - 2R_1 \rightarrow R_3\text{,}\) we have that \(\ell_{3,1} = 2\text{.}\)
     
  • \(R_3 - 5R_2 \rightarrow R_3\text{:}\)
    \begin{equation*} \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ 0 \amp 0 \amp -2 \amp 3 \\ 0 \amp 0 \amp -10 \amp 12 \end{bmatrix} \sim \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ 0 \amp 0 \amp -2 \amp 3 \\ 0 \amp 0 \amp 0 \amp -3\end{bmatrix} \end{equation*}
    Since we used the row operation \(R_3 - 5R_2 \rightarrow R_3\text{,}\) we have that \(\ell_{3,2} = 5\text{.}\) Since the matrix is in echelon form, we have found \(U\text{.}\)
 
Therefore, we have that
\begin{equation*} L = \begin{bmatrix} 1 \amp 0 \amp 0 \\ -4 \amp 1 \amp 0 \\ 2 \amp 5 \amp 1 \end{bmatrix} \quad \text{and} \quad U = \begin{bmatrix} 4 \amp -3 \amp -1 \amp 5 \\ 0 \amp 0 \amp -2 \amp 3 \\ 0 \amp 0 \amp 0 \amp -3\end{bmatrix}\text{.} \end{equation*}
We can check by verifying that \(A = LU\text{.}\)
 

Example 2.33.

Compute the LU factorization of:
\begin{equation*} A= \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 6 \amp 0 \amp 10 \amp 12\\ 3 \amp 3 \amp 8 \amp 10\\ 6 \amp 2 \amp 12 \amp 15\\ 3 \amp 3 \amp 8 \amp 9 \end{bmatrix}. \end{equation*}
Solution.
Note \(U\) will be a \(5 \times 4\) matrix and \(L\) will be a \(5 \times 5\) matrix. I start by writing the following near the bottom of the page:
\begin{equation*} L = \begin{bmatrix} 1 \amp 0 \amp 0 \amp 0 \amp 0 \\ \amp 1 \amp 0 \amp 0 \amp 0 \\ \amp \amp 1 \amp 0 \amp 0 \\ \amp \amp \amp 1 \amp 0 \\ \amp \amp \amp \amp 1 \end{bmatrix} \end{equation*}
 
Note the following row operations:
\begin{equation*} A = \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 6 \amp 0 \amp 10 \amp 12\\ 3 \amp 3 \amp 8 \amp 10\\ 6 \amp 2 \amp 12 \amp 15\\ 3 \amp 3 \amp 8 \amp 9 \end{bmatrix} \overset{\begin{array}{rr} R_2 - 2R_1 \rightarrow R_2 \\ R_3 - R_1 \rightarrow R_3 \end{array}}{\underset{\begin{array}{rr} R_4 - 2R_1 \rightarrow R_4 \\ R_5 - R_1 \rightarrow R_5 \end{array}}{\Huge \sim}} \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 0 \amp 2 \amp 2 \amp 2\\ 0 \amp 4 \amp 4 \amp 5\\ 0 \amp 4 \amp 4 \amp 5\\ 0 \amp 4 \amp 4 \amp 4 \end{bmatrix} \end{equation*}
Therefore, we have that \(\ell_{2,1} = 2, \ell_{3,1} = 1, \ell_{4,1} = 2,\) and \(\ell_{5,1} = 1\text{.}\)
 
Note that if we were only interested getting to echelon form, then we could do \(R_4 - R_3 \rightarrow R_4\text{;}\) however, since we are interested in \(LU\)-factorization we must pay attention to Caution 2. Therefore, we must use Row 2 next since it is the next pivot column. Note the following:
\begin{equation*} \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 0 \amp 2 \amp 2 \amp 2\\ 0 \amp 4 \amp 4 \amp 5\\ 0 \amp 4 \amp 4 \amp 5\\ 0 \amp 4 \amp 4 \amp 4 \end{bmatrix} \overset{\begin{array}{r} R_3 - 2R_2 \rightarrow R_3 \end{array}}{\underset{\begin{array}{r} R_4 - 2R_2 \rightarrow R_4 \\ R_5 - 2R_2 \rightarrow R_5 \end{array}}{\Huge \sim}} \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 0 \amp 2 \amp 2 \amp 2\\ 0 \amp 0 \amp 0 \amp 1\\ 0 \amp 0 \amp 0 \amp 1\\ 0 \amp 0 \amp 0 \amp 0 \end{bmatrix} \end{equation*}
Thus, we have that \(\ell_{3,2} = 2, \ell_{4,2} = 2,\) and \(\ell_{5,2} = 2\text{.}\)
 
We are still not in echelon form, so we must do one more row operation.
\begin{equation*} \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 0 \amp 2 \amp 2 \amp 2\\ 0 \amp 0 \amp 0 \amp 1\\ 0 \amp 0 \amp 0 \amp 1\\ 0 \amp 0 \amp 0 \amp 0 \end{bmatrix} \overset{R_4 - R_3 \rightarrow R_4}{\Huge \sim} \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 0 \amp 2 \amp 2 \amp 2\\ 0 \amp 0 \amp 0 \amp 1\\ 0 \amp 0 \amp 0 \amp 0\\ 0 \amp 0 \amp 0 \amp 0 \end{bmatrix} \end{equation*}
Therefore, \(\ell_{4,3} = 1\text{.}\) What about \(\ell_{5,3}\text{?}\) Well, we didn’t need to do anything to Row 5, so can we put anything in that position? NO. Notice, that if you put a non-zero number in that position, say \(\ell_{5,3} = 1\text{,}\) then we are saying that we need to do the row replacement \(R_5 - R_3 \rightarrow R_5\) which would mess up our row of zeros at the bottom of the matrix, thus undoing our work in getting to echelon form. Therefore, we need to make NO change to Row 5, which is doing the row replacement \(R_5 - 0R_3 \rightarrow R_5\text{,}\) so \(\ell_{5,3}\) must be \(0\text{.}\)
 
Since the final matrix is in echelon form, we have found \(U\text{.}\) However, there is still one more entry of \(L\) to determine, \(\ell_{5,4}\text{.}\) Do we have to put a zero there? NO. Note that whatever number is in that position, say \(\ell_{5,4} = 5\text{,}\) would result in the row replacement \(R_5 - 5R_4 \rightarrow R_5\) which actually doesn’t change anything since Row 4 is a row consisting entirely of \(0\)’s. Therefore, \(\ell_{5,4} = t\text{,}\) where \(t \in \mathbb{R}\) can be anything you wish. Thus, we have that
\begin{equation*} L = \begin{bmatrix} 1 \amp 0 \amp 0 \amp 0 \amp 0 \\ 2 \amp 1 \amp 0 \amp 0 \amp 0 \\ 1 \amp 2 \amp 1 \amp 0 \amp 0 \\ 2 \amp 2 \amp 1 \amp 1 \amp 0 \\ 1 \amp 2 \amp 0 \amp t \amp 1 \end{bmatrix} \quad \text{and} \quad U = \begin{bmatrix} 3 \amp -1 \amp 4 \amp 5\\ 0 \amp 2 \amp 2 \amp 2\\ 0 \amp 0 \amp 0 \amp 1\\ 0 \amp 0 \amp 0 \amp 0\\ 0 \amp 0 \amp 0 \amp 0 \end{bmatrix} \end{equation*}
One consequence of this is that the \(LU\) factorization of a matrix might not be unique. It is true that \(U\) is unique; however, there could be many possible \(L\)’s that work.
 
How do we use \(LU\) factorization to help us solve \(A\vec x=\vec b\text{?}\)
If \(A = LU\text{,}\) then the equation \(A \vec x = \vec b \) can be written \(LU \vec x = \vec b \text{,}\) or equivalently, \(L (U\vec x) = \vec b \text{.}\) Setting \(\vec y = U\vec x \text{,}\) we do the following:
  1. Forward substitution to solve for \(\vec y\) in \(L\vec y=\vec b\text{.}\)
  2. Backward substitution to solve for \(\vec x\) in \(U\vec x=\vec y\text{.}\)
Since \(L\) is an invertible matrix, the system \(L\vec y = \vec b\) will always be consistent and have a unique solution. The system \(U\vec x = \vec y\) may be inconsistent. Since \(U\) is in echelon form, we can quickly determine whether \(U\vec x = \vec y\) is inconsistent by checking if there is a nonzero entry in \(\vec y\) that corresponds to a row of zeros in \(U\text{.}\)

Example 2.34.

Solve the linear system with:
\begin{equation*} A=LU= \begin{bmatrix} 1 \amp 0 \amp 0 \amp 0\\ 1 \amp 1 \amp 0 \amp 0\\ 0 \amp 2 \amp 1 \amp 0\\ 0 \amp 0 \amp 1 \amp 1 \end{bmatrix} \begin{bmatrix} 1 \amp 0 \amp 0\\ 0 \amp 2 \amp 1\\ 0 \amp 0 \amp 2\\ 0 \amp 0 \amp 0 \end{bmatrix}, \qquad \vec b= \begin{bmatrix} 2\\5\\8\\2 \end{bmatrix}. \end{equation*}
Solution.
We start with determining a vector \(\vec y\) so that \(L\vec y = \vec b\text{.}\) Note that
\begin{equation*} \begin{bmatrix} L \amp \rvert \amp \vec b \end{bmatrix} = \begin{bmatrix} 1 \amp 0 \amp 0 \amp 0 \amp \rvert \amp 2 \\ 1 \amp 1 \amp 0 \amp 0 \amp \rvert \amp 5 \\ 0 \amp 2 \amp 1 \amp 0 \amp \rvert \amp 8 \\ 0 \amp 0 \amp 1 \amp 1 \amp \rvert \amp 2 \end{bmatrix} \sim \cdots \sim \begin{bmatrix} 1 \amp 0 \amp 0 \amp 0 \amp \rvert \amp 2 \\ 0 \amp 1 \amp 0 \amp 0 \amp \rvert \amp 3 \\ 0 \amp 0 \amp 1 \amp 0 \amp \rvert \amp 2 \\ 0 \amp 0 \amp 0 \amp 1 \amp \rvert \amp 0 \end{bmatrix} \end{equation*}
Therefore, we have that \(\vec y = \begin{bmatrix} 2 \\ 3 \\ 2 \\ 0 \end{bmatrix}\text{.}\)
 
Next, we determine \(\vec x\) so that \(U\vec x = \vec y\text{.}\) Note that
\begin{equation*} \begin{bmatrix} U \amp \rvert \amp \vec y \end{bmatrix} = \begin{bmatrix} 1 \amp 0 \amp 0 \amp \rvert \amp 2 \\ 0 \amp 2 \amp 1 \amp \rvert \amp 3 \\ 0 \amp 0 \amp 2 \amp \rvert \amp 2 \\ 0 \amp 0 \amp 0 \amp \rvert \amp 0 \end{bmatrix} \sim \cdots \sim \begin{bmatrix} 1 \amp 0 \amp 0 \amp \rvert \amp 2 \\ 0 \amp 1 \amp 0 \amp \rvert \amp 1\\ 0 \amp 0 \amp 1 \amp \rvert \amp 1 \\ 0 \amp 0 \amp 0 \amp \rvert \amp 0 \end{bmatrix} \end{equation*}
Therefore, we have that \(\vec x = \begin{bmatrix} 2 \\ 1 \\ 1 \end{bmatrix}\) is the solution to the system \(A \vec x = \vec b\text{.}\)
PrevTopNext