LU Factorization with MATLAB | MATLAB Tutorial

LU Factorization with MATLAB

LU Factorization

[A]{x}={b}                                                                         1

Equation 1 can be rearranged to give

[A]{x}-{b}=0                                                                     2

Suppose that above eq. 2 could be expressed as an upper triangular system. For example, a 3×3 system:

Eq. 3 can also be expressed in matrix notation and rearranged to give

[U]{x}-{d}=0                                                                     4

Now assume that there is a lower diagonal matrix with 1’s on the diagonal.

that has the property that when Eq. 4 is premultiplied by it, Eq. 2 is the result. That is,

[L]{[U]{x}-{d}}=[A]{x}-{b}                                            6

If this equation holds, it follows from the rules for matrix multiplication that

[L][U]=[A]                                                                          7

and

[L]{d}={b}                                                                          8

A two step strategy(Fig.1) for obtaining solutions can be based on Eqs. 3, 7, and 8:

  1. LU factorization step. [A] is factored or “decomposed” into lower [L] and upper [U] triangular matrices.
  2. Substitution step. [L] and [U] are used to determine a solution {x} for a right-hand side {b}. This step itself consists of two steps. First, Eq. 8 is used to generate an intermediate vector {d} by forward substitution. Then, the result is substituted into Eq. 3 which can be solved by back substitution for {x}.

Figure 1: The steps in LU factorization

Gauss Elimination as LU Factorization

Gauss elimination can be used to decompose [A] into [L] and [U]. This can be easily seen for [U], which is a direct product of the forward elimination. Recall that the forward-elimination step is intended to reduce the original coefficient matrix [A] to the form

which is in the desired upper triangular format, the matrix [L] is also produced during the step. This can be readily illustrated for a three-equation system,

The first step in Gauss elimination is to multiply row 1 by the factor

 

and subtract the result from the second row to eliminate a21 . Similarly, row 1 is multiplied by

and the result subtracted from the third row to eliminate a31 . The final step is to multiply the modified second row by

and subtract the result from the third row to eliminate .

After elimination, the [A] matrix can therefore be written as

 

This matrix, in fact, represents an efficient storage of the LU factorization of [A] [A] -> [L][U]

where

and

Problem 1:  Derive LU factorization based on the Gauss elimination,

Coefficient matrix:

After forward elimination, the following upper triangular matrix was obtained:

The elements a21 and a31 were eliminated by using the factors

and the element a32 was eliminated by using the factor

Thus, the lower triangular matrix is

Consequently, the LU factorization is

This result can be verified by performing the multiplication of [L][U] to give

After the matrix is decomposed, a solution can be generated for a particular right-hand side vector {b}. This is done in two steps. First, a forward-substitution step is executed by solving Eq. 8 for {d}. At the end of this step, the right-hand side will be in the same state that it would have been had we performed forward manipulation on [A] and {b} simultaneously.

The forward-substitution step can be represented concisely as

The second step then merely amounts to implementing back substitution to solve Eq. 3.

Problem 2: Complete the problem 1 by generating the final solution with forward and back substitution.

the forward – elimination phase of conventional Gauss elimination resulted in

The forward-substitution phase is implemented by applying Eq. 8:

We can solve the first equation for d1 = 7.85, which can be substituted into the second equation to solve for

d2=-19.3-0.03333(7.85)=-19.5617

Both d1 and d2 can be substituted into the third equation to give

d3=71.4-0.1(7.85)+0.02713(-19.5617)=70.0843

Thus,

This can be substituted into Eq. 3, [U]{x} = {d}:

Matlab m-file:

function x = LUGauss(A,b)
% input:
% A = coefficient matrix
% b = right hand side vector
% output:
% x = solution vector

[m,n]=size(A);
if m~=n
 error('Matrix A must be square'); 
end

[L,U] = lu(A); % LU factorization
disp('L = ');
disp(L)
disp('U = ');
disp(U)
d = L\b;
x = U\d; % solution
disp('x = ');
disp(x)

The M-file can be used to solve the problem 1 and 2.

>> A = [3 -.1 -.2;.1 7 -.3;.3 -.2 10];

b = [7.85; -19.3; 71.4];

>> LUGauss(A,b);

L =

1.0000         0               0

0.0333    1.0000         0

0.1000   -0.0271    1.0000

U =

3.0000   -0.1000   -0.2000

0           7.0033   -0.2933

0                0         10.0120

x =

3.0000

-2.5000

7.0000

LU Factorization with Pivoting

Just as for standard Gauss elimination, partial pivoting is necessary to obtain reliable solutions with LU factorization. One way to do this involves using a permutation matrix. The approach consists of the following steps:

  1. Elimination. The LU factorization with pivoting of a matrix [A] can be represented in matrix form as
[P][A] = [L][U]

The upper triangular matrix, [U], is generated by elimination with partial pivoting, while storing the multiplier factors in [L] and employing the permutation matrix, [P], to keep track of the row switches.

  1. Forward substitution. The matrices [L] and [P] are used to perform the elimination step with pivoting on {b} in order to generate the intermediate right-hand-side vector, {d}. This step can be represented concisely as the solution of the following matrix formulation:
[L]{d} = [P]{b}

  1. Back substitution. The final solution is generated in the same fashion as done previously for Gauss elimination. This step can also be represented concisely as the solution of the matrix formulation:
[U]{x} = {d}

Problem 3: Compute the LU factorization and find the solution for the system,

Before elimination, we set up the initial permutation matrix:

We immediately see that pivoting is necessary, so prior to elimination we switch the rows:

At the same time, we keep track of the pivot by switching the rows of the permutation matrix:

We then eliminate a21 by subtracting the factor l21 = a21/a11 = 0.0003/1 = 0.0003 from the second row of A. In so doing, we compute that the new value of a22 = 3 – 0.0003(1) = 2.9997. Thus, the elimination step is complete with the result: 

Before implementing forward substitution, the permutation matrix is used to reorder the right-hand-side vector to reflect the pivots as in

Then, forward substitution is applied as in

which can be solved for d1 = 1 and d2 = 2.0001 − 0.0003(1) = 1.9998. At this point, the system is

Applying back substitution gives the final result:

The LU factorization algorithm requires the same total flops as for Gauss elimination. The only difference is that a little less effort is expended in the factorization phase since the operations are not applied to the right-hand side. Conversely, the substitution phase takes a little more effort.

Matlab m-file:


function x = LUPivot(A,b)

% input:
% A = coefficient matrix
% b = right hand side vector
% output:
% x = solution vector

[m,n]=size(A);
if m~=n
 error('Matrix A must be square');
end
L = eye(n);
P = L;
U = A;
for k=1:n
 [pivot m]=max(abs(U(k:n,k)));
 m = m+k-1;
end
end

 

The M-file can be used to solve the problem 3.

>> A = [0.0003 3.000;1.0000 1.0000];

>> b = [2.0001 ; 1.0000];

>> LUPivot(A,b);

L =

1.0000         0

0.0003    1.0000

U =

1.0000    1.0000

0    2.9997

x =

0.3333

0.6667


I am a Certified MATLAB Associate. I love MATLAB Coding & Enjoy Programming. I can make anything & everything possible with MATLAB.