Linear Systems COS 323 Linear Systems a11x1 a12 x2 a13 x3 b1 a21x1 a22 x2 a23
Download ReportTranscript Linear Systems COS 323 Linear Systems a11x1 a12 x2 a13 x3 b1 a21x1 a22 x2 a23
Linear Systems COS 323 Linear Systems a11x1 a12 x2 a13 x3 b1 a21x1 a22 x2 a23 x3 b2 a31x1 a32 x2 a33 x3 b3 a11 a21 a 31 a12 a13 a22 a23 a32 a33 x1 b1 x2 b2 x3 b 3 Linear Systems • Solve Ax=b, where A is an nn matrix and b is an n1 column vector • Can also talk about non-square systems where A is mn, b is m1, and x is n1 – Overdetermined if m>n: “more equations than unknowns” – Underdetermined if n>m: “more unknowns than equations” Can look for best solution using least squares Singular Systems • A is singular if some row is linear combination of other rows • Singular systems can be underdetermined: 2 x1 3x2 5 4 x1 6 x2 10 or inconsistent: 2 x1 3x2 5 4 x1 6 x2 11 Inverting a Matrix • Usually not a good idea to compute x=A-1b – Inefficient – Prone to roundoff error • In fact, compute inverse using linear solver – Solve Axi=bi where bi are columns of identity, xi are columns of inverse – Many solvers can solve several R.H.S. at once Gauss-Jordan Elimination • Fundamental operations: 1. Replace one equation with linear combination of other equations 2. Interchange two equations 3. Re-label two variables • Combine to reduce to trivial system • Simplest variant only uses #1 operations, but get better stability by adding #2 (partial pivoting) or #2 and #3 (full pivoting) Gauss-Jordan Elimination • Solve: 2 x1 3x2 7 4 x1 5 x2 13 • Only care about numbers – form “tableau” or “augmented matrix”: 2 4 3 5 7 13 Gauss-Jordan Elimination • Given: 2 4 3 5 7 13 • Goal: reduce this to trivial system 1 0 0 1 ? ? and read off answer from right column Gauss-Jordan Elimination 2 4 3 5 7 13 • Basic operation 1: replace any row by linear combination with any other row • Here, replace row1 with 1/2 * row1 + 0 * row2 1 4 3 2 5 13 7 2 Gauss-Jordan Elimination 1 4 3 2 5 13 7 2 • Replace row2 with row2 – 4 * row1 1 0 3 2 1 1 7 2 • Negate row2 1 0 3 2 1 1 7 2 Gauss-Jordan Elimination 1 0 3 2 1 1 7 2 • Replace row1 with row1 – 3/2 * row2 1 0 0 1 2 1 • Read off solution: x1 = 2, x2 = 1 Gauss-Jordan Elimination • For each row i: – Multiply row i by 1/aii – For each other row j: • Add –aji times row i to row j • At the end, left part of matrix is identity, answer in right part • Can solve any number of R.H.S. simultaneously Pivoting • Consider this system: 0 2 1 3 2 8 • Immediately run into problem: algorithm wants us to divide by zero! • More subtle version: 0.001 1 3 2 2 8 Pivoting • Conclusion: small diagonal elements bad • Remedy: swap in larger element from somewhere else Partial Pivoting 0 2 1 3 2 8 • Swap rows 1 and 2: 2 0 3 1 8 2 • Now continue: 1 0 3 2 1 4 2 1 0 0 1 1 2 Full Pivoting 0 2 1 3 2 8 • Swap largest element onto diagonal by swapping rows 1 and 2 and columns 1 and 2: 3 1 2 0 8 2 • Critical: when swapping columns, must remember to swap results! Full Pivoting 3 1 1 0 2 0 2 3 23 1 0 0 1 8 2 2 3 8 * Swap results 1 and 2 3 2 1 • Full pivoting more stable, but only slightly Operation Count • For one R.H.S., how many operations? • For each of n rows: – Do n times: • For each of n+1 columns: – One add, one multiply • Total = n3+n2 multiplies, same # of adds • Asymptotic behavior: when n is large, dominated by n3 Faster Algorithms • Our goal is an algorithm that does this in 1/ n3 operations, and does not require 3 all R.H.S. to be known at beginning • Before we see that, let’s look at a few special cases that are even faster Tridiagonal Systems • Common special case: a11 a21 0 0 a12 0 0 a22 a23 0 a32 a33 a34 0 a43 a44 b1 b2 b3 b4 • Only main diagonal + 1 above and 1 below Solving Tridiagonal Systems • When solving using Gauss-Jordan: – Constant # of multiplies/adds in each row – Each row only affects 2 others a11 a21 0 0 a12 0 0 a22 a23 0 a32 a33 a34 0 a43 a44 b1 b2 b3 b4 Running Time • 2n loops, 4 multiply/adds per loop (assuming correct bookkeeping) • This running time has a fundamentally different dependence on n: linear instead of cubic – Can say that tridiagonal algorithm is O(n) while Gauss-Jordan is O(n3) Big-O Notation • Informally, O(n3) means that the dominant term for large n is cubic • More precisely, there exist a c and n0 such that running time c n3 if n > n0 • This type of asymptotic analysis is often used to characterize different algorithms Triangular Systems • Another special case: A is lower-triangular a11 a21 a31 a 41 0 0 0 a22 0 0 a32 a33 0 a42 a43 a44 b1 b2 b3 b4 Triangular Systems • Solve by forward substitution a11 a21 a31 a 41 0 0 0 a22 0 0 a32 a33 0 a42 a43 a44 b1 x1 a11 b1 b2 b3 b4 Triangular Systems • Solve by forward substitution a11 a21 a31 a 41 0 0 0 a22 0 0 a32 a33 0 a42 a43 a44 b2 a21 x1 x2 a22 b1 b2 b3 b4 Triangular Systems • Solve by forward substitution a11 a21 a31 a 41 0 0 0 a22 0 0 a32 a33 0 a42 a43 a44 b3 a31x1 a32 x2 x3 a33 b1 b2 b3 b4 Triangular Systems • If A is upper triangular, solve by backsubstitution a11 0 0 0 0 a12 a13 a14 a15 a22 a23 a24 a25 0 a33 a34 a35 0 0 a44 a45 0 0 0 a55 b5 x5 a55 b1 b2 b3 b4 b5 Triangular Systems • If A is upper triangular, solve by backsubstitution a11 0 0 0 0 a12 a13 a14 a15 a22 a23 a24 a25 0 a33 a34 a35 0 0 a44 a45 0 0 0 a55 b4 a45 x5 x4 a44 b1 b2 b3 b4 b5 Triangular Systems • Both of these special cases can be solved in O(n2) time • This motivates a factorization approach to solving arbitrary systems: – Find a way of writing A as LU, where L and U are both triangular – Ax=b LUx=b Ly=b Ux=y – Time for factoring matrix dominates computation Cholesky Decomposition • For symmetric matrices, choose U=LT • Perform decomposition a11 a12 a 13 a22 a23 a13 l11 a23 l21 a33 l31 • Ax=b LLTx=b a12 l22 l32 0 l11 0 0 l33 0 Ly=b 0 l22 l21 0 l31 l32 l33 LTx=y Cholesky Decomposition a11 a12 a 13 a13 l11 a23 l21 a33 l31 a12 a22 a23 0 l22 l32 0 l11 0 0 l33 0 l21 l22 0 l11 a11 l11 a11 2 l11l21 a12 l21 a12 l11 l11l31 a13 l31 a13 l11 l21 l22 a22 l22 a22 l21 2 2 l21l31 l22l32 a23 2 a23 l21l31 l32 l22 l31 l32 l33 Cholesky Decomposition a11 a12 a 13 a12 a22 a23 a13 l11 a23 l21 a33 l31 0 l22 l32 0 l11 0 0 l33 0 i 1 lii aii lik2 k 1 i 1 l ji aij lik l jk k 1 lii l21 l22 0 l31 l32 l33 Cholesky Decomposition • This fails if it requires taking square root of a negative number • Need another condition on A: positive definite For any v, vT A v > 0 (Equivalently, all positive eigenvalues) Cholesky Decomposition • Running time turns out to be 1/6n3 – Still cubic, but much lower constant • Result: this is preferred method for solving symmetric positive definite systems LU Decomposition • Again, factor A into LU, where L is lower triangular and U is upper triangular Ax=b LUx=b Ly=b Ux=y • Last 2 steps in O(n2) time, so total time dominated by decomposition Crout’s Method a11 a21 a 31 a12 a22 a32 a13 l11 a23 l21 a33 l31 0 l22 l32 0 u11 0 0 l33 0 • More unknowns than equations! • Let all lii=1 u12 u22 0 u13 u23 u33 Crout’s Method a11 a21 a 31 a12 a22 a32 a13 1 a23 l21 a33 l31 u11 a11 0 1 l32 0 u11 0 0 1 0 l21u11 a21 l21 a21 u11 l31u11 a31 l31 a31 u11 u12 u22 0 u12 a12 l21u12 u22 a22 u22 a22 l21u12 l31u12 l32u22 a32 a32 l31u12 l32 u22 u13 u23 u33 Crout’s Method a11 a21 a 31 a12 a22 a32 a13 1 a23 l21 a33 l31 0 1 l32 0 u11 0 0 1 0 • For i = 1..n – For j = 1..i j 1 u ji a ji l jk uki k 1 – For j = i+1..n i 1 l ji a ji l jk uki k 1 uii u12 u22 0 u13 u23 u33 Crout’s Method • Interesting note: # of outputs = # of inputs, algorithm only refers to elements not output yet – Can do this in-place! – Algorithm replaces A with matrix of l and u values, 1s are implied u11 l21 l 31 u12 u22 l32 u13 u23 u33 – Resulting matrix must be interpreted in a special way: not a regular matrix – Can rewrite forward/backsubstitution routines to use this “packed” l-u matrix LU Decomposition • Running time is 1/3n3 – Only a factor of 2 slower than symmetric case – This is the preferred general method for solving linear equations • Pivoting very important – Partial pivoting is sufficient, and widely implemented