MA/CS 375

1. Sit in your groups MA/CS 375 Fall 2002 Lecture 18 MA/CS 375 Fall 2002

2. Projects? • Demo from each group MA/CS 375 Fall 2002

3. Matrices and Systems MA/CS 375 Fall 2002

4. Recall • Given a matrix M then if it exists the inverse of a matrix M-1 is that matrix which satisfies: MA/CS 375 Fall 2002

5. Examples • If what is • If what is • If A is an NxN matrix how can we calculate its inverse ?  MA/CS 375 Fall 2002

6. More Examples • If what is • If what is • More generally when does an inverse exist MA/CS 375 Fall 2002

7. Example (Problems Calculating an Inverse) • Previously we saw how Matlab can get wrong answers when finite precision effects come into play (monsters) • The same type of behavior can be seen when Matlab calculates the inverse of a matrix. • Let’s consider the simple matrix: MA/CS 375 Fall 2002

8. Example cont • The exact inverse of A can be calculated as: • Now let’s see how well Matlab does: MA/CS 375 Fall 2002

9. Example cont • The inverse of A can be calculated as: • Now let’s see how well this exact solution works in Matlab: MA/CS 375 Fall 2002

10. Example cont With this formulation of the product of A and its inverse only satisfiesthe definition to 6 decimalplaces for delta=0.001 MA/CS 375 Fall 2002

11. Example cont Now let us compare theanswer Matlab gives ususing inv(A): MA/CS 375 Fall 2002

12. Example cont Now let us compare theanswer Matlab gives ususing inv(A): Looks pretty good,doesn’t it! MA/CS 375 Fall 2002

13. Example cont But not so quick. Now try A*inv(A) Now Matlab is only goodto 6 decimal places. MA/CS 375 Fall 2002

14. Example comments • So what it is it about thismatrix which is causing this odd behaviour?. • For small delta it is nearly singular. For those with some background the determinant of A is: • i.e. the inverse of A nearly does not exist • Matlab is actually solving for a matrix near to A MA/CS 375 Fall 2002

15. Example team exercise • Test the accuracy of Matlab’s inv(A) • Team 1 use delta=1e-2 • Team 2 use delta = 1e-4 • Team 3 use delta = 1e-6 • Team 4 use delta= 1e-9 • Team 5 use delta = 1e-11 • Try both A*inv(A) and inv(A)*A MA/CS 375 Fall 2002

16. Solving Triangular Systems • So let’s go back to basics. • There are some systems of equations which are “easy” to solve. • One example is a system with a triangular matrix structure. MA/CS 375 Fall 2002

17. Lower Triangular Systems equations matrix form MA/CS 375 Fall 2002

18. Lower Triangular Systems Volunteer to solve this? i.e. find x,y so that these two equations are satisfied equations MA/CS 375 Fall 2002

19. Lower Triangular Systems • note that we can easily figure out x • then knowing x we can easily figure out y. • i.e. 2x = 1  x = ½ • 3x+2y = 2  y = ½(2-3x) = ½(2-3/2) = 1/4 equations MA/CS 375 Fall 2002

20. Larger Problem • Volunteers?. • on the board • using Matlab MA/CS 375 Fall 2002

21. General Triangular System(matrix form) MA/CS 375 Fall 2002

22. General Triangular System We can write down an algorithm to solve this problem using the “Forward substitution algorithm”: MA/CS 375 Fall 2002

23. Forward Substitution Algorithm 1) Solve the 1st equation 2) Use the value from step 1 to solve the 2nd equation 3) Use the values from steps 1 and 2 to solve the 3rd equation 4) keep going until you solve the whole system. MA/CS 375 Fall 2002

24. Solving Triangular Systems in Matlab • Given a lower triangular matrix L • Given a vector b • Solve Lx = b for the unknown vector x MA/CS 375 Fall 2002

25. Solving Triangular Systems in MatlabIn action: • Given a lower triangular matrix L • Given a vector b • Solve Lx = b for the unknown vector x MA/CS 375 Fall 2002

26. Solving Triangular Systems in MatlabObservations • Given a lower triangular matrix L • Given a vector b • Solve Lx = b for the unknown vector x • The algorithm scales as O(N2) • We relied on the diagonal terms of L to be non-zero • What happens if there is a zero diagonal term ? MA/CS 375 Fall 2002

27. Class Exercise • Write the lower triangular system solver using only one loop • i.e. vectorize the inner loop • Test your algorithm against thebuilt in Matlab matrix divisionroutine. • You have 10 minutes from when I say go – hard limit (no hand in required) replace MA/CS 375 Fall 2002

28. Upper Triangular Systems MA/CS 375 Fall 2002

29. Upper Triangular Systems Volunteer: tell us how to solve this? (hint: this is not difficult) MA/CS 375 Fall 2002

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

31. 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

32. Matlab Code for Backward Substitution MA/CS 375 Fall 2002

33. Testing Backwards Substitution MA/CS 375 Fall 2002

34. Finite Precision Effects • Notice that x(1) depends on the values of x(2),x(3),…,x(N) • Remembering that round off rears its ugly head when we deal with (small+large) numbers • Volunteer to design a U matrix and b vector, which do not give good answers when used in this algorithm. • (goal is to give you some intuition for cases which will break an algorithm) MA/CS 375 Fall 2002

35. Summary of Lecture 18 • We have seen how finite precision can effect the calculation of an inverse matrix, by example. • We have discussed algorithms to solve Lx=b and Ux=b • We have seen how to vectorize part of the forward/backward substitution algorithms • We tried to break the forward/backward substitution algorithms MA/CS 375 Fall 2002

36. Next Lecture • Ok – so we know how to solve L or U matrices • We will see how to write fairly general matrices as a product of a lower and upper matrix. E.g. A = LU • We will investigate LU factorization, including sensitivity to the structure of A • Introduce the idea of pivoting MA/CS 375 Fall 2002