LU Decomposition

Mathematics

The LU decomposition is one of the most basic tenets of numerical linear algebra. It is used to solve the linear system of equations Ax=b

While it is possible to solve this system by using an explicit inverse of A, i.e. x = A^-1b, this is almost never a good idea for efficiency reasons. A numerical analyst would never do that. Instead, a factorization of A is more effective. One standard factorization is to factorize A=LU where L is a lower triangular matrix and U is an upper triangular matrix. Then Ax=b may be solved as follows:

Ax = b
(LU)x = b
L(Ux) = b

now let c denote Ux and solve the system

Lc = b

Single L is lower triangular, this is a simple forward substitution. Once we know c we then solve

Ux = c

via backward substitution and we are done. This method has the advantage that multiple right hand side solves can be done quite quickly.

Other important decompositions include the Choleski Decomposition and the QR Decomposition.

LU Decomposition is also sometimes called LU factorization. It seems to pop up frequently in college-level classes in math, numerical methods and computer science, but far fewer are examples sighted. And since there are more explanations of just the theory (as above) than you can shake a stick at, I'll only bore you with a sample work-through.

Sorry about the <PRE>s, folks, but this is matrix math, after all.

---

We start with the matrix A=

|   2    1    0   -4  |
|   1    0  0.25  -1  |
|  -2  -1.1 0.25  6.2 |
|   4   2.2  0.3 -2.4 |

with solutions b=

|  -3  |
| -1.5 |
|  5.6 |
|  2.2 |

Now, we're looking to create two matrices, U and L. To do so, we "zero out" entries below the diagonals using row operations to add a multiple of one row to another.

Blah blah blah. Here's what that means: We want an upper triangular matrix which looks like this:

| X X X X |
| 0 X X X |
| 0 0 X X |
| 0 0 0 X |

This is what U is supposed to look like when we're done. To have a matrix in this form allows back substitution which means you can find out what the 'bottom' one is and work backward from there. But getting back on track, we want to add a multiple of the first row to the second to make the first column of the second row a zero. We'll use the multiple of one-half, so we subtract one half of the first row from the second, giving us:

|   2    1    0   -4  |
|   0  -0.5 0.25   1  |
|  -2  -1.1 0.25  6.2 |
|   4   2.2  0.3 -2.4 |

for U and as such we insert 0.5 into our L matrix:

|  1   0   0   0  |
| 0.5  1   0   0  |
|  ?   ?   1   0  |
|  ?   ?   ?   1  |

We subtract -1 times the first row from the third (i.e. add it), and 2 times it from the fourth, giving us U₁=

|   2    1    0   -4  |
|   0  -0.5 0.25   1  |
|   0  -0.1 0.25  2.2 |
|   0   0.2  0.3  5.6 |

and L₁=

|  1   0   0   0  |
| 0.5  1   0   0  |
| -1   ?   1   0  |
|  2   ?   ?   1  |

Let's zero out the second column. To do this, we utilize the second row, subtracting 0.2 times it from the third and 0.4 times it from the fourth, giving us U₂=

|   2    1    0   -4  |
|   0  -0.5 0.25   1  |
|   0    0   0.2   2  |
|   0    0   0.4   6  |

and L₂=

|  1   0   0   0  |
| 0.5  1   0   0  |
| -1  0.2  1   0  |
|  2  0.4  ?   1  |

Leaving us with the third column to zero out, a multiple of twice the third row subtracted from the fourth, giving us a final U=

|   2    1    0   -4  |
|   0  -0.5 0.25   1  |
|   0    0   0.2   2  |
|   0    0    0    2  |

and L=

|  1   0   0   0  |
| 0.5  1   0   0  |
| -1  0.2  1   0  |
|  2  0.4  2   1  |

So what do you do with them now? We solve Ax=b by making it LUx=b (after all, LU=A) and then find a vector z from Lz=b and then solving Ux=z and (I'll save you the further matrix stuff) then we get x=

| 0.5 |
|  2  |
| -2  |
| 1.5 |

---

And that's it. If you've got more questions, if I've made a mistake, or anything like that, /msg me.

node your homework