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

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 ]

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

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

Operation Count LU Factorization Complexity of matrix algorithms Iterative Methods • flop counts • vector-vector operations • matrix-vector product • matrix-matrix product 4

Operation Count LU Factorization Flop counts Iterative Methods 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

Operation Count LU Factorization Vector-vector operations Iterative Methods • inner product of two n-vectors x T y = x 1 y 1 + x 2 y 2 + . . . + x n y n n multiplications and n − 1 additions = 2 n flops (2 n if n ≫ 1) • addition or subtraction of n -vectors: n flops • scalar multiplication of n -vector : n flops 6

Operation Count LU Factorization Matrix-vector product Iterative Methods matrix-vector product with m × n -matrix A : y = A x m elements in y ; each element requires an inner product of length n : ( 2 n − 1 ) m flops approximately 2 mn 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 ≪ 2 mn 7

Operation Count LU Factorization Matrix-matrix product Iterative Methods 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 ( 2 n − 1 ) flops approximately 2 mnp for large n . 8

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

Operation Count LU Factorization Overview Iterative Methods • 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

Operation Count LU Factorization Definitions Iterative Methods Definition (Triangular Matrices) An n × n matrix is said to be upper triangular if a ij = 0 for i > j and lower triangular if a ij = 0 for i < j . Also A is said to be triangular if it is either upper triangular or lower triangular. Example:  3 0 0   3 5 1  2 1 0 0 1 3     1 4 3 0 0 7 Definition (Diagonal Matrices) An n × n matrix is diagonal if a ij = 0 whenever i � = j . Example:   1 0 0 0 1 0   0 0 3 11

Operation Count LU Factorization Linear systems in practice Iterative Methods How much work in Gaussian elimination? • We eliminate x 1 , x 2 , . . . x n in this order and then find the values of x n , x n − 1 , . . . , x 1 by back substitution. • Elimination process: when the variables x 1 , x 2 , . . . , x k − 1 have been elminated, we are confronted with a system of n − k + 1 equations in the remaining n − k + 1 variables. • To eliminate x k 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: � n k = 1 ( n − k )( n − k + 1 ) + � n Additions: k = 1 ( n − k ) = n ( n − 1 )( 2 n + 5 ) / 6 Multiplications: � n k = 1 ( n − k )( n − k + 1 ) + � n k = 1 ( n − k ) = n ( n − 1 )( 2 n + 5 ) / 6 � n Divisions: k = 1 ( n − k + 1 ) = n ( n + 1 ) / 2 Hence: about n 3 / 3 multiplications that dominate 12

Operation Count LU Factorization Linear systems in practice Iterative Methods 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 (2 n 3 / 3 operations in the general case for LU decompositions) Moreover, factorization methods lend themselves well to solve sequences of linear systems. 13

Factor-solve approach to solve Ax = b , first write A as a product of ‘simple’ matrices A = A 1 A 2 · · · A k then solve ( A 1 A 2 · · · A k ) x = b by solving k equations A 1 z 1 = b , A 2 z 2 = z 1 , . . . , A k − 1 z k − 1 = z k − 2 , A k x = z k − 1 examples • Cholesky factorization (for positive definite A ) A = LL T k = 2 , • sparse Cholesky factorization (for sparse positive definite A ) A = PLL T P k = 4 , LU factorization 7-2

Complexity of factor-solve method # flops = f + s • f is cost of factoring A as A = A 1 A 2 · · · A k (factorization step) • s is cost of solving the k equations for z 1 , z 2 , . . . z k − 1 , x (solve step) • usually f ≫ s example: positive definite equations using the Cholesky factorization f = (1 / 3) n 3 , s = 2 n 2 LU factorization 7-3

Multiple right-hand sides two equations with the same matrix but different right-hand sides x = ˜ Ax = b, A ˜ b • factor A once ( f flops) • solve with right-hand side b ( s flops) • solve with right-hand side ˜ b ( s flops) cost: f + 2 s 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

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) n 3 if A has order n LU factorization 7-5

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) n 3 flops) 2. solve ( PLU ) x = b in three steps • permutation: z 1 = P T b (0 flops) • forward substitution: solve Lz 2 = z 1 ( n 2 flops) • back substitution: solve Ux = z 2 ( n 2 flops) total cost : (2 / 3) n 3 + 2 n 2 flops, or roughly (2 / 3) n 3 this is the standard method for solving Ax = b LU factorization 7-6

Multiple right-hand sides two equations with the same matrix A (nonsingular and n × n ): x = ˜ Ax = b, A ˜ b • factor A once • forward/back substitution to get x • forward/back substitution to get ˜ x cost: (2 / 3) n 3 + 4 n 2 or roughly (2 / 3) n 3 exercise : propose an efficient method for solving A T ˜ x = ˜ Ax = b, b LU factorization 7-7

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 − 1 L − 1 P T • gives interpretation of solve step: evaluate x = A − 1 b = U − 1 L − 1 P T b in three steps z 1 = P T b, z 2 = L − 1 z 1 , x = U − 1 z 2 LU factorization 7-8

Computing the inverse solve AX = I by solving n equations AX 1 = e 1 , AX 2 = e 2 , . . . , AX n = e n X i is the i th column of X ; e i is the i th unit vector of size n • one LU factorization of A : 2 n 3 / 3 flops • n solve steps: 2 n 3 flops total : (8 / 3) n 3 flops conclusion : do not solve Ax = b by multiplying A − 1 with b LU factorization 7-9

LU factorization without pivoting partition A , L , U as block matrices: � � � � � � a 11 A 12 1 0 u 11 U 12 A = , L = , U = A 21 A 22 L 21 L 22 0 U 22 • a 11 and u 11 are scalars • L 22 unit lower-triangular, U 22 upper triangular of order n − 1 determine L and U from A = LU , i.e. , � � � � � � a 11 A 12 1 0 u 11 U 12 = A 21 A 22 L 21 L 22 0 U 22 � � u 11 U 12 = u 11 L 21 L 21 U 12 + L 22 U 22 LU factorization 7-10

recursive algorithm : • determine first row of U and first column of L u 11 = a 11 , U 12 = A 12 , L 21 = (1 /a 11 ) A 21 • factor the ( n − 1) × ( n − 1) -matrix A 22 − L 21 U 12 as A 22 − L 21 U 12 = L 22 U 22 this is an LU factorization (without pivoting) of order n − 1 cost : (2 / 3) n 3 flops (no proof) LU factorization 7-11

Example LU factorization (without pivoting) of  8 2 9  A = 4 9 4   6 7 9 write as A = LU with L unit lower triangular, U upper triangular       8 2 9 1 0 0 u 11 u 12 u 13  = A = 4 9 4 l 21 1 0 0 u 22 u 23      6 7 9 l 31 l 32 1 0 0 u 33 LU factorization 7-12