Abstract

We present an algorithm for approximating an eigensubspace of a spectral component of an analytic Fredholm valued function. Our approach is based on numerical contour integration and the analytic Fredholm theorem. The presented method can be seen as a variant of the FEAST algorithm for infinite dimensional nonlinear eigenvalue problems. Numerical experiments illustrate the performance of the algorithm for polynomial and rational eigenvalue problems.

1. Introduction

In this paper we study analytic operator eigenvalue problems defined on an open connected subset in a separable Hilbert space . Throughout this paper we assume that is an analytic function with values in the space of bounded linear operators. A scalar is called an eigenvalue of if is not injective. Hence, the eigenvalue problem is to find and such thatSuch problems are, for example, used to study the dispersion and damping properties of waves [13]. Given a closed contour , we would like to approximate all the eigenvalues of inside to the sufficient degree of accuracy. In this paper the numerical method is based on contour integrals of the generalized resolvent . The state-of-the-art results in the contour integration based methods for solving nonlinear matrix eigenvalue problems are presented in [46] including the references therein. Results for contour integration based solution methods for Fredholm valued eigenvalue problems can be found in, for example, [7, 8].

The spectrum of the operator function is the set of all such that is not invertible in and the resolvent set is defined as the complement . For we define the operator norm by the expression , where denotes the spectral radius. We call an operator a Fredholm operator if the dimensions of its null space and of the orthogonal complement of its range are finite. By we denote the set of all Fredholm operators on and the number is called the index of . In what follows we will assume that for all . If in addition the resolvent set of such is nonempty, the analytic Fredholm theorem, for example, [9, Theorem 1.3.1], implies that the generalized resolvent is finitely meromorphic. This in turn implies that the spectrum is countable and the geometric multiplicity of , that is, , is finite. Moreover, the associated Jordan chains of generalized eigenvectors have finite length bounded by the algebraic multiplicity; see [9].

The results of this paper are a combination of matrix techniques from [6] with the specialization of the results from [7] to Hilbert spaces. In particular, we leverage the technique of block operator matrix representation of Fredholm valued operator functions and prove that the convergence rate for the inexact subspace iteration algorithm depends primarily on the spectral properties of the operator function. Also, we make the case for the problem dependent determination of the number of integration nodes depending on the clustering of eigenvalues towards the contour of integration (cf. [6]), where the use of nodes of Gauss-Legendre integration formula is recommended. To assess the quality of a computed eigenpair, rank one perturbations of the operator function are studied. In particular, we construct a perturbation based on the residual functional and estimate the approximation errors by estimating the norm of the residual by an auxiliary subspace technique. Our algorithm consists of the inexact subspace iteration for the zeroth moment of the resolvent to construct the approximate eigenspace for the eigenvalues contained inside a contour . We then use the moment method of Beyn et al. [4, 7] to extract information on eigenvalues from the computed approximate eigenspace. As the convergence criterion we use a hierarchical residual estimate. In the case in which the convergence criterion is not satisfied the procedure is repeated. This structure of the algorithm will be replicated directly in the structure of the paper and is presented as Algorithm 2.

The paper is organized as follows. In Section 2 we establish a criterion based on the residual norm for assessing approximations of simple eigenvalues. In Section 3 we present inexact subspace iteration algorithm based on contour integration and prove that its convergence rate essentially depends on the properties of the operator function. It is shown that the influence of the integration formula diminishes exponentially with the number of integration nodes. Finally, in Section 4 we present numerical experiments.

2. Notation and Basic Analytic Results

In this section we present the machinery of quasi-matrices from [1012]. In particular we present basic results on the angles between finite dimensional subspaces of a Hilbert space in terms of quasi-matrix notation. Finally, we will prove an error estimation result for simple eigenvalues of a Fredholm valued function.

A quasi-matrix is a bounded linear operator from the finite dimensional space to an (in general) infinite dimensional Hilbert space . Then, the product denotes the Gram matrix: which depends on the inner product of . Let , , , and be orthogonal projections such that and , where is the identity operator on . Furthermore, let , , , and identify and the two spaces and by isomorphism. Let and take Then, , , where is defined by restricting onto appropriate spaces. Hence, the bounded operator has a block operator matrix representation in terms of and is computed following the rules of matrix algebra: The multiplication of two block operator matrices also follows the rules of matrix algebra. Represent, for example, by , where and and . Then, the operator has the block representation: Here we have exemplarily taken the trivial partition of unity and and we assume that both and . This illustrates the flexibility we have in choosing the block operator matrix representations for bounded operators. For more details on this construction see a recent book by . Tretter [13]. To make the paper more readable we will use to denote the identity operator on the finite dimensional space and will denote the identity operator on .

Let be given such that . Let be a unitary operator such that and and . Note that in this setting we have . For the quasi-matrix such that we can write where and . With this notation we compute and so and . Furthermore, note that is the orthogonal projection onto and so from [14, Theorem 2] it follows thatThe last identity has been established by spectral calculus using the fact that ; see [15]. Let ; then we define the unique number such that We call the maximal canonical angle between the spaces and . Since (7) implies that is a positive definite matrix, therefore must be invertible. By direct computation using spectral calculus, as in [15], we establish that When we are considering several pairs of quasi-matrices and we will write to denote the maximal canonical angle between subspaces and .

This definition can be extended to subspaces of different dimensions using pairs of orthogonal projections and their singular value decomposition [14, 15]. In this case we call the norm of the difference of the two projections the gap distance between subspaces; see [14, 15].

2.1. Application of Abstract Results to Operators Defined by Sesquilinear Forms

The abstract Fredholm analytic theorem is stated for bounded operators between Hilbert spaces. Since we are interested in finite element computations, our problems will always be stated in the variational form which assumes working with unbounded sectorial operators. See [16] for definitions and the terminology relating to unbounded sectorial operators and forms. Below we will formulate the Fredholm analytic theorem in this setting.

Let be a densely defined closed symmetric form in which is semibounded from below with a strictly positive lower bound; see [16, Section VI.2]. We will call such quadratic forms positive definite forms and use to denote its domain of definition. To simplify the notation we will write in the rest of the paper. Additionally, let , , be a sequence of sesquilinear (not necessarily symmetric) forms which are relatively compact [16] with regard to . Moreover, assume that , , is a sequence of scalar analytic functions in . Define the family of sesquilinear forms: In a variationally posed eigenvalue problem we are seeking a scalar and a vector such that To this variational formulation we construct a representation by an operator-valued function and then apply the analytic Fredholm theorem to establish the structure of the spectrum.

Let , be a self-adjoint positive definite operator which represents the form in the sense of Kato’s second representation theorem; see [16, Theorem VI.2.23]. Let be defined by The operators , , are obviously compact and for the value of the generalized resolvent is the operator:Here is the unbounded sectorial operator with domain such that Obviously the operator-valued function satisfies the requirements of the analytic Fredholm theorem and . Let us note that we will use to define a derivative of an operator-valued analytic function. We can now define the notion of a semisimple eigenvalue.

Definition 1. Let be as in (15) and let be an eigenvalue. The eigenvalue is semisimple if for each there is a such that If , then is called a simple eigenvalue.

Note that in the quasi-matrix notation we will freely write To this end we identify the vectors with a mapping .

With these conventions we state, informally, the generalized argument principle proved by Gohberg and Sigal in [1719]. It states that for the closed contour the number satisfies and it equals the total multiplicity of the eigenvalues enclosed by . We also have the following consequence of [9, Theorem 1.3.1].

Proposition 2. Let us assume that we have a variational eigenvalue problem (12) with the operator representation given by (15). Then the spectrum consists of a countable collection of eigenvalues with finite multiplicity. Further, let the component of inside a contour consist only of semisimple eigenvalues , , such that . Then there are quasi-matrices , , such that and , , and an open once connected neighborhood containing , and and an operator-valued function which is analytic on and taking values in such that and , .

Proof. For write and define the Fredholm valued function: Recall that and that maps one to one on . Now apply [9, Theorem 1.3.1] on .

2.2. Estimating Eigenvalues inside a Contour

To count the semisimple eigenvalues inside a contour we will use the approach of [20]. We will limit our consideration on the case of semisimple eigenvalues and rank one perturbations. First we will present results for Fredholm valued operator functions and then formulate the result for operators defined by sesquilinear forms.

Lemma 3. Let be an analytic function and let be a bounded operator such that . Assume that is a simple closed contour such that the component of inside consists solely of a simple eigenvalue . If is invertible for all and all , then has a simple eigenvalue inside and where essentially depends on and and the length of the the integration curve .

Proof. Recall that and define the function: By the generalized argument principle—see [1719]—the value of equals the total multiplicity of the eigenvalues inside . In particular, by Proposition 2 we have that there are vectors and such that and where we used the circularity of the trace. By the assumptions of the theorem is invertible for all and is a rank one bounded operator. Therefore, there are vectors and so that and using Sherman-Morrison formula (see [2123] and the references therein), we write and so it follows thatAnd in particular is a smooth function. Since , due to Rouche’s theorem, conclude that takes values only in natural numbers, it must be constant for all . Let us denote this eigenvalue of by . Define the operators Proposition 2 implies , , where , , , and and so there exists a vector such that and . We now compute and so the second assertion follows by the following computation: The claim on the constant follows from the observation that and .

With this result we can now formulate the main result of this section. It will be used to assess the quality of a given approximation, regardless of its origin.

Proposition 4. Let denote the family of sesquilinear forms (11) with operator representation given by (15). Assume that a contour encloses only a simple eigenvalue and but no other points of . Let , with . Then there is such that The constant does not depend on and but on contour and the restriction of and to .

Proof. We will now construct a particular operator which will be used to assess the quality of the pair , . Define the sesquilinear form by the formula It is obviously relatively compact with respect to and so we can define the Fredholm valued operator function , where is the operator defined by the form Further, and so . We construct the operator as the operator defined by Now recall that (22) and (35) imply By construction the operator function satisfies the assumption of Lemma 3. Finally, note that and so from (36) it follows thatNow recall the definition of from Lemma 3 to conclude the proof.

2.3. A Sketch for a Practical Algorithm for Error Estimation

Proposition 4 will be the basis for practical error estimation. Let , be a family of finite dimensional spaces such that the orthogonal projections onto converge strongly to as . Let further for . Assume that and are given. For , define Obviously it holds, recalling the definition from (39), that However, in the case in which , satisfies the standard saturation property, for example, for fixed , and some , we can also prove The constant essentially depends on which in turn depends on but not on the magnitude of ; for more details see [24, 25].

3. Contour Integration Based Subspace Iteration

In the following section we will present the inexact subspace iteration algorithm based on contour integration and prove basic convergence results using quasi-matrix notation. We consider the spectral transform functions from [5, 6, 26, 27] in the context of eigenvalues of operator-valued analytic functions. Let be given and let be a closed curve which encloses either a set of simple eigenvalues or a single semisimple eigenvalue whose multiplicity is .

Let us consider the operators and their approximations: Here and , , are integration weights and integration nodes. Based on Theorem  4.7 in [4] we establish, for and defined by the node trapeze integration formula for the contour integral (44), the following estimates: Here is simple closed contour in such that , , and is the maximum order of poles for the inverse of . The constants in (46) are in general different and depend on the maximum of the integrand on the contour . For more details see [4, 7]. Subsequently we conclude that and , in the norm resolvent sense.

3.1. Extracting Information on Eigenvalues Enclosed by a Contour

Based on Proposition 2 we see that operators and are finite rank operators such that is the space spanned by linearly independent eigenvectors associated with semisimple eigenvalues which are enclosed by . Rather than providing a technical proof of these claims, which will be a subject of subsequent reports, we will present numerical evidence on two judiciously chosen examples. Further, we see that based on [7] we can establish the following technical result for the operators and .

Proposition 5. Let be the operator-valued function from (15) and let be the contour which encloses solely the , counting according to algebraic multiplicity, semisimple eigenvalues of . Define the matrices and , , as in (44) and let , , be a quasi-matrix such that and are invertible. Then the eigenvalues of the matrix pairs and are precisely the eigenvalues , , where we count according to multiplicity. Furthermore, if satisfy and , then and are eigenvectors of associated with , .

Proof. We will prove the statement for the matrix pair and note that the proof for the matrix pair is equivalent. Based on Proposition 2 we see that there are quasi-matrices and and a neighborhood of such that Here is an analytic operator-valued function. Now we compute Based on the assumptions of the theorem we conclude that and have to be invertible and so the conclusion readily follows.

Before we proceed, note that norm resolvent convergence of and implies the convergence of spectra and associated spaces to those of and , . In particular, this means that and , for large enough, will have two well separated components in the spectrum and so we are motivated to use subspace iteration on the operator level to improve the quality of the approximate eigenvalues. Note that for a vector the vectors and can be computed using formula (45) without ever forming a representation for the underlying operator.

3.2. Inexact Subspace Iteration Based on Quadrature Formula

We have the following algorithm for a generic bounded operator of type (44) which has a component of , dominant eigenvalues and whose action on a vector can be represented by a formula (45).

For Algorithm 1 we present standard convergence results in Theorems 6 and 8. First we consider a special case in which is a Hermitian (bounded self-adjoint) operator. We will see that a possible presence of singularities in the operator function will not pollute the convergence rate.

Pick a random quasi-matrix ,
Set
repeat
 Compute
 Compute an orthogonal quasi-matrix , using Gram-Schmidt so that .
 Compute  criterion =
until  criterion < tol
return  

Require  , , and .
 Compute as output from Algorithm 1 with tolerance and as starting value.
 Compute and .
 Form ,
 Compute and so that ,
 Compute eigenvectors , .
 Compute
 Form the index set
 Set
return   and ,

Theorem 6. Let be a bounded Hermitian operator such that , where , , are eigenvalues of finite multiplicity and , and , . Then for the sequence , where , , we have where , , is such that and , .

Proof. Let and represent with respect to the partition of unity by , with . Since is self-adjoint it can be represented by the block operator matrix: Let us represent, with respect to the same partition of unity, the quasi-matrix as Without reducing the level of generality we may assume that is invertible, since otherwise would contain eigenvectors of . This corresponds to the trivial situation and will not be further discussed.
Since it follows from (6) that The assumption of the separation of the spectra of implies that must also be invertible, and so for which implies thatAnd so from the optimality of the canonical angles (see [15]), it follows thatConsequently, it follows thatNow the conclusion of the theorem readily follows.

Remark 7. Note that from (8) and (9) it is clear that does not depend on the quasi-matrices and , but rather on the orthogonal projections onto and . In Theorem 6 we have stated the convergence result for the angle . Here we have ignored the choice of a basis of which is a very important part of Algorithm 1. Note that in practical computations the choice of the basis is crucial to achieve an efficient implementation; see [28].

In the case in which the operator is not Hermitian but has clearly separated invariant subspace associated with the finite collection, counting according to the algebraic multiplicity of dominant eigenvalues we have the theorem below.

Theorem 8. Let be a bounded operator such that its spectrum has two disjoint components and so that . Let further , where we have counted according to the algebraic multiplicity of an eigenvalue. Further let a partition of unity exist such that and and let further and . Then for a quasi-matrix we have Here is the solution of the Sylvester equation and .

Proof. For such that with satisfying the assumptions of the theorem choose as the orthogonal projection onto the maximal invariant subspace associated with and ; then has the block operator representation as claimed in the theorem. We note that where . Operator , such that exists and is unique, for example, [29, Lemma A.1], providing the spectra of and , are disjoint. The rest of the proof follows analogously as in Theorem 6.

Remark 9. In [29, Lemma A.1] there is also an estimate of in terms of the pseudospectra distances.

This result implies that the convergence rate for the subspace iteration essentially depends on the properties of the operator and not on the subspace which is iterated. Further, the dependence on the number of integration nodes diminishes exponentially. To extract eigenvalues and eigenvectors we use the method from [7] and apply Proposition 5 on the quasi-matrix which is returned by Algorithm 1. More precisely we have the following. For simplicity we write .

Theorem 10. Let the conditions of Theorem 8 be satisfied and also assume the same notation conventions. For a quasi-matrix , we set and apply Algorithm 1 with fixed . Then

Proof. Recall that From standard results on norm convergence of bounded operators (see [16]) we conclude that Here is the partition of unity in which has a block triangular form: Assumption on the spectral separation from Theorem 8 implies that is invertible and obviously Norm convergence (64) implies that for large enough must also be invertible; for details see [30]. For such we have Triangle inequality implies thatand so Let us now assume that is large enough so that Then we can apply Theorem 8 to and conclude that Equivalently since and in norm we conclude that We now apply the triangle inequality for the sine of the maximal canonical angle to conclude the proof.

Using Theorem 10 we will define the effective convergence rate of the inexact contour based subspace iteration. Recall that Starting from (63) we define the effective convergence rate for the inexact subspace iteration for as

Remark 11. We will project on carefully constructed finite dimensional spaces, for example, finite element spaces, to experimentally estimate in the following section. We will not further elaborate on this procedure but solely present the results of experiments for illustration purposes.

3.3. Auxiliary Subspace Error Estimates and the Convergence Criterion

In this paper we will use the inequality to filter the eigenpairs which have been extracted using Proposition 5 from the subspace returned by Algorithm 1. More precisely, eigenpairs for which is too large will be discarded. We will indicate the importance of this step of the algorithm in the experiments section.

In the case in which the number of eigenvalues returned by Algorithm 2 is not satisfactory, the procedure can be extended as an innerouter iteration scheme by setting and repeating the procedure. One further possibility to improve performance is to modify Algorithm 1 so that at the beginning of the repeat loop we remove all those directions from which are identified by Algorithm 2 as converged and then proceed as in the original versions of Algorithm 1.

4. Numerical Experiments

For the finite element discretizations, we use the space of piecewise linear and continuous finite elements on a given subdivision of an interval . We denote this space by . For our computations we set the tolerance tol so as to balance the error when solving the source problem by finite element approximation with the error in the integration. For the contour we chose a circle and as integration nodes and weights we use the trapeze formula from [7] and Gauss-Legendre nodes and weights from [5].

To estimate the approximation errors in our experiments we will use an auxiliary subspace error estimation technique. To this end we use to denote the space of piecewise quadratic and continuous function on the same subdivision of the interval which was used to define the space . To estimate the error from (39) for the function and we use the formula; see [24] for further information:

Example 12. We study the quadratic eigenvalue problem: for and . As a benchmark we use the eigenvalue computed with a spectral discretization using the chebfun system; see [10]. The timings are seconds for chebfun and seconds for subspace iteration. We used the trapeze formula to define the operator with and achieved the residuals as small as 1e-12. The results are presented in Figure 1. From Figure 1 we see that the spectrum of (75) is well separated, and so the sufficient separation of the spectra of and can be achieved with few integration nodes : for example, . Figure 2(a) depicts the dependence of the effective convergence rate for the inexact subspace iteration on the number of integration nodes for the trapeze formula. The contour was a circle which enclosed the five eigenvalues marked by crosses.

Remark 13. Note that in [5, 6] it was shown that Gauss-Legendre and trapeze integration require similar number of integration nodes to be applicable in inexact subspace iteration algorithm. A slight experimental advantage was noticed in favor of Gauss-Legendre, but both approaches can be used in competitive algorithms depending on the particular context of the application.

Example 14. We will now consider the following rational eigenvalue problem:whereis the characteristic function of the interval .

This eigenvalue problem has a singularity at and . The eigenvalues below the singularity and those in are discrete (finite multiplicity) and they accumulate at and at from below [31].

The results in Figures 3 and 4 are presented solely for benchmark purposes. They were computed using trapeze formula with nodes and validated using the approach from [8]. Looking at Figures 3 and 4 we see that in general there are two classes of eigenfunctions: those that are confined to the interval and those that are not. Further discussion of the qualitative properties of eigenfunctions is outside the scope of this paper. We will address the modeling aspects elsewhere.

As a first experiment in Table 1 we present the empirical study of the convergence of the residual measures and the relative eigenvalue error computed from the sequence , , where and and are Gauss-Legendre integration nodes and weights from [5]. We have used the contour which encloses only two well-separated eigenvalues, two left most crosses, in Figure 3.

The convergence rate for the well-separated eigenvalues and the standard choice of the integration nodes from [6] seem satisfactory. However, the eigenvalues of (76) cluster towards . In Figure 2(b) we see that many more integration nodes are necessary to achieve a similar estimate of the effective convergence rate compared to the quadratic eigenvalue problem in Figure 2(a). To this end we will compare below the measured convergence rates for the criterion in Algorithm 1 on Examples (76) and (75).

Remark 15. Results presented in Figure 5 show that, unlike what was suggested in [6], the number of integration nodes might have to be adaptively adjusted for some contours . Also we emphasize that the residual criterionin Algorithm 2 is particularly important in the case of the clustering of eigenvalues towards like in Figures 3 and 4. In comparison with benchmark results from Figures 3 and 4 three spurious eigenvalues, which would otherwise have been declared as converged, were discarded because their residuals were larger than a threshold. With this in mind we justify the application of the modified algorithm from [7]. The algorithm in [7] had residual error control but controlled the accuracy solely by increasing the number of integration nodes (e.g., we used nodes for benchmark results in Figures 3 and 4). In our modification we combine subspace iteration with contour integration with problem adapted number of nodes (e.g., with and steps of inexact subspace iteration acceleration we reach the same accuracy (see Table 1) as a priori predicted for based on [8]). Finally, we note that performance issues are outside the scope of this paper. As is indicated in [6] the algorithm offers large potential to leverage parallel processing. This will be the topic of further research.

4.1. Stability of Convergence Rate of Inexact Subspace Iteration

From Figure 2 we see that the convergence rate for the subspace iteration essentially depends on . We will now see that this convergence rate is relatively robust to perturbations as long as the distance from the curve and is not too small. To this end we consider the problem as a perturbation of (75). Its eigenvalues are presented on Figure 6. We see that both the eigenvalues of (75) and (79) are equally well separated from and so the convergence rates, as measured by the decay rate of in Algorithm 1, are similar. However, when comparing the results from Figure 7(c) with those from Figure 7(a) we see that we need more iterations to achieve a comparable reduction in the convergence criterion for Algorithm 1 in the case of eigenvalue clustering towards . On the other hand, Figure 7(b) indicates that the nonlinear perturbation did not change the convergence rate significantly. This justifies the consideration of the adaptivity both in the choice of integration nodes and in representing the operators.

5. Conclusion

We have shown that the subspace iteration nonlinear eigensolver based on spectral transformation of the analytic Fredholm valued function converges at rates that depend primarily on the problem and not on the discretization. Further, we have seen that the distance from the contour and the spectrum does limit the accuracy of the approximation based on numerical integration of the resolvent. Also, residual norm is a reliable estimator of the approximation error even in the presence of poles. In comparison, the convergence in the case of the quadratic eigenvalue problem which had a more pronounced spectral gap than the rational problem was much faster and more robust. This suggests that a subspace iteration algorithm that increases the number of integration nodes based on a posteriori computable criterion is promising for nonlinear eigenvalue problems.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Christian Engström gratefully acknowledges the support of the Swedish Research Council under Grant no. --. The work of Luka Grubišić has in part been supported by the Croatian Research Foundation Grant HRZZ- and the MZOS-NSF grant “Estimates for Finite Element Approximation Error by Auxiliary Subspace Method.” Luka Grubišić is thankful to J. Ovall, Portland State University, for the hospitality during the initial phase of the work on this paper as well as for many helpful discussions.