This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Outline Scientiﬁc Computing: An Introductory Survey
Chapter 4 – Eigenvalue Problems
1 Eigenvalue Problems Prof. Michael T. Heath
Department of Computer Science University of Illinois at UrbanaChampaign 2 Existence, Uniqueness, and Conditioning 3 Computing Eigenvalues and Eigenvectors Copyright c 2002. Reproduction permitted for noncommercial, educational use only. Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Eigenvalue Problems Eigenvalues and Eigenvectors Geometric Interpretation 1 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Eigenvalue Problems Eigenvalues and Eigenvectors Geometric Interpretation 2 / 87 Eigenvalue Problems Eigenvalues and Eigenvectors Eigenvalue problems occur in many areas of science and engineering, such as structural analysis Eigenvalues are also important in analyzing numerical methods Theory and algorithms apply to complex matrices as well as real matrices With complex matrices, we use conjugate transpose, AH , instead of usual transpose, AT Standard eigenvalue problem : Given n × n matrix A, ﬁnd scalar λ and nonzero vector x such that Ax = λx λ is eigenvalue, and x is corresponding eigenvector λ may be complex even if A is real Spectrum = λ(A) = set of eigenvalues of A Spectral radius = ρ(A) = max{λ : λ ∈ λ(A)} Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Eigenvalue Problems Eigenvalues and Eigenvectors Geometric Interpretation 3 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Eigenvalue Problems Eigenvalues and Eigenvectors Geometric Interpretation 4 / 87 Geometric Interpretation Examples: Eigenvalues and Eigenvectors
A= 10 1 : λ1 = 1, x1 = , 02 0 11 1 : λ1 = 1, x1 = , 02 0 λ2 = 2, x2 = λ2 = 2, x2 = 0 1 1 1 1 −1 −1 1 i 1 Matrix expands or shrinks any vector lying in direction of eigenvector by scalar factor Expansion or contraction factor is given by corresponding eigenvalue λ Eigenvalues and eigenvectors decompose complicated behavior of general linear transformation into simpler actions A= A= A= A= 3 −1 1 : λ1 = 2, x1 = , −1 3 1 1.5 0.5 1 : λ1 = 2, x1 = , 0.5 1.5 1 λ2 = 4, x2 = λ2 = 1, x2 = λ2 = −i, x2 = 01 1 : λ1 = i, x1 = , −1 0 i √ where i = −1
5 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 6 / 87 Characteristic Polynomial
Equation Ax = λx is equivalent to (A − λI )x = 0 which has nonzero solution x if, and only if, its matrix is singular Eigenvalues of A are roots λi of characteristic polynomial det(A − λI ) = 0 in λ of degree n Fundamental Theorem of Algebra implies that n × n matrix A always has n eigenvalues, but they may not be real nor distinct Complex eigenvalues of real matrix occur in complex conjugate pairs: if α + iβ is eigenvalue of real matrix, then √ so is α − iβ , where i = −1
Michael T. Heath Scientiﬁc Computing 7 / 87 Example: Characteristic Polynomial
Characteristic polynomial of previous example matrix is det 3 −1 10 −λ −1 3 01 det 3 − λ −1 −1 3 − λ = = (3 − λ)(3 − λ) − (−1)(−1) = λ2 − 6λ + 8 = 0 so eigenvalues are given by √ 6 ± 36 − 32 λ= , or λ1 = 2, 2 λ2 = 4 Michael T. Heath Scientiﬁc Computing 8 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Characteristic Polynomial Relevant Properties of Matrices Conditioning Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Characteristic Polynomial Relevant Properties of Matrices Conditioning Companion Matrix
Monic polynomial p(λ) = c0 + c1 λ + · · · + cn−1 λn−1 + λn is characteristic polynomial of companion matrix 0 0 · · · 0 −c0 1 0 · · · 0 −c1 Cn = 0 1 · · · 0 −c2 . . .. . . .. . . . . .. . 0 0 · · · 1 −cn−1 Roots of polynomial of degree > 4 cannot always computed in ﬁnite number of steps So in general, computation of eigenvalues of matrices of order > 4 requires (theoretically inﬁnite) iterative process
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 9 / 87 Characteristic Polynomial, continued Computing eigenvalues using characteristic polynomial is not recommended because of
work in computing coefﬁcients of characteristic polynomial sensitivity of coefﬁcients of characteristic polynomial work in solving for roots of characteristic polynomial Characteristic polynomial is powerful theoretical tool but usually not useful computationally Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 10 / 87 Example: Characteristic Polynomial
Consider A= 1 1 √
mach Multiplicity and Diagonalizability
Multiplicity is number of times root appears when polynomial is written as product of linear factors Eigenvalue of multiplicity 1 is simple Defective matrix has eigenvalue of multiplicity k > 1 with fewer than k linearly independent corresponding eigenvectors Nondefective matrix A has n linearly independent eigenvectors, so it is diagonalizable X −1 AX = D where X is nonsingular matrix of eigenvectors
11 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 12 / 87 where is positive number slightly smaller than Exact eigenvalues of A are 1 + and 1 − Computing characteristic polynomial in ﬂoatingpoint arithmetic, we obtain det(A − λI ) = λ2 − 2λ + (1 − which has 1 as double root Thus, eigenvalues cannot be resolved by this method even though they are distinct in working precision
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 2 ) = λ 2 − 2λ + 1 Eigenspaces and Invariant Subspaces
Eigenvectors can be scaled arbitrarily: if Ax = λx, then A(γ x) = λ(γ x) for any scalar γ , so γ x is also eigenvector corresponding to λ Eigenvectors are usually normalized by requiring some norm of eigenvector to be 1 Eigenspace = Sλ = {x : Ax = λx} Subspace S of Rn (or Cn ) is invariant if AS ⊆ S For eigenvectors x1 · · · xp , span([x1 · · · xp ]) is invariant subspace Relevant Properties of Matrices
Properties of matrix A relevant to eigenvalue problems Property diagonal tridiagonal triangular Hessenberg Deﬁnition aij = 0 for i = j aij = 0 for i − j  > 1 aij = 0 for i > j (upper) aij = 0 for i < j (lower) aij = 0 for i > j + 1 (upper) aij = 0 for i < j − 1 (lower) AT A = AAT = I AH A = AAH = I A = AT A = AH AH A = AAH
Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 14 / 87 orthogonal unitary symmetric Hermitian normal
13 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Examples: Matrix Properties
Transpose: 12 34
T Examples, continued
Orthogonal:
H = 13 24 1 + i 1 + 2i 2 − i 2 − 2i = 1−i 2+i 1 − 2i 2 + 2 i Conjugate transpose: Symmetric: 12 23 13 24 √ √ 01 −1 0 √2/2 √2/2 , , 10 0 −1 − 2/2 2/2 √ √ i√2/2 √ 2/2 Unitary: − 2/2 −i 2/2 11 12 120 Normal: 0 1 2 201 Nonorthogonal: Nonnormal: 11 01
Michael T. Heath Scientiﬁc Computing 16 / 87 Nonsymmetric: Hermitian: 1 1+i 1−i 2 1 1+i 1+i 2
Michael T. Heath Scientiﬁc Computing 15 / 87 NonHermitian: Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Characteristic Polynomial Relevant Properties of Matrices Conditioning Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Characteristic Polynomial Relevant Properties of Matrices Conditioning Properties of Eigenvalue Problems
Properties of eigenvalue problem affecting choice of algorithm and software Are all eigenvalues needed, or only a few? Are only eigenvalues needed, or are corresponding eigenvectors also needed? Is matrix real or complex? Is matrix relatively small and dense, or large and sparse? Does matrix have any special properties, such as symmetry, or is it general matrix? Conditioning of Eigenvalue Problems Condition of eigenvalue problem is sensitivity of eigenvalues and eigenvectors to changes in matrix Conditioning of eigenvalue problem is not same as conditioning of solution to linear system for same matrix Different eigenvalues and eigenvectors are not necessarily equally sensitive to perturbations in matrix Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 17 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Characteristic Polynomial Relevant Properties of Matrices Conditioning 18 / 87 Conditioning of Eigenvalues
If µ is eigenvalue of perturbation A + E of nondefective matrix A, then µ − λk  ≤ cond2 (X ) E
2 Conditioning of Eigenvalues
If (A + E )(x + ∆x) = (λ + ∆λ)(x + ∆x), where λ is simple eigenvalue of A, then ∆λ y
2· x y H x 2 E 2 = 1 E cos(θ) 2 where λk is closest eigenvalue of A to µ and X is nonsingular matrix of eigenvectors of A Absolute condition number of eigenvalues is condition number of matrix of eigenvectors with respect to solving linear equations Eigenvalues may be sensitive if eigenvectors are nearly linearly dependent (i.e., matrix is nearly defective) For normal matrix (AH A = AAH ), eigenvectors are orthogonal, so eigenvalues are wellconditioned
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 19 / 87 where x and y are corresponding right and left eigenvectors and θ is angle between them For symmetric or Hermitian matrix, right and left eigenvectors are same, so cos(θ) = 1 and eigenvalues are inherently wellconditioned Eigenvalues of nonnormal matrices may be sensitive For multiple or closely clustered eigenvalues, corresponding eigenvectors may be sensitive
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 20 / 87 Problem Transformations
Shift : If Ax = λx and σ is any scalar, then (A − σ I )x = (λ − σ )x, so eigenvalues of shifted matrix are shifted eigenvalues of original matrix Inversion : If A is nonsingular and Ax = λx with x = 0, then λ = 0 and A−1 x = (1/λ)x, so eigenvalues of inverse are reciprocals of eigenvalues of original matrix Powers : If Ax = λx, then Ak x = λk x, so eigenvalues of power of matrix are same power of eigenvalues of original matrix Polynomial : If Ax = λx and p(t) is polynomial, then p(A)x = p(λ)x, so eigenvalues of polynomial in matrix are values of polynomial evaluated at eigenvalues of original matrix
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 21 / 87 Similarity Transformation
B is similar to A if there is nonsingular matrix T such that B = T −1 A T Then By = λy ⇒ T −1 AT y = λy ⇒ A(T y ) = λ(T y ) so A and B have same eigenvalues, and if y is eigenvector of B , then x = T y is eigenvector of A Similarity transformations preserve eigenvalues and eigenvectors are easily recovered Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 22 / 87 Example: Similarity Transformation
From eigenvalues and eigenvectors for previous example, 3 −1 −1 3 and hence 0.5 0.5 0.5 −0.5 3 −1 −1 3 1 1 20 = 1 −1 04 1 1 1 1 = 1 −1 1 −1 20 04 Diagonal Form
Eigenvalues of diagonal matrix are diagonal entries, and eigenvectors are columns of identity matrix Diagonal form is desirable in simplifying eigenvalue problems for general matrices by similarity transformations But not all matrices are diagonalizable by similarity transformation Closest one can get, in general, is Jordan form, which is nearly diagonal but may have some nonzero entries on ﬁrst superdiagonal, corresponding to one or more multiple eigenvalues So original matrix is similar to diagonal matrix, and eigenvectors form columns of similarity transformation matrix Michael T. Heath Scientiﬁc Computing 23 / 87 Michael T. Heath Scientiﬁc Computing 24 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Triangular Form
Any matrix can be transformed into triangular (Schur ) form by similarity, and eigenvalues of triangular matrix are diagonal entries Eigenvectors of triangular matrix less obvious, but still straightforward to compute If U11 u U13 0 vT A − λI = 0 O 0 U33 is triangular, then U11 y = u can be solved for y , so that y x = −1 0 is corresponding eigenvector
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 25 / 87 Block Triangular Form
If A11 A12 · · · A22 · · · A= .. . A 1p A 2p . . . App with square diagonal blocks, then
p λ(A) =
j =1 λ(Ajj ) so eigenvalue problem breaks into p smaller eigenvalue problems Real Schur form has 1 × 1 diagonal blocks corresponding to real eigenvalues and 2 × 2 diagonal blocks corresponding to pairs of complex conjugate eigenvalues
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 26 / 87 Forms Attainable by Similarity
A distinct eigenvalues real symmetric complex Hermitian normal arbitrary real arbitrary arbitrary T nonsingular orthogonal unitary unitary orthogonal unitary nonsingular B diagonal real diagonal real diagonal diagonal real block triangular (real Schur) upper triangular (Schur) almost diagonal (Jordan) Power Iteration
Simplest method for computing one eigenvalueeigenvector pair is power iteration, which repeatedly multiplies matrix times initial starting vector Assume A has unique eigenvalue of maximum modulus, say λ1 , with corresponding eigenvector v1 Then, starting from nonzero vector x0 , iteration scheme xk = Axk−1 converges to multiple of eigenvector v1 corresponding to dominant eigenvalue λ1 Given matrix A with indicated property, matrices B and T exist with indicated properties such that B = T −1 AT If B is diagonal or triangular, eigenvalues are its diagonal entries If B is diagonal, eigenvectors are columns of T
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 27 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 28 / 87 Convergence of Power Iteration
To see why power iteration converges to dominant eigenvector, express starting vector x0 as linear combination
n Example: Power Iteration
Ratio of values of given component of xk from one iteration to next converges to dominant eigenvalue λ1 1.5 0.5 0 For example, if A = and x0 = , we obtain 0.5 1.5 1 T k xk ratio 0 0.0 1.0 1 0.5 1.5 1.500 2 1.5 2.5 1.667 3 3.5 4.5 1.800 4 7.5 8.5 1.889 5 15.5 16.5 1.941 6 31.5 32.5 1.970 7 63.5 64.5 1.985 8 127.5 128.5 1.992 Ratio is converging to dominant eigenvalue, which is 2
29 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 30 / 87 x0 =
i=1 αi vi where vi are eigenvectors of A Then xk = Axk−1 = A2 xk−2 = · · · = Ak x0 =
n n λk αi vi = λk i 1
i=1 α1 v1 +
i=2 (λi /λ1 )k αi vi Since λi /λ1  < 1 for i > 1, successively higher powers go to zero, leaving only component corresponding to v1
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods Limitations of Power Iteration
Power iteration can fail for various reasons Starting vector may have no component in dominant eigenvector v1 (i.e., α1 = 0) — not problem in practice because rounding error usually introduces such component in any case There may be more than one eigenvalue having same (maximum) modulus, in which case iteration may converge to linear combination of corresponding eigenvectors For real matrix and starting vector, iteration can never converge to complex vector Normalized Power Iteration Geometric growth of components at each iteration risks eventual overﬂow (or underﬂow if λ1 < 1) Approximate eigenvector should be normalized at each iteration, say, by requiring its largest component to be 1 in modulus, giving iteration scheme yk = Axk−1 xk = yk / yk With normalization, yk
∞ ∞ ∞ → λ1 , and xk → v1 / v1 Michael T. Heath Scientiﬁc Computing 31 / 87 Michael T. Heath Scientiﬁc Computing 32 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Example: Normalized Power Iteration
Repeating previous example with normalized scheme, k 0 1 2 3 4 5 6 7 8 xT k 0.000 0.333 0.600 0.778 0.882 0.939 0.969 0.984 0.992 yk 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
∞ Geometric Interpretation
Behavior of power iteration depicted geometrically 1.500 1.667 1.800 1.889 1.941 1.970 1.985 1.992 Initial vector x0 = v1 + v2 contains equal components in eigenvectors v1 and v2 (dashed arrows) Repeated multiplication by A causes component in v1 (corresponding to larger eigenvalue, 2) to dominate, so sequence of vectors xk converges to v1
33 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 34 / 87 < interactive example >
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods Power Iteration with Shift
Convergence rate of power iteration depends on ratio λ2 /λ1 , where λ2 is eigenvalue having second largest modulus May be possible to choose shift, A − σ I , such that λ2 − σ λ2 < λ1 − σ λ1 so convergence is accelerated Shift must then be added to result to obtain eigenvalue of original matrix Example: Power Iteration with Shift In earlier example, for instance, if we pick shift of σ = 1, (which is equal to other eigenvalue) then ratio becomes zero and method converges in one iteration In general, we would not be able to make such fortuitous choice, but shifts can still be extremely useful in some contexts, as we will see later Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 35 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 36 / 87 Inverse Iteration
If smallest eigenvalue of matrix required rather than largest, can make use of fact that eigenvalues of A−1 are reciprocals of those of A, so smallest eigenvalue of A is reciprocal of largest eigenvalue of A−1 This leads to inverse iteration scheme Ayk = xk−1 xk = yk / yk
∞ Inverse Iteration, continued Inverse iteration converges to eigenvector corresponding to smallest eigenvalue of A Eigenvalue obtained is dominant eigenvalue of A−1 , and hence its reciprocal is smallest eigenvalue of A in modulus which is equivalent to power iteration applied to A−1 Inverse of A not computed explicitly, but factorization of A used to solve system of linear equations at each iteration
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 37 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 38 / 87 Example: Inverse Iteration
Applying inverse iteration to previous example to compute smallest eigenvalue yields sequence k xT yk ∞ k 0 0.000 1.0 1 −0.333 1.0 0.750 2 −0.600 1.0 0.833 3 −0.778 1.0 0.900 4 −0.882 1.0 0.944 5 −0.939 1.0 0.971 6 −0.969 1.0 0.985 which is indeed converging to 1 (which is its own reciprocal in this case) < interactive example >
Michael T. Heath Scientiﬁc Computing 39 / 87 Inverse Iteration with Shift
As before, shifting strategy, working with A − σ I for some scalar σ , can greatly improve convergence Inverse iteration is particularly useful for computing eigenvector corresponding to approximate eigenvalue, since it converges rapidly when applied to shifted matrix A − λI , where λ is approximate eigenvalue Inverse iteration is also useful for computing eigenvalue closest to given value β , since if β is used as shift, then desired eigenvalue corresponds to smallest eigenvalue of shifted matrix Michael T. Heath Scientiﬁc Computing 40 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Rayleigh Quotient
Given approximate eigenvector x for real matrix A, determining best estimate for corresponding eigenvalue λ can be considered as n × 1 linear least squares approximation problem xλ ∼ Ax = From normal equation xT xλ = xT Ax, least squares solution is given by xT Ax λ= T xx This quantity, known as Rayleigh quotient, has many useful properties
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 41 / 87 Example: Rayleigh Quotient
Rayleigh quotient can accelerate convergence of iterative methods such as power iteration, since Rayleigh quotient xT Axk /xT xk gives better approximation to eigenvalue at k k iteration k than does basic method alone For previous example using power iteration, value of Rayleigh quotient at each iteration is shown below k 0 1 2 3 4 5 6 xT k 0.000 0.333 0.600 0.778 0.882 0.939 0.969 yk 1.0 1.0 1.0 1.0 1.0 1.0 1.0
∞ xT Axk /xT xk k k 1.500 1.800 1.941 1.985 1.996 1.999
42 / 87 1.500 1.667 1.800 1.889 1.941 1.970 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods Rayleigh Quotient Iteration
Given approximate eigenvector, Rayleigh quotient yields good estimate for corresponding eigenvalue Conversely, inverse iteration converges rapidly to eigenvector if approximate eigenvalue is used as shift, with one iteration often sufﬁcing These two ideas combined in Rayleigh quotient iteration σk = xT Axk /xT xk k k (A − σk I )yk+1 = xk xk+1 = yk+1 / yk+1 starting from given nonzero vector x0
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 43 / 87 Rayleigh Quotient Iteration, continued Rayleigh quotient iteration is especially effective for symmetric matrices and usually converges very rapidly Using different shift at each iteration means matrix must be refactored each time to solve linear system, so cost per iteration is high unless matrix has special form that makes factorization easy Same idea also works for complex matrices, for which transpose is replaced by conjugate transpose, so Rayleigh quotient becomes xH Ax/xH x ∞ Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 44 / 87 Example: Rayleigh Quotient Iteration Deﬂation
After eigenvalue λ1 and corresponding eigenvector x1 have been computed, then additional eigenvalues λ2 , . . . , λn of A can be computed by deﬂation, which effectively removes known eigenvalue Let H be any nonsingular matrix such that Hx1 = αe1 , scalar multiple of ﬁrst column of identity matrix (Householder transformation is good choice for H ) Then similarity transformation determined by H transforms A into form λ bT HAH −1 = 1 0B where B is matrix of order n − 1 having eigenvalues λ2 , . . . , λ n Using same matrix as previous examples and randomly chosen starting vector x0 , Rayleigh quotient iteration converges in two iterations k 0 1 2 xT k 0.807 0.397 0.924 1.000 1.000 1.000 σk 1.896 1.998 2.000 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 45 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 46 / 87 Deﬂation, continued Deﬂation, continued
Alternative approach lets u1 be any vector such that uT x1 = λ1 1 Then A − x1 uT has eigenvalues 0, λ2 , . . . , λn 1 Possible choices for u1 include
u1 = λ1 x1 , if A is symmetric and x1 is normalized so that x1 2 = 1 u1 = λ1 y1 , where y1 is corresponding left eigenvector (i.e., T AT y1 = λ1 y1 ) normalized so that y1 x1 = 1 u1 = AT ek , if x1 is normalized so that x1 ∞ = 1 and k th component of x1 is 1 Thus, we can work with B to compute next eigenvalue λ2 Moreover, if y2 is eigenvector of B corresponding to λ2 , then bT y2 α x2 = H −1 , where α = y2 λ2 − λ1 is eigenvector corresponding to λ2 for original matrix A, provided λ1 = λ2 Process can be repeated to ﬁnd additional eigenvalues and eigenvectors Michael T. Heath Scientiﬁc Computing 47 / 87 Michael T. Heath Scientiﬁc Computing 48 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Simultaneous Iteration
Simplest method for computing many eigenvalueeigenvector pairs is simultaneous iteration, which repeatedly multiplies matrix times matrix of initial starting vectors Starting from n × p matrix X0 of rank p, iteration scheme is Xk = AXk−1 span(Xk ) converges to invariant subspace determined by p largest eigenvalues of A, provided λp  > λp+1  Also called subspace iteration Orthogonal Iteration
As with power iteration, normalization is needed with simultaneous iteration Each column of Xk converges to dominant eigenvector, so columns of Xk become increasingly illconditioned basis for span(Xk ) Both issues can be addressed by computing QR factorization at each iteration ˆ Qk Rk = Xk−1 ˆ Xk = AQk ˆ where Qk Rk is reduced QR factorization of Xk−1 This orthogonal iteration converges to block triangular form, and leading block is triangular if moduli of consecutive eigenvalues are distinct
49 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 50 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods QR Iteration
For p = n and X0 = I , matrices ˆ ˆ Ak = QH AQk k generated by orthogonal iteration converge to triangular or block triangular form, yielding all eigenvalues of A QR iteration computes successive matrices Ak without forming above product explicitly Starting with A0 = A, at iteration k compute QR factorization Qk Rk = Ak−1 and form reverse product Ak = Rk Qk
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 51 / 87 QR Iteration, continued
Successive matrices Ak are unitarily similar to each other Ak = Rk Qk = QH Ak−1 Qk k Diagonal entries (or eigenvalues of diagonal blocks) of Ak converge to eigenvalues of A Product of orthogonal matrices Qk converges to matrix of corresponding eigenvectors If A is symmetric, then symmetry is preserved by QR iteration, so Ak converge to matrix that is both triangular and symmetric, hence diagonal Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 52 / 87 Example: QR Iteration
72 24 Compute QR factorization Let A0 = .962 −.275 A0 = Q1 R1 = .275 .962 and form reverse product A1 = R1 Q1 = 7.83 .906 .906 3.17 7.28 3.02 0 3.30 QR Iteration with Shifts Convergence rate of QR iteration can be accelerated by incorporating shifts Qk Rk = Ak−1 − σk I Ak = Rk Qk + σk I where σk is rough approximation to eigenvalue Good shift can be determined by computing eigenvalues of 2 × 2 submatrix in lower right corner of matrix Offdiagonal entries are now smaller, and diagonal entries closer to eigenvalues, 8 and 3 Process continues until matrix is within tolerance of being diagonal, and diagonal entries then closely approximate eigenvalues
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 53 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 54 / 87 Example: QR Iteration with Shifts
Repeat previous example, but with shift of σ1 = 4, which is lower right corner entry of matrix We compute QR factorization A0 − σ1 I = Q1 R1 = .832 .555 .555 −.832 3.61 1.66 0 1.11 Preliminary Reduction
Efﬁciency of QR iteration can be enhanced by ﬁrst transforming matrix as close to triangular form as possible before beginning iterations Hessenberg matrix is triangular except for one additional nonzero diagonal immediately adjacent to main diagonal Any matrix can be reduced to Hessenberg form in ﬁnite number of steps by orthogonal similarity transformation, for example using Householder transformations Symmetric Hessenberg matrix is tridiagonal Hessenberg or tridiagonal form is preserved during successive QR iterations
55 / 87 Michael T. Heath Scientiﬁc Computing 56 / 87 and form reverse product, adding back shift to obtain A1 = R1 Q1 + σ1 I = 7.92 .615 .615 3.08 After one iteration, offdiagonal entries smaller compared with unshifted algorithm, and eigenvalues closer approximations to eigenvalues < interactive example >
Michael T. Heath Scientiﬁc Computing Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Preliminary Reduction, continued
Advantages of initial reduction to upper Hessenberg or tridiagonal form Work per QR iteration is reduced from O(n3 ) to O(n2 ) for general matrix or O(n) for symmetric matrix Fewer QR iterations are required because matrix nearly triangular (or diagonal) already If any zero entries on ﬁrst subdiagonal, then matrix is block triangular and problem can be broken into two or more smaller subproblems Preliminary Reduction, continued
QR iteration is implemented in twostages symmetric −→ tridiagonal or general −→ Hessenberg −→ triangular −→ diagonal Preliminary reduction requires deﬁnite number of steps, whereas subsequent iterative stage continues until convergence In practice only modest number of iterations usually required, so much of work is in preliminary reduction Cost of accumulating eigenvectors, if needed, dominates total cost Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 57 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 58 / 87 Cost of QR Iteration Krylov Subspace Methods
Krylov subspace methods reduce matrix to Hessenberg (or tridiagonal) form using only matrixvector multiplication For arbitrary starting vector x0 , if Kk = x0 Ax0 · · · then
− Kn 1 AKn = Cn Approximate overall cost of preliminary reduction and QR iteration, counting both additions and multiplications Symmetric matrices
for eigenvalues only 9n for eigenvalues and eigenvectors
43 3n 3 Ak−1 x0 General matrices
10n3 for eigenvalues only 25n3 for eigenvalues and eigenvectors where Cn is upper Hessenberg (in fact, companion matrix) To obtain better conditioned basis for span(Kn ), compute QR factorization Qn Rn = Kn so that
− QH AQn = Rn Cn Rn 1 ≡ H n with H upper Hessenberg
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 59 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 60 / 87 Krylov Subspace Methods
Equating k th columns on each side of equation AQn = Qn H yields recurrence Aqk = h1k q1 + · · · + hkk qk + hk+1,k qk+1 relating qk+1 to preceding vectors q1 , . . . , qk
H Premultiplying by qj and using orthonormality, H hjk = qj Aqk , Arnoldi Iteration
x0 = arbitrary nonzero starting vector q1 = x0 / x0 2 for k = 1, 2, . . . uk = Aqk for j = 1 to k H hjk = qj uk uk = uk − hjk qj end hk+1,k = uk 2 if hk+1,k = 0 then stop qk+1 = uk /hk+1,k end j = 1, . . . , k These relationships yield Arnoldi iteration, which produces unitary matrix Qn and upper Hessenberg matrix Hn column by column using only matrixvector multiplication by A and inner products of vectors
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 61 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 62 / 87 Arnoldi Iteration, continued
If Qk = q1 · · · then Hk = QH AQk k is upper Hessenberg matrix Eigenvalues of Hk , called Ritz values, are approximate eigenvalues of A, and Ritz vectors given by Qk y , where y is eigenvector of Hk , are corresponding approximate eigenvectors of A Eigenvalues of Hk must be computed by another method, such as QR iteration, but this is easier problem if k n
Michael T. Heath Scientiﬁc Computing 63 / 87 Arnoldi Iteration, continued qk Arnoldi iteration fairly expensive in work and storage because each new vector qk must be orthogonalized against all previous columns of Qk , and all must be stored for that purpose. So Arnoldi process usually restarted periodically with carefully chosen starting vector Ritz values and vectors produced are often good approximations to eigenvalues and eigenvectors of A after relatively few iterations Michael T. Heath Scientiﬁc Computing 64 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Lanczos Iteration
Work and storage costs drop dramatically if matrix is symmetric or Hermitian, since recurrence then has only three terms and Hk is tridiagonal (so usually denoted Tk ) q0 = 0 β0 = 0 x0 = arbitrary nonzero starting vector q1 = x0 / x0 2 for k = 1, 2, . . . uk = Aqk H αk = qk uk uk = uk − βk−1 qk−1 − αk qk βk = uk 2 if βk = 0 then stop qk+1 = uk /βk end
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 65 / 87 Lanczos Iteration, continued
αk and βk are diagonal and subdiagonal entries of symmetric tridiagonal matrix Tk As with Arnoldi, Lanczos iteration does not produce eigenvalues and eigenvectors directly, but only tridiagonal matrix Tk , whose eigenvalues and eigenvectors must be computed by another method to obtain Ritz values and vectors If βk = 0, then algorithm appears to break down, but in that case invariant subspace has already been identiﬁed (i.e., Ritz values and vectors are already exact at that point) Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 66 / 87 Lanczos Iteration, continued
In principle, if Lanczos algorithm were run until k = n, resulting tridiagonal matrix would be orthogonally similar to A In practice, rounding error causes loss of orthogonality, invalidating this expectation Problem can be overcome by reorthogonalizing vectors as needed, but expense can be substantial Alternatively, can ignore problem, in which case algorithm still produces good eigenvalue approximations, but multiple copies of some eigenvalues may be generated Krylov Subspace Methods, continued
Great virtue of Arnoldi and Lanczos iterations is their ability to produce good approximations to extreme eigenvalues for k n Moreover, they require only one matrixvector multiplication by A per step and little auxiliary storage, so are ideally suited to large sparse matrices If eigenvalues are needed in middle of spectrum, say near σ , then algorithm can be applied to matrix (A − σ I )−1 , assuming it is practical to solve systems of form (A − σ I )x = y Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 67 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 68 / 87 Example: Lanczos Iteration
For 29 × 29 symmetric matrix with eigenvalues 1, . . . , 29, behavior of Lanczos iteration is shown below Jacobi Method
One of oldest methods for computing eigenvalues is Jacobi method, which uses similarity transformation based on plane rotations Sequence of plane rotations chosen to annihilate symmetric pairs of matrix entries, eventually converging to diagonal form Choice of plane rotation slightly more complicated than in Givens method for QR factorization To annihilate given offdiagonal pair, choose c and s so that J T AJ = = is diagonal c −s s c ab bd cs −s c c2 a − 2csb + s2 d c2 b + cs(a − d) − s2 b c2 b + cs(a − d) − s2 b c2 d + 2csb + s2 a
Michael T. Heath Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 70 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 69 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Jacobi Method, continued
Transformed matrix diagonal if c2 b + cs(a − d) − s2 b = 0 Dividing both sides by c2 b, we obtain s (a − d) s2 − 2 =0 c b c Making substitution t = s/c, we get quadratic equation 1+ (a − d) − t2 = 0 b for tangent t of √ angle of rotation, from which we can recover c = 1/( 1 + t2 ) and s = c · t 1+t Advantageous numerically to use root of smaller magnitude
Michael T. Heath Scientiﬁc Computing 71 / 87 Example: Plane Rotation
Consider 2 × 2 matrix A= 12 21 Quadratic equation for tangent reduces to t2 = 1, so t = ±1 Two roots of same magnitude, so we arbitrarily choose √ √ t = −1, which yields c = 1/ 2 and s = −1/ 2 Resulting plane rotation J gives √ √ √ √ 3 0 1/√2 1/√2 1 2 1/√2 −1/√2 = J T AJ = 0 −1 −1/ 2 1/ 2 2 1 1/ 2 1/ 2
Michael T. Heath Scientiﬁc Computing 72 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Jacobi Method, continued
Starting with symmetric matrix A0 = A, each iteration has form T Ak+1 = Jk Ak Jk where Jk is plane rotation chosen to annihilate a symmetric pair of entries in Ak Plane rotations repeatedly applied from both sides in systematic sweeps through matrix until offdiagonal mass of matrix is reduced to within some tolerance of zero Resulting diagonal matrix orthogonally similar to original matrix, so diagonal entries are eigenvalues, and eigenvectors are given by product of plane rotations
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 73 / 87 Jacobi Method, continued
Jacobi method is reliable, simple to program, and capable of high accuracy, but converges rather slowly and difﬁcult to generalize beyond symmetric matrices Except for small problems, more modern methods usually require 5 to 10 times less work than Jacobi One source of inefﬁciency is that previously annihilated entries can subsequently become nonzero again, thereby requiring repeated annihilation Newer methods such as QR iteration preserve zero entries introduced into matrix Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 74 / 87 Example: Jacobi Method 102 0 2 1 Let A0 = 211 First annihilate (1,3) and (3,1) entries using rotation 0.707 0 −0.707 0 1 0 J0 = 0.707 0 0.707 to obtain 3 0.707 0 T 2 0.707 A1 = J0 A0 J0 = 0.707 0 0.707 −1
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 75 / 87 Example, continued Next annihilate (1,2) and (2,1) entries using rotation 0.888 −0.460 0 0.888 0 J1 = 0.460 0 0 1 to obtain A2 =
T J 1 A 1 J1 3.366 0 0.325 1.634 0.628 = 0 0.325 0.628 −1 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 76 / 87 Example, continued Example, continued
Beginning new sweep, again annihilate (1,3) and (3,1) entries using rotation 0.998 0 −0.070 1 0 J3 = 0 0.070 0 0.998 to obtain A4 =
T J 3 A 3 J3 Next annihilate (2,3) and (3,2) entries using rotation 1 0 0 0 0.975 −0.221 J2 = 0 0.221 0.975 to obtain 3.366 0.072 0.317 T 0 A3 = J2 A2 J2 = 0.072 1.776 0.317 0 −1.142 3.388 0.072 0 0.072 1.776 −0.005 = 0 −0.005 −1.164 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 77 / 87 Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 78 / 87 Example, continued Bisection or SpectrumSlicing
For real symmetric matrix, can determine how many eigenvalues are less than given real number σ Process continues until offdiagonal entries reduced to as small as desired Result is diagonal matrix orthogonally similar to original matrix, with the orthogonal similarity transformation given by product of plane rotations < interactive example > By systematically choosing various values for σ (slicing spectrum at σ ) and monitoring resulting count, any eigenvalue can be isolated as accurately as desired For example, symmetric indeﬁnite factorization A = LDLT makes inertia (numbers of positive, negative, and zero eigenvalues) of symmetric matrix A easy to determine By applying factorization to matrix A − σ I for various values of σ , individual eigenvalues can be isolated as accurately as desired using interval bisection technique Michael T. Heath Scientiﬁc Computing 79 / 87 Michael T. Heath Scientiﬁc Computing 80 / 87 Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Problem Transformations Power Iteration and Variants Other Methods Sturm Sequence
Another spectrumslicing method for computing individual eigenvalues is based on Sturm sequence property of symmetric matrices Let A be symmetric matrix and let pr (σ ) denote determinant of leading principal minor of order r of A − σ I Then zeros of pr (σ ) strictly separate those of pr−1 (σ ), and number of agreements in sign of successive members of sequence pr (σ ), for r = 1, . . . , n, equals number of eigenvalues of A strictly greater than σ Determinants pr (σ ) are easy to compute if A is transformed to tridiagonal form before applying Sturm sequence technique
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 81 / 87 DivideandConquer Method
Express symmetric tridiagonal matrix T as T= T1 O + β uuT O T2 Can now compute eigenvalues and eigenvectors of smaller matrices T1 and T2 To relate these back to eigenvalues and eigenvectors of original matrix requires solution of secular equation, which can be done reliably and efﬁciently Applying this approach recursively yields divideandconquer algorithm for symmetric tridiagonal eigenproblems
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 82 / 87 Relatively Robust Representation
With conventional methods, cost of computing eigenvalues of symmetric tridiagonal matrix is O(n2 ), but if orthogonal eigenvectors are also computed, then cost rises to O(n3 ) Another possibility is to compute eigenvalues ﬁrst at O(n2 ) cost, and then compute corresponding eigenvectors separately using inverse iteration with computed eigenvalues as shifts Key to making this idea work is computing eigenvalues and corresponding eigenvectors to very high relative accuracy so that expensive explicit orthogonalization of eigenvectors is not needed RRR algorithm exploits this approach to produce eigenvalues and orthogonal eigenvectors at O(n2 ) cost
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 83 / 87 Generalized Eigenvalue Problems
Generalized eigenvalue problem has form Ax = λBx where A and B are given n × n matrices If either A or B is nonsingular, then generalized eigenvalue problem can be converted to standard eigenvalue problem, either (B −1 A)x = λx or (A−1 B )x = (1/λ)x This is not recommended, since it may cause
loss of accuracy due to rounding error loss of symmetry if A and B are symmetric Better alternative for generalized eigenvalue problems is QZ algorithm
Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 84 / 87 QZ Algorithm
If A and B are triangular, then eigenvalues are given by λi = aii /bii , for bii = 0 QZ algorithm reduces A and B simultaneously to upper triangular form by orthogonal transformations First, B is reduced to upper triangular form by orthogonal transformation from left, which is also applied to A Next, transformed A is reduced to upper Hessenberg form by orthogonal transformation from left, while maintaining triangular form of B , which requires additional transformations from right QZ Algorithm, continued Finally, analogous to QR iteration, A is reduced to triangular form while still maintaining triangular form of B , which again requires transformations from both sides Eigenvalues can now be determined from mutually triangular form, and eigenvectors can be recovered from products of left and right transformations, denoted by Q and Z Michael T. Heath Eigenvalue Problems Existence, Uniqueness, and Conditioning Computing Eigenvalues and Eigenvectors Scientiﬁc Computing Problem Transformations Power Iteration and Variants Other Methods 85 / 87 Michael T. Heath Scientiﬁc Computing 86 / 87 Computing SVD
Singular values of A are nonnegative square roots of eigenvalues of AT A, and columns of U and V are orthonormal eigenvectors of AAT and AT A, respectively Algorithms for computing SVD work directly with A, without forming AAT or AT A, to avoid loss of information associated with forming these matrix products explicitly SVD is usually computed by variant of QR iteration, with A ﬁrst reduced to bidiagonal form by orthogonal transformations, then remaining offdiagonal entries are annihilated iteratively SVD can also be computed by variant of Jacobi method, which can be useful on parallel computers or if matrix has special structure
Michael T. Heath Scientiﬁc Computing 87 / 87 ...
View
Full Document
 Fall '08
 Layton,A
 Eigenvectors, Vectors, Matrices, Michael T. Heath, Conditioning Computing Eigenvalues

Click to edit the document details