Effect of Perturbations on Singularity in Linear Algebra - Prof. Luke Olson, Study notes of Computer Science

Download Effect of Perturbations on Singularity in Linear Algebra - Prof. Luke Olson and more Study notes Computer Science in PDF only on Docsity! Lecture 9 Conditioning, Banded Systems L. Olson Department of Computer Science University of Illinois at Urbana-Champaign September 22, 2006 L. Olson (UIUC) CS 257 September 22, 2006 1 / 33 The Algorithms Listing 1: (forward) GE with SPP 1 Initialize ` = [1, . . . , n] 2 Set s to be the max of rows 3 for k = 1 to n 4 rmax = 0 5 for i = k to n 6 r = |a`ik/s`i | 7 if(r > rmax) 8 rmax = r 9 j = i 10 end 11 swap `j and `k 12 for i = k + 1 to n 13 xmult = a`ik/a`kk 14 a`ik = xmult 15 for j = k + 1 to n 16 a`ij = a`ij − xmult · a`k j 17 end 18 end 19 end See Gauss algorithm on page 285 of NMC. L. Olson (UIUC) CS 257 September 22, 2006 3 / 33 Geometric Interpretation of Singularity The equations for two parallel and coincident lines are[ 2 1 2 1 ] [ x1 x2 ] = [ 6 6 ] The equations for two nearly parallel lines are[ 2 1 2 + δ 1 ] [ x1 x2 ] = [ 6 6 + δ ] L. Olson (UIUC) CS 257 September 22, 2006 7 / 33 Geometric Interpretation of Singularity 0 1 2 3 4 0 2 4 6 8 A and b are consistent A is nonsingular 0 1 2 3 4 0 2 4 6 8 A and b are inconsistent A is singular 0 1 2 3 4 0 2 4 6 8 A and b are consistent A is singular 0 1 2 3 4 0 2 4 6 8 A and b are consistent A is ill conditioned L. Olson (UIUC) CS 257 September 22, 2006 8 / 33 Effect of Perturbations to b Consider the solution of a 2× 2 system where b = [ 1 2/3 ] One expects that the exact solutions to Ax = [ 1 2/3 ] and Ax = [ 1 0.6667 ] will be different. Should these solutions be a lot different or a little different? L. Olson (UIUC) CS 257 September 22, 2006 9 / 33 Effect of Perturbations to b (Proof, p. 2) Take norm of equation (2) ‖δxb‖ = ‖A−1 δb‖ Applying consistency requirement of matrix norms ‖δx‖ 6 ‖A−1‖‖δb‖ (3) Similarly, Ax = b gives ‖b‖ = ‖Ax‖, and ‖b‖ 6 ‖A‖‖x‖ (4) Rearrangement of equation (4) yields 1 ‖x‖ 6 ‖A‖ ‖b‖ (5) L. Olson (UIUC) CS 257 September 22, 2006 12 / 33 Effect of Perturbations to b (Proof) Multiply Equation (4) by Equation (3) to get ‖δxb‖ ‖x‖ 6 ‖A‖‖A−1‖‖δb‖ ‖b‖ (6) Summary: If x + δxb is the exact solution to the perturbed system A(x + δxb) = b + δb then ‖δxb‖ ‖x‖ 6 ‖A‖‖A−1‖‖δb‖ ‖b‖ L. Olson (UIUC) CS 257 September 22, 2006 13 / 33 Effect of Perturbations to A Perturb A with δA such that ‖δA‖ ‖A‖  1, The perturbed system is (A + δA)(x + δxA) = b Analysis shows that ‖δxA‖ ‖x + δxA‖ 6 ‖A‖‖A−1‖‖δA‖ ‖A‖ Thus, the effect of the perturbation is small only if ‖A‖‖A−1‖ is small. ‖δxA‖ ‖x + δxA‖  1 only if ‖A‖‖A−1‖ ∼ 1 L. Olson (UIUC) CS 257 September 22, 2006 14 / 33 Computational Stability In Practice, applying Gaussian elimination with partial pivoting and back substitution to Ax = b gives the exact solution, x̂, to the nearby problem (A + E)x̂ = b where ‖E‖∞ 6 εm‖A‖∞ Gaussian elimination with partial pivoting and back substitution “gives exactly the right answer to nearly the right question.” — Trefethen and Bau L. Olson (UIUC) CS 257 September 22, 2006 17 / 33 Computational Stability An algorithm that gives the exact answer to a problem that is near to the original problem is said to be backward stable. Algorithms that are not backward stable will tend to amplify roundoff errors present in the original data. As a result, the solution produced by an algorithm that is not backward stable will not necessarily be the solution to a problem that is close to the original problem. Gaussian elimination without partial pivoting is not backward stable for arbitrary A. If A is symmetric and positive definite, then Gaussian elimination without pivoting in backward stable. L. Olson (UIUC) CS 257 September 22, 2006 18 / 33 The Residual Let x̂ be the numerical solution to Ax = b. x̂ , x (x is the exact solution) because of roundoff. The residual measures how close x̂ is to satisfying the original equation r = b − Ax̂ It is not hard to show that ‖x̂ − x‖ ‖x̂‖ 6 κ(A) ‖r‖ ‖b‖ Small ‖r‖ does not guarantee a small ‖x̂ − x‖. If κ(A) is large the x̂ returned by Gaussian elimination and back substitution (or any other solution method) is not guaranteed to be anywhere near the true solution to Ax = b. L. Olson (UIUC) CS 257 September 22, 2006 19 / 33 Summary of Limits to Numerical Solution of Ax = b 1 κ(A) indicates how close A is to being numerically singular 2 If κ(A) is “large”, A is ill-conditioned and even the best numerical algorithms will produce a solution, x̂ that cannot be guaranteed to be close to the true solution, x 3 In practice, Gaussian elimination with partial pivoting and back substitution produces a solution with a small residual r = b − Ax̂ even if κ(A) is large. L. Olson (UIUC) CS 257 September 22, 2006 22 / 33 The Backslash Operator Consider the scalar equation 5x = 20 =⇒ x = (5)−120 The extension to a system of equations is, of course Ax = b =⇒ x = A−1b where A−1b is the formal solution to Ax = b In MATLAB notation the system is solved with 1 x = A\b L. Olson (UIUC) CS 257 September 22, 2006 24 / 33 The Backslash Operator Given an n× n matrix A, and an n× 1 vector b the \ operator performs a sequence of tests on the A matrix. MATLAB attempts to solve the system with the method that gives the least roundoff and the fewest operations. When A is an n× n matrix: 1 MATLAB examines A to see if it is a permutation of a triangular system If so, the appropriate triangular solve is used. 2 MATLAB examines A to see if it appears to be symmetric and positive definite. If so, MATLAB attempts a Cholesky factorization and two triangular solves. 3 If the Cholesky factorization fails, or if A does not appear to be symmetric, MATLAB attempts an LU factorization and two triangular solves. L. Olson (UIUC) CS 257 September 22, 2006 26 / 33 Tridiagonal  d1 c1 a1 d2 c2 a2 d3 c3 . . . . . . . . . ai−1 di ci . . . . . . . . . . . . . . . . . . an−1 dn  Start forward elimination (without any special pivoting) 1 subtract a1/d1 times row 1 from row 2 2 this eliminates a1, changes d2 and does not touch c2 3 continuing: di = di − ( ai−1 di−1 ci−1 ) for i = 2 . . . n L. Olson (UIUC) CS 257 September 22, 2006 29 / 33 Tridiagonal  d̃1 c1 d̃2 c2 d̃3 c3 . . . . . . d̃i ci . . . . . . . . . . . . d̃n  This leaves an upper triangular (2-band). With back substitution: 1 xn = b̃n/d̃n 2 xn−1 = (1/d̃n−1)(b̃n−1 − cn−1xn) 3 xi = (1/d̃i)(b̃i − cixi+1) L. Olson (UIUC) CS 257 September 22, 2006 30 / 33 Tridiagonal Algorithm 1 input: n, a, d, c, b 2 for i = 2 to n 3 xmult = ai−1/di−1 4 di = di − xmult · ci−1 5 bi = bi − xmult · bi−1 6 end 7 x_{n} = b_{n}/d_{n} 8 for i = n − 1 down to 1 9 xi = (bi − cixi+1)/di 10 end L. Olson (UIUC) CS 257 September 22, 2006 32 / 33