Direct methods for symmetric eigenvalue Algorithms problems - PowerPoint PPT Presentation

CES 703/5 Imre P´ olik Outline Introduction Direct methods for symmetric eigenvalue Algorithms problems Eigenvalues of a tridiagonal matrix Conclusion Literature Imre P´ olik, PhD McMaster University School of Computational Engineering and Science February 4, 2008

CES 703/5 Outline Imre P´ olik Introduction 1 Outline Theoretical background Introduction Posing the question Algorithms Perturbation theory Eigenvalues of a tridiagonal matrix Algorithms 2 Conclusion Direct and iterative methods Literature Transformation to tridiagonal form Counting the eigenvalues Finding the eigenvectors Eigenvalues of a tridiagonal matrix 3 QR method Divide and conquer Bisection with inverse iteration Jacobi method Conclusion 4

CES 703/5 Theoretical background ( A is real, symmetric) Imre P´ olik Eigenvalue: Ax = λx Outline A − λI is singular Introduction Theoretical n eigenvalues (with multiplicities) background Posing the question Perturbation theory Characteristic polynomial: p ( λ ) = det( A − λI ) Algorithms If Q T Q = I then A and Q T AQ have the same Eigenvalues of a tridiagonal matrix eigenvalues Conclusion If P is nonsingular then A and P T AP have the same Literature inertia Eigendecomposition ( Av i = λ i v i ): n � λ i v i v T A = i i =1 The v i s should be orthogonal

CES 703/5 Posing the question Imre P´ olik Find all the eigenvalues Outline Find some eigenvalues Introduction Theoretical Find largest/smallest eigenvalue(s) background Posing the question Perturbation theory Find largest/smallest magnitude eigenvalue(s) Algorithms Find eigenvalue(s) around a given number Eigenvalues of a tridiagonal matrix Find the corresponding eigenvectors Conclusion Find the eigenvalues in [ α, β ] Literature Count the eigenvalues in [ α, β ]

CES 703/5 Perturbation theory Imre P´ olik Eigenvalues are well-conditioned (Weyl) Outline Introduction | λ i ( A ) − λ i ( A + E ) | ≤ � E � 2 Theoretical background Posing the question Perturbation theory Eigenvectors can be ill-conditioned Algorithms     Eigenvalues of a 2 0 0 2 0 0 tridiagonal matrix 0 1 0 1 0 ε     Conclusion 0 ε 1 0 0 1 + ε Literature Problem: distance between eigenvalues Problem: small eigenvalues Shifting: λ ( A − αI ) = λ ( A ) − α

CES 703/5 Direct and iterative methods Imre P´ olik Eigenvalues are roots of polynomials Outline Every method has to iterate Introduction Direct Algorithms Direct and iterative transforms the matrix to a simpler form (tridiagonal) methods Transformation to finds all the eigenvalues and eigenvectors tridiagonal form Counting the (almost) always works eigenvalues Finding the eigenvectors converges quickly Eigenvalues of a Iterative tridiagonal matrix only multiplication by A or A T Conclusion does not access the elements of A Literature only a few eigenvectors/values

CES 703/5 Tridiagonalization (Hessenberg reduction) Imre P´ olik Householder reflections Outline     ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ Introduction � 1 � 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗     Algorithms =     Direct and iterative I − 2 uu T 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗     methods Transformation to � �� � ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ tridiagonal form Q Counting the � �� � eigenvalues A Finding the eigenvectors With this Q , due to symmetry: Eigenvalues of a tridiagonal matrix   Conclusion ∗ ∗ 0 0 Literature ∗ ∗ ∗ ∗   QAQ =   0 ∗ ∗ ∗   0 ∗ ∗ ∗ 3 n 3 if Q is needed for the eigenvectors Cost 4 3 n 3 , or 8 LAPACK: ssytrd , Matlab: hess

CES 703/5 Count the eigenvalues in [ α, β ) Imre P´ olik Factorize A − αI = L α D α L T α Outline Number of negative eigenvalues in D α gives the number Introduction of eigenvalues of A that are smaller than α Algorithms Direct and iterative Factorize A − βI = L β D β L T methods β Transformation to tridiagonal form Number of negative eigenvalues in D β gives the number Counting the eigenvalues Finding the of eigenvalues of A that are smaller than β eigenvectors Eigenvalues of a The difference gives the number of eigenvalues in [ α, β ) tridiagonal matrix Cost: two LDL factorizations (can exploit structure) Conclusion Literature

CES 703/5 Finding the eigenvectors Imre P´ olik Solve ( A − λI ) v = 0 Outline λ is only approximate Introduction extremely ill-conditioned system Algorithms only if λ is an exact eigenvalue Direct and iterative methods Inverse iteration Transformation to tridiagonal form Counting the repeat eigenvalues Finding the eigenvectors w i +1 = ( A − αI ) − 1 v i Eigenvalues of a tridiagonal matrix w i +1 v i +1 = Conclusion � w i +1 � 2 Literature fast convergence can be modified to find the eigenvalues, too

CES 703/5 Eigenvalues of a tridiagonal matrix Imre P´ olik Outline Introduction 1 Introduction Algorithms Algorithms 2 Eigenvalues of a tridiagonal matrix QR method Divide and conquer Eigenvalues of a tridiagonal matrix 3 Bisection with inverse iteration QR method Jacobi method Conclusion Divide and conquer Literature Bisection with inverse iteration Jacobi method Conclusion 4

CES 703/5 QR method Imre P´ olik Factor A i = Q i R i , set A i +1 = R i Q i , repeat Outline Shifting: A i − α i I = Q i R i , A i +1 = R i Q i + α i I Introduction If α i is close to an eigenvalue then convergence is fast Algorithms Eigenvalues of a (quadratic) tridiagonal matrix QR method Fastest method for eigenvalues Divide and conquer Bisection with inverse iteration Fastest method for eigenvectors if n ≤ 25 Jacobi method 3 n 3 + O ( n 2 ) , cubic! convergence Total cost: 4 Conclusion Literature LAPACK: xSTEQR , xSTERF , Matlab: eig

CES 703/5 QR method – Choosing the shift Imre P´ olik Outline   Introduction a 1 c 1 Algorithms c 1 a 2 c 2     Eigenvalues of a ... ... ...   A = tridiagonal matrix   QR method   c n − 2 a n − 1 c n − 1 Divide and conquer   Bisection with inverse iteration c n − 1 a n Jacobi method Conclusion Literature Single shift: α = a n Wilkinson shift: � a n − 1 � c n − 1 α is the eigenvalue of closest to a n c n − 1 a n

CES 703/5 Divide and conquer Imre P´ olik Partition the tridiagonal matrix into two halves Outline Solve the smaller problems recursively Introduction Assemble the results Algorithms Eigenvalues of a Cost: O ( n 2 . 3 ) on average tridiagonal matrix QR method Fastest method for eigenvectors if n ≥ 25 Divide and conquer Bisection with inverse iteration Not easy to implement in a stable way Jacobi method Conclusion Efficiently parallelizable Literature LAPACK: xSTEVD

CES 703/5 Bisection with inverse iteration Imre P´ olik Count the eigenvalues in [ α, β ) Outline Bisect the interval, count the eigenvalues in the Introduction subintervals, then repeat Algorithms LDL T is too expensive in general ⇒ Eigenvalues of a tridiagonal matrix QR method Transform A to tridiagonal form first Divide and conquer Bisection with inverse iteration Cost: O ( nk ) for k eigenvalues Jacobi method Excellent parallelization Conclusion Literature Only linear convergence (due to bisection) Compute eigenvectors with inverse iteration LAPACK: xSTEBZ , xSTEIN

CES 703/5 Jacobi method Imre P´ olik Transform A to (almost) diagonal form Outline Zero out a jk , j � = k with Givens rotation Introduction Algorithms � T � � � � � � � cos θ − sin θ A jj A jk cos θ − sin θ λ 1 0 = Eigenvalues of a sin θ cos θ A kj A kk sin θ cos θ 0 λ 2 tridiagonal matrix QR method Costs O ( n ) flops per rotation Divide and conquer Bisection with inverse iteration How to choose A jk ? Jacobi method Conclusion �� A 2 Literature off( A ) := jk j<k off( A ′ ) 2 = off( A ) 2 − A 2 jk Always the largest A jk : too expensive Row-wise: cheap but still efficient (5-10 sweeps) In general Jacobi is too slow, but very accurate

CES 703/5 Conclusion Imre P´ olik Transform A to tridiagonal form Outline Only eigenvalues: use QR Introduction Eigenvalues and eigenvectors: use D&C Algorithms Eigenvalues of a Small eigenvalues with high accuracy: use Jacobi tridiagonal matrix Conclusion Eigenvalues in [ α, β ) : use bisection Literature Most algorithms are available in LAPACK

Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der CES 703/5 Vorst, editors. Imre P´ olik Templates for the Solution of Algebraic Eigenvalue Problems: Outline A Practical Guide . Introduction SIAM, Philadelphia, 2000. Algorithms Eigenvalues of a James W. Demmel. tridiagonal matrix Applied Numerical Linear Algebra . Conclusion Literature SIAM, Philadelphia, 1997.