Numerical Linear Algebra: Inverse & Rayleigh Iteration, Condition & Perturbation, Slides of Advanced Computer Architecture

Download Numerical Linear Algebra: Inverse & Rayleigh Iteration, Condition & Perturbation and more Slides Advanced Computer Architecture in PDF only on Docsity! Lecture 15 1 / 21Docsity.com Overview Inverse iteration Rayleigh quotient iteration Condition Perturbation Stability 2 / 21Docsity.com Inverse iteration and Rayleigh quotient iteration Inverse iteration: I one of the most valuable tools of numerical linear algebra I a standard method of calculating one or more eigenvectors of a matrix if the eigenvalues are already known Obtaining an eigenvalue estimate from an eigenvector estimate (the Rayleigh quotient) Obtaining an eigenvector estimate from an eigenvalue estimate (inverse iteration) Rayleigh quotient algorithm: Initialize v(0) randomly with ‖v(0)‖ = 1 λ(0) = (v(0))>Av(0) = corresponding Rayleigh quotient for k = 1, 2, . . . do Solve w = (A− λ(k−1)I )−1v(k−1) // apply (A− µ(k−1)I )−1 v(k) = w‖w‖ // normalize λ(k) = (v(k))>Av(k) // Rayleigh quotient end for 6 / 21Docsity.com Rayleigh quotient iteration Theorem Rayleigh quotient iteration converges to an eigenvalue/eigenvector pair for all except a set of measure zero of starting vectors v(0). When it converges, the convergence is ultimately cubic in the sense that if λJ is an eigenvalue of A and v(0) is sufficiently close to the eigenvector qJ , then ‖v(k+1) − (±qJ)‖ = O(‖v(k) − (±qJ)‖3) , and |λ(k+1) − λJ | = O(|λ(k) − λJ |3) as k →∞. The ± sign means that at each step k, one or the other choice of sign is to be taken, and then the indicated bound holds. 7 / 21Docsity.com Example Consider the symmetric matrix A = 2 1 11 3 1 1 1 4  and let v(0) = (1, 1, 1)>/ √ 3 the initial eigenvector estimate When Rayleigh quotient iteration is applied to A, the following values λ(k) are computed by the first 3 iterations: λ(0) = 5, λ(1) = 5.2131 . . . , λ(2) = 5.214319743184 . . . The actual value of the eigenvalue corresponding to the eigenvector closest to v(0) is λ = 5.214319743377 After three iterations, Rayleigh quotient iteration has produced a result accurate to 10 digits With three more iterations, the solution increases this figure to about 270 digits, if our machine precision were high enough 8 / 21Docsity.com Absolute condition number Let δx denote a small perturbation of x and write δf = f (x + δx)− f (x) Absolute condition number: κ̂ = κ̂(x) = lim δ→0 sup ‖δx‖≤δ ‖δf ‖ ‖δx‖ For most problems, it can be interpreted as a supremum over all infinitesimal perturbations δx, and can write simply as κ̂ = sup δx ‖δf ‖ ‖δx‖ If f is differentiable, we can evaluate the condition number by means of the derivative of f Let J(x) be the matrix whose i , j entry is the partial derivative ∂fi/∂xj evaluated at x, known as the Jacobian of f at x By the definition, δf ≈ J(x)δx with ‖δx‖ → 0, and the absolute condition number becomes κ̂ = ‖J(x)‖ where ‖J(x)‖ represents the norm of J(x) 11 / 21Docsity.com Relative condition number Relative condition number: κ = lim δ→0 sup ‖δx‖≤δ ( ‖δf ‖ ‖f (x)‖ / ‖δx‖ ‖x‖ ) Again, assume δx and δf are infinitesimal, κ = sup δx ( ‖δf ‖ ‖f (x)‖ / ‖δx‖ ‖x‖ ) If f is differentiable, we can express this in terms of Jacobian κ = ‖J(x)‖ ‖f (x)‖/‖x‖ Relative condition number is more important in numerical analysis A problem is well-conditioned if κ is small (e.g., 1, 10, 102), and ill-conditioned if κ is large (e.g., 106, 1016) 12 / 21Docsity.com Well-conditioned and ill-conditioned problems Consider the problem of obtaining the scalar x/2 from x ∈ C. The Jacobian of the function f : x → x/2 is the derivative J = f ′ = 1/2 κ = ‖J‖ ‖f (x)‖/‖x‖ = 1/2 (x/2)/x = 1. This problem is well-conditioned Consider x2 − 2x + 1 = (x − 1)2, with a double root x = 1. A small perturbation in the coefficient may lead to large change in the roots x2 − 2x + 0.9999 = (x − 0.99)(x − 1.01). In fact, the roots can change in proportion to the square root of the change in the coefficients, so in this case the Jacobian is infinite (the problem is not differentiable), and κ =∞ 13 / 21Docsity.com Condition of matrix-vector multiplication (cont’d) For certain choices of x, we have α = 1, and consequently κ = ‖A‖‖A−1‖ If we use 2-norm, this will occur when x is a multiple of a minimal right singular vector of A If A ∈ Cm×n with m ≥ n has full rank, the above equations hold with A−1 replaced by the pseudo-inverse A† = (AHA)−1AH ∈ Cn×m Given A, compute A−1b from input b? Mathematically, identical to the problem just considered except that A is repalced by A−1 16 / 21Docsity.com Condition of matrix-vector multiplication (cont’d) Theorem Let A ∈ Cm×n be nonsingular and consider the equation Ax = b. The problem of computing b given x, has condition number κ = ‖A‖‖x‖ ‖b‖ ≤ ‖A‖‖A−1‖ (1) with respect to perturbation of x. The problem of computing x given b, has condition number κ = ‖A−1‖‖b‖ ‖x‖ ≤ ‖A‖‖A−1‖ (2) with respect to perturbation of b. If we use 2-norm, then the equality holds in (1) if x is a multiple of a right singular vector of A corresponding to the minimal singular value σm. Likewise, the equality hold in (2) if b is a multiple of a left singular vector of A corresponding to the maximal singular value σ1. 17 / 21Docsity.com Condition number of a matrix The product ‖A‖‖A−1‖ is the condition number of A (relative to the norm ‖ · ‖), denoted by κ(A) κ(A) = ‖A‖‖A−1‖ The condition number is attached to a matrix, not a problem If κ(A) is small, A is said to be well-conditioned; if κ(A) is large, A is ill-conditioned. If A is singular, it is customary to write κ(A) =∞ If ‖ · ‖ = ‖ · ‖2, then ‖A‖ = σ1, and ‖A−1‖ = 1/σm, and thus κ(A) = σ1 σm in the 2-norm, which is the formula for computing 2-norm condition numbers of matrices. The ratio σ1/σm can be interpreted as the eccentricity of the hyperellipse that is the image of the unit sphere of Cm under A 18 / 21Docsity.com