Notes Block Approach to LU Assignment 1 is out (due October 5) - PDF document

Notes Block Approach to LU � Assignment 1 is out (due October 5) � Rather than get bogged down in details of GE (hard to see forest for trees) � Matrix storage: usually column-major � Partition the equation A=LU � Gives natural formulas for algorithms � Extends to block algorithms cs542g-term1-2006 1 cs542g-term1-2006 2 Cholesky Factorization Pivoting � LU and Modified Cholesky can fail � If A is symmetric positive definite, can cut • Example: if A 11 =0 work in half: A=LL T � Go back to Gaussian Elimination ideas: reorder • L is lower triangular the equations (rows) to get a nonzero entry � In fact, nearly zero entries still a problem • Possibly due to cancellation error => few significant � If A is symmetric but indefinite, possibly digits still have the Modified Cholesky • Dividing through will taint rest of calculation factorization: A=LDL T � Pivoting strategy: reorder to get the biggest entry • L is unit lower triangular on the diagonal • Partial pivoting: just reorder rows • D is diagonal • Complete pivoting: reorder rows and columns (expensive) cs542g-term1-2006 3 cs542g-term1-2006 4 Pivoting in LU Symmetric Pivoting � Can express it as a factorization: � Problem: partial (or complete) pivoting destroys symmetry A=PLU � How can we factor a symmetric indefinite matrix • P is a permutation matrix: just the identity with reliably but twice as fast as unsymmetric its rows (or columns) permuted matrices? • Store the permutation, not P! � One idea: symmetric pivoting PAP T =LDL T • Swap the rows the same as the columns � But let D have 2x2 as well as 1x1 blocks on the diagonal • Partial pivoting: Bunch-Kaufman (LAPACK) • Complete pivoting: Bunch-Parlett (safer) cs542g-term1-2006 5 cs542g-term1-2006 6

Reconsidering RBF Gibbs � RBF interpolation has advantages: � Globally smooth • Mesh-free calculation also • Optimal in some sense makes for • Exponential convergence (each point extra overshoot/ data point improves fit everywhere) undershoot • Defined everywhere (Gibbs phenomena) � But some disadvantages: around • It � s a global calculation discontinuities (even with compactly supported functions) � Can � t easily control • Big dense matrix to form and solve effect (though later we � ll revisit that… cs542g-term1-2006 7 cs542g-term1-2006 8 Noise Linear Least Squares � If data contains noise (errors), RBF strictly � Idea: instead of interpolating data + noise, interpolates them approximate � If the errors aren � t spatially correlated, lots � Pick our approximation from a space of of discontinuities: RBF interpolant functions we expect (e.g. not wiggly -- becomes wiggly maybe low degree polynomials) to filter out the noise � Standard way of defining it: k � f ( x ) = � i � i ( x ) i = 1 ( ) n � 2 � = argmin f j � f ( x j ) � j = 1 cs542g-term1-2006 9 cs542g-term1-2006 10 Rewriting Normal Equations � Write it in matrix-vector form: � First attempt at finding minimum: set the gradient equal to zero 2 � � n k � � f i � � j � j ( x i ) = b � Ax 2 2 (called “the normal equations”) � � � � � 2 = 0 i = 1 j = 1 � x b � Ax 2 ( ) b = f 1 T � f 2 f n � ( ) = 0 � x ( b � Ax ) T ( b � Ax ) ( ) T x = � 1 � k � � A ij = � i ( x j ) (a rectangular n � k matrix) ( ) = 0 � x b T b � 2 x T A T b + x T A T Ax � 2 A T b + 2 A T Ax = 0 A T Ax = A T b cs542g-term1-2006 11 cs542g-term1-2006 12

Good Normal Equations Bad Normal Equations � A T A is a square k � k matrix � What if k=n? (k probably much smaller than n) At least for 2-norm condition number, k(A T A)=k(A) 2 � Symmetric positive (semi-)definite • Accuracy could be a problem… � In general, can we avoid squaring the errors? cs542g-term1-2006 13 cs542g-term1-2006 14