FEAST for differential eigenvalue problems

TL;DR

This paper introduces an operator analogue of the FEAST eigensolver for differential operators, enabling precise computation of targeted spectral regions, especially high-frequency modes, with robustness and structure preservation.

Contribution

It develops a novel operator-based FEAST method with a rational filter for unbounded regions, preserving operator structure and achieving high-precision eigenvalue computation.

Findings

Accurately computes eigenvalues in specified spectral regions.

Handles unbounded search regions with a new rational filter.

Preserves operator structure for normal or self-adjoint operators.

Abstract

An operator analogue of the FEAST matrix eigensolver is developed to compute the discrete part of the spectrum of a differential operator in a region of interest in the complex plane. Unbounded search regions are handled with a novel rational filter for the right half-plane. If the differential operator is normal or self-adjoint, then the operator analogue preserves that structure and robustly computes eigenvalues to near machine precision accuracy. The algorithm is particularly adept at computing high-frequency modes of differential operators that possess self-adjoint structure with respect to weighted Hilbert spaces.

Equations

$$\mathcal{L}u=\lambda u,u(\pm 1)=\cdots =u^{(N/2)}(\pm 1)=0.$$

$$-\frac{d^{2}u}{d{x}^2}=\lambda u,u(\pm 1)=0.$$

$$P_{\mathcal{V}}=\frac{1}{2\pi i}{\int }_{\partial \Omega }(zI-A{)}^{-1}dz.$$

$$P_{\mathcal{V}}Y\approx \frac{1}{2\pi i}\sum _{k=1}^{\ell }{w}_k(z_{k}I-A{)}^{-1}Y.$$

$${\mathcal{P}}_{\mathcal{V}}=\frac{1}{2\pi i}{\int }_{\partial \Omega }(z\mathcal{I}-\mathcal{L}{)}^{-1}dz.$$

$${\kappa }_{\mathcal{H}}(\lambda )=\frac{\Vert u{\Vert }_{\mathcal{H}}\Vert w{\Vert }_{\mathcal{H}}}{(w,u{)}_{\mathcal{H}}}.$$

$$\Vert u{\Vert }_{\mathcal{H}}\Vert w{\Vert }_{\mathcal{H}}=\Vert x{\Vert }_{{\mathbb{C}}^m}\left(\Vert y{\Vert }_{{\mathbb{C}}^m}+\Vert v{\Vert }_{\mathcal{H}}\right)\ge \Vert x{\Vert }_{{\mathbb{C}}^m}\Vert y{\Vert }_{{\mathbb{C}}^m}.$$

$${\kappa }_{{\mathbb{C}}^m}(\lambda )=\frac{\Vert x{\Vert }_{{\mathbb{C}}^m}\Vert y{\Vert }_{{\mathbb{C}}^m}}{(y,x{)}_{{\mathbb{C}}^m}}\le \frac{\Vert u{\Vert }_{\mathcal{H}}\Vert w{\Vert }_{\mathcal{H}}}{(w,u{)}_{\mathcal{H}}}={\kappa }_{\mathcal{H}}(\lambda ),$$

$${\lambda }_{ϵ}(\mathcal{L})=\{z\in \mathbb{C}:\Vert (z\mathcal{I}-\mathcal{L}{)}^{-1}{\Vert }_{\mathcal{H}}>1/ϵ\}.$$

$$\bigcup _{\Vert \mathcal{E}{\Vert }_{\mathcal{H}}<ϵ}\lambda (\mathcal{L}+\mathcal{E})={\lambda }_{ϵ}(\mathcal{L}).$$

$$\Vert (zI-\mathcal{L}{)}^{-1}{\Vert }_{\mathcal{H}}=\underset{f\in \mathcal{H},\Vert f{\Vert }_{\mathcal{H}}=1}{\sup }\Vert (z\mathcal{I}-\mathcal{L}{)}^{-1}f{\Vert }_{\mathcal{H}}\ge \underset{x\in {\mathbb{C}}^m,\Vert x{\Vert }_{{\mathbb{C}}^m}=1}{\max }\Vert (z\mathcal{I}-\mathcal{L}{)}^{-1}Qx{\Vert }_{\mathcal{H}}.$$

$$\Vert (zI-\mathcal{L}{)}^{-1}{\Vert }_{\mathcal{H}}\ge \Vert (Q^{\ast }(z\mathcal{I}-\mathcal{L})Q{)}^{-1}{\Vert }_{{\mathbb{C}}^m}=\Vert (zI-Q^{\ast }\mathcal{L}Q{)}^{-1}{\Vert }_{{\mathbb{C}}^m}.$$

$$\hat{\lambda }\in {\lambda }_{ϵ}(Q^{\ast }\mathcal{L}Q)\subset {\lambda }_{ϵ}(\mathcal{L}).$$

$${\hat{\mathcal{P}}}_{\mathcal{V}}=\frac{1}{2\pi i}\sum _{k=1}^{\ell }{w}_k(z_k\mathcal{I}-\mathcal{L}{)}^{-1}.$$

$$(z_k\mathcal{I}-\mathcal{L})g_{i,k}=f_i,g_{i,k}(\pm 1)=\cdots =g_{i,k}^{(N/2)}(\pm 1)=0,1\le i\le m.$$

$${\hat{L}}_{ij}=({\hat{q}}_i,\mathcal{L}{\hat{q}}_j{)}_{\mathcal{H}},1\le i,j\le m,$$

$$(f,g{)}_{L^2([-1,1])}={\int }_{-1}^1\overline{f(x)}g(x)dx.$$

$$(f,g{)}_{L^2([-1,1])}\approx \sum _{k=1}^p{w}_k\overline{f(x_k)}g(x_k).$$

$$-\frac{d^{2}u}{d{x}^2}+x^{2}u=\lambda \cosh (x)u,u(\pm 1)=0.$$

$$(v,u{)}_w={\int }_{-1}^1\overline{v}u\cosh (x)dx.$$

$$\sqrt{{\lambda }_n}\sim \frac{n\pi }{{\int }_{-1}^1\sqrt{\cosh (x)}dx},n\to \infty .$$

$$-\frac{d^{2}u}{d{x}^2}=\lambda x^{3}u,u(\pm 1)=0,$$

$$\sqrt{{\lambda }_n}\sim \frac{(n-1/4)\pi }{{\int }_0^1{x}^{3/2}dx},{\lambda }_n>0,n\to \infty .$$

$$(v,u{)}_{|w|}={\int }_{-1}^1\overline{v}u|x{|}^{3}dx.$$

$${\mathcal{P}}_{\mathcal{V}}=\frac{1}{2\pi i}{\int }_{\partial \Omega }(z{\mathcal{L}}_2-{\mathcal{L}}_1{)}^{-1}{\mathcal{L}}_{2}dz.$$

$$s(z)=\sum _{k=1}^{\ell }\frac{w_k}{z_k-z},z\in \mathbb{C}\setminus \{z_1,\dots ,z_l\}.$$

$$\chi (z)=\frac{1}{2\pi i}{\int }_{\partial \Omega }\frac{dw}{w-z},z\in \mathbb{C}\setminus \partial \Omega .$$

$$\rho (\mathcal{B})=\max \{|z|:z\in \lambda (\mathcal{B})\}.$$

$$\rho ({(\mathcal{I}-\mathcal{P}}_{\mathcal{V}})\mathcal{B})<|{\lambda }_m|.$$

$$\Vert {\hat{u}}_k-u{\Vert }_{\mathcal{H}}\le {\left(\left|\rho /\lambda \right|+{ϵ}_k\right)}^k\Vert (I-{\mathcal{P}}_{\mathcal{V}})Fx{\Vert }_{\mathcal{H}},k=1,2,3,\dots ,$$

Full text

{AMS}

34L16, 65F15

1 Introduction

In this paper, we consider differential eigenvalue problems posed on the interval $[-1,1]$, i.e.,

[TABLE]

Here, $\mathcal{L}$ is a linear, ordinary differential operator of even order $N$. A complex number $\lambda$ and a function $u$ satisfying Eq. 1 are called an eigenvalue and eigenfunction of $\mathcal{L}$, respectively. We focus on computing the eigenvalues of $\mathcal{L}$ contained in a simply connected region $\Omega\subset\mathbb{C}$. Throughout, we assume that the boundary $\partial\Omega$ is a rectifiable, simple closed curve and that the spectrum $\lambda(\mathcal{L})$ of $\mathcal{L}$ is discrete, does not intersect $\partial\Omega$, and only finitely many eigenvalues counting multiplicities are in $\Omega$. To simplify discussion about the eigenfunctions of Eq. 1, we assume that there are eigenfunctions of $\mathcal{L}$ that form a basis for the invariant subspace of $\mathcal{L}$ associated with $\Omega$.

Since the development of the QR algorithm in the 1960s, the standard methods for solving Eq. 1 have adopted a “discretize-then-solve” paradigm. These algorithms first discretize $\mathcal{L}$ to obtain a finite matrix eigenvalue problem and then solve the matrix eigenvalue problem with algorithms from numerical linear algebra [18, 22, 26, 40]. Motivated by mathematical software for highly adaptive computations with functions [19], we propose an alternative strategy: an algorithm that solves Eq. 1 by directly manipulating $\mathcal{L}$ at the continuous level and only discretizes functions, not operators. By designing an eigensolver for $\mathcal{L}$ rather than intermediate discretizations, we are able to leverage spectrally accurate approximation schemes for functions while avoiding several pitfalls that plague spectral discretizations of Eq. 1 (see [61], [23, Ch. 2], and [60, Ch. 30]). For this reason, we view our proposed algorithms as adopting a “solve-then-discretize” paradigm. This paradigm has been applied to Krylov methods [24], iterative eigensolvers [29], and contour integral projection eigensolvers [8] for differential operators. Related techniques for computing with operators on infinite dimensional spaces have been proposed and studied in [15, 38].

As an example of the advantages of our methodology, consider the simplest possible differential eigenvalue problem given by

[TABLE]

The eigenvalues of Eq. 2 are $\lambda_{k}=(k\pi/2)^{2}$, for $k\geq 1$, and are well-conditioned due to the fact that the eigenfunctions form a complete orthonormal set in the Hilbert space $L^{2}([-1,1])$ [32, p. 382]. However, spectral discretizations of Eq. 2 lead to highly non-normal matrices with eigenvalues that are far more ill-conditioned than expected. Due to this ill-conditioning, the accuracy in the computed eigenvalues can be extremely variable and difficult to predict, ranging from a few digits to nearly full precision (see Fig. 1).

It is possible to use structure-preserving spectral discretizations to solve Eq. 2 accurately [14, 50]. However, there is a lack of literature on designing spectral discretizations of Eq. 1 when $\mathcal{L}$ is self-adjoint or normal with respect to an inner product other than $L^{2}([-1,1])$. On the other hand, our solve-then-discretize methodology automatically preserves the normality or self-adjointness of $\mathcal{L}$ with respect to a relevant Hilbert space $\mathcal{H}$, provided that the inner product $(\cdot,\cdot)_{\mathcal{H}}$ can be evaluated.

At the heart of our approach is an operator analogue of the FEAST matrix eigensolver, which we briefly outline:

• (1)

We construct a basis for the eigenspace $\mathcal{V}$ corresponding to $\Omega$ by sampling the range of the associated spectral projector $\mathcal{P}_{\mathcal{V}}$.

• (2)

We extract an $\mathcal{H}$-orthonormal basis for $\mathcal{V}$ with a continuous analogue of the QR factorization [58].

• (3)

We perform a Rayleigh–Ritz projection [47, p. 98] of $\mathcal{L}$ onto $\mathcal{V}$ with the orthonormal basis in (2). We solve the resulting matrix eigenvalue problem to obtain approximations to the eigenvalues of $\mathcal{L}$ in $\Omega$.

As with the FEAST matrix eigensolver, the spectral projector $\mathcal{P}_{\mathcal{V}}$ is applied approximately via a quadrature rule approximation. For matrices, this involves solving shifted linear systems, while for differential operators one needs to solve shifted linear differential equations. We solve these differential equations with the ultraspherical spectral method, which is a well-conditioned spectral method that is capable of resolving solutions that exhibit layers, rapid oscillations, and weak corner singularities [37].

Critically, we discretize basis functions for $\mathcal{V}$ as opposed to discretizing the differential operator $\mathcal{L}$ when solving Eq. 1. While discretizations of a normal operator $\mathcal{L}$ can lead to non-normal matrices, the Rayleigh–Ritz projection described in (3) always leads to a normal matrix eigenvalue problem when $\mathcal{L}$ is normal (see Theorems 3.1 and 3.3). In fact, we prove that using a sufficiently good approximate basis for $\mathcal{V}$ does not significantly increase the sensitivity of the eigenvalues when $\mathcal{L}$ is normal (see Section 5.2 for a precise statement). The result is a highly accurate eigensolver for normal differential operators $\mathcal{L}$, requiring $\mathcal{O}(mMN\log(N)+m^{2}N+m^{3})$ floating point operations, where $m={\rm dim}(\mathcal{V})$ and $M$ and $N$ are the polynomial degrees used to resolve the variable coefficients in $\mathcal{L}$ and the eigenfunctions in $\mathcal{V}$, respectively.

The eigensolver we develop is competitive in the high-frequency regime because it efficiently resolves oscillatory basis functions in $\mathcal{V}$. Furthermore, it handles operators that are self-adjoint or normal with respect to non-standard Hilbert spaces. Finally, our algorithm is parallelizable like the FEAST matrix eigensolver [44]. This work is a step towards closing the gap between the frequency regimes that are accessible to computational techniques and asymptotic methods for differential eigenvalue problems posed on higher dimensional domains [10, 11].

The paper is organized as follows. We begin in Section 2 by reviewing FEAST for matrix eigenvalue problems. In Section 3 we introduce an analogue of FEAST for differential operators and show that the operator analogue preserves eigenvalue sensitivity. In Section 4 we discuss a practical implementation of the operator analogue and provide two examples from Sturm–Liouville theory to illustrate its capabilities in the high-frequency regime. We analyze the convergence and stability of this implementation in Section 5. Sections 6 and 7 develop further applications of the solve-then-discretize paradigm, including an operator analogue of the Rayleigh Quotient iteration and an extension of FEAST to unbounded search regions.

2 The FEAST matrix eigensolver

The FEAST matrix eigensolver uses approximate spectral projection to compute the eigenvalues of a matrix $A\in\mathbb{C}^{n\times n}$ in a region of interest $\Omega\subset\mathbb{C}$ [33] (see Fig. 2). It is usually more computationally efficient than standard eigensolvers when the number of eigenvalues in $\Omega$ is much smaller than $n$. The dominating computational cost of FEAST is solving several independent shifted linear systems, but these can be performed in parallel [33].

There are three essential ingredients to FEAST:

• (i)

Spectral projector. Let $\lambda_{1},\dots,\lambda_{m}$ be the eigenvalues of $A$ in $\Omega$ and let $\mathcal{V}$ be the associated invariant subspace of $A$, i.e., $A\mathcal{V}=\mathcal{V}$. The spectral projector onto $\mathcal{V}$ is defined as

[TABLE]

The important fact here is that ${\rm range}(P_{\mathcal{V}})=\mathcal{V}$ and so $P_{\mathcal{V}}$ is a projection onto the invariant subspace of $A$ [32].

• (ii)

Basis for $\mathcal{V}$. FEAST uses the spectral projector to construct a basis for $\mathcal{V}$. It begins with a matrix $Y\in\mathbb{C}^{n\times m}$ with linearly independent columns that are not in ${\rm ker}(P_{\mathcal{V}})$, then it computes $Z=P_{\mathcal{V}}Y$. The columns of $Z$ span $\mathcal{V}$ and a QR factorization of $Z$ provides an orthonormal basis, $Q$, for $\mathcal{V}$.

• (iii)

Rayleigh–Ritz projection. Having obtained an orthonormal basis for $\mathcal{V}$, FEAST solves $A_{Q}x=\lambda x$ using a dense eigensolver [44], where $A_{Q}=Q^{*}AQ$. Since ${\rm range}(Q)=\mathcal{V}$, the eigenvalues of $A_{Q}$ are the eigenvalues of $A$ that lie inside $\Omega$. When $A_{Q}$ is diagonalizable, the eigenvectors of $A$ are given by $u_{i}=Qx_{i}$, for $i=1,\dots,m$, where $x_{1},\dots,x_{m}$ are the corresponding eigenvectors of $A_{Q}$.

For practical computation, FEAST approximates the contour integral in Eq. 3 with a quadrature rule. Given a quadrature rule with nodes $z_{1},\ldots,z_{\ell}$ and weights $w_{1},\ldots,w_{\ell}$, one can approximate $P_{\mathcal{V}}Y$ by

[TABLE]

In this case, the eigenpairs of $A_{Q}$ provide approximations to the eigenpairs of $A$, known as Ritz values and vectors [54]. To refine the accuracy of the Ritz values and vectors, a more accurate quadrature rule can be used to compute $P_{\mathcal{V}}Y$ [54]. FEAST also refines the approximate eigenvalues and eigenvectors by applying $P_{\mathcal{V}}$ to the $n\times m$ block of approximate eigenvectors using a quadrature rule and iterating (ii) and (iii) until convergence. To fully understand this refinement process, one must examine FEAST through the lens of rational subspace iteration [54].

When the dimension $m$ of the invariant subspace $\mathcal{V}$ is unknown, there are several techniques for estimating $m$ and selecting an appropriate value [33, 36, 54]. Most of these techniques can be incorporated into the operator analogue of FEAST. Consequently, we assume that $m$ is known throughout the paper and focus on the algorithmic and theoretical aspects of FEAST that are relevant in the operator setting.

Curiously, the originally proposed FEAST algorithm does not compute an orthonormal basis for $\mathcal{V}$ before performing the Rayleigh–Ritz projection [44, 33].111The FEAST algorithm for non-Hermitian matrices utilizes dual bases for the left and right eigenspaces to improve stability [33]. However, when $Q$ has orthonormal columns and ${\rm range}(Q)$ is an invariant subspace of $A$, then the eigenvalues of the small matrix $Q^{*}AQ$ are no more sensitive to perturbations than the original eigenvalues of $A$. This highly desirable property follows from an examination of the structure of the left and right invariant subspaces of $Q^{*}AQ$ or, alternatively, from the $\epsilon$-pseudospectra of $Q^{*}AQ$ [60, p. 382].

3 An operator analogue of FEAST

The FEAST matrix algorithm provides a natural starting point for an operator analogue because it provides a recipe to construct a small matrix $Q^{*}AQ$ whose eigenvalues coincide with those of $A$ inside $\Omega$ and have related invariant subspaces. Moreover, the eigenstructure of $Q^{*}AQ$ reflects the eigenstructure of $A$ when the columns of $Q$ are orthonormal. As the sensitivity of the eigenvalues of $A$ depends intimately on the structure of the associated eigenvectors, $Q^{*}AQ$ may be used to compute the desired eigenvalues of $A$ efficiently without sacrificing accuracy. Here, we generalize FEAST so that it constructs a matrix whose eigenvalues coincide with those of a differential operator inside $\Omega$.

3.1 FEAST for differential operators

In place of a matrix $A$ acting on vectors from $\mathbb{C}^{n}$, we now consider a differential operator $\mathcal{L}$ acting on functions from a Hilbert space $\mathcal{H}$. As described in Section 2, the FEAST recipe prescribes a spectral projection to compute a basis for $\mathcal{V}$, which is then used for the Rayleigh–Ritz projection to construct a matrix representation on $\mathcal{V}$. Throughout, we require that $\mathcal{L}$ be a closed operator222An operator $\mathcal{A}:\mathcal{D}(\mathcal{A})\rightarrow\mathcal{H}$ is closed if its graph is a closed linear subspace of $\mathcal{H}\times\mathcal{H}$ [32, p. 165]. and that its domain $\mathcal{D}(\mathcal{L})$ is dense in $\mathcal{H}$.

• (i)

Spectral projector. Although $\mathcal{L}$ is unbounded, the resolvent $(z\mathcal{I-L})^{-1}$ is bounded when $z\not\in\lambda(\mathcal{L})$ and the spectral projector onto $\mathcal{V}$ may be defined via contour integral [32, p. 178]. It is given by

[TABLE]

• (ii)

Basis for $\mathcal{V}$. With the spectral projector at our disposal, we apply $\mathcal{P}_{\mathcal{V}}$ to functions $f_{1},\dots,f_{m}$ in $\mathcal{H}\setminus{\rm ker}(\mathcal{P}_{\mathcal{V}})$ to obtain a basis of functions $v_{1},\dots,v_{m}$ for $\mathcal{V}$. Orthonormalizing $v_{1},\dots,v_{m}$ with respect to the inner product $(\cdot,\cdot)_{\mathcal{H}}$ on $\mathcal{H}$ gives us an $\mathcal{H}$-orthonormal basis $q_{1},\dots,q_{m}$ for $\mathcal{V}$.

• (iii)

Rayleigh–Ritz projection. To compute a matrix representation $L$ of $\mathcal{L}$ on $\mathcal{V}$, the Rayleigh–Ritz projection is performed using the inner product on $\mathcal{H}$. The elements of $L$ are given by $L_{ij}=(q_{i},\mathcal{L}q_{j})_{\mathcal{H}}$ for $1\leq i,j\leq m$. The eigenvalues of $L$ are precisely the eigenvalues $\lambda_{1},\dots,\lambda_{m}$ of $\mathcal{L}$ that lie inside $\Omega$. The eigenfunctions of $\mathcal{L}$ are recovered from the eigenvectors $x_{1},\dots,x_{m}$ of $L$ by computing $u_{i}=\sum_{k=1}^{m}x_{i}^{(k)}q_{k}$, for $i=1,\dots,m$, where $x_{i}^{(k)}$ is the $k^{th}$ component of $x_{i}$.

To avoid a clutter of indices, we employ the notation of quasimatrices.333A quasimatrix is a matrix whose columns (or rows) are functions defined on an interval $[a,b]$, in contrast to matrices whose columns (or rows) are vectors [19, Ch.6]. If $Q$ is the quasimatrix with columns $q_{1},\dots,q_{m}$, then the matrix $L$ whose elements are $L_{ij}=(q_{i},\mathcal{L}q_{j})_{\mathcal{H}}$ in (iii) is expressed compactly in quasimatrix notation as $L=Q^{*}\mathcal{L}Q$. Here, $Q^{*}$ is the conjugate transpose of the quasimatrix $Q$ so its rows are complex conjugates of the functions $q_{1},\dots,q_{m}$.

The analogue of FEAST for differential operators is summarized in Algorithm 2 using quasimatrix notation so that it resembles its matrix counterpart.

Keep in mind that Algorithm 2 is a formal algorithm. In general, we cannot apply the spectral projector exactly, nor represent the basis $\mathcal{V}$ exactly with finite memory. A practical implementation is discussed in Section 4.

3.2 Condition number of the Ritz values

As illustrated in Fig. 1, the eigenvalues of matrix discretizations of $\mathcal{L}$ can be more sensitive to perturbations than the eigenvalues of $\mathcal{L}$. The advantage of our FEAST approach in Section 3.1 is that the Ritz values, i.e., the eigenvalues of $Q^{*}\mathcal{L}Q$, are no more sensitive to perturbations than the original eigenvalues of $\mathcal{L}$ when ${\rm range}(Q)$ is an invariant subspace of $\mathcal{L}$.

To see this, let $\lambda$ be a simple eigenvalue of a differential operator $\mathcal{L}$. Let $u,w\in\mathcal{H}$ satisfy $\mathcal{L}u=\lambda u$ and $\mathcal{L}^{*}w=\overline{\lambda}w$, where $\overline{\lambda}$ denotes the complex conjugate of $\lambda$. The condition number444Although this formula is usually associated with the condition number for a simple eigenvalue of a matrix, its proof extends to our general setting [55, Theorem 5]. of $\lambda$ is given by [6, Theorem 2.3]

[TABLE]

The condition number $\kappa_{\mathcal{H}}(\lambda)$ quantifies the worst-case first-order sensitivity of $\lambda$ to perturbations of $\mathcal{L}$. For instance, if we compute $\lambda$ using a backward stable algorithm in floating point arithmetic, we expect to achieve an accuracy of at least $\kappa_{\mathcal{H}}(\lambda)\epsilon_{\text{mach}}$, where $\epsilon_{\text{mach}}$ is machine precision [59, Theorem 15.1].

Theorem 3.1**.**

Let $\mathcal{L}:\mathcal{D}(\mathcal{L})\rightarrow\mathcal{H}$ be a closed and densely defined operator on a Hilbert space $\mathcal{H}$, $Q:\mathbb{C}^{m}\rightarrow\mathcal{H}$ be an invariant subspace of $\mathcal{L}$ satisfying $Q^{*}Q=I$, and $L=Q^{*}\mathcal{L}Q$. Suppose that $u\in{\rm range}(Q)$ satisfies $\mathcal{L}u=\lambda u$ and $w$ satisfies $\mathcal{L}^{*}w=\overline{\lambda}w$, where $\mathcal{L}^{*}$ denotes the adjoint of $\mathcal{L}$ and $\lambda$ is a simple eigenvalue of $\mathcal{L}$ with condition number $\kappa_{\mathcal{H}}(\lambda)$. Then,

$LQ^{*}u=\lambda Q^{*}u$* and $L^{*}Q^{*}w=\overline{\lambda}Q^{*}w$,*

• 2)

$(Q^{*}w,Q^{*}u)_{\mathbb{C}^{m}}=(w,u)_{\mathcal{H}}$, and

• 3)

$\kappa_{\mathbb{C}^{m}}(\lambda)\leq\kappa_{\mathcal{H}}(\lambda)$.

Proof 3.2**.**

Denote $x=Q^{*}u$ and $y=Q^{*}w$. We prove the statements of the theorem in order. 1) Since $u\in{\rm range}(Q)$, we can write $u=Qx$. Then, $\mathcal{L}(Qx)=\lambda(Qx)$ implies that $Q^{*}\mathcal{L}Qx=\lambda x$ using the fact that $Q^{*}Q=I$. For the left eigenvector, we write $w=Qy+v$ for some $v\in{\rm range}(Q)^{\perp}$. Rewriting the adjoint equation for $w$, we find that $\mathcal{L}^{*}(Qy+v)=\overline{\lambda}(Qy+v)$ and multiplying by $Q^{*}$ on both sides yields $Q^{*}\mathcal{L}^{*}Qy=\overline{\lambda}y$. Here, we have used the fact that $Q^{*}\mathcal{L}^{*}v=0$, which holds because $v^{*}\mathcal{L}Q=0$. 2) By calculating $(w,u)_{\mathcal{H}}=(Qy+v,Qx)_{\mathcal{H}}$, we find that $(w,u)_{\mathcal{H}}=(Qy,Qx)_{\mathcal{H}}$ because $v\in{\rm range}(Q)^{\perp}$. Moreover, since $Q^{*}Q=I$ we conclude that $(w,u)_{\mathcal{H}}=(y,Q^{*}Qx)_{\mathbb{C}^{m}}=(y,x)_{\mathbb{C}^{m}}$. 3) We know that $\lVert u\rVert_{\mathcal{H}}=(Qx,Qx)_{\mathcal{H}}=(x,x)_{\mathbb{C}^{m}}=\lVert x\rVert_{\mathbb{C}^{m}}$ and $\lVert w\rVert_{\mathcal{H}}=(Qy+v,Qy+v)_{\mathcal{H}}=\lVert y\rVert_{\mathbb{C}^{m}}+\lVert v\rVert_{\mathcal{H}}$. Therefore,

[TABLE]

Referring to 2) for equality of the inner products in the denominator, we have

[TABLE]

*which concludes the proof. *

Theorem 3.1 shows that if $\mathcal{L}$ is a normal operator, then $u=w$ and we have $\kappa_{\mathbb{C}^{m}}(\lambda)=\kappa_{\mathcal{H}}(\lambda)=1$. For non-normal operators, item 3) of Theorem 3.1 may seem to erroneously indicate that ill-conditioning in the eigenvalues of $\mathcal{L}$ can be overcome by a Rayleigh–Ritz projection. However, when $\mathcal{L}$ is non-normal the spectral projector $\mathcal{P}_{\mathcal{V}}$ is an oblique projection and computing the basis $Q$ may be itself an ill-conditioned problem. Theorem 3.1 also illustrates why the operator analogue of FEAST leads to a well-conditioned matrix eigenvalue problem when the differential eigenvalue problem is well-conditioned. By computing an $\mathcal{H}$-orthonormal basis for the Rayleigh–Ritz projection, the relevant structure in the eigenspaces of $\mathcal{L}$ and $\mathcal{L}^{*}$ is preserved. However, the first-order analysis above is limited to simple eigenvalues.

3.3 Pseudospectra of $\mathbf{Q^{*}\mathcal{L}Q}$

To go beyond first-order sensitivity analysis, we compare the $\epsilon$-pseudospectra of $\mathcal{L}$ and $Q^{*}\mathcal{L}Q$. Fix any $\epsilon>0$ and let $\mathcal{L}:\mathcal{D}(\mathcal{L})\rightarrow\mathcal{H}$ be a closed operator with a domain $\mathcal{D}(\mathcal{L})$ that is dense in $\mathcal{H}$. The $\epsilon$-pseudospectrum of $\mathcal{L}$ is defined as the set [60, p. 31]

[TABLE]

Here, we adopt the usual convention that $\lVert(z\mathcal{I-L})^{-1}\rVert_{\mathcal{H}}=\infty$ when $z\in\lambda(\mathcal{L})$ so that $\lambda(\mathcal{L})\subset\lambda_{\epsilon}(\mathcal{L})$. The $\epsilon$-pseudospectrum set of $\mathcal{L}$ bounds the region in which the eigenvalues of the perturbed operator $\mathcal{L+E}$ with $\lVert\mathcal{E}\rVert_{\mathcal{H}}<\epsilon$ can be found [60, p. 31]. This means that $\lambda(\mathcal{L+E})\subset\lambda_{\epsilon}(\mathcal{L})$. In fact, there is an equivalence so that [60, p. 31]

[TABLE]

This allows us to relate the sensitivity of the eigenvalues of $\mathcal{L}$ and $Q^{*}\mathcal{L}Q$ by comparing the resolvent norms $\lVert(z\mathcal{I-L})^{-1}\rVert_{\mathcal{H}}$ and $\lVert(zI-Q^{*}\mathcal{L}Q)^{-1}\rVert_{\mathbb{C}^{m}}$, respectively.

A useful generalization of Theorem 3.1 is that the $\epsilon$-pseudospectrum of $Q^{*}\mathcal{L}Q$ is contained in the $\epsilon$-pseudospectrum of $\mathcal{L}$. Since this holds for any $\epsilon>0$, it demonstrates that the eigenvalues (even those with multiplicity) of $Q^{*}\mathcal{L}Q$ are no more sensitive to perturbations than those of $\mathcal{L}$. This inclusion result is well-known in the matrix case where projection methods are a popular method for approximating the $\epsilon$-pseudospectra of large data-sparse matrices [60, p. 381].

Theorem 3.3**.**

*Let $\mathcal{L}:\mathcal{D}(\mathcal{L})\rightarrow\mathcal{H}$ be a closed and densely defined operator on a Hilbert space $\mathcal{H}$. For a fixed $\epsilon>0$, suppose that $Q:\mathbb{C}^{m}\rightarrow\mathcal{H}$ satisfies $Q^{*}Q=I$ and that ${\rm range}(Q)$ is an invariant subspace of $\mathcal{L}$. Then, $\lambda_{\epsilon}(Q^{*}\mathcal{L}Q)\subset\lambda_{\epsilon}(\mathcal{L})$. *

Proof 3.4**.**

We follow the proof of Proposition 40.1 in [60, p. 382] for matrices, but with a closed operator. Since $Qx\in\mathcal{H}$ for any $x\in\mathbb{C}^{m}$, we have that

[TABLE]

Now, when ${\rm range}(Q)$ is an invariant subspace of $\mathcal{L}$ and $Q^{*}Q=I$, we have that $\lVert(z\mathcal{I-L})^{-1}Qx\rVert_{\mathcal{H}}=\lVert Q^{*}(z\mathcal{I-L})^{-1}Qx\rVert_{\mathbb{C}^{m}}$. Because $QQ^{*}f=f$ for all $f\in\mathcal{V}$, we can check that $Q^{*}(z\mathcal{I-L})^{-1}Q=(Q^{*}(z\mathcal{I-L})Q)^{-1}.$ Since $Q^{*}Q=I$, it follows that

[TABLE]

*Therefore, $z\in\lambda_{\epsilon}(\mathcal{L})$ whenever $z\in\lambda_{\epsilon}(Q^{*}\mathcal{L}Q)$. *

The inclusion in Theorem 3.3 may be strict, indicating that the eigenvalues of $Q^{*}\mathcal{L}Q$ are less sensitive than those of $\mathcal{L}$. For example, this may occur when the projection onto ${\rm range}(Q)$ targets a subset of well-conditioned eigenvalues of $\mathcal{L}$. However, we emphasize that ill-conditioning in the eigenvalues of $\mathcal{L}$ cannot be overcome by a Rayleigh–Ritz projection: in general, the situation is complicated [60, Ch. 40].

Theorem 3.3 is useful for studying the stability of Algorithm 2. If an approximate eigenvalue $\hat{\lambda}$ of $Q^{*}\mathcal{L}Q$ is computed with an error tolerance of $\epsilon>0$, then

[TABLE]

From this, we know by Eq. 8 that $\hat{\lambda}$ is an eigenvalue of a perturbed operator $\mathcal{L+E}$ with $\lVert\mathcal{E}\rVert_{\mathcal{H}}<\epsilon$. In other words, the operator analogue of FEAST, Algorithm 2, is backward stable. As we see in Section 5, Theorem 3.3 is also the starting point for a stability analysis when the spectral projection is no longer exact and the Rayleigh–Ritz projection is performed with a matrix $\hat{Q}$ that only approximates a basis for an invariant subspace of $\mathcal{L}$.

4 A practical differential eigensolver based on an operator analogue of FEAST

The operator analogue of FEAST requires the manipulation of objects such as differential operators, functions, and contour integrals (see Algorithm 2). For a practical implementation, these objects must be discretized; however, we avoid discretizing $\mathcal{L}$ directly. Instead, we construct polynomial approximants to the basis for $\mathcal{V}$ by approximately solving shifted linear ODEs. These polynomial approximants are used in the Rayleigh–Ritz projection to compute the eigenvalues of $\mathcal{L}$ in $\Omega$.

Let $z_{1},\dots,z_{\ell}$ and $w_{1},\dots,w_{\ell}$ be a set of quadrature nodes and weights to approximate the integral in Eq. 5. As FEAST does in the matrix case, we approximate $\mathcal{P}_{\mathcal{V}}$ in Eq. 5 with a quadrature rule as follows:

[TABLE]

If $F$ is a quasimatrix with columns $f_{1},\dots,f_{m}\in\mathcal{H}$, then $\mathcal{P}_{\mathcal{V}}F$ is replaced by the approximation $\mathcal{\hat{P}}_{\mathcal{V}}F=\frac{1}{2\pi i}\sum_{k=1}^{\ell}w_{k}(z_{k}\mathcal{I-L})^{-1}F$. Therefore, to compute $\mathcal{\hat{P}}_{\mathcal{V}}F$ we need to solve $\ell$ shifted linear ODEs, each with $m$ righthand sides, i.e.,

[TABLE]

If the quasimatrix with columns $g_{1,k},\dots,g_{m,k}$ is denoted by $G_{k}$ for $k=1,\dots,\ell$, then we have $\mathcal{\hat{P}}_{\mathcal{V}}F=\sum_{k=1}^{\ell}w_{k}G_{k}$.

To construct a basis for $\mathcal{V}$, it is important to choose $F$ so that the columns of $\hat{V}=\mathcal{\hat{P}}_{\mathcal{V}}F$ are linearly independent and, if possible, well-conditioned. In analogy with the implementation of matrix FEAST [44, 33], we obtain the columns of $F$ by selecting $m$ band-limited random functions555A periodic band-limited random function on $[-L,L]$ is a periodic function defined by a truncated Fourier series with random (e.g. standard Gaussian distributed) coefficients. In the non-periodic setting, the Fourier series is defined on a larger interval $[-L^{\prime},L^{\prime}]$ and the domain is then truncated [21]. on $[-1,1]$ [21]. When $\mathcal{L}$ is a normal operator, this typically yields a well-conditioned basis $\hat{V}$.

We now outline the key implementation details of our differential eigensolver:

• (i)

Approximate spectral projection. To compute $\hat{V}=\mathcal{\hat{P}}_{\mathcal{V}}F$, we solve the shifted linear ODEs in Eq. 10 using the ultraspherical spectral method [37]. The ultraspherical spectral method leads to well-conditioned linear systems and is capable of accurately resolving the functions $g_{i,k}$ even when they are highly oscillatory or have boundary layers. Moreover, an adaptive QR factorization automatically determines the degree of the polynomial interpolants needed to approximate the functions $g_{i,k}$ to near machine precision [39, 38]. After accurately resolving the functions $g_{i,k}$, we can accurately compute a basis for $\mathcal{V}$ provided that both the spectral projector is well-conditioned (i.e., $\mathcal{L}$ is not highly non-normal) and the quadrature rule is sufficiently accurate.

• (ii)

Orthonormal basis. To compute an orthonormal basis $\hat{Q}$ for the columns of $\hat{V}$, we compute a QR factorization of the quasimatrix $\hat{V}$ by Householder triangularization [58]. The Householder reflectors are constructed with respect to the inner product $(\cdot,\cdot)_{\mathcal{H}}$ so that the columns of $\hat{Q}$ are $\mathcal{H}$-orthonormal.

• (iii)

Computing $\mathbf{\hat{Q}^{*}\mathcal{L}\hat{Q}}$. To construct the matrix $\hat{L}=\hat{Q}^{*}\mathcal{L}\hat{Q}$, we apply $\mathcal{L}$ to the columns of $\hat{Q}$ and then evaluate the action of $\hat{Q}^{*}$ on $\mathcal{L}\hat{Q}$. Multiplying $\hat{Q}^{*}$ with $\mathcal{L}\hat{Q}$ involves taking the inner products

[TABLE]

where $\hat{q}_{i}$ denotes the $i$th column of $\hat{Q}$. The eigenvalues $\hat{\lambda}_{1},\dots,\hat{\lambda}_{m}$ and eigenvectors $\hat{x}_{1},\dots,\hat{x}_{m}$ of the matrix $\hat{L}$ are computed using the QR algorithm [25, p. 385].

Critically, the inner product $(\cdot,\cdot)_{\mathcal{H}}$ used in the QR factorization of $\hat{V}$ and the construction of $\hat{Q}^{*}\mathcal{L}\hat{Q}$ depends on the choice of the Hilbert space $\mathcal{H}$. As long as we are able to evaluate $(\cdot,\cdot)_{\mathcal{H}}$, we can exploit the fact that $\mathcal{L}$ is self-adjoint or a normal operator with respect to $(\cdot,\cdot)_{\mathcal{H}}$ so that we can accurately compute the eigenvalues of $\mathcal{L}$ in $\Omega$ (see Theorem 5.2). For this reason, our algorithm is able to accurately compute the eigenvalues and eigenfunctions of differential operators that are self-adjoint with respect to non-standard Hilbert spaces (see Section 4.1).

Evaluating the inner product $(\cdot,\cdot)_{\mathcal{H}}$ usually means computing an integral, which we approximate with a quadrature rule. For example, if $\mathcal{H}=L^{2}([-1,1])$,

[TABLE]

Given the Gauss–Legendre quadrature nodes $x_{1},\dots,x_{p}$ and weights $w_{1},\dots,w_{p}$ on $[-1,1]$, then one uses the approximation [12]

[TABLE]

A practical implementation of the operator analogue of FEAST is presented in Algorithm 3. As with matrix FEAST, there are two approaches for improving the accuracy of the Ritz values $\hat{\lambda}_{1},\dots,\hat{\lambda}_{m}$ and vectors $\hat{Q}\hat{x}_{1},\dots,\hat{Q}\hat{x}_{m}$. The first is to improve the accuracy of the quadrature rule in Eq. 9. The second is to iterate the algorithm by replacing $F$ by the quasimatrix $\hat{U}$ with columns $\hat{u}_{i}=\hat{Q}\hat{x}_{i}$ for $1\leq i\leq m$, repeating the process if necessary.666When $\mathcal{L}$ is non-normal the Ritz vectors $\hat{Q}\hat{x}_{1},\dots,\hat{Q}\hat{x}_{m}$ may become numerically linearly dependent, which can lead to an ill-conditioned basis $\hat{V}$ in susbequent iterations. The robustness of Algorithm 3 may be improved by computing the Schur vectors $v_{1},\dots,v_{m}$ of $\hat{L}$ and using the orthonormal basis $\hat{Q}\hat{v}_{1},\dots,\hat{Q}\hat{v}_{m}$ to seed the next iteration [53]. For normal operators, this iteration generates a sequence of quasimatrices $\hat{Q}_{k}$ with $\mathcal{H}$-orthonormal columns that converge to an $\mathcal{H}$-orthonormal basis for the invariant subspace $\mathcal{V}$ as $k\rightarrow\infty$. This can be viewed as a rational subspace iteration and geometric convergence of the Ritz pairs is typical (see Section 5).

With either refinement strategy, the accuracy of the Ritz pairs may be monitored using the residual norm (see step 6 of Algorithm 3) as a proxy, just as in the matrix case. For normal operators, the error in the eigenvalues and eigenvectors computed by Algorithm 3 is typically $\mathcal{O}(\epsilon)$, where $\epsilon$ is the threshold for the residual norm in step 6. We defer a discussion of the convergence and stability of Algorithm 3 to Section 5. Additional resources on residual norm bounds for eigenvalues and eigenvectors of matrices and extensions to closed linear operators are found in [9, 13, 52].

In practice, when $\mathcal{L}$ is non-normal it may be beneficial to use a dual Rayleigh–Ritz projection $\hat{Q}_{L}^{*}\mathcal{L}\hat{Q}_{R}$, where the columns of $\hat{Q}_{R}$ approximate an orthonormal basis for the target eigenspace of $\mathcal{L}$ and the columns of $\hat{Q}_{L}$ approximate an orthonormal basis for the associated eigenspace of the adjoint $\mathcal{L}^{*}$. In the case of matrix FEAST, the use of the dual projection leads to a non-normal matrix eigensolver with improved robustness [33]. Although it is not difficult to adapt Algorithm 3 to an operator analogue of FEAST that uses dual projection, we focus on the implementation and analysis of the one-sided iteration.

Typically, solving the ODEs in (10) dominates the computational cost of Algorithm 3. With the ultraspherical spectral method, the computational complexity of solving the linear ODEs with $m$ distinct right hand sides is $\mathcal{O}(mMN\log(N))$ floating point operations (flops) [37]. Here, $N$ and $M$ are, respectively, the degrees of the truncated Chebyshev series needed to resolve the columns of $G_{k}$ and the variable coefficients in $\mathcal{L}$ to within the tolerance $\epsilon$ specified in Algorithm 3. In addition to the ODE solve, the QR factorization in (ii) requires $\mathcal{O}(m^{2}N)$ flops [58], while the dense eigenvalue computation with a small $m\times m$ matrix in (iii) takes $\mathcal{O}(m^{3})$ flops [25, p. 391]. The complexity of one iteration of Algorithm 3 is therefore $\mathcal{O}(mMN\log(N)+m^{2}N+m^{3})$ flops. In practice, convergence to machine precision usually occurs within two or three iterations.

4.1 Computing high-frequency eigenmodes

Algorithm 3 adaptively and accurately resolves basis functions for highly oscillatory eigenmodes and preserves the sensitivity of the eigenvalues of the differential operator $\mathcal{L}$, so it is well-suited to computing high-frequency eigenmodes when $\mathcal{L}$ is self-adjoint or normal with respect to $(\cdot,\cdot)_{\mathcal{H}}$. We provide two examples from Sturm–Liouville theory to illustrate the effectiveness of the solve-then-discretize methodology in the high-frequency regime.

4.1.1 A regular Sturm–Liouville eigenvalue problem

First consider a regular Sturm–Liouville eigenvalue problem (SLEP) given by

[TABLE]

This defines a self-adjoint differential operator with respect to the inner product

[TABLE]

Consequently, Eq. 12 possesses a complete $(\cdot,\cdot)_{w}$-orthonormal basis of eigenfunctions $u_{1},u_{2},u_{3},\ldots$ for the weighted Hilbert space $\mathcal{H}_{w}=\{u:\lVert u\rVert_{w}=\sqrt{(u,u)_{w}}<\infty\}$ and an unbounded set of real eigenvalues $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\cdots$.

Asymptotics for the large eigenvalues of Eq. 12 are given by [2]

[TABLE]

To accurately compute the large eigenvalues of Eq. 12 with Algorithm 3, we prescribe circular search regions with unit radius centered at the values given by the asymptotic formula in Eq. 14 (see Fig. 3). Each search region contains one eigenvalue.

4.1.2 An indefinite Sturm–Liouville eigenvalue problem

Next, we consider the following indefinite SLEP:

[TABLE]

which is closely related to models of light propagation in a nonhomogeneous material [2, 62]. Since the weight function $x^{3}$ changes sign at $x=0$, Eq. 15 has a bi-infinite sequence of eigenvalues [7]. We index them in order as $\cdots\leq\lambda_{-2}\leq\lambda_{-1}<0<\lambda_{1}\leq\lambda_{2}\leq\cdots$. The asymptotics for the positive eigenvalues are given by [3]

[TABLE]

A similar expansion holds for the negative eigenvalues [3].

In contrast to the previous example, the indefinite weight function $x^{3}$ means that Eq. 15 is not immediately associated with a self-adjoint operator on a Hilbert space. Instead, Eq. 15 is usually studied through the lens of a Krein space and the eigenfunctions form a Riesz basis for the Hilbert space with the inner product [16]

[TABLE]

We use the leading order asymptotics in Eq. 16 to identify search regions that are likely to contain an eigenvalue of Eq. 15. Because the ultraspherical spectral method used to solve the ODEs in step 2 of Algorithm 3 is efficient when applied to ODEs with smooth variable coefficients, it is convenient to treat Eq. 15 as a generalized eigenvalue problem, i.e., as $\mathcal{L}_{1}u=\lambda\mathcal{L}_{2}u$, where $\mathcal{L}_{1}u=-\frac{d^{2}u}{dx^{2}}$ and $\mathcal{L}_{2}u=x^{3}u$. The eigenvalues of the pencil $z\mathcal{L}_{2}-\mathcal{L}_{1}$ are then computed with a straightforward generalization of Algorithm 3 that is based on the spectral projector for the generalized eigenvalue problem, i.e.,

[TABLE]

The eigenvalues and eigenfunctions are automatically resolved to essentially machine precision because of the use of the adaptive QR solver (see Fig. 4).

5 Convergence and stability

The primary consequence of the approximations introduced in Algorithm 3 is that the spectral projector is no longer applied exactly. Therefore, the basis $\hat{Q}$ computed for the Rayleigh–Ritz projection is not an exact basis for the invariant subspace $\mathcal{V}$ of $\mathcal{L}$ and may require further refinement. Here, we view the iterative refinement procedure used in Algorithm 3 as a rational subspace iteration applied to a normal differential operator $\mathcal{L}$ in order to provide a preliminary analysis of the stability of the iteration and the sensitivity of the Ritz values. The main results may be summarized as follows.

• (i)

Algorithm 3 yields a sequence of quasimatrices $\hat{Q}_{1},\dots,\hat{Q}_{k}$ that (generically) converge geometrically to an orthonormal basis for the eigenspace $\mathcal{V}$ (see Theorem 5.1).

• (ii)

If $\hat{Q}_{k}$ is a sufficiently good approximation to an orthonormal basis for $\mathcal{V}$, then the $\epsilon$-psuedospectrum of $\hat{Q}_{k}^{*}\mathcal{L}\hat{Q}_{k}$ is contained in the $2\epsilon$-psuedospectrum of $\mathcal{L}$ itself (see Theorem 5.2).

• (iii)

Under mild conditions on the initial quasimatrix $F$ in Algorithm 3, the sequence $\lVert\mathcal{L}(\hat{Q}_{k}-Q)\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}$ is uniformly bounded as $k\rightarrow\infty$ (see Lemma 5.4).

Taken together, these results demonstrate that each iteration of Algorithm 3 yields uniformly consistent Ritz pairs that converge linearly to the desired eigenpair and that the unboundedness of $\mathcal{L}$ does not lead to instability. Note that this analysis does not take into account the impact of finite-precision arithmetic or the fact that the shifted differential equations in Eq. 10 are not solved exactly at each iteration (see the discussion at the end of Section 5.1). However, (ii) ensures that the eigenvalues of the small matrix $\hat{Q}_{k}^{*}\mathcal{L}\hat{Q}_{k}$ are not much more sensitive than the eigenvalues of $\mathcal{L}$. Therefore, provided that the eigenvalue problem for $\mathcal{L}$ is well-conditioned and we compute a sufficiently accurate approximation to a basis for $\mathcal{V}$, then we expect that the eigenvalues computed with Algorithm 3 provide an accurate approximation to the desired eigenvalues of $\mathcal{L}$.

5.1 Rational subspace iteration for differential operators

In analogy to the matrix case [54], Algorithm 3 may be interpreted as a filtered subspace iteration. Filtered subspace iteration is a variant of standard subspace iteration for computing a target subset of eigenvalues of a matrix $A$ [47, Ch. 5]. The main idea is to choose a filter function $s(\cdot)$ that is large on the targeted eigenvalues of $A$ and small on the unwanted eigenvalues of $A$. Applying the spectral transformation $s(A)$,777A spectral transformation $s(\cdot)$ may be applied to $A$ via the eigendecomposition of $A$, or more generally the Jordan decomposition. For example, if $A$ has eigendecomposition $A=X\Lambda X^{-1}$ with $\Lambda={\rm diag}(\lambda_{1},\dots,\lambda_{n})$, then $s(A)=X\,s(\Lambda)X^{-1}$, where $s(\Lambda)={\rm diag}(s(\lambda_{1}),\dots,s(\lambda_{n}))$. one uses standard subspace iteration to compute a basis for the eigenspace of $s(A)$ corresponding to its largest eigenvalues, i.e., the targeted eigenvalues of $A$. With an approximate basis for the eigenspace available, the eigenvalues and eigenvectors can be extracted with a Rayleigh–Ritz step.

From this perspective, Algorithm 3 computes the eigenvalues of $\mathcal{L}$ in $\Omega$ with the aid of a rational filter function induced by the quadrature rule in Eq. 9, i.e., $s(\cdot)$ is given by

[TABLE]

The functional calculus for unbounded normal operators ensures that if $\lambda_{i}$ is an eigenvalue of $\mathcal{L}$ with eigenfunction $u_{i}$, then $s(\lambda_{i})$ is an eigenvalue of $s(\mathcal{L})$ with eigenfunction $u_{i}$ [45, VIII.5].888The result [45, VIII.5] is stated for closed self-adjoint operators on $\mathcal{H}$, however, it extends immediately to closed normal operators on $\mathcal{H}$ if the spectral decomposition of a closed normal operator [46, Theorem 13.33] is used in place of the spectral decomposition of a self-adjoint operator. For information on the spectral decomposition of unbounded normal operators and the associated functional calculus, see [46, Ch. 13] and [20]. As the degree of the quadrature rule is increased, the rational function becomes an increasingly good approximation to the Cauchy integral

[TABLE]

Therefore, the eigenvalues of $\mathcal{L}$ in $\Omega$ are usually $\mathcal{O}(1)$ in size under the spectral transformation $s(\cdot)$ while the eigenvalues outside of $\Omega$ are much smaller.

We now turn to the convergence of the iteration described in Algorithm 3, which we interpret as a subspace iteration applied to the bounded linear operator $\mathcal{\hat{P}}_{\mathcal{V}}=s(\mathcal{L})$. It is helpful to introduce the notions of the spectral radius and a dominant eigenspace of a bounded linear operator $\mathcal{B}$. The spectral radius of a bounded linear operator $\mathcal{B}$ on a Hilbert space $\mathcal{H}$ is defined as [17, p. 99]

[TABLE]

The spectral radius is useful because it characterizes the asymptotic behavior of $\lVert\mathcal{B}^{k}\rVert_{\mathcal{H}}$, in the sense that $\rho(\mathcal{B})=\lim_{k\rightarrow\infty}\lVert\mathcal{B}^{k}\rVert_{\mathcal{H}}^{1/k}$ [17, Theorem 4.1.3]. Let $\mathcal{V}$ be an invariant subspace of $\mathcal{B}$ associated with eigenvalues $\lambda_{1}\geq\dots\geq\lambda_{m}$ and a spectral projector $\mathcal{P}_{\mathcal{V}}$. We say that $\mathcal{B}$ has dominant eigenspace $\mathcal{V}$ if

[TABLE]

The following theorem is an extension of a convergence analysis [47, p. 119] for matrix subspace iteration to the setting of bounded linear operators with a dominant eigenspace. We omit the details of the proof, as they are identical to those found in the proof of Lemma 3.1 and Lemma 3.2 of [27].

Theorem 5.1**.**

Let $\mathcal{B}$ be a bounded linear operator on a Hilbert space $\mathcal{H}$ with dominant eigenspace $\mathcal{V}$, defined in Eq. 22, having ${\rm dim}(\mathcal{V})=m$. Select a quasimatrix $F:\mathbb{C}^{m}\rightarrow\mathcal{H}$ such that the columns of $\mathcal{P}_{\mathcal{V}}F$ are linearly independent and suppose the columns of the quasimatrix $\hat{Q}_{k}:\mathbb{C}^{m}\rightarrow\mathcal{H}$ form an orthonormal basis for ${\rm range}(\mathcal{B}^{k}F)$, for $k=1,2,3,\dots$. If $u\in\mathcal{V}$ is an eigenvector of $\mathcal{B}$ with eigenvalue $\lambda$, then there is a function $\hat{u}_{k}\in{\rm range}(\hat{Q}_{k})$ such that

[TABLE]

*where $\rho=\rho((\mathcal{I-P}_{\mathcal{V}})\mathcal{B})$, $\epsilon_{k}\rightarrow 0$ as $k\rightarrow\infty$, and $u=\mathcal{P}_{\mathcal{V}}Fx$. *

Although we have neglected the effects of approximately solving the ODEs in Eq. 10 and the impact of round-off errors in our brief analysis of rational subspace iteration for normal differential operators, we mention two recent results for rational subspace iteration with matrices [48] and self-adjoint differential operators [27, 28].

• •

For matrices, small errors made during application of the spectral projector generally do not alter the convergence behavior of subspace iteration [48]. In this case, the sequence $\hat{Q}_{k}$ no longer converges to an exact basis for $\mathcal{V}$. However, the matrices $\hat{Q}_{k}$ approximate a basis for $\mathcal{V}$ and the approximation error converges geometrically to a constant determined by the sizes of the errors introduced at each iteration [48].

• •

For self-adjoint differential operators (closed and densely defined on $\mathcal{H}$), rational subspace iteration converges to a subspace even when the resolvent operator is discretized to solve the ODEs in Eq. 10 [27, 28]. The distance between the computed subspace and the target eigenspace (in a distance metric between subspaces) is proportional to the approximation error in the discretized resolvent [27, 28].

We expect that similar statements hold for normal operators on $\mathcal{H}$, but a rigorous and detailed convergence analysis is more subtle and beyond the scope of this paper.

5.2 A pseudospectral inclusion theorem

As ${\rm range}(\hat{Q})$ is not an invariant subspace of $\mathcal{L}$, the $\epsilon$-pseudospectrum of $\hat{Q}^{*}\mathcal{L}\hat{Q}$ is not, in general, contained in the $\epsilon$-pseudospectrum of $\mathcal{L}$. However, if $\lVert\hat{Q}-Q\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}$ is sufficiently small, then the $\epsilon$-pseudospectrum of $\hat{Q}^{*}\mathcal{L}\hat{Q}$ is contained in the $2\epsilon$-pseudospectrum of $\mathcal{L}$.

Theorem 5.2**.**

Consider a closed operator $\mathcal{L}$ with domain $\mathcal{D}(\mathcal{L})$ that is densely defined on a Hilbert space $\mathcal{H}$ and fix $\epsilon>0$. Let $Q:\mathbb{C}^{m}\rightarrow\mathcal{D}(\mathcal{L})\bigcap\mathcal{D}(\mathcal{L}^{*})$ satisfy $Q^{*}Q=I$, and let ${\rm range}(Q)$ be an $m$-dimensional invariant subspace of $\mathcal{L}$. If $\hat{Q}:\mathbb{C}^{m}\rightarrow\mathcal{D}(\mathcal{L})$ satisfies

[TABLE]

*then $\lambda_{\epsilon}(\hat{Q}^{*}\mathcal{L}\hat{Q})\subset\lambda_{2\epsilon}(\mathcal{L})$. *

Proof 5.3**.**

Consider $z\in\lambda_{\epsilon}(\hat{Q}^{*}\mathcal{L}\hat{Q})$. If $z\in\lambda_{\epsilon}(\hat{Q}^{*}\mathcal{L}\hat{Q})\cap\lambda_{\epsilon}(\mathcal{L})$, there is nothing to prove, so assume without loss of generality that $z\not\in\lambda_{\epsilon}(\mathcal{L})$. If we denote $R_{Q}(z)=(zI-Q^{*}\mathcal{L}Q)^{-1}$, $R_{\hat{Q}}(z)=(zI-\hat{Q}^{*}\mathcal{L}\hat{Q})^{-1}$, and $E=\hat{Q}-Q$, then we have that $R_{\hat{Q}}(z)=[R_{Q}(z)^{-1}-B]^{-1}$, where $B=Q^{*}\mathcal{L}E+E^{*}\mathcal{L}Q+E^{*}\mathcal{L}E$. Employing a formula for the inverse of the sum of two matrices, we obtain $R_{\hat{Q}}(z)=R_{Q}(z)+R_{Q}(z)[I-BR_{Q}(z)]^{-1}BR_{Q}(z)$ [31].

Now, $\lVert B\rVert_{\mathbb{C}^{m}}\leq\lVert Q^{*}\mathcal{L}E\rVert_{\mathbb{C}^{m}}+\lVert E^{*}\mathcal{L}Q\rVert_{\mathbb{C}^{m}}+\lVert E^{*}\mathcal{L}E\rVert_{\mathbb{C}^{m}}$. Since $\lVert E^{*}\rVert_{\mathcal{H}\rightarrow\mathbb{C}^{m}}=\lVert E\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}$ and $\lVert Q^{*}\mathcal{L}\rVert_{\mathcal{H}\rightarrow\mathbb{C}^{m}}=\lVert\mathcal{L}^{*}Q\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}$ [32, p. 256], our hypothesis indicates that the sum of the three terms comprising $B$ is bounded in norm by

[TABLE]

Moreover, since $z\not\in\lambda_{\epsilon}(\mathcal{L})$, we have that $\lVert R_{Q}(z)\rVert_{\mathbb{C}^{m}}\leq 1/\epsilon$ by Theorem 3.3. Therefore, $\lVert BR_{Q}(z)\rVert_{\mathbb{C}^{m}}\leq 1/2$.

Because $\lVert BR_{Q}(z)\rVert_{\mathbb{C}^{m}}\leq 1/2$, we may use the Neumann series to compute $(I-BR_{Q}(z))^{-1}=\sum_{k=0}^{\infty}(BR_{Q}(z))^{k}$. We see that $R_{\hat{Q}}(z)=R_{Q}(z)\left(I+\sum_{k=1}^{\infty}(BR_{Q}(z))^{k}\right)$ and therefore,

[TABLE]

*Now, if $z\in\lambda_{\epsilon}(\hat{Q}^{*}\mathcal{L}\hat{Q})$, then $\lVert R_{Q}(z)\rVert_{\mathbb{C}^{m}}\geq\lVert R_{\hat{Q}}(z)\rVert_{\mathbb{C}^{m}}/2>1/(2\epsilon)$. By Theorem 3.3, we have that $\lVert(zI-\mathcal{L})^{-1}\rVert_{\mathcal{H}}\geq\lVert R_{Q}(z)\rVert_{\mathbb{C}^{m}}$. Collecting inequalities yields the result $\lVert(zI-\mathcal{L})^{-1}\rVert_{\mathcal{H}}>1/(2\epsilon)$, i.e., $z\in\lambda_{2\epsilon}(\mathcal{L})$. *

A consequence of Theorem 5.2 is that Algorithm 3 possesses a type of stability provided that $\mathcal{L}$ is uniformly bounded on the sequence $E_{1},E_{2},E_{3},\dots$, where $E_{k}=\hat{Q}_{k}-Q$ for $k\geq 1$. If $\mathcal{L}$ is uniformly bounded on $\{E_{k}\}_{k=1}^{\infty}$, then there is a $\Lambda\geq 0$ such that ${\rm sup}_{k\geq 1}\lVert\mathcal{L}E_{k}\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}\leq\Lambda$. Applying Theorem 5.2, we see that Algorithm 3 computes elements in the $2\epsilon$-pseudospectrum of $\mathcal{L}$ provided that a basis for $\mathcal{V}$ is resolved to within $\epsilon/(2(\lVert\mathcal{L}^{*}Q\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}+\lVert\mathcal{L}Q\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}+\Lambda))$.

We now verify, with two mild constraints placed on the choice of the initial quasimatrix $F$, that $\mathcal{L}$ is uniformly bounded on the sequence $\{\hat{Q}_{k}\}_{k=1}^{\infty}$ generated by Algorithm 3. Note that this implies that $\mathcal{L}$ is uniformly bounded on $\{E_{k}\}_{k=1}^{\infty}$ because $E_{k}=\hat{Q}_{k}-Q$ and ${\rm range}(Q)\subset\mathcal{D}(\mathcal{L})$. The constraints on $F$ are generically satisfied when $F$ is selected as in Section 4. In the statement of the bound on $\lVert\mathcal{L}\hat{Q}_{k}\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}$, we use the notation $\sigma_{\min}(\mathcal{P}_{\mathcal{V}}F)$ and $\sigma_{\min}((\mathcal{I-P}_{\mathcal{V}})F)$ to denote the smallest singular values of the quasimatrices $\mathcal{P}_{\mathcal{V}}F$ and $(\mathcal{I-P}_{\mathcal{V}})F$, respectively.999The singular value decomposition of a quasimatrix $A:\mathbb{C}^{m}\rightarrow\mathcal{H}$ is the decomposition $A=U\Sigma V^{*}$, where $U:\mathbb{C}^{m}\rightarrow\mathcal{H}$ is a quasimatrix with $\mathcal{H}$-orthonormal columns, $\Sigma\in\mathbb{C}^{m\times m}$ is a diagonal matrix with non-negative entries $\sigma_{1}\geq\dots\geq\sigma_{m}$, and $V\in\mathbb{C}^{m\times m}$ is a unitary matrix [57].

Lemma 5.4**.**

Consider a closed, normal operator $\mathcal{L}$ with domain $\mathcal{D}(\mathcal{L})$ that is densely defined on a Hilbert space $\mathcal{H}$. Let $\mathcal{\hat{P}}_{\mathcal{V}}$ be the bounded operator on $\mathcal{H}$ defined in Eq. 9 and suppose that $\mathcal{\hat{P}}_{\mathcal{V}}$ has a dominant eigenspace of $\mathcal{V}$ (see Eq. 22) with ${\rm dim}(\mathcal{V})=m$. Let $F$, $\mathcal{P}_{\mathcal{V}}$, and $\{\hat{Q}_{k}\}_{k=1}^{\infty}$ be as in Theorem 5.1 with $\mathcal{B}=\mathcal{\hat{P}}_{\mathcal{V}}$. Suppose that $\mathcal{\hat{P}}_{\mathcal{V}}^{k}F$ (for each $k\geq 1$) and $(\mathcal{I-P}_{\mathcal{V}})F$ each have linearly independent columns and that ${\rm range}(F)\subset\mathcal{D}(\mathcal{L})$. Then, we have that

[TABLE]

*where $M=\max\left\{1/\sigma_{\min}(\mathcal{P}_{\mathcal{V}}F),1/\sigma_{\min}((\mathcal{I-P}_{\mathcal{V}})F)\right\}$. *

Proof 5.5**.**

Since $\hat{Q}_{k}$ is an orthonormal basis for $\mathcal{\hat{P}}_{\mathcal{V}}^{k}F$, there is a matrix $R_{k}\in\mathbb{C}^{m\times m}$ such that $\mathcal{\hat{P}}_{\mathcal{V}}^{k}F=\hat{Q}_{k}R_{k}$. By the assumption that $\mathcal{\hat{P}}_{\mathcal{V}}^{k}F$ has linearly independent columns, we know that $R_{k}$ is invertible. We obtain that

[TABLE]

We use the spectral projector $\mathcal{P}_{\mathcal{V}}$ to rewrite Eq. 23 as

[TABLE]

Now, ${\rm range}(\mathcal{P}_{\mathcal{V}}F)$ and ${\rm range}((\mathcal{I-P}_{\mathcal{V}})F)$ are invariant under $\mathcal{\hat{P}}_{\mathcal{V}}$ [32, p. 178]. Consequently, there are matrices $D_{1},D_{2}\in\mathbb{C}^{m\times m}$ such that

[TABLE]

Substituting Eq. 25 into Eq. 24 yields the following useful equation for $\hat{Q}_{k}$:

[TABLE]

Applying $\mathcal{L}$ to both sides of Eq. 26 and commuting with the spectral projectors $\mathcal{P}_{\mathcal{V}}$ and $\mathcal{I-P}_{\mathcal{V}}$ [32, p. 179], we obtain

[TABLE]

Since ${\rm range}(F)\subset\mathcal{D}(\mathcal{L})$, we have that $\lVert\mathcal{L}F\rVert_{\mathbb{C}^{m}\rightarrow\mathcal{H}}<\infty$. Additionally, since $\mathcal{L}$ is normal, the spectral projectors have norms equal to $1$ [32, p. 277]. Therefore, it remains to find a uniform bound for $\lVert D_{1}^{k}R_{k}^{-1}\rVert_{\mathbb{C}^{m}}$ and $\lVert D_{2}^{k}R_{k}^{-1}\rVert_{\mathbb{C}^{m}}$ as $k\rightarrow\infty$.

For brevity, we prove uniform boundedness of $\lVert D_{1}^{k}R_{k}^{-1}\rVert_{\mathbb{C}^{m}}$ and note that the proof for $\lVert D_{2}^{k}R_{k}^{-1}\rVert_{\mathbb{C}^{m}}$ is essentially identical. We begin by commuting $\mathcal{\hat{P}}_{\mathcal{V}}$ with the spectral projectors in Eq. 25 and substituting the QR factorization of $\mathcal{\hat{P}}_{\mathcal{V}}^{k}F$ to see that

[TABLE]

Using the psuedoinverse $(\mathcal{P}_{\mathcal{V}}F)^{+}$ of the quasimatrix101010The pseudoinverse of a quasimatrix $A:\mathbb{C}^{m}\rightarrow\mathcal{H}$ may be defined via the SVD as $A^{+}=V\Sigma^{+}U^{*}$, where $\Sigma^{+}$ is the diagonal matrix with entries $\Sigma^{+}_{ii}=1/\sigma_{i}$ if $\sigma_{i}\neq 0$ and [math] otherwise. It is easy to verify familiar properties from the matrix case [25, p. 290], i.e., if $A$ has linearly independent columns, then $A^{+}A=I$ and $\lVert A^{+}\rVert_{\mathcal{H}\rightarrow\mathbb{C}^{m}}=1/\sigma_{\min}(A)$. $\mathcal{P}_{\mathcal{V}}F$ and noting that $\mathcal{P}_{\mathcal{V}}F$ has linearly independent columns, Eq. 28 implies that

[TABLE]

Now, we know that $\lVert\mathcal{P}_{\mathcal{V}}\hat{Q}_{k}\rVert_{\mathcal{H}}\leq 1$, because $\lVert\mathcal{P}_{\mathcal{V}}\rVert_{\mathcal{H}}=1$ and $\hat{Q}_{k}$ has orthonormal columns. We conclude that

[TABLE]

A similar argument shows that

[TABLE]

Taking norms in Eq. 27 and substituting the bounds from Eq. 30 and Eq. 31, we find

[TABLE]

*The lemma follows immediately from Eq. 32. *

Theorem 5.1, Theorem 5.2, and Lemma 5.4 provide a preliminary analysis to explain why Algorithm 3 accurately computes the eigenvalues of normal operators with a dominant eigenspace $\mathcal{V}$. Theorem 5.1 allows us to accurately resolve an orthonormal basis $Q$ for $\mathcal{V}$ by refining the quasimatrix $\hat{Q}_{k}$ with subspace iteration. Lemma 5.4 confirms that $\mathcal{L}\hat{Q}_{k}$ does not grow without bound as $\hat{Q}_{k}$ is refined. Finally, Theorem 5.2 demonstrates that the eigenvalues are computed to the expected accuracy, provided that the basis for $\mathcal{V}$ has been resolved.

6 An operator analogue of the Rayleigh Quotient Iteration

It is useful to have operator analogues for other eigensolvers too; particularly, when the eigenvalues of interest are difficult to target with a pre-selected search region $\Omega\subset\mathbb{C}$. The Rayleigh Quotient Iteration (RQI) is a generalization of the inverse iteration that incorporates dynamic shifting to obtain cubic (for Hermitian problems) or quadratic (non-Hermitian problems) convergence [43]. Given a matrix $A\in\mathbb{C}^{n\times n}$ and an initial vector $\tilde{y}_{0}\in\mathbb{C}^{n}$, RQI computes the iterates

[TABLE]

The vectors $y_{k}$ typically converge to a nearby eigenvector of $A$, while the sequence $\beta_{k}$ converges to the associated eigenvalue of $A$ [41]. In the matrix setting, Eq. 33 is often used to compute interior eigenvalues or refine an estimate of an invariant subspace [42, 43].

Replacing a matrix $A$ by a differential operator $\mathcal{L}:\mathcal{D}(\mathcal{L})\rightarrow\mathcal{H}$, as in Eq. 1, and the vectors $\tilde{y}_{k}$ by functions $f_{k}\in\mathcal{D}(\mathcal{L})$, we obtain an operator analogue of RQI. One needs to select an initial function $f_{0}\in\mathcal{D}(\mathcal{L})$ and solve a sequence of ODEs, i.e.,

[TABLE]

At each iteration, the shift $\beta_{k}$ is computed from the Rayleigh Quotient $(f_{k},\mathcal{L}f_{k})_{\mathcal{H}}$ (in strong form) and the solution $f_{k+1}$ is normalized after each iteration. Analogous to the matrix setting, we observe that the operator analogue of the Rayleigh Quotient Iteration converges cubically for self-adjoint operators and quadratically otherwise [29].

We note that block generalizations of RQI (RSQR and GRQI [1]) are also easily extended to the differential operator setting. In this case, a sequence of quasimatrices $\hat{Q}_{k}$ with $\mathcal{H}$-orthonormal columns are generated to approximate an invariant subspace of $\mathcal{L}$ and a Rayleigh–Ritz projection is performed to compute approximate eigenvalues and eigenvectors. As with the operator analogue of FEAST, Theorem 5.2 implies that the iteration Eq. 34 accurately computes eigenvalues of normal differential operators when the basis for the target eigenspace is sufficiently resolved.

6.1 Free vibrations of an airplane wing

The improved convergence rate of RQI can offer much faster computation time than subspace iteration, often requiring only $3$ or $4$ ODE solves to reach an accuracy of essentially machine precision [29]. We now employ Eq. 34 for the rapid computation of vibrational modes of an airplane wing.

An airplane wing may be crudely modeled as a thin, cantilevered beam of length $L$ with a linear taper. The governing equation for free vibrations is [30]

[TABLE]

The variable coefficient $1+x$ accounts for the linear taper of the wing, while the boundary conditions on $u^{\prime\prime}$ and $u^{\prime\prime\prime}$ at $x=1$ express the natural requirement that the bending moment and shear force vanish at the endpoint.

To compute a few of the smoothest modes of Eq. 35 we use the eigenfunctions $w_{n}$ of the cantilevered beam equation with constant coefficients, given in closed form by [30]

[TABLE]

Here $\beta_{n}$ is the $n$th root of $g(\beta)=\cosh(\beta L)\cos(\beta L)+1$ [30]. We target a mode of Eq. 35 by setting $f_{0}(x)=w_{n}(x)$. Figure 5 shows the modes that are computed using initial guesses $w_{1},\dots,w_{4}$, corresponding to the smallest four positive roots of $g(\beta)$.

7 Computing eigenvalues in unbounded regions

The stability analysis of solutions to time-dependent partial differential equations (PDEs) provides an abundant source of differential eigenvalue problems. Consider the initial boundary value problem (IBVP) with periodic boundary conditions

[TABLE]

Here, $\mathcal{L}$ and $\mathcal{N}$ are linear and nonlinear ordinary differential operators (with respect to the variable $x$), respectively. In many instances, Eq. 37 supports steady-states, traveling wave states, or other phenomena whose stability is of critical importance in the physical problem under study [35, 4, 49]. When $\mathcal{L}$ is self-adjoint or normal, the stability analysis often reduces to determining whether or not the eigenvalues of $\mathcal{L}$ are contained in one half-plane [34, 4, 49, 60]. We now show how to modify the spectral projector in Eq. 5 to derive a practical rational filter to compute (finitely many) eigenvalues of $\mathcal{L}$ in the right half-plane.

7.1 A rational filter for the half-plane

Let $\mathcal{L}$ be a closed linear operator that is densely defined on a Hilbert space $\mathcal{H}$. Suppose that $\mathcal{L}$ is a normal operator with a spectrum in the left half-plane ${\rm Re}(z)<0$ except for finitely many eigenvalues $\lambda_{1},\dots,\lambda_{m}$ (including multiplicities) such that ${\rm Re}(\lambda_{i})>0$ for $1\leq i\leq m$. Denote the eigenspace associated with $\lambda_{1},\dots,\lambda_{m}$ by $\mathcal{V}$ and consider search regions that are semi-circles of radius $R$, i.e.,

[TABLE]

To construct a computable spectral projector onto the right half-plane we consider taking $R\rightarrow\infty$. We adopt the following strategy:

• (i)

Introduce a $1/R$ decay into the integrand of the spectral projector Eq. 5 as $R\rightarrow\infty$, while preserving the projection onto $\mathcal{V}$.

• (ii)

Split the projector into an integral over the vertical part of $\partial\Omega_{R}$ and an integral over the circular arc of $\partial\Omega_{R}$. By taking $R\rightarrow\infty$, we observe that the contribution from the circular arc goes to [math] due to the additional $1/R$ decay in the integrand.

• (iii)

Map the imaginary axis to the interval $[-1,1]$ and approximate the spectral projector by a quadrature rule.

Select $a\in\mathbb{R}^{+}$ and consider the family of functions that are analytic in the right half-plane defined by

[TABLE]

By Cauchy’s Integral Formula, we know that $\mathcal{P}_{R}(\lambda)=(\lambda+a)^{-1}$ if $\lambda\in\Omega_{R}$ and is zero otherwise [51]. Taking the limit $R\rightarrow\infty$, we obtain

[TABLE]

Using functional calculus for unbounded normal operators we can extend $\mathcal{P}(\lambda)$ to an operator-valued function $\mathcal{P}(\mathcal{L})$ [46, Theorem 13.24].111111In [46, Theorem 13.24], $E_{x,y}$ is the spectral measure of $\mathcal{L}$ [46, Theorem 13.33]. Moreover, we have that $\mathcal{P}(\mathcal{L})u=\mathcal{P}(\lambda)u$ when $\mathcal{L}u=\lambda u$ [45, VIII.5]. Consequently, ${\rm range}(\mathcal{P}(\mathcal{L}))=\mathcal{V}$.

Now, take the change-of-variables $x=\frac{2}{\pi}\tan^{-1}y$ in Eq. 40 to obtain

[TABLE]

Using Gauss–Legendre quadrature nodes $x_{1},\dots,x_{\ell}$ and weights $w_{1},\dots,w_{\ell}$ on $[-1,1]$, we can approximate $\mathcal{P}(\mathcal{L})$ by

[TABLE]

Figure 6 (right) shows the derived rational filter $\mathcal{\hat{P}}(\lambda)$ in the complex plane.

7.2 Stability of thin fluid films

To demonstrate the utility of the filter in Eq. 42, we assess the stability of the steady-state solutions to a PDE governing the motion of a thin film of fluid supported below by a flat substrate. The PDE is

[TABLE]

where $u(x,t)$ is a positive, periodic function representing the thickness of the fluid [35]. The nonlinear term models gravitational effects and substrate-fluid interactions [35].

A droplet steady-state $u_{\text{ss}}(x)$ of Eq. 43, rescaled so that it is supported on $[0,l]$ with contact angle $\pi/4$, is stable if all the eigenvalues of a fourth-order differential operator are in the left half-plane. The associated differential eigenproblem is [34]

[TABLE]

We compute the steady-state $u_{\text{ss}}(x)$ by solving the second order nonlinear ODE [35]

[TABLE]

Here, $\delta$ is a dimensionless quantity relating the rescaled problem to the original contact angle [35]. The length $l$ of the droplet’s base and $\delta$ may be calculated analytically [35].

In Fig. 7, we show an approximation to the rescaled steady-state $u_{\text{ss}}$ along with the right-most eigenvalues of Eq. 44. Using the rational filter in Eq. 42 with $\ell=20$ (the degree of the quadrature rule defining the filter) to perform the approximate spectral projection in Algorithm 3, we are able to identify an eigenvalue of Eq. 44 in the right half-plane, which indicates that the droplet (see Fig. 7 (left)) is unstable.

Techniques for selecting the dimension $m$ of the subspace $\mathcal{V}$ [54, 33] are important in stability analysis as one is trying to determine the number of eigenvalues in the right half-plane. To select $m$, we monitor the singular values of the matrix $\hat{V}^{*}\hat{V}$ after each iteration and adjust the number of basis functions by removing columns of $\hat{V}$ associated with singular values that are close to machine precision (relative to the largest singular value) [54]. This procedure usually allows us to capture the dominant eigenspace of the filtered operator $\mathcal{\hat{P}}(\mathcal{L})$ that includes the target eigenspace as well as any eigenvalues clustered near the imaginary axis. We then determine whether there are any eigenvalues in the right half-plane by sorting through the computed eigenvalues. However, this strategy may break down, for instance, if there is an eigenvalue close to a quadrature node. Additionally, the sharp decay of the filter Eq. 42 across the imaginary axis is softened as $\lvert{\rm Im}(z)\rvert\rightarrow\infty$, which can lead to difficulties when there are clusters of eigenvalues near the imaginary axis with large imaginary part. In this case, one may need to take a large number of basis functions to accurately resolve the dominant eigenvalues of $\mathcal{\hat{P}}(\mathcal{L})$.

Conclusions

An operator analogue of the FEAST matrix eigensolver is derived to solve differential eigenvalue problems without discretizing the operator. This approach leads to an algorithm that can exploit spectrally accurate techinques for computing with functions while preserving the structure of $\mathcal{L}$. The result is an efficient, automated, and accurate eigensolver for normal and self-adjoint differential operators. This eigensolver is adept in the high-frequency regime and may provide a new direction towards robust high-frequency eigenvalue computations.

The implementation described in Section 4 extends to higher dimensions in a straightforward way for simple geometries where spectral methods apply [56]. In the case of more complicated geometries, one may still benefit from the advantages of the “solve-then-discretize” framework outlined in Algorithm 3 provided that one is able to accurately compute inner products and solutions to the shifted linear differential equations. Although we have focused on the strong form of the eigenvalue problem in Eq. 1, it may be necessary to work with the weak form. In this case, one can still follow Algorithm 3 provided that the shifted linear systems are solved in weak form and the Rayleigh–Ritz projection is performed with the associated sesquilinear form.

Acknowledgements

We thank Anthony Austin for providing detailed feedback on the results in Section 3 and on our rational filter in Section 7. We thank David Bindel for helpful discussions about spectral discretizations of differential eigenvalue problems. Finally, we thank Marcus Webb, Sheehan Olver, Yuji Nakatsukasa, Heather Wilber, and Matthew Colbrook for providing comments on a draft.