3.2: LU Decomposition
( \newcommand{\kernel}{\mathrm{null}\,}\)
The process of Gaussian Elimination also results in the factoring of the matrix A to
A=LU,
where L is a lower triangular matrix and U is an upper triangular matrix. Using the same matrix A as in the last section, we show how this factorization is realized. We have
(−32−16−673−44)→(−32−10−250−23)=M1 A
where
M1 A=(100210101)(−32−16−673−44)=(−32−10−250−23)
Note that the matrix M1 performs row elimination on the first column. Two times the first row is added to the second row and one times the first row is added to the third row. The entries of the column of M1 come from 2=−(6/−3) and 1=−(3/−3) as required for row elimination. The number −3 is called the pivot.
The next step is
(−32−10−250−23)→(−32−10−2500−2)=M2(M1 A)
where
M2(M1 A)=(1000100−11)(−32−10−250−23)=(−32−10−2500−2)
Here, M2 multiplies the second row by −1=−(−2/−2) and adds it to the third row. The pivot is −2.
We now have
M2M1 A=U
or
A=M−11M−12U
The inverse matrices are easy to find. The matrix M1 multiples the first row by 2 and adds it to the second row, and multiplies the first row by 1 and adds it to the third row. To invert these operations, we need to multiply the first row by −2 and add it to the second row, and multiply the first row by −1 and add it to the third row. To check, with
M1M−11=I
we have
(100210101)(100−210−101)=(100010001)
Similarly,
M−12=(100010011)
Therefore,
L=M−11M−12
is given by
L=(100−210−101)(100010011)=(100−210−111)
which is lower triangular. The off-diagonal elements of M−11 and M−12 are simply combined to form L. Our LU decomposition is therefore
(−32−16−673−44)=(100−210−111)(−32−10−2500−2).
Another nice feature of the LU decomposition is that it can be done by overwriting A, therefore saving memory if the matrix A is very large.
The LU decomposition is useful when one needs to solve Ax=b for x when A is fixed and there are many different b’s. First one determines L and U using Gaussian elimination. Then one writes
(LU)x=L(Ux)=b.
We let
y=Ux
and first solve
Ly=b
for y by forward substitution. We then solve
Ux=y
for x by backward substitution. When we count operations, we will see that solving (LU)x=b is significantly faster once L and U are in hand than solving Ax=b directly by Gaussian elimination.
We now illustrate the solution of LUx=b using our previous example, where
L=(100−210−111),U=(−32−10−2500−2),b=(−1−7−6)
With y=Ux, we first solve Ly=b, that is
(100−210−111)(y1y2y3)=(−1−7−6)
Using forward substitution
y1=−1y2=−7+2y1=−9y3=−6+y1−y2=2
We now solve Ux=y, that is
(−32−10−2500−2)(x1x2x3)=(−1−92)
Using backward substitution,
−2x3=2→x3=−1−2x2=−9−5x3=−4→x2=2−3x1=−1−2x2+x3=−6→x1=2
and we have once again determined
(x1x2x3)=(22−1)