A randomized FEAST algorithm for generalized eigenvalue problems
aa r X i v : . [ m a t h . NA ] D ec A RANDOMIZED FEAST ALGORITHM FOR GENERALIZEDEIGENVALUE PROBLEMS
GUOJIAN YIN ∗ Abstract.
The FEAST algorithm, due to Polizzi, is a typical contour-integral based eigensolverfor computing the eigenvalues, along with their eigenvectors, inside a given region in the complexplane. It was formulated under the circumstance that the considered eigenproblem is Hermitian. TheFEAST algorithm is stable and accurate, and has attracted much attention in recent years. However,it was observed that the FEAST algorithm may fail to find the target eigenpairs when applying itto the non-Hermitian problems. Efforts have been made to adapt the FEAST algorithm to non-Hermitian cases. In this work, we develop a new non-Hermitian scheme for the FEAST algorithm.The mathematical framework will be established, and the convergence analysis of our new methodwill be studied. Numerical experiments are reported to demonstrate the effectiveness of our methodand to validate the convergence properties.
Key words. generalized eigenvalue problems, contour integral, spectral projection
AMS subject classifications.
1. Introduction.
Large-scale non-Hermitian eigenvalue problems arise in vari-ous areas of science and engineering, such as dynamic analysis of structures [14], linearstability analysis of the Navier-Stokes equation in fluid dynamics [6], the electron en-ergy and position problems in quantum chemistry [10], and resonant state calculation[22]. In most practical applications, it is not the whole spectrum but rather a signif-icant part of it is of interest to the users [20]. For example, in the model reductionof a linear dynamical system, one only needs to know the response over a range offrequencies, see [3, 11].Consider the generalized eigenvalue problem A x = λB x , (1.1)where A, B ∈ C n × n . The scalars λ ∈ C and the associated vectors x ∈ C n , x = 0, arecalled the eigenvalues and their associated (right) eigenvectors, respectively [3, 8, 13].When B is the identity matrix, then (1.1) becomes a standard eigenvalue problem.In this work, our objective is to compute the eigenvalues of (1.1) inside a given regionin the complex plane, along with the corresponding eigenvectors.Computing the partial spectrum of a large-scale problem is very difficult in prac-tice. Maybe the most straightforward method is first using the well-known QZ method[19] to compute the whole spectrum and then selecting the target eigenvalues. This di-rect method costs about O ( n ) [13], consequently, it is prohibitively expensive whenthe size of considered problem is large. In the past decades, the most successfulmethods for solving the partial spectrum of a large eigenproblem are based on theprojection techniques [3, 4, 29], of which perhaps the Krylov subspace approaches arethe most wildly used [23, 24]. However, the existing projection methods mainly focuson computing the extreme eigenvalues [28] or the eigenvalues close to a given shift[14].Recently, a class of eigensolvers based on contour integrals were proposed forcomputing the eigenvalues inside a given region in the complex plane [2, 5, 9, 21,25, 26, 27]. Unlike the well-known Krylov subspace methods, these new methods use ∗ School of Mathematics, Sun Yat-sen University, Guangzhou, P. R. China( [email protected] ). 1 pecifically defined contour integrals to generate subspaces to contain the eigenspacecorresponding to the target eigenvalues. Then the projection techniques are used toextract the target eigenpairs. Two typical examples of these contour-integral basedeigensolvers are the Sakurai-Sugiura (SS) method [26] and the FEAST algorithmdeveloped by Polizzi in [21]. By noticing that the SS method always suffers fromnumerical instability [2, 17], Sakurai et al. turned to use the Rayleigh-Ritz procedureto extract the target eigenpairs, and leaded to a more stable contour-integral basedeigensolver, called CIRR [16, 27].The derivation of both CIRR and FEAST is under the assumptions that A and B are Hermitian matrices and B is positive definite, i.e., (1.1) is a Hermitian problem.It was shown in [32] that CIRR and FEAST may fail to find the target eigenpairswhen (1.1) is a non-Hermitian problem. Motivated by this observation, the authorsin [32] developed a non-Hermitian FEAST algorithm to make the FEAST algorithmalso applicable for the non-Hermitian problems. Instead of the orthogonal projectiontechnique used in the FEAST algorithm, the non-Hermitian FEAST algorithm pro-posed in [32] uses to the oblique projection technique with appropriately chosen leftsubspace to extract desired eigenpairs.In this work, we would like to formulate another non-Hermitian scheme for theFEAST algorithm. We find that the FEAST algorithm can deal with the non-Hermitian problems if the left subspace spanned by a random matrix. The theoreticalanalysis will be given to justify our findings. The convergence properties also will bestudied to show the effectiveness of our method.The paper is organised as follows. In Section 2, we briefly describe the FEASTalgorithm [21]. In Section 3, we review the non-Hermitian variant of the FEASTalgorithm proposed in [32]. We formulate our new non-Hermitian FEAST algorithmand give convergence analysis in Section 4. In Section 5, numerical experiments arereported to illustrate the numerical performance of our method.Throughout the paper, we use the following notation and terminology. The sub-space spanned by the columns of a matrix X is denoted by span { X } . The rankand conjugate transpose of X are denoted by rank( X ) and X ∗ respectively. Thealgorithms are presented in Matlab style.
2. Introduction to FEAST.
In this section, we provide a brief review of theFEAST algorithm [21]. The algorithm was formulated under the assumptions that A and B are Hermitian and B is positive definite, in which case the eigenvalues of (1.1)are real-valued [8]. The FEAST algorithm was developed for finding all eigenvaluesof (1.1) within a specified interval, say [ σ , σ ], and their associated eigenvectors.Without loss of generality, assume that the eigenvalues inside [ σ , σ ] are λ ≤ λ ≤ . . . ≤ λ s . Therefore, there are s eigenvalues inside [ σ , σ ].Essentially, the FEAST algorithm belongs to the family of subspace iteration withorthogonal projection [30]. Unlike the better known Krylov subspace methods, theFEAST algorithm constructs a subspace that envelops the desired eigenspace via thecontour integral defined as V := 12 π √− I Γ ( zB − A ) − dzY, (2.1)where Γ is any contour that contains { λ i } si =1 inside, and Y is an n × s random matrix.To formulate the FEAST algorithm, we need the following theorem. Theorem 2.1 ( [29] ). Let A and B be n × n Hermitian matrices and that B is ositive definite. Then there exists an n × n matrix X = [ x , x , . . . , x n ] for which X ∗ BX = I n and X ∗ AX = Λ = diag( λ , λ · · · , λ n ) , (2.2) where I n is the n × n identity matrix, { λ i } ni =1 are the eigenvalues of the matrix pencil zB − A , and the columns { x i } ni =1 of X are their associated eigenvectors. By (2.2) and the residue theorem in complex analysis [1], we have V = 12 π √− X I Γ ( zI n − Λ) − dzY ( BX ) − = X (: , s ) ( X (: , s ) ) ∗ Y. (2.3)Then the columns of V form a basis for the eigenspace span { X (: , s ) } , if ( X (: , s ) ) ∗ Y is full-rank. Forming the s × s matrices ˆ A = V ∗ AV and ˆ B = V ∗ BV , solving theproblem (1.1) now is reduced to computing the eigenpairs of the projected eigenvalueproblem ˆ A y = λ ˆ B y , (2.4)according to the Rayleigh-Ritz procedure [21, 29].To generate the projected eigenproblem (2.4), the most important task is to com-pute the basis vectors V . In view of (2.1) and (2.3), we know that V has to becomputed numerically by a quadrature scheme. Let Γ be the circle with center at γ = ( σ + σ ) / ρ = ( σ − σ ) /
2, applying the q -point Gauss-Legendrequadrature [7] to compute V numerically yields V = 12 π √− I Γ ( zB − A ) − dzY ≈ q X j =1 ω j ( z j − γ )( z j B − A ) − Y, (2.5)where z j = γ + ρe √− θ j , θ j = (1 + t j ) π , and t j is the j th Gaussian node withassociated weight ω j . From (2.5), one can see that the dominant computational workof the FEAST algorithm is solving the linear systems of the form( z j B − A ) X j = Y, j = 1 , . . . , q. (2.6)The complete FEAST algorithm is given as follows.
Algorithm 1.
Input Hermitian matrices A and B with B being positive definite,a uniformly-distributed random matrix Y ∈ R n × t , where t ≥ s , the circle Γ enclosingthe interval [ σ , σ ] , and a convergence tolerance ǫ . The function “ Feast ” computeseigenpairs (ˆ λ i , ˆ x i ) of (1.1) that satisfy ˆ λ i ∈ [ σ , σ ] and s X i =1 ˆ λ i < ǫ, (2.7) and they are output in the vector Λ s and the matrix X s . Function [Λ s , X s ] = Feast ( A, B, Y, Γ , ǫ )1. Compute V approximately by (2.5).2. Set ˆ A = V ∗ AV and ˆ B = V ∗ BV .3. Solve the generalized eigenproblem of size t : ˆ A y = ˆ λ ˆ B y , to obtain theeigenpairs { (ˆ λ i , y i ) } ti =1 .4. Compute ˆ x i = V y i , i = 1 , , . . . t .5. Check if { (ˆ λ i , ˆ x i ) } ti =1 satisfy the convergence criteria (2.7). If s eigenpairssatisfy (2.7), stop. Otherwise, set X t = [ˆ x , ˆ x , . . . , ˆ x t ] and Y = BX t , thengo back to Step 1. he FEAST algorithm is an accurate and reliable technique [18, 30]. It transformsthe difficulty of solving the eigenproblem (1.1) to that of solving linear systems (2.6).Since the quadrature nodes z j and the columns of the right-hand sides in (2.6) areindependent, the FEAST algorithm can be easily implemented on parallel machines.Due to these appealing features, the FEAST algorithm attracts much attention re-cently.
3. A non-Hermitian FEAST algorithm.
The FEAST algorithm was for-mulated when (1.1) is a Hermitian problem. However, when it comes to the non-Hermitian problem, it was found in [32] that the FEAST algorithm may fail tocompute the desired eigenpairs; a simple example was given to illustrate this fact.Motivated by this observation, the authors in [32] developed a non-Hermitian FEASTalgorithm so as to adapt FEAST to the non-Hermitian cases. The key to the successof their non-Hermitian FEAST algorithm is that the oblique projection technique,instead of the orthogonal projection technique used in FEAST, with appropriatelychosen left subspace is used to extract the desired eigenpairs.The only requirement for the non-Hermitian FEAST algorithm proposed in [32] isthat the matrix pencil zB − A is regular, which means the method is able to deal withthe most common generalized eigenproblems [3]. Recall that a matrix pencil zB − A is regular if det( zB − A ) is not identically zero for all z ∈ C . As with the Jordancanonical form for a matrix, there exists a canonical form for the regular matrix pencil zB − A . Theorem 3.1 (
The Weierstrass canonical form [12, 32] ). Let zB − A be a regularmatrix pencil of order n . Then there exist nonsingular matrices S and T ∈ C n × n suchthat T AS = (cid:20) J d I n − d (cid:21) and T BS = (cid:20) I d N n − d (cid:21) , (3.1) where J d is a d × d matrix in Jordan canonical form with its diagonal entries corre-sponding to the eigenvalues of zB − A , N n − d is an ( n − d ) × ( n − d ) nilpotent matrixalso in Jordan canonical form, and I d denotes the identity matrix of order d . Let J d be of the form J d = J d ( λ ) 0 · · · J d ( λ ) · · · · · · J d m ( λ m ) (3.2)where P mi =1 d i = d and J d i ( λ i ) are d i × d i matrices of the form J d i ( λ i ) = λ i · · · λ i · · · λ i , i = 1 , , . . . , m with λ i being the eigenvalues. Here the λ i are not necessarily distinct and can berepeated according to their multiplicities. et N n − d be of the form N n − d = N d ′ · · · N d ′ · · · · · · N d ′ m ′ , where P m ′ i =1 d ′ i = n − d and N d ′ i are d ′ i × d ′ i matrices of the form N d ′ i = · · ·
00 0 1 .... . . . . . . . . 0... . . . . . . 10 · · · , i = 1 , , . . . , m ′ . Partition S into block form S = [ S , S , . . . , S m , S m +1 ], where each S i ∈ C n × d i ,1 ≤ i ≤ m , and S i into S i = [ s i , s i , . . . , s id i ] with s ij ∈ C n , 1 j d i . It was verifiedin [32] that for any eigenvalue λ i , 1 ≤ i ≤ m ,( λ i B − A ) S (cid:20) I d N n − d (cid:21) = BS (cid:20) λ i I d − J d λ i N n − d − I n − d (cid:21) . (3.3)By comparing the first d columns on both sides above, we get( λ i B − A ) s ij = B s ij − , ≤ j ≤ d i , ≤ i ≤ m, (3.4)with s i ≡ . We can see that s i are the eigenvectors corresponding to the eigenvalues λ i for all 1 ≤ i ≤ m .Let Γ be a positively oriented simple closed curve enclosing the desired eigenval-ues. Again without loss of generality, we let the eigenvalues of (1.1) enclosed by Γ be { λ , . . . , λ l } , and s := d + d + · · · + d l be the number of eigenvalues inside Γ withmultiplicity taken into account. Define the contour integral Q := 12 π √− I Γ ( zB − A ) − Bdz. (3.5)According to the residue theorem in complex analysis [1], it was verified in [32] that Q = 12 π √− I Γ ( zB − A ) − Bdz = S (cid:20) I s
00 0 (cid:21) S − = S (: , s ) ( S − ) (1: s, :) . (3.6)One can show that Q = Q , which means Q is a projector onto subspace K =span { S (: , s ) } . Define U := QY = S (: , s ) ( S − ) (1: s, :) Y, (3.7)where Y is an n × s random matrix. Therefore, U is the projection of Y onto thesubspace K . Now we would like to show that the columns of U form a basis for thesubspace K . We begin with emma 3.2 ( [32] ). Let Y ∈ R n × s . If the entries of Y are random numbers froma continuous distribution and that they are independent and identically distributed(i.i.d.), then with probability 1, the matrix ( S − ) (1: s, :) Y is nonsingular. According to (3.7) and Lemma 3.2, we can conclude that the columns of U forma basis for the subspace K . Note that K contains the eigenspace corresponding to thedesired eigenvalues (see (3.4) for details), it is natural to take K as the right subspace.The FEAST algorithm takes advantage of the often used orthogonal projection tech-nique to extract desired eigenpairs. The authors in [32] found that this extractionapproach may fail to compute the desired eigenpairs for the FEAST algorithm when(1.1) is a non-Hermitian problem. To address this deficiency, they resorted to theoblique projection method, and developed a non-Hermitian FEAST algorithm. Intheir method, the left subspace is taken as B K ; the approximate eigenpairs (¯ λ, ¯ x ) areobtained by imposing the Petrov-Galerkin condition [3, 24]:( A ¯ x − ¯ λB ¯ x ) ⊥ B K , (3.8)where ¯ λ ∈ C and ¯ x ∈ K . It was shown in [32] that the columns of BU form a basisfor B K . Therefore (3.8) can be written in matrix form( BU ) ∗ ( AU y − ¯ λBU y ) = 0 , (3.9)where y ∈ C s satisfying ¯ x = U y . Accordingly, solving the eigenvalues of (1.1) insideΓ now is reduced to solve the projected eigenproblem A y = ¯ λB y , (3.10)with A = ( BU ) ∗ AU and B = ( BU ) ∗ BU. (3.11)The key to the success of the non-Hermitian FEAST algorithm proposed in [32] isthat the left subspace is taken as B K , instead of K used in the FEAST algorithm.Due to this, below we call this non-Hermitian FEAST algorithm BFEAST for theease of reference. The following theorem justifies their choice of the left subspace. Theorem 3.3.
Let { (¯ λ i , y i ) } si =1 be the eigenpairs of the projected eigenproblem(3.10). Then { (¯ λ i , U y i ) } si =1 are the eigenpairs of (1.1) located inside Γ . In order to generate the projected eigenproblem (3.10), the most important taskis to compute the projection U (see (3.7)). In practice, we have to compute U ap-proximately by a quadrature rule: U = QY = 12 π √− I Γ ( zB − A ) − BdzY ≈ π √− q X j =1 ω j ( z j B − A ) − BY, (3.12)where z j are the quadrature nodes on Γ associated with weights ω j . From (3.12), weknow that the dominant work is solving q linear systems of the form( z i B − A ) X i = BY. (3.13)The non-Hermitian FEAST algorithm (BFEAST) can be described as follows.
Algorithm 2.
Input
A, B ∈ C n × n , an i.i.d. random matrix Y ∈ R n × t where t ≥ s , a closed curve Γ , a convergence tolerance ǫ , and “ max iter ” to control the maximum umber of iterations. The function “ BFEAST ” computes eigenpairs (¯ λ i , ¯ x i ) of (1.1)that satisfies ¯ λ i inside Γ and k A ¯ x i − ¯ λ i B ¯ x i k k A ¯ x i k + k B ¯ x i k < ǫ. (3.14) The results are stored in the vector Λ s and the matrix X s . Function [Λ s , X s ] = BFEAST ( A, B, Y, Γ , ǫ, max iter )1. For k = 1 , · · · , max iter
2. Compute U approximately by the quadrature rule (3.12).3. Compute QR decompositions: U = U R and BU = U R .
4. Form A = U ∗ AU and B = U ∗ BU .5. Solve the projected eigenproblem A y = ¯ λB y of size t to obtain eigenpairs { (¯ λ i , y i ) } ti =1 . Set ¯ x i = U y i , i = 1 , , . . . , t .6. Set Λ s = [ ] and X s = [ ].7. For i = 1 : t
8. If (¯ λ i , ¯ x i ) satisfies (3.14), then Λ s = [Λ s , ¯ λ i ] and X s = [ X s , ¯ x i ].9. End10. If there are s eigenpairs satisfying (3.14), stop. Otherwise, set Y = U .11. End.
4. A randomized FEAST algorithms.
In the previous two sections, we re-viewed the FEAST algorithm, as well as its non-Hermitian variation, i.e., the BFEASTalgorithm. In this part, we first formulate another non-Hermitian scheme for theFEAST algorithm. After that, the convergence analysis will be given to illustrate theeffectiveness of our new method.
In [32], the authors used the obliqueprojection technique, rather than the wildly used orthogonal projection technique, toextend FEAST to the non-Hermitian problems. The key step is that they take the leftsubspace to B K instead of K used in the original FEAST algorithm. Here we presentanother scheme for the non-Hermitian FEAST algorithm. The intuition behind ournew method is inspired by Lemma 3.2. The following theorem validates our intuition. Theorem 4.1.
Let R be an n × s random matrix, whose entries are independentand identically distributed (i.i.d.). Define e A = R ∗ AU and e B = R ∗ BU. (4.1)
Let { (˜ λ i , y i ) } si =1 be the eigenpairs of the projected eigenproblem e A y = ˜ λ e B y . (4.2) Then { (˜ λ i , U y i ) } si =1 are the eigenpairs of (1.1) located inside Γ .Proof . By (3.1), one can verify that(˜ λB − A ) S (: , s ) = ( T − ) (: , s ) (˜ λI (1: s, s ) − J (1: s, s ) ) . (4.3)By (3.7), (4.1) and (4.3), we have˜ λ e B − e A = R ∗ (˜ λB − A ) U (4.4)= R ∗ (˜ λB − A ) S (: , s ) ( S − ) (1: s, :) Y = R ∗ ( T − ) (: , s ) (˜ λI (1: s, s ) − J (1: s, s ) )( S − ) (1: s, :) Y. herefore, the characteristic polynomial of the eigenproblem (4.2) isdet(˜ λ e B − e A ) = det( R ∗ ( T − ) (: , s ) ) det(˜ λI (1: s, s ) − J (1: s, s ) ) det(( S − ) (1: s, :) Y ) . By Lemma 3.2, we know that R ∗ ( T − ) (: , s ) and ( S − ) (1: s, :) Y are nonsingular. As aresult, due to the special structure of J (1: s, s ) , the roots of the characteristic polyno-mial det(˜ λ e B − e A ) are λ , λ , . . . , λ l with multiplicities d , d , . . . , d l respectively.Recall that λ , . . . , λ l are not necessary distinct. Without loss of generality, let usconsider the case where ˜ λ = λ = λ and ˜ λ = λ i for 3 i l . Since (˜ λ e B − e A ) y = 0and R ∗ ( T − ) (: , s ) is nonsingular, by (4.4) we have(˜ λI (1: s, s ) − J (1: s, s ) )( S − ) (1: s, :) Y y = 0 . The above equation in turn implies that ( S − ) (1: s, :) Y y = α e + β e d +1 for somescalars α and β not both zero due to the special structure of J , where e and e d +1 are the first and the ( d + 1)th columns of the s × s identity matrix, respectively.Therefore, U y = S (: , s ) ( S − ) (1: s, :) Y y = αS (: , s ) e + βS (: , s ) e d +1 . Note that the first and the ( d + 1)th columns of S are the eigenvector of (1.1)corresponding to the eigenvalues λ and λ respectively by (3.4). Thus, their linearcombinations are the eigenvectors associated with ˜ λ . The proof is completed.Theorem 4.1 tells us that the FEAST algorithm can deal with the non-Hermitianeigenproblems if we take the left subspace spanned by a random matrix. Due to theusage of random matrix, we call our new non-Hermitian FEAST algorithm RFEASTfor ease of reference. Algorithm 3.
Input
A, B ∈ C n × n , an i.i.d. random matrix Y ∈ R n × t where t ≥ s , a closed curve Γ , a convergence tolerance ǫ , and “ max iter ” to control the maximumnumber of iterations. The function “ RFEAST ” computes eigenpairs (˜ λ i , ˜ x i ) of (1.1)that satisfies ˜ λ i inside Γ and k A ˜ x i − ˜ λ i B ˜ x i k k A ˜ x i k + k B ˜ x i k < ǫ. (4.5) The results are stored in the vector Λ s and the matrix X s . Function [Λ s , X s ] = RFEAST ( A, B, Y, Γ , ǫ, max iter )1. For k = 1 , · · · , max iter
2. Compute U approximately by the quadrature rule (3.12).3. Generate an n × t random matrix R , and compute QR decompositions: U = U R and R = U R .
4. Form ˜ A = U ∗ AU and ˜ B = U ∗ BU .5. Solve the projected eigenproblem ˜ A y = ˜ λ ˜ B y of size t to obtain eigenpairs { (˜ λ i , y i ) } ti =1 . Set ˜ x i = U y i , i = 1 , , . . . , t .6. Set Λ s = [ ] and X s = [ ].7. For i = 1 : t
8. If (˜ λ i , ˜ x i ) satisfies (4.5), then Λ s = [Λ s , ˜ λ i ] and X s = [ X s , ˜ x i ].9. End10. If there are s eigenpairs satisfying (4.5), stop. Otherwise, set Y = U .11. End. ur method can make the FEAST algorithm applicable to the non-Hermitianproblems. Obviously, in each iteration, the dominant work in our method is to com-pute the projection U by a quadrature scheme, see (3.12) for details. As with othercontour-integral based methods, our algorithm replaces the difficulty of solving theeigenvalue problem (1.1) by the difficulty of solving the linear systems (3.13), and hasa good potential to be parallelized. In this part, we study the convergence propertiesof our new method (Algorithm 3) to show its effectiveness.An important quantity for the convergence properties of projection methods isthe distance of the exact eigenvector from the search subspace [24]. We begin theconvergence analysis of our method from this perspective. For notational convenience,we represent the approximate projection computed in the k th iteration in Algorithm3 by U ( k ) .Compute the contour integral f ( µ ) = 12 π √− I Γ z − µ dz (4.6)by a q -point quadrature rule on Γ: f ( µ ) ≈ ˜ f ( µ ) = 12 π √− q X i =1 ω i z i − µ , (4.7)where z j are the quadrature nodes on Γ associated with weights ω j . Let Γ be anunit circle with center at γ and µ = γ + re √− θ , θ ∈ ( − π, π ]. Therefore, µ is locatedinside Γ when 0 ≤ r < r >
1. In theory, accordingto the residue theorem, we have that f ( µ ) = 1 if µ is located inside Γ and f ( µ ) = 0if µ is located outside Γ [1]. Compute the approximation ˜ f ( µ ) of f ( µ ) by the Gauss-Legendre quadrature with 16 integration points on Γ. Fig | ˜ f ( µ ) | . We can see that | ˜ f ( µ ) | is close 1 when µ is contained inside Γ and is closeto 0 when µ is outside Γ. Without loss of generality, assume that | ˜ f ( λ ) | ≥ | ˜ f ( λ ) | ≥ · · · ≥ | ˜ f ( λ l ) | > | ˜ f ( λ l +1 ) | ≥ · · · ≥ | ˜ f ( λ m ) | . (4.8) Theorem 4.2.
Let d = 0 , then S (: , d + ... + d j − ) is an eigenvector correspondingto λ j . Let t be the size of starting vectors Y satisfying t > s , then there exists aninteger l ′ , l ′ > l , such that P l ′ − i d i < t ≤ P l ′ i d i . Suppose the eigenvalues outside Γ are simple, which implies that d i = 1 for i = l + 1 , . . . , m . There exists a vector v ( k ) j ∈ span { U ( k ) } , j = 1 , . . . , l , such that k S (: , d + ... + d j − ) − v ( k ) j k ≤ τ j (cid:0) | ˜ f ( λ l ′ ) || ˜ f ( λ j ) | (cid:1) k , (4.9) where τ j is a constant. In particular, k ( I n − Q ( k ) ) S (: , d + ... + d j − ) k ≤ τ j (cid:0) | ˜ f ( λ l ′ ) || ˜ f ( λ j ) | (cid:1) k , (4.10) where Q ( k ) is the orthogonal projector onto the subspace span { U ( k ) } . | ˜ f ( µ ) | -0.500.511.5 θ = 0 θ = π /6 θ = π r l og | ˜ f ( µ ) | -10-8-6-4-202 θ = 0 θ = π /6 θ = π Fig. 4.1 . Here Γ is a unit circle with center at γ . The approximation ˜ f ( µ ) , where µ = γ + re √− θ , θ ∈ ( − π, π ] , is computed by the Gauss-Legendre quadrature with quadrature nodes.Therefore, µ is located inside Γ when ≤ r < and is located outside Γ when r > . The leftpicture shows the general shape of | ˜ f ( µ ) | , while the right one shows the logarithmic scale shape ofthe function. Proof . Let U (0) = Y be an n × t random matrix and Z = ( S − ) (1: t, :) U (0) . ByLemma 3.2, we know Z is nonsingular, then U (0) = SS − U (0) = (cid:2) S (: , t ) ( S − ) (1: t, :) + S (: ,t +1: n ) ( S − ) ( t +1: n, :) (cid:3) U (0) = h S (: , t ) + S (: ,t +1: n ) ( S − ) ( t +1: n, :) U (0) Z − i Z = (cid:2) ( S (: , s ) + S (: ,t +1: n ) E (0) ) , V (0) (cid:3) Z, (4.11)where E (0) is the first s columns of ( S − ) ( t +1: n, :) U (0) Z − , and V (0) is the last ( t − s )columns of matrix S (: , t ) + S (: ,t +1: n ) ( S − ) ( t +1: n, :) U (0) Z − .Under the assumption that the eigenvalues outside Γ are simple, the matrix N n − d in (3.1) is a zero matrix. Let D = 12 π √− q X j =1 ω j (cid:20) ( z j I d − J d ) −
00 0 (cid:21) . (4.12)It was shown in [15, 31] that | D ( i,i ) | >
0, for i = 1 , . . . , s . According to (2.5), (3.12),(4.11) and (4.12), we have U (1) = 12 π √− q X j =1 ω j ( z j B − A ) − BU (0) = SDS − U (0) = (cid:2) ( S (: , s ) + S (: ,t +1: n ) E (1) ) D (1: s, s ) , V (1) (cid:3) Z, (4.13)where E (1) = D ( t +1: n,t +1: n ) E (0) ( D (1: s, s ) ) − and V (1) = SDS − V (0) . Denote theQR decomposition of U ( k ) by U ( k ) = U k R k . By induction, we have the followingrelationship U ( k ) = SDS − U k − = SDS − U ( k − ( R k − ) − = (cid:2) ( S (: , s ) + S (: ,t +1: n ) E ( k ) )( D (1: s, s ) ) k , V ( k ) (cid:3) Z ( R k − . . . R ) − , (4.14) here E ( k ) = ( D ( t +1: n,t +1: n ) ) k E (0) (( D (1: s, s ) ) − ) k and V ( k ) = ( SDS − ) k V (0) .Since Z ( R k − . . . R ) − is nonsingular, from (4.14) we conclude that the columnsof S (: , s ) + S (: ,t +1: n ) E ( k ) are in the subspace span { U ( k ) } . In particular, vectors v ( k ) j = S (: , d + ... + d j − ) + S (: ,t +1: n ) ( E ( k ) ) (: , d + ... + d j − ) ∈ span { U ( k ) } for j = 1 , . . . , l . Bythe special structure of ( D (1: s, s ) ) − ) k , we have k S (: , d + ... + d j − ) − v ( k ) j k = k S (: ,t +1: n ) ( E ( k ) ) (: , d + ... + d j − ) k = (cid:0) | ˜ f ( λ j ) | (cid:1) k k S (: ,t +1: n ) ( D ( t +1: n,t +1: n ) ) k ( E (0) ) (: , d + ... + d j − ) k ≤ τ j (cid:0) | ˜ f ( λ l ′ ) || ˜ f ( λ j ) | (cid:1) k , (4.15)where τ j = k S (: ,t +1: n ) k k ( E (0) ) (: , d + ... + d j − ) k .Moreover, k ( I n − Q ( k ) ) S (: , d + ... + d j − ) k = min v ∈ span { U ( k ) } k S (: , d + ... + d j − ) − v k ≤ τ j (cid:0) | ˜ f ( λ l ′ ) || ˜ f ( λ j ) | (cid:1) k . (4.16)Note that l ′ > l , which means λ l ′ is located outside Γ. Suppose the | ˜ f ( λ l ′ ) | isabout 1 . × − , we can expect that there exists a vector v ( k ) j in span { U ( k ) } suchthat k S (: , d + ... + d j − ) − v ( k ) j k → − k .Let P ( k ) be the oblique projector onto the subspace span { U ( k ) } and orthogonal tothe left subspace generated by a random matrix in the k th iteration in Algorithm 3.Define approximate operators A k = P ( k ) AQ ( k ) and B k = P ( k ) BQ ( k ) . The followingtheorem gives an upper bound for the residual norm of the exact eigenpair with respectto the approximate operator pair ( A k , B k ). Theorem 4.3.
Let σ ( k ) = k P ( k ) ( A − λB )( I n − Q ( k ) ) k . Then the residual norms of ( λ j , S (: , d + ... + d j − ) ) , j = 1 , . . . , l , for the approximate operator pair ( A k , B k ) saftisfy k ( A k − λ j B k ) S (: , (: , d + ... + d j − )) k ≤ σ ( k ) τ j (cid:0) | ˜ f ( λ l ′ ) || ˜ f ( λ j ) | (cid:1) k (4.17) Proof . Similar to the proof of Lemma 2 in [15] and Theorem 4.7 in [24], we have k ( A k − λ j B k ) S (: , d + ... + d j − ) k = k P ( k ) ( A − λ j B )( I n − Q ( k ) ) S (: , d + ... + d j − ) k = k P ( k ) ( A − λ j B )( I n − Q ( k ) )( I n − Q ( k ) ) S (: , d + ... + d j − ) k ≤ σ ( k ) k ( I n − Q ( k ) ) S (: , d + ... + d j − ) k . By (4.10), we establish our result (4.17).Theorem 4.3 says that the residual norm associated with the exact eigenpair( λ j , S (: , d + ... + d j − ) ) converges at a rate of (cid:0) | ˜ f ( λ l ′ ) | / | ˜ f ( λ j ) | (cid:1) with respect to theiteration counts. able 5.1 Test problems from Matrix Market that are used in our experiments.
No. Problem Type n Region: ( γ, ρ ) s − . × , . × ) 2302 DWG961 gen. 961 (5 . × , . × ) 1573 UTM1700 gen. 1700 (4 . , .
0) 964 MHD4800 gen. 4800 ( − . , .
0) 1695 OLM5000 stand. 5000 ( − . × , . × ) 2046 DW8192 stand. 8192 (1 . , .
2) 270
5. Numerical Experiments.
In this section, we present some numerical exper-iments to demonstrate the numerical performance of our new non-Hermitian FEASTalgorithm (RFEAST). The experiments are organized into threes sets. The first setaims at demonstrating the convergence behavior of our new method. The secondset is devoted to comparing our technique with another non-Hermitian variant ofFEAST, that is the BFEAST algorithm (Algorithm 2). In the last set, we would liketo compare our RFEAST method with the
Matlab built-in function eig . For theapproximation eigenpairs (˜ λ i , ˜ x i ), define the relative residual norms r i = k A ˜ x i − ˜ λ i B ˜ x i k k A ˜ x i k + k B ˜ x i k . (5.1)We use the maximum relative residual norm defined as Res = max ≤ i ≤ s r i to assessthe accuracy achieved by the test methods. All computations are carried out in Matlab version R2014b on a MacBook with an Intel Core i5 2.5 GHz processor and8 GB RAM.The test matrices presented in
Table . They are the real-world problems from scientific and engineering ap-plications. All test problems are non-Hermitian. The first four test problems aregeneralized eigenvalue problems and the last two test problems are standard. Theregion of interest for each test problem is a circle with center at γ and radius ρ . Thevalue of s is the number of eigenvalues inside the target region. In all experiments,we use the Gauss-Legendre quadrature with 16 quadrature nodes to compute theapproximate projection U (see (3.12)). The generalized shifted linear systems (see(3.13)) involved are computed by direct method. We first use the Matlab function lu to compute the LU decomposition of A − z j B, j = 1 , , . . . , q , and then performthe triangular substitutions to get the corresponding solutions. The FEAST algorithm is a stable and fasttechnique [18, 30]. It was formulated for the Hermitian problems [21]. The goal ofour work is to adapt FEAST to the non-Hermitian cases. Meanwhile, we hope thatour method retains the effectiveness of the FEAST algorithm. The objective of thisexperiment is two-fold. First, we would like to validate the convergence propertiesanalysed in Section 4. Second, we would like to demonstrate the influence of the sizeof starting vectors on our new method.In each iteration there are t − s spurious eigenvalues. The spurious eigenvaluesoutside the target region can be easily detected according to their coordinates. Forthe spurious eigenvalues inside the target region, in [32] the authors introduced atolerance η to filter them. The idea behind is that the spurious eigenvalues can not http://math.nist.gov/MatrixMarket/ 12 teration l og R es -16-14-12-10-8-6-4-2 Problem 1 t = ⌈ . ⌉ t = ⌈ . ⌉ Iteration l og R es -14-12-10-8-6-4-2 Problem 2 t = ⌈ . ⌉ t = ⌈ . ⌉ Iteration l og R es -14-12-10-8-6-4-2 Problem 3 t = ⌈ . ⌉ t = ⌈ . ⌉ Iteration l og R es -7-6-5-4-3-2-1 Problem 4 t = ⌈ . ⌉ t = ⌈ . ⌉ Iteration l og R es -12-10-8-6-4-2 Problem 5 t = ⌈ . ⌉ t = ⌈ . ⌉ Iteration l og R es -12-10-8-6-4-2 Problem 6 t = ⌈ . ⌉ t = ⌈ . ⌉ Fig. 5.1 . The maximum relative residual norms in different iterations. achieve high accuracy; as the iteration process proceeds, there will be a gap in accuracybetween the desired eigenpairs and the spurious ones. If the relative residual norm ofan eigenpair is less that η , then the eigenpair is viewed as desired one and referred asfiltered eigenpair. In the experiment, we set the filtering tolerance η = 1 . × − . In Fig
Res ’s from the iteration that the number of filtered eigenpairsattains s for the first time to the 10th iteration for the cases t = ⌈ . s ⌉ and t = ⌈ . s ⌉ ,respectively. Here, we assume that the number s of eigenvalues inside the region ofinterest is known. Theorem 4.3 tells us the residual norm will converge with the factor | ˜ f ( λ l ′ ) | / | ˜ f ( λ i ) | for the exact eigenpair ( λ i , x i ) with respect to the iteration counts. Fig
Res decreases monotonically, asexpected, until the accuracy can not be further improved. On the other hand, alarger subspace size t leads to a smaller | ˜ f ( λ l ′ ) | , and then leads to faster convergence.Taking the Problem 2 as an example, it is clear to see that our method convergesalmost linearly with a factor for both t = ⌈ . s ⌉ and t = ⌈ . s ⌉ . Precisely, theconvergence rate is ablout 1 . × − for the former case and is about 1 . × − forthe latter. To converge to the minimum residual norm, which is about 1 . × − ,it needs 4 iterations when we take the size t to ⌈ . s ⌉ , but 10 iterations are requiredfor the case t = ⌈ . s ⌉ .Increasing the value of t will lead to a faster convergence rate, however, it alsoresults in a considerable increase in computational cost in each iteration since t rep-resents the number of the right-hand sides in each shifted linear system involved (see(3.13)). Both RFEAST and BFEAST [32] aim tomake the FEAST algorithm applicable for the non-Hermitian problems. The onlydifference between the two non-Hermitian FEAST methods is the choice of the leftsubspace. In BFEAST, the left subspace is spanned by BU (see (3.7)), while in ourmethod the left subspace is spanned by a random matrix. The dominant work ofboth methods is computing the approximate projection U (see (3.12)). Thus thecomputational cost required by both non-Hermitian FEAST algorithms in each iter-ation is almost the same. Due to this, in this experiment we compare the numericalperformance of the two methods through the accuracy achieved in each iteration.In [32], the authors presented a technique to select a suitable size of the starting teration l og R es -16-14-12-10-8-6-4-2 Problem 1
BFEASTRFEAST
Iteration l og R es -14-12-10-8-6-4-2 Problem 2
BFEASTRFEAST
Iteration l og R es -8-7-6-5-4-3-2-1 Problem 3
BFEASTRFEAST
Iteration l og R es -8-7-6-5-4-3-2-1 Problem 4
BFEASTRFEAST
Iteration l og R es -14-12-10-8-6-4 Problem 5
BFEASTRFEAST
Iteration l og R es -14-12-10-8-6-4-2 Problem 6
BFEASTRFEAST
Fig. 5.2 . The convergence behavior of two non-Hermitian FEAST algorithms. vectors Y for the BFEAST algorithm. To facilitate the comparisons, here we alsouse this technique to start our method. We depict Res ’s computed by the two testmethods from the iteration that the number of filtered eigenpairs attains s for the firsttime to the 10th iteration in Fig η is also taken to 1 . × − . The convergence curves of two non-HermitianFEAST methods are almost parallel, which means the two methods converge withalmost the same rate. We have shown in Theorem 4.3 that the upper bound forthe residual norms of exact eigenpairs ( λ j , x j ) are σ ( k ) τ j (cid:0) | ˜ f ( λ l ′ ) | / | ˜ f ( λ j ) | (cid:1) k in the k th iteration (see (4.17)). Recall that our method shares the same right subspacewith the BFEAST algorithm. In view of the proof of Theorem 4.3, we are able toestablish a similar upper bound for the BFEAST algorithm simply via replacing theoblique projector P ( k ) in the expression of σ ( k ) with Z ( k ) , where Z ( k ) is the obliqueprojector onto span { U ( k ) } and orthogonal to the left subspace B K . More precisely,we can write the upper bound for BFEAST as κ ( k ) τ j (cid:0) | ˜ f ( λ l ′ ) | / | ˜ f ( λ j ) | (cid:1) k , where κ ( k ) = k Z ( k ) ( A − λB )( I n − Q ( k ) ) k . Therefore, the two methods have almost the sameconvergence rate, which is | ˜ f ( λ l ′ ) | / | ˜ f ( λ j ) | for the eigenpair ( λ j , x j ). This can interpretwhy two non-Hermitian FEAST algorithms exhibit essentially the same convergencebehavior.On the other hand, it can be seen from Fig κ ( k ) and σ ( k ) , due to thedifferent choices of the left subspaces. The BFEAST algorithm works better than ourmethod possibly because the constant κ ( k ) in the BFEAST algorithm is smaller than σ ( k ) in our method, and therefore the upper bound in BFEAST is sharper than theone in our our method. eig function. In this experiment, we com-pare our method with the
Matlab built-in function eig in terms of timing. Sincethe target eigenvalues are the interior ones of non-Hermitian problems, when using eig to compute the eigenvalues inside the regions presented in
Tab ǫ to able 5.2 Comparison of eig and our method in terms of timing.
No. eig
Our method1 13.87 3.222 14.13 10.943 68.40 35.964 2719 .
36 180.075 2397.71 27.586 77653.86 398.82 . × − and take the parameter max iter = 10.The amount of time, which is measured in seconds, required by eig and ourRFEAST algorithm is reported in Table
Matlab function eig , although the parallelism offered by ourmethod is not used in the tests. The difference in CPU times is more obvious whenthe size of test problem grows larger. Therefore, our method is much more efficientthan the
Matlab ’s eig function.
6. Conclusions.
In this work, we have developed a new scheme to make theFEAST algorithm applicable for the non-Hermitian problems. The key step is thatthe left subspace used to extract the desired eigenpairs in our method is spannedby a random matrix. Theoretical analysis shown that our method can deal withthe non-Hermitian cases. The resulting method retains the feature of parallelismoffered by the original FEAST algorithm and does not increase the computationalcost. The convergence properties of our new method were also investigated. Numericalexperiments were reported to demonstrate the numerical performance of our newmethod and to validated the convergence analysis.
7. Acknowledgment.
I would like to thank Professor Raymond H. Chan, mythesis advisor, at The Chinese University of Hong Kong and Professor Man-ChungYeung at University of Wyoming for their help and fruitful discussions in the prepa-ration of this paper.
REFERENCES[1]
L. Ahlfors , Complex Analysis , 3rd Edition, McGraw-Hill, Inc., 1979.[2]
A. P. Austin and L. N. Trefethen , Computing eigenvalues of real real symmetric matriceswith rational filters in real arithmetic , preprint.[3]
Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. Van Der Vorst , Templates for theSolution of Algebraic Eigenvalue Problems: A Practical Guide , SIAM, Philadelphia, 2000.[4]
Z. Bai, J. Demmel, and M. Gu , An inverse free parallel spectral divide and conquer algorithmfor nonsymmetric eigenproblem , Numer. Math., 76 (1997), pp. 279–308.[5]
W.-J. Beyn , An integral method for solving nonlinear eigenvalue problems , Lin. Alg. Appl.,436 (2012), pp. 3839–3863[6]
K. A. Cliffe, A. Spence, and S. J. Tavener , The numerical analysis of bifurcation problemswith application to fluid mechanics , Acta Numer., 9 (2000), pp. 39–131.[7]
P. J. Davis and P. Rabinowitz , Methods of numerical integration , 2nd Edition, AcademicPress, Orlando, FL, 1984.[8]
J. Demmel , Applied Numerical Linear Algebra , SIAM, Philadelphia, 1997.[9]
E. Di Napoli, E. Polizzi, and Y. Saad , Efficient estimation of eigenvalue counts in aninterval , http://arxiv.org/abs/1308.4275.[10]
B. Ford and G. Hall , The generalized eigenvalue problem in quantum chemistry , Comput.Phys. Commun., 8 (1974), pp. 337–348.[11]
K. Gallivan, E. Grimme, and P. Van Dooren , A rational Lanczos algorithm for modelreduction , Numer. Algorithms, 12 (1996), pp. 33–64.1512]
F. R. Gantmacher , The Theory of Matrices , Chelsea, New York, 1959.[13]
G. H. Golub and C. F. Van Loan , Matrix Computations , 3rd Edition, Johns Hopkins Uni-versity Press, Baltimore, MD, 1996.[14]
R. G. Grimes, J. D. Lewis, and H. D. Simon , A shifted block Lanczos algorithm for solvingsparse symmetric generalized eigenproblems , SIAM J. Matrix Anal. Appl., 15 (1994), pp.228–272.[15]
A. Imakura, L. Du, and T. Sakurai , Error bounds of Rayleigh-Ritz type contour integral-based eigensolver for solving generalized eigenvalue problems , Numer. Algor., (accepted).[16]
T. Ikegami and T. Sakurai , Contour integral eigensolver for non-Hermitian systems: aRayleigh-Ritz-type approach , Taiwanese J. Math., 14 (2010), pp. 825–837.[17]
T. Ikegami, T. Sakurai, and U. Nagashima , A filter diagonalization for generalized eigen-value problems based on the Sakurai-Sugiura projection method , J. Comp. Appl. Math.,233 (2010), pp. 1927–1936.[18]
L. Kr¨amer, E. Di Napoli, M. Galgon, B. Lang, and P. Bientinesi , Dissecting the FEASTalgorithm for generalized eigenproblems , J. Comput. Appl. Math., 244 (2013), pp. 1–9.[19]
C. B. Moler and G. W. Stewart , An algorithm for generalized matrix eigenvalue problems ,SIAM J. Numer. Anal., 10 (1973), pp. 241–256.[20]
Y. Nakatsukasa and N. J. Higham , Stable and efficient spectral divide and conquer algo-rithms for the symmetric eigenvalue decomposition and the SVD , SIAM J. Sci. Comput.,35 (2013), pp. A1325–A1349.[21]
E. Polizzi , Density-matrix-based algorithm for solving eigenvalue problems , Phys. Rev. B 79(2009), 115112.[22]
W. P. Reinhardt , Complex coordinates in the theory of atomic and molecular structure anddynamics , Annu. Rev. Phys. Chem., 33 (1982), pp. 223–255.[23]
A. Ruhe , Rational Krylov: A practical algorithm for large sparse nonsymmetric matrix pencils ,SIAM J. Sci. Comput., 19 (1998), pp. 1535–1551.[24]
Y. Saad , Numerical Methods for Large Eigenvalue Problems , SIAM, Philadelphia, 2011.[25]
T. Sakurai, Y. Futamura, and H. Tadano , Efficient parameter estimation and implementa-tion of a contour integral-based eigensolver , J. Alg. Comput. Tech., 7 (2013), pp. 249–269.[26]
T. Sakurai and H. Sugiura , A projection method for generalized eigenvalue problems usingnumerical integration , J. comput. Appl. Math., 159 (2003), pp. 119–128.[27]
T. Sakurai and H. Tadano , CIRR: A Rayleigh–Ritz type method with contour integral forgeneralized eigenvalue problems , Hokkaido Math. J. 36 (2007), pp. 745–757.[28]
B. K. Sriperumbudur, D. A. Torres, and G. R. G. Lanckriet , A majorization-minimizationapproach to the sparse generalized eigenvalue problem , Mach. Learn., 85 (2011), pp. 3–39.[29]
G. W. Stewart , Matrix Algorithms, Vol. II, Eigensystems , SIAM, Philadelphia, 2001.[30]
P. T. P. Tang and E. Polizzi , FEAST as a subspace iteration eigensolver accelerated byapproximate spectral projection , SIAM J. Matrix Anal. Appl., 35 (2014), pp. 354–390.[31]
G. Yin , A Contour-integral Based Method for Counting the Eigenvalues Inside a Region inthe Complex Plane , https://arxiv.org/abs/1503.05035.[32]