 # FEAST Eigensolver for Nonlinear Eigenvalue Problems

The linear FEAST algorithm is a method for solving linear eigenvalue problems. It uses complex contour integration to calculate the eigenvectors whose eigenvalues that are located inside some user-defined region in the complex plane. This makes it possible to parallelize the process of solving eigenvalue problems by simply dividing the complex plane into a collection of disjoint regions and calculating the eigenpairs in each region independently of the eigenpairs in the other regions. In this paper we present a generalization of the linear FEAST algorithm that can be used to solve nonlinear eigenvalue problems. Like its linear progenitor, the nonlinear FEAST algorithm can be used to solve nonlinear eigenvalue problems for the eigenpairs whose eigenvalues lie in a user-defined region in the complex plane, thereby allowing for the calculation of large numbers of eigenpairs in parallel. We describe the nonlinear FEAST algorithm, and use several physically-motivated examples to demonstrate its properties.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The nonlinear eigenvalue problem (NLEVP) consists of finding vectors

and scalars that satisfy

 T(λ)x=0, (1)

where is some matrix-valued function of , which we call the eigenvector residual function. The linear generalized eigenvalue problem (GEP) is a special case of (1), with

 T(λ)=λB−A,   A,B∈Cn×n. (2)

Another well-known example of a NLEVP is a quadratic eigenvalue problem

 T(λ)=λ2A2+λA1+A0,   A2,A1,A0∈Cn×n. (3)

Quadratic eigenvalue problems can arise, for example, in models of physical systems undergoing damped oscillations, e.g., HolGL05 ; LiC15 . In general can take a variety of different forms, depending on the physical model. Higher-degree polynomials are possible, as are non-polynomial functions of , e.g.  MehW02 ; Kau06 ; Sol06 ; Lia07 .

The usefulness of any physical model that produces the problem (1) is limited by the ease with which that problem can be solved. General nonlinear eigenvalue problems carry with them some unique challenges, e.g. their eigenvectors do not form a basis for , and the particular form of can have an impact on which solution methods are most effective. However, one of the core challenges for nonlinear eigenvalue problems is the same as for linear eigenvalue problems: solution methods that are effective for small values of do not necessarily scale well to very large-dimensional problems, where one is usually interested in finding a small number of specific, physically-important eigenvalue-eigenvector pairs.

Most methods for solving such large-dimensional problems require generating an initial guess that is close enough to the desired eigenpairs to ensure convergence, as well as a method for calculating successive eigenpairs whose eigenvalues are close to each other, but without converging repeatedly to the same eigenpair Vos14 . These challenges are exacerbated in the situation that large numbers of eigenpairs are desired. The difficulty of calculating large numbers of eigenpairs is compounded further for NLEVP algorithms that work by approximating the desired eigenvectors in some subspace of dimension ; in this case, calculating large numbers of eigenpairs requires a large value of , leading to bad scaling for the solution of the reduced-dimension eigenvalue problems that are solved inside that subspace. For linear eigenvalue problems the algorithmic complexity of solving the reduced dimension problem of dimension is , and solving nonlinear eigenvalue problems always requires an amount of work that is at least equal to solving a single linear eigenvalue problem; keeping the dimension of the reduced-size problem small is thus of paramount importance for ensuring good scaling to larger problem sizes.

A class of algorithms that address these challenges is contour integral methods, e.g., AsaSTIK09 ; AsaSTIK10 ; Bey12 . Contour integral methods exploit the Cauchy integral formula in order to calculate only the eigenpairs whose eigenvalues lie in a specific, user-defined region of the complex plane. In doing so they eliminate the need to calculate initial guesses and separate closely-positioned eigenvalues for the original, large problem of dimension . Moreover, contour integral approaches allow for the efficient, parallel calculations of large numbers of eigenpairs by dividing the region of interest in the complex plane into disjoint subregions, solving for the eigenvalues in each subregion independently of the eigenvalues elsewhere.

The linear FEAST algorithm is an example of a contour integral method for linear eigenvalue problems. It has the distinction of being able to solve large, generalized linear eigenvalue problems robustly with excellent parallel scalability. In the following sections, we describe a generalization of the linear FEAST algorithm that can be applied to solving nonlinear eigenvalue problems. We describe the resulting nonlinear FEAST (NLFEAST) algorithm and solve some physically-motivated model problems to illustrate its properties.

## 2 Nonlinear Eigenvalue Problems

Given a non-empty open set and a matrix-valued function , we consider a nonlinear eigenvalue problem (NLEVP): Find a scalar and nonzero vectors and such that

 T(λ)x=0,y∗T(λ)=0. (4)

In this paper, we discuss only the nonlinear dependence of on the eigenparameter , leaving the nonlinear dependence on eigenvectors and multi-parameters outside of our focus. Analogously to the linear case, is called a eigenvalue and the corresponding right and left eigenvectors of , respectively. We call the set of all eigenvalues of the spectrum of and denote i by , i.e.,

 σ(T):={λ∈D:det(T(λ))=0}.

We refer to as an eigentriple of , and either or as an eigenpair. Although for the choice and , (4) reduces to the standard eigenvalue problem , and the generalized eigenvalue problem , respectively, the characteristics of a general nonlinear eigenvalue problem differ significantly from those of its linear counterparts, e.g. (even if regular, i.e., det) may have infinitely many eigenvalues, the eigenvectors belonging to distinct eigenvalues do not have to be linearly independent, and the algebraic multiplicity of an isolated eigenvalue, although finite, is not bounded by the problem size  TisG17 . All these particularities make general nonlinear eigenvalue problems much harder to solve. For a more detailed description we refer the readers to the survey papers by Tisseur and Meerbergen TisM01 for quadratic eigenvalue problems, Mackey, Mackey and Tisseur MacMT15 for general polynomial eigenvalue problems, and by Voss Vos14 , and Tisseur and Güttel TisG17 for general nonlinear eigenvalue problems. A variety of applications in science and engineering, e.g., dynamic analysis of structures, vibrations of fluid-solid structures, the electronic behavior of quantum dots are discussed by Mehrmann and Voss MehV04 .

For the sake of simplictly, our main focus in this paper is on a class of square polynomial eigenvalue problems (PEPs): Find a scalar and nonzero vectors such that

 P(λ)x=0,y∗P(λ)=0, (5)

where is an matrix polynomial

 P(λ)=k∑i=0λiAi,Ai∈Cn×n. (6)

If is regular, i.e., det, (5) has finite eigenvalues (roots of det) and infinite eigenvalues. Using linearization, it is easy to see that the matrix polynomial of degree has eigenvalues (finite or infinite) and up to right and left associated eigenvectors. As mentioned before, even in the case of distinct eigenvalues, the associated eigenvectors are not necessarily linearly independent.

For example, let us consider the quadratic eigenvalue problem (QEP): Find a scalar and nonzero vectors such that

 Q(λ)x=(λ2A2+λA1+A0)x=0, (7) y∗Q(λ)=y∗(λ2A2+λA1+A0)=0, (8)

with complex matrices , and . Then the eigenvalues of are equal to those of the companion linearization (GohLR09, , Ch. 1), e.g.,

 L(λ)=λ[A2OOI]+[A1A0−IO](the first companion form)

or

 L(λ)=λ[A2OO−A0]+[A1A0A0O].

The right eigenvectors of are , where are the right eigenvectors of . Similarly, the left eigenvectors of are , where are the left eigenvectors of .

## 3 The FEAST Algorithm

The original FEAST algorithm, introduced for solving the generalized eigenvalue problems, i.e.  Pol09 , is a subspace iteration method that uses the Rayleigh-Ritz projection and an approximate spectral projector as a filter. FEAST calculates the eigenpairs whose eigenvalues lie in a specific, user-defined region in the complex plane. It can therefore be used to calculate a large number of eigenpairs in parallel by dividing the complex plane into a collection of disjoint regions and solving for the eigenpairs in each region independently of the eigenpairs in other regions.

FEAST selects the eigenvectors and eigenvalues of interest by first projecting any random initial subspace onto the subspace spanned by the eigenvectors of interest, and then using the Rayleigh-Ritz procedure in this subspace to extract eigenvalue/eigenvector approximations. Projection onto the subspace of interest is accomplished by using complex contour integration, i.e.

 Q=ρ(A,B)X(0)=12πi∮C(zB−A)−1Bdz X(0), (9)

where is the new, projected subspace, and the integration is performed in a closed contour around the region of the complex plane where we would like to find eigenvalues. If the integration is performed exactly then the filter becomes a spectral projector with , where and are the exact left and right eigenvector subspaces associated with the eigenvalues of interest (i.e. the eigenvalues that are located within the closed curve ). In practice, however, it is impossible to perform the integration in (9) exactly, so instead it is approximated by using a quadrature rule. The subspace that is generated by the practical FEAST algorithm is then

 Q=nc∑j=1ωj(zjB−A)−1BX(0)=nc∑jωjQj, (10)

where and are integration nodes and weights, respectively, and the summation is performed by solving linear systems

. Because each of these linear systems is independent of the others, they can all be solved simultaneously in parallel. FEAST iteratively refines the estimates for the eigenvectors of interest by repeatedly applying the summation (

10), orthogonalizing the resulting subspace and using it in the Rayleigh-Ritz procedure.

In addition to parallelizing the solution of eigenvalue problems by dividing up the eigenspectrum, FEAST has the additional benefit of being able to systematically improve its rate of convergence. The rate of convergence of the FEAST algorithm is related to the accuracy of the approximation for the integral (9) TanP14 . If that integral is evaluated exactly, then the exact solution to the eigenvalue problem is found in a single iteration. Approximating (9) by using (10) generally means that the solution to the eigenvalue problem will be found in some number of iterations larger than one; reducing the number of required iterations is as simple as increasing the number of terms in the quadrature rule summation in (10), which in turn can be accomplished by solving a larger number of linear systems in parallel. The speed with which an eigenvalue problem can be solved by using FEAST is thus limited only by the amount of parallel processing power that is available.

The FEAST numerical library package was first released under the free BSD license in September 2009 (v1.0) to address the Hermitian linear eigenvalue problems. The current FEAST version (v.3) released in June 2015 allows to solve the non-Hermitian linear eigenvalue problems, see KesPT16 .

In the following section, we discuss an extension of the FEAST algorithm to tackle nonlinear eigenvalue problems. Our goal is to solve NLEVPs by taking advantage of the attractive features of the FEAST framework: we would like to solve for eigenpairs whose eigenvalues are in a user-defined region of the complex plane, by iteratively refining a fixed-dimension subspace using contour integration in such a way that the rate of convergence can be enhanced by using parallel computing to improve the accuracy of the contour integration.

## 4 FEAST for Nonlinear Eigenvalue Problems

### 4.1 Previous and Related Work

Other researchers have previously investigated the use of contour integration (using Cauchy’s integral in particular) in the complex plane for solving nonlinear eigenvalue problems.

In a series of papers Asakura at al. introduced the nonlinear variants of the Sakurai-Sugiura (SS) method with block Hankel matrices (SS-H method) for the polynomial eigenvalue problems AsaSTIK09 and eigenvalue problems of analytic matrix functions  AsaSTIK10 . Although they are cost efficient and highly scalable, the accuracy of obtained solutions is relatively low. Almost at the same time Beyn Bey12

introduced a highly accurate algorithm that uses the zeroth and first-order moments matrices to reduce an NLEVP with

eigenvalues inside to a linear eigenvalue problem of dimension . The main idea of Beyn’s integral method is to probe a Jordan decomposition of the matrix following the Keldysh’s theorem (MennickenMoller, , Theorem 1.6.5), which is conceptually very simple but known to be highly sensitive to perturbations. Moreover, since the value of parameter (the number of linearly independent eigenvectors) is not known in advance, the practical realization of Beyn’s algorithm requires various adaptations which makes its overall computational cost relatively high.

Recently, Yokota and Sakurai YokS13 addressed the problem of low accuracy in the nonlinear SS-H method. Their method and the nonlinear variant of the Sakurai-Sugiura method with block Hankel matrices (SS-H) AsaSTIK09 use the same contour integrals of the SS-H method, however, they differ in the way the approximate eigenpairs are extracted from the underlying subspaces. The method of Yokota and Sakurai is a projection-based method which uses the Rayleigh-Ritz procedure to obtain the approximate eigenpairs from a subspace. It does not require any fixed point iterations and gives better accuracy than the methods of Asakura et al. AsaSTIK09 and Beyn Bey12 .

Both the various SS methods and the Beyn method use the moments of the Cauchy integral of , i.e.,

 12πi∮CzkT−1(z)dz,  k≥0. (11)

Beyn’s method uses the and moments, and the SS-type methods use as many moments as are necessary to reach convergence. The original FEAST algorithm uses only the moment, which, in the case when is linear, turns out to be a spectral projector associated with the eigenvalues of interest. Unfortunately, subspace iteration using only the moment in (11) does not work in the case of the nonlinear . When is nonlinear, taking an initial set of approximate eigenvectors and refining it using the quadrature approximation of the moment of (11

), as is done in the linear FEAST algorithm, does not bring the resulting subspace closer to the desired eigenspace. Although Beyn’s method resolves this issue by continuously refining the accuracy of the contour integration for a single multiplication, it requires performing a new matrix factorization for each newly added quadrature node. The SS-type methods, on the other hand, refine their solutions by increasing the dimension of their search subspace by calculating additional moments of (

11) using the same quadrature rule each time; the dimension of the search subspace is increased iteratively until the desired solution is sufficiently accurate.

We propose, instead, to modify the contour integral (11) such that the subspace iteration framework of the original FEAST algorithm can be used, thereby making it possible to solve nonlinear eigenvalue problems using a fixed-dimension subspace that is refined by solving linear systems at a constant number of fixed shifts in the complex plane.

### 4.2 The Nonlinear FEAST Algorithm

To develop a nonlinear FEAST algorithm that can use contour integration to iteratively refine a subspace of a fixed dimension by using fixed-location shifts in the complex plane, we propose to study the following (modified) form of the contour integral

 12πi∮CI−T−1(z)T(λ)z−λdz, (12)

where is the current Ritz estimation of an eigenvalue that is inside the contour . For the implementation of nonlinear FEAST, we use a block version of (12) in order to generate a refined a subspace from an initial set of approximate eigenvectors of a nonlinear eigenvalue problem, i.e.,

 Q=12πi∮C(X(0)−T−1(z)T(X(0),Λ(0)))(zI−Λ(0))−1dz, (13)

where is the block eigenvector residual function for the initial estimate for the eigenvectors , and the diagonal matrix of the corresponding Ritz values is . For example, for the polynomial eigenvalue problem in equation (6), the block form of the residual function is

 T(X,Λ)=k∑i=0AiXΛi. (14)

The nonlinear FEAST algorithm follows the same essential steps as the linear FEAST algorithm: a subspace is formed by refining an initial set of estimated eigenvectors by using the contour integration in (13), then a new set of approximate eigenvectors is found in the resulting subspace by using the Rayleigh-Ritz procedure to solve a projected eigenvalue problem (which is still nonlinear, but with a substantially smaller dimension). As in the linear case we use a numerical integration scheme to evaluate the integral in (13), and we solve a linear system whenever a multiplication by is required.

The contour integral in (13) is mathematically equivalent to (9) for the linear, generalized eigenvalue problem. However, the two are different when depends nonlinearly on . The relationship between the contour integrals for linear and nonlinear FEAST is much like the relationship between the shift-and-invert iteration for the linear eigenvalue problem and the residual inverse iteration for the nonlinear eigenvalue problem. The residual inverse iteration Neu85 is a modification of the shift-and-invert iteration which allows it to be applied to to nonlinear problems while using a constant shift . The original FEAST algorithm can be interpreted as a generalization of the shift-and-invert iteration that uses contour integration in order to efficiently use multiple shifts simultaneously. We note that the nonlinear FEAST algorithm, in turn, can be interpreted as a generalization of residual inverse iteration that uses contour integration in order to efficiently use multiple shifts simultaneously.

In this paper, we will present the algorithmic framework and various example applications of the NLFEAST algorithm, leaving the mathematical details to a forthcoming paper. The full NLFEAST algorithm is described in Algorithm 1.

The Rayleigh-Ritz step in the nonlinear FEAST works the same as in its linear counterpart: the residual function is projected onto the subspace , and the resulting nonlinear eigenvalue problem of reduced-dimension is solved using any suitable method, here we use a simple linearization. The resulting Ritz values and associated Ritz vectors are used as the new estimates of the desired eigenpairs.

Careful attention must be paid while selecting the desired eigenpairs from the solutions of the nonlinear eigenvalue problem of reduced-dimension. In many cases, such as those considered in this paper, it is possible to find all the solutions of the reduced-dimension eigenvalue problem; for the nonlinear FEAST algorithm, however, we want the Rayleigh-Ritz procedure to return only a number of eigenpairs that is equal to the dimension of the FEAST subspace . The desired eigenpairs are thus selected by choosing those whose eigenvalues are closest to being inside the FEAST search contour

. This heuristic appears to ensure convergence to only the eigenpairs whose eigenvalues are inside the region of interest in the complex plane.

The subspace that is used for the Rayleigh-Ritz procedure, i.e., the one that is generated by performing the contour integration, should be orthonormalized, for example by using the decomposition. Orthonormalization improves the numerical stability of the FEAST iterations by preventing the norms of the desired eigenvectors from diverging and by preventing the Rayleigh-Ritz procedure from producing spurious eigenpairs, i.e., the Ritz pairs that do not correspond to any eigenpairs of the full-size eigenvalue problem.

In the following section we describe the results of using Algorithm 1 to solve several example nonlinear eigenvalue problems.

## 5 Numerical Examples

In this section, we illustrate the behavior of the nonlinear FEAST algorithm on several polynomial eigenvalue problems of practical relevance.

###### Example 1.

A nonoverdamped mass-spring system

Let us consider the Hermitian (self-adjoint) quadratic eigenvalue problem (HQEP), i.e., , associated with the mass-spring system discussed in (TisM01, , §3.9), i.e.,

 T(λ)x=(λ2A2+λA1+A0)x=0,

with

 A2=In,A1=τ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣3−100…0−13−10…00−13−1…0⋮⋮⋮⋱⋮⋮0…0−13−10…00−13⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,

and

 A0=κ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣3−100…0−13−10…00−13−1…0⋮⋮⋮⋱⋮⋮0…0−13−10…00−13⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

The eigenvalues of a Hermitian are either real or they come in complex conjugate pairs . Since matrices and are real, the right and left eigenvectors coincide.

Following (LiC15, , §4.2), we first study the nonoverdamped system. We set and choose , . All the eigenvalues of are plotted in Figure 1. One possible problem of interest is to calculate the eigenvectors whose eigenvalues have no imaginary part; the exponentially increasing or decaying solutions to the original equations of motion for the mass-spring system are linear combinations of these eigenvectors. In this case, we seek the real eigenvalues lying within the interval , see Figure 2.

We use the nonlinear FEAST approach presented in Algorithm 1 to compute approximations of those real eigenvalues. The contour is chosen as an ellipse centered at the midpoint of the interval with radius on the real and on the imaginary axis using integration nodes, and the convergence criteria for the residuum is . With the subspace of size , the nonlinear FEAST algorithm finds all eigenvalues within the interval in three FEAST iterations, see Figure 2. The eigenvalue approximations obtained using the linearization approach and the nonlinear FEAST algorithm are listed in Table 1. For completeness, we also list all the final residuals, i.e., .

Table 2 shows the number of NLFEAST iterations that are required for convergence when the search interval is enlarged. When the number of quadrature nodes is kept constant at , as it is here, then enlarging the search interval has the same effects as placing the quadrature nodes farther away from the eigenvalues of interest, making the contour integration itself less accurate. Modest changes to the size of the interval result in only a slightly larger amount of work being required for convergence; larger changes in the interval size tend to require proportionally more NLFEAST iterations. In practice, decreased performance due to the size of the search interval relative to the distribution of the desired eigenvalues can be mitigated simply by using a larger number of quadrature nodes which, in turn, requires only that one use more parallel processing power in order to solve more linear systems simultaneously. Figure 1: All eigenvalues of the nonoverdamped mass-spring system for n=1000, τ=0.6202 and κ=0.4807 obtained via linearization. Figure 2: Zoom towards all the real eigenvalues of the nonoverdamped mass-spring system for τ=0.6202, κ=0.4807. The NLFEAST search contour is shown as a green dashed line, with the locations of the quadrature shifts indicated by green squares. Exact solutions to the nonlinear eigenvalue problem, obtained by linearization, are shown as blue circles, and the approximations obtained by using NLFEAST are indicated by red triangles. Using a subspace of dimension m0=22, NLEAST correctly identifies the 20 eigenvalues inside of the search interval, in addition to the two eigenvalues that are closest to the contour without being inside of it.
###### Example 2.

An overdamped mass-spring system

We now consider the same mass-spring system as in Example 1, but with new parameter values and . This choice of parameters results in a hyperbolic QEP with real and non-positive eigenvalues (TisM01, , §3.9). For all eigenvalues of obtained via linearization are plotted in Figure 3. Half of the eigenvalues are clustered very close to zero, and the other half are distributed throughout the interval . The spectral gap between the lowest and the highest eigenvalues is a characteristic of this class of problems.

In the case of an overdamped system there are no oscillations; the dynamics of the system consist entirely of decreasing exponentials. The fastest dynamics of the system are determined by the eigenvectors whose eigenvalues come before the eigenvalue gap, i.e., the ones with the largest magnitudes. As an example, we calculate the real eigenvalues lying within the interval . The contour is chosen as a circle centered at the midpoint of the interval with radius , using quadrature nodes, and the convergence criteria for the residual is . With the subspace of size , the nonlinear FEAST algorithm finds all eigenvalues within the interval in ten () FEAST iterations.

Figure 4 illustrates the the eigenvalues of the problem, the FEAST contour and quadrature points, and the resulting eigenvalue estimates that are calculated by NLFEAST. Figure 3: All eigenvalues of the overdamped mass-spring system with τ=10, κ=5 and n=50 obtained via linearization. The eigenvalues we seek with NLFEAST are highlighted in red. A large number of eigenvalues is clustered relatively close to zero; these are separated from the rest by a spectral gap. Figure 4: The nonlinear FEAST approximations of m=19 eigenvalues inside interval (−30,−11) of the overdamped mass-spring system for n=50, τ=10 and κ=5.
###### Example 3.

Scattering resonances in 1D

Another application in which nonlinear eigenvalue problems arise is open boundary quantum transmission problems. With open boundary quantum transmission problems, one seeks to calculate the quasi-bound states of a quantum potential where particles can either enter or leave the system, preferably without having to explicitly model the external sources or sinks to which that quantum potential is connected. This can be done by solving a nonlinear eigenvalue problem ShaPL95 , in contrast to many other problems in quantum mechanics that can be approached by solving a linear eigenvalue problem with the system Hamiltonian.

As an example of open boundary quantum transmission problems, let us consider the problem of scattering resonances in D with the following compactly supported finite square model potential

 V(r)={−V0,r∈[−L,L]0,otherwise, (19)

with and width  CanN17 .

We are interested to determine the Siegert states Sie39 and the associated resonances such that for all

 L∫−Lu′v′+Vu¯u′dx =ık(u(L)¯v(L)+u(−L)¯v(−L)) (20) + k2L∫−Lu¯vdx.

In the case of being bounded with compact support in , the set of discrete solutions of (20) can be approximated using the finite element space associated with the grid , . This discretization approach results in the following quadratic eigenvalue problem

 T(kh)uh=(k2hAh+ıkhBh−Ch)uh=0, (21)

where

 Ah=h6⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣2100…01410…00141…0⋮⋮⋮⋱⋮⋮0…01410…0012⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ ∈ Rn+2×n+2,
 Bh=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣10…0000…00⋮⋮⋮⋮⋮00…0000…01⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ ∈ Rn+2×n+2,

and

 Ch=1h⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1−100…0−12−10…00−12−1…0⋮⋮⋮⋱⋮⋮0…0−12−10…00−11⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦+V0Ah ∈ Rn+2×n+2.

The associated linear eigenvalue problem has a form

 [−ıBhChIO][khuhuh]=kh[AhOOI][khuhuh]. (22)

A finite element discretization of (20) over the grid of points results in a quadratic eigenvalue problem of size with scattering resonances plotted in Figure 5. We are interested in all twenty-two complex scattering resonances lying inside the circle centered at with radius , see Figure 5. The nonlinear FEAST computes accurate approximations of those twenty-two () scattering resonances in four iterations using integration nodes and the subspace of size , see Figure 6. Figure 5: All scattering resonances for potential V0=10 obtained via linearization (22). Figure 6: The nonlinear FEAST approximations of m0=30 eigenvalues for the scattering resonance problem for potential V0=10. NLFEAST captures all of the 22 eigenvalues inside of the integration contour, plus the eight (8) eigenvalues that are closest to the integration contour while still being outside of it.
###### Example 4.

Quartic Eigenvalue Problem. As an example of a polynomial eigenvalue problem with degree larger than two, let us consider the following quartic problem

 P(λ)x=(λ4A4+λ3A3+λ2A2+λA1+A0)x=0.

Eigenvalue problems of this form can come from, for example, discretizations of the Orr-Sommerfeld equation BriM84 ; TisH01 . The Orr-Sommerfeld equation arises from a linearization of the incompressible Navier-Stokes equation in which the perturbations of the pressure and velocity are assumed to be periodic in time.

To illustrate the behavior of the nonlinear FEAST algorithm, we use a simple example of a quartic eigenvalue problem provided by the NLEVP collection nlevp_collection : the so-called butterfly problem (distribution of eigenvalues in the complex plane resembles a shape of a butterfly). The butterfly problem is a structured quartic matrix pencil with eigenvalues, the construction of which is described in MehW02 . We use the NLFEAST algorithm to calculate the eigenvalues that are located inside of some arbitrarily chosen region in the complex plane. This problem is illustrated in Figure 7. Figure 7: An illustration of the locations of the eigenvalues of the butterfly problem in the complex plane, as well as the search region C in which we calculate eigenvalues using NLFEAST.

We calculate 13 eigenvalues inside of the indicated region by using a subspace of dimension , using several different numbers of quadrature nodes . The largest (at each iteration) eigenvector residual associated with the eigenvectors whose eigenvalues are inside the search region is plotted in Figure 8. Using quadrature nodes, NLFEAST does not converge at all. We can achieve steady convergence by using quadrature nodes, and the rate of convergence increases with increasing values of . For , convergence to the desired tolerance of occurs in only five NLFEAST iterations. Because the linear system for each individual quadrature node is independent of the linear systems associated with all the other quadrature nodes, the rate of convergence of NLFEAST can be systematically improved by using additional parallel processing power to solve a larger number of linear systems simultaneously in parallel. Figure 8: Convergence of the NLFEAST algorithm for the eigenvectors whose eigenvalues are inside the search region when solving the butterfly problem. Convergence trajectories are shown for different numbers of integration quadratute points nc. NLFEAST is not able to converge effectively when using only 8 quadrature points, and converges rapidly when using 128 points. When using enough parallel processing power, a single iteration takes the same amount of time regardless of the number of quadrature points that are used.

Figure 9 shows the NLFEAST-estimated eigenvalues for the experiments from Figure 8 that use and quadrature points in the numerical integration. The case is not able to converge because the integration is not sufficiently accurate to achieve convergence of the two eigenvalues that are farthest-separated from the main cluster of eigenvalues; using allows the accurate convergence of NLFEAST for all of eigenvalues inside the search region . Figure 9: Illustrations of the NLFEAST-estimated eigenvalues for the butterfly problem. The results when using nc=8 quadrature nodes (top), and the results when using nc=32 quadrature nodes (bottom). Increasing the number of quadrature points improves the quality of the integration, allowing NLFEAST to better-identify eigenvalues that are farther away from the main cluster.

## 6 Summary and conclusions

In this paper we have described an extension of the linear FEAST eigenvalue algorithm for solving nonlinear eigenvalue problems. The resulting nonlinear FEAST algorithm (NLFEAST) uses the same series of operations as the linear FEAST algorithm, but with a modified contour integral that allows for using a fixed collection of shifts and a fixed subspace dimension in solving nonlinear eigenvalue problems. Where the linear FEAST algorithm can be interpretted as a generalization of shift-and-invert iterations that can use multiple shifts, the nonlinear FEAST algorithm can be interpretted as a generalization of residual inverse iterations that can use multiple shifts.

Like the linear version, the NLFEAST algorithm can be used to calculate eigenvectors corresponding to the eigenvalues located in a specific, user-defined region in the complex plane, allowing for the parallel calculation of large numbers of eigenpairs. Moreover, analogously to the linear FEAST, the convergence rate of the NLFEAST can be systematically improved by solving additional linear systems in parallel to refine the numerical contour integration.

In this paper, for the sake of simplicity, we have treated only polynomial eigenvalue problems. This makes the implementation of NLFEAST relatively straight-forward, in the sense that the reduced, nonlinear eigenvalue problem in the NLFEAST algorithm can be solved easily via linearization. We expect, however, that the NLFEAST algorithm as described in this paper will prove to be a general method of solution for nonlinear eigenvalue problems of any form, rather than just polynomial eigenvalue problems. In the case of more general (non-polynomial) nonlinear eigenvalue problems, the reduced nonlinear eigenvalue problem will be itself of an arbitrary form, and therefore will require using a Newton- or a projection-type solution method instead of a simple linearization technique.

## Acknowledgements

The authors would like to thank Eric Cances and Boris Nectoux for providing the scattering problem example. The work of Agnieszka Miedlar was partially supported by the Simons Foundation Collaboration Grant for Mathematicians Eigenvalue Computations in Modern Applications. This work was also supported by the National Science Foundation, under grant #CCF-1510010.

## References

• (1) U. B. Holz, G. H. Golub, K. H. Law, A subspace approximation method for the quadratic eigenvalue problem, SIAM J. Matrix Anal. Appl. 26 (2) (2005) 498–521.
• (2) H. Li, Y. Cai, Solving the real eigenvalues of Hermitian quadratic eigenvalue problems via bisection, Electron. J. Linear Algebra 30 (2015) 721–743.
• (3) V. Mehrmann, D. Watkins, Polynomial eigenvalue problems with Hamiltonian structure, Electron. Trans. Numer. Anal. 13 (2002) 106–118.
• (4) L. L. Kaufman, Eigenvalue problems in fiber optic design, SIAM J. Matrix Anal. Appl. 28 (1) (2006) 105–117.
• (5) S. I. Solov’ëv, Preconditioned iterative methods for a class of nonlinear eigenvalue problems, Linear Algebra Appl. 415 (1) (2006) 210–229.
• (6) B.-S. Liao, Subspace projection methods for model order reduction and nonlinear eigenvalue computation, Ph.D. thesis, Dept. of Mathematics, University of California at Davis (2007).
• (7) H. Voss, Nonlinear Eigenvalue Problems - Chapter 60, in: L. Hogben (Ed.), Handbook of Linear Algebra, Second Edition, Discrete Mathematics and its Applications, Chapman & Hall/CRC, Boca Raton, FL, 2013, pp. 1063–1086.
• (8) J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, K. Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM Lett. 1 (2009) 52–55.
• (9) J. A. T., Sakurai, H. Tadano, T. Ikegami, K. Kimura, A numerical method for polynomial eigenvalue problems using contour integral, Jpn. J. Ind. Appl. Math. 27 (1) (2010) 73–90.
• (10) W. J. Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra Appl. 436 (10) (2012) 3839–3863.
• (11) S. Güttel, F. Tisseur, The nonlinear eigenvalue problem, Acta Numer. 26 (2017) 1–94.
• (12) F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev. 43 (2) (2001) 235–286.
• (13) D. S. Mackey, N. Mackey, F. Tisseur, Polynomial eigenvalue problems: theory, computation, and structure, in: Numerical algebra, matrix theory, differential-algebraic equations and control theory, Springer, Cham, 2015, pp. 319–348.
• (14) V. Mehrmann, H. Voss, Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods, GAMM Mitt. Ges. Angew. Math. Mech. 27 (2) (2004) 121–152 (2005).
• (15) I. Gohberg, P. Lancaster, L. Rodman, Matrix polynomials, Vol. 58 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2009, reprint of the 1982 original [MR0662418].
• (16) E. Polizzi, Density-matrix-based algorithm for solving eigenvalue problems, Phys. Rev. B (79) (2009) 115112.
• (17) P. T. P. Tang, E. Polizzi, FEAST as a subspace iteration eigensolver accelerated by approximate spectral projection, SIAM J. Matrix Anal. Appl. 35 (2) (2014) 354–390.
• (18) J. Kestyn, E. Polizzi, P. T. P. Tang, FEAST eigensolver for non-Hermitian problems, Tech. Rep. 5 (2016).
• (19) R. Mennicken, M. Möller, Non-self-adjoint boundary eigenvalue problems, Vol. 192 of North-Holland Mathematics Studies, North-Holland Publishing Co., Amsterdam, 2003.
• (20) S. Yokota, T. Sakurai, A projection method for nonlinear eigenvalue problems using contour integrals, JSIAM Lett. 5 (2013) 41–44.
• (21) A. Neumaier, Residual inverse iteration for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 22 (5) (1985) 914–923.
• (22) Z.-a. Shao, W. Porod, C. S. Lent, D. J. Kirkner, An eigenvalue method for open-boundary quantum transmission problems, Journal of applied physics 78 (4) (1995) 2177–2186.
• (23) E. Cancès, B. Nectoux, Private Communication (2017).
• (24) A. J. F. Siegert, On the derivation of the dispersion formula for nuclear reactions, Phys. Rev. 56 (1939) 750–752.
• (25) T. J. Bridges, P. Morris, Differential eigenvalue problems in which the parameter appears nonlinearly, J. Comput. Phys. 55 (3) (1984) 437–460.
• (26) F. Tisseur, N. J. Higham, Structured pseudospectra for polynomial eigenvalue problems, with applications, SIAM J. Matrix Anal. Appl. 23 (1) (2001) 187–208.
• (27) T. Betcke, N. J. Higham, V. Mehrmann, C. Schröder, F. Tisseur, NLEVP: A collection of nonlinear eigenvalue problems, http://www.mims.manchester.ac.uk/research/numerical-analysis/nlevp.html.