Solving large scale eigenvalue problems Lecture 10, May 2, 2018: - PowerPoint PPT Presentation

Solving large scale eigenvalue problems Solving large scale eigenvalue problems Lecture 10, May 2, 2018: More on Lanczos and Arnoldi http://people.inf.ethz.ch/arbenz/ewp/ Peter Arbenz Computer Science Department, ETH Z¨ urich E-mail: arbenz@inf.ethz.ch Large scale eigenvalue problems, Lecture 10, May 2, 2018 1/25

Solving large scale eigenvalue problems Survey Survey of today’s lecture ◮ Restarting Arnoldi ◮ Thick restarts for Lanczos ◮ Krylov-Schur algorithm ◮ Rational Krylov space methods ◮ Polynomial filtering Large scale eigenvalue problems, Lecture 10, May 2, 2018 2/25

Solving large scale eigenvalue problems Survey When is a basis generating a Krylov space? Let v 1 , . . . , v k be linearly independent n -vectors. Is the subspace V := span { v 1 , . . . , v k } a Krylov space, i.e., is there a vector q ∈ V such that V = K k ( A , q )? Theorem V = span { v 1 , . . . , v k } is a Krylov space if and only if there is a k-by-k matrix M such that R := AV − VM , V = [ v 1 , . . . , v k ] , (1) has rank one and span { v 1 , . . . , v k , R ( R ) } has dimension k + 1 . For a proof see the lecture notes or paper by Stewart [1]. Large scale eigenvalue problems, Lecture 10, May 2, 2018 3/25

Solving large scale eigenvalue problems Survey The Lanczos algorithm with thick restarts Apply theorem to the case where a subspace is spanned by some Ritz vectors. For A = A ∗ a Lanczos relation is given by AQ k − Q k T k = β k +1 q k +1 e T k . (2) Let spectral decomposition of tridiagonal T k be S k = [ s ( k ) 1 , . . . , s ( k ) T k S k = S k Θ k , k ] , Θ k = diag( ϑ 1 , . . . , ϑ k ) . Then, for all i , the Ritz vector y i = Q k s ( k ) ∈ K k ( A , q ) i gives rise to the residual r i = A y i − y i ϑ i = β k +1 q k +1 s ( k ) ∈ K k +1 ( A , q ) ⊖ K k ( A , q ) . ki Large scale eigenvalue problems, Lecture 10, May 2, 2018 4/25

Solving large scale eigenvalue problems Survey The Lanczos algorithm with thick restarts (cont.) So, for any set of indices 1 ≤ i 1 < · · · < i j ≤ k we have A [ y i 1 , y i 2 , . . . , y i j ] − [ y i 1 , y i 2 , . . . , y i j ] diag( ϑ i 1 , . . . , ϑ i j ) = β k +1 q k +1 [ s ( k ) k , i 1 , s ( k ) k , i 2 , . . . , s ( k ) k , i j ] . Theorem = ⇒ any set [ y i 1 , y i 2 , . . . , y i j ] of Ritz vectors forms a Krylov space. The generating vector differs for each set { i 1 , · · · , i j } . Split the indices 1 , . . . , k in two sets. ◮ First set: ‘good’ Ritz vectors that we want to keep and that we collect in Y 1 ◮ Second set: ‘bad’ Ritz vectors that we want to remove. They go into Y 2 . Large scale eigenvalue problems, Lecture 10, May 2, 2018 5/25

Solving large scale eigenvalue problems Survey The Lanczos algorithm with thick restarts (cont.) � Θ 1 � = β k +1 q k +1 [ s ∗ 1 , s ∗ A [ Y 1 , Y 2 ] − [ Y 1 , Y 2 ] 2 ] . (3) Θ 2 Keeping the first set of Ritz vectors and purging (deflating) the rest yields AY 1 − Y 1 Θ 1 = β k +1 q k +1 s ∗ 1 . We now can restart a Lanczos procedure by orthogonalizing A q k +1 against Y 1 =: [ y ∗ 1 , . . . , y ∗ j ] and q k +1 . Large scale eigenvalue problems, Lecture 10, May 2, 2018 6/25

Solving large scale eigenvalue problems Survey The Lanczos algorithm with thick restarts (cont.) From the equation k s ( k ) σ i = β k +1 e ∗ A y i − y i ϑ i = q k +1 σ i , i we get q ∗ k +1 A y i = σ i , whence j � r k +1 = A q k +1 − β k q k +1 − σ i y i ⊥ K k +1 ( A , q . ) i =1 From this point on the Lanczos algorithm proceeds with the ordinary three-term recurrence. Large scale eigenvalue problems, Lecture 10, May 2, 2018 7/25

Solving large scale eigenvalue problems Survey The Lanczos algorithm with thick restarts (cont.) We finally arrive at relation AQ m − Q m T m = β m +1 q m +1 e T m with Q m = [ y 1 , . . . , y j , q k +1 , . . . , q m + k − j ] and   ϑ 1 σ 1 . ... .   .     ϑ j σ j   T m =   ...   σ 1 · · · σ j α k +1    ... ...    β m + k − j − 1   β m + k − j − 1 α m + k − j Large scale eigenvalue problems, Lecture 10, May 2, 2018 8/25

Solving large scale eigenvalue problems Survey The Lanczos algorithm with thick restarts (cont.) ◮ This procedure, called thick restart, has been suggested by Wu & Simon [2]. ◮ It allows to restart with any number of Ritz vectors. ◮ In contrast to the implicitly restarted Lanczos procedure, we need the spectral decomposition of T m . (Its computation is not an essential overhead.) ◮ The spectral decomposition admits a simple sorting of Ritz values. ◮ We could further split the first set of Ritz pairs into converged and unconverged ones, depending on the value β m +1 | s k , i | . If this quantity is below a given threshold we set the value to zero and lock (deflate) the corresponding Ritz vector, i.e., accept it as an eigenvector. Large scale eigenvalue problems, Lecture 10, May 2, 2018 9/25

Solving large scale eigenvalue problems Survey Thick restart Lanczos algorithm 1: Let us be given k Ritz vectors y i and a residual vector r k such that A y i = ϑ i y i + σ i r k +1 , i = 1 , . . . , k . The value k may be zero in which case r 0 is the initial guess. This algorithm computes an orthonormal basis y 1 , . . . , y k , q k +1 , . . . , q m that spans a m -dimensional Krylov space whose generating vector is not known unless k = 0. 2: z := A q k +1 ; 3: α k +1 := q ∗ k +1 z ; 4: r k +1 = z − α k +1 q k +1 − � k i =1 σ i y i 5: β k +1 := � r k +1 � 6: for i = k + 2 , . . . , m do q i := r i − 1 /β i − 1 . 7: z := A q i ; 8: Large scale eigenvalue problems, Lecture 10, May 2, 2018 10/25

Solving large scale eigenvalue problems Survey Thick restart Lanczos algorithm (cont.) α i := q ∗ i z ; 9: r i = z − α i q i − β i − 1 q i − 1 10: β i = � r i � 11: 12: end for ◮ Problem of losing orthogonality is similar to plain Lanczos. ◮ Wu–Simon [2] investigate various reorthogonalizing strategies known from plain Lanczos (full, selective, partial). ◮ In their numerical experiments the simplest procedure, full reorthogonalization, performs similarly or even faster than the more sophisticated reorthogonalization procedures. Large scale eigenvalue problems, Lecture 10, May 2, 2018 11/25

Solving large scale eigenvalue problems Krylov–Schur algorithm Krylov–Schur algorithm Krylov–Schur algorithm (Stewart [1]) is generalization of the thick-restart procedure for non-Hermitian problems. The Arnoldi algorithm constructs the Arnoldi relation AQ m = Q m H m + r m e ∗ m , where H m is Hessenberg and [ Q m , r m ] has full rank. Let H m = S m T m S ∗ m be a Schur decomposition of H m with unitary S m and upper triangular T m . Similarly as before we have s ∗ = e ∗ AY m = Y m T m + r m s ∗ , Y m = Q m S m , m S m . (4) Large scale eigenvalue problems, Lecture 10, May 2, 2018 12/25

Solving large scale eigenvalue problems Krylov–Schur algorithm Krylov–Schur algorithm (cont.) The upper triangular form of T m eases analysis of Ritz pairs. It admits moving unwanted Ritz values to the lower-right of T m . We collect ‘good’ and ‘bad’ Ritz vectors in matrices Y 1 and Y 2 : � T 11 � T 12 = β k +1 q k +1 [ s ∗ 1 , s ∗ A [ Y 1 , Y 2 ] − [ Y 1 , Y 2 ] 2 ] . (5) T 22 Keeping the first set of Ritz vectors and purging the rest yields AY 1 − Y 1 T 11 = β k +1 q k +1 s ∗ 1 . Large scale eigenvalue problems, Lecture 10, May 2, 2018 13/25

Solving large scale eigenvalue problems Krylov–Schur algorithm Krylov–Schur algorithm (cont.) Thick-restart Lanczos: eigenpair considered found if β k +1 | s ik | sufficiently small. Determination of converged subspace not so easy with general Krylov–Schur. However, if we manage to bring s 1 into the form � � � � s ′ 0 1 s 1 = = s ′′ s ′′ 1 1 then we found an invariant subspace: � T ′ T ′ � ∗ ] A [ Y ′ 1 , Y ′′ 1 ] − [ Y ′ 1 , Y ′′ 11 12 = β k +1 q k +1 [ 0 T , s ′′ 1 ] 1 T ′ 22 i.e., AY ′ 1 = Y ′ 1 T ′ 11 Large scale eigenvalue problems, Lecture 10, May 2, 2018 14/25

Solving large scale eigenvalue problems Krylov–Schur algorithm Krylov–Schur algorithm (cont.) ◮ In most cases s ′ 1 consists of a single small element or of two small elements in the case of a complex-conjugate eigenpair of a real nonsymmetric matrix [1]. ◮ These small elements are then declared zero and the columns in Y ′ 1 are locked, i.e., they are not altered anymore in the future computations. ◮ Orthogonality against them has to be enforced in the continuation of the eigenvalue computation though. Large scale eigenvalue problems, Lecture 10, May 2, 2018 15/25

Solving large scale eigenvalue problems Polynomial filtering Polynomial filtering ◮ SI-Lanczos/Arnoldi are the standard way to compute eigenvalues close to some target/shift σ . ◮ This works as long as the factorization of A − σ I (or of A − σ B ) is feasible. ◮ An alternative that has emerged in the last few years is polynomial filtering: Replace the matrix A by p ( A ) where p ∈ P d for some d . A u i = λ i u i = ⇒ p ( A ) u i = p ( λ i ) u i . Eigenvectors remain unchanged (at least if λ i � = λ j implies p ( λ i ) � = p ( λ j )). Large scale eigenvalue problems, Lecture 10, May 2, 2018 16/25