# A Subspace Iteration Algorithm for Fredholm Valued Functions

Authors: Christian Engström; Luka Grubišić
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 1. 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 S:Ω→L(H) is an analytic function with values in the space L(H) of bounded linear operators. A scalar λ∈Ω is called an eigenvalue of S if S(λ) is not injective. Hence, the eigenvalue problem is to find λ∈Ω and u∈H∖{0} such that(1)Sλu=0.Such problems are, for example, used to study the dispersion and damping properties of waves [1–3]. Given a closed contour Γ, we would like to approximate all the eigenvalues of S inside Γ to the sufficient degree of accuracy. In this paper the numerical method is based on contour integrals of the generalized resolvent S-1. The state-of-the-art results in the contour integration based methods for solving nonlinear matrix eigenvalue problems are presented in [4–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σ(S) of the operator function S is the set of all λ∈Ω such that S(λ) is not invertible in L(H) and the resolvent set is defined as the complement ρ(S)=Ω∖σ(S). For F∈L(H) we define the operator norm by the expression F=spr(F∗F), where spr(·) denotes the spectral radius. We call an operator F∈L(H) a Fredholm operator if the dimensions of its null space Ker⁡(F) and of the orthogonal complement of its range CoKer(F)=Ran(F)⊥ are finite. By Φ(H) we denote the set of all Fredholm operators on H and the number ind(F)=dim⁡Ker(F)-dim⁡CoKer(F) is called the index of F∈Φ(H). In what follows we will assume that S(λ)∈Φ(H) for all λ∈Ω. If in addition the resolvent set of such S is nonempty, the analytic Fredholm theorem, for example, [9, Theorem 1.3.1], implies that the generalized resolvent z↦S(z)-1 is finitely meromorphic. This in turn implies that the spectrum σ(S) is countable and the geometric multiplicity of λ, that is, dim(Ker⁡S(λ)), 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 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 Section2 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 [10–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 operatorV from the finite dimensional space Cr to an (in general) infinite dimensional Hilbert space H. Then, the product V∗V denotes the Gram matrix: (2)V∗Vij=Vei,Vej,i,j=1,2,…,r,which depends on the inner product (·,·) of H. Let P1, P2, Q1, and Q2 be orthogonal projections such that P1⊕P2=I and Q1⊕Q2=I, where I is the identity operator on H. Furthermore, let Hi≔PiH, H^i≔QiH, i=1,2, and identify H and the two spaces H1⊕H2 and H^1⊕H^2 by isomorphism. Let B∈L(H) and take (3)u=u1u2∈H1⊕H2≔H.Then, QiBu=Bi1u1+Bi2u2, i=1,2, where Bij:Hj→H^i is defined by restricting Bij≔QiBPj onto appropriate spaces. Hence, the bounded operator B has a block operator matrix representation B:H1⊕H2→H^1⊕H^2 in terms of Bij and Bu is computed following the rules of matrix algebra: (4)Bu=B11B12B21B22u1u2=B11u1+B12u2B21u1+B22u2.The multiplication of two block operator matrices also follows the rules of matrix algebra. Represent, for example, A:H→H by A=A1A2, where H=H1⊕H2 and A1:H1→H and A2:H2→H~. Then, the operator AB:H→H has the block representation: (5)AB=A1B11+A2B21A1B12+A2B22.Here we have exemplarily taken the trivial partition of unity Q1=I and Q0=0 and we assume that both P1≠0 and P2≠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 C. Tretter [13]. To make the paper more readable we will use Ir to denote the identity operator on the finite dimensional space Cr and I will denote the identity operator on H.LetP1⊕P2=I be given such that dim⁡Ran(P1)=r. Let V be a unitary operator such that V=V1V2 and Ran(P1)=Ran(V1) and Ran(P2)=Ran(V2). Note that in this setting we have P1=V1V1∗. For the quasi-matrix X:Cr→H such that X∗X=Ir we can write (6)X=V1C+V2S,where C=V1∗X and S=V2∗X. With this notation we compute (7)Ir=X∗X=C∗C+S∗S,and so S≤1 and C≤1. Furthermore, note that Q=XX∗ is the orthogonal projection onto Ran(X) and so from [14, Theorem 2] it follows that(8)P1-Q=V1V1∗-XX∗=Ir-CC∗-CS∗-SC∗-SS∗=S.The last identity has been established by spectral calculus using the fact that dimRan(P1)=r; see [15]. Let S<1; then we define the unique number θ∈[0,π/2〉 such that (9)sin⁡θ=S.We call θ the maximal canonical angle between the spaces Ran(X) and Ran(P1)=Ran(V1). Since S<1 (7) implies that C∗C is a positive definite matrix, therefore C must be invertible. By direct computation using spectral calculus, as in [15], we establish that (10)cos⁡θ=C,tan⁡θ=SC-1.When we are considering several pairs of quasi-matrices V1 and X we will write θ(V1,X) to denote the maximal canonical angle between subspaces Ran(V1) and Ran(X).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.Lett0 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(t0)⊂H to denote its domain of definition. To simplify the notation we will write V≔Dom(t0) in the rest of the paper. Additionally, let ci, i=1,…,d, be a sequence of sesquilinear (not necessarily symmetric) forms which are relatively compact [16] with regard to t0. Moreover, assume that fi·, i=1,…,d, is a sequence of scalar analytic functions in Ω. Define the family of sesquilinear forms: (11)tzϕ,ψ≔t0ϕ,ψ+∑i=1dfizciϕ,ψ,ϕ,ψ∈V,z∈Ω.In a variationally posed eigenvalue problem we are seeking a scalar λ∈Ω and a vector ϕ∈V∖{0} such that (12)tλϕ,ψ=0,ψ∈V.To this variational formulation we construct a representation by an operator-valued function T(·) and then apply the analytic Fredholm theorem to establish the structure of the spectrum.LetT0, Dom(T01/2)=V be a self-adjoint positive definite operator which represents the form t0 in the sense of Kato’s second representation theorem; see [16, Theorem VI.2.23]. Let Ki be defined by (13)Kiu,v=ciT0-1/2u,T0-1/2v,u,v∈H.The operators Ki, i=1,…,d, are obviously compact and for z∈Ω∖σ(T) the value of the generalized resolvent is the operator:(14)T-1z=T0-1/2I+∑j=1dfjzKj-1T0-1/2.Here T(z) is the unbounded sectorial operator with domain Dom(T(z)) such that (15)Tzϕ,ψ=t0ϕ,ψ+∑i=1dfizciϕ,ψ,ϕ∈DomTz,ψ∈V.Obviously the operator-valued function (16)Sz=I+∑j=1dfjzKj,z∈Ωsatisfies the requirements of the analytic Fredholm theorem and σ(S)=σ(T). Let us note that we will use (17)T′zϕ,ψ=∑i=1dfi′zciϕ,ψ,ϕ∈DomT′z,ψ∈V,to define a derivative of an operator-valued analytic function. We can now define the notion of a semisimple eigenvalue.Definition 1.
LetT be as in (15) and let μ∈σ(T) be an eigenvalue. The eigenvalue μ is semisimple if for each ψ∈Ker(T(μ))∖{0} there is a ψ∈Ker(T∗(μ)) such that (18)T′μϕ,ψ≠0.If dim⁡Ker⁡(T(μ))=1, then μ is called a simple eigenvalue.Note that in the quasi-matrix notation we will freely write(19)T′μϕ,ψ=ψ∗T′μϕ.To this end we identify the vectors with a mapping ψ:C→H.With these conventions we state, informally, the generalized argument principle proved by Gohberg and Sigal in [17–19]. It states that for the closed contour Γ⊂ρ(T) the number (20)MT,Γ=tr∫ΓTz-1T′zdz=tr∫ΓT′zTz-1dzsatisfies M(T,Γ)∈N 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 T:Ω→H given by (15). Then the spectrum σ(T) consists of a countable collection of eigenvalues with finite multiplicity. Further, let the component of σ(T) inside a contour Γ consist only of semisimple eigenvalues λi, i=1,…,r, such that nλi=dim⁡Ker(T(λi)). Then there are quasi-matrices Xi,Yi:Cnλi→H, i=1,…,r, such that Ran(Xi)⊂Ker(T(λi)) and Ran(Yi)⊂Ker(T∗(λi)), i=1,…,r, and an open once connected neighborhood U containing λi, i=1,…,r and Γ and an operator-valued function H which is analytic on U and taking values in L(H) such that (21)T-1z=∑i=1r1z-λiXiYi∗+Hz,z∈U,and Yi∗(T′(λi)Xi)=Inλi, i=1,…,r.Proof.
Forz∈Ω∖σ(T) write (22)T-1z=T0-1/2I+∑j=1dfjzKj-1T0-1/2and define the Fredholm valued function: (23)Sz=I+∑j=1dfjzKj,z∈Ω.Recall that Dom(T01/2)=V and that T01/2 maps V one to one on H. Now apply [9, Theorem 1.3.1] on S.


### 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.
LetS:Ω→Φ(H) be an analytic function and let E be a bounded operator such that dim⁡Ran(E)=1. Assume that Γ⊂ρ(S) is a simple closed contour such that the component of σ(S) inside Γ consists solely of a simple eigenvalue λ. If S(z)+τE is invertible for all z∈Γ and all τ∈[0,1], then T+E has a simple eigenvalue λ~ inside Γ and (24)λ-λ~≤CE,where C essentially depends on maxz∈ΓS(z)-1|Ran(E) and maxz∈ΓS(z)-∗|Ran(E∗) and the length of the the integration curve Γ.Proof.
Recall that(Sz+E)′=S′(z) and define the function: (25)fτ=12πitr∫ΓS′zSz+τE-1dz.By the generalized argument principle—see [17–19]—the value of f(τ) equals the total multiplicity of the eigenvalues inside Γ. In particular, by Proposition 2 we have that there are vectors x and y such that (S′(λ)x,y)=y∗S′(λ)x=1 and (26)f0=trS′λxy∗=try∗S′λx=1,where we used the circularity of the trace. By the assumptions of the theorem S(z) is invertible for all z∈Γ and E is a rank one bounded operator. Therefore, there are vectors v and u so that E=uv∗ and using Sherman-Morrison formula (see [21–23] and the references therein), we write (27)Sz+τE-1=Sz-1-τ1+v∗S-1zuSz-1uv∗Sz-1,and so it follows that(28)Sz+τ1E-1-Sz+τ2E-1=τ2-τ11+v∗S-1zuSz-1uv∗Sz-1=τ2-τ11+v∗S-1zuSz-1uSz-∗v∗.And in particular f is a smooth function. Since f, due to Rouche’s theorem, conclude that f takes values only in natural numbers, it must be constant for all τ∈[0,1]. Let us denote this eigenvalue of S~=S+E by λ~. Define the operators (29)Aq=∫ΓzqS-1zdz,A~q=∫ΓzqS~-1zdz,q=0,1.Proposition 2 implies A(1)=λxy∗, A~(1)=λ~x~y~∗, where S(λ)x=0, S(λ)∗y=0, S~(λ~)x~=0, and S~(λ)∗y~=0 and so there exists a vector ϕ such that ϕ∗A(0)ϕ≠0 and ϕ∗A~(0)ϕ≠0. We now compute (30)λ=ϕ∗∫ΓzS-1zdzϕϕ∗∫ΓS-1zdzϕ,λ~=ϕ∗∫ΓzS~-1zdzϕϕ∗∫ΓS~-1zdzϕand so the second assertion follows by the following computation: (31)λ-λ~=ϕ∗∫ΓzS-1zdzϕϕ∗∫ΓS-1zdzϕ-ϕ∗∫ΓzS~-1zdzϕϕ∗∫ΓS~-1zdzϕ≤1ϕ∗∫ΓS-1zdzϕϕ∗∫ΓzS-1zdzϕ-ϕ∗∫ΓzS~-1zdzϕ+ϕ∗∫ΓzS~-1zdzϕϕ∗∫ΓS~-1zdzϕϕ∗∫ΓS~-1zdzϕ-ϕ∗∫ΓS-1zdzϕ≤1ϕ∗∫ΓS-1zdzϕϕ∗∫ΓzS-1z-zS~-1zdzϕ+λ~ϕ∗∫ΓS~-1z-S-1zdzϕ=1ϕ∗∫ΓS-1zdzϕϕ∗∫Γz1+v∗S-1zuSz-1ESz-1dzϕ+λ~ϕ∗∫Γ11+v∗S-1zuSz-1ESz-1dzϕ≤CE.
The claim on the constantC follows from the observation that Ran(E)=span{u} and Ran(E∗)=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.
Lett(·) denote the family of sesquilinear forms (11) with operator representation T(·):Ω→H given by (15). Assume that a contour Γ⊂ρ(T) encloses only a simple eigenvalue λ and μ∈C but no other points of σ(T). Let u∈V, with u=1. Then there is λ∈σ(T) such that (32)λ-μ≤Csupϕ∈V∖0⁡tμu,ϕt0ϕ,ϕ.The constant C does not depend on μ and u but on contour Γ and the restriction of T-1(·) and T-∗(·) to Γ.Proof.
We will now construct a particular operatorEμ,u which will be used to assess the quality of the pair μ, u. Define the sesquilinear form eμ,u:V×V→C by the formula (33)eμ,uψ,ϕ=-tμu,ϕψ,u,ψ,ϕ∈V.It is obviously relatively compact with respect to t0 and so we can define the Fredholm valued operator function Ω∋z→T~(z), where T~(z) is the operator defined by the form (34)t~λψ,ϕ=tλψ,ϕ+eμ,uψ,ϕ.Further, (35)T~μu,ϕ=tμu,ϕ+eμ,uu,ϕ=tμu,ϕ-tμu,ϕu,u=tμu,ϕ-tμu,ϕ=0,ϕ∈V,and so μ∈σ(T~). We construct the operator Eμ,u as the operator defined by (36)Eμ,uϕ,ψ=eμ,uT0-1/2ϕ,T0-1/2ψ,ϕ,ψ∈H.Now recall that (22) and (35) imply (37)T~-1z=T0-1/2I+∑j=1dfjzKj+Eμ,u-1T0-1/2.By construction the operator function S~:Ω→L(H)(38)S~z=I+∑j=1dfjzKj+Eμ,u=Sz+Eμ,usatisfies the assumption of Lemma 3. Finally, note that (39)supϕ∈V∖0⁡eμ,uμu,ϕt0ϕ,ϕ=supϕ∈V∖0⁡tμu,ϕt0ϕ,ϕ≕tμ,uμu,·t0,-1,and so from (36) it follows that(40)Eμ,u=tμ,uμu,·t0,-1.Now recall the definition of C from Lemma 3 to conclude the proof.


### 2.3. A Sketch for a Practical Algorithm for Error Estimation

Proposition4 will be the basis for practical error estimation. Let V⊃Qh, h>0 be a family of finite dimensional spaces such that the orthogonal projections Qh onto Qh converge strongly to I as h→∞. Let further Qh1⊂Qh2 for h1≤h2. Assume that ϕ∈Qh1 and μ∈C are given. For h2, h2>h1 define (41)tμu,·Qh2,t0,-1=supϕ∈Qh2∖0⁡tμu,ϕt0ϕ,ϕ.Obviously it holds, recalling the definition from (39), that (42)tμu,·Qh2,t0,-1≤tμ,uμu,·t0,-1.However, in the case in which Qh, h>0 satisfies the standard saturation property, for example, (I-Qh2)ϕ≤q(I-Qh1)ϕ for fixed q, 0<q<1 and some ϕ, we can also prove (43)tμu,·Qh2,t0,-1≤tμu,·t0,-1≤Ctμu,·Qh2,t0,-1.The constant C essentially depends on q which in turn depends on h2-h1 but not on the magnitude of h1; for more details see [24, 25].


## 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.Lett0 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(t0)⊂H to denote its domain of definition. To simplify the notation we will write V≔Dom(t0) in the rest of the paper. Additionally, let ci, i=1,…,d, be a sequence of sesquilinear (not necessarily symmetric) forms which are relatively compact [16] with regard to t0. Moreover, assume that fi·, i=1,…,d, is a sequence of scalar analytic functions in Ω. Define the family of sesquilinear forms: (11)tzϕ,ψ≔t0ϕ,ψ+∑i=1dfizciϕ,ψ,ϕ,ψ∈V,z∈Ω.In a variationally posed eigenvalue problem we are seeking a scalar λ∈Ω and a vector ϕ∈V∖{0} such that (12)tλϕ,ψ=0,ψ∈V.To this variational formulation we construct a representation by an operator-valued function T(·) and then apply the analytic Fredholm theorem to establish the structure of the spectrum.LetT0, Dom(T01/2)=V be a self-adjoint positive definite operator which represents the form t0 in the sense of Kato’s second representation theorem; see [16, Theorem VI.2.23]. Let Ki be defined by (13)Kiu,v=ciT0-1/2u,T0-1/2v,u,v∈H.The operators Ki, i=1,…,d, are obviously compact and for z∈Ω∖σ(T) the value of the generalized resolvent is the operator:(14)T-1z=T0-1/2I+∑j=1dfjzKj-1T0-1/2.Here T(z) is the unbounded sectorial operator with domain Dom(T(z)) such that (15)Tzϕ,ψ=t0ϕ,ψ+∑i=1dfizciϕ,ψ,ϕ∈DomTz,ψ∈V.Obviously the operator-valued function (16)Sz=I+∑j=1dfjzKj,z∈Ωsatisfies the requirements of the analytic Fredholm theorem and σ(S)=σ(T). Let us note that we will use (17)T′zϕ,ψ=∑i=1dfi′zciϕ,ψ,ϕ∈DomT′z,ψ∈V,to define a derivative of an operator-valued analytic function. We can now define the notion of a semisimple eigenvalue.Definition 1.
LetT be as in (15) and let μ∈σ(T) be an eigenvalue. The eigenvalue μ is semisimple if for each ψ∈Ker(T(μ))∖{0} there is a ψ∈Ker(T∗(μ)) such that (18)T′μϕ,ψ≠0.If dim⁡Ker⁡(T(μ))=1, then μ is called a simple eigenvalue.Note that in the quasi-matrix notation we will freely write(19)T′μϕ,ψ=ψ∗T′μϕ.To this end we identify the vectors with a mapping ψ:C→H.With these conventions we state, informally, the generalized argument principle proved by Gohberg and Sigal in [17–19]. It states that for the closed contour Γ⊂ρ(T) the number (20)MT,Γ=tr∫ΓTz-1T′zdz=tr∫ΓT′zTz-1dzsatisfies M(T,Γ)∈N 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 T:Ω→H given by (15). Then the spectrum σ(T) consists of a countable collection of eigenvalues with finite multiplicity. Further, let the component of σ(T) inside a contour Γ consist only of semisimple eigenvalues λi, i=1,…,r, such that nλi=dim⁡Ker(T(λi)). Then there are quasi-matrices Xi,Yi:Cnλi→H, i=1,…,r, such that Ran(Xi)⊂Ker(T(λi)) and Ran(Yi)⊂Ker(T∗(λi)), i=1,…,r, and an open once connected neighborhood U containing λi, i=1,…,r and Γ and an operator-valued function H which is analytic on U and taking values in L(H) such that (21)T-1z=∑i=1r1z-λiXiYi∗+Hz,z∈U,and Yi∗(T′(λi)Xi)=Inλi, i=1,…,r.Proof.
Forz∈Ω∖σ(T) write (22)T-1z=T0-1/2I+∑j=1dfjzKj-1T0-1/2and define the Fredholm valued function: (23)Sz=I+∑j=1dfjzKj,z∈Ω.Recall that Dom(T01/2)=V and that T01/2 maps V one to one on H. Now apply [9, Theorem 1.3.1] on S.


## 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.
LetS:Ω→Φ(H) be an analytic function and let E be a bounded operator such that dim⁡Ran(E)=1. Assume that Γ⊂ρ(S) is a simple closed contour such that the component of σ(S) inside Γ consists solely of a simple eigenvalue λ. If S(z)+τE is invertible for all z∈Γ and all τ∈[0,1], then T+E has a simple eigenvalue λ~ inside Γ and (24)λ-λ~≤CE,where C essentially depends on maxz∈ΓS(z)-1|Ran(E) and maxz∈ΓS(z)-∗|Ran(E∗) and the length of the the integration curve Γ.Proof.
Recall that(Sz+E)′=S′(z) and define the function: (25)fτ=12πitr∫ΓS′zSz+τE-1dz.By the generalized argument principle—see [17–19]—the value of f(τ) equals the total multiplicity of the eigenvalues inside Γ. In particular, by Proposition 2 we have that there are vectors x and y such that (S′(λ)x,y)=y∗S′(λ)x=1 and (26)f0=trS′λxy∗=try∗S′λx=1,where we used the circularity of the trace. By the assumptions of the theorem S(z) is invertible for all z∈Γ and E is a rank one bounded operator. Therefore, there are vectors v and u so that E=uv∗ and using Sherman-Morrison formula (see [21–23] and the references therein), we write (27)Sz+τE-1=Sz-1-τ1+v∗S-1zuSz-1uv∗Sz-1,and so it follows that(28)Sz+τ1E-1-Sz+τ2E-1=τ2-τ11+v∗S-1zuSz-1uv∗Sz-1=τ2-τ11+v∗S-1zuSz-1uSz-∗v∗.And in particular f is a smooth function. Since f, due to Rouche’s theorem, conclude that f takes values only in natural numbers, it must be constant for all τ∈[0,1]. Let us denote this eigenvalue of S~=S+E by λ~. Define the operators (29)Aq=∫ΓzqS-1zdz,A~q=∫ΓzqS~-1zdz,q=0,1.Proposition 2 implies A(1)=λxy∗, A~(1)=λ~x~y~∗, where S(λ)x=0, S(λ)∗y=0, S~(λ~)x~=0, and S~(λ)∗y~=0 and so there exists a vector ϕ such that ϕ∗A(0)ϕ≠0 and ϕ∗A~(0)ϕ≠0. We now compute (30)λ=ϕ∗∫ΓzS-1zdzϕϕ∗∫ΓS-1zdzϕ,λ~=ϕ∗∫ΓzS~-1zdzϕϕ∗∫ΓS~-1zdzϕand so the second assertion follows by the following computation: (31)λ-λ~=ϕ∗∫ΓzS-1zdzϕϕ∗∫ΓS-1zdzϕ-ϕ∗∫ΓzS~-1zdzϕϕ∗∫ΓS~-1zdzϕ≤1ϕ∗∫ΓS-1zdzϕϕ∗∫ΓzS-1zdzϕ-ϕ∗∫ΓzS~-1zdzϕ+ϕ∗∫ΓzS~-1zdzϕϕ∗∫ΓS~-1zdzϕϕ∗∫ΓS~-1zdzϕ-ϕ∗∫ΓS-1zdzϕ≤1ϕ∗∫ΓS-1zdzϕϕ∗∫ΓzS-1z-zS~-1zdzϕ+λ~ϕ∗∫ΓS~-1z-S-1zdzϕ=1ϕ∗∫ΓS-1zdzϕϕ∗∫Γz1+v∗S-1zuSz-1ESz-1dzϕ+λ~ϕ∗∫Γ11+v∗S-1zuSz-1ESz-1dzϕ≤CE.
The claim on the constantC follows from the observation that Ran(E)=span{u} and Ran(E∗)=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.
Lett(·) denote the family of sesquilinear forms (11) with operator representation T(·):Ω→H given by (15). Assume that a contour Γ⊂ρ(T) encloses only a simple eigenvalue λ and μ∈C but no other points of σ(T). Let u∈V, with u=1. Then there is λ∈σ(T) such that (32)λ-μ≤Csupϕ∈V∖0⁡tμu,ϕt0ϕ,ϕ.The constant C does not depend on μ and u but on contour Γ and the restriction of T-1(·) and T-∗(·) to Γ.Proof.
We will now construct a particular operatorEμ,u which will be used to assess the quality of the pair μ, u. Define the sesquilinear form eμ,u:V×V→C by the formula (33)eμ,uψ,ϕ=-tμu,ϕψ,u,ψ,ϕ∈V.It is obviously relatively compact with respect to t0 and so we can define the Fredholm valued operator function Ω∋z→T~(z), where T~(z) is the operator defined by the form (34)t~λψ,ϕ=tλψ,ϕ+eμ,uψ,ϕ.Further, (35)T~μu,ϕ=tμu,ϕ+eμ,uu,ϕ=tμu,ϕ-tμu,ϕu,u=tμu,ϕ-tμu,ϕ=0,ϕ∈V,and so μ∈σ(T~). We construct the operator Eμ,u as the operator defined by (36)Eμ,uϕ,ψ=eμ,uT0-1/2ϕ,T0-1/2ψ,ϕ,ψ∈H.Now recall that (22) and (35) imply (37)T~-1z=T0-1/2I+∑j=1dfjzKj+Eμ,u-1T0-1/2.By construction the operator function S~:Ω→L(H)(38)S~z=I+∑j=1dfjzKj+Eμ,u=Sz+Eμ,usatisfies the assumption of Lemma 3. Finally, note that (39)supϕ∈V∖0⁡eμ,uμu,ϕt0ϕ,ϕ=supϕ∈V∖0⁡tμu,ϕt0ϕ,ϕ≕tμ,uμu,·t0,-1,and so from (36) it follows that(40)Eμ,u=tμ,uμu,·t0,-1.Now recall the definition of C from Lemma 3 to conclude the proof.


## 2.3. A Sketch for a Practical Algorithm for Error Estimation

Proposition4 will be the basis for practical error estimation. Let V⊃Qh, h>0 be a family of finite dimensional spaces such that the orthogonal projections Qh onto Qh converge strongly to I as h→∞. Let further Qh1⊂Qh2 for h1≤h2. Assume that ϕ∈Qh1 and μ∈C are given. For h2, h2>h1 define (41)tμu,·Qh2,t0,-1=supϕ∈Qh2∖0⁡tμu,ϕt0ϕ,ϕ.Obviously it holds, recalling the definition from (39), that (42)tμu,·Qh2,t0,-1≤tμ,uμu,·t0,-1.However, in the case in which Qh, h>0 satisfies the standard saturation property, for example, (I-Qh2)ϕ≤q(I-Qh1)ϕ for fixed q, 0<q<1 and some ϕ, we can also prove (43)tμu,·Qh2,t0,-1≤tμu,·t0,-1≤Ctμu,·Qh2,t0,-1.The constant C essentially depends on q which in turn depends on h2-h1 but not on the magnitude of h1; 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 T(·) be given and let Γ⊂Ω be a closed curve which encloses either a set of r simple eigenvalues or a single semisimple eigenvalue whose multiplicity is r.Let us consider the operators(44)Aq=12πi∫ΓzqTz-1dz,Bq=12πi∫ΓzqTz-1T′zdz,q=0,1,and their approximations: (45)ANq=∑k=1NωkzkqTzk-1,BNq=∑k=1NωkzkqTzk-1T′zk,q=0,1.Here ωi and zi∈Γ, i=1,…,N, are integration weights and integration nodes. Based on Theorem  4.7 in [4] we establish, for ωi and zi∈Γ defined by the N node trapeze integration formula for the contour integral (44), the following estimates: (46)Aq-ANq≤C1dT-κe-C2NdT,Bq-BNq≤C1dT-κe-C2NdT.Here Γ is simple closed contour in Ω such that σ(T)∩Γ=∅, d(T)=minλ∈σ(T)dist(λ,Γ), and κ is the maximum order of poles for the inverse of T. 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 AN(q)→A(q) and BN(q)→B(q), q=0,1 in the norm resolvent sense.


### 3.1. Extracting Information on Eigenvalues Enclosed by a Contour

Based on Proposition2 we see that operators A(0) and B(0) are finite rank operators such that Ran(A(0))=Ran(B(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 A(1) and B(1).Proposition 5.
LetT be the operator-valued function from (15) and let Γ be the contour which encloses solely the r, counting according to algebraic multiplicity, semisimple eigenvalues of T. Define the matrices A(q) and B(q), q=0,1, as in (44) and let Q:Cr→H~, Q∗Q=I, be a quasi-matrix such that Q∗A(0)Q and Q∗B(0)Q are invertible. Then the eigenvalues of the matrix pairs (Q∗A(1)Q,Q∗A(0)Q) and (Q∗B(1)Q,Q∗B(0)Q) are precisely the eigenvalues λi, i=1,…,r, where we count according to multiplicity. Furthermore, if ai,bi∈Cr satisfy Q∗A(1)Qai=λiQ∗A(0)Qbi and Q∗B(1)Qai=λiQ∗B(0)Qbi, then Qai and Qbi are eigenvectors of T associated with λi, i=1,…,r.Proof.
We will prove the statement for the matrix pair(Q∗B(1)Q,Q∗B(0)Q) and note that the proof for the matrix pair (Q∗A(1)Q,Q∗A(0)Q) is equivalent. Based on Proposition 2 we see that there are quasi-matrices X:Cr→H~ and Y:Cr→H~ and a neighborhood U of Γ such that (47)Tz-1=Xz-λ1-10⋯00z-λ2-1⋯0⋯⋯⋯⋯0⋯0z-λr-1Y∗+Hz,z∈U.Here H is an analytic operator-valued function. Now we compute (48)Q∗B1Q=Q∗Xλ10⋯00λ2⋯0⋯⋯⋯⋯0⋯0λrY∗QQ∗B0Q=Q∗XY∗Q.Based on the assumptions of the theorem we conclude that (Q∗X) and (Y∗Q) have to be invertible and so the conclusion readily follows.Before we proceed, note that norm resolvent convergence ofAN(q) and BN(q) implies the convergence of spectra and associated spaces to those of A(q) and B(q), q=0,1. In particular, this means that AN(0) and BN(0), for N 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 AN(0)v and BN(0)v 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 operatorB of type (44) which has a component of r, r∈N dominant eigenvalues and whose action on a vector can be represented by a formula (45).For Algorithm1 we present standard convergence results in Theorems 6 and 8. First we consider a special case in which B 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.Algorithm 1:Contour based inexact subspace iteration for B, tol, and N.
Pick a random quasi-matrixQ1=q1q2⋯qr, Q1:Cr→L2(Ω)
Setn=1
repeat
ComputeYn+1=BNQn
Compute an orthogonal quasi-matrixQn+1, Qn+1∗Qn+1=I using Gram-Schmidt so that Ran(Qn+1)=Ran(Yn+1).
Computecriterion = ‖Qn+1-Qn(Qn∗Qn+1)‖
n≔n+1
until  criterion < tol
return  Qn+1Algorithm 2:Approximating eigenvectors and eigenvalues with tolerance tol2.
Require  Q1:Cr→H~, Q1∗Q1=Ir, tol1 and N.
ComputeQ as output from Algorithm 1 with tolerance tol1 and Q1 as starting value.
ComputeZ1=BN(1)Q and Z0=BN(0)Q.
FormW(q)=Q∗Zq, q=0,1
Computevi and λi so that W(1)ui=λiW(0)ui, i=1,…,r
Compute eigenvectorsvi=Qn+1ui, i=1,…,r.
Compute‖t(λi)[vi,⋅]‖Qh2,t0,-1
Form the index seti1,…,ir2=i:‖t(λi)[vi,⋅]‖Qh2,t0,-1≤tol2
SetV=vi1⋯vir2
return  V:Cr2→H~ and λij, j=1,…,r2Theorem 6.
LetB be a bounded Hermitian operator such that σ(B)={λ1,…,λr}∪Σ, where λi, i=1,…,r, are eigenvalues of finite multiplicity and |λ1|≥⋯≥|λr|>τ>0, and τ>|λ|, λ∈Σ. Then for the sequence V(n+1)≔BV(n), where V(0):Cr→H, (V(0))∗V(0)=Ir, we have (49)tan⁡θVn+1,U≤τnλrntan⁡θV0,U,where U:Cr→H, U∗U=I, is such that U=u1u2⋯ur and Bui=λiui, i=1,…,r.Proof.
LetP1=UU∗ and represent U:Cr→H with respect to the partition of unity P1⊕P2=I by U=U1U2, with Ran(P1)=Ran(U1). Since B is self-adjoint it can be represented by the block operator matrix: (50)B=B1100B22.Let us represent, with respect to the same partition of unity, the quasi-matrix V(0) as (51)V0=V10V20.Without reducing the level of generality we may assume that V1(0) is invertible, since otherwise Ran(V(0)) would contain eigenvectors of B. This corresponds to the trivial situation and will not be further discussed.
Since(V(0))∗V(0)=Ir it follows from (6) that (52)tan⁡θU1,V0=V20V10-1.The assumption of the separation of the spectra of B implies that B11 must also be invertible, and so (53)Vn≔BnV0=B11nV10B22kV20=IFnB11nV10for F(n)=B22nV2(0)(V1(0))-1B11-n which implies that(54)VnV10-1B11-n=IFn=U1+U2Fn.And so from the optimality of the canonical angles (see [15]), it follows that(55)tan⁡θVn,U≤Fn.Consequently, it follows that(56)tan⁡θVn,U≤Fn≤B22nF0B11-n≤τnλrntan⁡θV0,U1.Now the conclusion of the theorem readily follows.Remark 7.
Note that from (8) and (9) it is clear that θ(V(n),U) does not depend on the quasi-matrices V(n) and U, but rather on the orthogonal projections onto Ran(V(n)) and Ran(U)=Ran(P1). In Theorem 6 we have stated the convergence result for the angle θ(V(n),U). Here we have ignored the choice of a basis of Ran(V(n)) 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 operatorB 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.
LetB be a bounded operator such that its spectrum has two disjoint components Σ1 and Σ2 so that σ(B)=Σ1∪Σ2. Let further Σ1={λ1,λ2,…,λr}, where we have counted λi∈σ(B) according to the algebraic multiplicity of an eigenvalue. Further let a partition of unity P1⊕P2=I exist such that r=dim⁡(Ran(P1)) and (57)B=B11B120B22,and let further sep(B11,B22)>1 and B22B11-1<1. Then for a quasi-matrix V:Cr→H we have (58)tan⁡θBnV,P1≤B22B11-1n1+Gtan⁡θV,P1.Here G is the solution of the Sylvester equation GB22-B11G=B12 and sep(B11,B22)=infG=1GB22-B11G.Proof.
ForB such that σ(B)=Σ1∪Σ2 with Σ1={λ1,λ2,…,λr} satisfying the assumptions of the theorem choose P1 as the orthogonal projection onto the maximal invariant subspace associated with {λ1,λ2,…,λr} and P2=I-P1; then B has the block operator representation as claimed in the theorem. We note that (59)B=B11B120B22=I-G0IB1100B22IG0I,where GB22-B11G=B12. Operator G, such that GB22-B11G=B12 exists and is unique, for example, [29, Lemma A.1], providing the spectra of B11 and B22, 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 sep(B11,B22) 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 nodesN diminishes exponentially. To extract eigenvalues and eigenvectors we use the method from [7] and apply Proposition 5 on the quasi-matrix Qn+1 which is returned by Algorithm 1. More precisely we have the following. For simplicity we write B=B(0).Theorem 10.
Let the conditions of Theorem8 be satisfied and also assume the same notation conventions. For a quasi-matrix V:Cr→H, we set Q1≔V and apply Algorithm 1 with fixed N∈N. Then (60)sin⁡θRanQn,P1≤B22B11-1+C3dT-κe-C2NdTn1+G+C4dT-κe-C2NdTtan⁡θQ1,P1+C5dT-κe-C2NdT.Proof.
Recall that(61)B-BN≤C1dT-κe-C2NdT.From standard results on norm convergence of bounded operators (see [16]) we conclude that (62)sin⁡θP1,PN1=P1-PN1≤C5dT-κe-C2NdT.Here (PN)1+(PN)2=I is the partition of unity in which BN has a block triangular form: (63)BN=BN11BN120BN22.Assumption on the spectral separation from Theorem 8 implies that B11 is invertible and obviously (64)B11-BN11≤C1dT-κe-C2NdT(65)B22-BN22≤C1dT-κe-C2NdT.Norm convergence (64) implies that for N large enough (BN)11 must also be invertible; for details see [30]. For such N we have (66)B11-1-BN11-1≤B11-1B11-BN11BN11-1≤B11-1BN11-1B11-BN11≤C3dT-κe-C2NdT.Triangle inequality implies that(67)BN22≤B22+BN22-B22,BN11-1≤B11-1+BN11-1-B11-1and so (68)BN22BN11-1≤B22B11-1+C4dT-κe-C2NdT.
Let us now assume thatN is large enough so that (69)B22B11-1+C4dT-κe-C2NdT<1.Then we can apply Theorem 8 to BN and conclude that (70)sin⁡θBNnV,PN1≤tan⁡θBNnV,PN1≤BN22BN11-1n1+GNtan⁡θV,PN1.
Equivalently sincesep((BN)22,(BN)11)→sep(B22,B11) and (BN)12→B12 in norm we conclude that (71)GN≤G+C5dT-κe-C2NdT.We now apply the triangle inequality for the sine of the maximal canonical angle to conclude the proof.Using Theorem10 we will define theeffective convergence rate of the inexact contour based subspace iteration. Recall that (72)B=12πi∫ΓTz-1T′zdz,BN=∑k=1NωkzkqTzk-1T′zk.Starting from (63) we define theeffective convergence rate for the inexact subspace iteration for B as(73)ηΓ,N≔BN22BN11-1.Remark 11.
We will projectBN on carefully constructed finite dimensional spaces, for example, finite element spaces, to experimentally estimate η(Γ,N) 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 inequalityt(μ)[u,·]Qh2,t0,-1≤t(μ)[u,·]t0,-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(μ)[u,·]Qh2,t0,-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 Algorithm2 is not satisfactory, the procedure can be extended as an innerouter iteration scheme by setting Q1=V 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 Q1 which are identified by Algorithm 2 as converged and then proceed as in the original versions of Algorithm 1.


## 3.1. Extracting Information on Eigenvalues Enclosed by a Contour

Based on Proposition2 we see that operators A(0) and B(0) are finite rank operators such that Ran(A(0))=Ran(B(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 A(1) and B(1).Proposition 5.
LetT be the operator-valued function from (15) and let Γ be the contour which encloses solely the r, counting according to algebraic multiplicity, semisimple eigenvalues of T. Define the matrices A(q) and B(q), q=0,1, as in (44) and let Q:Cr→H~, Q∗Q=I, be a quasi-matrix such that Q∗A(0)Q and Q∗B(0)Q are invertible. Then the eigenvalues of the matrix pairs (Q∗A(1)Q,Q∗A(0)Q) and (Q∗B(1)Q,Q∗B(0)Q) are precisely the eigenvalues λi, i=1,…,r, where we count according to multiplicity. Furthermore, if ai,bi∈Cr satisfy Q∗A(1)Qai=λiQ∗A(0)Qbi and Q∗B(1)Qai=λiQ∗B(0)Qbi, then Qai and Qbi are eigenvectors of T associated with λi, i=1,…,r.Proof.
We will prove the statement for the matrix pair(Q∗B(1)Q,Q∗B(0)Q) and note that the proof for the matrix pair (Q∗A(1)Q,Q∗A(0)Q) is equivalent. Based on Proposition 2 we see that there are quasi-matrices X:Cr→H~ and Y:Cr→H~ and a neighborhood U of Γ such that (47)Tz-1=Xz-λ1-10⋯00z-λ2-1⋯0⋯⋯⋯⋯0⋯0z-λr-1Y∗+Hz,z∈U.Here H is an analytic operator-valued function. Now we compute (48)Q∗B1Q=Q∗Xλ10⋯00λ2⋯0⋯⋯⋯⋯0⋯0λrY∗QQ∗B0Q=Q∗XY∗Q.Based on the assumptions of the theorem we conclude that (Q∗X) and (Y∗Q) have to be invertible and so the conclusion readily follows.Before we proceed, note that norm resolvent convergence ofAN(q) and BN(q) implies the convergence of spectra and associated spaces to those of A(q) and B(q), q=0,1. In particular, this means that AN(0) and BN(0), for N 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 AN(0)v and BN(0)v 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 operatorB of type (44) which has a component of r, r∈N dominant eigenvalues and whose action on a vector can be represented by a formula (45).For Algorithm1 we present standard convergence results in Theorems 6 and 8. First we consider a special case in which B 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.Algorithm 1:Contour based inexact subspace iteration for B, tol, and N.
Pick a random quasi-matrixQ1=q1q2⋯qr, Q1:Cr→L2(Ω)
Setn=1
repeat
ComputeYn+1=BNQn
Compute an orthogonal quasi-matrixQn+1, Qn+1∗Qn+1=I using Gram-Schmidt so that Ran(Qn+1)=Ran(Yn+1).
Computecriterion = ‖Qn+1-Qn(Qn∗Qn+1)‖
n≔n+1
until  criterion < tol
return  Qn+1Algorithm 2:Approximating eigenvectors and eigenvalues with tolerance tol2.
Require  Q1:Cr→H~, Q1∗Q1=Ir, tol1 and N.
ComputeQ as output from Algorithm 1 with tolerance tol1 and Q1 as starting value.
ComputeZ1=BN(1)Q and Z0=BN(0)Q.
FormW(q)=Q∗Zq, q=0,1
Computevi and λi so that W(1)ui=λiW(0)ui, i=1,…,r
Compute eigenvectorsvi=Qn+1ui, i=1,…,r.
Compute‖t(λi)[vi,⋅]‖Qh2,t0,-1
Form the index seti1,…,ir2=i:‖t(λi)[vi,⋅]‖Qh2,t0,-1≤tol2
SetV=vi1⋯vir2
return  V:Cr2→H~ and λij, j=1,…,r2Theorem 6.
LetB be a bounded Hermitian operator such that σ(B)={λ1,…,λr}∪Σ, where λi, i=1,…,r, are eigenvalues of finite multiplicity and |λ1|≥⋯≥|λr|>τ>0, and τ>|λ|, λ∈Σ. Then for the sequence V(n+1)≔BV(n), where V(0):Cr→H, (V(0))∗V(0)=Ir, we have (49)tan⁡θVn+1,U≤τnλrntan⁡θV0,U,where U:Cr→H, U∗U=I, is such that U=u1u2⋯ur and Bui=λiui, i=1,…,r.Proof.
LetP1=UU∗ and represent U:Cr→H with respect to the partition of unity P1⊕P2=I by U=U1U2, with Ran(P1)=Ran(U1). Since B is self-adjoint it can be represented by the block operator matrix: (50)B=B1100B22.Let us represent, with respect to the same partition of unity, the quasi-matrix V(0) as (51)V0=V10V20.Without reducing the level of generality we may assume that V1(0) is invertible, since otherwise Ran(V(0)) would contain eigenvectors of B. This corresponds to the trivial situation and will not be further discussed.
Since(V(0))∗V(0)=Ir it follows from (6) that (52)tan⁡θU1,V0=V20V10-1.The assumption of the separation of the spectra of B implies that B11 must also be invertible, and so (53)Vn≔BnV0=B11nV10B22kV20=IFnB11nV10for F(n)=B22nV2(0)(V1(0))-1B11-n which implies that(54)VnV10-1B11-n=IFn=U1+U2Fn.And so from the optimality of the canonical angles (see [15]), it follows that(55)tan⁡θVn,U≤Fn.Consequently, it follows that(56)tan⁡θVn,U≤Fn≤B22nF0B11-n≤τnλrntan⁡θV0,U1.Now the conclusion of the theorem readily follows.Remark 7.
Note that from (8) and (9) it is clear that θ(V(n),U) does not depend on the quasi-matrices V(n) and U, but rather on the orthogonal projections onto Ran(V(n)) and Ran(U)=Ran(P1). In Theorem 6 we have stated the convergence result for the angle θ(V(n),U). Here we have ignored the choice of a basis of Ran(V(n)) 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 operatorB 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.
LetB be a bounded operator such that its spectrum has two disjoint components Σ1 and Σ2 so that σ(B)=Σ1∪Σ2. Let further Σ1={λ1,λ2,…,λr}, where we have counted λi∈σ(B) according to the algebraic multiplicity of an eigenvalue. Further let a partition of unity P1⊕P2=I exist such that r=dim⁡(Ran(P1)) and (57)B=B11B120B22,and let further sep(B11,B22)>1 and B22B11-1<1. Then for a quasi-matrix V:Cr→H we have (58)tan⁡θBnV,P1≤B22B11-1n1+Gtan⁡θV,P1.Here G is the solution of the Sylvester equation GB22-B11G=B12 and sep(B11,B22)=infG=1GB22-B11G.Proof.
ForB such that σ(B)=Σ1∪Σ2 with Σ1={λ1,λ2,…,λr} satisfying the assumptions of the theorem choose P1 as the orthogonal projection onto the maximal invariant subspace associated with {λ1,λ2,…,λr} and P2=I-P1; then B has the block operator representation as claimed in the theorem. We note that (59)B=B11B120B22=I-G0IB1100B22IG0I,where GB22-B11G=B12. Operator G, such that GB22-B11G=B12 exists and is unique, for example, [29, Lemma A.1], providing the spectra of B11 and B22, 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 sep(B11,B22) 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 nodesN diminishes exponentially. To extract eigenvalues and eigenvectors we use the method from [7] and apply Proposition 5 on the quasi-matrix Qn+1 which is returned by Algorithm 1. More precisely we have the following. For simplicity we write B=B(0).Theorem 10.
Let the conditions of Theorem8 be satisfied and also assume the same notation conventions. For a quasi-matrix V:Cr→H, we set Q1≔V and apply Algorithm 1 with fixed N∈N. Then (60)sin⁡θRanQn,P1≤B22B11-1+C3dT-κe-C2NdTn1+G+C4dT-κe-C2NdTtan⁡θQ1,P1+C5dT-κe-C2NdT.Proof.
Recall that(61)B-BN≤C1dT-κe-C2NdT.From standard results on norm convergence of bounded operators (see [16]) we conclude that (62)sin⁡θP1,PN1=P1-PN1≤C5dT-κe-C2NdT.Here (PN)1+(PN)2=I is the partition of unity in which BN has a block triangular form: (63)BN=BN11BN120BN22.Assumption on the spectral separation from Theorem 8 implies that B11 is invertible and obviously (64)B11-BN11≤C1dT-κe-C2NdT(65)B22-BN22≤C1dT-κe-C2NdT.Norm convergence (64) implies that for N large enough (BN)11 must also be invertible; for details see [30]. For such N we have (66)B11-1-BN11-1≤B11-1B11-BN11BN11-1≤B11-1BN11-1B11-BN11≤C3dT-κe-C2NdT.Triangle inequality implies that(67)BN22≤B22+BN22-B22,BN11-1≤B11-1+BN11-1-B11-1and so (68)BN22BN11-1≤B22B11-1+C4dT-κe-C2NdT.
Let us now assume thatN is large enough so that (69)B22B11-1+C4dT-κe-C2NdT<1.Then we can apply Theorem 8 to BN and conclude that (70)sin⁡θBNnV,PN1≤tan⁡θBNnV,PN1≤BN22BN11-1n1+GNtan⁡θV,PN1.
Equivalently sincesep((BN)22,(BN)11)→sep(B22,B11) and (BN)12→B12 in norm we conclude that (71)GN≤G+C5dT-κe-C2NdT.We now apply the triangle inequality for the sine of the maximal canonical angle to conclude the proof.Using Theorem10 we will define theeffective convergence rate of the inexact contour based subspace iteration. Recall that (72)B=12πi∫ΓTz-1T′zdz,BN=∑k=1NωkzkqTzk-1T′zk.Starting from (63) we define theeffective convergence rate for the inexact subspace iteration for B as(73)ηΓ,N≔BN22BN11-1.Remark 11.
We will projectBN on carefully constructed finite dimensional spaces, for example, finite element spaces, to experimentally estimate η(Γ,N) 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 inequalityt(μ)[u,·]Qh2,t0,-1≤t(μ)[u,·]t0,-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(μ)[u,·]Qh2,t0,-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 Algorithm2 is not satisfactory, the procedure can be extended as an innerouter iteration scheme by setting Q1=V 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 Q1 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δh of an interval [a,b]. We denote this space by Vh(a,b). 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 useQh(a,b) to denote the space of piecewise quadratic and continuous function on the same subdivision of the interval [a,b] which was used to define the space Vh(a,b). To estimate the error Eμ,u from (39) for the function u∈Vh(a,b) and μ∈C we use the formula; see [24] for further information: (74)tμu,·Qha,b,t0,-1=supϕ∈Qha,b∖0⁡tμu,ϕt0ϕ,ϕ.Example 12.
We study the quadratic eigenvalue problem:(75)∂xxu=λγ∂xxu+δu+λ2uu-2=u2=0for γ=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 BN with N=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 (BN)11 and (BN)22 can be achieved with few integration nodes N: for example, N=25. Figure 2(a) depicts the dependence of the effective convergence rate for the inexact subspace iteration η(Γ,N)=(BN)22(BN)11-1 on the number of integration nodes N for the trapeze formula. The contour was a circle which enclosed the five eigenvalues marked by crosses.Figure 1
Eigenvalues for the quadratic eigenvalue problem.Figure 2
The dependence of the effective convergence rateη(Γ,N) as defined in (73) on the number of integration nodes N. This illustrates the negative effect of the distance to Γ on the convergence rate.
(a)
Quadratic eigenvalue problem with well-separated spectrum from Example12 and Γ as in Figure 1
(b)
Rational eigenvalue problem with tightly clustered eigenvalues from Example14 and Γ as in Figure 4Remark 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:(76)d2dx2u+λ1-λχx≤9/10·u+λ3-λχx≤9/10·u+λχx≤9/10u+10χx>9/10u=0u∈H1-1,1with  periodic  b.c.,where(77)
χ
x
≤
9
/
10
:
x
⟼
0
:
x
∉
-
9
10
,
9
10
1
:
x
∈
-
9
10
,
9
10is the characteristic function of the interval -9/10,9/10.This eigenvalue problem has a singularity atλ=1 and λ=3. The eigenvalues below the singularity λ=1 and those in 〈1,3〉 are discrete (finite multiplicity) and they accumulate at λ=1 and at λ=3 from below [31].The results in Figures3 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.Figure 3
Eigenvectors for eigenvalues in[0,1〉.Figure 4
Eigenvectors for eigenvalues in〈1,3〉.As a first experiment in Table1 we present the empirical study of the convergence of the residual measures and the relative eigenvalue error computed from the sequence V(n+1)=BNV(n), V0:C2→V, where N=16 and ωi and zi 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.Table 1
Convergence history for the rational eigenvalue problem Example14 and the two leftmost eigenvalues from Figure 3.
Iteration numbern
Maximal eigenvalue error
Maximal residual estimate
1
2.0880
e
-
008
1.2607
e
-
009
2
4.9292
e
-
009
1.3579
e
-
010
3
4.9282
e
-
009
1.3580
e
-
010The 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 z=1. In Figure 2(b) we see that many more integration nodes are necessary to achieve a similar estimate of the effective convergence rate η(Γ,N) 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 Figure5 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(78)tμu,·Qha,b,t0,-1≤tol2in 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 N=16 and 3 steps of inexact subspace iteration acceleration we reach the same accuracy (see Table 1) as a priori predicted for N=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.Figure 5
We displayQn+1-Qn(Qn∗Qn+1) for n=1,…,10.
(a)
Eigenvalues clustering towardsΓ
(b)
Decay ofQn+1-Qn(Qn∗Qn+1) for the rational eigenvalue problem and N=350


### 4.1. Stability of Convergence Rate of Inexact Subspace Iteration

From Figure2 we see that the convergence rate for the subspace iteration essentially depends on B22B11-1≈η(Γ,N). We will now see that this convergence rate is relatively robust to perturbations as long as the distance from the curve Γ and σ(T) is not too small. To this end we consider the problem (79)∂xxu=λγ∂xxu+δu+λ2u+14exp⁡-λu-2=u2=0as 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 Qn+1-Qn(Qn∗Qn+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 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.Figure 6
Eigenvalues for the eigenvalue problem (79).Figure 7
We displayQn+1-Qn(Qn∗Qn+1) for n=1,…,10. For Γ we have used the contours as in Figures 1, 6, and 5(a), respectively.
(a)
Problem (75) and N=25
(b)
Problem (79) and N=25
(c)
Problem (76) and N=350


## 4.1. Stability of Convergence Rate of Inexact Subspace Iteration

From Figure2 we see that the convergence rate for the subspace iteration essentially depends on B22B11-1≈η(Γ,N). We will now see that this convergence rate is relatively robust to perturbations as long as the distance from the curve Γ and σ(T) is not too small. To this end we consider the problem (79)∂xxu=λγ∂xxu+δu+λ2u+14exp⁡-λu-2=u2=0as 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 Qn+1-Qn(Qn∗Qn+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 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.Figure 6
Eigenvalues for the eigenvalue problem (79).Figure 7
We displayQn+1-Qn(Qn∗Qn+1) for n=1,…,10. For Γ we have used the contours as in Figures 1, 6, and 5(a), respectively.
(a)
Problem (75) and N=25
(b)
Problem (79) and N=25
(c)
Problem (76) and N=350


## 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.



---

*Source: 459895-2015-11-02.xml*