Linear Systems II CS3220 Summer 2008 Jonathan Kaldor
Revisiting the LU Factorization • Goal: Solve Ax = b • Represents n linear equations, with n unknowns • Assume A is full rank (invertible, nonsingular) • Idea: Simplify without changing the solution • What form is simpler to solve?
Triangular Matrices
Triangular Matrices • Upper Triangular Matrix a 11 a 12 .... a 1n x 1 b 1 0 a 22 .... a 2n x 2 b 2 . . . . . . = . . . . . . . . a n-1,n-1 a n-1,n x n-1 b n-1 0 0 .... a nn x n b n
Triangular Matrices • Upper Triangular Matrix a 11 a 12 .... a 1n x 1 b 1 0 a 22 .... a 2n x 2 b 2 . . . . . . = . . . . . . . . a n-1,n-1 a n-1,n x n-1 b n-1 0 0 .... a nn x n b n Solve for x n
Triangular Matrices • Upper Triangular Matrix a 11 a 12 .... a 1n x 1 b 1 0 a 22 .... a 2n x 2 b 2 Substitute x n , . . . . . . = solve for x n-1 . . . . . . . . a n-1,n-1 a n-1,n x n-1 b n-1 0 0 .... a nn x n b n Solve for x n
Triangular Matrices • Upper Triangular Matrix Repeat for each variable a 11 a 12 .... a 1n x 1 b 1 0 a 22 .... a 2n x 2 b 2 Substitute x n , . . . . . . = solve for x n-1 . . . . . . . . a n-1,n-1 a n-1,n x n-1 b n-1 0 0 .... a nn x n b n Solve for x n
Triangular Matrices • Upper Triangular Matrix Repeat for each variable a 11 a 12 .... a 1n x 1 b 1 0 a 22 .... a 2n x 2 b 2 Substitute x n , . . . . . . = solve for x n-1 . . . . . . . . a n-1,n-1 a n-1,n x n-1 b n-1 0 0 .... a nn x n b n Solve for x n • “Backward Substitution”
LU Factorization • Convert Ax = b into MA x= Mb , where MA is upper triangular • M = M n-1 ... M 2 M 1 • M 1 - matrix that eliminates first variable from equations 2...n, M 2 is... • Property: If C and D are both upper (or both lower) triangular, CD is also upper (or lower) triangular (so M is lower tri)
LU Factorization • So after n-1 steps, have: M n-1 ... M 2 M 1 Ax = M n-1 ... M 2 M 1 b • Define U = M n-1 ... M 2 M 1 A . • Upper triangular by construction • Multiply by L = M 1-1 M 2-1 ... M n-1-1 • Then LU = A
Gaussian Elimination • What do the update matrices M i look like? 1 0 ... 0 0 ... 0 0 1 ... 0 0 ... 0 Ones along ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ the diagonal 0 0 ... 1 0 ... 0 0 0 ... -a i+1,i /a i,i 1 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 0 0 ... -a n,i /a i,i 0 ... 1 Entries in the i-th column, below the diagonal, chosen to cancel i-th variable in equations i+1:n
Gaussian Elimination • What do the inverse matrices M i-1 look like?
Gaussian Elimination • What do the inverse matrices M i-1 look like? 1 0 ... 0 0 ... 0 0 1 ... 0 0 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ Note: entries below diagonal 0 0 ... 1 0 ... 0 have flipped signs 0 0 ... a i+1,i /a i,i 1 ... 0 ⋮ ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ 0 0 ... a n,i /a i,i 0 ... 1
The L in LU • So each M i-1 is lower triangular as well • In fact, unit lower triangular (ones on diagonal) • How about the product M 1-1 M 2-1 ... M n-1-1 ?
Solving with LU factorization • Have LUx = b • Replace Ux with y • Solve Ly = b (forward substitution) • Then solve Ux = y (backward substitution)
What does L look like?
What does L look like? • We know: L is lower triangular, consists of product of M i-1 matrices
What does L look like? • We know: L is lower triangular, consists of product of M i-1 matrices • Do we need to multiply M i-1 together to form L ?
What does L look like? • We know: L is lower triangular, consists of product of M i-1 matrices • Do we need to multiply M i-1 together to form L ? • No!
What does L look like? • We know: L is lower triangular, consists of product of M i-1 matrices • Do we need to multiply M i-1 together to form L ? • No! • Do we need to scale/update rows like U ?
What does L look like? • We know: L is lower triangular, consists of product of M i-1 matrices • Do we need to multiply M i-1 together to form L ? • No! • Do we need to scale/update rows like U ? • No! (chalkboard)
A Vectorized U update • How can we quickly apply M i to partially computed U matrix?
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 -1 -2 0 0 1 -1
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 1 2 1 -1 -2 0 0 0 1 0 1 -1 0 1 -1 Scale first equation by -(-1/1), add to second equation
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 0 0 1 0 1 -1 Scale first equation by -(-1/1), add to second equation
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 0 0 1 0 1 -1
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 0 0 1 0 1 -1 Scale second equation by -(2/0), add to first equation
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 ? 0 0 1 0 1 -1 Scale second equation by -(2/0), add to first equation
Pivoting • Gaussian (and Gauss-Jordan) elimination both have failure cases • Consider the following matrix: 1 2 1 ? 0 0 1 Division by zero! 0 1 -1 Scale second equation by -(2/0), add to first equation
Pivoting 1 2 1 0 0 1 0 1 -1
Pivoting • Does that mean that the system is singular? 1 2 1 0 0 1 0 1 -1
Pivoting • Does that mean that the system is singular? • Does the order we specify equations matter? 1 2 1 0 0 1 0 1 -1
Pivoting • Observe: if we swap the second and third equations, Gaussian elimination works as expected 1 2 1 1 2 1 0 0 1 0 1 -1 0 1 -1 0 0 -1 • We can reorder equations as-needed in order to get “good” pivots (pivot: A jj )
What makes a good pivot? • Again, back to our example. Suppose A(2,2) = 1e-15 • Instead of dividing by zero, dividing by a very small number • Almost as bad • Observation: Want to swap equations so that diagonal element is as large as possible in absolute value
Partial Pivoting • At i-th step • Look at elements in i-th column below the diagonal , find largest in absolute value (say, j-th) • Swap i-th and j-th rows. Largest entry now in A(i,i) • Advantage: all multipliers in L ≤ 1 (recall, multipliers are A(j,i)/A(i,i) )
Permutation Matrices • Want to represent these row swaps using matrix notation • Suppose we want to swap rows i and j in A • Take identity matrix • Swap rows i and j, call this matrix P • Then PA is the same as A with rows i and j interchanged
Permutation Matrices • Note: P -1 = P T • PA : swaps rows • AP : swaps columns • Like M i matrices, never form them explicitly • Oftentimes, don’t even need to explicitly swap rows
Pivoting in Matrix Formulation • Introduce P i matrices during computation M n-1 P n-1 ... M 2 P 2 M 1 P 1 Ax = M n-1 P n-1 ... M 2 P 2 M 1 P 1 b • Job of each P i : move row with largest element in i-th column to the i-th row
Pivoting in Matrix Formulation • U = M n-1 P n-1 ... M 2 P 2 M 1 P 1 A • L is still formed from M i , but multipliers below diagonal get permutated • P = P n-1 ... P 2 P 1 • End up with LU = PA
An Example 1 2 2 4 4 2 4 6 4
LU with Partial Pivoting • Factor A into PA = LU • To solve Ax = b : • Note: A = P T LU • So we have P T LUx = b • Multiply by P to get LUx = Pb
LU with Partial Pivoting • Now have LUx = Pb • Solve Ly = Pb for the vector y using forward substitution • Solve Ux = y for the vector x using backward substitution
LU with Partial Pivoting • What happens if all the elements in a column below the diagonal are zero?
Solving in MATLAB • MATLAB provides an operator to solve linear systems • A \ b solves equation Ax = b • Can also specify matrix B • Be careful: also / operator (solves AB -1 ) • Solves more than just n x n systems (more to come!)
Solving in MATLAB • To compute LU factorization of matrix A , use function lu(A) • Syntax: [L, U] = lu(A) • L is not lower triangular... ‘psychologically lower triangular’ • Or: [L,U,P] = lu(A) • P: permutation matrix
Wait! ‘Partial’ Pivoting?! • We reordered equations to get a good pivot element... could we also relabel our unknowns to get a better element? • Short answer: Yes • Longer answer: Yes, but typically partial pivoting is good enough
LU Factorization: A Summary • Generalized solution method for n x n systems of linear equations in n unknowns • Requires pivoting in order to avoid certain failure cases
Recommend
More recommend
