LU Decomposition

What is LU Decomposition?

LU Decomposition is simply the process of factoring a matrix A into an upper triangular matrix U, and a lower triangular matrix L so that L×U = A.

Ok, well, why should I care?

Many engineering applications involve solving large sets of linear equations, of the form Ax = b. Typically, the matrix A describes the physical system, the vector b describes the stimulus applied to the system, and the vector x is the response of the system.

In situations where we need to apply many different stimuli to a system, LU Decomposition allows us to essentially solve half of the problem. By cleverly creating L and U, we can later solve the system of equation using O(n2) operations compared to the O(n3) required for Gaussian Elimination.

How do I do this black magic?!

There are two common ways to perform the decomposition, but the most common are Doolittle Decomposition and Crout Decomposition. We'll illustrate these methods using this matrix:

A = \begin{bmatrix} 3 & 2 & 9 \\ 11 & 5 & 1 \\ 7 & 9 & 3 \end{bmatrix}
1. Doolittle Decomposition

Doolittle Decomposition is based on Guassian Elimination. As we place the matrix A in row echelon form, we save the factors that we used to eliminate matrix elements in L. Once A is in row echelon form, it is the pertinent upper triangular matrix, U.

1. We automatically know several of elements in both L and U.

A = \begin{bmatrix} 3 & 2 & 9 \\ 11 & 5 & 1 \\ 7 & 9 & 3 \end{bmatrix}\textbf{ }U = \begin{bmatrix} 3 & 2 & 9 \\ 0 & \textbf{-} & \textbf{-} \\ 0 & 0 & \textbf{-} \end{bmatrix}\textbf{ }L = \begin{bmatrix} 1 & 0 & 0 \\ \textbf{-} & 1 & 0 \\ \textbf{-} & \textbf{-} & 1 \end{bmatrix}
2. Next, we perform gaussian elimination to place the second row in row echelon form:

R_2 - \left ( \frac{a_{21}}{a_{11}} \right )R_1 \rightarrow R_2\textbf{, }l_{21} = \frac{a_{21}}{a_{11}}
R_3 - \left ( \frac{a_{31}}{a_{11}} \right )R_1 \rightarrow R_3\textbf{, }l_{31} = \frac{a_{21}}{a_{11}}
A = \begin{bmatrix} 3 & 2 & 9 \\ 0 & -2.333 & -32 \\ 0 & -4.333 & -18 \end{bmatrix}\textbf{ }
U = \begin{bmatrix} 3 & 2 & 9 \\ 0 & -2.333 & -32 \\ 0 & 0 & \textbf{-} \end{bmatrix}\textbf{ }L = \begin{bmatrix} 1 & 0 & 0 \\ 3.667 & 1 & 0 \\ 2.333 & \textbf{-} & 1 \end{bmatrix}
3. Finally, we place the third row in row echelon form:

R_3 - \left ( \frac{a_{32}}{a_{22}} \right )R_2 \rightarrow R_3\textbf{, }l_{32} = \frac{a_{32}}{a_{22}}
A = \begin{bmatrix} 3 & 2 & 9 \\ 0 & -2.333 & -32 \\ 0 & 0 & -77.428 \end{bmatrix}\textbf{ }
U = \begin{bmatrix} 3 & 2 & 9 \\ 0 & -2.333 & -32 \\ 0 & 0 & -77.428 \end{bmatrix}\textbf{ }L = \begin{bmatrix} 1 & 0 & 0 \\ 3.667 & 1 & 0 \\ 2.333 & -1.857 & 1 \end{bmatrix}

Which yields our final decomposition. We can verify the decomposition by calculating L×U and comparing it to A (they should be equivalent).

2. Crout Decomposition

Doolittle Decomposition is a by-product of a method that we're familiar with: Gaussian Elimination. On the other hand, Crout Decomposition is rooted in the definition of the factorization itself. For this method, our decomposition takes the form:

\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix}=\begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \times \begin{bmatrix} 1 & u_{12} & u_{13} \\ 0 & 1 & u_{23} \\ 0 & 0 & 1 \end{bmatrix}

If we carry out this matrix multiplication, we end up with 9 equations for the 9 unknowns in the U and L matrices. These equations are:

 a_{11} = l_{11} a_{12} = l_{11}u_{12} a_{13} = l_{11}u_{13} a_{21} = l_{21} a_{22} = l_{21}u_{12} + l_{22} a_{23} = l_{21}u_{13} + l_{22}u_{23} a_{31} = l_{31} a_{32} = l_{31}u_{12} + l_{32} a_{33} = l_{31}u_{13} + l_{32}u_{33} + l_{33}

This set of equations can be solved for the unknowns, yielding the factorization:

L = \begin{bmatrix} 3 & 0 & 0 \\ 11 & -2.333 & 0 \\ 7 & 4.333 & -77.429 \end{bmatrix}\textbf{ }U = \begin{bmatrix} 1 & 0.667 & 3 \\ 0 & 1 & 13.714 \\ 0 & 0 & 1 \end{bmatrix}

Using either of these decompositions, we can finish solving the system Ax = b by calculating:

x = U^{-1}L^{-1}b

This is typically broken down into two steps. First we solve for an intermediate vector, D using backwards substitution:

D = L^{-1}b

And then we solve for x using forward substitution:

x = U^{-1}D

Both of these methods can easily be extended to larger matrices. Here is a MATLAB function that implements Doolittle Decomposition and another that implements Crout Decomposition.

by Chris McComb