A harmonic FEAST algorithm for non-Hermitian generalized eigenvalue problems

Linear Algebra and its Applications 578 (2019) 75–94

Contents lists available at ScienceDirect

Linear Algebra and its Applications www.elsevier.com/locate/laa

A harmonic FEAST algorithm for non-Hermitian generalized eigenvalue problems Guojian Yin School of Mathematics, Sun Yat-sen University, Guangzhou, PR China

a r t i c l e

i n f o

Article history: Received 1 January 2018 Accepted 25 April 2019 Available online 2 May 2019 Submitted by V. Mehrmann MSC: 15A18 58C40 65F15 Keywords: Generalized eigenvalue problems Generalized Schur form Contour integral Spectral projection

a b s t r a c t The FEAST algorithm, a contour-integral based eigensolver, was developed for computing the eigenvalues inside a given interval, along with their eigenvectors, of a Hermitian generalized eigenproblem. The FEAST algorithm is a fast and stable technique, and is easily parallelizable. However, it was shown that this algorithm may fail to find the desired eigenpairs when applied to non-Hermitian problems. Efforts have been made to adapt FEAST to non-Hermitian problems. In this work, we aim at formulating a new non-Hermitian scheme for the FEAST algorithm. Our new scheme is based on the partial generalized Schur decomposition. Unlike the existing nonHermitian FEAST methods, our method seeks approximation to the generalized Schur vectors instead of the potentially ill-conditioned eigenvectors. Our new method can compute not only the eigenvectors but also the right and left generalized Schur vectors corresponding to target eigenvalues. We present standard benchmark numerical experiments in which our method achieves better stability and efficiency comparing with the existing non-Hermitian FEAST algorithms. © 2019 Elsevier Inc. All rights reserved.

E-mail address: [email protected]. https://doi.org/10.1016/j.laa.2019.04.036 0024-3795/© 2019 Elsevier Inc. All rights reserved.

76

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

1. Introduction Large-scale eigenvalue problems arise in various areas of science and engineering, such as dynamic analysis of structures [15], linear stability analysis of the Navier-Stokes equation in fluid dynamics [6], the electron energy and position problems in quantum chemistry [12], and resonant state calculation [25]; see [35] for more examples. We consider the generalized eigenproblem Ax = λBx,

(1.1)

where A and B are general square matrices of order n, the scalars λ and associated nonzero vectors x are the eigenvalues and their associated (right) eigenvectors, respectively [3,8,14]. (1.1) is the most general form of eigenvalue problem. If B is an identity matrix, then (1.1) becomes a standard eigenproblem. In this work, we are concerned with computing eigenvalues located inside a given region in the complex region, along with their eigenvectors. In practical applications, it is not the whole spectrum, but rather one of its components that is of interest to the users [3,23,32]. When (1.1) is a dense eigenproblem, we can use the QZ method [22], which is an analog of the well-known QR method, to solve all the eigenvalues and then determine the target eigenvalues belonging to the desired spectral component. However, the QZ method requires O(n3 ) floating point operations, consequently, it is prohibitively expensive when (1.1) is a large-scale problem. In the practical situations, the problems encountered are large-scale and sparse. In the past decades, the most successful methods for approximating a component of the spectrum of a large sparse eigenproblem are based on projection techniques [3,4,32], of which perhaps the Krylov subspace approaches are the most wildly used [26,27]. However, the existing projection methods mainly focus on computing the extreme eigenvalues [31] or the eigenvalues close to a given shift [15]. Computing interior eigenvalues of a large-scale eigenproblem is still one of the most difficult problems in computational linear algebra today [10]. Recently, a class of contour-integral based eigensolvers were proposed for computing the eigenvalues inside a given region in the complex plane [2,5,9,24,28–30]. These novel eigensolvers have two appealing features: (i) replacing the difficulties of solving eigenproblem (1.1) by the difficulties of solving linear systems [24]; and (ii) possessing good potential to be computed on modern scalable parallel machines [24,29]. Two typical examples of the contour-integral based eigensolvers are the Sakurai-Sugiura (SS) method [29] and the FEAST algorithm developed by Polizzi in [24]. The SS method suffers from numerical instability [2,18] because it needs to compute the eigenvalues of ill-conditioned Hankel matrices. Later Sakurai et al. turned to use the Rayleigh-Ritz procedure to extract the target eigenpairs, which yielded a stable contour-integral based eigensolver, that is, the so-called CIRR method [17,30]. The FEAST algorithm and the CIRR method belong to the family of subspace iteration methods [33]. They can be briefly summarized as follows: first using specifically defined contour integrals to gen-

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

77

erate the invariant subspaces corresponding to the target eigenvalues, then using the orthogonal projection techniques to extract the target eigenpairs. The formulation of both CIRR and FEAST requires that matrices A and B are Hermitian matrices and B is positive definite, i.e., (1.1) is a Hermitian problem. It was pointed out in [37] that FEAST may fail to find the target eigenpairs when (1.1) is a nonHermitian problem. By noting this fact, the authors in [37] proposed a non-Hermitian scheme for the FEAST algorithm so as to extend it to the non-Hermitian cases. The key idea in their method is that using the oblique projection technique with an appropriately chosen left subspace, rather than the orthogonal projection technique used in FEAST, to extract the desired eigenpairs. In [19,21], a dual-subspace FEAST algorithm was developed for adapting the FEAST algorithm to the non-Hermitian problems. In the dual-subspace FEAST algorithm, two subspaces are needed to construct through two contour integrals at each iteration to contain the left and right eigenspaces, respectively, corresponding to the desired eigenvalues. As with the original FEAST algorithm [24], the dominant work of the non-Hermitian FEAST algorithm proposed in [37] at each iteration is computing one contour integral approximately by a quadrature rule. Consequently, if only the (right) eigenvalues or the (right) eigenvalues and their associated eigenvectors are needed, the computational cost required by the scheme proposed in [37] is about half of that required by the dual-subspace FEAST algorithm. In this work, we formulate another non-Hermitian scheme for the FEAST algorithm. Our motivation is to avoid the potential ill-conditioning of the eigenbasis. To this end, we extract the desired eigenpairs from the subspace spanned by generalized Schur vectors corresponding to the desired eigenpairs. We borrow this idea from the one used to formulate the well-known JDQZ method [11]. The resulting extraction approach is a harmonic Schur-Rayleigh-Ritz procedure [34]. The experiments demonstrate that our method is more stable and efficient than the existing non-Hermitian FEAST algorithms. In addition, our method can compute not only the eigenvalues inside a given region and the associated eigenvectors, but also the associated right and left generalized Schur vectors. The paper is organised as follows. In Section 2, we briefly describe the FEAST algorithm [24]. In Section 3, we review the non-Hermitian variant of the FEAST algorithm that is proposed in [37]. We formulate our new non-Hermitian FEAST algorithm in Section 4. In Section 5, numerical experiments are reported to illustrate the numerical performance of our method. Throughout the paper, we use the following notation and terminology. The subspace spanned by the columns of a matrix X is denoted by span{X}. The rank and conjugate transpose of X are denoted by rank(X) and X ∗ respectively. The algorithms are presented in Matlab style. We use the Matlab conventions when we refer to entries in matrices and vectors.

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

78

2. Introduction to FEAST We want to develop a new scheme to adapt the FEAST algorithm [24] to nonHermitian problems, so it is necessary for us to provide a brief introduction to the FEAST algorithm. Polizzi formulated the FEAST algorithm in [24] with the aim of computing the eigenvalues, along with the corresponding eigenvectors, inside a given interval of a Hermitian eigenvalue problem. Assume that the interval of interest is [σ1 , σ2 ] and the eigenvalues inside [σ1 , σ2 ] are {λi }si=1 , which are not necessarily distinct. Below we briefly describe how to use the FEAST algorithm to compute the target eigenpairs. Essentially, the FEAST algorithm belongs to the family of subspace iteration method [33]. However, unlike the well-known Krylov subspace methods [27], it constructs the invariant subspace corresponding to the desired eigenvalues via the contour integral defined as  1 √ F := (zB − A)−1 dzY, (2.1) 2π −1 Γ

where Γ is any Jordan contour that contains {λi }si=1 inside, and Y is an n × s random matrix. Below we show that the columns of F form a basis of the invariant subspace corresponding to the desired eigenvalues. We begin with the following theorem. Theorem 2.1 ([32]). Let A and B be n × n Hermitian matrices and that B is positive definite. Then there exists an n × n matrix X = [x1 , x2 , . . . , xn ] for which X ∗ BX = In

and

X ∗ AX = Λ = diag(λ1 , λ2 · · · , λn ),

(2.2)

where In is the n × n identity matrix, {λi }ni=1 are the eigenvalues of the matrix pencil zB − A, and the columns {xi }ni=1 of X are their associated eigenvectors. By (2.2) and the residue theorem in complex analysis [1], it is easy to verify that F =

1 √

2π −1

 X

(zIn − Λ)−1 dzY (BX)−1 = X(:,1:s) (X(:,1:s) )∗ Y.

(2.3)

Γ

Therefore, if (X(:,1:s) )∗ Y is full-rank, we can conclude that the columns of F form a basis for span{X(:,1:s) }, which is exactly the desired eigenspace. With this knowledge at hand, the desired eigenpairs can be extracted from subspace span{F } by a projection technique [3,27]. In the FEAST algorithm, the orthogonal projection technique is used [24]. The desired eigenpairs (λ, x) are computed by imposing the Galerkin condition: (Ax − λBx) ⊥ span{F }, which can be transformed into solving the projected eigenproblem

(2.4)

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

ˆ = λBy, ˆ Ay

79

(2.5)

with x = F y and Aˆ = F ∗ AF

ˆ = F ∗ BF. and B

(2.6)

Now the task of solving large-scale eigenproblem (1.1) is reduced to the one of solving the s × s eigenproblem (2.5). To get the projected eigenproblem (2.5), the most important task is forming the basis vectors F . In practice, F (see (2.1)) has to be computed numerically by a quadrature scheme. In [24], Polizzi used the q-point Gauss-Legendre quadrature [7] to compute F approximately with Γ being taken to the circle with center at γ = (σ1 + σ2 )/2 and radius ρ = (σ2 − σ1 )/2, which leads to F =

1 √ 2π −1



1 ωj (zj − γ)(zj B − A)−1 Y, 2 j=1 q

(zB − A)−1 dzY ≈

Γ

(2.7)

√

where zj = γ + ρe −1θj , θj = (1 + tj )π, and tj is the jth Gaussian node with associated weight ωj . From (2.7) we see that the dominant computational work of the FEAST algorithm at each iteration is solving q linear systems of the form (zj B − A)Xj = Y,

j = 1, . . . , q.

(2.8)

According to the description above, the FEAST algorithm can be summarised as follows. Algorithm 1: The FEAST Algorithm. Input:

Hermitian matrices A and B with B being positive definite, a uniformly-distributed random matrix Y ∈ Rn×t with t ≥ s, and the circle Γ enclosing the interval [σ1 , σ2 ]. ˆi , x ˆ i are located inside [σ1 , σ2 ]. ˆ i )}si=1 , λ Output: The approximate eigenpairs {(λ 1. Compute F approximately by (2.7). ˆ = F ∗ BF . ˆ = F ∗ AF and B 2. Set A ˆ By, ˆ i , yi )}t . ˆ =λ ˆ to obtain the eigenpairs {(λ 3. Solve the eigenproblem: Ay i=1 ˆ i = F yi , i = 1, 2, . . . t. 4. Compute x ˆ t ] and 5. If s approximate eigenpairs converge, stop. Otherwise, set Xt = [ˆ x1 , . . . , x Y = BXt , then go back to Step 1.

The FEAST algorithm is an accurate and stable technique [20,33] that recently attracted much attention. The dominant computational cost is solving the linear systems of the form (2.8), which means that the FEAST algorithm replaces the difficulties of solving the eigenproblem (1.1) by those of solving linear systems (2.8) [24]. The linear systems can be solved by any method of choice. On the other hand, since the quadrature nodes zj and the columns of the right-hand sides in (2.8) are independent of each other, the involved linear systems can be solved in parallel.

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

80

3. A non-Hermitian FEAST algorithm Initially, Polizzi formulated the FEAST algorithm under the assumptions that A and B are Hermitian and B is positive definite, that is, the considered eigenproblem (1.1) is Hermitian [24]. When applying the FEAST algorithm to non-Hermitian problems, it was observed in [37] that it may fail to find the desired eigenpairs. The authors in [37] gave a simple example to illustrate this fact and developed a non-Hermitian version for the FEAST algorithm. The key to the success of their method is that the oblique projection technique, rather than the orthogonal projection technique used in the FEAST algorithm, with an appropriately chosen left subspace is used to extract desired eigenpairs. We are now in a position to derive the non-Hermitian FEAST method proposed in [37]. We begin with the introduction of the canonical form for a regular matrix pencil. A matrix pencil zB − A is regular if det(zB − A) is not identically zero for all z ∈ C, which is the most common case in the generalized eigenvalue problems [3]. Theorem 3.1 (The Weierstrass canonical form [13,37]). Let zB − A be a regular matrix pencil of order n. Then there exist nonsingular matrices S and T ∈ C n×n such that 

J T AS = d 0

0





and

In−d

I T BS = d 0

0 Nn−d

 ,

(3.1)

where Jd is a d × d matrix in Jordan canonical form with its diagonal entries corresponding to the eigenvalues of zB − A, Nn−d is an (n − d) × (n − d) nilpotent matrix also in Jordan canonical form, and Id denotes the identity matrix of order d. Let Jd be of the form ⎡

Jd1 (λ1 ) ⎢ 0 ⎢ Jd = ⎢ .. ⎢ . ⎣ 0 where

m i=1

0 Jd2 (λ2 ) .. . 0

··· ··· .. . ···

⎤ 0 ⎥ 0 ⎥ ⎥ .. ⎥ . ⎦ Jdm (λm )

(3.2)

di = d and Jdi (λi ) are di × di matrices of the form ⎡

λi

⎢ ⎢ 0 ⎢ ⎢ Jdi (λi ) = ⎢ ⎢ ⎢ . ⎢ . ⎣ . 0

1

0

···

λi .. .

1 .. .

..

.

..

..

.

···

.

0

⎤ 0 .. ⎥ . ⎥ ⎥ ⎥ ⎥ 0 ⎥, ⎥ ⎥ 1 ⎦ λi

i = 1, 2, . . . , m

with λi being the eigenvalues. Here the λi are not necessarily distinct and can be repeated according to their multiplicities.

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

81

Let Nn−d be of the form ⎡

Nn−d

where

m i=1

0 Nd2 .. . 0

Nd1 ⎢ 0 ⎢ =⎢ ⎢ .. ⎣ . 0

··· ··· .. . ···

⎤

0 0 .. .

⎥ ⎥ ⎥, ⎥ ⎦

Ndm

di = n − d and Ndi are di × di matrices of the form ⎡

Ndi

0

1

⎢ ⎢0 0 ⎢ ⎢ .. =⎢ . ⎢ ⎢. ⎢. ⎣. 0 ···

0

···

1 .. .

..

.

..

..

.

.

⎤ 0 .. ⎥ .⎥ ⎥ ⎥ ⎥ 0⎥, ⎥ ⎥ 1⎦ 0

0

i = 1, 2, . . . , m .

Let us partition S into block form S = [S1 , S2 , . . . , Sm , Sm+1 ], where each Si ∈ C n×di , 1 ≤ i ≤ m, and Si into Si = [si1 , si2 , . . . , sidi ] with sij ∈ C n , 1  j  di . It was verified in [37] that for any eigenvalue λi , 1 ≤ i ≤ m, 

I (λi B − A)S d 0

0 Nn−d





λ I − Jd = BS i d 0

0

λi Nn−d − In−d

 .

(3.3)

By comparing the first d columns on both sides above, we get (λi B − A)sij = Bsij−1 ,

1 ≤ j ≤ di ,

1 ≤ i ≤ m,

(3.4)

with si0 ≡ 0. It can be seen that si1 are the eigenvectors corresponding to the eigenvalues λi for all 1 ≤ i ≤ m. Let Γ be a positively oriented simple closed curve enclosing the desired eigenvalues. Again without loss of generality, we let the eigenvalues of (1.1) enclosed by Γ be {λ1 , . . . , λl }, and s := d1 + d2 + · · · + dl be the number of eigenvalues inside Γ with multiplicity taken into account. Define the contour integral 1 √ Q := 2π −1



(zB − A)−1 Bdz.

(3.5)

Γ

According to the residue theorem in complex analysis [1], it was verified in [37] that 1 √ Q= 2π −1

 Γ

(zB − A)−1 Bdz = S

Is 0

 0 S −1 = S(:,1:s) (S −1 )(1:s,:) . 0

(3.6)

82

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

One can easily verify that Q2 = Q, which means Q is a projector onto the subspace K := span{S(:,1:s) }. Define U := QY = S(:,1:s) (S −1 )(1:s,:) Y,

(3.7)

where Y is an n × s random matrix. Obviously, U is the projection of Y onto the subspace K. In order to show that the columns of U form a basis for the subspace K, we need the following result. Lemma 3.2 ([37]). Let Y ∈ Rn×s be a sample of a random matrix whose entries are independent identically distributed (i.i.d.) continuous distributions, then with probability 1, the matrix (S −1 )(1:s,:) Y is nonsingular. According to (3.7) and Lemma 3.2, we can conclude that the columns of U form a basis for the subspace K almost surely if Y is a random matrix. Here we refer the reader to [16], in which the authors discussed detailedly what random matrix should we use and how much sampling do we need. Note that K contains the eigenspace corresponding to the desired eigenvalues (see (3.4) for details), it is natural to take K as the right subspace (search subspace). Recall that the FEAST algorithm takes advantage of the often used orthogonal projection technique to extract the desired eigenpairs, in which case the left subspace (test subspace) is the same as the right one. The authors in [37] observed that this extraction approach may result in the FEAST algorithm failing to find the desired eigenpairs when the problem involved is non-Hermitian. To address this deficiency, they resorted to the oblique projection method and developed a non-Hermitian scheme for the FEAST algorithm in [37]. In their scheme the left subspace is taken as BK, rather ¯ x ¯ ) are obtained by imposing the than K. Consequently, the approximate eigenpairs (λ, Petrov-Galerkin condition [3,27]: ¯ x ¯ ) ⊥ BK, (A¯ x − λB

(3.8)

¯ ∈ C and x ¯ ∈ K. It was shown in [37] that the columns of BU form a basis of where λ BK. Therefore (3.8) can be written in matrix form ¯ (BU )∗ (AU y − λBU y) = 0,

(3.9)

¯ = U y. Accordingly, solving the eigenvalues of (1.1) inside Γ where y ∈ C s satisfying x now is reduced to solving the smaller projected eigenproblem ¯ By, ¯ =λ ¯ Ay

(3.10)

with A¯ = (BU )∗ AU

¯ = (BU )∗ BU. and B

(3.11)

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

83

The following theorem justifies their choice of the left subspace. ¯ i , yi )}s be the eigenpairs of the projected eigenproblem Theorem 3.3 ([37]). Let {(λ i=1 s ¯ i , U yi )} (3.10). Then {(λ i=1 are the eigenpairs of (1.1) located inside Γ. In order to generate the projected eigenproblem (3.10), the most important task is to compute the projection U (see (3.7)). In practice, we have to compute U approximately by a quadrature rule:

U=

1 √

2π −1



˜= (zB − A)−1 BdzY ≈ U

Γ

1 √

2π −1

q 

ωj (zj B − A)−1 BY,

(3.12)

j=1

where zj are the quadrature nodes on Γ associated with weights ωj . Let  1 √ ωj (zj B − A)−1 B. 2π −1 j=1 q

˜= Q

(3.13)

˜ is an approximation of spectral projector Q (see (3.5)), which can be viewed as Then Q ˜ = QY ˜ . a rational filter [18,33,36], and U From (3.12), we see that the dominant work at each iteration is solving q linear systems of the form (zi B − A)Xi = BY,

j = 1, . . . , q.

(3.14)

The non-Hermitian FEAST algorithm described above uses the oblique projection technique to extract the desired eigenvalues, the right subspace is different from the left one. However, since the left subspace is spanned by BK, the dominant computational work at each iteration is solving q linear systems, as with the FEAST algorithm. Below we call the non-Hermitian FEAST algorithm described above BFEAST for the ease of reference. Algorithm 2: The BFEAST Algorithm. The matrices A and B, a uniformly-distributed random matrix Y ∈ Rn×t , where t ≥ s, and the Jordan curve Γ. ¯ i are located inside Γ. ¯i , x ¯ i )}si=1 , λ Output: The approximate eigenpairs {(λ 1. Compute U approximately by the quadrature rule (3.12). 2. Compute QR decompositions: U = U1 R1 and BU = U2 R2 . ¯ = U ∗ AU1 and B ¯ = U ∗ BU1 . 3. Form A 2 2 ¯ By ¯ =λ ¯ of size t to obtain eigenpairs 4. Solve the projected eigenproblem Ay ¯ i , yi )}t . Set x ¯ i = U1 yi , i = 1, 2, . . . , t. {(λ i=1 5. If s approximate eigenpairs converge, stop. Otherwise, set Y = U1 and go to Step 1. Input:

84

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

4. A new non-Hermitian FEAST algorithm In the preceding section, we reviewed a non-Hermitian variation of the FEAST algorithm, i.e., BFEAST. By the Weierstrass canonical form (3.1), we have AS(:,1:s) = BS(:,1:s) J(1:s,1:s) .

(4.1)

The diagonal entries of J(1:s,1:s) are the desired eigenvalues and S(:,d0 +...,di−1 +1) , d0 = 0, are the eigenvectors corresponding to the desired eigenvalues λi , i = 1, . . . , l. Especially, when the desired eigenvalues are nondefective, which means their algebraic and geometric multiplicities coincide, J(1:s,1:s) is a diagonal matrix and the columns S(:,i) are the eigenvectors corresponding to λi , i = 1, . . . , s. In this case, (4.1) is the partial eigenvalue decomposition of (A, B) [32,34]. The BFEAST algorithm computes a basis for K = span{S(:,1:s) }, which contains the desired eigenspace, through the contour integral defined as (3.7). It uses the oblique projection technique to extract the desired eigenpairs from the right subspace K with left subspace taken to BK. However, it is well-known that the eigenbasis can potentially be ill-conditioned when the problem under consideration is defective [32,34], in which case the algebraic multiplicities of some eigenvalues exceeds their geometric multiplicities [14]. To overcome the difficulties related to ill-conditioning of eigenbasis, we want to formulate a new non-Hermitian FEAST algorithm, with the hope of formulating a more stable scheme. To this end, we generate the generalized Schur vectors instead of potentially ill-conditioned eigenbasis via the contour-integral based technique. It should be pointed out that this idea is borrowed from the formulation of the well-known JDQZ method [11] and the GPLHR method developed in [34]. In what follows, we derive our new scheme detailedly. 4.1. The derivation of new method We begin our derivation with the definition of partial generalized Schur decomposition. Definition 4.1 ([11]). A partial generalized Schur form of dimension s for a matrix pair (A, B) is the decomposition AQs = Zs GA ,

BQs = Zs GB ,

(4.2)

where Qs and Zs are orthogonal n × s matrices, and GA and GB are upper triangular s × s matrices. A column (Qs )(:,i) is referred to as a generalized Schur vector, and we refer to a pair ((Qs )(:,i) , (GA )(i,i) /(GB )(i,i) ) as a generalized Schur pair. In the preceding section, we have shown that K = span{U } (see (2.1)) and K contains the eigenspace corresponding to the eigenvalues inside Γ. It is natural to take span{U } as the right subspace. For the time being, we still consider using oblique projection technique to extract the desired eigenpairs from span{U }. Let V be an orthogonal matrix

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

85

whose columns form an orthonormal basis of span{U }. Now the main task is to seek an appropriate left subspace. Suppose such a left subspace is already available, and the columns of W form its orthonormal basis. Then the desired approximate eigenpairs (λ, x) can be computed by imposing the Petrov-Galerkin condition [27]: (Ax − λBx) ⊥ span{W },

(4.3)

which can be transformed into solving the projected eigenproblem: W ∗ AV y = λW ∗ BV y,

(4.4)

where y ∈ C s satisfies x = V y. Applying the QZ decomposition to pair (W ∗ AV, W ∗ BV ) to yield generalized Schur form [14]: (PL )∗ (W ∗ AV )PR = HA

and (PL )∗ (W ∗ BV )PR = HB ,

(4.5)

where PR and PL are orthogonal s × s matrices, HA and HB are upper triangular s × s matrices. The eigenvalues of pair (W ∗ AV, W ∗ BV ) are {(HA )(i,i) /(HB )(i,i) }si=1 [14,22]. The formulation (4.2) is equivalent to (Zs )∗ AQs = GA ,

(Zs )∗ BQs = GB ,

(4.6)

from which we know that (GA )(i,i) /(GB )(i,i) are the eigenvalues of (GA , GB ). Let yi be the eigenvectors of pair (GA , GB ) associated with (GA )(i,i) /(GB )(i,i) , then ((GA )(i,i) /(GB )(i,i) , Qs yi ) are the eigenpairs of (A, B) [11,22]. Comparing (4.6) with (4.5), we see that we have constructed a partial generalized Schur form in (4.5) for matrix pair (A, B): V PR and W PL are the right and left generalized Schur vectors, respectively, along with the upper triangular matrices HA and HB . Note that V = orth(U ), to get the partial generalized Schur vectors V PR and W PL , our task is to compute the basis matrix W . To this end, we borrow an idea used to formulate the well-known JDQZ algorithm [11]. In view of (4.2), it was shown in [11] that span{Zs } = span{AQs } = span{BQs }.

(4.7)

Since V PR forms a Qs and W PL forms a Zs , we have span{W PL } = span{AV PR } = span{BV PR }.

(4.8)

On the other hand, since PL and PR are orthogonal s × s matrices, we have span{W } = span{AV } = span{BV }.

(4.9)

86

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

Motivated by (4.9), we set the left subspace span{W } as span{AV } ∪ span{BV }. Now we are ready to state the main result of our work. Theorem 4.2. Let V = orth(U ) and W = orth(AV − σBV ), where σ is a scalar and is different from the eigenvalues of (A, B). Applying the QZ decomposition to projected pair (W ∗ AV, W ∗ BV ) to generate generalized Schur form (4.5). Let yi be the eigenvectors of (HA , HB ) corresponding to eigenvalues (HA )(i,i) /(HB )(i,i) , then {((HA )(i,i) /(HB )(i,i) , V PR yi )}si=1 are exactly the eigenpairs of (1.1) inside Γ, that is, the desired eigenpairs. Proof. We have shown that (4.5) is a partial generalized Schur form of pair (A, B), therefore {((HA )(i,i) /(HB )(i,i) , V PR yi )}si=1 are the eigenpairs of (A, B). Note that (4.5) is generated by applying QZ decomposition to pairs (W ∗ AV, W ∗ BV ). As a result, to complete our proof, we only need to prove that the eigenvalues of (W ∗AV, W ∗ BV ) are exactly the eigenvalues of (1.1) inside Γ. By (4.9) we know that there is an s × s nonsingular matrix Δ satisfying AV = BV Δ.

(4.10)

Since W = orth(AV − σBV ), there exists a matrix D such that W = (AV − σBV )D = BV (Δ − σIs )D.

(4.11)

By (4.10) and (4.11), we know that the characteristic polynomial of the eigenproblem (W ∗ AV, W ∗ BV ) is det(λW ∗ BV − W ∗ AV ) = det[D∗ (Δ − σIs )∗ (BV )∗ BV (λIs − Δ)].

(4.12)

In view of (4.12), to prove our result, first we want to show that D is a nonsingular matrix. Since σ is not an eigenvalue of (A, B), A −σB is nonsingular. Since U is full-rank and V is an orthonormal basis of span{U }, we can verify that span{(A − σB)V } is an s-dimensional subspace. Note that W is an orthonormal basis of span{(A − σB)V }, we can conclude D is an s × s nonsingular matrix. Next, we want to show that (BV )∗ BV is nonsingular. By (3.1) and (3.7), we have BU = BS(:,1:s) (S −1 )(1:s,:) Y = (T −1 )(:,1:s) (S −1 )(1:s,:) Y.

(4.13)

Therefore, according to Lemma 3.2, we can conclude that BU is full-rank. Consequently, it is easy to conclude that (BV )∗ BV is nonsingular. Finally, we want to show that Δ − σIs is nonsingular. Suppose Δ has eigendecomposition Δ = XΔ ΛΔ (XΔ )−1 . Substituting the eigendecomposition into (4.10) we get A(V XΔ ) = B(V XΔ )ΛΔ , by which we can conclude that the eigenvalues of Δ are exactly the eigenvalues of (A, B) inside Γ. Recall that σ is not an eigenvalue of (A, B), thus Δ − σIs should be nonsingular.

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

87

According to above analysis, now we know that the roots of the characteristic polynomial det(λW ∗ BV − W ∗ AV ) are the ones of det(λIs − Δ), which are exactly the desired eigenvalues. The proof is completed. 2 According to Theorem 4.2, we know that the target eigenpairs can be obtained by applying the QZ algorithm to the projected pair (W ∗ AV, W ∗ BV ). Unlike the BFEAST algorithm (Algorithm 2), the desired eigenpairs are extracted from the subspace spanned by V PR , whose columns are the generalized Schur vectors corresponding to the target eigenvalues. Note that there is a parameter, i.e., σ, needed to select when constructing the subspace span{W }. The parameter σ can be viewed as a shift in the sense of seeking the eigenvalues close to σ [3,15,26]. Recall that we want to find the eigenvalues inside a Jordan contour. Therefore, it makes sense to select σ to the “center” of the target region and to be different from the desired eigenvalues. Then, the extraction approach in our scheme is the so-called harmonic Schur-Rayleigh-Ritz procedure [34]. For ease of reference, we call our new non-Hermitian FEAST algorithm HFEAST and summarize it as follows. Algorithm 3: The HFEAST Algorithm. The matrices A and B, a uniformly-distributed random matrix Y ∈ Rn×t , where t ≥ s, the Jordan curve Γ, a constant σ is chosen to be both inside Γ and in the resolvent set of the pencil (A, B). ˜ i are located inside Γ. ˜i , x ˜ i )}si=1 , λ Output: The approximate eigenpairs {(λ ˜ by (3.12). 1. Compute the approximate projection U ˜ ) and W = orth(AV − σBV ). 2. Compute orthogonalization: V = orth(U 3. Compute [HA , HB , PL , PR , VL , VR ] = qz(W ∗ AV, W ∗ BV ). ˜ i = (HA )(i,i) /(HB )(i,i) and x ˜ i = V PR (VR )(:,i) . 4. Compute λ 5. If s approximate eigenpairs converge, stop. Otherwise, set Y = V and go back to Step 1. Input:

Below we give some remarks on Algorithm 3 (HFEAST). 1. As with the BFEAST algorithm, the dominant work at each iteration in the HFEAST ˜ algorithm is computing q linear systems to get the approximation projection matrix U (see (3.12)) in Step 1. 2. In Step 3, the Matlab command qz is used to compute the QZ decomposition of projected pair (W ∗ AV, W ∗ BV ). 3. In Step 3, the right and left generalized Schur vectors V PR and W PL are generated. Meanwhile, with the left generalized Schur vectors W PL at hand, the left eigenvectors associated with the target eigenvalues can be solved. 4. Since the number s of the eigenvalues inside the target region is generally unknown in advance, one problem arises when applying Algorithm 3, that is, how to guarantee all desired eigenvalues are found when the algorithm stops. The authors in [37] proposed an empirical approach to determine the number s and then designed stopping criteria for the BFEAST algorithm, which are also applicable for our HFEAST algorithm.

88

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

In their empirical approach, a thresholding tolerance, say η, was introduced to select ˜ i, x ˜i) the desired approximation eigenvalues. To be more precise, for the eigenpairs (λ inside the target region (detected by checking the values of their coordinates), define the relative residual norms

ri =

˜iBx ˜ i 2 A˜ xi − λ , ˜ i 2 A˜ xi 2 + B x

(4.14)

˜i, x ˜ i ) are viewed as desired if ri are less than the thresholding tolerance η, then (λ approximation eigenpairs and referred to as filtered eigenpairs. As iterations proceed, the accuracy of desired approximation eigenpairs will be improved while the spurious eigenpairs inside the target region do not. As a result, there will exist a gap of accuracy between the desired approximation eigenpairs and the spurious eigenpairs inside the target region. Therefore, the number of filtered eigenvalues will remain unchanged and be equal to s from a certain iteration onwards. Based on this observation, the authors in [37] designed the stopping criteria for BFEAST as follows. If the numbers of filtered eigenpairs remain unchanged, meanwhile, all filtered eigenpairs attain the user-specified accuracy, then we stop the iteration process. ˜ } is refined by the filter Q ˜ (see (3.13)) 5. In Algorithm 3, the search subspace span{U iteratively. We do not deflate the directions corresponding to the spurious eigenvalues from the search subspace in the iteration procedure. This is because these directions may contain very useful information about the desired eigenspace. Actually, the filter ˜ helps to remove the redundant information contained in the search subspace during Q the iteration process. 6. We borrow an idea from the JDQZ method to formulate the resulting algorithm. However, unlike the JDQZ method, which expands the search subspace by solving a correction equation at each iteration, our method constructs a search subspace, with fixed dimension t, to contain the desired eigenspace directly by contour integrals (3.12). Then the desired approximation eigenpairs are extracted simultaneously by the harmonic Schur-Rayleigh-Ritz procedure. For the class of contour-integral based eigensolvers, which our method belongs to, the convergence rate of desired approximation eigenpair is related to the size of starting vectors t and the location of desired eigenvalue in the target region [18,28,33]. So deflating the converged eigenvalues in the iteration procedure will not improve the convergence rates of the unconverged ones. 5. Numerical experiments In this section, some numerical experiments are reported to demonstrate the numerical performance of our new non-Hermitian FEAST algorithm, that is, HFEAST. The

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

89

Table 5.1 Test problems from Matrix Market that are used in our experiments. No.

Problem

Type

n

Region: (γ, ρ)

s

1 2 3 4 5 6 7 8 9

RBS480 BFW782 DWG961 UTM1700 MHD4800 OLM5000 DW8192 CRY10000 BAURU5727

gen. gen. gen. gen. gen. stand. stand. stand. stand.

480 782 961 1700 4800 5000 8192 10000 40366

(0.2, 0.5) (−5.0 × 105 , 2.0 × 105 ) (−5.0 × 105 , 2.0 × 105 ) (4.0, 1.0) (−6.0, 3.0) (−1.0 × 104 , 0.6 × 104 ) (0.5, 0.03) (−5.0 × 102 , 0.3 × 102 ) (6.0 × 102 , 1.0 × 102 )

110 192 86 96 169 204 125 134 127

test problems are downloaded from the Matrix Market collection.1 They originate from the real-world applications. The present work aims at developing a more stable nonHermitian FEAST algorithm, so the test problems we select are non-Hermitian in the experiments. Their names are listed in Table 5.1. The first five test problems are generalized eigenproblems, while the remaining four test problems are standard ones. The target region for each test problem is a disk enclosed by a circle with center at γ and radius ρ. s is the number of eigenvalues inside the target region. In all experiments, we use the Gauss-Legendre quadrature with q = 16 quadrature nodes to compute the approximate ˜ (see (3.12)). The linear systems (see (3.14)) involved are computed by diprojection U rect method. We first use the Matlab function lu to compute the LU decompositions of A − zj B, j = 1, 2, . . . , q, and then perform the triangular substitutions to get the corresponding solutions. All computations are carried out in Matlab version R2014b on a MacBook with an Intel Core i5 2.5 GHz processor and 8 GB RAM. We use the maximum relative residual norm defined as Res = max1≤i≤s ri (see (4.14)) to assess the accuracy achieved by the test methods involved in the experiments. As with the BFEAST algorithm [37], we introduce a thresholding tolerance η to select the desired approximation eigenpairs. The approximation eigenpairs whose relative residual norms are less than η are viewed as desired ones and referred to as filtered eigenpairs. In the experiments, we set η = 1.0 ×10−2 , and the number of filtered eigenpairs at each iteration is denoted by sη . 5.1. Convergence behaviour of HFEAST The objective of the work is to develop a stable non-Hermitian FEAST algorithm. The FEAST algorithm is a fast and effective technique [20,33]. We hope our HFEAST algorithm retains the effectiveness and robustness of the FEAST algorithm. The aim of this experiment is two-fold. First, we want to show the convergence behaviour of our HFEAST algorithm for demonstrating its effectiveness. Second, we want to show the influence of the sizes of starting vectors on our new method. 1

http://math.nist.gov/MatrixMarket/.

90

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

Fig. 5.1. The convergence history associated with different sizes of starting vectors.

Our method rotates the search space, without reducing its dimension, towards the ˜ rather than deflating spurious directions desired algebraic eigenspace using the filter Q from it. Thus in our method t eigenvalues are computed at each iteration and t −s of them are spurious eigenvalues. The spurious eigenvalues which are located outside the target region can be easily detected by the values of their coordinates. We use thresholding tolerance η to select the desired approximation eigenpairs. The approximation eigenpairs whose relative residual norms are less than η are viewed as desired ones. However, the relative residual norms of some desired approximation eigenpairs may be larger than η at the first few iterations, consequently, it is difficult for us to distinguish them from the spurious eigenpairs inside the target region, which means that we are unable to obtain the values of Res at these iterations. So in Fig. 5.1 we plot Res at the iterations, where the number of filtered eigenvalues sη attains the actual number s, when the size t of starting vectors is taken to 1.2s, 1.4s and 1.6s, respectively. We can see from Fig. 5.1 that Res decreases monotonically and dramatically before accuracy can not be improved, which shows that the convergence of our HFEAST algorithm is fast. On the other hand, it can be observed that the larger the size t is, the faster we method converges. Note that t is also the size of right-hand sides. Therefore, larger value of t means that we need to compute more linear systems at each iteration. As a result, more computational costs are incurred at each iteration. 5.2. Comparison with the BFEAST algorithm The BFEAST algorithm extends the FEAST algorithm to the non-Hermitian cases successfully [37]. In this work, we aim at developing a more stable non-Hermitian scheme

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

91

Fig. 5.2. The maximum relative residual norms at different iterations.

for the FEAST algorithm. To this end, our HFEAST algorithm extracts the desired eigenpairs from the subspace spanned by generalized Schur vectors, rather than the eigenspace, corresponding to the desired eigenpairs. In this experiment, we compare our HFEAST algorithm with another non-Hermitian FEAST method, BFEAST, proposed in [37]. In this experiment, we assume that the number s of eigenvalues inside the target region is known, the size of starting vectors is taken to 1.3s for both non-Hermitian FEAST algorithms. As with the previous experiment, in Fig. 5.2, we depict the maximum relative residual norms Res at the iterations where the number of filtered eigenpairs sη attains the actual number s for both BFEAST and HFEAST. From Fig. 5.2, we can see that generally speaking our HFEAST algorithm has almost the same convergence behaviour with the BFEAST algorithm, they converge linearly before attain the accuracy the methods can achieve. The two non-Hermitian FEAST algorithms can achieve almost the same accuracy, although our method is slightly faster than its counterpart. Remarkably, we can see that our HFEAST is more stable than the BFEAST algorithm, for example, from Problem 4 and Problem 5. Here we should point out that the dominant computational costs at each iteration are almost the same for the two non-Hermitian FEAST algorithms, that is solving q linear systems of the form (3.14). 5.3. Comparison with Matlab function eig In this experiment, we compare our HFEAST algorithm with Matlab built-in function eig in terms of timing. Note that the target eigenvalues are the interior ones of non-Hermitian problems, when using eig to compute the eigenvalues inside the regions presented in Table 5.1, we have to first compute all eigenvalues in dense format and then

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

92

Table 5.2 Comparison of eig and HFEAST in terms of timing (in seconds). No.

eig

HFEAST

1 2 3 4 5 6 7 8 9

1.60 11.05 14.43 69.42 2409.82 2700.68 >10000 >10000 >10000

3.01 2.29 2.95 13.96 136.82 19.41 35.86 53.27 160.30

select the target eigenvalues according to the values of their coordinates. We set the convergence tolerance as 1.0 × 10−8 and take maximum iteration count max_iter = 10 for the HFEAST algorithm. The amount of time, which is measured in seconds, required by eig and HFEAST is reported in Table 5.2. It is clear to see that HFEAST is much faster than Matlab function eig, although the parallelism offered by HFEAST is not used in the tests. The difference in CPU times is more obvious when the size of test problem grows larger. 6. Conclusions In this work, we developed a new non-Hermitian FEAST algorithm. Our motivation was to overcome the difficulties of the potential ill-conditioning of eigenbasis and to formulate a more stable non-Hermitian FEAST algorithm. To this end, our HFEAST algorithm constructs the subspace spanned by (right) generalized Schur vectors, rather than the eigenspace, corresponding to the desired eigenpairs by the contour-integral technique. Our new method can also be used to compute the right generalized Schur vectors and left generalized Schur vectors. The convergence analysis of the proposed method will be our future work. Declaration of Competing Interest There is no competing interest. Acknowledgements I would like to thank Professor Raymond H. Chan, my thesis advisor, at The Chinese University of Hong Kong and Professor Man-Chung Yeung at University of Wyoming for their help and fruitful discussions in the preparation of this paper. We are grateful to the referees for many constructive comments and suggestions that led to substantial improvement of this paper. This research is supported by National Natural Science Foundation of China under grant number 11701593.

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

93

References [1] L. Ahlfors, Complex Analysis, 3rd edition, McGraw-Hill, Inc., 1979. [2] A.P. Austin, L.N. Trefethen, Computing eigenvalues of real symmetric matrices with rational filters in real arithmetic, SIAM J. Sci. Comput. 37 (2015) A1365–A1387. [3] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, H. Van Der Vorst, Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide, SIAM, Philadelphia, 2000. [4] Z. Bai, J. Demmel, M. Gu, An inverse free parallel spectral divide and conquer algorithm for nonsymmetric eigenproblem, Numer. Math. 76 (1997) 279–308. [5] W.-J. Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra Appl. 436 (2012) 3839–3863. [6] K.A. Cliffe, A. Spence, S.J. Tavener, The numerical analysis of bifurcation problems with application to fluid mechanics, Acta Numer. 9 (2000) 39–131. [7] P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, 2nd edition, Academic Press, Orlando, FL, 1984. [8] J. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997. [9] E. Di Napoli, E. Polizzi, Y. Saad, Efficient estimation of eigenvalue counts in an interval, http:// arxiv.org/abs/1308.4275. [10] H.-R. Fang, Y. Saad, A filtered Lanczos procedure for extreme and interior eigenvalue problems, SIAM J. Sci. Comput. 34 (2012) A2220–A2246. [11] D.R. Fokkema, G.L.G. Sleijpen, H.A. Van Der Vorst, Jacobi-Davidson style QR and QZ algorithms for the reduction of matrix pencils, SIAM J. Sci. Comput. 20 (1998) 94–125. [12] B. Ford, G. Hall, The generalized eigenvalue problem in quantum chemistry, Comput. Phys. Commun. 8 (1974) 337–348. [13] F.R. Gantmacher, The Theory of Matrices, Chelsea, New York, 1959. [14] G.H. Golub, C.F. Van Loan, Matrix Computations, 3rd edition, Johns Hopkins University Press, Baltimore, MD, 1996. [15] R.G. Grimes, J.D. Lewis, H.D. Simon, A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems, SIAM J. Matrix Anal. Appl. 15 (1994) 228–272. [16] N. Halko, P.G. Martinsson, J.A. Tropp, Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions, SIAM Rev. 53 (2011) 217–288. [17] T. Ikegami, T. Sakurai, Contour integral eigensolver for non-Hermitian systems: a Rayleigh-Ritztype approach, Taiwanese J. Math. 14 (2010) 825–837. [18] T. Ikegami, T. Sakurai, U. Nagashima, A filter diagonalization for generalized eigenvalue problems based on the Sakurai-Sugiura projection method, J. Comput. Appl. Math. 233 (2010) 1927–1936. [19] J. Kestyn, E. Polizzi, P. Tang, FEAST eigensolver for non-Hermitian problems, http://arxiv.org/ abs/1506.04463. [20] L. Krämer, E. Di Napoli, M. Galgon, B. Lang, P. Bientinesi, Dissecting the FEAST algorithm for generalized eigenproblems, J. Comput. Appl. Math. 244 (2013) 1–9. [21] S.E. Laux, Solving complex band structure problems with the FEAST eigenvalue algorithm, Phys. Rev. B 86 (2012) 075103. [22] C.B. Moler, G.W. Stewart, An algorithm for generalized matrix eigenvalue problems, SIAM J. Numer. Anal. 10 (1973) 241–256. [23] Y. Nakatsukasa, N.J. Higham, Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue decomposition and the SVD, SIAM J. Sci. Comput. 35 (2013) A1325–A1349. [24] E. Polizzi, Density-matrix-based algorithm for solving eigenvalue problems, Phys. Rev. B 79 (2009) 115112. [25] W.P. Reinhardt, Complex coordinates in the theory of atomic and molecular structure and dynamics, Annu. Rev. Phys. Chem. 33 (1982) 223–255. [26] A. Ruhe, Rational Krylov: a practical algorithm for large sparse nonsymmetric matrix pencils, SIAM J. Sci. Comput. 19 (1998) 1535–1551. [27] Y. Saad, Numerical Methods for Large Eigenvalue Problems, SIAM, Philadelphia, 2011. [28] T. Sakurai, Y. Futamura, H. Tadano, Efficient parameter estimation and implementation of a contour integral-based eigensolver, J. Algorithms Comput. Technol. 7 (2013) 249–269. [29] T. Sakurai, H. Sugiura, A projection method for generalized eigenvalue problems using numerical integration, J. Comput. Appl. Math. 159 (2003) 119–128. [30] T. Sakurai, H. Tadano, CIRR: a Rayleigh–Ritz type method with contour integral for generalized eigenvalue problems, Hokkaido Math. J. 36 (2007) 745–757.

94

G. Yin / Linear Algebra and its Applications 578 (2019) 75–94

[31] B.K. Sriperumbudur, D.A. Torres, G.R.G. Lanckriet, A majorization-minimization approach to the sparse generalized eigenvalue problem, Mach. Learn. 85 (2011) 3–39. [32] G.W. Stewart, Matrix Algorithms, Vol. II, Eigensystems, SIAM, Philadelphia, 2001. [33] P.T.P. Tang, E. Polizzi, FEAST as a subspace iteration eigensolver accelerated by approximate spectral projection, SIAM J. Matrix Anal. Appl. 35 (2014) 354–390. [34] E. Vecharynski, C. Yang, F. Xue, Generalized preconditioned locally harmonic residual method for non-Hermitian eigenproblems, SIAM J. Sci. Comput. 38 (2016) A500–A527. [35] L.N. Trefethen, M. Embree, Spectra and Pseudospectra, Princeton University Press, Princeton, NJ, 2005. [36] Y. Xi, Y. Saad, Computing partial spectra with least-squares rational filters, SIAM J. Sci. Comput. 38 (2016) A3020–A3045. [37] G. Yin, R.H. Chan, M-C. Yeung, A FEAST algorithm with oblique projection for generalized eigenvalue problems, Numer. Linear Algebra Appl. 24 (2017) e2092.