LU Decomposition and Matrix Inversion

The Islamic University of Gaza
Faculty of Engineering
Civil Engineering Department
Numerical Analysis
ECIV 3306
Chapter 10
LU Decomposition and Matrix Inversion

Introduction
• Gauss elimination solves [A] {x} ={B}
• It becomes insufficient when solving these

equations for different values of {B}
• LU decomposition works on the matrix [A] and the

vector {B} separately.
• LU decomposition is very useful when the vector of

variables {x} is estimated for different parameter
vectors {B} since the forward elimination process is
not performed on {B}.
LU Decomposition
If:
L: lower triangular matrix
U: upper triangular matrix
Then,
[A]{X}={B} can be decomposed into two
matrices [L] and [U] such that:
1. [L][U] = [A]  ([L][U]){X} = {B}
LU Decomposition
Consider:
[U]{X} = {D}
So, [L]{D} = {B}
2. [L]{D} = {B} is used to generate an
intermediate vector {D} by forward
substitution.
3. Then, [U]{X}={D} is used to get {X} by back
substitution.
Summary of LU Decomposition
LU Decomposition
As in Gauss elimination, LU decomposition must
employ pivoting to avoid division by zero and to
minimize round off errors. The pivoting is done
immediately after computing each column.

LU Decomposition
System of linear equations [A]{x}={B}
 a11 a12 a13   x1  b1 

a    
 21 a22 a23   x2   b2 
 a31 a32 a33   x3  b3 
Step 1: Decomposition [ L][U ]  [ A]
a11 a12 a13  1 0 0

 / 
[U ]   0 a22 a23 
/
[L ]  l 21 1 0 
 0 0 a33//   l 31 l 32 1 
LU Decomposition
Step 2: Generate an intermediate vector {D} by forward
substitution
1 0 0  d1  b1 
l     
 21 1 0 d 2   b2 
l31 l32 1 d 3  b3 
Step 3: Get {X} by back substitution.
a11 a12 a13   x1   d1 

0     
 a '22 a '23   x2   d 2 
 0 0 a ' '33   x3  d 3 
LU Decomposition-Example
 3  0.1  0.2   3 0.1 0.2 

 A  0.1 7  0.3
0.3  0.2 10 

 0 7.003 0.293
 
  0 0.19 10.02 
0.1 0.3
l 21   0.03333; l 31   0.1000
3 3
a32\ 0.19
l 32  \   0.02713
a22 7.003
3  0.1  0.2   1 0 0
 
U   0 7.003  0.293 [L ]  0.03333 1 0
0 0 10 .012   0.1000 .02713 1 
 
LU Decomposition-Example (cont’d)
Use previous L and D matrices to solve the system:
 3  0.1  0.2   x1   7.85 
 0.1 7  0 .31  x    19.3
  2   
0.3  0.2 10   x3   71.4 
Step 2: Find the intermediate vector {D} by forward substitution
 1 0 0  d1   7.85  d1   7.85 

0.0333         
 1 0 d 2    19.3  d 2    19.5617
0.1000  0.02713 1 d 3   71.4  d 3   70.0843 
LU Decomposition-Example (cont’d)
Step 3: Get {X} by back substitution.
3  0.1  0.2   x1   7.85   x1   3 

0 7.0033  0.2933  x    19.5617   x     2.5 
  2     2  
0 0 10.012   x3   70.0843   x3  7.00003
Decomposition Step
• % Decomposition Step
for k=1:n-1
[a,o]= pivot(a,o,k,n);
for i = k+1:n
a(i,k) = a(i,k)/a(k,k);
a(i,k+1:n)= a(i,k+1:n)-a(i,k)*a(k,k+1:n);
end
end
Partial Pivoting
• %Partial Pivoting
function [a,o] = pivot(a,o,k,n)
[big piv]=max(abs(a(k:n,k)));
piv=piv+(k-1);
if piv ~= k
temp = a(piv,:);
a(piv,:)= a(k,:);
a(k,:)=temp;
temp = o(piv);
o(piv)=o(k);
o(k)=temp;
end
Substitution Steps
%Forward Substitution
d(1)=bn(1);
for i=2:n
d(i)=bn(i)-a(i,1:i-1)*(d(1:i-1))';
end
• % Back Substitution
x(n)=d(n)/a(n,n);
for i=n-1:-1:1
x(i)=(d(i)-a(i,i+1:n)*x(i+1:n)')/a(i,i);
end
Matrix Inverse Using the LU
Decomposition
• LU decomposition can be used to obtain the

inverse of the original coefficient matrix [A].
• Each column j of the inverse is determined by
using a unit vector (with 1 in the jth raw ).
Matrix Inverse: LU Decomposition
[A] [A]-1 = [A]-1[A] = I
 A x1  {b}1  A x 2  {b}2  A x 3  {b}3

1  0 0
     
 A x1  0  A x 2  1   A x 3  0
0 0 1
     
1st column 2nd column 3rd column

of [A]-1 of [A]-1 of [A]-1
 A 1
 {x}1 {x}2 {x}3 
Matrix inverse using LU decomposition Example
 3  0.1  0.2   1 0 0 3  0.1  0. 2 

 A  0.1 7  0.3 [L ]  0.03333

1 0
 U   0 7.003  0.293
0.3  0.2 10   0.1000 .02713 1  0 0 10.012 
  
1A. [L]{d}1 = {b}1
 1 0 0  d1  1 d1   1 
0.0333         
 1 0 d 2   0  d 2    0.03333
0.1000  0.02713 1 d 3  0 d 3    0.1009  st
1 column
1B. Then, [U]{X}1={d}1 of [A]-1
3  0.1  0.2   x1   1   x1   0.33249 

0 7.0033  0.2933  x    0.03333   x    0.00518
  2     2  
0 0 10.012   x3    0.1009   x3   0.01008
Matrix inverse using LU decomposition
Example (cont’d)
2A. [L]{d}2 = {b}2
 1 0 0  d1  0 d1   0 
0.0333 1 0  d   1  d    1 
  2     2   
0.1000  0.02713 1 d 3  0 d 3  0.02713
2nd column
2B. Then, [U]{X}2={d}2 of [A]-1
3  0.1  0.2   x1   0   x1  0.004944

0 7.0033  0.2933  x    1    x   0.142903
  2     2  
0 0 10.012   x3  0.02713  x3   0.00271 
Example (cont’d)
3A. [L]{d}3 = {b}3
 1 0 0  d1  0 d1  0

0.0333 1 0  d   0  d   0
  2     2   
0.1000  0.02713 1 d 3  1 d 3  1
3B. Then, [U]{X}3={d}3 3rd column

of [A]-1
3  0.1  0.2   x1  0  x1  0.006798

0 7.0033  0.2933  x   0   x   0.004183
  2     2   
0 0 10.012   x3  1  x3   0.09988 
Example (cont’d)
 0.33249 0.004944 0.006798

[ A]1   0.00518 0.142903 0.004183
 0.01008 0.00271 0.09988 
Vector and Matrix Norms
Norm is a real-valued function that provides a measure of

size or “length” of vectors and matrices.
Norms are useful in studying the error behavior of
algorithms.
A simple example is a vector in three - dimensional Euclidean space
that can be represented as
 F    a b c
where a, b, and c are the distances along x, y, and z axes, repectively.
Vector and Matrix Norms (cont’d)
• The length of this vector can be simply computed as
F e  a2  b2  c2
Length or Euclidean norm of [F]
• For an n dimensional vector

 X    x1 x2  xn 
a Euclidean norm is computed as
n
Xe  i
x 2
i 1
For a matrix [A]

n n
Frobenius norm
Ae 
i 1 j 1
ai2, j
P-Norm
1/ p
 n 
p
X p   xi 
 i 1 
 1-Norm
 For Vector
n
X 1   xi
i 1
 For Matrix (Column Sum Norm)

n
A 1  max  ai , j
1 j  n i 1
• Uniform vector norm

X 
 max xi
1 i  n
• Uniform matrix norm (row sum Norm)

n
A 
 max  ai , j
1 i  n j 1
• Matrix Condition umber Defined as:
Cond  A  A  A1
• For a matrix [A], this number will be greater than or

equal to 1.
• If the coefficients of [A] are known to t-digit
precision (rounding errors~10-t) and Cond [A]=10c,
the solution [X] may be valid to only t-c digits
(rounding errors~10c-t).
Iterative Refinement
• Round-off errors can be reduced by the following procedure:
a11 x1  a12 x2  a13 x3  b1

a21 x1  a22 x2  a23 x3  b2
a31 x1  a32 x2  a33 x3  b3
•
 
Suppose an approximate solution vectors given by
~ T ~ ~ ~
X  [ x1 x2 x3 ]
• Substitute the result in the original system
~ ~ ~ ~
a11 x1  a12 x2  a13 x3  b1
~ ~ ~ ~
a21 x1  a22 x2  a23 x3 ....(Eq.1)
b2
~ ~ ~ ~
a31 x1  a32 x2  a33 x3  b3
Iterative Refinement (cont’d)
• Now, assume the exact solution is:

x1  ~
x1  x1
x ~
2 x  x
2 2
x3  ~
x3  x3
• Then:
a11 ( ~
x1  x1 )  a12 ( ~
x2  x2 )  a13 ( ~
x3  x3 )  b1
a21 ( ~
x1  x1 )  a22 ( ~
x2  x2 )  a23 ( ~
x3  …
x3 )  b2 …….(Eq.2)
a31 ( ~
x1  x1 )  a32 ( ~
x2  x2 )  a33 ( ~
x3  x3 )  b3
Iterative Refinement (cont’d)
Subtract Eq.2 from Eq.1, the result is a system of linear equations
that can be solved to obtain the correction factors  x
~
a11x1  a12 x2  a13x3  b1  b1  E1
~
a21x1  a22 x2  a23x3  b2  b2  E2
~
a31x1  a32 x2  a33x3  b3  b3  E3
The factors then can be applied to improve the solution as

specified by the equation: x1  ~
x1  x1
x ~
2 x  x
2 2
x3  ~
x3  x3
Iterative Refinement- Example (cont’d)
Solve: 4.23x1  1.06 x2  2.11x3  5.28

 2.53x1  6.77 x2  0.98 x3  5.22
1.85 x1  2.11x2  2.32 x3  2.58
The exact solution is ………  X  T  [1.00 1.00 1.00]
1- Solve the equations using [A]-1, such as {x}=[A]-1{c}
 
~ T
X  [0.991 0.997 1.00]
Iterative Refinement - Example
(cont’d)
2- Substitute the result in the original system [A]{x}={c}
~ ~ 5.24511 
~
[ A]{ X }  {C}  {C }  5.22246 
 
 2.59032
5.28  5.24511  0.0348 

~ 
E  {C}  {C}  5.22  5.22246    0.00246
 2.58  2.59032 0.0103 
Iterative Refinement - Example
(cont’d)
~
3- Solve the system of linear equations [ A]{X }  {E}
using [A]-1 to find the correction factors  x to yield
 X  T  [0.00822 0.00300  0.00000757]
x1  ~
x1  x1  0.999
x ~
2 x  x  1.00
2 2
x3  ~
x3  x3  1.00

LU Decomposition and Matrix Inversion

Uploaded by

Copyright:

Available Formats

LU Decomposition and Matrix Inversion

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

LU Decomposition and Matrix Inversion

Uploaded by

Copyright:

Available Formats

The Islamic University of Gaza

LU Decomposition and Matrix Inversion

• It becomes insufficient when solving these

• LU decomposition works on the matrix [A] and the

• LU decomposition is very useful when the vector of

As in Gauss elimination, LU decomposition must

employ pivoting to avoid division by zero and to

minimize round off errors. The pivoting is done

immediately after computing each column.

 a11 a12 a13   x1  b1 

Step 1: Decomposition [ L][U ]  [ A]

a11 a12 a13  1 0 0

a11 a12 a13   x1   d1 

Step 2: Find the intermediate vector {D} by forward substitution

 1 0 0  d1   7.85  d1   7.85 

3  0.1  0.2   x1   7.85   x1   3 

• LU decomposition can be used to obtain the

 A x1  {b}1  A x 2  {b}2  A x 3  {b}3

1st column 2nd column 3rd column

 3  0.1  0.2   1 0 0 3  0.1  0. 2 

3  0.1  0.2   x1   1   x1   0.33249 

3  0.1  0.2   x1   0   x1  0.004944

 1 0 0  d1  0 d1  0

3B. Then, [U]{X}3={d}3 3rd column

3  0.1  0.2   x1  0  x1  0.006798

 0.33249 0.004944 0.006798

Norm is a real-valued function that provides a measure of

• For an n dimensional vector

For a matrix [A]

 For Matrix (Column Sum Norm)

• Uniform vector norm

• Uniform matrix norm (row sum Norm)

• Matrix Condition umber Defined as:

• For a matrix [A], this number will be greater than or

a11 x1  a12 x2  a13 x3  b1

• Now, assume the exact solution is:

The factors then can be applied to improve the solution as

Solve: 4.23x1  1.06 x2  2.11x3  5.28

The exact solution is ………  X  T  [1.00 1.00 1.00]

1- Solve the equations using [A]-1, such as {x}=[A]-1{c}

5.28  5.24511  0.0348 

 X  T  [0.00822 0.00300  0.00000757]

You might also like