MA/CS 375

1. MA/CS 375 Fall 2002 Lecture 21 MA/CS 375 Fall 2002

2. This Lecture • Introduce the idea of pivoting • This allows us to deal with the case where the algorithm we discussed broke down • i.e. what to do to avoid the diagonal term we need to invert became very small or even zero. MA/CS 375 Fall 2002

3. Matrices and Systems MA/CS 375 Fall 2002

4. Review Example: Upper Triangular Systems MA/CS 375 Fall 2002

5. Review Upper Triangular Systems 1) Solve the Nth equation 2) Use the value from step 1 to solve the (N-1)th equation 3) Use the values from steps 1 and 2 to solve the (N-2)th equation 4) keep going until you solve the whole system. MA/CS 375 Fall 2002

6. Review Backwards Substitution MA/CS 375 Fall 2002

7. Review LU Factorization • By now you should be persuaded that it is not difficult to solve a problem of the form Lx=b or Ux=b using forward/backward elimination • Suppose we are given a matrix A and we are able to find a lower triangular matrix L and an upper triangular matrix U such that: A = LU • Then to solve Ax=b we can solve two simple problems: • Ly = b • Ux = y • (sanity check LUx = L(Ux) = Ly =b) MA/CS 375 Fall 2002

8. Review LU Factorization (in words) • Construct a sequence of multiplier matrices (M1, M2, M3, M4,…., MN-1) such that: • (MN-1…. M4 M3 M2 M1)A is upper triangle. • A suitable candidate for M1 is the matrix ( 4 by 4 case): MA/CS 375 Fall 2002

9. Review LU Factorization Sequence cont. MA/CS 375 Fall 2002

10. Review LU Factorization cont. What does M look like: MA/CS 375 Fall 2002

11. Review LU Factorization cont. What does M-1look like: Notice: M and its inverse are very closely related. MA/CS 375 Fall 2002

12. Review LU Factorization cont. • In summary: • set L=M-1 • set U=MA • the L is upper triangular • the U is lower triangular • and A = LU, we are done. MA/CS 375 Fall 2002

13. Review Matlab Routine ForLU Factorization • Build L and U • A is modified in place • Total cost O(N3) ( we can get this by noticing the triply nested loop) MA/CS 375 Fall 2002

14. Review Matlab Routine ForLU FactorizationAnd Solution of Ax=b • Build rhs • Solve Ly=b • Solve Ux=y MA/CS 375 Fall 2002

15. Class Example • Try using your LU solver on the following matrix: • The inverse is: • What L and U do you get? MA/CS 375 Fall 2002

16. New stuff: When Basic LU Does Not Work • Remember that we are required to compute 1/A(i,i) for all the successive diagonal terms • What happens when we get to the i’th diagonal term and this turns out to be 0 or very small ?? • At any point we can swap the i’th row with one lower down which does not have a zero (or small entry) in the i’th column. • This is known as partial pivoting. MA/CS 375 Fall 2002

17. Partial Pivoting Remember before that we constructed a setof M matrices so that is upper triangle Now we construct a sequence of M matricesand P matrices so that:is upper triangle. (MN-1…. M4 M3 M2 M1)A (MN-1 PN-1 …. M4 P4 M3 P3 M2 P2 M1 P1)A MA/CS 375 Fall 2002

18. Example of a Permutation Matrix P What does P do to a matrix A under multiplication ? MA/CS 375 Fall 2002

19. Example of a Permutation Matrix P What does P do to a matrix A under multiplication ? MA/CS 375 Fall 2002

20. Example of a Permutation Matrix P In this case P has swapped thetwo bottom rows of A. MA/CS 375 Fall 2002

21. LU with Partial Pivoting MA/CS 375 Fall 2002

22. LU with Partial Pivoting • find pivot • swap rows • Do elimination MA/CS 375 Fall 2002

23. LU with Partial Pivoting • build matrix A • call our LU routine • check that: • A(piv,:) = L*U • Looks good  MA/CS 375 Fall 2002

24. Team Exercise • Make sure you all have a copy of the .m files on your laptops (courtesy of Coutsias) • Locate GEpiv • For N=100,200,400,600 • build a random NxN matrix A • time how long GEpiv takes to factorize A • end • Plot the cputime taken per test against N using loglog • On the same graph plot N^3 using loglog • Observation MA/CS 375 Fall 2002

25. Notes on LU Factorization • Big fuss over nothing, right ? • What happens if I want to solve 40000 problems with Ax=b for 40000 different b vectors?. • Cost is O(40000*N2) for the solves and O(N3) for the LU factorization. • i.e. the cost of the factorization is negligable compared to the cost of the solves !!. MA/CS 375 Fall 2002

26. From Now On • When you need to solve a square non-symmetric system of equations (Ax=b) with lots of different b vectors use the Matlab built in routine lu: • [L,U,P] = lu(A); • y = L\(Pb); • x = U\y; • For a one time solve use: x=A\b MA/CS 375 Fall 2002

27. Interlude on Norms • We all know that the “size” of 2 is the same as the size of -2, namely: • ||2|| = sqrt(2*2) = 2 • ||-2|| = sqrt((-2)*(-2)) = 2 • The ||..|| notation is known as the norm function. MA/CS 375 Fall 2002

28. Norm of A Vector • We can generalize the scalar norm to a vector norm, however now there are more choices for a vector : MA/CS 375 Fall 2002

29. Norm of A Vector • Matlab command: • norm(x,1) • norm(x,2) • norm(x,inf) MA/CS 375 Fall 2002

30. Norm of A Matrix • We can generalize the vector norm to a matrix norm, however now there are more choices for a matrix : MA/CS 375 Fall 2002

31. Norm of A Matrix • Matlab command: • norm(A,1) • norm(A,2) • norm(A,inf) • norm(A,fro) MA/CS 375 Fall 2002

32. Matrix Norm in Action • Theorem 5: (page 185 van Loan) • If is the stored version of then where and i.e. an error of order ||A||1eps occurs when a real matrix is stored in floating point. MA/CS 375 Fall 2002

33. Next Lecture • We will see matrix norms in action… • Those of you who own a digital camera, please bring it in. • Otherwise bring a cd with some picture files on (jpeg, gif, tiff,…) – I have the previously taken pictures too.. • We will start talking about interpolation MA/CS 375 Fall 2002