Generalized preconditioned locally harmonic residual method for non-Hermitian eigenproblems Version generated on July 4, 2019. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, SciDAC program (the first and second authors) and National Science Foundation under Grant DMS-1419100 (the third author).

Generalized preconditioned locally harmonic residual method for non-Hermitian eigenproblems thanks: Version generated on July 4, 2019. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, SciDAC program (the first and second authors) and National Science Foundation under Grant DMS-1419100 (the third author).

Eugene Vecharynski    Chao Yang Computational Research Division, Lawrence Berkeley National Laboratory, One Cyclotron Road, MS 50F-1620L, Berkeley, CA 94720, USA ({evecharynski,cyang}@lbl.gov)    Fei Xue Department of Mathematics, University of Louisiana at Lafayette, P.O. Box 41010, Lafayette, LA 70504-1010, USA (fxue@louisiana.edu)
Abstract

We introduce the Generalized Preconditioned Locally Harmonic Residual (GPLHR) method for solving standard and generalized non-Hermitian eigenproblems. The method is particularly useful for computing a subset of eigenvalues, and their eigen- or Schur vectors, closest to a given shift. The proposed method is based on block iterations and can take advantage of a preconditioner if it is available. It does not need to perform exact shift-and-invert transformation. Standard and generalized eigenproblems are handled in a unified framework. Our numerical experiments demonstrate that GPLHR is generally more robust and efficient than existing methods, especially if the available memory is limited.

Key words. Eigenvalue, eigenvector, non-Hermitian, preconditioned eigensolver

1 Introduction

Large non-Hermitian eigenproblems arise in a variety of important applications, including resonant state calculation [3, 16, 33] or excited state analysis for equation-of-motion coupled-cluster (EOM-CC) approaches [14, 15] in quantum chemistry, linear stability analysis of the Navier-Stokes equation in fluid dynamics [5, 6], crystal growth simulation [21, 22], magnetohydrodynamics [18], power systems design [17], and many others; see, e.g., [45] for more examples.

In their most general form, these problems can be written as

\hb@xt@.01(1.1)

where and are general square matrices. We are particularly interested in the case in which both and are very large and sparse, or available only implicitly through a matrix-vector multiplication procedure. If is the identity matrix, then (LABEL:eq:gevp) becomes a standard eigenproblem.

The spectrum of (LABEL:eq:gevp), denoted by , is given by a set of numbers that make singular. The value of can be infinity in the case of a singular . Given a scalar , which we refer to as a shift, we seek to find a subset of eigenvalues that are closest to , and their associated (right) eigenvectors . These eigenvalues can be either extreme eigenvalues (e.g., eigenvalues with the largest magnitude) or interior eigenvalues that are inside the convex hull of the spectrum.

Our focus in this paper is on algorithms for computing interior eigenvalues and their corresponding eigenvectors of a non-Hermitian pencil. It is well known that these eigenpairs are often difficult to compute in practice. Traditional methods, such as the inverse subspace iteration or variants of the shift-invert Arnoldi algorithm, see, e.g., [2, 36], rely on using spectral transformations that require performing LU factorizations of and computing solutions to triangular linear systems. Such an approach can be prohibitively expensive, especially when the problem size is large and the LU factors of are not sparse.

There has been some effort in recent years to develop methods that are factorization free. Examples of such methods include the inexact inverse subspace iteration [34, 48] and inexact shift-invert Arnoldi methods [8, 49] in which linear systems of the form are solved iteratively. While these schemes can be considerably less expensive per iteration, the overall convergence of these methods is less predictable. There is often a tradeoff between the accuracy of the solution to linear systems and the convergence rate of the subspace or Arnoldi iterations. Setting an appropriate convergence criterion for an iterative solution of is not straightforward.

Another class of factorization-free methods include the Generalized Davidson (GD) method [27, 28] and the Jacobi–Davidson (JD) methods [7]. The GD methods generally do not rely on solution of linear systems. The JD style QR and QZ algorithms (JDQR and JDQZ) presented in [7] can be viewed as preconditioned eigensolvers in which a Newton-like correction equation is solved approximately by a preconditioned iterative solver.

The GD method can be easily extended into a block method in which several eigenpairs can be approximated simultaneously [38, 50]. Block GD is widely used in quantum chemistry. A block GD method is particularly well suited for modern high performance computers with a large number of processing units. This is because in a block method several (sparse) matrix vector multiplications, which often constitute the major cost of the algorithm, can be carried out simultaneously. Furthermore, one can easily implement data blocking, exploit data reuse and take advantage of BLAS3 in a block method. These techniques can lead to significant speedups on modern high performance computers as reported in [1]. However, the convergence of the block GD method depends to a large degree on the effectiveness of the preconditioner. In quantum chemistry applications, the matrix is often diagonal dominant. The diagonal part of the matrix can be used to construct an effective preconditioner. In other applications, it is less clear how to construct a good preconditioner when is not diagonal dominant [46].

The JDQR and JDQZ methods tend to be more robust with respect to the choice of preconditioners. However, extending JDQR and JDQZ to block methods is not straightforward. Because JDQR and JDQZ methods compute one eigenpair at a time, concurrency can only be exploited within a single matrix-vector multiplication procedure, which limits its parallel scalability to a large number of processors.

In this paper, we present a new family of methods for solving non-Hermitian eigenproblems, called the Generalized Preconditioned Locally Harmonic Residual (GPLHR) methods. The proposed scheme is a block method that constructs approximations to a partial (generalized) Schur form of the matrix pair associated with (LABEL:eq:gevp) iteratively. It does not require performing an LU factorization of . It allows a preconditioner to be used to construct a search space from which approximations to the Schur vectors can be extracted. We demonstrate that the GPLHR methods are generally more robust and efficient than JDQR/JDQZ as well as the block GD method when the same preconditioner is used, and when dimension of the search spaces are comparable. In addition, because GPLHR is a block method, it is well suited for high performance computers that consist of many cores or processors.

The GPLHR algorithm is a generalization of the recently introduced Preconditioned Locally Harmonic Residual (PLHR) method for computing interior eigenpairs of Hermitian eigenproblems [46] to the non-Hermitian case (hence the name). The GPLHR scheme can also be viewed as an extension of the well known LOBPCG [19] method and block inverse-free preconditioned (BIFP) Krylov subspace methods [9, 32] for computing extreme eigenpairs of Hermitian problems. At each step, all these methods construct trial subspaces of a fixed dimension and use them to extract the desired approximations by means of a subspace projection. The key difference is that GPLHR systematically uses a harmonic Rayleigh–Ritz procedure [26, 29], which allows for a robust eigen- or Schur pair extraction independent of the location of the targeted eigenvalues. Moreover, the proposed method explicitly utilizes properly chosen projectors to stabilize convergence, and is capable of performing iterations using the approximations to Schur vectors instead of potentially ill-conditioned eigenvectors.

This paper is organized as follows. In Section LABEL:sec:eig, we review standard eigenvalue-revealing decompositions for non-Hermitian eigenproblems and propose a new generalized Schur form which emphasizes only the right Schur vectors. This form is then used to derive the GPLHR algorithms for standard and generalized eigenproblems in Section LABEL:sec:pqr. The preconditioning options are discussed in Section LABEL:sec:prec. Section LABEL:sec:num presents numerical results. Our conclusions can be found in Section LABEL:sec:con. Appendix A briefly describes the approximate eigenbasis based GPLHR iteration (GPLHR-EIG).

2 Eigenvalue-revealing decompositions

For non-Hermitian eigenproblems, it is common (see, e.g., [24, 36, 44]) to replace (LABEL:eq:gevp) by an equivalent problem

\hb@xt@.01(2.1)

where the pair , called a generalized eigenvalue, defines an eigenvalue of the matrix pair . This formulation is advantageous from the numerical perspective. For example, it allows revealing the infinite or indeterminate eigenvalues of (LABEL:eq:gevp). The former corresponds to the case where is singular and , , whereas the latter is indicated by and takes place if both and are singular with a non-trivial intersection of their null spaces. Additionally, the separate treatment of and allows one to handle situations where either of the quantities is close to zero, which leads to underflow or overflow for . Note that a generalized eigenvalue is invariant with respect to multiplication by a scalar, i.e., for any nonzero , is also a generalized eigenvalue. In order to fix a representative pair, a normalizing condition should be imposed, e.g., .

In this paper, we assume that the pair is regular, i.e., is not identically zero for all and . The violation of this assumption gives a singular pair, which corresponds to an ill-posed problem. In this case, a regularization procedure should be applied to obtain a nearby regular pair, e.g., [2, Chapter 8.7], for which the eigenproblem is solved. In most cases, however, is regular, even when or or both are singular as long as the intersection of their null spaces is zero. Thus, regular pairs can have infinite eigenvalues, whereas the possibility of indeterminate eigenvalues is excluded.

2.1 Partial eigenvalue decomposition and generalized Schur form

Let us assume that the targeted eigenvalues of (LABEL:eq:gevp) are non-deficient, i.e., their algebraic and geometric multiplicities coincide. Then, given formulation (LABEL:eq:gevp_mod), the partial eigenvalue decomposition of can be written in the form

\hb@xt@.01(2.2)

where is a matrix whose columns correspond to the eigenvectors associated with the wanted eigenvalues . The -by- matrices and are diagonal, with the diagonal entries given by and , respectively. In our context, it is assumed that are the eigenvalues closest to .

Since, by our assumption, ’s are non-deficient, there exists associated linearly independent eigenvectors, so that is full-rank and (LABEL:eq:partial) is well defined. For a vast majority of applications, the partial eigenvalue decomposition (LABEL:eq:partial) does exist, i.e., the desired eigenvalues are non-deficient. Nevertheless, pursuing decomposition (LABEL:eq:partial) is generally not recommended, as it aims at computing a basis of eigenvectors which can potentially be ill-conditioned.

In order to circumvent the difficulties related to ill-conditioning or deficiency of the eigenbasis, we consider instead a partial generalized Schur decomposition

\hb@xt@.01(2.3)

see, e.g., [7, 36, 43]. Here, and are the -by- matrices of the right and left generalized Schur vectors, respectively. The factors and are -by- upper triangular matrices, such that the ratios of their diagonal entries correspond to the desired eigenvalues.

The advantage of the partial Schur form (LABEL:eq:gschur) is that it is defined for any and both and are orthonormal, so that they can be computed in a numerically stable manner. The matrix then gives an orthonormal basis of the invariant subspace associated with the eigenvalues of interest. If individual eigenvectors are needed, the right generalized Schur vectors can be transformed into eigenvector block by the means of the standard Rayleigh–Ritz procedure for the pair , performed over the subspace spanned by the columns of . Similarly,  contains an orthonormal basis of the left invariant subspace and can be used to retrieve the left eigenvectors.

2.2 The “-free” partial generalized Schur form

We start by addressing the question of whether it is possible to eliminate from decomposition (LABEL:eq:gschur) and obtain an eigenvalue-revealing Schur type form that contains only the right Schur vectors and triangular factors. Such a form would provide a numerically stable alternative to the eigenvalue decomposition (LABEL:eq:partial). It would also give an opportunity to handle the standard and generalized eigenproblems in a uniform manner. Most importantly, it would allow us to define a numerical scheme that emphasizes computation only of the generalized right Schur vectors, which are often the ones needed in a real application. The left Schur basis can be obtained as a by-product of computation.

For example, if is nonsingular, then can be eliminated in a straightforward manner by multiplying the bottom equation in (LABEL:eq:gschur) by and substituting the resulting expression for into the top equation. This yields the desirable decomposition , where is triangular with eigenvalues on the diagonal. Similarly, assuming that is nonsingular, one can obtain , where is also triangular, with inverted eigenvalues on the diagonal.

However, the above treatment involves inversions of or , and consequently could cause numerical instabilities whenever the triangular matrices to be inverted have very small diagonal entries. Moreover, in the presence of both zero and infinite eigenvalues, and are singular and cannot be inverted at all. Therefore, it would be ideal to construct a -free Schur form for any regular that does not rely on the inverse of the triangular matrices and .

In the rest of the section, we show that such a “-free” partial Schur form is possible. We start with the following lemma.

Lemma 2.1

Let and be upper triangular matrices, such that for all . Then for any upper triangular and , and any , such that

\hb@xt@.01(2.4)

and

\hb@xt@.01(2.5)

the matrix is upper triangular with for all .

Proof. The matrix is a combination of upper triangular matrices and, hence, is upper triangular. For each diagonal entry of , we have , where and do not equal zero simultaneously. It is then easy to check that for any , (LABEL:eq:G1) and (LABEL:eq:G2) satisfy for all .     

The following theorem gives the existence of a “-free” generalized Schur form.

Theorem 2.2

For any regular pair there exists a matrix with orthonormal columns and upper triangular matrices , such that

\hb@xt@.01(2.6)

and are (possibly infinite) eigenvalues of .

Proof. We start from the generalized Schur decomposition (LABEL:eq:gschur). Let and be upper triangular with diagonal elements defined by (LABEL:eq:G1) and (LABEL:eq:G2). Right-multiplying both sides of the top and bottom equalities in (LABEL:eq:gschur) by and , respectively, and then summing up the results, gives an equivalent system

\hb@xt@.01(2.7)

where

Since is a regular pair, and do not equal zero simultaneously for the same . Thus, by Lemma LABEL:lem:triu, is upper triangular with for all . In particular, the latter implies that the inverse of exists. Hence, we can eliminate from (LABEL:eq:gschur_tmp).

It follows from the bottom equality in (LABEL:eq:gschur_tmp) that . Substituting it to the top identity gives , which implies (LABEL:eq:gschur2) with

\hb@xt@.01(2.8)

Furthermore, since the diagonal entries of are all ones, for all . Thus, from (LABEL:eq:M), the diagonal elements of and are related to those of and by and . In particular, using definition (LABEL:eq:G1) and (LABEL:eq:G2) of and , if , we obtain

\hb@xt@.01(2.9)

which implies that , provided that is chosen such that .

Similarly, if ,

\hb@xt@.01(2.10)

Hence, by choosing , we get .     

Given the traditional generalized Schur decomposition (LABEL:eq:gschur), Theorem LABEL:thm:gschur2 suggests a simple approach for obtaining the “Q-free” form (LABEL:eq:gschur2). Namely, one starts with setting up the upper triangular matrices , , and defined according to Lemma LABEL:lem:triu, and then applies (LABEL:eq:M) to evaluate and . The right Schur basis is the same in (LABEL:eq:gschur) and (LABEL:eq:gschur2). Note that the definition (LABEL:eq:M) of the factors and  does not rely on inverting or and only requires an inverse of a triangular matrix which has all ones on the diagonal. Hence, it avoids inverting (nearly) singular triangular matrices. Additionally, the proposed formulas (LABEL:eq:G1) and (LABEL:eq:G2), as well as (LABEL:eq:Mj1), always have the largest modulus number between and in the denominator, which mitigates potential numerical issues introduced by the small diagonal elements of and .

Clearly, the choice of and in (LABEL:eq:G1) and (LABEL:eq:G2), and hence decomposition (LABEL:eq:gschur2), is not unique, as it depends on the value of the parameter and allows freedom in the choice of entries above diagonal in and . Therefore, in order to fix particular and , in practice, we choose and to be diagonal and set if and , otherwise. This yields the diagonal matrices

\hb@xt@.01(2.11)

and

\hb@xt@.01(2.12)

for which, by (LABEL:eq:Mj1), , if , and, by (LABEL:eq:Mj2), , , otherwise. Note that it is also possible to choose that give and , such that , but our implementation follows the simple formulas (LABEL:eq:G1_fix) and (LABEL:eq:G2_fix).

3 The GPLHR algorithm for computing a partial Schur form

For Hermitian problems, a class of powerful eigensolvers is based on preconditioned block iterations; e.g., [19, 32, 46]. Such methods fit into a unified framework consisting of two key ingredients: generation of the trial (or search) subspace and extraction of the approximate eigenvectors. Namely, given a number of eigenvector approximations and a preconditioner, at each step , a low-dimensional trial subspace is constructed. If properly chosen, contains improved approximations to the wanted eigenvectors, which are extracted from the subspace using a suitable projection procedure and are then used to start a new iteration. In this section, we extend this framework to the non-Hermitian case.

The results of Section LABEL:subsec:gschur_new are central to our derivations. In particular, instead of the standard partial generalized Schur decomposition (LABEL:eq:gschur), we focus on the “-free” form (LABEL:eq:gschur2), and seek approximation of the right Schur vectors and the associated triangular factors and . As we shall see, the approximation to the left Schur vectors , as well as to triangular matrices and , can be easily obtained as a by product of the proposed scheme.

3.1 Construction of the trial subspace

Let , , and be approximations to , , and in (LABEL:subsec:gschur_new) at an iteration , such that has orthonormal columns, and and are -by- upper triangular. In order to define a trial subspace that can be used for searching a new approximation , we adopt the following subspace orthogonal correction viewpoint.

We are interested in constructing a subspace that provides a good representation of the correction , such that

\hb@xt@.01(3.1)

where ; and , are corrections of the triangular factors and , respectively. Once the “correction subspace” is determined, the trial subspace can be defined as the subspace sum , where denotes the column space of .

After rearranging the terms in (LABEL:eq:schur_cor), we get

where is the matrix of the exact right Schur vectors. Then, using (LABEL:eq:gschur), we arrive at the identity

\hb@xt@.01(3.2)

where is the matrix of the exact left Schur vectors. Let be an orthogonal projector onto , where is a matrix whose orthonormal columns represent a basis of for some scalars , such that . Then, applying to both sides of (LABEL:eq:schur_cor1) and neglecting the high-order term gives the equation

\hb@xt@.01(3.3)

whose solution , constrained to satisfy , provides an approximation to the desired correction of (LABEL:eq:schur_cor). Here, is an orthogonal projector onto and, hence, .

Equation (LABEL:eq:sylvester) represents a generalized Sylvester equation [39] and can be viewed as a block extension of the standard single vector Jacobi–Davidson (JD) correction equation [7, 40], where the right-hand side is given by the projected residual of problem (LABEL:eq:gschur2). Note that a connection between the subspace correction and solution of a Sylvester equation was also observed in [31], though not in the context of the generalized Schur form computation. Thus, a possible option is to (approximately) solve (LABEL:eq:sylvester) for a block correction and then set .

However, with this approach, it is not straightforward to solve the matrix equation (LABEL:eq:sylvester) efficiently. Under certain assumptions, one can reduce (LABEL:eq:sylvester) to the standard Sylvester equation, but this is generally expensive for large problems.

Instead, we do not seek the solution of the generalized Sylvester equation, and use a different strategy to generate the correction subspace . To this end, we consider (LABEL:eq:sylvester) as an equation of the general form

\hb@xt@.01(3.4)

where and . Since is a linear operator, standard subspace projection techniques for solving linear systems, e.g., [10, 35], can be formally applied to (LABEL:eq:operator). Thus, we can define a preconditioned Krylov-like subspace for (LABEL:eq:operator). However, instead of using this subspace to solve (LABEL:eq:operator), we place it directly in an eigensolver’s trial subspace.

More precisely, let be a preconditioner for the operator . Then, following the analogy with the Krylov subspace methods for (block) linear systems [10, 11, 35], one can expect that the solution of (LABEL:eq:operator) can be well represented in the preconditioned block Krylov subspace

generated by the operator and the starting block , which contains all blocks of the form , with  [11]. In particular, this implies that each column of should be searched in

\hb@xt@.01(3.5)

where the blocks represent what we call the Krylov–Arnoldi sequence, generated by a preconditioned Arnoldi-like procedure for problem (LABEL:eq:operator), or equivalently, for the generalized Sylvester equation (LABEL:eq:sylvester). This procedure is stated in Algorithm LABEL:alg:barnoldi, where the input parameter determines the subspace size.

Input: preconditioning operator , the parameter . Output: the Krylov–Arnoldi basis 1:  ; ; 2:  for  do 3:     ; ; 4:     for  do 5:         ; 6:     end for 7:     ; 8:  end for 9:  Return .
Algorithm 1 Preconditioned block Arnoldi type procedure for equation (LABEL:eq:operator)

Thus, given the subspace (LABEL:eq:krylov_sep_block), capturing the orthogonal correction , the eigensolver’s trial subspace can be defined as , i.e.,

\hb@xt@.01(3.6)

Clearly, the Krylov–Arnoldi sequence vectors , generated by Algorithm LABEL:alg:barnoldi, represents a system of orthonormal columns that are orthogonal to due to the choice of the preconditioner . Hence, the trial subspace in (LABEL:eq:trial_partial) is spanned by orthonormal vectors, which gives a stable basis.

A desirable feature of the preconditioner for problem (LABEL:eq:operator) (or (LABEL:eq:sylvester)) is that it should approximate the inverse of the operator . One way to construct such an approximation is to replace the upper triangular matrices and in the expression for by diagonals and , respectively, where . This results in an approximation of that is given by the matrix . Thus, one can choose the preconditioner as

\hb@xt@.01(3.7)

where . Throughout, we consider as a preconditioner for eigenvalue problem (LABEL:eq:gevp), which should be distinguished from the preconditioner for the generalized Sylvester equation. Note that, in practice, assuming that is a finite shift that is different from an eigenvalue, we can set and , so that is an approximate shift-and-invert operator.

Finally, given the definition (LABEL:eq:prec_gen) of the preconditioner , we can derive explicit expressions for the blocks in the Krylov–Arnoldi sequence generated from Algorithm LABEL:alg:barnoldi. It is easy to check that

\hb@xt@.01(3.8)

and, for ,

\hb@xt@.01(3.9)

where and is the orthogonal projector onto ,

\hb@xt@.01(3.10)

We remark that the above construction of the trial subspace (LABEL:eq:trial_partial) is similar in spirit to the one used to devise the BIFP Krylov subspace methods in [9, 32] for computing extreme eigenpairs of Hermitian problems. The conceptual difference, however, is that we put the block correction equation (LABEL:eq:sylvester) in the center stage, whereas the BIFP methods are based on Krylov(-like) subspaces delivered by a simultaneous solution of problems for different shifts . In particular, this allows us to discover the projectors and in (LABEL:eq:prec_gen). These projectors are then blended into formulas (LABEL:eq:WA_gschur)–(LABEL:eq:projS) that define the trial subspace. Their presence has a strong effect on eigensolver’s robustness, as we shall demonstrate in the numerical experiments of Section LABEL:sec:num.

3.2 The trial subspace for standard eigenvalue problem

If , instead of the generalized form (LABEL:eq:gschur), we seek the standard partial Schur decomposition

\hb@xt@.01(3.11)

where contains orthonormal Schur vectors and is upper triangular with the wanted eigenvalues of on its diagonal. It is clear that (LABEL:eq:schur) is a special case of the “-free” form (LABEL:eq:gschur2) with , , and . Therefore, the derivation of the trial subspace, described in the previous section, is directly applicable here.

In particular, let be an approximate Schur basis, and let and be the associated upper triangular matrices that approximate and , respectively. Then, following the arguments of Section LABEL:subsec:trial_gschur, from (LABEL:eq:schur_cor) with , we get

\hb@xt@.01(3.12)

This equality is the analogue to (LABEL:eq:schur_cor1) for a standard eigenvalue problem. Here, in order to approximate the correction , it is natural to apply the projector to both sides of (LABEL:eq:schur_cor1_stand) and neglect the term . As a result, we obtain the equation

\hb@xt@.01(3.13)

where the solution , which is orthogonal to , approximates the desired . Note that in contrast to the correction equation (LABEL:eq:sylvester) for the generalized eigenvalue problem, which is based on two projectors and , (LABEL:eq:sylvester_st) involves only one projector . This is expected, as both and approximate the same vectors, .

Considering (LABEL:eq:sylvester_st) as an equation of the general form (LABEL:eq:operator), where and , we apply steps of Algorithm LABEL:alg:barnoldi with the preconditioner

\hb@xt@.01(3.14)

This yields the trial subspace (LABEL:eq:trial_partial), where

\hb@xt@.01(3.15)

and, for ,

\hb@xt@.01(3.16)

As before, we assume that in (LABEL:eq:SA_schur) and that is the orthogonal projector defined in (LABEL:eq:projS); the preconditioner in (LABEL:eq:prec) is chosen as an approximation of the shift-and-invert operator .

3.3 The harmonic Schur–Rayleigh–Ritz procedure

Let be a trial subspace of dimension at iteration defined by (LABEL:eq:trial_partial) with (LABEL:eq:WA_gschur)–(LABEL:eq:projS). We try to find a new orthonormal approximate Schur basis and the corresponding -by- upper triangular matrices and , such that each column of belongs to  and approximate the triangular factors in (LABEL:eq:gschur2).

To fulfill this task, we use a harmonic Rayleigh–Ritz projection [26, 29], adapted to the case of the “Q-free” Schur form (LABEL:eq:gschur2). Specifically, we let be a test subspace, where we assume that the target shift is not an eigenvalue. Then the approximations , , and can be determined by requiring each column of the residual of (LABEL:eq:gschur2) to be orthogonal to this test subspace, i.e.,

\hb@xt@.01(3.17)

This yields the projected problem

\hb@xt@.01(3.18)

where , are unknown; and contain orthonormal columns spanning and , respectively. Once we have computed , , and , approximations to the Schur vectors, determined by (LABEL:eq:pg), are given by .

Let us now consider the solution of the projected problem (LABEL:eq:proj_schur). We first show that the projected matrix pair is regular.

Proposition 3.1

Let contain orthonormal basis vectors of the trial and test subspaces, and assume that is different from any eigenvalue of . Then the projected pair is regular.

Proof. Assume, on the contrary, that is singular. Then, by definition of a singular matrix pair, regardless of the choice of ,

\hb@xt@.01(3.19)

Since is an orthonormal basis of , we can write , where is an -by- square matrix. This matrix is nonsingular, i.e., , because is nonsingular. Hence, , and

Thus, from the above equality and (LABEL:eq:det0), we have

Since , this identity implies that . But the matrix is Hermitian positive semidefinite and its singularity implies that , which contradicts the assumption that .     

Thus, problem (LABEL:eq:proj_schur) is of the form (LABEL:eq:gschur2) and the matrix pair is regular, whenever is not an eigenvalue of . Therefore, (LABEL:eq:proj_schur) can be solved by the approach described in Section LABEL:subsec:gschur_new.

Specifically, one first employs the standard sorted QZ algorithm [24] to compute the full generalized Schur form of . This will produce two unitary matrices of the left and right generalized Schur vectors, along with the upper triangular , ordered in such a way that the ratios of the first diagonal elements are closest to , and the corresponding Schur vectors are located at the leading columns of and . Then, we let , and obtain the desired Schur basis in (LABEL:eq:proj_schur) as . The triangular factors are given as and , where , and , and , are diagonal matrices defined by (LABEL:eq:G1_fix) and (LABEL:eq:G2_fix) with , replaced with , , respectively.

As a result, the described extraction approach, which we call a harmonic Schur–Rayleigh–Ritz (SRR) procedure, constructs a new basis of the approximate Schur vectors (the harmonic Schur–Ritz vectors) and the upper triangular matrices , , such that the ratios of their diagonal entries are approximations to the desired eigenvalues. Note that in the case of a standard eigenvalue problem, the proposed use of formulas (LABEL:eq:G1_fix) and (LABEL:eq:G2_fix) for evaluating the matrices and indeed guarantees that they converge to and , as and get closer to and , respectively.

Although both the construction of the trial subspace and the presented harmonic SRR procedure aim at finding the factors , , and of the “Q-free” form (LABEL:eq:gschur2), it is easy to see that, as a by-product, the scheme can also produce approximations to the left Schur vectors and the upper triangular factors , of the conventional generalized Schur decomposition (LABEL:eq:gschur) if . Specifically, the former is given by , whereas the latter correspond to and .

Note that the computed can be conveniently exploited at iteration to construct the new test subspace and, if , define the projector , because represents an orthonormal basis of .

Proposition 3.2

Let <