*Engineering, Tech, TeX, and Numbers. Lots of Numbers.*

- 2014.01.31: Latent Semantic Analysis
- 2014.01.25: On Cutting Firewood
- 2014.01.08: Dragons!
- 2013.10.13: A Haiku For Winter
- 2013.09.28: ShareLaTeX: TeX from anywhere!
- 2013.09.23: Google Trends Scipt

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*.

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(n ^{2})* operations compared to the

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}

**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*.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}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}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).

**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