Theory

Introduction

An LU factorization of a matrix $A$ is simply a lower triangular matrix $L$ and an upper triangular matrix $U$ such that $A=LU$.

For example

1. $$A=\left(\begin{array}{rrr}1& 2& 3\\ 2& 5& 10\\ 3& 2& -6\end{array}\right)=\left(\begin{array}{rrr}1& 0& 0\\ 2& 1& 0\\ 3& -4& 1\end{array}\right)\left(\begin{array}{rrr}1& 2& 3\\ 0& 1& 4\\ 0& 0& 1\end{array}\right)=LU$$ or
2. $$A=\left(\begin{array}{rrrrr}2& 4& -1& 5& -2\\ -4& -5& 3& -8& 1\\ 2& -5& -4& 1& 8\\ -6& 0& 7& -3& 1\end{array}\right)=\left(\begin{array}{rrrr}1& 0& 0& 0\\ -2& 1& 0& 0\\ 1& -3& 1& 0\\ -3& 4& 2& 1\end{array}\right)\left(\begin{array}{rrrrr}2& 4& -1& 5& -2\\ 0& 3& 1& 2& -3\\ 0& 0& 0& 2& 1\\ 0& 0& 0& 0& 5\end{array}\right)=LU.$$ Observe:

• In 1., both $L$ and $U$ have ones on the diagonal, and all matrices are square.
• In 2., $A$ and $U$ are not square, $L$ has ones on the diagonal, while $U$ has other values: $2,3,0,0$.

We require the $L$ matrix to have ones on the diagonal, while the $U$ matrix may have other diagonal values.

A more general rule for $U$ called “row-echelon form” will be explained below; it is applicable when zeros on the diagonal crop up in the process; in essence it means the entries should be mashed “high and to the right.”

Keys to LU factorization

The principle of elimination is to use row-add matrices acting on the left in order to eliminate nonzero entries of $A$ that are below the main diagonal entries, clearing a column. These row-add matrices are then collected into a single lower-triangular matrix that combines the actions of all of them.

Creating columns:
Observe first that we can combine row-add matrices which add from the same row into a single matrix having a single column of nonzero subdiagonal entries, for example:

$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 2& 1& 0& 0\\ 3& 0& 1& 0\\ 4& 0& 0& 1\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 2& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 3& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 4& 0& 0& 1\end{array}\right).$$

The order in the product on the RHS does not matter!

Keys to LU factorization (part 2)

Combining columns:
Furthermore, we can break up any lower-triangular matrix $L$ into a product of single-column matrices like the one above:

$$L=\left(\begin{array}{cccc}1& 0& 0& 0\\ 2& 1& 0& 0\\ 3& 5& 1& 0\\ 4& 6& 7& 1\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 2& 1& 0& 0\\ 3& 0& 1& 0\\ 4& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 5& 1& 0\\ 0& 6& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 7& 1\end{array}\right)=L_1\cdot L_2\cdot L_3.$$

Combining columns – order matters

The three factors in the RHS product cannot be reordered without changing the result!

By combining the above two factorizations, we can write any lower-triangular matrix $L$ (with ones on the diagonal) as a product of row-add matrices, and these row-add matrices can be grouped into single-column row-add combo matrices. If $L$ is lower-triangular with $n$ rows, and we write $L_i$ for the single-column combo with nonzeros in the $i^{\text{th}}$ column, then we have $L=L_1\cdot L_2\cdot \cdots \cdot L_{n-1}$. When you see $L_i$ you should think “single-column, row-add combo” and remember that the column of nonzeros lies below the main diagonal of ones.

Inverse of $L_i$:
The last critical fact is that the inverse of any single-column row-add combo is the same matrix but with opposite signs on the subdiagonal entries:

$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 2& 1& 0& 0\\ 3& 0& 1& 0\\ 4& 0& 0& 1\end{array}\right)\left(\begin{array}{rrrr}1& 0& 0& 0\\ -2& 1& 0& 0\\ -3& 0& 1& 0\\ -4& 0& 0& 1\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).$$

Computing LU factorizations

Here is the plan:

1. Postulate hypothetically that $A=LU=(L_1\cdots L_{n-1})U$, where $L_i$ are single-column row-add combos which combine to form the hypothetical $L$. We know $A$ and we seek $L$ and $U$.
2. Act upon $A$ and $LU$ on the left simultaneously by $L_1^{-1}$, then by $L_2^{-1}$, etc., through $L_{n-1}^{-1}$. We must invent $L_1^{-1}$, $L_2^{-1}$ etc. as we go along:
   1. Choose $L_1^{-1}$ such that $L_1^{-1}A$ has all zeros in the first column on the subdiagonal: In each nonzero subdiagonal entry of $L_1^{-1}$, put $\lambda$ such that row-add of the top entry into entry $2$ eliminates entry $2$, and so on to create the entire first column of $L_1^{-1}$ which acts to clear the first column of $A$.
   2. Continue with $L_2^{-1}$ chosen such that $L_2^{-1}(L_1^{-1}A)$ now has zeros in the first and second columns, and $L_2^{-1}$ has subdiagonal entries $\lambda$ in the second column.
   3. Continue through $L_{n-1}^{-1}$ and obtain $(L_1^{-1}\cdots L_{n-1}^{-1})A=U$. By flipping signs on $L_i$ for $i=1,\dots ,n-1$, we now know $L_i$ and simply combine the subdiagonal columns of all $L_i$ to create the whole $L$.
   4. Notice! The matrix $L_i$ is exactly the unique single-column row-add combo matrix which has $i^{\text{th}}$ column cleared by whatever clears the $i^{\text{th}}$ column of $A$. In other words, we cook up $L_i$ so that the same elimination works both on $A$ and on $L_i$, the difference being that $L_i$ has a $1$ at the top entry by definition, while the corresponding entry of $A$ may have other values.

In other words, we find the matrices $L_1,\dots ,L_{n-1}$ by backwards reasoning: row-add eliminations performed on $A$ should simultaneously effect eliminations on $L_1,\dots ,L_{n-1}$, converting them each to the identity matrix $I_n$. If this works, we know we have the $L_i$ ingredients inside the factorization $A=LU=(L_1\cdots L_{n-1})U$.

Pivoting

The algorithm presented above for LU factorization can fail down because of zeros. For example, consider $A$:

$$A=\left(\begin{array}{rrrr}0& 1& 1& 0\\ 2& -3& 5& -8\\ 6& -4& 17& -26\\ 8& -6& 5& -45\end{array}\right).$$

We cannot even begin to eliminate $2,6,8$ using row-add steps. In this situation, it is best to pivot to a different top row by using a row-swap permutation. The new top row is sometimes called the pivot row, and the new entry used to eliminate others by row-adding may be called the pivot entry or simply the pivot.

Write $P$ for the row-swap matrix that swaps rows $1$ and $2$. After pivoting, we have a new matrix:

$$PA=\left(\begin{array}{rrrr}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{rrrr}0& 1& 1& 0\\ 2& -3& 5& -8\\ 6& -4& 17& -26\\ 8& -6& 5& -45\end{array}\right)=\left(\begin{array}{rrrr}2& -3& 5& -8\\ 0& 1& 1& 0\\ 6& -4& 17& -26\\ 8& -6& 5& -45\end{array}\right).$$

The new matrix has an LU factorization now:

$$\left(\begin{array}{rrrr}2& -3& 5& -8\\ 0& 1& 1& 0\\ 6& -4& 17& -26\\ 8& -6& 5& -45\end{array}\right)=\left(\begin{array}{rrrr}1& 0& 0& 0\\ 0& 1& 0& 0\\ 3& 5& 1& 0\\ 4& 6& 7& 1\end{array}\right)\left(\begin{array}{rrrr}2& -3& 5& -8\\ 0& 1& 1& 0\\ 0& 0& -3& -2\\ 0& 0& 0& 1\end{array}\right).$$

It is appropriate to say that “$A$ has an LU factorization with pivots”, or to say that “we have a factorization $PA=LU$.” We do not say simply that “$A$ has an LU factorization” if pivots are required.

Every matrix $A$ has a factorization of the form $PA=LU$.

Numerical methods

For numerical stability, it is preferable to always pivot to the largest entry in a given column before performing elimination. This technique is called “partial pivoting.”

(“Complete pivoting” refers to a technique where elimination is performed to clear any uncleared column by choosing the largest remaining entry across the matrix. It is computationally expensive and therefore rarely used.)

The general procedure for computing $PA=LU$ factorizations is fairly complex because the natural algorithm interlaces elimination (using factors $L_i^{-1}$) with permutations for the pivots (using row-swap matrices $P_i$). You are not required to master this procedure for exams.

The general process of computing $U$ by successive elimination of terms with pivots, clearing out one column at a time, is called row reduction (see below). You are required to master row-reduction for exams.

The difference between row-reduction and the procedure for $PA=LU$ is that in the latter we do not keep track of $L_i$ and $P_i$ and do not calculate $L$ and $P$, we simply perform the row-add operations on $A$ until the $L$ factor inside $A$ is gone (by clearing columns) and the result is $U$.

Here is the theory for producing $PA=LU$ in general:

Factorization $PA=LU$

For $i=1,\dots ,n$, work through row $i$ of $A$:

1. Choose a row of $A$ to pivot into the first row, defining the corresponding row-swap matrix $P_1$. After swapping we have $P_{1}A$. (For numerical stability, choose the row of $A$ with largest entry in the first column; but in principle, can choose any row with nonzero entry in the first column.)
2. Construct $L_1^{-1}$ to eliminate entries in row $2$ and below, clearing column $1$ of $P_{1}A$. Obtain $L_1^{-1}{P}_{1}A$.
3. Proceed to row $2$: choose a pivot row and define row-swap $P_2$ and obtain $P_2{L}_1^{-1}{P}_{1}A$.
4. Construct $L_2^{-1}$ to clear column $2$ of $P_2{L}_1^{-1}{P}_{1}A$.
5. Repeat with $P_i$ a row-swap putting a chosen pivot from below row $i$ into row $i$.
6. Repeat with $L_i^{-1}$ clearing the $i^{\text{th}}$ column below the $i^{\text{th}}$ row.
7. End when $i=n$.

In this way we build a sequence:

$$L_{n-1}^{-1}{P}_{n-1}{L}_{n-2}^{-1}\dots P_2{L}_1^{-1}{P}_{1}A=U.$$

If we didn’t have all the $P_i$ matrices interlaced between the $L_i$’s, this would be $L_{n-1}^{-1}\dots L_1^{-1}A=U$, and therefore $A=L_1\dots L_{n}U$ as in the ordinary LU factorization. So what we want to do here is to extract the $P_i$’s out to the right and group them into $P$. This is the hardest part!

There is a trick we can use to make this happen. It turns out that $P_j{L}_i^{-1}{P}_j^{-1}$ (for $i<j$) is the same as $L_i^{-1}$, except that the nonzero stack of entries undergoes the row-swap of $P_j$. Using this fact, together with the fact that $P_j^{-1}=P_j$ because $P_j$ swaps two rows, we can rewrite our sequence in terms of other single-column matrices called $L_i^{\prime }$. For example, in the case $n=4$, we could write:

$$L_3^{\prime -1}=L_3,L_2^{\prime -1}=P_3{L}_2^{-1}{P}_3,L_1^{\prime -1}=P_3{P}_2{L}_1^{-1}{P}_2{P}_3.$$

From these definitions, it follows that

$$\left(L_3^{\prime -1}{L}_2^{\prime -1}{L}_1^{\prime -1}\right)\left(P_3{P}_2{P}_1\right)A=L_3^{-1}{P}_3{L}_2^{-1}{P}_2{L}_1^{-1}{P}_{1}A=U$$

and thus

$$\begin{array}{rl}\left(P_3{P}_2{P}_1\right)A& =\left(L_1^{\prime }{L}_2^{\prime }{L}_3^{\prime }\right)U\\ PA& =LU,\end{array}$$

where we define $P=P_3{P}_2{P}_1$ and $L=L_1^{\prime }{L}_2^{\prime }{L}_3^{\prime }$.

Here is a concise computational algorithm to perform these steps on a computer:

Algorithm on entries for $PA=LU$ factorization

Initialize: $P=L=I_n,U=A$.
For $j=1$ to $n-1$ :

• Choose $i\ge j$ with maximum (or simply nonzero) $|u_{ij}|$.
  – That is: $i$ is the new pivot row index.
• Exchange $u_{j,j:n}\leftrightarrow u_{i,j:n}$.
  – That is: swap current row $j$ with new pivot row $i$ on future columns of $U$.
• Exchange $p_{j,:}\leftrightarrow p_{i,:}$.
  – That is: swap current row $j$ with new pivot row $i$ on all of $P$ to create $P_j{P}_{j-1}\cdots P_1$.
• Exchange ${\ell }_{j,1:j-1}\leftrightarrow {\ell }_{i,1:j-1}$.
  – That is: swap current row $j$ with new pivot row $i$ on the subdiagonal entries of $L$ to create $P_j{L}_1^{\prime }\cdots L_{j-1}^{\prime }{P}_j$.
• For $k=j+1$ to $n$ :
  • Set ${\ell }_{kj}=u_{kj}/u_{jj}$
    – That is: write next subdiagonal stack of $L$: add terms to $P_j{L}_1^{\prime }\cdots L_{j-1}^{\prime }{P}_j$ to create $L_1^{\prime }\cdots L_j^{\prime }$.
  • Set $u_{k,j:n}=u_{k,j:n}-{\ell }_{kj}{u}_{j,j:n}$
    – That is: clear subdiagonals of column $j$ of $U$ by row-adding.

To understand the step of setting ${\ell }_{kj}$ for $k>j$, it is useful to make the following observation:

$$(P_j{L}_1^{\prime }\cdots L_{j-1}^{\prime }{P}_j)L_j=L_1^{\prime }\cdots L_j^{\prime }.$$

This notation is a little misleading though, because $L_1^{\prime }$ on the RHS has subdiagonals permuted by $P_j$, whereas $L_1^{\prime }$ on the LHS does not. The factor $P_j{L}_1^{\prime }{P}_j$ will, though. Therefore, the meaning of $L_i^{\prime }$ changes at each step in the comments on the above algorithm.

Row-echelon form

When $U$ has a number of zeros cropping up, we can sometimes use pivoting to clear out more columns. For example, consider the matrices:

$$A=\left(\begin{array}{rccc}0& 0& 1& 0\\ 0& 2& 3& 4\\ 0& 0& 4& 5\\ 0& 0& 0& 1\end{array}\right),B=\left(\begin{array}{rccc}0& 2& 3& 4\\ 0& 0& 1& 0\\ 0& 0& 4& 5\\ 0& 0& 0& 1\end{array}\right),C=\left(\begin{array}{rccc}0& 2& 3& 4\\ 0& 0& 1& 0\\ 0& 0& 0& 5\\ 0& 0& 0& 0\end{array}\right).$$

Here $A$ is technically upper-triangular and satisfies the rules for $U$ in an LU factorization. However, it is possible to perform a pivot and obtain $B$, which can be further reduced by two elimination steps to obtain $C$.

The form of $C$ is called row-echelon form, or sometimes just echelon form. To be precise: in row-echelon form, the leading entry (left-most nonzero entry) in each row is somewhere to the right of the leading entries in the rows above, and rows of only zeros are at the bottom.

The combined process of creating $U$ as in LU decomposition and then continuing pivots and eliminations to put $U$ in row-echelon form is called row reduction.

• Row-echelon form for $U$ is very useful for inversion problems (see below).
  – For inversion problems, we will perform row reduction simultaneously on $A$ and on a vector $\mathbf{b}$ or matrix $B$ upon which we are inverting the action of $A$.
• Row-echelon form for $U$ is not needed for finding determinants (see below).
  – Indeed, row-echelon $U$ can only differ from other triangular $U$ when the determinant is zero.