\(LU\) factorization

Consider the equation:
\[\begin{bmatrix}1&0&0\\3&1&0\\-2&2&1\end{bmatrix}\begin{bmatrix}2&-1&2\\0&4&-1\\0&0&3\end{bmatrix}\vec {x}=\begin{bmatrix}1\\8\\5\end{bmatrix}\]
This equation is unusual in that it involves two matrices on the left-hand side. If we multiply the matrices together, we get
\[\begin{bmatrix}2&-1&2\\6&1&5\\-4&10&-3\end{bmatrix}\vec {x}=\begin{bmatrix}1\\8\\5\end{bmatrix}\]
Gaussian elimination yields:
\[\left [\begin{array}{ccc|c} 2&-1&2&1\\6&1&5&8\\-4&10&-3&5 \end{array}\right ] \quad \rightsquigarrow \quad \left [\begin{array}{ccc|c} 1&0&0&2\\0&1&0&1\\0&0&1&-1 \end{array}\right ] \]
We conclude that
\[\vec {x}=\begin{bmatrix}2\\1\\-1\end{bmatrix}.\]

Now that we know the answer, we will turn to the original question and consider the advantages of the format of the original problem. Observe that the two matrices have a distinct pattern. The matrix on the left has zeros above the main diagonal, while the matrix on the right has zeros below the main diagonal. Matrices that follow such a pattern are called lower triangular and upper triangular matrices, respectively.

Let \(\vec {y}=\begin{bmatrix}2&-1&2\\0&4&-1\\0&0&3\end{bmatrix}\vec {x}\). Then we can write the original equation as

\[\begin{bmatrix}1&0&0\\3&1&0\\-2&2&1\end{bmatrix}\left (\begin{bmatrix}2&-1&2\\0&4&-1\\0&0&3\end{bmatrix}\vec {x}\right )=\begin{bmatrix}1&0&0\\3&1&0\\-2&2&1\end{bmatrix}\vec {y}=\begin{bmatrix}1\\8\\5\end{bmatrix}\]
The equation
\[\begin{bmatrix}1&0&0\\3&1&0\\-2&2&1\end{bmatrix}\vec {y}=\begin{bmatrix}1&0&0\\3&1&0\\-2&2&1\end{bmatrix}\begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix}=\begin{bmatrix}1\\8\\5\end{bmatrix}\]
corresponds to the system
\[\begin{matrix} y_1&&&=&1\\ 3y_1 & +y_2&&= &8 \\ -2y_1&+2y_2&+y_3&=&5 \end{matrix}\]

This system can be easily solved by substitution, giving us

\[y_1=1\]
\[3(1)+y_2=8 \rightsquigarrow y_2=5\]
\[-2(1)+2(5)+y_3=5 \rightsquigarrow y_3=-3\]
So,
\[\vec {y}=\begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix}=\begin{bmatrix}1\\5\\-3\end{bmatrix}\]

Recall that we let \(\vec {y}=\begin{bmatrix}2&-1&2\\0&4&-1\\0&0&3\end{bmatrix}\vec {x}\). The equation

\[\begin{bmatrix}2&-1&2\\0&4&-1\\0&0&3\end{bmatrix}\vec {x}=\begin{bmatrix}1\\5\\-3\end{bmatrix}\]
corresponds to the system
\[\begin{matrix} 2x_1& -&x_2&+&2x_3&=&1\\ & &4x_2&-&x_3&= &5 \\ & &&&3x_3&=&-3 \end{matrix}\]
which can be easily solved by back-substitution:
\[x_3=\answer {-1}\]
\[x_2=\answer {1}\]
\[x_1=\answer {2}\]

Thus, \(\vec {x}=\begin{bmatrix}2\\1\\-1\end{bmatrix}\). Observe that this is the same answer that we obtained in the beginning of the problem using Gaussian elimination.

Given a matrix equation such as

\[\begin{bmatrix}2&-1&2\\6&1&5\\-4&10&-3\end{bmatrix}\vec {x}=\begin{bmatrix}1\\8\\5\end{bmatrix}\]

of Exploration init:LUfactorization, often in practice we express the matrix on the left as a product of an upper triangular matrix \(U\) and a lower triangular matrix \(L\) and use substitution to solve the equation instead of using Gaussian elimination to find the solution. In Exploration init:LUfactorization,

\[\begin{bmatrix}2&-1&2\\6&1&5\\-4&10&-3\end{bmatrix}=LU=\begin{bmatrix}1&0&0\\3&1&0\\-2&2&1\end{bmatrix}\begin{bmatrix}2&-1&2\\0&4&-1\\0&0&3\end{bmatrix}\]

The process of taking a matrix \(A\) and expressing it as a product \(A=LU\) of a lower triangular matrix \(L\) and an upper triangular matrix \(U\) is called \(LU\) factorization. We also adopt the (common) convention that \(L\) is a matrix with a diagonal consisting entirely of \(1\)’s, which is often called a lower unit triangular matrix.

Not every matrix has an \(LU\) factorization, but we will see that we can correct for that.

Finding an \(LU\) factorization by Inspection

One way to find the \(LU\) factorizationis to simply look for it directly. We need

\[ \begin{bmatrix} 1 & 2 & 0 & 2 \\ 1 & 3 & 2 & 1 \\ 2 & 3 & 4 & 0 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ x & 1 & 0 \\ y & z & 1 \end{bmatrix} \begin{bmatrix} a & d & h & j \\ 0 & b & e & i \\ 0 & 0 & c & f \end{bmatrix} \]

By multiplying the latter matrices we get

\[ \begin{bmatrix} a & d & h & j \\ xa & xd+b & xh+e & xj+i \\ ya & yd+zb & yh+ze+c & yj+iz+f \end{bmatrix} \]

and from this, it is a simple exercise to determine each of the unknowns. For example, from the first column, we get \(a=1\) and then \(x=1,y=2.\)

See if you can continue to determine the entire \(LU\) factorization of \(A\).

\[ \begin{bmatrix} 1 & 2 & 0 & 2 \\ 1 & 3 & 2 & 1 \\ 2 & 3 & 4 & 0 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 2 & \answer {-1} & 1 \end{bmatrix} \begin{bmatrix} 1 & \answer {2} & \answer {0} & \answer {2} \\ 0 & \answer {1} & \answer {2} & \answer {-1} \\ 0 & 0 & \answer {6} & \answer {-5} \end{bmatrix} \]

Finding an \(LU\) factorization by the Multiplier Method

In the following example we demonstrate a method for computing an \(LU\) factorization known as the multiplier method.

Notice that when we perform the multiplier method, we are making repeated use of a certain type of elementary row operation, namely, we are adding a scalar multiple of one row to a row below it. The reason this helps to create an \(LU\) factorization depends upon the fact that the elementary matrices corresponding to such operations are lower triangular. To understand how this works, we begin with a lemma.

Consider the usual setup for finding the inverse \(\begin{bmatrix} L &|& I \end{bmatrix}\). Each row operation used on \(L\) to transform this matrix to reduced row echelon form changes only the entries in \(I\) below the main diagonal. Whether we have the special case of \(L\) given in 1 where the nonzero nondiagonal entries are in the left column, or if the single nonzero column is in another position, it is clear that multiplication by \(-1\) as described in the lemma gives us \(L^{-1}\).

For a simple illustration of the lemma, observe:

\begin{equation*} \left [\begin{array}{ccc|ccc} 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & a & 1 & 0 & 0 & 1 \end{array}\right ] \rightsquigarrow \left [\begin{array}{ccc|ccc} 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & -a & 1 \end{array}\right ] \end{equation*}

Now let \(A\) be an \(m\times n\) matrix, say

\begin{equation*} A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix}, \end{equation*}

and assume \(A\) can be row reduced to an upper triangular form using only row operation 3. Thus, in particular, \(a_{11}\neq 0\). Multiply on the left by

\begin{equation*} E_{1}= \begin{bmatrix} 1 & 0 & \cdots & 0 \\ - \frac {a_{21}}{a_{11}} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ -\frac {a_{m1}}{a_{11}} & 0 & \cdots & 1 \end{bmatrix} \end{equation*}

This is the product of elementary matrices which make modifications in the first column only. It is equivalent to taking \(-a_{21}/a_{11}\) times the first row and adding to the second. Then taking \(-a_{31}/a_{11}\) times the first row and adding to the third and so forth. The quotients in the first column of the above matrix are the multipliers. Thus the result is of the form

\begin{equation*} E_{1}A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n}^{\prime } \\ 0 & a_{22}^{\prime } & \cdots & a_{2n}^{\prime } \\ \vdots & \vdots & & \vdots \\ 0 & a_{m2}^{\prime } & \cdots & a_{mn}^{\prime } \end{bmatrix} \end{equation*}

By assumption, \(a_{22}^{\prime }\neq 0\) and so it is possible to use this entry to zero out all the entries below it in the matrix on the right by multiplication by a matrix of the form

\begin{equation*} E_{2}= \begin{bmatrix} 1 & \vec {0} \\ \vec {0} & E \end{bmatrix} \end{equation*}

where \(E\) is an \(( m-1)\times ( m-1)\) matrix of the form

\begin{equation*} E=\begin{bmatrix} 1 & 0 & \cdots & 0 \\ -\frac {a_{32}^{\prime }}{a_{22}^{\prime }} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ -\frac {a_{m2}^{\prime }}{a_{22}^{\prime }} & 0 & \cdots & 1 \end{bmatrix} \end{equation*}

Again, the entries in the first column below the 1 are the multipliers. Continuing this way, zeroing out the entries below the diagonal entries, finally leads us to

\begin{equation*} E_{m-1}E_{m-2}\cdots E_{1}A=U \end{equation*}

where \(U\) is upper triangular. Each \(E_{j}\) has all ones down the main diagonal and is lower triangular. Now we can multiply both sides by the inverses of the \(E_{j}\) in the reverse order. This yields

\begin{equation*} A=E_{1}^{-1}E_{2}^{-1}\cdots E_{m-1}^{-1}U \end{equation*}

By Lemma 5, this implies that the product of those \(E_{j}^{-1}\) is a lower triangular matrix having all ones down the main diagonal.

The above discussion and lemma explain how the multiplier method works. The expressions

\begin{equation*} -\frac {a_{21}}{a_{11}},-\frac {a_{31}}{a_{11}},\cdots , -\frac {a_{m1}}{a_{11}}, \end{equation*}

which we obtained in building \(E_{1}\), and which we denote respectively by \(M_{21},\cdots ,M_{m1}\) to save notation, are the multipliers. Therefore, according to the lemma, to find \(E_{1}^{-1}\) we simply write

\begin{equation*} \begin{bmatrix} 1 & 0 & \cdots & 0 \\ -M_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ -M_{m1} & 0 & \cdots & 1 \end{bmatrix} \end{equation*}

Similar considerations apply to the other \(E_{j}^{-1}.\) Thus \(L\) is a product of the form

\begin{equation*} \begin{bmatrix} 1 & 0 & \cdots & 0 \\ -M_{21} & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ -M_{m1} & 0 & \cdots & 1 \end{bmatrix} \cdots \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & 0 & \ddots & \vdots \\ 0 & \cdots & -M_{m(m-1)} & 1 \end{bmatrix}, \end{equation*}

where each factor has at most one nonzero column, the position of which moves from left to right as we scan the above product of matrices from left to right. It follows from what we know about the effect of multiplying on the left by an elementary matrix that the above product is of the form

\begin{equation*} \begin{bmatrix} 1 & 0 & \cdots & 0 & 0 \\ -M_{21} & 1 & \cdots & 0 & 0 \\ \vdots & -M_{32} & \ddots & \vdots & \vdots \\ -M_{(M-1)1} & \vdots & \cdots & 1 & 0 \\ -M_{M1} & -M_{M2} & \cdots & -M_{MM-1} & 1 \end{bmatrix} \end{equation*}

To sum up the procedure, beginning at the left column and moving toward the right, you simply insert, into the corresponding position in the identity matrix, \(-1\) times the multiplier which was used to zero out an entry in that position below the main diagonal in \(A,\) while retaining the main diagonal which consists entirely of ones. This will give us \(L\) as we create \(U\) using row operation 3.

We now return to Example 1, to understand why we could not use the multiplier method to find an \(LU\) factorization of the coefficient matrix. Suppose we had written

\[\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}\begin{bmatrix}0&12&-18\\1&2&-3\\1&-2&1\end{bmatrix}.\]

Our first step of Gaussian elimination would require a row swap to get a nonzero entry into the top left corner (we swapped rows 1 and 2 in Example 1). Unfortunately, the elementary matrix that accomplishes this, \(P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}\) is NOT lower triangular. So we cannot apply Lemma 5 to generate a lower triangular \(L\).

If the elementary matrices used to reduce our matrix to row-echelon form are all lower triangular, then we can find an \(LU\) factorization. But what about the general case?

For a proof of the above theorem, the reader is referred to [Nicholson].

The matrix \(P\) in the above theorem is called a permutation matrix. These matrices have other important applications, as we will see later.

Practice Problems

Use the given \(LU\) factorization of the coefficient matrix to solve the system of equations.
\[\begin{array}{ccccccc} x & +&2y&+&3z&= &5 \\ 2x&+&3y&+&z&=&6\\ 3x& +&5y&+&4z&=&11 \end{array}\]
Observe that an \(LU\) factorization of the coefficient matrix is
\[ \begin{bmatrix} 1 & 2 & 3 \\ 2 & 3 & 1 \\ 3 & 5 & 4\end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \end{bmatrix} \begin{bmatrix} 1 & 2 & 3 \\ 0 & -1 & -5 \\ 0 & 0 & 0 \end{bmatrix} \]

Problems -

Find the \(LU\) factorization of the coefficient matrix and use it to solve the system of equations.

\[ \begin{array}{ccccc} x& +&2y&=&5\\ 2x & +&3y&= &6 \end{array} \]
\[L=\begin{bmatrix}\answer {1}&0\\\answer {2}&\answer {1}\end{bmatrix},\quad U=\begin{bmatrix}\answer {1}&\answer {2}\\0&\answer {-1}\end{bmatrix}\]
\[x=\answer {-3}, y=\answer {4}\]
\[\begin{array}{ccccccc} x & +&2y&+&z&= &1 \\ & &y&+&3z&=&2\\ 2x& +&3y&&&=&6 \end{array}\]
\[L=\begin{bmatrix}\answer {1}&0&0\\\answer {0}&\answer {1}&0\\\answer {2}&\answer {-1}&\answer {1}\end{bmatrix},\quad U=\begin{bmatrix} \answer {1}&\answer {2}&\answer {1}\\0&\answer {1}&\answer {3}\\0&0&\answer {1}\end{bmatrix}\]
\[x=\answer {27},y=\answer {-16}, z=\answer {6}\]
Is there only one \(LU\) factorization for a given matrix?
Consider the equation
\[ \begin{bmatrix}0 & 1 \\0 & 1\end{bmatrix}=\begin{bmatrix}1& 0 \\1 & 1\end{bmatrix} \begin{bmatrix}0 & 1 \\0 & 0\end{bmatrix} \]

Look for all possible \(LU\) factorizations.

Can you show that every permutation matrix is invertible? If so, What does the inverse of a permutation matrix look like? (Recall that a permutation matrix is matrix \(P\) of Theorem 7.)

Text Source

The text in this section is an adaptation of Section 2.2 of Ken Kuttler’s A First Course in Linear Algebra. (CC-BY)

Ken Kuttler, A First Course in Linear Algebra, Lyryx 2017, Open Edition, p. 99-106.

Bibliography

[Nicholson] W. Keith Nicholson, Linear Algebra with Applications, Lyryx 2018, Open Edition, pp. 123–127.