MATH/APPM 4650
Class notes, April 20
Section 6.5: LU Factorization


Recall that it is easy (O(n²) computations) to solve Ux = b for x when U is an upper-triangular matrix.

% Solve Ux=b for x, when U is an upper triangular matrix with no zeros on the diagonal.
x[n] = b[n] / U[n,n]
for i from 1 to n-1 do
   x[n-i] = b[n-i]
   for j from n-i+1 to n do
      x[n-i] = x[n-i] - U[n-i,j]*x[j]
   end do
   x[n-i] = x[n-i] / U[n-i,n-i]
end do

Similarly it is easy to solve Lx = b for x if L is a lower-triangular matrix.

% Solve Lx=b for x, when L is a lower triangular matrix with no zeros on the diagonal.
x[1] = b[1] / L[1,1]
for i from 1 to n-1 do
   x[i] = b[i]
   for j from 1 to i-1 do
      x[i] = x[i] - L[i,j]*x[j]
   end do
   x[i] = x[i] / L[i,i]
end do

In general, most matrices can be written in the form A = LU, for some lower-triangular matrix L and upper-triangular matrix U.

So if we want to solve Ax = b, we can write LUx = b.

Letting y = Ux, we can think of this as two matrix equations: Ly = b and Ux = y.

First solve Ly = b for y. Then solve Ux = y for x.

Both solutions take O(n²) operations, so we have solved a matrix using only O(n²) operations.

Of course, the problem is how to find L and U given A. This takes O(n³) operations.

No free lunch!



If we look at exactly how we did the row reduction process, we can think of each step of Gaussian elimination as being a matrix multiplication.

1. "Multiply the second row of A by 5" is the same as "Replace A with M₁A, where M₁ is the following matrix."
   

   / 1	0 \
   \ 0	5 /

2. "Add 7 times row one to row two" is the same as "Replace A with M₂A, where M₂ is the following matrix."
   

   / 1	0 \
   \ 7	1 /

So for example, with a 3x3 matrix, to obtain the upper triangular form, we would do three of the second kind of operation.
U = M₃M₂M₁A where

	/   1	0	0 \
M1   =	| -m12	1	0 |
	\   0	0	1 /

	/     1	0	0 \
M2   =	|      0	1	0 |
	\ -m13	0	1 /

	/ 1	0	0 \
M3   =	| 0	1	0 |
	\ 0	-m23	1 /

(Important: notice that m₂₃ is the multiplier for the reduced matrix; that is, you clear out the first column first, then compute a₃₂/a₂₂ using the new values of a.

This means that A = M₁^-1M₂^-1M₃^-1U

But notice that

	/ 1	0	0 \
M1-1M2-1M3-1    =	| m12	1	0 |
	\ m13	m23	1 /

So we have factored A = LU, where L is a lower-triangular matrix (with 1s on the diagonal) and U is an upper-triangular matrix.

This can always be done, as long as we never have to switch rows.

Switching rows is a matrix like

/ 0	1 \
\ 1	0 /

and this messes up the lower-triangular property.

So the algorithm is Doolittle's method, slightly different from the book for clarity:
% Given a matrix A, factor into LU form.
% First do Gaussian elimination, saving all the multipliers to put in the off-diagonal components of L.
L = identity matrix
for k from 1 to n-1 do
   for i from k+1 to n do
      L[i,k] = a[i,k] / a[k,k]
      for j from k+1 to n do
         a[i,j] = a[i,j] - L[i,k]*a[k,j]
      end do
   end do
end do

% Save the new entries of A as U.
U = zero matrix
for i from 1 to n do
   for j from i to n do
      U[i,j] = a[i,j]
   end do
end do

Why are we doing this?

If we want to solve the equation Ax=b multiple times, for the same A but different b, we should not use the full Gaussian elimination algorithm each time.
Do the LU factorization once, then solve Ly=b, Ux=y repeatedly.

Note: A factorization A = LU is never unique. We have to specify the diagonals elements on one of the elements. (Either L or U but not both.) It's usually easiest to demand that one of them have 1s on the diagonal, since Gaussian elimination gives you that automatically.

If we had wanted to specify that the diagonal elements of U are all 1s, we would have done Gaussian elimination on all the terms above the diagonal.

---

What if we have to switch rows?

A matrix to switch row one and row two of a 3x3 matrix looks like

	/ 0	1	0 \
P =	| 1	0	0 |
	\ 0	0	1 /

In general, a permutation matrix is a matrix that has exactly one '1' in each row and each column.

There are n! permutation matrices of size nxn.

They all satisfy P^T = P^-1.

If you do Gaussian elimination using the algorithm above, and you ever need to switch a row, you can incorporate that switch into a permutation matrix. (See board for example.)