Transcript ENGINEERING 212 98A

Chapter 10
LU Decomposition
L and U Matrices
Lower Triangular Matrix
 l 11
l
 L    21
l 31

 l 41
0
l 22
l 32
l 42
0

0
0
l 33
0

l 34 l 44 
0
 u11 u12 u13 u13 
 0 u

u
u
22
23
24 
U  
 0
0 u 33 u 34 


0
0 u44 
 0
Upper Triangular Matrix
LU Decomposition
 Another method for solving matrix equations
Idea behind the LU Decomposition - start with
 A  x   b
We know (because we did it in Gauss Elimination)
we can write
 u11

U  x   d   



u12
u13
u22
u23
u33
u14   x1  d 1 





u24   x 2  d 2 
  
u34   x 3  d 3 
   
u44   x 4  d 4 
LU Decomposition
Assume there exists [L]
 Such that
 l 11
l
 L    21
l 31

 l 41
l 22
l 32
l 33
l 42
l 34





l 44 
L  U  x  d    A x  b
This implies
LU   A
&
L d   b
The Steps of LU Decomposition
 L U    A
 LU  x  B
LU Decomposition
 LU Decomposition
* Based on Gauss elimination
*More efficient
 Decomposition Methods (not unique)
* Doolittle decomposition lii = 1
* Crout decomposition
uii = 1 (omitted)
* Cholesky decomposition (for symmetric
matrices) uii = lii
LU Decomposition
Three Basic Steps
(1) Factor (decompose) [A] into [L] and [U]
(2) given {b}, determine {d} from [L]{d} = {b}
(3) using [U]{x} = {d} and back-substitution,
solve for {x}
Advantage: Once we have [L] and [U], we can
use many different {b}’s without repeating the
decomposition process
LU Decomposition
LU decomposition / factorization
[A]{x}=[L][U]{x}={b}
Forward substitution
[L]{d}={b}
Back substitution
[U]{x}={d}
Forward substitutions are more efficient
than elimination
Simple Truss
F45
4
F14
5
F35
F24
F25
1 
H1
V1

F12


2
W
F23
3
V2
Exampe: Forces in a Simple Truss
1
0

0

0
0

0
0

0
0

0
0
1
0
0
0
0
0
1
0
0
0 1
sin 
cos
0
0
0
0
0
0
sin 
cos 
sin 
cos
0
0
0
0 1
0
 sin 
0  sin 
0
0
0
0 1
0 0
0 0
0
0
0
0 0
0 0
0 0
0  cos
0
0
0
0
0
0
0
1
0
0
0
cos 
0
0  sin 
0
cos
0   V1   0 
 
0   H 1   0 
0   V 3  100
  

0   F12   0 
sin 
0   F14   0 
   

 cos
0   F23   0 
0
0   F24   0 
  

0
1  F25   0 
sin 
0   F35   0 
  

cos  1  F45   0 
0
0
0
0
[A] depends on geometry only; but {b} varies with applied load
LU Factorization
Gauss elimination [A]{x}={b}
elimination steps need to be repeated for
each different {b}
LU factorization / decomposition
[A]=[L][U]
does not involve the RHS {b} !!!
Doolittle LU Decomposition
Doolittle Algorithm
Get 1’s on diagonal of [L] (lii =1)
Operate on rows and columns sequentially,
narrowing down to single element
Identical to LU decomposition based on
Gauss elimination (see pp238-239), but
different approach
Doolittle LU Decomposition
 a 11
a
A   21
a 31

 a 41
 u11
l u
A   21 11
l 31 u11

 l 41 u11
a 12
a 22
a 13
a 23
a 32
a 33
a 42
a 43
a 14   1 0
0
a 24   l 21
1 0


 l 31 l 32
a 34
1
 
a 44   l 41 l 42 l 43
0   u11 u12 u13 u14 
0   0 u 22 u 23 u 24 
0  0
0 u 33 u 34 


1  0
0
0 u44 
u12
u13
l 21 u12  u22
l 21 u13  u23
l 31 u12  l 32 u 22
l 31 u13  l 32 u 23  u 33
l 41 u12  l 42 u 22
l 41 u13  l 42 u 23  l 43 u 33


l 21 u14  u 24


l 31 u14  l 32 u 24  u 34

l 41 u14  l 42 u 24  l 43 u 34  u44 
u14
1st row : u11  a11 ; u12  a12 ; u13  a13 ; u14  a14
1st column: l21  a21 / u11; l31  a31 / u11; l41  a41 / u11
Doolittle LU Decomposition
 u11
l u
A   21 11
l 31 u11

 l 41 u11
u12
u13
l 21 u12  u22
l 21 u13  u23
l 31 u12  l 32 u 22
l 31 u13  l 32 u 23  u 33
l 41 u12  l 42 u 22
l 41 u13  l 42 u 23  l 43 u 33


l 21 u14  u 24


l 31 u14  l 32 u 24  u 34

l 41 u14  l 42 u 24  l 43 u 34  u44 
u14
2nd row : l21u12  u22  a22 ; l21u13  u23  a23 ; l21u14  u24  a24
u22  a22  l21u12

u23  a23  l21u13
u  a  l u
 24
24
21 14
2nd column: l31u12  l32u22  a32 ; l41u12  l42u22  a42
l32  (a32  l31u12 ) / u22

l42  (a42  l41u12 ) / u22
Doolittle LU Decomposition
 u11
l u
A   21 11
l 31 u11

 l 41 u11
u12
u13
l 21 u12  u22
l 21 u13  u23
l 31 u12  l 32 u 22
l 31 u13  l 32 u 23  u 33
l 41 u12  l 42 u 22
l 41 u13  l 42 u 23  l 43 u 33


l 21 u14  u 24


l 31 u14  l 32 u 24  u 34

l 41 u14  l 42 u 24  l 43 u 34  u44 
u14
3rd row : l31u13  l32u23  u33  a33 ; l31u14  l32u24  u34  a34
u33  a33  l31u13  l32u23

u34  a34  l31u14  l32u24
3rdcolumn: l43  (a43  l41u13  l42u23 )/ u33
4th row: u44  a44  l41u14  l42u24  l43u34
Algorithm of Doolittle Decomposition
[ A]  [ L][U ]
(N by N matrix)
Staring the first row of U  ,
u1,i  a1,i , for i  1, 2,......, N ;
then the first column of  L  , l j ,1  a j ,1 / u1,1 , for j  2,......, N ;
Then alternatively determine the 2nd row of U  ,
u2,i  a2,i  l2,1u1,i for i  2,3,......, N ;
and 2nd column of  L  ;
l j ,2  (a j ,2  l j ,1u1,2 ) / u2,2 ,
and n
th
for j  3,......, N ; then
n 1
row of U  , un,i  an,i   ln, k uk ,i for i  n,..., N ;
k 1
and n
th
n 1


column of  L  , l j , n   a j , n   l j , k uk , n  / un, n ,


k 1
for j  n  1,..., N ; ..........until N th row of [U ].
Forward Substitution
Once [L] is formed, we can use forward substitution
instead of forward elimination for different {b}’s
 l11
l
L d    21
 l 31

 l 41
d 1
d
 2

d 3
d 4
0
l 22
l 32
l 42
0   d 1   l11 d 1


0 0  d 2   l 21 d 1  l 22 d 2

 

l 33
0  d 3   l 31 d 1  l 32 d 2  l 33 d 3

  
l 43 l 44  d 4   l 41 d 1  l 42 d 2  l 43 d 3  l 44 d 4 
0
 b1 / l11
 ( b2  l 21 d 1 ) / l 22
 ( b3  l 31 d 1  l 32 d 2 ) / l 33
 b1 
b 
 2
 
b3 
b4 
Very efficient for
large matrices !
 ( b4  l 41 d 1  l 42 d 2  l 43 d 3 ) / l 44
Back Substitution
 u11
0
U  x  
0

0
 x4
x
 3

 x2
 x1
u14   x 1 
u24   x 2 
 
u34   x 3 

u44   x 4 
 u11 x 1  u12 x 2  u13 x 3  u14 x 4 
u x  u x  u x

 22 2

23 3
24 4


 u33 x 3  u34 x 4

 u44 x 4

 ( d 3  u34 x 4 ) / u33
Identical to
Gauss
elimination
u12
u13
u22
u23
0
u33
0
0
 d 4 / u44
 ( d 2  u23 x 3  u24 x 4 ) / u22
 ( d 1  u12 x 2  u13 x 3  u14 x 4 ) / u11
d1 
d 
 2
 
d 3 
d 4 
Forward Substitution
Example:
0 0
 1
 1
1 0

L d  
 0 1/2 1

1 14
 6
d 1
d
 2

d 3
d 4
1
0 d1 



d
0  2
  
0  d 3 
 
1 d 4 
 1
  1
 
    b
 2
 1
 1
 0
 1  d 1  1  1  0


 d   
 2  (1/2)d2  2
 2
  33
 1  6d1  d 2  14d3  1  6  14(2) 33
Back-Substitution
1
0
U x  
0

0
 x4

 x3

 x2
 x 1
3  x 1 
0   x 2 
 
0 1
4  x3 
 
0
0  70  x 4 
0
2
2
4
 33/  70  33/70
 4 x 4  2  4/35
 2 x 3  8/35
 1  2 x 3  3 x 4  13/70
 1
 0




 2
 33
  13/70
 8/35

{ x}  
  4/35


 33/70
Forward and Back Substitutions
Forward-substitution
i l
d i  bi   lij d j
for i  1,2, , n
j 1
Error in textbook
(p. 165)
Back-substitution (identical to Gauss elimination)
x n  d n / ann
di 
xi 
n
u x
ij
j i1
uii
j
for i  n  1, n  2,  , 3, 2,1
Forward and Back Substitutions
Example: Forward and Back Substitutions
»
»
»
L
A=[1 0 2 3; -1 2 2 -3; 0 1 1 4; 6 2 2 4];
b=[1 -1 2 1]';
[L,U] = LU_factor(A); (without
=
1.0000
0
0
-1.0000
1.0000
0
0
0.5000
1.0000
6.0000
1.0000
14.0000
U =
1
0
2
3
0
2
4
0
0
0
-1
4
0
0
0
-70
» x=LU_solve(L,U,b)
x =
-0.1857
0.2286
-0.1143
0.4714
pivoting)
0
0
0
1.0000
Forward and back
substitution
MATLAB
M-file
LU_factor
Pivoting in LU Decomposition
Still need pivoting in LU decomposition
Messes up order of [L]
What to do?
Need to pivot both [L] and a permutation
matrix [P]
Initialize [P] as identity matrix and pivot
when [A] is pivoted  Also pivot [L]
LU Decomposition with Pivoting
Permutation matrix [ P ]
- permutation of identity matrix [ I ]
Permutation matrix performs “bookkeeping”
associated with the row exchanges
Permuted matrix [ P ] [ A ]
LU factorization of the permuted matrix
[P][A]=[L][U]
Solution
[ L ] [ U ] {x} = [ P ] {b}
Permutation Matrix
Bookkeeping for row exchanges
Example: [ P1] interchanges row 1 and 3
0
0

1

0
0 1 0   a 11

1 0 0  a 21
0 0 0  a 31

0 0 1 a 41
a 12
a 22
a 13
a 23
a 32
a 33
a 42
a 43
a 14  a 31
a 24  a 21

a 34   a 11
 
a 44  a 41
a 32
a 22
a 33
a 23
a 12
a 13
a 42
a 43
a 32
a 22
a 33
a 23
a 42
a 43
a 12
a 13
a 34 
a 24 
a 14 

a 44 
Multiple permutations [ P ]
0
0

0

1
0 1 0   a 11

1 0 0  a 21
0 0 1 a 31

0 0 0  a 41
a 12
a 22
a 13
a 23
a 32
a 33
a 42
a 43
a 14  a 31
a 24  a 21

a 34  a 41
 
a 44   a 11
a 34 
a 24 
a 44 

a 14 
LU Decomposition with Pivoting
x1
 x1
6 x1
 2 x2
x2
 2 x2
 2 x3
 2 x3
 x3
 2 x3
 3 x4
 3 x4
 4 x4
 4 x4
1
 1
2
1
 1
 1
A b  
 0

 6
0
2
1
2
2
3
1
2  3  1
1
4
2

2
4
1
Start with
 1
 1
A  
 0

 6
0
2
1
2
2
3
2  3
1
4

2
4
1

 1


L  

1 


1


1

 1


P   

1 


1


No need to consider {b} in decomposition
Forward Elimination
 6
 1
A  
 0

 1
2 2
4
2 2  3

1 1
4

0 2
3
f 21  1/6
f 31  0
Interchange rows 1 & 4
f 41  1/6
Gauss elimination of first column
2
2
4
6
0
7/3 7/3  7/3

A  
0
1
1
4


0

1/3
5/3
7/3


1


  1/6 1


L  

0
1 


1/6
1


Save fi1 in the first column of [L]
0
0
P   
0

1
0 0
1 0
0 1
0 0
1
0

0

0
Forward Elimination
2
2
4
6
0
7/3 7/3  7/3

A  
0
1
1
4


0

1/3
5/3
7/3


f 32  3/7
No interchange required
f 42  1/7
Gauss elimination of second column
2
2
4
6
0 7/3 7/3  7/3

A  
0
0
0
5


0
0
2
2


1


  1/6

1

L  

0
3/7 1 


1/6

1/7
1


Save fi2 in the second column of [L]
0
0
P   
0

1
0
1
0
0
0
0
1
0
1
0

0

0
Forward Elimination
2
2
4
6
0 7/3 7/3  7/3

A  
0
0
0
5


0
2
2
0
Interchange rows 3 &4
Partial pivoting for [L] and [P]
2
2
4
6
0 7/3 7/3  7/3

U   
0
0
2
2


0
0
0
5


1


  1/6

1

L  
 1/6  1/7 1 


0
3/7
1


0
0
P   
1

0
0 0 1
1 0 0 
0 0 0

0 1 0
Forward Elimination
2
2
4
6
 0 7/3 7/3  7/3

A  
0
0
2
2


0
0
0
5


f 43  0
Gauss elimination of third column
2
2
4
6
 0 7/3 7/3  7/3

U   
0
0
2
2


0
0
0
5


1


  1/6

1

L  
 1/6  1/7 1 


0
3/7
0
1


Save fi3 in third column of [L]
0
0
P   
1

0
0 0 1
1 0 0

0 0 0

0 1 0
LU Decomposition with Pivoting
 1
 1
A  
 0

 6
0
2
1
2
2
2
1
2
3
3
;
4

4
0
0
P   
1

0
0
1
0
0
0
0
0
1
1
0

0

0
Gauss elimination with partial pivoting
0
 1
  1/6
1

L 
 1/6  1/7

3/7
 0
0
0
1
0
0
2
2
4
6
 1
0 7/3 7/3  7/3
 1
0 
 


U  
P b   
0
0
0
2
2
 1



 2 
1
0
0
5
0
LU Decomposition with Pivoting
0
 1
  1/6
1

LU  
 1/6  1/7

3/7
 0
0
0
1
0
0
0 
0

1
2
2
4  6
6
0 7/3 7/3  7/3   1


0
0
2
2  1

 
0
0
0
5

  0
1
0

  1/6
1

Forward L  d  
 1/6  1/7
substitution

0
3/7

0 0
0 0

1 0

0 1
 d1 
d 
 2
 
d 3 
d 4 
2
2
4   x1 
6
0 7/3 7/3  7/3  x 

  2  




U
x

Back
0
0
2
2  x3 

 
substitution
0
0
0
5

  x4 
2
2
0
1
2
4
2  3
 P A

2
3

1
4
 d1 
 1
d 
  1
 
 2

 
 
 1
d 3 
 2 
d 4 
1
 x1 

x 
  5/6


 2



 
 5/7 
 x3 
 33/14
 x 4 
1

  5/6




5/7


 33/14
  13/70
 8/35





4/35


 33/70
LU Decomposition with Pivoting
partial pivoting
LU Decomposition with Pivoting
» A=[1 0 2 3; -1 2 2 -3; 0 1 1 4; 6 2 2 4];
» b=[1 -1 2 1]';
» [L,U,P]=LU_pivot(A);
L =
1.0000
0
0
-0.1667
1.0000
0
0.1667
-0.1429
1.0000
0
0.4286
0
U =
6.0000
2.0000
2.0000
0
2.3333
2.3333
0
0
2.0000
0
0
0
T1 =
6
2
2
4
-1
2
2
-3
1
0
2
3
0
1
1
4
T2 =
6
2
2
4
-1
2
2
-3
1
0
2
3
0
1
1
4
0
0
0
1.0000
4.0000
-2.3333
2.0000
5.0000
T1 = [L][U]
T2 = [P][A]
Verify
T1= T2
LU Decomposition with Pivoting
» A=[1 0 2 3; -1 2 2 -3; 0 1 1 4; 6 2 2 4];
» b=[1; -1; 2; 1];
» [L,U,P]=LU_pivot(A);
» Pb=P*b
Pb =
1
-1
1
2
» x=LU_Solve(L,U,Pb)
d =
1.0000
-0.8333
0
0
d =
1.0000
-0.8333
0.7143
0
d =
1.0000
-0.8333
0.7143
2.3571
LU Decomposition
x =
0
0
-0.1143
0.4714
x =
0
0.2286
-0.1143
0.4714
[L][U]{x} = [P]{b}
Forward Substitution
x =
-0.1857
0.2286
-0.1143
0.4714
x =
-0.1857
0.2286
-0.1143
0.4714
Back Substitution
MATLAB’s Methods
LU factorization [L,U] = lu (A)
-- returns [L] and [U]
[L, U, P] = lu (A)
-- returns also the permutation matrix [P]
Cholesky factorization: R = chol(A)
-- [A] = [R][R]
Determinant: det (A)
MATLAB function lu
» A=[1 0 2 3; -1 2 2 -3; 0 1 1 4; 6 2 2 4];
» b=[1; -1; 2; 1];
» [L,U]=lu(A)
L =
0.1667
-0.1429
-0.1667
1.0000
0
0.4286
1.0000
0
U =
6.0000
2.0000
0
2.3333
0
0
0
0
» L*U
ans =
1
0
2
-1
2
2
0
1
1
6
2
2
1.0000
0
0
0
0
0
1.0000
0
2.0000
2.3333
2.0000
0
4.0000
-2.3333
2.0000
5.0000
3
-3
4
4
LU Decomposition
without Pivoting
» d=L\b
d =
1.0000
-0.8333
0.7143
2.3571
» x=U\d
x =
-0.1857
0.2286
-0.1143
0.4714
Forward and
Back
Substitutions
Cholesky LU Factorization
If [A] is symmetric and positive definite, it is
convenient to use Cholesky decomposition.
[A] = [L][L]T= [U]T[U]
Cholesky factorization can be used for either
symmetric or non-symmetric matrices
No pivoting or scaling needed if [A] is symmetric
and positive definite (all eigenvalues are positive)
If [A] is not positive definite, the procedure may
encounter the square root of a negative number
Cholesky LU Factorization
For a general non-symmetric matrix
 a11
a
A   21
 a31

 a41
a12
a13
a22
a32
a23
a33
a42
a43
a14   d 11
a24   l 21

a34   l 31
 
a44   l 41
0
d 22
l 32
l 42
0   d 11 u12 u13 u14 
0
0   0 d 22 u23 u24 


d 33
0  0
0 d 33 u34 


l 43 d 44   0
0
0 d 44 
0
For symmetric and positive definite matrices
 a11
a
A   12
 a13

 a14
a12
a13
a22
a23
a23
a33
a24
a34
a14   u11
a24   u12

a34   u13
 
a44   u14
0
u22
u23
u24
0   u11 u12 u13 u14 
0
0   0 u22 u23 u24 


u33
0  0
0 u33 u34 


u34 u44   0
0
0 u44 
0
Cholesky LU Decomposition
2
 u11

A  u11 u12
u11 u13

 u11 u14
u11 u12
u11 u13
2
2
u12
 u22
u13 u12  u23 u22
u13 u12  u23 u22
u14 u12  u24 u22
2
2
2
u13
 u23
 u33
u14 u13  u24 u23  u34 u33


u14 u12  u24 u22

u14 u13  u24 u23  u34 u33 

2
2
2
2
u14
 u24
 u34
 u44

u11 u14
1st column / row: u11  a11 ; u12  a12 /u11; u13  a13/u11; u14  a14 /u11
2
2
2 nd column / row : u12
 u22
 a22 ; u13 u12  u23 u22  a23 ; u14 u12  u24 u22  a24
2
u22  a22  u12
; u23  ( a23  u13 u12 ) / u22 ; u24  ( a24  u14 u12 ) / u22
2
2
2
3 rd column / row : u13
 u23
 u33
 a33 ; u13 u14  u23 u24  u33 u34  a34
2
2
u33  a33  u13
 u23
; u34  ( a34  u14 u13  u24 u23 ) / u33
2
2
2
2
4 th row : u14
 u24
 u34
 u44
 a44
2
2
2
 u44  a44  u14
 u24
 u34
Example: Cholesky LU
2
 9  6 12  3  u11

 6
5 9
2   u11 u12

A  
 12  9 21
0   u11 u13

 

3
2
0
6

  u11 u14
u11 u12
u11 u13
2
2
u12
 u22
u13 u12  u23 u22
u13 u12  u23 u22
u14 u12  u24 u22
2
2
2
u13
 u23
 u33
u14 u13  u24 u23  u34 u33


u14 u12  u24 u22

u14 u13  u24 u23  u34 u33 

2
2
2
2
u14
 u24
 u34
 u44

u11 u14
1st column/row: u11  9  3; u12  6 / 3  2; u13  12/ 3  4 ; u14  3 / 3  1
2
2
2 nd column/row: u12
 u22
 5; u13 u12  u23 u22  9; u14 u12  u24 u22  2
u22  5  ( 2 ) 2  1; u23  ( 9  4 ( 2)) / 1  1; u24  ( 2  ( 1)( 2)) / 1  0
2
2
2
3 rd column/row: u13
 u23
 u33
 21; u13 u14  u23 u24  u33 u34  0
u33  21  ( 4 ) 2  ( 1) 2  2; u34  (0  ( 1)(4 )  (0 )(1)) / 2  2
2
2
2
2
4 th row : u14
 u24
 u34
 u44
6
 u44  6  ( 1) 2  (0 ) 2  ( 2 ) 2  1
Cholesky LU Factorization
[A] = [U]T[U] = [U] [U]
Recurrence relations
i 1
uii  aii   uki2
k 1
i 1
uij 
aij   uki ukj
k 1
uii
for j  i  1, , n
Script file for Cholesky Decomposition
Symmetric Matrix L = U
Compute only the upper triangular elements
» A=[9 -6
A =
9
-6
12
-3
» [L,U] =
L =
3
-2
4
-1
U =
3
0
0
0
12 -3; -6 5 -9 2; 12 -9 21 0; -3 2 0 6]
-6
12
-3
5
-9
2
-9
21
0
2
0
6
Cholesky(A)
0
1
-1
0
-2
1
0
0
0
0
2
2
4
-1
2
0
0
0
0
1
-1
0
2
1
 9  6 12  3
 6

5

9
2

A  
 12  9 21
0



3
2
0
6


Symmetric
[L] = [U]
MATLAB Function: chol
»
»
»
U
A=[9 -6 12 -3; -6 5 -9 2; 12 -9 21 0; -3 2 0 6];
b=[24; -19; 51; 11];
U = chol(A)
=
3
-2
4
-1
0
1
-1
0
0
0
2
2
0
0
0
1
» U'*U
ans =
9
-6
12
-3
-6
5
-9
2
12
-9
21
0
-3
2
0
6
» d = U'\b
Forward substitution
d =
8
-3
8
3
Back substitution
» x = U\d
x =
1
-2
1
3