Computing several eigenvalues of nonlinear eigenvalue problems by selection

Computing more than one eigenvalue for (large sparse) one-parameter polynomial and general nonlinear eigenproblems, as well as for multiparameter linear and nonlinear eigenproblems, is a much harder task than for standard eigenvalue problems. We present simple but efficient selection methods based on divided differences to do this. Selection means that the approximate eigenpair is picked from candidate pairs that satisfy a certain suitable criterion. The goal of this procedure is to steer the process away from already detected pairs. In contrast to locking techniques, it is not necessary to keep converged eigenvectors in the search space, so that the entire search space may be devoted to new information. The selection techniques are applicable to many types of matrix eigenvalue problems; standard deflation is feasible only for linear one-parameter problems. The methods are easy to understand and implement. Although the use of divided differences is well known in the context of nonlinear eigenproblems, the proposed selection techniques are new for one-parameter problems. For multiparameter problems, we improve on and generalize our previous work. We also show how to use divided differences in the framework of homogeneous coordinates, which may be appropriate for generalized eigenvalue problems with infinite eigenvalues. While the approaches are valuable alternatives for one-parameter nonlinear eigenproblems, they seem the only option for multiparameter problems.


Introduction
In large sparse matrix eigenvalue problems, a common task is to compute a few eigenvalues closest to a given target, largest in magnitude, or rightmost in the complex plane.If we have already (approximately) computed a number of eigenpairs, and would like to compute a new pair, it is of importance to avoid convergence to one of the previously computed pairs.
Let A ∈ C n×n .For the standard eigenvalue problem avoidance of previous vectors can be conveniently achieved by computing Schur vectors instead of eigenvectors.This technique is based on the Schur decomposition AQ = QR for a matrix A, where Q is unitary and R is upper triangular.If (λ 1 , q 1 ), . . ., (λ d , q d ) are Schur pairs computed earlier in the process and has the same Schur pairs as A, except for λ 1 , . . ., λ d which are replaced by zero eigenvalues.In a subspace method, the search space may then be kept orthogonal to q 1 , . . ., q d so that the subspace method does not notice these zero eigenvalues.For the generalized eigenvalue problem (GEP) Ax = λBx, where B ∈ C n×n , the generalized Schur decomposition for matrix pencils may be exploited in a similar way.
The Jacobi-Davidson QR (JDQR [28]) method for (1) and Jacobi-Davidson QZ (JDQZ [28]) method for GEP are both examples of methods that are based on the described strategies.Modified deflation techniques are available for other types of linear eigenvalue problems such as the (generalized) singular value problem [8,9].
In this paper we discuss a new approach to find several eigenvalues for the regular nonlinear eigenvalue problem (NEP) where F (λ) is an n × n matrix, whose elements are analytic functions of the complex parameter λ.The regularity means that det(F (λ)) does not vanish identically.As a special case, we will consider the polynomial eigenvalue problem (PEP) where all matrices are n × n.In some applications, the leading matrix Am may be singular, and eigenvalues may be infinite.For this reason, it may be beneficial to consider homogeneous coordinates, which we will do in Section 3.Moreover, because of the practical importance, as well as for convenience of presentation, we will focus in particular on the quadratic eigenvalue problem (QEP) Already for the QEP no deflation procedure comparable to those for the standard and generalized eigenvalue problems is known, which is natural in view of the existence of 2n eigenpairs.This implies that if we compute eigenvalues with, for instance, the Jacobi-Davidson method [28], we may find the same eigenvalue again without preventive measures.It is possible to linearize the QEP or PEP to a GEP, for which a deflation procedure is possible.However, a clear drawback of this is that linearization increases the dimension of the problem by a factor m − 1 (see also the comments in Section 5).Also, the mathematical properties of linearizations are interesting but not straightforward; see, e.g., [10].For the NEP, linearizations are even more involved.It is therefore relevant to study techniques that can be directly applied to the problem at hand.Several strategies have been mentioned to compute several eigenvalues of nonlinear eigenproblems.One option to find more than one eigenvalue is locking, which was studied for the QEP by Meerbergen [21].The essence of this approach is to carry out no further computations on sufficiently converged Schur vectors, and retain them in the search space.While this method may be effective, a disadvantage is that the size of the search space grows steadily.In [5,6,18], a technique called nonequivalence deflation has been proposed, which replaces the original problem by another.Kressner [20] has developed a block method, while Effenberger [4] proposes a deflation strategy for nonlinear eigenproblems.Several methods have been discussed and compared in [4] and by Güttel and Tisseur [7].Some further comparisons can be found in Section 5.In this paper, we propose an alternative simple and elegant strategy: computing several eigenvalues by selection.In each iteration we pick an approximate eigenpair from (Ritz) pairs that meet a certain selection criterion.This new strategy is particularly simple to comprehend and implement.We present several selection methods to compute more than one eigenpair of linear and nonlinear, and one-parameter and multiparameter eigenvalue problems (MEPs).We present various variants for the QEP and PEP, and for linear and nonlinear multiparameter eigenvalue problems.The selection criteria can be used for all types of eigenvalue problems, provided that expressions for the divided difference and derivative of the problem are available; see (7).
This work contains contributions in three directions.Although we have already successfully used some of these selection techniques in our work on linear multi-parameter eigenvalue problems ( [11,15], followed by nonlinear two-parameter eigenproblems [13,25] and linear three-parameter eigenvalue problems [12] very recently), we will present an improvement on these criteria for these problems.Secondly, to the best of our knowledge, the use of selection techniques to compute several eigenvalues in the form as described in this paper is new for one-parameter nonlinear eigenproblems: the QEP (4), PEP (3), and general NEP (2).Thirdly, these problems may have infinite eigenvalues when the leading coefficient matrix is singular.Therefore, methods exploiting homogeneous coordinates may be attractive for these problems (cf., e.g., [3,14]).We therefore introduce divided differences in homogeneous coordinates, which is nontrivial, and new to the best of our knowledge.
Finally, this paper is also meant to serve as an overview paper on selection techniques.We hope that it will inform about, popularize, and facilitate the use of these effective and easy-to-implement methods for eigenvalue problems.
There are various subspace expansion methods for eigenvalue problems.In this paper we will focus on the Jacobi-Davidson method.However, we want to stress that there are several other options to perform a subspace expansion, such as nonlinear Arnoldi [31] for NEPs, and Krylov type methods for MEPs [22].The selection techniques operate independently of the expansion of the subspace.
The rest of this paper has been organized as follows.In Section 2 we introduce a new selection criterion for nonlinear one-parameter eigenvalue problems.Eigenvalue problems involving matrices which are not full rank may have infinite eigenvalues.To deal with these in a consistent framework, homogeneous coordinates are the proper viewpoint; this is studied in Section 3. Section 4 focuses on the use of selection criteria for linear and nonlinear MEPs, our original motivation to study these techniques.We end with some numerical experiments and conclusions in Sections 6 and 7.

Selection for nonlinear one-parameter eigenvalue problems
We first introduce some basic notation and facts.The pair (λ, x) is an eigenpair if F (λ)x = 0 for a nonzero vector x.If y * F (λ) = 0 for a nonzero y, then y is a left eigenvector for the eigenvalue λ.We assume that both x and y have unit norm.We say that λ 0 is a simple eigenvalue when f (λ) := det(F (λ)) has a simple zero at λ = λ 0 .Neumaier [24] proves the following proposition about the left and the right eigenvectors of a simple eigenvalue.The same result with an alternative proof is presented in [27].
The divided difference for F is defined as which is continuous in both variables λ and µ.An alternative way to denote this without distinction of cases is . We will use similar expressions for convenience and brevity in Section 4.
This divided difference enjoys the following two key properties that can be used in the selection process.First, it is easy to see that when x j is the right eigenvector for λ j , and y i is the left eigenvector for a different eigenvalue λ i = λ j .Secondly, when λ i is a simple eigenvalue, then it follows from Proposition 1 that y * i F [λ i , λ i ] x i = 0. Based on these two observations, we develop a new selection criterion for NEPs using this F [•, •]orthogonality of right and left eigenvectors.Suppose that we have already computed eigentriplets (λ 1 , x 1 , y 1 ), . . ., (λ d , x d , y d ) for (2).We assume that all computed eigenvalues λ 1 , . . ., λ d are simple; the problem is allowed to have multiple eigenvalues, as long as the computed ones are simple.Our selection criteria are not suitable for multiple eigenvalues.
Suppose that (θ, v) is a candidate approximation for the next eigenpair, where also v has unit norm.To steer convergence to a pair different from the previously detected eigenpairs, and in view of (6), we only consider approximate eigenpairs for which |y * i F [λ i , θ] v| is sufficiently small for i = 1, . . ., d.To be precise, in the selection of the candidate approximate eigenpairs we require that max i=1,...,d where 0 < η < 1 is a fixed constant, which controls the strictness of the selection.Note that the denominator is a well-known quantity that arises, e.g., in the eigenvalue condition number (cf.[7,Th. 2.20]).
The next proposition explains the behavior of the criterion close to an eigenpair.Suppose that we have already computed the eigenpair (λ 1 , x 1 ), and now approximate a pair (λ 2 , x 2 ).Let us assume that (λ 2 + εφ, x 2 + εw), for small ε and vectors of unit norm, is an approximation for (λ 2 , x 2 ).This is a realistic assumption in this situation; see [17] for options to extract an approximate eigenvalue from an approximate eigenvector for the QEP and PEP.
Proposition 2 Let y 1 = 0 be a left eigenvector for a simple eigenvalue λ 1 of the nonlinear eigenvalue problem F (λ)x = 0. Let x 2 = 0 be a right eigenvector for a simple eigenvalue λ 2 = λ 1 and let (θ, v), where θ = λ 2 + εφ and v = x 2 + εw, be a candidate for the next eigenpair.Then where Proof The proposition follows from the Taylor series expansion where we take into account that y * 1 F (λ 1 ) = 0 and F (λ 2 )x 2 = 0.
Proposition 2 indicates that selection based on divided differences may be difficult for eigenvalues which are very close.This is in line with the earlier remark that the selection methods are not suited for multiple eigenvalues.We will again see the expression for C in Section 3.
Under the assumptions of Proposition 2, x 2 ) and converges to 1 when (θ, v) → (λ 1 , x 1 ).Although this shows that the selection criteria (7) may work for any value η in the interval (0, 1), we generally recommend to use η ≈ 0.1, to avoid the already computed eigenvalues while simultaneously avoiding the selection process to be unnecessarily strict.We will use this value in the experiments in Section 6.
The divided difference for the PEP (3) has the following simple explicit expression.
Proposition 3 For the polynomial eigenvalue problem (3) we have In particular, the divided difference for the QEP (4) is Proof This follows from an easy calculation.
The selection criterion for a candidate eigenpair (θ, v) for the QEP using divided differences then becomes max i=1,...,d We note that an important key to the success and simplicity of the selection criteria is the fact that divided differences can be elegantly generalized for many types of eigenproblems.More specifically, we will describe a homogeneous variant of divided differences and a selection criteria in Section 3, and a multiparameter variant in Section 4.
We now proceed to discuss some practical matters of the proposed techniques.A main disadvantage of our selection criterion is that one needs the left eigenvectors corresponding to converged eigenvalues during the process.For symmetric and Hermitian eigenvalue problems and right-definite multiparameter eigenvalue problems (see Section 4), these left eigenvectors come without extra computations.In other cases, we generally need extra work to compute them.For large-scale problems, this may comprise several additional matrix-vector products to solve y from F (λ) * y = 0. Pseudocode for this task is given in Algorithm 1.
Algorithm 1: An iterative approach for computing a null vector of a singular matrix Input: (Almost) singular Z, (random) nonzero starting vector y 0 , tolerance ε.Output: An approximate null vector y of Z with Zy ≤ ε.
1: Compute b = Zy 0 / Zy 0 2: Solve approximately Zx = b with an iterative method, e.g., (preconditioned) GMRES with tolerance ε 3: Set y = (x − y 0 ) / x − y 0 We note that in some situations we may have a reasonable approximation to the left eigenvector y in the process.Moreover, and more importantly, any available preconditioner will usually be of great help.For instance, for eigencomputations, often a preconditioner M ≈ F (τ ), where τ ∈ C is the target of interest, is at our disposal.
For problems that are not truly large-scale, we can solve the system in step 2 in Algorithm 1 with a direct instead of an iterative method.In particular, this holds for applications of the selection techniques for multiparameter eigenvalue problems, where the problem size is often relatively small: here, p vectors of length n are sought while the total problem size is of size n p , where p is the number of parameters.
In either case, the extra work to compute the left eigenvectors will often be relatively small compared with the overall effort of computing the eigenpairs.In addition, the left eigenvectors are very useful information to determine the condition numbers of the eigenvalues, an important quantity to assess the reliability of the eigenvalue.For instance, for the QEP an absolute eigenvalue condition number is given by Tisseur [30]: where we assume that x and y are of unit norm.For general PEPs there is a straightforward analogous expression [30].We will study the extra costs of the selection in more detail later in this section.
As already mentioned in the introduction, the selection criteria presented in this paper can be combined with appropriate subspace methods, which expand given subspaces, and extract potential eigenpairs from the subspaces.As an example, we now give a pseudocode for a Jacobi-Davidson method to compute several eigenpairs for the polynomial eigenvalue problem in Algorithm 2 (cf.also [16]).
Algorithm 2: Jacobi-Davidson type method for computing several eigenpairs for the polynomial eigenvalue problem (3) using selection.
Input: Matrix polynomial P (λ) = m i=0 λ i A i , desired number of eigenpairs d, starting vector v, tolerance ε, minimum and maximum dimensions of search space mindim, maxdim, threshold η for (7).Output: d approximate eigenpairs. 1: Compute kth rows and columns of H

(possibly with preconditioner)
Let us comment on Algorithm 2. The key steps for the selection process are indicated in boldface.In step 3, rgs denotes repeated Gram-Schmidt orthogonalization or any other method to compute an orthonormal basis in a stable way.There are several strategies available for the extraction of an approximate eigenpair from the search space V k in step 5.For the extraction of approximate eigenvectors, standard and harmonic Rayleigh-Ritz are default options [16], while [17] discusses choices for an approximate eigenvalue.We can use the selection criterion independently of these extraction choices.The key selection approach described in this paper is used in steps 5 and 9.
In step 5, it might happen that none of the Ritz vectors satisfies the selection criterion.In such a case we pick the best Ritz pair given the target regardless of the selection criteria and skip step 8 so that the method does not find the same eigenpair twice.The idea is that as the subspace expands, Ritz approximations for other eigenvalues should appear in the extraction.
We now give some more details about the costs of Algorithm 2 with selection, compared to the original version without it.Consider the criterion for the polynomial eigenvalue problem of degree m.We look at (3) and (11) to understand the costs.The selection criterion does not require any extra matrix-vector products (MVs), except for the m + 1 MVs y * A j per detected left eigenvector y, and any MVs necessary to compute this left eigenvector.The costs of computing y depends on the problem.For problems with certain structure, for instance symmetry, the left vector may come for free.For nearly symmetric problems, x may be a good approximation to y.For general problems, several steps with an iterative solver may be necessary to solve for y.A preconditioner that we already have for the eigenvalue problem will generally be of great help.As an alternative for not-too-large problems, an exact solve with P (θ) is often an efficient method to find a very good approximation to y.
Next, we study the costs of the divided differences.Per Ritz vector that we choose to test, we need a product of a "primitive Ritz vector" c with the n × k matrix V k (cost 2nk), some vector additions (cost 2nm) and one inner product for each of the d already detected eigenpairs (cost 2nd).This is summarized in the following table.

Computation
Approx.costs Note y (Variable) Only when eigentriplet found Only when eigentriplet found Every step, per pair to test (7) The dominant costs of lines 2 and 3 in this table will usually be ≈ 2nk per Ritz pair that does not satisfy the criterion (7), so that we have to take the next candidate pair.We have to compute and store the detected left eigenvectors y i , but in contrast to locking, no orthogonalization costs with respect to converged vectors is required as they do not form part of the search space.Note that the extra costs of the selection procedure may be considered low compared to the complexity O(nk 2 ) necessary for orthonormalizing the basis.Also, note that we get valuable extra information about the left eigenvalues and the condition numbers of the eigenvalues.
Example 1 To avoid convergence to the same eigenpairs, one may wonder if the proposed selection criterion based on divided differences ( 5) is really needed: could it be an alternative to simply compare the angles of a current approximate eigenvector v with the already detected eigenvectors x 1 , . . ., x d instead, and make sure that v is sufficiently different?This simple example for the standard eigenvalue problem shows that this is generally not a stable approach; the one based on the divided differences is preferable.
Consider the 2 × 2 matrix for small δ and ε, and suppose that the eigenpair (0, e 1 ) has already been found, where e 1 is the first standard basis vector.The (nonnnormalized) corresponding left eigenvector is y = [δ, −ε] T .Since δ ≈ 0, λ = 0 is close to being a multiple eigenvalue, and therefore the second eigenvector is numerically ill defined: all vectors v with unit norm and associated Rayleigh quotient θ = v * Av have a small residual r = Av − θv.
Comparing the angle of v with e 1 only rules out v that are close to e 1 .With the divided difference requirement y * v ≈ 0 the vector v is forced to be close to the true (non-normalized) second eigenvector [ε, δ] T .

Selection for polynomial eigenproblems in homogeneous coordinates
The context of this section is restricted to polynomial eigenvalue problems (3) rather than the general nonlinear eigenvalue problem.The PEP (3) may have infinite eigenvalues when the leading matrix Am is singular.We therefore study these problems in homogeneous coordinates, and propose a new notion for divided differences in this setting.We consider the QEP (4) for the ease of the presentation, but the techniques carry over to the general polynomial eigenvalue problem (3).Already at this place, we would like to stress the fact that a homogeneous approach may be very valuable in itself, apart from the use for infinite eigenvalues that may be present.As we will see in this section, homogeneous coordinates lead to a different selection criterion (see Experiment 6.1 for a numerical example), that may be seen as a mediator between the selection criterion for the standard QEP and the reverse QEP.Homogeneous techniques are mathematically elegant and pleasing, and account for the fact that eigenvalue problems may be scaled and transformed in various ways.
In the context of divided differences, we have two pairs (α 1 , β 1 ) and (α 2 , β 2 ): a detected eigenvalue λ = α 1 /β 1 , and a new approximation to an eigenvalue θ = α 2 /β 2 .It is common to normalize (or scale) which can be done without loss of generality.Moreover, we might also assume that the β's are real, nonnegative, and in the interval [0, 1]; however, this turns out to be sometimes undesirable in our context of divided differences, for the following reason.When (α 2 , β 2 ) → (α 1 , β 1 ) as projective coordinates, we want the componentwise convergence α 1 → α 2 and β 1 → β 2 in what is to follow.This is not satisfied for, e.g., the pairs (i, ε) and (1, 0), for ε → 0. The first pair converges to the second pair, the infinite eigenvalue 1/0, but there is no component-wise convergence.
Therefore, instead, we scale the pair (α 1 , β 1 ) such that the coordinate with the maximal absolute value is in [0, 1].For convenience of notation, we assume that this maximum coordinate is β 1 , but this is not a restriction.As a result, we know that The other pair is then scaled accordingly to the first pair: β 2 ∈ [0, 1], so that (say) β 2 > 0.7 when (α 2 , β 2 ) is close enough to (α 1 , β 1 ).This scaling ensures that convergence also implies componentwise convergence.
We note that a vector of the form DQ(α, β) x has also been exploited in the context of homogeneous Jacobi-Davidson [14].
Divided differences make a distinction between eigenvectors belonging to the same, and belonging to different eigenvalues.We will use these differences to find out whether approximate eigenpairs are likely to converge to already detected eigenpairs which need to be avoided, or to new pairs.Similar to the nonhomogeneous divided difference Q[λ, µ], we now would like to derive an expression for a divided difference Q[(α 1 , β 1 ), (α 2 , β 2 )] for the homogeneous problem (12), which may be used in the selection process.This is done in the following definition, which is new to the best of our knowledge.
Note that this expression is both elegant and related to the chordal distance.We recall that the chordal distance of two homogeneous numbers (α 1 , β 1 ) and (α 2 , β 2 ) is (see, e.g., [29, p. 139]) Note that this is the sine of the angle between the two projective numbers interpreted as vectors; cf.[3, p. 75].With the standard scaling |α| 2 + |β| 2 = 1 of projective numbers, this implies that the absolute value of the denominator in the first case of Definition 1 is equal to the chordal distance (14).
The next result is a justification of Definition 1: ] is a continuous function of its two variables.
In view of the restriction that the numbers are on the complex unit circle, the following orthogonality condition holds (cf., e.g., [3, Eq. ( 3)]) For the numerator of the divided differences we have where , and that such a procedure also extends to general polynomial eigenvalue problems (3).
Assuming α 1 = 0, we now have, using (15) and the definitions of dα and dβ, and, similarly, In the case that α 1 = 0, it is easy to check that the first expression equals β −1 which is equal to β 1 in this case, and the second expression equals dβ dα −1 β −1 1 = 0 and therefore is equal to −α 1 , since dβ should vanish in view of (15).Suppose (α 1 , β 1 ) is an eigenvalue with right eigenvector x and (α 2 , β 2 ) is an eigenvalue with left eigenvector y.Proposition 5 shows that the following two desirable properties are satisfied: These two properties will help us in the selection process to avoid convergence towards an already detected eigenvalue: we would like to only select candidate approximate eigenpairs (in homogeneous form) ((θ, η), v) for which the divided differences y * i Q[(α i , β i ), (θ, η)] v are small enough for all previously detected eigenvalues (α i , β i ) with corresponding left eigenvectors y i .For this reason, and in line with (7), during the selection process we require that max i=1,...,d where, for instance, η = 0.1.Elegantly, the denominator in this criterion also appears in the denominator of the condition number [3,Thm. 4.2].
The above selection criteria may be exploited for polynomial eigenvalue problems (3), especially if they are expected to have infinite eigenvalues.
We can show that for a PEP the relation between the standard selection criteria and its homogeneous counterpart depends only on the magnitude of the eigenvalues.As a result we see that, if the magnitudes of the new candidate for the eigenvalue and the computed eigenvalue do not differ much, then ( 7) and ( 16) return very close values and both criteria should give the same decision.
Lemma 1 Let x be the eigenvector for a simple eigenvalue λ of the polynomial eigenvalue problem (3) of degree m.If P (α, β) is the homogeneous variant of (3), i.e., P (α, β) = β m P (α/β) for β = 0, then Proof From the partial derivatives DαP (α, β) = β m−1 P (α/β), In the next result, we assume for convenience of presentation that the eigenvalues and approximation are real, as this greatly simplifies the expressions.
Proposition 6 Let y 1 be the left eigenvector of a simple real eigenvalue λ 1 of the polynomial eigenvalue problem (3) of degree m and let P (α, β) be the homogeneous variant of (3).Let x 2 be the right eigenvector for a simple real eigenvalue λ 2 = λ 1 and let (θ, v) = (λ 2 + εφ, x 2 + εw) be a candidate for the next eigenpair, where θ ∈ R. Then where C is given by (9), , and Proof It follows from an expansion that up to O(ε 2 )-terms For the numerator of ( 17) a multivariate Taylor series expansion gives, omitting the O(ε 2 )-terms, where we have applied Lemma 1.If we also use Lemma 1 for the denominator of ( 17) and combine the results, we obtain (17).
Finally, we give a result that elegantly connects the homogeneous divided differences with divided difference of the standard QEP and of the reverse QEP defined by λ 2 Q(λ −1 ) = A + λB + λ 2 C. Let λ and θ be given in homogeneous coordinates by (α 1 , β 1 ) and (α 2 , β 2 ), respectively.By Proposition 3, the divided difference expression for the QEP is (λ + θ)A + B, which in homogeneous coordinates corresponds to A divided difference expression for the reverse QEP is B + (λ −1 + θ −1 )C, or in homogeneous coordinates.The numerator of the homogeneous divided differences as defined in Definition 1 is

It may be checked that
Therefore, the homogeneous approach may be viewed as a mediator between the divided differences of the QEP and reversed QEP.(We note that in homogeneous Jacobi-Davidson [14], the homogeneous vector used in the subspace expansion has also been shown to be a mediator between those arising in the standard and reverse QEP.)In fact, it may be shown that a similar expression also holds for the PEP (3).This result illustrates a mathematically pleasant property of homogeneous coordinates.

Multiparameter eigenvalue problems
As discussed in the introduction, linear two-parameter eigenvalue problems have been the origin of our interest in selection criteria [11,15].We briefly review various previous results, whereby we also improve on our previously proposed criteria.We will keep the discussion as concise as possible, referring to the given references for more information.
Consider the linear two-parameter eigenvalue problem where the task is to find one or more eigenvalues (λ, µ) together with their eigenvectors of the form We first briefly follow [11].
, which we assume to be nonsingular.A left eigenvector y 1 ⊗ y 2 and right eigenvector x 1 ⊗ x 2 corresponding to different simple eigenvalues (λ 1 , µ 1 ) and (λ 2 , µ 2 ), respectively, are ∆ 0 -orthogonal: For a selection criterion, we would like an approximate eigenvector v 1 ⊗ v 2 to be sufficiently ∆ 0 -orthogonal to already detected left eigenvectors y In our previous criterion, as was proposed in [11], we required potential v 1 ⊗ v 2 to satisfy max i=1,...,d (y (y While criterion (18) has turned out to perform satisfactorily in the numerical tests in [11,15], it may be unnecessarily strict: if one eigenvalue has been detected with right and left eigenvector x 1 ⊗ x 2 and y 1 ⊗ y 2 for which the right-hand side of ( 18) is small, the selection procedure may reject many or all candidate Ritz pairs.Therefore, instead of ( 18), we propose the new modified criterion (cf.( 7)) max i=1,...,d (y with, e.g., η = 0.1.This criterion has been successfully used very recently in [12]. In [15] the special but important right-definite case has been treated, where all matrices A i , B i , and C i are Hermitian, and ∆ 0 is positive definite.In this situation, the right and left eigenvectors coincide, and therefore eigenvectors x 1 ⊗ x 2 and x 1 ⊗ x 2 corresponding to different eigenvalues are ∆ 0 -orthogonal: (x 1 ⊗x 2 ) * ∆ 0 ( x 1 ⊗ x 2 ) = 0. We note that [15] has been the first paper where a selection criterion to compute several eigenvalues has been proposed and used, in the context of linear two-parameter eigenproblems.
Besides being simple and easy to implement, selection criteria can be elegantly extended to other types of (multiparameter) eigenvalue problems.The ∆ 0 -orthogonality can be nicely extended to nonlinear two-parameter eigenvalue problems as follows.
For the polynomial and general nonlinear two-parameter eigenvalue problem stands for the operator determinant A ⊗ D − B ⊗ C; see also [25].In these papers, it has been shown that this divided difference has the desired property that the quantity (y is nonzero when (θ, η) converges to (λ i , µ i ), while it is 0 when the pair converges to another eigenvalue.
We now illustrate the adaptivity and flexibility of the selection criterion by the following generalization for the differentiable nonlinear three-parameter eigenvalue problem T 3 (λ, µ, ν) x 3 = 0, While the special case of a linear case of this problem has been treated recently in [12,Lem. 4.2], we now define a divided difference for the nonlinear case (20).
Proposition 7 The quantity ) is a simple eigenvalue with right and left eigenvector and y 1 ⊗ y 2 ⊗ y 3 , respectively; it equals 0 when Proof A rather straightforward generalization of [23,Prop. 3.2] shows that = 0, since the sum of the columns is the zero vector.Finally, when some, but not all, of the coordinates of (λ 1 , µ 1 , ν 1 ) are equal to (λ 2 , µ 2 , ν 2 ), the determinant vanishes as well.Indeed, it can be checked that: when (λ 1 , µ 1 , ν 1 ) and (λ 2 , µ 2 , ν 2 ) agree in one of three coordinates then the columns where the coordinates do not agree differ by a factor −1; -when (λ 1 , µ 1 , ν 1 ) and (λ 2 , µ 2 , ν 2 ) agree in two of three coordinates then the column where the coordinates do not agree is zero.
This result implies that selection criteria in the line of ( 19) can be exploited.
For multiparameter eigenvalue problems, locking becomes less and less attractive as the number of parameters increases.For instance, when we are prepared to solve projected eigenvalue problems of dimension approximately 100 at the subspace extraction step, the search spaces are limited to dimension 10 for two parameters, and even to dimension 5 for three parameters.As locking keeps the converged vectors in the search space, this technique is generally not an option for MEPs.
Therefore, selection effectively creates more space in the subspaces to contain new information.However, even when using selection criteria to compute several eigenvalues, already detected vectors may sometimes turn up in the search space in practice.Therefore, we will be limited by the size of the search space at some point, and we cannot expect to compute arbitrarily many eigenvalues.

Comparison with other approaches
A good comparison of various approaches has already been carried out in [4].Here we briefly discuss differences of selection criteria compared to other methods for computing several eigenvalues of one-parameter eigenvalue problems.
Besides our selection criteria, there are several alternatives for the computation of several eigenvalues for nonlinear eigenvalue problems and linear and nonlinear multiparameter eigenvalue problems.Nonequivalence deflation [5,6,18] has an elegant mathematical foundation, but changes the original problem, and might suffer from instabilities.Block methods may be used to compute several eigenvalues simultaneously [20], but also has some drawbacks as indicated in [4].Locking, which keeps the eigenvectors in the search space [21] [4,Ch. 6], leaves less space for new vectors; to find new eigenvectors, the search spaces have to grow.Especially for multiparameter eigenvalue problems, where the dimension of the projected problems grows as n p , with p the number of parameters, locking is not a realistic alternative.
We will now discuss differences with the method by Effenberger [4], which we consider state-of-the-art and of particular importance, in more detail.This method, as our approach, also computes the eigenvalues successively while preventing convergence to the same eigenpairs.However, this method and the one proposed here are still of very different nature.First, the method in [4] is far from trivial to implement.During the computations, it modifies the original problem by adding rows and columns so that the problem size steadily increases.It is unable to deal with infinite eigenvalues, as it does not use homogeneous coordinates.Moreover, it is an open question if the approach can be generalized to multiparameter eigenvalue problems.As a big advantage, Effenberger's method has been designed with the aim of also computing multiple and clustered eigenvalues in a stable way.Our proposed approach, on the other hand, is (much) simpler, both conceptually and with respect to implementation (just a few lines of codes on top of an existing code).Our method is designed to handle infinite eigenvalues by homogeneous coordinates, and the problem remains unmodified during the iterations.Also, the techniques are elegantly generalizable to various types of eigenproblems.On the other hand, as stated before, the method is not suitable to compute multiple eigenvalues.
We note that standard deflation methods (see Section 1) for the generalized eigenvalue problems can be used for polynomial one-parameter and multiparameter eigenvalue problems when one is prepared to linearize the problem into a (much) larger problem.For instance, in [22], a Krylov-Schur type method has been proposed for the linear two-parameter eigenvalue problem, which works on the operators A main disadvantage of this approach is that it works on vectors of length n 2 , instead of n for a direct approach.The action with the ∆ −1 0 ∆ i operators can be done in O(n 3 ) effort instead of the expected O(n 6 ) by solving a Sylvester equation.Therefore, when n is small enough, this method may still be worthwhile for two-parameter problems.For three-parameter problems, the situation looks far less favorable [12].

Numerical examples
We present some numerical examples obtained with Matlab.Several successful experiments with several types of multiparameter eigenvalue problems have been carried out and described in [11][12][13]15,25].Therefore, we concentrate ourselves mostly on the new use for polynomial eigenvalue problems.
Experiment 1 We consider the QEP utrecht1331 with target τ = −70 − 2000i as in [16], and an exact LU preconditioner based on this target.This is a quite challenging interior eigenvalue problem, due to the difficult spectrum, the interior location of the target, and the different scales of the real and imaginary parts; see Fig. 1(a).Approximate eigenpairs (θ, v) are computed to relative tolerance 10 −6 , meaning At first, we take selection threshold η = 0.1 in (11).We use the Jacobi-Davidson method with harmonic extraction [16] and 10 steps of bicgstab to solve the correction equations.The left eigenvectors are solved by an exact solve with Q(θ) * when θ has sufficiently converged.For the value extraction we use the onedimensional Galerkin gal1 approach from [17].With minimum and maximum subspace sizes of 20 and 40, we find 12 eigentriplets in 200 iterations; the convergence history is displayed in Fig. 1(b).The eigenvalues are detected after 10,12,14,16,78,96,105,118,133,147,165, and 178 iterations.Elegantly, when we sort the eigenvalues with respect to distance to the target, these are eigenvalues number 1 through 12, in this order!The longer "hiccup" after several eigenpairs (here the 5th) may occur in many problems, and is likely due to the fact that new information needs to be inserted in the search space.We note that it seems important that the search spaces are allowed to be sufficiently large; otherwise, at some point, the convergence may stop altogether.For instance, using minimal subspace size 15 and maximal subspace size 25, only 4 eigenvalues are detected in 200 iterations, with indices 5, 7, 6, and 10.Favorably, the process seems to be not very sensitive with respect to the precise threshold value of η: the choices of η = 0.01, 0.2, and 0.5 result in 9, 13, and 13 found eigenpairs, respectively.Although the problem does not have infinite eigenvalues, we may also use the homogeneous divided differences of Section 3. Note that this method is different from the standard divided difference.In this case, we also find 12 eigenpairs in 200 iterations, after 10, 12, 15, 70, 81, 91, 107, 124, 143, 154, 170, and 190 iterations, respectively.
Experiment 2 We consider a popular challenge: the problem gyroscopic, a model of a gyroscopic dynamical system, of size n = 10000; cf.[1, p. 654].Here, A is diagonal with elements uniformly from [0, 1] with additionally a 11 = 0, B is tridiagonal with −1s on the subdiagonal and 1s on the superdiagonal, and C is diagonal with elements uniformly from (−1, 0).Therefore, A is symmetric positive semidefinite, B is skew-symmetric, and C is symmetric negative definite, which is typical for this type of system.The matrix A is singular and the QEP has infinite eigenvalues.Therefore, it seems appealing to exploit the homogeneous technique of Section 3. We take target τ = 80i, and an exact LU preconditioner based on this target.An eigenpair is considered converged if the residual is below 10 −4 .All other parameters are as in Experiment 6.1.We find 10 eigenpairs in 800 iterations, after 108, 109, 111, 115, 118, 143, 162, 176, 625, and 777 iterations, respectively.Here, we see again the same pattern of first spending several iterations to obtain a good subspace, then the quick detection of a number of eigenvalues, followed by a new period of enriching the subspace before new eigenpairs are found.
Experiment 3 For the next experiment, we take the largest cubic polynomial eigenvalue problem (λ 3 A 3 + λ 2 A 2 + λA 1 + A 0 ) x = 0 of the nlevp toolbox [2]: the problem plasma drift, with coefficient matrices of size 512; see Figure 2. We note that this spectrum is quite challenging, with close eigenvalue and eigenvalues of high multiplicity.Our target is τ = 0, and as in the previous experiment we use an exact LU preconditioner based on this target, so LU = A 0 .For the value extraction we use the two-dimensional minimum residual mr2 approach from [17].The other settings are the same as in Experiment 1.With η = 0.1, the Jacobi-Davidson method finds 19 eigenvalues in 200 outer iterations; cf.Fig. 1(c).With respect to distance to the target, these are approximations to eigenvalues with index 1 through 12, 511, 512, 514, 14, 16, 515, and 510, respectively.This "alternating" behavior is quite typical for iterative eigensolvers; cf. also [11,15].The high indices can be explained by the fact that there are several eigenvalues of high multiplicity close to the origin.This illustrates that the selection method may work fine for problems with multiple eigenvalues, as long as the computed eigenvalues are simple.Other choices for η result in 14 (η = 0.01), 10 (η = 0.2), and 11 (η = 0.5) eigenvalues.
where we seek (λ, µ, η) such that there exists a nonzero solution y(x).This problem can be decomposed into a 3-parameter eigenvalue problem that consists of three 2-point boundary value problems of the form for i = 1, 2, 3.A smooth function y(x) that satisfies ( 21) can be constructed from the functions y 1 (x 1 ), y 2 (x 2 ), y 3 (x 3 ).The 3-parameter eigenvalue problem ( 22) has the Klein oscillation property, which means that for each triple of nonnegative integers (m 1 , m 2 , m 3 ) there exist a triple of values (λ, µ, η) such that (21) has a solution y(x) that has m 1 zeros on interval (0, 1), m 2 zeros on (1, 2), and m 3 zeros on (2, 3).We discretize ( 22) using the Chebyshev collocation on 200 points (cf.[12]), which leads to an algebraic 3-parameter eigenvalue problem of the form The solutions with indices (j 1 , j 2 , j 3 ) such that j 1 +j 2 +j 3 is small correspond to eigenvalues (λ, µ, η) close to (0, 0, 0).To find eigenvalues close to the origin, we apply the Jacobi-Davidson method, for details see [12].We restrict the subspace dimensions between 5 and 10 and solve the corresponding correction equations approximately by 10 steps of GMRES, where we use the exact LU preconditioner based on the target, i.e., A j = L j U j for j = 1, 2, 3.The Jacobi-Davidson method returns 20 eigenvalues after performing 40 subspace updates.The first nine eigenvalues converged are provided in Table 1 together with their indices, while the corresponding solutions y(x) of ( 21) are illustrated in Figure 3.Note that the indices in Table 1 confirm that the eigenvalues converged are indeed the ones closest to the origin.We have presented several selection criteria for computing several eigenvalues for nonlinear one-parameter, and linear and nonlinear multiparameter eigenvalue problems.Selection means that an approximate eigen-Fig.3 First 9 solutions of ( 21) corresponding to the eigenvalues listed in Table 1.
pair is picked from candidate pairs that satisfy a certain suitable criterion.The goal of this process is to steer the process away from already previously found pairs.These criteria are easy to understand and implement, and also elegantly extend to various types of eigenproblems.We have also developed a divided difference and selection criterion in homogeneous coordinates.This not only has the potential to handle infinite eigenvalues, but also is a valuable alternative approach in itself.The methods work directly on the original problem; no linearizations (as for instance discussed in [10]) are necessary.They require the computation of the left eigenvector, which implies some extra costs for nonsymmetric problems.However, these additional costs are often relatively small compared to the total costs.For certain problems with structure, such as symmetric problems, the left eigenvectors come for free.Also, left eigenvectors provide valuable information on the condition number and reliability of the computed eigenvalues.
A main advantage of the selection techniques is that the search spaces effectively may contain more useful vectors for the computations of new eigenvectors.Instead of locking, which keeps the converged vectors in the search space, the search spaces can now be more fully used for new information.Moreover, and also important for practical use, the selection criteria are relatively easy to understand and implement compared with several existing approaches.
For the quadratic and polynomial (one-parameter) eigenvalue problem, the presented methods are new, and a valuable alternative to other methods such as locking or block methods (cf.[4]); a more detailed comparison can be found in Section 5.For linear and nonlinear multiparameter eigenvalue problems, we would like to stress the fact that the presented selection techniques seem to be the only realistic option.While for multiparameter eigenvalue problems we already proposed selection criteria in the past, in this paper we propose updated and less strict criteria of the type (19) instead of (18).
The approach can also be applied to general nonlinear eigenproblems F (λ)x = 0, as long as we can evaluate the derivative F (λ) and the divided difference F [λ, µ].
We note that for challenging problems, it sometimes is not easy to find more than about 10 eigenpairs with the selection criterion.Reasons for this may be that a larger part of the search space is occupied by already detected eigenvectors, or that the preconditioner is of lower quality for the new eigenvalues.In this case, it may be a good idea to start a new process with a modified target and preconditioner.

Experiment 4
We consider the 4-point boundary value problem y

Table 1
(21)first 9 eigenvalues of the 4-point boundary value problem(21)retrieved by the Jacobi-Davidson method with the origin as the target point.