Finding Eigenvalues: Arnoldi Iteration and the QR Algorithm Will - PowerPoint PPT Presentation

Finding Eigenvalues: Arnoldi Iteration and the QR Algorithm Will Wheeler Algorithms Interest Group Jul 13, 2016

Outline Finding the largest eigenvalue ◮ Largest eigenvalue ← power method ◮ So much work for one eigenvector. What about others? ◮ More eigenvectors ← orthogonalize Krylov matrix ◮ build orthogonal list ← Arnoldi iteration Solving the eigenvalue problem ◮ Eigenvalues ← diagonal of triangular matrix ◮ Triangular matrix ← QR algorithm ◮ QR algorithm ← QR decomposition ◮ Better QR decomposition ← Hessenberg matrix ◮ Hessenberg matrix ← Arnoldi algorithm

Power Method How do we find the largest eigenvalue of an m × m matrix A ? ◮ Start with a vector b and make a power sequence: b , Ab , A 2 b , . . . ◮ Higher eigenvalues dominate: A n ( v 1 + v 2 ) = λ n 1 v 1 + λ n 2 v 2 ≈ λ n 1 v 1 Vector sequence (normalized) converged? ◮ Yes: Eigenvector ◮ No: Iterate some more

Krylov Matrix Power method throws away information along the way. ◮ Put the sequence into a matrix K = [ b , Ab , A 2 b , . . . , A n − 1 b ] (1) = [ x n , x n − 1 , x n − 2 , . . . , x 1 ] (2) ◮ Gram-Schmidt approximates first n eigenvectors v 1 = x 1 (3) v 2 = x 2 − x 2 · x 1 x 1 (4) x 1 · x 1 v 3 = x 3 − x 3 · x 1 x 1 − x 3 · x 2 x 2 (5) x 1 · x 1 x 2 · x 2 . . . (6) ◮ Why not orthogonalize as we build the list?

Arnoldi Iteration ◮ Start with a normalized vector q 1 x k = Aq k − 1 (next power) (7) k − 1 � y k = x k − ( q j · x k ) q j (Gram-Schmidt) (8) j =1 q k = y k / | y k | (normalize) (9) ◮ Orthonormal vectors { q 1 , . . . , q n } span Krylov subspace ( K n = span { q 1 , Aq 1 , . . . , A n − 1 q 1 } ) ◮ These are not the eigenvectors of A , but make a similarity (unitary) transformation H n = Q † n AQ n

Arnoldi Iteration ◮ Start with a normalized vector q 1 x k = Aq k − 1 (next power) (10) k − 1 � y k = x k − ( q j · x k ) q j (Gram-Schmidt) (11) j =1 q k = y k / | y k | (normalize) (12) ◮ Orthonormal vectors { q 1 , . . . , q n } span Krylov subspace ( K n = span { q 1 , Aq 1 , . . . , A n − 1 q 1 } ) ◮ These are not the eigenvectors of A , but make a similarity (unitary) transformation H n = Q † n AQ n

Arnoldi Iteration - Construct H n We can construct the matrix H n along the way. ◮ For k ≥ j : h j , k − 1 = q j · x k = q † j A q k − 1 ◮ For k = j : h j , k − 1 = | y k | (equivalent) ◮ For k < j : h j , k − 1 = 0 ◮ ⇒ H n is upper Hessenberg   h 1 , 1 h 1 , 2 h 1 , 3 h 1 , 4 h 1 , 5 h 1 , 6   h 2 , 1 h 2 , 2 h 2 , 3 h 2 , 4 h 2 , 5 h 2 , 6     0 h 3 , 2 h 3 , 3 h 3 , 4 h 3 , 5 h 3 , 6   H n = (13)   0 0 h 4 , 3 h 4 , 4 h 4 , 5 h 4 , 6     0 0 0 h 5 , 4 h 5 , 5 h 5 , 6 0 0 0 0 h 6 , 5 h 6 , 6 ◮ Why is this useful? H n = Q † n AQ n → same eigenvalues as A

Finding Eigenvalues How do we find eigenvalues of a matrix? ◮ Start with a triangular matrix A ◮ The diagonal elements are the eigenvalues   a 11 a 12 a 13 a 14 a 15 a 16   0 a 22 a 23 a 24 a 25 a 26     0 0 a 33 a 34 a 35 a 36   A = (14)   0 0 0 a 44 a 45 a 46     0 0 0 0 a 55 a 56 0 0 0 0 0 a 66 What if A isn’t triangular?

QR algorithm Factorize A = QR ◮ Q is orthonormal (sorry this is a different Q ) ◮ R is upper triangular (aka right triangular) Algorithm: 1. Start A 1 = A 2. Factorize A k = Q k R k 3. Construct A k +1 = R k Q k = Q − 1 k Q k R k Q k = Q − 1 k A k Q k ← similarity transform 4. A k converges to upper triangular

QR Algorithm Example (numpy.linalg.qr)   0 . 313 0 . 106 0 . 899   A 1 = 0 . 381 0 . 979 0 . 375 0 . 399 0 . 488 0 . 876   1 . 649 0 . 05 − 0 . 256   A 2 = 0 . 035 0 . 502 0 . 534 − 0 . 014 0 . 004 0 . 017   1 . 653 e + 00 2 . 290 e − 02 2 . 305 e − 01   A 3 = 5 . 966 e − 03 5 . 063 e − 01 − 5 . 342 e − 01 8 . 325 e − 05 − 9 . 429 e − 05 9 . 666 e − 03   1 . 653 e + 00 1 . 872 e − 02 − 2 . 285 e − 01   A 4 = 1 . 800 e − 03 5 . 063 e − 01 5 . 349 e − 01 − 4 . 812 e − 07 1 . 803 e − 06 9 . 554 e − 03

QR Decomposition - Gram-Schmidt How do we get the matrices Q and R ? http://www.seas.ucla.edu/ vandenbe/103/lectures/qr.pdf Recursively solve for the first column of Q and the first row of R : ◮ A = QR � � r 11 � ◮ � � � R 12 = a 1 A 2 q 1 Q 2 0 R 22 ◮ � � � � = q 1 R 12 + Q 2 R 22 a 1 A 2 r 11 q 1 ◮ r 11 = || a 1 || , q 1 = a 1 / r 11 ◮ Note q ⊺ 1 q 1 = 1 and q ⊺ 1 Q 2 = 0 ◮ R 12 = q ⊺ 1 A 2 ◮ Solve A 2 − q 1 R 12 = Q 2 R 22

QR Decomposition of Hessenberg Matrix What if we have a matrix in upper Hessenberg form?   h 1 , 1 h 1 , 2 h 1 , 3 h 1 , 4 h 1 , 5 h 1 , 6   h 2 , 1 h 2 , 2 h 2 , 3 h 2 , 4 h 2 , 5 h 2 , 6     0 h 3 , 2 h 3 , 3 h 3 , 4 h 3 , 5 h 3 , 6   H n = (15)   0 0 h 4 , 3 h 4 , 4 h 4 , 5 h 4 , 6     0 0 0 h 5 , 4 h 5 , 5 h 5 , 6 0 0 0 0 h 6 , 5 h 6 , 6 We just need to remove the numbers below the diagonal by combining rows → Givens rotation.

Givens Rotation A Givens rotation acts only on two rows and leaves the others unchanged.   γ − σ 0 0 0   σ γ 0 0 0     G 1 = 0 0 1 0 0 (16)     0 0 0 1 0 0 0 0 0 1 � a � � r � √ a 2 + b 2 Want G = , r = 0 b √ a 2 + b 2 ◮ γ = a / √ a 2 + b 2 ◮ σ = − b /

QR Decomposition of Hessenberg Matrix   � � � � � h 11 h 12 h 13 h 14 h 15   � � � � 0  h 22 h 23 h 24 h 25    G 2 G 1 H n = � � � (17)   0 0 h 33 h 34 h 35     0 0 h 43 h 44 h 45 0 0 0 h 54 h 55 So H n = G ⊺ 1 G ⊺ 2 . . . G ⊺ n − 1 R = QR

QR Decomposition of Hessenberg Matrix QR algorithm is iterative; is RQ also upper Hessenberg? Yes: Acting on the right, the Givens rotations mix two columns instead of rows, but change the same zeros.   � � � r 11 r 12 r 13 r 14 r 15   � r 21 r 22 � r 23 � r 24 r 25     G 2 G 1 H n = 0 � � (18) r 32 r 33 r 34 r 35     0 0 0 r 44 r 45 0 0 0 0 r 55

Householder Reflections ◮ Arnoldi finds a few eigenvalues of a large matrix ◮ In practice, QR uses Householder reflections ◮ www.math.usm.edu/lambers/mat610/sum10/lecture9.pdf Reflection across the plane perpendicular to a unit vector v : P = I − 2 vv ⊺ (19) Px = x − 2 vv ⊺ x (20) Px − x = − v (2 v ⊺ x ) (21) Want to arrange first column into zeros: find v so that Px = α e 1 .

Householder Reflections Want to arrange first column into zeros: find v so that Px = α e 1 . || x || = α (22) x − 2 vv ⊺ x = α e 1 (23) 1 2( x − α e 1 ) = v ( v ⊺ x ) (24) So we know ◮ v ∝ x − α e 1 ◮ v is a unit vector x − α e 1 v = (25) || x − α e 1 ||

Lanczos What if we apply the Arnoldi algorithm to a Hermitian matrix? ◮ Start with a normalized vector q 1 x k = Aq k − 1 (next power) (26) k − 1 � y k = x k − ( q j · x k ) q j (Gram-Schmidt) (27) j = k − 2 q k = y k / | y k | (normalize) (28) Since H n = Q † n AQ n , ◮ H n is also Hermitian ◮ Entries of H n above the first diagonal are 0 (tridiagonal) ◮ We don’t need to calculate them!

Lanczos What if we apply the Arnoldi algorithm to a Hermitian matrix? Since H n = Q † n AQ n , ◮ H n is also Hermitian ◮ Entries of H n above the first diagonal are 0 (tridiagonal) ◮ We don’t need to calculate them!   h ∗ 0 0 0 h 11 21  h ∗  h 21 h 22 0 0   32   h ∗ H n = 0 0 (29) h 32 h 33   43   h ∗ 0 0 h 43 h 44 54 0 0 0 h 54 h 55