[PPT] - LU Factorization Marco Chiarandini Department of Mathematics & PowerPoint Presentation

SLIDE 1

DM559 Linear and Integer Programming

LU Factorization

Marco Chiarandini

Department of Mathematics & Computer Science University of Southern Denmark [Based on slides by Lieven Vandenberghe, UCLA]

SLIDE 2

Operation Count LU Factorization Iterative Methods

Outline

1. Operation Count
2. LU Factorization
3. Iterative Methods

2

SLIDE 3

Operation Count LU Factorization Iterative Methods

Outline

1. Operation Count
2. LU Factorization
3. Iterative Methods

3

SLIDE 4

Operation Count LU Factorization Iterative Methods

Complexity of matrix algorithms

flop counts
vector-vector operations
matrix-vector product
matrix-matrix product

4

SLIDE 5

Operation Count LU Factorization Iterative Methods

Flop counts

floating-point operation (flop)

one floating-point addition, subtraction, multiplication, or division
other common definition: one multiplication followed by one addition

flop counts of matrix algorithm

total number of flops is typically a polynomial of the problem dimensions
usually simplified by ignoring lower-order terms

applications

a simple, machine-independent measure of algorithm complexity
not an accurate predictor of computation time on modern computers

5

SLIDE 6

Operation Count LU Factorization Iterative Methods

Vector-vector operations

inner product of two n-vectors

xTy = x1y1 + x2y2 + . . . + xnyn n multiplications and n − 1 additions = 2n flops (2n if n ≫ 1)

addition or subtraction of n-vectors: n flops
scalar multiplication of n-vector : n flops

6

SLIDE 7

Operation Count LU Factorization Iterative Methods

Matrix-vector product

matrix-vector product with m × n-matrix A: y = Ax m elements in y; each element requires an inner product of length n: (2n − 1)m flops approximately 2mn for large n special cases

m = n, A diagonal: n flops
m = n, A lower triangular: n(n + 1) flops
A very sparse (lots of zero coefficients): #flops ≪ 2mn

7

SLIDE 8

Operation Count LU Factorization Iterative Methods

Matrix-matrix product

product of m × n-matrix A and n × p-matrix B: C = AB mp elements in C; each element requires an inner product of length n: mp(2n − 1)flops approximately 2mnp for large n.

8

SLIDE 9

Operation Count LU Factorization Iterative Methods

Outline

1. Operation Count
2. LU Factorization
3. Iterative Methods

9

SLIDE 10

Operation Count LU Factorization Iterative Methods

Overview

factor-solve method
LU factorization
solving Ax = b with A nonsingular
the inverse of a nonsingular matrix
LU factorization algorithm
effect of rounding error
sparse LU factorization

10

SLIDE 11

Operation Count LU Factorization Iterative Methods

Definitions

Definition (Triangular Matrices) An n × n matrix is said to be upper triangular if aij = 0 for i > j and lower triangular if aij = 0 for i < j. Also A is said to be triangular if it is either upper triangular or lower triangular. Example:   3 0 0 2 1 0 1 4 3     3 5 1 0 1 3 0 0 7   Definition (Diagonal Matrices) An n × n matrix is diagonal if aij = 0 whenever i = j. Example:   1 0 0 0 1 0 0 0 3  

11

SLIDE 12

Operation Count LU Factorization Iterative Methods

Linear systems in practice

How much work in Gaussian elimination?

We eliminate x1, x2, . . . xn in this order and then find the values of

xn, xn−1, . . . , x1 by back substitution.

Elimination process: when the variables x1, x2, . . . , xk−1 have been

elminated, we are confronted with a system of n − k + 1 equations in the remaining n − k + 1 variables.

To eliminate xk amounts to n − k + 1 divisions, followed by

(n − k)(n − k + 1) multiplications and (n − k)(n − k + 1) additions.

In total for the whole method:

Additions: n

k=1(n − k)(n − k + 1) + n k=1(n − k) = n(n − 1)(2n + 5)/6

Multiplications: n

k=1(n − k)(n − k + 1) + n k=1(n − k) = n(n − 1)(2n + 5)/6

Divisions: n

k=1(n − k + 1) = n(n + 1)/2

Hence: about n3/3 multiplications that dominate

12

SLIDE 13

Operation Count LU Factorization Iterative Methods

Linear systems in practice

Gaussian elimination is not suitable for large-scale systems because of:

computer roundoff errors
memory usage
speed

Computer methods are based on factorizations or iterative methods. Factorization methods need much less computations in practice for sparse matrices (2n3/3 operations in the general case for LU decompositions) Moreover, factorization methods lend themselves well to solve sequences of linear systems.

13

SLIDE 14

Factor-solve approach

to solve Ax = b, first write A as a product of ‘simple’ matrices A = A1A2 · · · Ak then solve (A1A2 · · · Ak)x = b by solving k equations A1z1 = b, A2z2 = z1, . . . , Ak−1zk−1 = zk−2, Akx = zk−1 examples

Cholesky factorization (for positive definite A)

k = 2, A = LLT

sparse Cholesky factorization (for sparse positive definite A)

k = 4, A = PLLTP

LU factorization 7-2

SLIDE 15

Complexity of factor-solve method

#flops = f + s

f is cost of factoring A as A = A1A2 · · · Ak (factorization step)
s is cost of solving the k equations for z1, z2, . . . zk−1, x (solve step)
usually f ≫ s

example: positive definite equations using the Cholesky factorization f = (1/3)n3, s = 2n2

LU factorization 7-3

SLIDE 16

Multiple right-hand sides

two equations with the same matrix but different right-hand sides Ax = b, A˜ x = ˜ b

factor A once (f flops)
solve with right-hand side b (s flops)
solve with right-hand side ˜

b (s flops) cost: f + 2s instead of 2(f + s) if we solve second equation from scratch conclusion: if f ≫ s, we can solve the two equations at the cost of one

LU factorization 7-4

SLIDE 17

LU factorization

LU factorization without pivoting A = LU

L unit lower triangular, U upper triangular
does not always exist (even if A is nonsingular)

LU factorization (with row pivoting) A = PLU

P permutation matrix, L unit lower triangular, U upper triangular
exists if and only if A is nonsingular (see later)

cost: (2/3)n3 if A has order n

LU factorization 7-5

SLIDE 18

Solving linear equations by LU factorization

solve Ax = b with A nonsingular of order n factor-solve method using LU factorization

1. factor A as A = PLU ((2/3)n3 flops)
2. solve (PLU)x = b in three steps
permutation: z1 = P Tb (0 flops)
forward substitution: solve Lz2 = z1 (n2 flops)
back substitution: solve Ux = z2 (n2 flops)

total cost: (2/3)n3 + 2n2 flops, or roughly (2/3)n3 this is the standard method for solving Ax = b

LU factorization 7-6

SLIDE 19

Multiple right-hand sides

two equations with the same matrix A (nonsingular and n × n): Ax = b, A˜ x = ˜ b

factor A once
forward/back substitution to get x
forward/back substitution to get ˜

x cost: (2/3)n3 + 4n2 or roughly (2/3)n3 exercise: propose an efficient method for solving Ax = b, AT ˜ x = ˜ b

LU factorization 7-7

SLIDE 20

Inverse of a nonsingular matrix

suppose A is nonsingular of order n, with LU factorization A = PLU

inverse from LU factorization

A−1 = (PLU)−1 = U −1L−1P T

gives interpretation of solve step: evaluate

x = A−1b = U −1L−1P Tb in three steps z1 = P Tb, z2 = L−1z1, x = U −1z2

LU factorization 7-8

SLIDE 21

Computing the inverse

solve AX = I by solving n equations AX1 = e1, AX2 = e2, . . . , AXn = en Xi is the ith column of X; ei is the ith unit vector of size n

one LU factorization of A: 2n3/3 flops
n solve steps: 2n3 flops

total: (8/3)n3 flops conclusion: do not solve Ax = b by multiplying A−1 with b

LU factorization 7-9

SLIDE 22

LU factorization without pivoting

partition A, L, U as block matrices: A =

a11

A12 A21 A22

,

L =

1

L21 L22

,

U =

u11

U12 U22

a11 and u11 are scalars
L22 unit lower-triangular, U22 upper triangular of order n − 1

determine L and U from A = LU, i.e.,

a11

A12 A21 A22

=
1

L21 L22 u11 U12 U22

=
u11

U12 u11L21 L21U12 + L22U22

LU factorization

7-10

SLIDE 23

recursive algorithm:

determine first row of U and first column of L

u11 = a11, U12 = A12, L21 = (1/a11)A21

factor the (n − 1) × (n − 1)-matrix A22 − L21U12 as

A22 − L21U12 = L22U22 this is an LU factorization (without pivoting) of order n − 1 cost: (2/3)n3 flops (no proof)

LU factorization 7-11

SLIDE 24

Example

LU factorization (without pivoting) of A =   8 2 9 4 9 4 6 7 9   write as A = LU with L unit lower triangular, U upper triangular A =   8 2 9 4 9 4 6 7 9   =   1 l21 1 l31 l32 1     u11 u12 u13 u22 u23 u33  

LU factorization 7-12

SLIDE 25

first row of U, first column of L:

  8 2 9 4 9 4 6 7 9   =   1 1/2 1 3/4 l32 1     8 2 9 u22 u23 u33  

second row of U, second column of L:
9

4 7 9

−
1/2

3/4 2 9

=
1

l32 1 u22 u23 u33

8

−1/2 11/2 9/4

=
1

11/16 1 8 −1/2 u33

third row of U: u33 = 9/4 + 11/32 = 83/32

conclusion: A =   8 2 9 4 9 4 6 7 9   =   1 1/2 1 3/4 11/16 1     8 2 9 8 −1/2 83/32  

LU factorization 7-13

SLIDE 26

Not every nonsingular A can be factored as A = LU

A =   1 2 1 −1   =   1 l21 1 l31 l32 1     u11 u12 u13 u22 u23 u33  

first row of U, first column of L:

  1 2 1 −1   =   1 1 l32 1     1 u22 u23 u33  

second row of U, second column of L:
2

1 −1

=
1

l32 1 u22 u23 u33

u22 = 0, u23 = 2, l32 · 0 = 1 ?

LU factorization 7-14

SLIDE 27

LU factorization (with row pivoting)

if A is n × n and nonsingular, then it can be factored as A = PLU P is a permutation matrix, L is unit lower triangular, U is upper triangular

not unique; there may be several possible choices for P, L, U
interpretation: permute the rows of A and factor P TA as P TA = LU
also known as Gaussian elimination with partial pivoting (GEPP)
cost: (2/3)n3 flops

we will skip the details of calculating P, L, U

LU factorization 7-15

SLIDE 28

Example

  5 5 2 9 6 8 8   =   1 1 1     1 1/3 1 15/19 1     6 8 8 19/3 −8/3 135/19   the factorization is not unique; the same matrix can be factored as   5 5 2 9 6 8 8   =   1 1 1     1 1 3 −19/5 1     2 9 5 5 27  

LU factorization 7-16

SLIDE 29

Effect of rounding error

10−5

1 1 1 x1 x2

=
1
exact solution:

x1 = −1 1 − 10−5, x2 = 1 1 − 10−5 let us solve the equations using LU factorization, rounding intermediate results to 4 significant decimal digits we will do this for the two possible permutation matrices: P =

1

1

r

P =

1

1

LU factorization

7-17

SLIDE 30

first choice of P: P = I (no pivoting)

10−5

1 1 1

=
1

105 1 10−5 1 1 − 105

L, U rounded to 4 decimal significant digits

L =

1

105 1

,

U =

10−5

1 −105

forward substitution
1

105 1 z1 z2

=
1
=

⇒ z1 = 1, z2 = −105 back substitution

10−5

1 −105 x1 x2

=
1

−105

=

⇒ x1 = 0, x2 = 1 error in x1 is 100%

LU factorization 7-18

SLIDE 31

second choice of P: interchange rows

1

1 10−5 1

=
1

10−5 1 1 1 1 − 10−5

L, U rounded to 4 decimal significant digits

L =

1

10−5 1

,

U =

1

1 1

forward substitution
1

10−5 1 z1 z2

=
1
=

⇒ z1 = 0, z2 = 1 backward substitution

1

1 1 x1 x2

=
1
=

⇒ x1 = −1, x2 = 1 error in x1, x2 is about 10−5

LU factorization 7-19

SLIDE 32

conclusion:

for some choices of P, small rounding errors in the algorithm cause very

large errors in the solution

this is called numerical instability: for the first choice of P, the

algorithm is unstable; for the second choice of P, it is stable

from numerical analysis: there is a simple rule for selecting a good

(stable) permutation (we’ll skip the details, since we skipped the details

f the factorization algorithm)
in the example, the second permutation is better because it permutes

the largest element (in absolute value) of the first column of A to the 1,1-position

LU factorization 7-20

SLIDE 33

Sparse linear equations

if A is sparse, it is usually factored as A = P1LUP2 P1 and P2 are permutation matrices

interpretation: permute rows and columns of A and factor ˜

A = P T

1 AP T 2

˜ A = LU

choice of P1 and P2 greatly affects the sparsity of L and U: many

heuristic methods exist for selecting good permutations

in practice: #flops ≪ (2/3)n3; exact value is a complicated function of

n, number of nonzero elements, sparsity pattern

LU factorization 7-21

SLIDE 34

Conclusion

different levels of understanding how linear equation solvers work: highest level: x = A\b costs (2/3)n3; more efficient than x = inv(A)*b intermediate level: factorization step A = PLU followed by solve step lowest level: details of factorization A = PLU

for most applications, level 1 is sufficient
in some situations (e.g., multiple right-hand sides) level 2 is useful
level 3 is important only for experts who write numerical libraries

LU factorization 7-22

SLIDE 35

Operation Count LU Factorization Iterative Methods

Outline

1. Operation Count
2. LU Factorization
3. Iterative Methods

14

SLIDE 36

Operation Count LU Factorization Iterative Methods

Numerical Solutions

A matrix A is said to be ill conditioned if relatively small changes in the

entries of A can cause rlatively large changes in the solutions of Ax = b.

A is said to be well conditioned if relatively small changes in the entries
f A result in relatively small changes in the solutions of Ax = b.
reaching RREF as in Gauss-Jordan requires more computation and more

numerical instability hence disadvantageous.

Gauss elimination is a direct method: the amount of operations can be

specified in advance. Indirect or Iterative methods work by iteratively improving approximate solutions until a desired accuracy is reached. Amount of operations depend on the accuracy required. (way to go if the matrix is sparse)

15

SLIDE 37

Operation Count LU Factorization Iterative Methods

Gauss-Seidel Iterative Method

Example

+ 25.00 = 64.062

17