LU Factorization

Problem Statement:

Given a n-by-n matrix A. The LU Factorization of A refers to writing A = L * U, where L and U are lower-triangular matrix and upper-triangular matrices respectively.

The product L * U equals the matrix A shown below.

How does one find matrices L and U from A?

Example:

l21 = a21/a11 = ¼

l31 = a31/a11 = ½

l41 = a41/a11 = ¾

l32 = a32/a22 = ¾
l42 = a42/a22 = ½

l43 = a43/a33 = ¼

L * U = A (check this)

How do we program this?

Modify Gaussian Elimination to develop an algorithm to compute triangular matrices L and U such that A = L * U. (see board).

Vectorized Gauss:

for p = 1:n-1

m = -A(p+1:n,p)/A(p,p);

A(p+1:n,p+1:n) = A(p+1:n,p+1:n) + m*A(p,p+1:n);

b(p+1:n) = b(p+1:n) + m*b(p);

end

Example:

>> A

A =

4             12              8              4

1              7             18              9

2              9             20             20

3             11             15             14

>> [L,U]=lu(A)

L =

1              0              0              0

1/4            1              0              0

1/2            3/4            1              0

3/4            1/2            1/4            1

U =

4             12              8              4

0              4             16              8

0              0              4             12

0              0              0              4

>> L*U

ans =

4             12              8              4

1              7             18              9

2              9             20             20

3             11             15             14

Solving the system Ax = b using LU-Factorization

Matrix A is a known n-by-n matrix and b is an n-by-1 vector.

Write A*x = b as

L*U*x = b    (1)

Let

y = U*x      (2)

Substitute equation (2) into equation (1) above to get

L*y = b      (3)

Solve equation (3) for y using Forward Substitution. Finally, solve equation (2) for x using Back Substitution.

Example:

Call this equation 1:

Let

Call this equation 2. Substitute equation 2 into equation 1.

y1 = 24

y2 = 16

y3 = 12

y4 = 4

Solve for x

X1 = -1

X2 =  2

X3 =  0

X4 =  1

Gaussian Elimination requires O(n3) arithmetic operations and Back and Forward Substitution require O(n2) arithmetic operations. Suppose you want to solve m≥1 matrix equations of the form

Ax = bj, where bj is an n-by-1 vector for j = 1, 2, …, m and A is a fixed n-by-n coefficient matrix. The direct method (Gaussian Elimination followed by Back Substitution) to solve these m systems would require O(mn3) arithmetic operations. Using LU Factorization to solve this same system requires n3 + mn2 arithmetic operations which is a great improvement for large m.