Computational Methods CMSC AMSC MAPL 460 LU Decomposition Sensitivity Ramani Duraiswami Dept of Computer Science Solving Linear Systems One idea compute inverse Not usually a good idea unless inverse is computable easily and accurately using some matrix property Leads to increased errors and is more expensive usually Easy systems to solve Diagonal system Triangular system On board and then matlab Cost of diagonal solve is O n x zeros n 1 for k 1 n x k b k A k k end Solving a triangular system Cost of solving a triangular system Loop of size n Each loop has a cost of k or n k So total cost is n 1 n 2 n n n2 Gaussian Elimination Zero elements of first column below 1st row multiplying 1st row by 0 3 and add to 2nd row multiplying 1st row by 0 5 and add to 3rd row Results in Zero elements of first column below 2nd row Swap rows Multiply 2nd row by 0 04 and add to 3rd Representing linear systems as matrix vector equations Represent it as a matrix vector equation linear system We will apply the familiar elimination technique and then see its matrix equivalent Solution Start from last equation which can be solved by division Next substitute in the previous line and continue This describes the way to do the algorithm by hand How to represent it using matrices Also how do we solve another system that has the same matrix Upper triangular matrix we end up with will be the same but the sequence of operations on the r h s needs to be repeated Gaussian Elimination LU Matrix decomposition It turns out that Gaussian elimination corresponds to a particular matrix decomposition Product of permutation lower triangular and upper triangular matrices What is a permutation matrix It rearranges a system of equations and changes the order Multiplying by it swaps the order of rows in a matrix Essentially a rearrangement of the identity Nice property transpose is its inverse PPT I LU Decomposition What is an upper triangular matrix Elements below diagonal are zero Lower triangular matrix Elements above diagonal are zero Unit lower triangular matrix Elements along diagonal are one Upper triangular part of Gauss Elimination is clear final matrix we end up with What about lower triangular and permutation LU PA Identify the elements of L and P L has the multipliers we used in the elimination steps P has a record of the row swaps we did to avoid dividing by small numbers In fact we can write each step of Gaussian elimination in matrix form Solving a system with the LU decomposition Ax b LU PA PT LUx b L Ux Pb Solve Ly Pb Then Ux y Solving a system with the LU decomposition Ax b LU PA PT LUx b L Ux Pb Solve Ly Pb Then Ux y Look at LU code Initialize Matrix size Permutation vector Second output argument to max is index of max element If max element is zero then we need not eliminate Exchange rows update permutation vector Look at LU code Multipliers for each row below diagonal Note multipliers are stored in the lower triangular part of A Vectorized update A i k A k j multiplies column vector by row vector to produce a square rank 1 matrix of order n k matrix is then subtracted from the submatrix of the same size in the bottom right corner of A In a programming language without vector and matrix operations this update of a portion of A would be done with doubly nested loops on i and j Cost is n2 and done n times for a total cost of n3 Computes decomposition in the matrix A itself Here they are separated but when memory is important it can be left there Code to solve linear system using LU In Matlab the backslash operator can be used to solve linear systems For square matrices it employs LU or special variants Lower triangular Upper triangular symmetric Symmetric LU is called Cholesky decomposition A LLT Upper and lower triangular are equal transposes If matrix not positivedefinite go to regular solution Code continues Call LU Solve y Lb Solve x Uy LU Wrap up Is pivoting necessary in LU Consider Exact solution is Let 0 5 Solution without pivoting gives Is pivoting necessary With pivoting Elimination gives With answers Close to exact Another example from the book Another example from the book Correct answer is 0 1 1 T How accurate are answers from LU We solve the equation Ax b Let true solution be x Let obtained solution be x Then error is e x x Error is not computable also called Forward error New concept residual also called Backward error Residual is the difference between the original right hand side and the right hand side obtained with the obtained solution r b Ax Guarantee LU produces answers with small residuals on computers with IEEE floating point Do small residuals mean small errors Return to our example Compute residual We have exactly solved a nearby problem Another example Solution has small residual but very large error In fact signs of the solution are opposite Why Can condition numbers tell us what is going on Condition numbers Recall definition of condition number Condition Number of a Matrix A measure of how close a matrix is to singular cond A A A A 1 i maximum stretch max i maximum shrink min i i cond I 1 cond singular matrix Norm can be any norm One norm is easy to compute Relation between condition number and error In words relative error is smaller than norm of residual divided by norm of rhs times condition number So larger condition number means larger error Properties of the condition number Closing remarks Never compute matrix inverse Use a stable algorithm Check residual and condition number of problem If condition number is large do not trust solution Can problem be reformulated somehow