A Subspace Iteration Algorithm for Fredholm Valued Functions

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.


Introduction
In this paper we study analytic operator eigenvalue problems defined on an open connected subset Ω ⊆ C in a separable Hilbert space H. Throughout this paper we assume that  : Ω → L(H) is an analytic function with values in the space L(H) of bounded linear operators.A scalar  ∈ Ω is called an eigenvalue of  if () is not injective.Hence, the eigenvalue problem is to find  ∈ Ω and  ∈ H \ {0} such that  ()  = 0. ( Such problems are, for example, used to study the dispersion and damping properties of waves [1][2][3].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  −1 .The state-of-the-art results in the contour integration based methods for solving nonlinear matrix eigenvalue problems are presented in [4][5][6] 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 L(H) and the resolvent set is defined as the complement () = Ω \ ().For  ∈ L(H) we define the operator norm by the expression ‖‖ = √spr( * ), where spr(⋅) denotes the spectral radius.We call an operator  ∈ L(H) a Fredholm operator if the dimensions of its null space Ker() and of the orthogonal complement of its range CoKer() = Ran() ⊥ are finite.By Φ(H) we denote the set of all Fredholm operators on H and the number ind() = dim Ker() − dim CoKer() is called the index of  ∈ Φ (H).In what follows we will assume that () ∈ Φ(H) 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   → () −1 is finitely meromorphic.This in turn implies that the spectrum () is countable and the geometric multiplicity of , that is, dim(Ker ()), 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 16 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 2 Mathematical Problems in Engineering 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.

Notation and Basic Analytic Results
In this section we present the machinery of quasi-matrices from [10][11][12].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 C  to an (in general) infinite dimensional Hilbert space H.Then, the product  *  denotes the Gram matrix: which depends on the inner product (⋅, ⋅) of H. Let Then, The multiplication of two block operator matrices also follows the rules of matrix algebra.
Here we have exemplarily taken the trivial partition of unity  1 =  and  0 = 0 and we assume that both  1 ̸ = 0 and  2 ̸ = 0.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 C  and  will denote the identity operator on H.
We call  the maximal canonical angle between the spaces Ran() and Ran( 1 ) = Ran( 1 ).Since ‖‖ < 1 (7) implies that  *  is a positive definite matrix, therefore  must be invertible.By direct computation using spectral calculus, as in [15], we establish that cos  = ‖‖ , When we are considering several pairs of quasi-matrices  1 and  we will write ( 1 , ) to denote the maximal canonical angle between subspaces Ran( 1 ) and Ran().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].Mathematical Problems in Engineering 3

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 t 0 be a densely defined closed symmetric form in H 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 Dom(t 0 ) ⊂ H to denote its domain of definition.To simplify the notation we will write V := Dom(t 0 ) in the rest of the paper.Additionally, let c  ,  = 1, . . ., , be a sequence of sesquilinear (not necessarily symmetric) forms which are relatively compact [16] with regard to t 0 .Moreover, assume that   (⋅),  = 1, . . ., , 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  ∈ V \ {0} 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  0 , Dom( 1/2 0 ) = V be a self-adjoint positive definite operator which represents the form t 0 in the sense of Kato's second representation theorem; see [16,Theorem VI.2.23].Let   be defined by The operators   ,  = 1, . . ., , are obviously compact and for  ∈ Ω \ () the value of the generalized resolvent is the operator: Here () is the unbounded sectorial operator with domain Dom(()) 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.
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 [17][18][19].It states that for the closed contour Γ ⊂ () the number satisfies (, Γ) ∈ N and it equals the total multiplicity of the eigenvalues enclosed by Γ.We also have the following consequence of [9, Theorem 1.3.1].(12) with the operator representation  : Ω → H 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   ,  = 1, . . ., , such that    = dim Ker((  )).Then there are quasi-matrices   ,   :

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  : Ω → Φ(H) be an analytic function and let  be a bounded operator such that dim Ran() = 1.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  ∈ [0, 1], then  +  has a simple eigenvalue λ inside Γ and where  essentially depends on max ∈Γ ‖() −1 | Ran() ‖ and max ∈Γ ‖() − * | Ran( * ) ‖ and the length of the the integration curve Γ.
Proof.Recall that (() + )  =   () and define the function: By the generalized argument principle-see [17][18][19]-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 (  (), ) =  *   () = 1 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 V and  so that  = V * and using Sherman-Morrison formula (see [21][22][23] and the references therein), we write and so it follows that And 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  ∈ [0, 1].Let us denote this eigenvalue of S =  +  by λ.Define the operators Proposition 2 implies  (1) =  * , Ã(1) = λ x ỹ * , where () = 0, () *  = 0, S( λ) x = 0, and S() * ỹ = 0 and so there exists a vector  such that  *  (0)  ̸ = 0 and  * Ã(0)  ̸ = 0. We now compute and so the second assertion follows by the following computation: The claim on the constant  follows from the observation that Ran() = span{} and Ran( * ) = span{V}.
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 t(⋅) denote the family of sesquilinear forms (11) with operator representation (⋅) : Ω → H given by (15).Assume that a contour Γ ⊂ () encloses only a simple eigenvalue  and  ∈ C but no other points of ().Let  ∈ V, with ‖‖ = 1.Then there is  ∈ () such that The constant  does not depend on  and  but on contour Γ and the restriction of ‖ −1 (⋅)‖ 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 e , : V × V → C by the formula It is obviously relatively compact with respect to t 0 and so we can define the Fredholm valued operator function Ω ∋  → T(), where T() is the operator defined by the form Further, and so  ∈ ( T).We construct the operator  , as the operator defined by Now recall that ( 22) and (35) imply By construction the operator function S : Ω → L(H) satisfies the assumption of Lemma and so from (36) it follows that Now recall the definition of  from Lemma 3 to conclude the proof.

A Sketch for a Practical Algorithm for Error Estimation.
Proposition 4 will be the basis for practical error estimation.Let V ⊃ Q ℎ , ℎ > 0 be a family of finite dimensional spaces such that the orthogonal projections  ℎ onto Q ℎ converge strongly to  as ℎ → ∞.Let further Obviously it holds, recalling the definition from (39), that However, in the case in which Q ℎ , ℎ > 0 satisfies the standard saturation property, for example, ‖(− ℎ 2 )‖ ≤ ‖(− ℎ 1 )‖ for fixed , 0 <  < 1 and some , we can also prove The constant  essentially depends on  which in turn depends on ℎ 2 − ℎ 1 but not on the magnitude of ℎ 1 ; for more details see [24,25].

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   ∈ Γ,  = 1, . . ., , 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 () ∩ Γ = 0, () = min ∈() dist(, Γ), 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  ()  →  () ,  = 0, 1 in the norm resolvent sense.

Extracting Information on Eigenvalues Enclosed by a
Contour.Based on Proposition 2 we see that operators  (0) and  (0) are finite rank operators such that Ran( (0) ) = Ran( (0) ) 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  (1) and  (1) .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  () ,  = 0, 1, as in (44) and let  : C  → H,  *  = , be a quasi-matrix such that  *  (0)  and  *  (0)  are invertible.Then the eigenvalues of the matrix pairs ( *  (1) ,  *  (0) ) and ( *  (1) ,  *  (0) ) are precisely the eigenvalues   ,  = 1, . . ., , where we count according to multiplicity.Furthermore, if   ,   ∈ C  satisfy  *  (1)   =    *  (0)   and  *  (1)   =    *  (0)   , then   and   are eigenvectors of  associated with   ,  = 1, . . ., .
Proof.We will prove the statement for the matrix pair ( *  (1) ,  *  (0) ) and note that the proof for the matrix pair ( *  (1) ,  *  (0) ) is equivalent.Based on Proposition 2 we see that there are quasi-matrices  : C  → H and  : C  → H and a neighborhood U 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  () ,  = 0,1.In particular, this means that  (0)   and  (0)  , 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 V the vectors  (0)  V and  (0)  V can be computed using formula (45) without ever forming a representation for the underlying operator.
operator  of type (44) which has a component of ,  ∈ N dominant eigenvalues and whose action on a vector can be represented by a formula (45).
Proof.Let  1 =  * and represent  : C  → H with respect to the partition of unity with Ran( 1 ) = Ran( 1 ).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  (0) as Without reducing the level of generality we may assume that  (0) 1 is invertible, since otherwise Ran( (0) ) would contain eigenvectors of .This corresponds to the trivial situation and will not be further discussed.
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 Ran( () ) and Ran() = Ran( 1 ).In Theorem 6 we have stated the convergence result for the angle ( () , ).Here we have ignored the choice of a basis of Ran( () ) 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 Σ 1 and Σ 2 so that () = Σ 1 ∪ Σ 2 .Let further Σ 1 = { 1 ,  2 , . . .,   }, where we have counted   ∈ () according to the algebraic multiplicity of an eigenvalue.Further let a partition of unity  1 ⊕  2 =  exist such that  = dim(Ran( 1 )) and and let further ( 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  +1 which is returned by Algorithm 1.More precisely we have the following.For simplicity we write  =  (0) .Theorem 10.Let the conditions of Theorem 8 be satisfied and also assume the same notation conventions.For a quasi-matrix  : C  → H, we set  1 :=  and apply Algorithm 1 with fixed  ∈ N. Then Proof.Recall that From standard results on norm convergence of bounded operators (see [16]) we conclude that sin  ( 1 , ( Here (  ) 1 + (  ) 2 =  is the partition of unity in which   has a block triangular form: Assumption on the spectral separation from Theorem 8 implies that  11 is invertible and obviously Norm convergence (64) implies that for  large enough (  ) 11 must also be invertible; for details see [30].
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 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.

Auxiliary Subspace Error Estimates and the Convergence Criterion.
In this paper we will use the inequality ‖t()[, ⋅]‖ Q ℎ 2 ,t 0 ,−1 ≤ ‖t()[, ⋅]‖ t 0 ,−1 to filter the eigenpairs which have been extracted using Proposition 5 from the subspace returned by Algorithm 1.More precisely, eigenpairs for which ‖t()[, ⋅]‖ Q ℎ 2 ,t 0 ,−1 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  1 =  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  1 which are identified by Algorithm 2 as converged and then proceed as in the original versions of Algorithm 1.

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 V ℎ (, ).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 Q ℎ (, ) to denote the space of piecewise quadratic and continuous function on the same subdivision of the interval [, ] which was used to define the space V ℎ (, ).To estimate the error  , from (39) for the function  ∈ V ℎ (, ) and  ∈ C we use the formula; see [24] Example 12.We study the quadratic eigenvalue problem: for  = 0.05 and  = 0.3.As a benchmark we use the eigenvalue computed with a spectral discretization using the chebfun system; see [10].The timings are 2.5 seconds for chebfun and 1.6 seconds for subspace iteration.We used the trapeze formula to define the operator   with  = 25 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 (  ) 11 and (  ) 22 can be achieved with few integration nodes : for example,  = 25.Figure 2(a) depicts the dependence of the effective convergence rate for the inexact subspace iteration (Γ, ) = ‖(  ) 22 ‖‖(  ) −1 11 ‖ 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 The results in Figures 3 and 4 are presented solely for benchmark purposes.They were computed using trapeze formula with 4500 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 [−9/10, 9/10] 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  (+1) =    () ,  0 : C 2 → V, where  = 16 and   and   are Gauss-Legendre integration nodes and weights from [5].We have used the contour which encloses only two wellseparated 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  = 1.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 criterion     t () [, ⋅]    Q ℎ (,),t 0 ,−1 ≤ tol 2 (78) in 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 4500 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  = 16 and 3 steps of inexact subspace iteration acceleration we reach the same accuracy (see Table 1) as a priori predicted for  = 4500 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.

Stability of Convergence Rate of Inexact Subspace Iteration.
From Figure 2 we see that the convergence rate for the subspace iteration essentially depends on ‖ 22 ‖‖ −1 11 ‖ ≈ (Γ, ).We will now see that this convergence rate is relatively robust 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 ‖ +1 −   ( *   +1 )‖ 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 by the Croatian Research Foundation Grant HRZZ-9345 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.

Table 1 :
Convergence history for the rational eigenvalue problem Example 14 and the two leftmost eigenvalues from Figure3.