Fully computable a posteriori error bounds for eigenfunctions

TL;DR

This paper introduces two fully computable a posteriori error bounds for eigenfunction approximation of compact self-adjoint operators, applicable to clustered and multiple eigenvalues, with improved accuracy demonstrated through numerical examples.

Contribution

It proposes two novel algorithms for a posteriori error estimation in eigenfunction approximation, extending existing methods to handle complex eigenvalue problems with enhanced precision.

Findings

Algorithms effectively estimate errors in eigenfunction approximations.

The methods perform well on matrix and Laplace operator problems.

Numerical results confirm the efficiency and accuracy of the proposed algorithms.

Abstract

For compact self-adjoint operators in Hilbert spaces, two algorithms are proposed to provide fully computable a posteriori error estimate for eigenfunction approximation. Both algorithms apply well to the case of tight clusters and multiple eigenvalues, under the settings of target eigenvalue problems. Algorithm I is based on the Rayleigh quotient and the min-max principle that characterizes the eigenvalue problems. The formula for the error estimate provided by Algorithm I is easy to compute and applies to problems with limited information of Rayleigh quotients. Algorithm II, as an extension of the Davis--Kahan method, takes advantage of the dual formulation of differential operators along with the Prager--Synge technique and provides greatly improved accuracy of the estimate, especially for the finite element approximations of eigenfunctions. Numerical examples of eigenvalue problems of matrices and the Laplace operators over convex and non-convex domains illustrate the efficiency of the proposed algorithms.

Tables

Table 1. Table 1: Error estimates of eigenvalues and eigenvectors for the generalized matrix eigenvalue problem

cluster	two-sided bounds	eigenvalue	bound (10)	bound (11)	bound (23)
	on eigenvalues	enclosure size	on ${\delta }_a$	on ${\delta }_b$	on ${\delta }_a$
$\left\{{\lambda }_1\right\}$	${19.78151183}_{16}^{25}$	$8.43\times {10}^{-10}$	$4.10\times {10}^{-7}$	$2.60\times {10}^{-7}$	$6.23\times {10}^{-6}$
$\{{\lambda }_2,{\lambda }_3\}$	${49.5769652}_{58}^{60}$	$1.00\times {10}^{-9}$	$1.13\times {10}^{-6}$	$6.12\times {10}^{-7}$	$8.89\times {10}^{-7}$
$\left\{{\lambda }_4\right\}$	${79.63484906}_{38}^{50}$	$1.10\times {10}^{-9}$	$5.33\times {10}^{-6}$	$1.98\times {10}^{-6}$	$1.20\times {10}^{-6}$
$\{{\lambda }_5,{\lambda }_6\}$	${99.50227701}_{43}^{55}$	$1.16\times {10}^{-9}$	$1.06\times {10}^{-5}$	$3.46\times {10}^{-6}$	$2.04\times {10}^{-6}$

Table 2. Table 2: The leading 8 8 8 clusters for the eigenvalue problem A ​ x = λ ​ x 𝐴 𝑥 𝜆 𝑥 Ax=\lambda x

$E_1$	$E_2$	$E_3$	$E_4$	$E_5$	$E_6$	$E_7$	$E_8$
$\left\{{\lambda }_1\right\}$	$\{{\lambda }_2,{\lambda }_3\}$	$\left\{{\lambda }_4\right\}$	$\{{\lambda }_5,{\lambda }_6\}$	$\{{\lambda }_7,{\lambda }_8\}$	$\{{\lambda }_9,{\lambda }_{10}\}$	$\left\{{\lambda }_{11}\right\}$	$\{{\lambda }_{12},{\lambda }_{13}\}$

Table 3. Table 3: The directed distance δ b ​ ( E k , E ^ k ) subscript 𝛿 𝑏 subscript 𝐸 𝑘 subscript ^ 𝐸 𝑘 \delta_{b}(E_{k},\widehat{E}_{k}) , its bound ( 11 ), and Davis–Kahan’s estimate ( 46 ) for clusters k = 1 , 2 , … , 8 𝑘 1 2 … 8 k=1,2,\dots,8 .

Cluster	$\left\{{\lambda }_1\right\}$	$\{{\lambda }_2,{\lambda }_3\}$	$\left\{{\lambda }_4\right\}$	$\{{\lambda }_5,{\lambda }_6\}$
${\delta }_b(E_k,{\hat{E}}_k)$	1.26 $\times {10}^{-3}$	2.39 $\times {10}^{-3}$	1.78 $\times {10}^{-3}$	2.16 $\times {10}^{-3}$
Bound (11)	$1.08\times {10}^{-2}$	$2.12\times {10}^{-2}$	$4.36\times {10}^{-2}$	$7.24\times {10}^{-2}$
Davis–Kahan (46)	$1.36\times {10}^{-1}$	$1.41\times {10}^{-1}$	$2.04\times {10}^{-1}$	$2.14\times {10}^{-1}$

Table 4. Table 4: The four leading clusters for the unit square.

Cluster	Eigenvalues
1	${\lambda }_1=2{\pi }^2$
2	${\lambda }_2={\lambda }_3=5{\pi }^2$
3	${\lambda }_4=8{\pi }^2$
4	${\lambda }_5={\lambda }_6=10{\pi }^2$

Table 5. Table 5: Effectivity indices of estimates ( 10 ), ( 23 ) and ( 11 ) (without the residual error estimation; the unit square)

mesh size	${\delta }_a(E_1,{\hat{E}}_1)$		${\delta }_a(E_2,{\hat{E}}_2)$		${\delta }_a(E_3,{\hat{E}}_3)$		${\delta }_a(E_4,{\hat{E}}_4)$
	(10)	(23)	(10)	(23)	(10)	(23)	(10)	(23)
h=1/32	1.29	1.29	2.53	1.81	5.08	2.98	9.41	4.80
h=1/64	1.29	1.29	2.55	1.81	5.08	2.96	9.42	4.72
h=1/128	1.29	1.29	2.55	1.81	5.08	2.96	9.42	4.70

Table 6. Table 6: Effectivity indices of estimation from Theorem 5.3 (with the residual error estimation; the unit square)

mesh size	cluster 1		cluster 2		cluster 3		cluster 4
	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$
h=1/32	1.00045	1.27	1.0065	3.12	1.078	8.189	1.12	10.385
h=1/64	1.00011	1.27	1.0016	3.12	1.018	7.76	1.026	9.53
h=1/128	1.00003	1.27	1.0004	3.12	1.0043	7.66	1.0064	9.36

Table 7. Table 7: Comparison of the effectivity index for different approaches applied to δ ¯ ​ ( E 1 , E ^ 1 ) ¯ 𝛿 subscript 𝐸 1 subscript ^ 𝐸 1 \overline{\delta}(E_{1},\widehat{E}_{1}) (the same type uniform mesh as CanDusMadStaVoh2017 for the unit square domain )

$h$	$\overline{\delta }$	$ℐ\left(\overline{\delta }\right)$ from Thm 3.1	$ℐ\left(\overline{\delta }\right)$ from Thm 5.3	$ℐ\left(\overline{\delta }\right)$ from (45)	$ℐ\left(\overline{\delta }\right)$ of CanDusMadStaVoh2017
1/20	$0.3494$	$1.29$	1.00035	1.67	1.13
1/40	$0.1745$	1.29	1.00008	1.67	1.07

Table 8. Table 8: Lower bounds on the leading eigenvalues for the L-shaped domain.

${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\lambda }_4$	${\lambda }_5$
${9.6397}_1^3$	${15.1972}_5^6$	${19.7392}_0^1$	${29.5214}_7^9$	${31.9126}_2^4$

Table 9. Table 9: Effectivity index of estimation on δ a subscript 𝛿 𝑎 \delta_{a} and δ b subscript 𝛿 𝑏 \delta_{b} (standard FEM; L-shaped domain)

(a) Estimate of ${\delta }_a$ and ${\delta }_b$ by (23) and (11) (no residual error estimation)
mesh size	cluster 1		cluster 2		cluster 3		cluster 4
	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$
h=1/32	1.67	24.99	3.57	159.94	4.48	178.45	12.74	394.83
h=1/64	1.69	38.71	4.03	365.49	5.056	407.52	13.83	898.22

Table 10. Table 10: Comparison of the effectivity index for different approaches (L-shaped domain)

Approaches	Mesh type	#(DOF)	$ℐ\left(\overline{\delta }(E_1,{\hat{E}}_1)\right)$
Theorem 3.1	uniform mesh ($h=1/64$)	12545	1.69
Theorem 5.3	uniform mesh ($h=1/64$)	12545	1.024
Estimate of CanDusMadStaVoh2017	unstructured mesh	24925	2.51

Table 11. Table 11: Lower and upper bounds of eigenvalues for the dumbbell shaped domain. Two times refined initial mesh and third order finite element spaces were used.

cluster	lower and upper eigenvalue bounds
1	${\lambda }_1={19.736}_{634}^{729},{\lambda }_2={19.736}_{635}^{729}$
2	${\lambda }_3={49.33}_{761}^{809},{\lambda }_4={49.33}_{761}^{809},{\lambda }_5={49.348020}_5^8,{\lambda }_6={49.348020}_5^8$
3	${\lambda }_7={78.9568}_{290}^{301},{\lambda }_8={78.9568}_{290}^{301}$
4	${\lambda }_9={98.6}_{69041}^{71154},{\lambda }_{10}={98.6}_{69041}^{71154},{\lambda }_{11}={98.69604}_{39}^{41},{\lambda }_{12}={98.69604}_{39}^{41}$

Table 12. Table 12: Effectivity index of estimation on δ a subscript 𝛿 𝑎 \delta_{a} and δ b subscript 𝛿 𝑏 \delta_{b} (Dumbbell-shaped domain)

(a) Estimate of ${\delta }_a$ and ${\delta }_b$ by (23) and (11) (no residual error estimation)
mesh size	cluster 1		cluster 2		cluster 3		cluster 4
	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$	$ℐ\left({\delta }_a\right)$	$ℐ\left({\delta }_b\right)$
h=0.068	1.29	29.94	2.41	29.64	5.10	40.27	8.32	42.96

Equations

$$a(u_i,v)={\lambda }_{i}b(\gamma u_i,\gamma v)\forall v\in V.$$

$$a(Tf,v)=b(f,\gamma v)\forall v\in V.$$

$$b(\gamma u_i,\gamma v)={\mu }_{i}a(u_i,v)\forall v\in V.$$

$$0\le \underset{i\to d}{\lim }{\mu }_i\le \cdots \le {\mu }_2\le {\mu }_1,$$

$$0<{\lambda }_1\le {\lambda }_2\le \cdots ,\text{and}a(u_i,v)={\lambda }_{i}b(\gamma u_i,\gamma v)\forall v\in V,$$

$$a(u_i,u_j)={\delta }_{ij},i,j=1,2,\dots ,$$

$$v=\sum _{i=1}^{d}a(v,u_i)u_i+v_0,\text{where }{v}_0\in \mathrm{Ker}(T\gamma ).$$

$$\Vert v{\Vert }_a^2=\sum _{i=1}^d|a(v,u_i){|}^2+\Vert v_0{\Vert }_a^2\text{and}\Vert \gamma v{\Vert }_b^2=\sum _{i=1}^d|a(v,u_i){|}^2/{\lambda }_i\forall v\in V.$$

$$\left\{\begin{array}[]{l}\displaystyle V=H_{0}^{1}(\Omega),\quad a(u,v)=\int_{\Omega}\nabla u\cdot\nabla v\,\mbox{d}x;\\ \displaystyle W=L^{2}(\Omega),\quad b(u,v)=\int_{\Omega}uv\,\mbox{d}x;\\ \gamma\mbox{ being the identity operator}.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle V=H^{1}(\Omega),\quad a(u,v)=\int_{\Omega}\nabla u\cdot\nabla v+uv\,\mbox{d}x;\\ \displaystyle W=L^{2}(\partial\Omega),\quad b(u,v)=\int_{\partial\Omega}uv\,\mbox{d}s;\\ \gamma\mbox{ being the trace operator}.\end{array}\right.$$

$${\delta }_a(E,\hat{E})=\underset{\begin{array}{c}v\in E\\ \Vert v{\Vert }_a=1\end{array}}{\max }\underset{\hat{v}\in \hat{E}}{\min }\Vert v-\hat{v}{\Vert }_a\text{and}{\delta }_b(E,\hat{E})=\underset{\begin{array}{c}v\in E\\ \Vert \gamma v{\Vert }_b=1\end{array}}{\max }\underset{\hat{v}\in \hat{E}}{\min }\Vert \gamma (v-\hat{v}){\Vert }_b.$$

$${\delta }_a(E,\hat{E})$$

$${\delta }_a^2(E,\hat{E})$$

$$\Vert u-\hat{u}{\Vert }_a^2=\Vert u{\Vert }_a^2+\Vert \hat{u}{\Vert }_a^2-2\Vert u{\Vert }_a\Vert \hat{u}{\Vert }_a\sqrt{1-{\delta }_a^2(E,\hat{E})}.$$

$${\hat{\epsilon }}_a(E,E^{\prime })=\underset{\begin{array}{c}v\in E\\ \Vert v{\Vert }_a=1\end{array}}{\max }\underset{\begin{array}{c}{v}^{\prime }\in E^{\prime }\\ \Vert v^{\prime }{\Vert }_a=1\end{array}}{\max }a(v,v^{\prime })\text{and}{\hat{\epsilon }}_b(E,E^{\prime })=\underset{\begin{array}{c}v\in E\\ \Vert \gamma v{\Vert }_b=1\end{array}}{\max }\underset{\begin{array}{c}{v}^{\prime }\in E^{\prime }\\ \Vert \gamma v^{\prime }{\Vert }_b=1\end{array}}{\max }b(\gamma v,\gamma v^{\prime }).$$

$$F={\left(a(v_i,v_j^{\prime })\right)}_{m\times m^{\prime }},G={\left(a(v_i,v_j)\right)}_{m\times m},H={\left(a(v_i^{\prime },v_j^{\prime })\right)}_{m^{\prime }\times m^{\prime }}.$$

$${\hat{\epsilon }}_a^2(E,E^{\prime })={\lambda }_{max}(F^T{G}^{-1}F,H)={\lambda }_{max}(F{H}^{-1}{F}^T,G),$$

$${\hat{\epsilon }}_a^2(E,E^{\prime })\le \frac{{\eta }_F}{(1-{\eta }_G)(1-{\eta }_H)}.$$

$$a(v,v^{\prime })={\mathbf{c}}^{T}F{\mathbf{c}}^{\prime },\Vert v{\Vert }_a^2={\mathbf{c}}^{T}G\mathbf{c},\text{and}\Vert v^{\prime }{\Vert }_a^2=({\mathbf{c}}^{\prime }{)}^{T}H{\mathbf{c}}^{\prime }.$$

$${\hat{\epsilon }}_a(E,E^{\prime })=\underset{{\mathbf{c}}^{T}G\mathbf{c}=1}{\max }\underset{({\mathbf{c}}^{\prime }{)}^{T}H{\mathbf{c}}^{\prime }=1}{\max }{\mathbf{c}}^{T}F{\mathbf{c}}^{\prime }=\underset{{\mathbf{c}}^{T}G\mathbf{c}=1}{\max }\underset{|{\tilde{\mathbf{c}}}^{\prime }|=1}{\max }{\mathbf{c}}^{T}F{L}^{-T}{\tilde{\mathbf{c}}}^{\prime }=\underset{{\mathbf{c}}^{T}G\mathbf{c}=1}{\max }|{\mathbf{c}}^{T}F{L}^{-T}|,$$

$${\hat{\epsilon }}_a^2(E,E^{\prime })=\underset{0\ne \mathbf{c}\in {\mathbb{R}}^m}{\max }\frac{{\mathbf{c}}^{T}F{L}^{-T}{L}^{-1}{F}^T\mathbf{c}}{{\mathbf{c}}^{T}G\mathbf{c}}={\lambda }_{max}(F{H}^{-1}{F}^T,G).$$

$${\lambda }_{max}(F{H}^{-1}{F}^T,G)=\Vert Q^{-1}F{H}^{-1}{F}^T{Q}^{-T}{\Vert }_2\le \Vert G^{-1}{\Vert }_2\Vert H^{-1}{\Vert }_2\Vert F^{T}F{\Vert }_2.$$

$$\Vert G^{-1}{\Vert }_2=\frac{1}{{\lambda }_{min}(G)}=\frac{1}{1-{\lambda }_{max}(I-G)}\le \frac{1}{1-\Vert I-G{\Vert }_2}\le \frac{1}{1-{\eta }_G}.$$

$${\delta }_a^2(E_K,{\hat{E}}_K)$$

$${\delta }_b^2(E_K,{\hat{E}}_K)$$

$${\hat{\lambda }}_N^{(K)}=\underset{\hat{v}\in {\hat{E}}_K}{\max }\frac{\Vert \hat{v}{\Vert }_a^2}{\Vert \gamma \hat{v}{\Vert }_b^2},{\theta }_a^{(K)}=\sum _{k=1}^{K-1}\frac{\rho -{\lambda }_{n_k}}{{\lambda }_{n_k}}{\left[{\hat{\epsilon }}_a({\hat{E}}_k,{\hat{E}}_K)+{\delta }_a(E_k,{\hat{E}}_k)\right]}^2,\text{and}$$

$${\theta }_b^{(K)}=\sum _{k=1}^{K-1}\left(\rho -{\lambda }_{n_k}\right){\left[{\hat{\epsilon }}_b({\hat{E}}_k,{\hat{E}}_K)+{\delta }_b(E_k,{\hat{E}}_k)\right]}^2.$$

$$\rho \Vert \gamma \hat{u}{\Vert }_b^2-\Vert \hat{u}{\Vert }_a^2+\Vert {\hat{u}}_0{\Vert }_a^2=\sum _{i=1}^d\frac{\rho -{\lambda }_i}{{\lambda }_i}a(\hat{u},u_i{)}^2$$

$$\vartheta (\hat{u})=\sum _{i=1}^{n-1}\frac{\rho -{\lambda }_i}{{\lambda }_i}a(\hat{u},u_i{)}^2,$$

$$\rho \Vert \gamma \hat{u}{\Vert }_b^2-\Vert \hat{u}{\Vert }_a^2-\vartheta (\hat{u})\le \sum _{i=n}^N\frac{\rho -{\lambda }_i}{{\lambda }_i}a(\hat{u},u_i{)}^2.$$

$$\Vert P_k\hat{u}{\Vert }_a^2=\sum _{i=n_k}^{N_k}a(\hat{u},u_i{)}^2,\Vert \gamma P_k\hat{u}{\Vert }_b^2=\sum _{i=n_k}^{N_k}\frac{a(\hat{u},u_i{)}^2}{{\lambda }_i},$$

$$\rho \Vert \gamma \hat{u}{\Vert }_b^2-\Vert \hat{u}{\Vert }_a^2-\vartheta (\hat{u})\le (\rho -{\lambda }_n)\Vert \gamma P_K\hat{u}{\Vert }_b^2\le \frac{\rho -{\lambda }_n}{{\lambda }_n}\Vert P_K\hat{u}{\Vert }_a^2.$$

Full text

††journal: Numerische Mathematik

1 Introduction

The weak form of eigenvalue problems for linear elliptic partial differential operators motivates us to consider an abstract setting of a pair of Hilbert spaces $V$ and $W$ with inner products $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$, respectively, and a compact operator $\gamma:V\to W$ between them. We seek eigenvalues $0<\lambda_{1}\leq\lambda_{2}\leq\cdots$, repeated according to their multiplicity, and corresponding eigenfunctions $u_{i}\in V\setminus\{0\}$, $i=1,2,\dots$ such that

[TABLE]

Section 2 provides more details about the well-posedness of this problem and gives examples how the (generalized) matrix, Laplace, and Steklov eigenvalue problems fit into this setting.

The problem to determine eigenvalues $\lambda_{i}$ is well posed in the sense that small perturbations of the data lead to small perturbations of eigenvalues. However, the variation of eigenfunctions $u_{i}$ upon a small perturbation of the data is not necessarily small, and can even be discontinuous. For example, if two close and simple eigenvalues merge to one multiple eigenvalue then the two corresponding orthogonal eigenfunctions abruptly change into a two dimensional eigenspace. Thus, eigenfunction determination in case of tightly clustered or multiple eigenvalues is an ill-conditioned problem. Any attempt to estimate the error of approximate eigenfunctions has to take into the account this ill-conditioning.

Our approach is to consider the space spanned by eigenfunctions corresponding to all eigenvalues within a cluster. This space is well conditioned provided the cluster is well separated from the rest of the spectrum. We propose error estimators that bound the directed distance (Meyer:2000, , §5.15) between the approximate and the exact space spaces of eigenfunctions in norms induced by both inner products $a$ and $b$, see estimates (10) and (11) below. Further, we present a bound on the distance in the $a$-sense obtained from the distance in the $b$-sense; see (23) below.

The error estimators in our approach provide rigorous upper bounds on these distances without any a priori information about the approximate eigenfunctions. These estimators are fully computable in terms of approximate eigenfunctions and two-sided bounds on eigenvalues and their quality depends on the width of clusters and spectral gaps between them. For the Laplace eigenvalue problem, the proposed estimators generalize the idea from BirBooSwaWen:1966 ; see Remark 5 below.

To compute the needed two-sided bounds on individual eigenvalues, we use the recently developed approach based on the finite element method with explicit error estimation Liu2015 ; Liu2020 (see also, LiuOis2013 ; CarGal2014 ; CarGed2014 ) for the lower bounds on eigenvalues and the Lehmann–Goerisch method Lehmann1949 ; Lehmann1950 ; GoeHau1985 for their high-precision improvements. Note that the Lehmann–Goerisch method should be attributed to T. Kato as well, because his independently developed method Kato1949 , gives essentially the same bounds as Lehmann’s method. In the current paper, we focus on the estimation of eigenfunctions and the two-sided bounds of eigenvalues are assumed to be known.

Error estimates for symmetric elliptic partial differential eigenvalue problems are widely studied in the literature. We refer to classical works Chatelin1983 ; BabOsb:1991 ; Boffi:2010 for the fundamental theories. The majority of the existing literature concerns error estimates valid asymptotically or containing unknown constants; see, e.g., ArmDur2004 ; DarDurPad2012 ; DurGasPad1999 ; giani2018posteriori ; GiaHal2012 ; HuHuaLin2014 ; JiaCheXie2013 ; MehMie2011 ; Yang2010 . Recently, fully computable (containing no unknown constants) and guaranteed (bounding the error from above on all meshes, not only asymptotically) error estimates for eigenvalue problems appeared. Papers CarGal2014 ; CarGed2014 ; Liu2015 ; Liu2020 ; LiuOis2013 ; SebVej2014 ; Vejchodsky2018b ; Vejchodsky2018 concern the eigenvalues. Particularly, as a general framework, the method proposed in Liu2015 ; Liu2020 has been applied to eigenvalue problems of various differential operators, including the Stokes operator Xie2LIU-2018 , the Steklov operator you-xie-liu-2019 , and biharmonic operators related to the quadratic interpolation error constants liu-you:2018 ; LiaoYuLiu2019 .

Concerning eigenfunctions, Davis and Kahan davis1970rotation provide a fundamental result for the estimation of the distance of eigenspaces utilizing the strong residual $R=A\hat{x}-\hat{\lambda}\hat{x}$, where $\{\hat{x},\hat{\lambda}\}$ stands for an approximate eigenpair of a matrix (or operator) $A$. Recent paper Nakatsukasa2020 extends Davis–Kahan’s approach by further orthogonal decomposition of $A$ with respect to the trial eigenspaces.

Similar to Davis–Kahan’s approach, papers CanDusMadStaVoh2017 ; CanDusMadStaVoh2018 ; CanDusMadStaVoh2019 ; CarGed2014 ; HongXieYueZhang2018 consider the weak residual. These methods provide guaranteed, robust, and optimally convergent a posteriori bounds for eigenvalues and corresponding eigenfunctions for finite element approximations assuming an a priori knowledge of bounds on eigenvalues. All these existing approaches introduce dual functions from $H(\mbox{div};\Omega)$ to approximate the gradient of eigenfuctions and need to solve an auxiliary problem to find them. Alternative approach toyonaga2002verified formulates the error estimation problem as a fixed-point problem. Section 6 contains the comparison of the efficiency of our proposed method and the one in CanDusMadStaVoh2019 .

In contrast to the existing literature, we propose two algorithms to estimate error of approximate eigenfunctions. Proposed error estimates have the following features.

• •

Algorithm I in §3 is based on the Rayleigh quotient and the min-max principle that characterizes the eigenvalue problem. This algorithm only utilizes Rayleigh quotients of approximate eigenfunctions and thus the estimator is easy to compute, especially when compared to Algorithm II and the existing literature. A defect of this algorithm is that the estimate for the concerned cluster depends on previous clusters and the width of the concerned cluster.

• •

Algorithm II in §5 utilizes the residual error for variational representation of differential operators. Here, Davis–Kahan’s method, originally proposed for strongly formulated eigenvalue problems, is successfully extended to weakly formulated eigenvalue problems. In case of the Laplace eigenvalue problem, we obtain sharper estimate for the error of approximate eigenfunctions by further utilizing an auxiliary flux $p_{h}\in H(\mbox{div})$.

The estimates by Algorithm II are independent from the cluster index, its width (difference of the largest and smallest eigenvalue in the cluster), and its degree (number of eigenvalues in the cluster). The quality of the estimate depends solely on the residual error and the spectral gap (distance of the cluster from its neighboring clusters). Numerical examples presented below show that this estimate is very accurate, especially for finite element approximations.

The rest of the paper is organized as follows. Section 2 briefly recalls the main properties of the abstract eigenvalue problem. Section 3 presents the a posteriori error bounds for eigenfunctions in both $a$- and $b$-sense. Section 4 derives an error bound for the directed distance measured in the $a$-sense based on the error bound in the $b$-sense. Section 5 utilizes the residual error of the differential operators and extended Davis–Kahan’s theorem to obtain the error estimation. Section 6 provides numerical examples for matrix and Laplace eigenvalue problems. Finally, Section 7 draws the conclusions.

2 Eigenvalue problem for a compact self-adjoint operator

As we mentioned above, weak formulations of elliptic partial differential eigenvalue problems motivate us to consider Hilbert spaces $V$ and $W$ with inner products $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$, respectively, and a compact operator $\gamma:V\to W$ between them. Let $T:W\to V$ be the solution operator such that for given $f\in W$, the element $Tf\in V$ is uniquely determined by the identity

[TABLE]

It is easy to see that $T\gamma:V\to V$ is a compact self-adjoint operator. The eigenvalue problem for the operator $T\gamma$ is equivalent to finding eigenfunctions $u_{i}\in V\setminus\{0\}$ and corresponding eigenvalues $\mu_{i}\geq 0$ such that

[TABLE]

The Hilbert–Schmidt Theorem (see, e.g., (debnath2005introduction, , Theorem 4.10.1)) provides the existence of a countable sequence of (finite) non-zero eigenvalues

[TABLE]

of the operator $T\gamma$, where $d=\operatorname{dim}V$ can be infinite. Note that these eigenvalues are repeated according to their multiplicities and that if $\mbox{dim}(\operatorname{Ker}(T\gamma))>0$ then there are also zero eigenvalues of $T\gamma$.

Eigenvalue problems for linear elliptic partial differential operators are usually formulated for the reciprocals of $\mu_{i}$. Denoting $\lambda_{i}=1/\mu_{i}$ for $\mu_{i}>0$, we clearly have

[TABLE]

namely the eigenvalue problem (1).

The eigenfunctions $u_{i}\in V$ are assumed to be normalized such that

[TABLE]

where $\delta_{ij}$ stands for the Kronecker delta. With this normalization, the Hilbert–Schmidt Theorem also implies that any $v\in V$ can be expressed as

[TABLE]

Noticing that $b(\gamma v_{0},\gamma v)=a(T\gamma v_{0},v)=0$ for all $v\in V$, we have

[TABLE]

Remark 1

The Laplace eigenvalue problem in a domain $\Omega$ with homogeneous Dirichlet boundary conditions fits the above setting for

[TABLE]

The Steklov eigenvalue problem fits this settings for

[TABLE]

The generalized matrix eigenvalue problem $Au_{i}=\lambda_{i}Bu_{i}$ for a symmetric positive definite matrix $A\in\mathbb{R}^{d\times d}$ and a symmetric positive semidefinite matrix $B\in\mathbb{R}^{d\times d}$ fits the above setting for $d<\infty$, $V=\mathbb{R}^{d}$, $W=\operatorname{range}B$, $\gamma=B^{1/2}$, $a(u,v)=v^{T}Au$, and $b(u,v)=v^{T}u$.

In order to formulate the bound on eigenfunctions, a notation for clusters of eigenvalues has to be introduced. Let us focus on the leading $K$ clusters. Let $n_{k}$ and $N_{k}$ stand for indices of the first and the last eigenvalue in the $k$-th cluster, $k=1,2,\dots,K$, respectively. In particular, $n_{1}=1$, $n_{k+1}=N_{k}+1$, and the $k$-th cluster is formed of $N_{k}-n_{k}+1$ eigenvalues $\lambda_{n_{k}}$, $\lambda_{n_{k}+1}$, …, $\lambda_{N_{k}}$; see Figure 1. Notice that the eigenvalues in a cluster do not necessarily equal to each other. To simplify the notation, we set $n=n_{K}$ and $N=N_{K}$.

Each cluster is associated with the space $E_{k}=\operatorname{span}\{u_{n_{k}},u_{n_{k}+1},\dots,u_{N_{k}}\}$ of exact eigenfunctions. Similarly, arbitrary linearly independent approximations $\hat{u}_{i}\in V$ of exact eigenfunctions $u_{i}$, $i=1,2,\dots,N_{K}$, form the corresponding approximate spaces $\widehat{E}_{k}=\operatorname{span}\{\hat{u}_{n_{k}},\hat{u}_{n_{k}+1},\dots,\hat{u}_{N_{k}}\}$. Thus $\mbox{dim}(\widehat{E}_{k})=N_{k}-n_{k}+1$. Spaces $\widehat{E}_{k}$, $k=1,2,\dots,K$, of approximate eigenfunctions need not be orthogonal to each other.

Remark 2

Given an eigenvalue problem, the a priori information about the distribution of eigenvalues is usually unknown. Therefore, in practical problems, the partition of the spectrum into clusters is determined numerically using two-sided bounds on eigenvalues.

If the enclosing intervals of eigenvalues overlap, it is natural to define a cluster for such eigenvalues. Consequently, the internal distribution of eigenvalues inside the cluster is regarded as uncertain and only the lower and upper bound on the whole cluster will be used for analysis.

3 Algorithm I: Estimation for eigenfunctions based on the Rayleigh quotients

The goal of this section is to derive an estimate of the directed distance between spaces $E_{K}$ and $\widehat{E}_{K}$ of exact and approximate eigenfunctions for the $K$-th cluster. The directed distance between two subspaces $E$ and $\widehat{E}$ of $V$ is defined both in the $a$- and $b$-sense as follows; see e.g., (Meyer:2000, , §5.15),

[TABLE]

The directed distance is not symmetric, in general, but if the two subspaces are of the same finite dimension, which is the case in this paper, then the directed distance is symmetric and coincides with the gap between spaces. The following characterization of the directed distance will be used in the proof of the main theorem; see also (Meyer:2000, , §5.15).

Note that in CanDusMadStaVoh2019 , the Hilbert–Schmidt norm for the cluster residual error is utilized to estimate the error of approximate eigenfunction.

Lemma 1

Let $E$ and $\widehat{E}$ be the two finite dimensional subspaces of the Hilbert space $V$. Let $\widehat{P}:V\to\widehat{E}$ be the orthogonal projector with respect to $a(\cdot,\cdot)$. Then, the directed distance defined in (6) satisfies

[TABLE]

Proof

The first statement in (7) follows immediately from (6), because the orthogonal projection $\widehat{P}v$ is the closest element to $v$ in the subspace $\widehat{E}$. For the second statement, we apply the same idea to the orthogonal projector $\widehat{\Pi}:\gamma V\to\gamma\widehat{E}$ with respect to the scalar product $b(\cdot,\cdot)$. This projector satisfies $\widehat{\Pi}\gamma=\gamma\widehat{P}$ and, hence, the second statement holds true. Identities (8) follow from (7) due to the equality $\|\widehat{P}v\|_{a}^{2}+\|v-\widehat{P}v\|_{a}^{2}=\|v\|_{a}^{2}$ and a similar equality for $\widehat{\Pi}$.

Remark 3

For one dimensional spaces $E=\operatorname{span}\{u\}$ and $\widehat{E}=\operatorname{span}\{\hat{u}\}$, we have $\delta_{a}(E,\widehat{E})=\sin\alpha$, where $\alpha$ is the angle between $u$ and $\hat{u}$. Consequently, the difference of $u$ and $\hat{u}$ in the $a$-norm can be expressed by the cosine theorem as

[TABLE]

Hence, an upper bound on $\delta_{a}(E,\widehat{E})$ yields immediately an upper bound on the error $\|u-\hat{u}\|_{a}$. Moreover, if $\|u\|_{a}=\|\hat{u}\|_{a}=1$ then $\|u-\hat{u}\|_{a}^{2}=\delta_{a}^{2}(E,\widehat{E})+O(\delta_{a}^{4}(E,\widehat{E}))$ and the difference between the directed distance $\delta_{a}(E,\widehat{E})$ and the $a$-norm of the error is negligible for small $\delta_{a}(E,\widehat{E})$. The same conclusions clearly hold also for the norm induced by $b$.

The main result of this section is presented in Theorem 3.1 below. In order to formulate it, we introduce measures of non-orthogonality between finite dimensional subspaces $E$ and $E^{\prime}$ of $V$ both in the $a$- and $b$-sense as

[TABLE]

Both these measures of non-orthogonality can be easily computed or estimated by using the following lemma. To simplify the exposition, we present the result for $\hat{\varepsilon}_{a}$ only, because the statement for $\hat{\varepsilon}_{b}$ is completely analogous.

Lemma 2

Let $v_{1},v_{2},\dots,v_{m}$ and $v^{\prime}_{1},v^{\prime}_{2},\dots,v^{\prime}_{m^{\prime}}$ form bases of subspaces $E$ and $E^{\prime}$ of $V$, respectively. Define matrices $F$, $G$, $H$ as follows,

[TABLE]

Then

[TABLE]

where $\lambda_{max}(A,B)$ denotes the maximum eigenvalue of eigen-problem $Ax=\lambda Bx$.

Further, suppose $\|F^{T}F\|_{2}\leq\eta_{F}$, $\|I-G\|_{2}\leq\eta_{G}$, $\|I-H\|_{2}\leq\eta_{H}$. If $\eta_{G},\eta_{H}<1$, then

[TABLE]

Proof

Expand $v\in E$ and $v^{\prime}\in E^{\prime}$ as $v=\sum_{i=1}^{m}c_{i}v_{i}$ and $v^{\prime}=\sum_{j=1}^{m^{\prime}}c^{\prime}_{j}v^{\prime}_{j}$ and consider vectors $\boldsymbol{c}\in\mathbb{R}^{m}$ and $\boldsymbol{c^{\prime}}\in\mathbb{R}^{m^{\prime}}$ of coefficients $c_{i}$ and $c^{\prime}_{j}$, respectively. Then

[TABLE]

Thus, definition (9) gives

[TABLE]

where $\boldsymbol{\tilde{c}^{\prime}}=L^{T}\boldsymbol{c^{\prime}}$, $H=LL^{T}$ is the Cholesky decomposition of matrix $H$, and $|\cdot|$ stands for the Euclidean norm. Consequently,

[TABLE]

Expression $\hat{\varepsilon}_{a}^{2}(E,E^{\prime})=\lambda_{max}(F^{T}G^{-1}F,H)$ can be proved analogously.

To prove the upper bound on $\hat{\varepsilon}_{a}$, we use Cholesky decomposition $G=QQ^{T}$. Noticing that $\|A\|_{2}=\|A^{T}\|_{2}=\sqrt{\|A^{T}A\|_{2}}$ holds for a general matrix $A$, we have

[TABLE]

Finally, we estimate $\|G^{-1}\|_{2}$ and $\|H^{-1}\|_{2}$. If $\eta_{G}<1$, then

[TABLE]

With the same argument for $H^{-1}$, we easily draw the conclusion.

Remark 4

This lemma is used for spaces $\widehat{E}_{k}$ and $\widehat{E}_{k^{\prime}}$ of approximate eigenfunctions with $k\neq k^{\prime}$, therefore matrices $F$, $G$, and $H$ are available and $\lambda_{max}(F^{T}G^{-1}F,H)$ as well as $\lambda_{max}(FH^{-1}F^{T},G)$ can be computed. Alternatively, guaranteed estimates $\eta_{F}$, $\eta_{H}$, and $\eta_{G}$ can be obtained by the Gershgorin circle theorem. These estimates are expected to be good for $\hat{\varepsilon}_{a}(\widehat{E}_{k},\widehat{E}_{k^{\prime}})$, because if the approximate eigenfunctions in $\widehat{E}_{k}$ and the ones in $\widehat{E}_{k^{\prime}}$ are appropriately orthonormalized, then $F^{T}F\approx 0$, $G\approx I_{m}$, and $H\approx I_{m^{\prime}}$.

The following theorem provides the desired estimates of the directed distances $\delta_{a}(E_{K},\widehat{E}_{K})$ and $\delta_{b}(E_{K},\widehat{E}_{K})$ defined in (6).

Theorem 3.1

Let the above specified partition of the spectrum into $K$ clusters be arbitrary. Let $\hat{u}_{i}\in V$ for $i=1,2,\dots,N$ with $N<d=\operatorname{dim}V$ be such that $\operatorname{dim}\widehat{E}_{k}=N_{k}-n_{k}+1$ for all $k=1,2,\dots,K$. Let $\lambda_{n}<\rho\leq\lambda_{N+1}$. Then

[TABLE]

where

[TABLE]

Proof

Following (5), the proof is based on the equality

[TABLE]

for an arbitrary and fixed $\hat{u}\in\widehat{E}_{K}$, where $\hat{u}_{0}$ is the component of $\hat{u}$ in $\operatorname{Ker}(T\gamma)$ and $d$ can be infinity. Recalling that $n=n_{K}$ and $N=N_{K}$, using definition

[TABLE]

and inequalities $\|\hat{u}_{0}\|_{a}^{2}\geq 0$ and $\rho\leq\lambda_{N+1}$, we easily derive bound

[TABLE]

Introducing orthogonal projectors $P_{k}:V\rightarrow E_{k}$, $k=1,2,\dots,K$, with respect to the inner product $a(\cdot,\cdot)$, recalling their properties

[TABLE]

and using $\lambda_{n}\leq\lambda_{i}$ for $i=n,\dots,N$, we immediately obtain estimates

[TABLE]

Similarly, the fact that $\lambda_{n_{k}}\leq\lambda_{i}$ for $i=n_{k},\dots,N_{k}$ implies the following estimate on $\vartheta(\hat{u})$:

[TABLE]

Now, we bound $\|P_{k}\hat{u}\|_{a}$ and $\|\gamma P_{k}\hat{u}\|_{b}$ for $k=1,2,\dots,K-1$. Introducing $z_{k}=P_{k}\hat{u}/\|P_{k}\hat{u}\|_{a}\in E_{k}$ and the $a$-orthogonal projector $\widehat{P}_{k}:V\rightarrow\widehat{E}_{k}$, definition (9) and relation (7) imply

[TABLE]

Since $P_{k}\hat{u}=a(\hat{u},z_{k})z_{k}$ and $\|\widehat{P}_{k}z_{k}\|_{a}\leq\|z_{k}\|_{a}=1$, these estimates provide the bound

[TABLE]

Analogous steps yield the bound

[TABLE]

Consequently, estimates (16) provide

[TABLE]

The desired lower bounds on $\|P_{K}\hat{u}\|_{a}$ and $\|\gamma P_{K}\hat{u}\|_{b}$ follow from (15) and (17):

[TABLE]

The final step is to apply these two estimates in expressions (8) for the distances $\delta_{a}^{2}(\widehat{E}_{K},E_{K})$ and $\delta_{b}^{2}(\widehat{E}_{K},E_{K})$. Indeed,

[TABLE]

and the statement (10) follows by elementary manipulations utilizing the definition of $\hat{\lambda}^{(K)}_{N}$. Statement (11) follows analogously. Note that $\delta_{a}(E_{K},\widehat{E}_{K})=\delta_{a}(\widehat{E}_{K},E_{K})$ and $\delta_{b}(E_{K},\widehat{E}_{K})=\delta_{b}(\widehat{E}_{K},E_{K})$, because $\operatorname{dim}E_{K}=\operatorname{dim}\widehat{E}_{K}=N-n+1$.

Remark 5

Theorem 3.1 is a direct and nontrivial generalization of (BirBooSwaWen:1966, , Corollary 1). Indeed, if all eigenvalues $\lambda_{i}$ are simple and well separated (forming clusters of size one), $E_{i}=\operatorname{span}\{u_{i}\}$ and $\widehat{E}_{i}=\operatorname{span}\{\hat{u}_{i}\}$ stand for the corresponding one-dimensional exact and approximate eigenspaces, $\hat{\lambda}_{i}=\|\hat{u}_{i}\|_{a}^{2}/\|\hat{u}_{i}\|_{b}^{2}$, and the corresponding approximate eigenfunctions $\hat{u}_{i}$, $i=1,2,\dots,k$, are mutually orthogonal, then bounds (10) and (11) become

[TABLE]

and

[TABLE]

where the estimate for $\delta_{b}$ coincides with the statement in (BirBooSwaWen:1966, , Corollary 1).

Remark 6

In this remark, we show the necessity of terms $\rho-\lambda_{n}$, $\hat{\lambda}_{N}^{(K)}-\lambda_{n}$ and $\theta_{a}^{(K)}$ or $\theta_{b}^{(K)}$ appearing in bounds (10) and (11).

• (i)

The difference $\rho-\lambda_{n}$ is determined by the spectral gap between the last cluster and the following eigenvalues. To have a sharp bound of the approximation error, the gap should not be too small.

• (ii)

Quantity $\hat{\lambda}_{N}^{(K)}-\lambda_{n}$ corresponds to the width of the last cluster. Note that in case of non-degenerated cluster, i.e., $\lambda_{N}>\lambda_{n}$, this difference causes bounds (10) and (11) not to converge to zero even if $\hat{\lambda}_{N}^{(K)}\to\lambda_{N}$. In fact, under the assumption that only the two-side bound of the whole cluster is available, the distribution of exact eigenvalues is uncertain and so are the corresponding exact eigenvectors. To illustrate this uncertainty in eigenvectors, let us consider a matrix eigenvalue problem $Ax=\lambda x$ with $A$ selected as

[TABLE]

Eigenvalues of $A$ are given by

[TABLE]

and the corresponding the eigenvectors are

[TABLE]

Let $a(u,v):=u^{T}Av$, $b(u,v):=u^{T}v$ for column vector $u,v\in\mathbb{R}^{4}$. Considering the approximate eigenvectors $\widehat{u}_{i}=e_{i}$ for $i=1,2,3$ and spaces $E_{1}=\operatorname{span}\{u_{1},u_{2},u_{3}\}$, $\widehat{E}_{1}=\operatorname{span}\{\widehat{u}_{1},\widehat{u}_{2},\widehat{u}_{3}\}$, we compute the distance $\delta_{b}^{2}(E_{1},\widehat{E}_{1})$ exactly and apply bound (11) with $\rho=\lambda_{4}$ to obtain

[TABLE]

Clearly, the quality of approximate vectors measured by the distance $\delta_{b}$ depends on the position of $\lambda_{2}$ (determined by $\alpha$). For such a system with uncertain $\alpha$, one cannot expect a sharper estimation for $\delta_{b}$.

If further information about eigenvalues inside a tight cluster is known, we can (theoretically) split the cluster into smaller clusters consisting of a single or a multiple eigenvalue and for these clusters bounds (10) and (11) do converge. Bounds (10) and (11) are naturally computed iteratively starting from the first cluster. Accuracy of this procedure is illustrated on numerical examples in Section 6.

• (iii)

Values of $\theta_{a}^{(K)}$ and $\theta_{b}^{(K)}$ measure errors in all previous clusters. Notice that quantity $\theta_{a}^{(K)}$ (and similarly $\theta_{b}^{(K)}$) depends on $\delta_{a}(E_{k},\widehat{E}_{k})$ and $\hat{\varepsilon}_{a}(\widehat{E}_{k},\widehat{E}_{K})$ for $k=1,2,\dots,K-1$. The distance $\delta_{a}(E_{k},\widehat{E}_{k})$ accounts for errors in spaces of eigenfunctions of previous clusters and $\hat{\varepsilon}_{a}(\widehat{E}_{k},\widehat{E}_{K})$ for possible non-orthogonality of approximate eigenfunctions.

Note that the dependence of bounds (10) and (11) on errors in previous clusters is necessary, because these bounds utilize only the information about eigenvalue gaps and Rayleigh quotients. Indeed, consider the one-dimensional Dirichlet eigenvalue problem:

[TABLE]

and approximate eigenfunctions

[TABLE]

Clearly, $\hat{u}_{2}$ is a bad approximation of the exact eigenfunction $u_{2}(x)=\sin(2x)$. Taking naturally $E_{2}=\operatorname{span}\{u_{2}\}$, $\widehat{E}_{2}=\operatorname{span}\{\hat{u}_{2}\}$, distances $\delta_{a}^{2}(E_{2},\widehat{E}_{2})$ and $\delta_{b}^{2}(E_{2},\widehat{E}_{2})$ are relatively large. On the other hand, the corresponding Rayleigh quotient gives the exact eigenvalue, namely $R(\hat{u}_{2})=\hat{\lambda}_{2}^{(2)}=\lambda_{2}=4$. The difference $\hat{\lambda}_{2}^{(2)}-\lambda_{2}$ in (10) and (11) is zero and values of these upper bounds are determined by terms $\theta_{a}^{(2)}$ and $\theta_{b}^{(2)}$ only, reflecting the non-orthogonality of $\hat{u}_{1}$ and $\hat{u}_{2}$.

Terms $\theta_{a}^{(2)}$ and $\theta_{b}^{(2)}$ could be avoided if, for example, guaranteed estimates of the residual are employed. This approach will be introduced in §5.1 as Algorithm II.

The following theorem addresses the question of efficiency of bounds presented in Theorem 3.1. To formulate it, we denote the right-hand sides of estimates (10) and (11) by $\overline{\delta}_{a}(E_{K},\widehat{E}_{K})$ and $\overline{\delta}_{b}(E_{K},\widehat{E}_{K})$, respectively.

Theorem 3.2

If eigenspaces $\widehat{E}_{k}$ for $k=1,2,\dots,K$ are mutually orthogonal then there exists a generic constant $C>0$ determined by the leading ${N_{K}}$ exact and approximate eigenvalues and independent of the width of clusters such that

[TABLE]

Proof

Take $\hat{u}\in\widehat{E}_{K}$ such that $\|\hat{u}\|^{2}_{a}=\hat{\lambda}^{(K)}_{N}\|\gamma\hat{u}\|_{b}^{2}$. Let $u=P_{K}\hat{u}\in E_{K}$ be the best approximation to $\hat{u}$ under the norm $\|\cdot\|_{a}$. By Lemma 1, $\|\hat{u}-u\|_{a}\leq\delta_{a}(E_{K},\widehat{E}_{K})$. Suppose $u=\sum_{i=n}^{N}c_{i}u_{i}$ and take $\tilde{\lambda}\in[\lambda_{n},\lambda_{N}]$ such that $\tilde{\lambda}\sum_{i=n}^{N}c_{i}^{2}=\sum_{i=n}^{N}c_{i}^{2}\lambda_{i}$. By expanding the squared terms $\|\hat{u}-u\|_{a}^{2}$ and $\|\gamma(\hat{u}-u)\|_{b}^{2}$, we have

[TABLE]

Noticing that $\|\gamma(\hat{u}-u)\|_{b}^{2}\leq 1/\lambda_{1}\|\hat{u}-u\|_{a}^{2}$ and $\lambda_{n}\leq\tilde{\lambda}\leq\lambda_{N}$, the equation above tells

[TABLE]

Therefore, for the term $(\hat{\lambda}_{N}^{(K)}-\lambda_{n})$ in (10) and (11), we have

[TABLE]

The mutual orthogonality of $\widehat{E}_{k}$ for $k=1,2,\dots,K$ implies $\hat{\epsilon}_{a}(E_{k},E_{K})=\hat{\epsilon}_{b}(E_{k},E_{K})=0$. Consequently, terms $\theta_{a}^{(K)}$ and $\theta_{b}^{(K)}$ can be estimated recursively, leading to the inequality (20).

As pointed out in (ii) of Remark 6, the term $\lambda_{N_{k}}-\lambda_{n_{k}}$ is necessary for clusters of non-zero width. If all clusters are of zero size and errors in previous clusters are comparable to the error in the cluster of interest, i.e.,

[TABLE]

then inequality (20) immediately yields the following efficiency result

[TABLE]

Remark 7

In Section 6, we apply the general bounds (10) and (11) in the particular context of the finite element method. This method, specifically, computes the best approximation ($a$-orthogonal projection) of the exact solution of the underlying boundary value problem in the finite element space. As the numerical results in Section 6 reveal, the bound (10) is optimal in the sense that it has the same rate of convergence as $\delta_{a}(\widehat{E}_{K},E_{K})$. However, the bound (11) has a lower rate of convergence than $\delta_{b}(\widehat{E}_{K},E_{K})$.

The Aubin–Nitsche technique together with explicit a priori error estimates LiuOis2013 will provide fully computable and optimal bounds under $b$-norm. However, since the application of the Aubin–Nitsche technique utilizes special properties about FEM approximation, such an approach is beyond the problem setting of this paper and will be addressed in a subsequent publication; partial results of this approach can be found in liu-Vejchodsky-arxiv-v1 .

4 Sharp $a$-norm estimates based on the $b$-norm bounds

This section provides an estimate of the distance $\delta_{a}(E,\widehat{E})$ by the distance $\delta_{b}(E,\widehat{E})$ between the space of exact eigenfunctions $E$ and the space of approximate eigenfunctions $\widehat{E}$. The idea is motivated by the following well known formula (see e.g. (Boffi:2010, , page 55))

[TABLE]

for the exact eigenpair $\lambda_{i}>0$, $u_{i}\in V$, arbitrary approximate eigenfunction $\hat{u}_{i}\in V$ and approximate eigenvalue $\hat{\lambda}_{i}=\|\hat{u}_{i}\|_{a}^{2}/\|\gamma\hat{u}_{i}\|_{b}^{2}$. In the context of the finite element method for the Laplace eigenvalue problem, this identity essentially says that the error $\|u_{i}-\hat{u}_{i}\|_{a}$ in the energy norm is dominated by the error of the approximate eigenvalue itself, because the error $\|\gamma(u_{i}-\hat{u}_{i})\|_{b}$ in the $L^{2}$ norm has a higher order of convergence.

The following estimate is theoretically independent of the partition of eigenvalues into clusters, but its natural usage is to bound $\delta_{a}(E_{k},\widehat{E}_{k})$ by $\delta_{b}(E_{k},\widehat{E}_{k})$, where $k$ is the index of a cluster as it is introduced at the end of Section 2.

Theorem 4.1

Let $u_{n},\dots,u_{N}$ be the exact eigenfunctions of (1) and $0<n\leq N$ the corresponding indices. Let $\hat{u}_{n},\dots,\hat{u}_{N}\in V$ be linearly independent. Let $E=\operatorname{span}\{u_{n},\dots,u_{N}\}$ and $\widehat{E}=\operatorname{span}\{\hat{u}_{n},\dots,\hat{u}_{N}\}$. Then

[TABLE]

where $\lambda_{n}$ and $\lambda_{N}$ are exact eigenvalues corresponding to $u_{n}$ and $u_{N}$ and

[TABLE]

Proof

Consider the linear mapping $\tau:E\to E$ defined by

[TABLE]

Since $\lambda_{i}>0$ for all $i=n,\dots,N$, $\tau$ is a bijection. Given arbitrary $u\in E$ and $\hat{u}\in\widehat{E}$, we clearly have

[TABLE]

This enables us to estimate the distance (6) between $E$ and $\widehat{E}$ as follows

[TABLE]

where the last inequality follows from the fact that

[TABLE]

Since $\tau$ is a bijection, it is easy to show that

[TABLE]

This equality together with bounds $\|u\|_{a}^{2}\leq\lambda_{N}\|\gamma u\|_{b}^{2}$ for all $u\in E$ and $\|\hat{u}\|_{a}^{2}\leq\hat{\lambda}_{N}\|\gamma\hat{u}\|_{b}^{2}$ for all $\hat{u}\in\widehat{E}$ imply

[TABLE]

where we note that $\max_{\hat{u}\in\widehat{E},\ \|\gamma\hat{u}\|_{b}=1}b(\gamma u,\gamma\hat{u})$ is non-negative and the last equality follows from (8) using characterization

[TABLE]

The proof is finished by substituting (25) to (24).

Remark 8

The estimate (23) will provide a sharper bound on $\delta_{a}(E,\widehat{E})$ if the eigenvector approximation has a higher convergence rate under $\|\cdot\|_{b}$ norm compared with $\|\cdot\|_{a}$ norm, which usually happens for the finite element approximation of elliptic eigenvalue problems.

For a multiple eigenvalue, we have $\lambda_{n}=\lambda_{N}$ and the right-hand side of (23) is close to

[TABLE]

where we use the approximation $(1-\delta_{b}^{2})^{1/2}\approx 1-\delta_{b}^{2}/2$. This shows that for example in the context of the finite element method for the Laplace eigenvalue problem, the bound (23) has the optimal rate of convergence. Numerical experiments indicate that the bound (23) combined with the sub-optimal estimate (11) can still provide a sharper bound on $\delta_{a}$ than (10); see Section 6.

Define the following quantity to measure the distance between $E$ and $\widehat{E}$:

[TABLE]

This quantity can be easily bounded by $\delta_{a}$ as follows

[TABLE]

The following lemma relates $\tilde{\delta}$ and $\delta_{b}$.

Lemma 3

By using $\delta_{b}(E,\widehat{E})$, we have the following estimate for $\tilde{\delta}(E,\widehat{E})$,

[TABLE]

Proof

We utilize the same mapping $\tau:E\to E$ as in the proof of Theorem 4.1 to obtain the estimate for $\tilde{\delta}$. Note that for any $u\in E$, $\|\gamma u\|_{b}=1$, we have $\|\gamma\tau(u)\|_{b}\geq\lambda_{n}$ and consequently

[TABLE]

Remark 9

If approximate eigenvalues converge to the exact ones, i.e. if $|\lambda_{k,h}-\lambda_{k}|\to 0$, then the right-hand sides of estimates (23) and (28) do not converge to zero, due to the existence of the term corresponding to the cluster width $|\lambda_{N}-\lambda_{n}|$. The necessity of this term is confirmed by the discussion of the eigenvalue problem for matrix (18), where $\delta_{a}\to\sqrt{2\beta+\beta^{2}}$ and $\delta_{b}\to 0$ when $\alpha\to 0$. Note that $\beta$ determines the cluster width for the problem in (18). Despite the inefficiency caused by the cluster width, the advantage of utilizing (28) is that for clusters with narrow width, an easy-to-obtain sharp bound on $\delta_{b}$ will greatly improve the precision for $\hat{\delta}$ and thus $\delta_{a}$.

5 Algorithm II: Estimation based on the residual error

Estimates in Theorem 3.1 of §3 only utilize the information about eigenvalues distribution and Rayleigh quotients for approximate eigenfunctions. As discussed in Remark 6, such estimates have the advantage to provide guaranteed error bound for problems with uncertainty information about the objective eigenvalue problem. As a result, the bounds in (10) and (11) depend on the quality of previous clusters and the width of the cluster of interest.

In this section, we utilize the residual error estimation of approximate eigenfunctions to obtain sharper bounds of the distances between the exact and the approximate eigenspaces. For this purpose, we introduce the approach based on the Davis–Kahan’s method davis1970rotation . Note that the original Davis–Kahan’s method only deals with strongly formulated differential operators and we extend it to weakly formulated eigenvalue problems. It is worth to point out that, compared to very cheap estimates given in Theorem 3.1, the approach utilizing the residual will require more computational effort to compute or bound its suitable norm. The complexity of this effort varies and depends on the problem.

The Davis–Kahan’s $\sin\theta$ theorem in davis1970rotation bounds the error of an approximate eigenspace by a norm of the residual and a spectral gap. The bound requires to apply the operator to an approximate eigenfunction to estimate the norm of the residual. This can be easily done for strongly defined operators and sufficiently smooth approximations. However, for weakly defined operators and non-smooth approximations the original Davis–Kahan’s $\sin\theta$ theorem cannot be directly used. Therefore, we extend the Davis–Kahan’s $\sin\theta$ theorem to weakly defined eigenvalue problems in the form (3).

For the reader’s convenience, we quote Davis–Kahan’s result davis1970rotation using the notation of this paper.

Theorem 5.1 (Davis-Kahan’s $\sin\theta$ theorem)

Let $U$ be a Hilbert space with inner product $\langle\cdot,\cdot\rangle$. Let $\|\cdot\|$ be a unitary invariant norm of $U$. Let $A:U\to U$ be a self-adjoint operator and $E$ be the subspace spanned by the eigenvectors of $A$ in the cluster of interest. That is,

[TABLE]

Let $A_{h}$ be a perturbation of the operator $A$, which has an invariant subspace $\widehat{E}$ such that $\widehat{E}=\operatorname{span}\{\hat{u}_{n},\dots,\hat{u}_{N}\}$ and

[TABLE]

Let $P:U\to E$ and $\widehat{P}:U\to\widehat{E}$ be the two projection operators with respect to $\langle\cdot,\cdot\rangle$. Suppose the eigenvalues of $A$ corresponding to $E^{\perp}$ are excluded from $(\hat{\lambda}_{n}-\sigma,\hat{\lambda}_{N}+\sigma)$. Then

[TABLE]

Since the proof in the original paper is divided into several parts for more general discussion, here we provide a reformulation of the proof in davis1970rotation in a concise way.

Proof

Define $c=(\hat{\lambda}_{n}+\hat{\lambda}_{N})/2$, $r=(\hat{\lambda}_{N}-\hat{\lambda}_{n})/2$. Since $AP=PA$, we have

[TABLE]

Since the eigenvalues of $A-cI$ over $E^{\perp}$ are excluded from $(-r-\sigma,r+\sigma)$, we have

[TABLE]

Also, since the eigenvalue of $A_{h}-cI$ over $\widehat{E}$ belongs to $[-r,r]$, we have

[TABLE]

Thus, we draw the conclusion from (29), (30), (31).

Remark 10

If the operator $A$ is strongly defined as a differential operator, then for any given smooth test function $u_{h}$ from $U$, we can easily evaluate $Au_{h}$. However, if the problem is weakly formulated, one cannot evaluate $Au_{h}$ directly. In the succeeding subsection, we will introduce the extension of Davis–Kahan’s result to weakly formulated eigenvalue problems.

5.1 Extension of Davis–Kahan’s $\sin\theta$ theorem to weakly formulated problems

Let us recall the weakly formulated eigenvalue problem (3) and the composition $T\gamma:V\to V$ of the compact operator $\gamma$ and the solution operator $T$. We set $U=\mbox{Ker}(T\gamma)^{\perp}=\mbox{Ker}(\gamma)^{\perp}(\subset V)$ and notice that the operator $T\gamma$ is invertible on $U$. Thus, the operator $A=\left(T\gamma|_{U}\right)^{-1}:U\to U$ is well defined and satisfies

[TABLE]

Since $T\gamma$ is self-adjoint, the operator $A$ is self-adjoint as well and we may apply the Davis–Kahan’s method for it.

As above, we consider an eigenvalue cluster of interest and denote the space of exact eigenfunctions corresponding to eigenvalues in the cluster by $E=\operatorname{span}\{u_{n},\dots,u_{N}\}$. These eigenfunctions are approximated by the one from the finite dimensional space $\widehat{E}=\operatorname{span}\{\hat{u}_{n},\dots,\hat{u}_{N}\}$. Note that $\widehat{E}$ need not be necessarily a subspace of $U$. This setting enables us to define an approximate operator $A_{h}:\widehat{E}\to\widehat{E}$ by the identity

[TABLE]

Thus, the discrete eigenvalue problem to find $\hat{\lambda}_{i}$ and $\hat{u}_{i}\in\widehat{E}\setminus\{0\}$ such that

[TABLE]

is equivalent to the eigenvalue problem

[TABLE]

Recall that $P$ and $\widehat{P}$ are projections that map $V$ to $E$ and $\widehat{E}$ with respect to the inner product $a(\cdot,\cdot)$, respectively. To measure the residual error of $(A-A_{h})$, let us introduce the quantity $\epsilon$ by

[TABLE]

Here, the norm $\|\cdot\|_{-1}$ is defined as

[TABLE]

Note that Definitions 3.4, 3.5 from (CanDusMadStaVoh2019, , §3.2) and Definition 2.1 from CanDusMadStaVoh2017 use a similar quantity to $\epsilon$ to derive the estimator.

Lemma 4

Given a normed space $B$ with a norm $\|\cdot\|_{B}$, let $L:B\to B$ be a linear mapping such that $\|L\|_{B}\leq r$. Let $F:B\to\mathbb{R}$ be an absolutely homogeneous functional (i.e., $F(\alpha v)=|\alpha|F(v)$ for all $\alpha\in\mathbb{R}$ and $v\in B$). Then

[TABLE]

Proof

The proof is easily done by noting the following inequality,

[TABLE]

where the first inequality follows from the absolute homogeneity of $F$ and the second one from the fact that $\operatorname{Range}(L)\subset B$ and from the definition of the operator norm. The final equality uses again the absolute homogeneity of $F$.

We formulate the extension of the Davis–Kahan’s $\sin\theta$ theorem below in Theorem 5.2.

Theorem 5.2

Suppose the eigenvalues of $A$ corresponding to $E^{\perp}$ are excluded from $(\hat{\lambda}_{n}-\sigma,\hat{\lambda}_{N}+\sigma)$ i.e., assume that there exists $\sigma>0$ such that

[TABLE]

Then

[TABLE]

Proof

Define $c:=(\hat{\lambda}_{n}+\hat{\lambda}_{N})/2$, $r:=(\hat{\lambda}_{N}-\hat{\lambda}_{n})/2$. The proof is based on the triangle inequality

[TABLE]

where $u\in U$. Note for any $\widehat{u}$ in $\widehat{E}$, $(A_{h}-cI)\widehat{u}\in\widehat{E}$ and $\|\gamma(A_{h}-cI)\widehat{u}\|_{b}\leq r\|\gamma\widehat{u}\|_{b}$. Lemma 4 leads to the following inequality,

[TABLE]

Note that for any $v\in E^{\perp}$, $\|\gamma(A-cI)v\|_{b}\geq(r+\sigma)\|\gamma v\|_{b}$. For any $u\in E^{\perp}$, by taking $v:=(A-cI)|_{E^{\perp}}^{-1}u\in E^{\perp}$, we have

[TABLE]

That is,

[TABLE]

Noticing that $(I-P)\widehat{P}=(A-cI)|_{E^{\perp}}^{-1}(A-cI)(I-P)\widehat{P}$, Lemma 4 along with (37) implies

[TABLE]

Inequalities (35), (36), and (38) lead to the following result.

[TABLE]

Next, we apply the above estimation to bound $\delta_{b}(E,\widehat{E})$. Note that

[TABLE]

Therefore, we obtain the following estimation for $\delta_{b}(E,\widehat{E})$:

[TABLE]

Applying relation (28) between $\delta_{b}$ and $\tilde{\delta}$ to the statement of Theorem 5.2, we obtain the following bound on $\delta_{b}(E,\widehat{E})$.

[TABLE]

Inequality (39) can be reformulated as a quadratic inequality for $t:=\delta_{b}^{2}$,

[TABLE]

with $\alpha_{1}=2\epsilon^{2}\lambda_{n}$ and $\alpha_{2}=\epsilon^{2}\cdot(\lambda_{N}+\hat{\lambda}_{N})$. By solving this quadratic inequality, we derive the following theorem.

Theorem 5.3

Let the two solution of $g(t)=0$ be $t_{1},t_{2}$. In case $\alpha_{1}<\alpha_{2}$ and $0<t_{1}\leq\alpha_{2}<t_{2}$, then we have $\delta_{b}^{2}\leq t_{1}$. That is,

[TABLE]

Estimates for $\delta_{a}$ and $\tilde{\delta}$ are available through the relation (23) and (28).

Remark 11

Note that for small value of $\delta_{b}$, we have $\sqrt{1-\delta_{b}^{2}}\approx(1-\delta_{b}^{2}/2)$. Hence

[TABLE]

The quality of the estimator for $\delta_{b}$ depends on the residual error, the width of the specified cluster and the quality of the eigenvalue estimation. As we will see in next sub-section, the quantity $\epsilon$ does not depend on the width of the cluster. In solving practical problems, if the width of the specified cluster is large, one can try to split the cluster to reduce the cluster width.

In the case of solving the Laplace eigenvalue problem by the conforming linear finite element method, for which it is expected that $O(|\hat{\lambda}_{N}-\lambda_{N}|)=O(h^{2})$, $\epsilon=O(h)$ and $\delta_{b}=O(h^{2})$. Our proposed estimation coincides with this fact, if the cluster width is zero or small enough compared to $O(|\hat{\lambda}_{N}-\lambda_{N}|)$.

Remark 12

The following identity for any $u_{h}\in\widehat{E}$ and $v\in U$, will be useful for computing a tight bound on $\epsilon$:

[TABLE]

5.2 Weakly formulated residual error estimation

In this subsection, we consider a concrete eigenvalue problem of the Laplacian and utilize the Prager–Synge technique PraSyn:1947 to compute a tight bound on the residual error term $\epsilon$. The technique relies on a flux reconstruction, which we explain for the case of the Laplace eigenvalue problem with homogeneous Dirichlet boundary condition. The Prager–Synge technique has wide applications in both a posteriori and a priori error estimation; see for example Liu-2022 for a recent successful application in validating the solution to the Navier–Stokes equation in 3D domains.

Let $V=H^{1}_{0}(\Omega)$, $W=L^{2}(\Omega)$, $\gamma$ be the identity, $a(u,v)=(\nabla u,\nabla v)$, and $b(\gamma u,\gamma v)=(u,v)$, where the parenthesis denote the $L^{2}(\Omega)$ scalar product. Further, let $\|\cdot\|_{0}$ stands for the $L^{2}(\Omega)$ norm.

Lemma 5

Let $u_{h}\in H^{1}_{0}(\Omega)$ and $p\in H(\operatorname{div},\Omega)$ be arbitrary. Then

[TABLE]

Proof

Using the particular forms of $a$ and $b$ for the Laplace eigenvalue problem, the negative norm in the definition of $\epsilon$ can be expressed with the aid of (41) as

[TABLE]

Denoting $\varphi=(I-P)v$, choosing arbitrary $p\in H(\operatorname{div},\Omega)$, applying the divergence theorem, Cauchy–Schwarz inequality, and the Friedrichs inequality $\|\varphi\|_{0}\leq\|\nabla\varphi\|_{0}/\sqrt{\lambda_{1}}$, we obtain

[TABLE]

Inserting this estimate to (42) together with inequality $\|\nabla\varphi\|_{0}\leq\|\nabla v\|_{0}$ finishes the proof.

Remark 13

The residual error $\|(I-P)(A-A_{h})u_{h}\|_{-1}$ is the same as the one in (CanDusMadStaVoh2017, , Definition 2.1), up to the unimportant factor $I-P$. However, the assumptions and error bounds for approximate eigenfunctions obtained in CanDusMadStaVoh2017 are different from our approach. We assume high-precision bounds for eigenvalues in the cluster of interest, which can be obtained by applying, for example, the Lehmann–Goerisch theorem with the FEM over refined meshes; while CanDusMadStaVoh2017 considers bounds for eigenvalues and eigenvectors simultaneously. Numerical results in §6 compare the accuracy of these two approaches and show that our approach is capable of high precision.

Lemma 6

Let $W_{h}$ be a finite dimensional subspace of $H(\operatorname{div},\Omega)$ and $\widehat{E}$ the space of approximate eigenvalues corresponding to the cluster of interest, then

[TABLE]

Further, if $A_{h}\widehat{E}\subset\operatorname{div}W_{h}$, we have

[TABLE]

Proof

Since $\widehat{P}V=\widehat{E}$ and using Lemma 5, we have

[TABLE]

In case that $A_{h}u_{h}\in\operatorname{div}W_{h}$, we can thus force $\operatorname{div}p_{h}+A_{h}u_{h}=0$ in seeking the minimizer $p_{h}$. Inserting this bound to the definition of $\epsilon$ finishes the proof.

Note that in the case of the finite element method the approximate eigenvectors $u_{h}$ forming $\widehat{E}$ are piecewise polynomial. A natural choice for $W_{h}$ in Lemma 6 is a suitable Raviart–Thomas space of higher order such that we can always find $p_{h}\in W_{h}$ satisfying $\operatorname{div}p_{h}+A_{h}u_{h}=0$. There are alternative methods to obtain this flux; see, e.g., CanDusMadStaVoh2017 ; Rep-2008-book . Assuming this choice, we consider linear operator $R_{h}:\widehat{E}\to W_{h}$ mapping $u_{h}$ to $p_{h}=R_{h}u_{h}$, which solves the following minimization problem:

[TABLE]

With the operator $R_{h}$, the evaluation of the residual error term $\epsilon$ reduces to solving a maximization problem

[TABLE]

The right-hand side of the above inequality corresponds to an eigenvalue problem for matrices with dimension $\operatorname{dim}(\widehat{E})$.

Remark 14

The bound in Theorem 5.3 utilizes relations (23) and (28) to obtain estimates for $\delta_{a}$ and $\tilde{\delta}$. Notice that a large cluster width $\lambda_{N}-\lambda_{n}$ will cause a drop of precision of this bound. Therefore, the following subsection derives bounds on $\delta_{a}$ and $\tilde{\delta}$ independent from the cluster width. We achieve it by introducing another application of Davis–Kahan’s theorem and estimating $\delta_{a}$ and $\tilde{\delta}$ directly.

5.3 Direct estimation of $\delta_{a}$: another application of Davis–Kahan’s theorem

Corresponding to the operator $T\gamma$ defined in (2), let us introduce a local discrete operator $\widehat{T}{\gamma}:\widehat{E}\to\widehat{E}$. For given $\hat{u}\in\widehat{E}$, the image $\widehat{T}{\gamma}\hat{u}$ is uniquely determined as the solution to problem

[TABLE]

Note that $\widehat{T}{\gamma}\hat{u}_{i}=\hat{\mu}_{i}\hat{u}_{i}$ for $i=n,\dots,N$, where $\hat{u}_{i}$ and $\hat{\mu}_{i}$ ($\hat{\mu}_{i}\geq\hat{\mu}_{i+1}$) are the discrete eigenpairs over $\widehat{E}$, i.e., $b(\gamma\hat{u}_{i},\gamma\hat{v})=\hat{\mu}_{i}a(\hat{u}_{i},\hat{v})$ for all $\hat{v}\in\widehat{E}$. Suppose the eigenvalues of $T$ corresponding to $E^{\perp}$ are excluded from $(\hat{\mu}_{N}-\sigma,\hat{\mu}_{n}+\sigma)$, then from Davis–Kahan’s theorem we have

[TABLE]

Note the following variational equation holds for any $\hat{u}\in\widehat{E}$ and $v\in V$,

[TABLE]

Then, we have

[TABLE]

In case of the Laplace eigenvalue problem, suppose $\widehat{E}$ is the approximation by using a conforming FEM space and $W_{h}$ the FEM approximation to $H(\mbox{div})$, for example, the Raviart–Thomas space. Then, with analogous argument as in Lemma 6, we have

[TABLE]

Compared with the approach utilizing (40) and (23) for estimating $\delta_{a}$, the estimate (45) is independent from the cluster width $|\lambda_{N}-\lambda_{n}|$. Such an estimate is expected to deal with clusters of relative large width. Note that, by further utilizing the relation in (27) and (34), the estimate (45) also applies to $\tilde{\delta}$ and thus for $\delta_{b}$. That is,

[TABLE]

This estimate of $\delta_{b}$ has the convergence rate as $O(h^{2})$ for liner conforming FEM approximation if solution $u$ has the $H^{2}$ regularity.

Remark 15

In practical computations with narrow cluster widths, the estimate (45) gives relatively inaccurate bounds on $\delta_{a}$ in comparison with Theorem 5.3. This is because the estimate (45) introduces to $\delta_{a}$ an additional overestimating factor of convergence rate $O(h)$. On the other hand, (43) introduces an overestimating factor of rate $O(h)$ for $\epsilon$ and the approach through Theorem 5.3, which depends on (40) and (23), finally leads to an overestimating factor for $\delta_{a}$ with the convergence rate $O(h^{2})$. Note that the estimator for $\delta_{a}$ itself is still of convergence rate $O(h)$.

To summarize, we extended the Davis–Kahan theorem for weakly formulated eigenvalue problems and derived bounds (40) for narrow clusters and (45) for wide clusters. These estimates are independent from the cluster index, its width and degree. The quality of these bounds is determined by the residual error and the gap between the concerned eigenvalue cluster and its neighbouring clusters. Thus, these estimates enjoy the same favourable properties as those in CanDusMadStaVoh2019 . As demonstrated by numerical examples in the following section, the estimates using the residual information provide sharper bound for $\delta_{a}$ and $\delta_{b}$ than the results in CanDusMadStaVoh2019 , especially for finite element approximations for eigenvalue problems of differential operators.

6 Numerical examples

This section numerically illustrates the accuracy of proposed bounds on the directed distances of spaces of exact and approximate eigenfunctions for matrix and Laplace eigenvalue problems. Refer to Remark 1 for how these problems fit to the considered abstract context.111See https://ganjin.online/xfliu/EigenVecEstimation for source codes and demonstrations of all presented examples.

6.1 The matrix eigenvalue problem

First, we apply estimates of the directed distance between spaces of eigenvectors to the generalized matrix eigenvalue problem $Ax=\lambda Bx$, where $A,B\in\mathbb{R}^{961\times 961}$ are the stiffness and mass matrices of the Laplacian in the unit square discretized by linear conforming finite element method on a uniform mesh with $h=\sqrt{2}/32$ (cf. Figure 2). In agreement with the abstract eigenvalue problem (1), the $a$- and $b$- norms of a vector $v$ are given as

[TABLE]

The approximate eigenvectors are computed by the MATLAB command eigs. The error of these approximations is estimated by Theorem 3.1. The needed two-sided bounds on eigenvalues are computed by using the interval arithmetic providing guaranteed enclosing intervals of exact eigenvalues.

The upper bounds on eigenvalues are easily obtained by the Rayleigh-Ritz method, while the lower bounds need more effort by applying Sylvester’s law of inertia; see the detailed implementation in Behnke1991 .

Table 1 lists the resulting two-sided bounds of eigenvalues and estimates of the directed distance of the corresponding subspaces of eigenvectors. Notice that since the approximate eigenvectors do no have better precision in the $b$-norm, the bound (23) has no advantage over (11).

6.2 Comparison with the $\sin\theta$ theorem of Davis and Kahan

Since the proposed bounds and Davis–Kahan’s estimates can be easily compared in the case of the standard matrix eigenvalue problem, we consider the problem $Ax=\lambda x$, where $A$ is the same matrix as in Subsection 6.1. Its eigenvalues and eigenvectors can be easily computed with high accuracy. The leading $13$ eigenvalues of $A$ are naturally clustered as in Table 2.

To generate approximate eigenvectors $\hat{u}_{i}$, we perturb the exact eigenvectors $u_{i}$ by uniform random errors of magnitude $10^{-4}$ and apply the Gram–Schmidt process to satisfy the orthonormality requirement of the Davis–Kahan’s method. Then, we compute approximate eigenvalues by the Rayleigh quotient as $\hat{\lambda}_{i}=\hat{u}_{i}^{T}A\hat{u}_{i}/\hat{u}_{i}^{T}\hat{u}_{i}$ and observe that differences $|\lambda_{i}-\hat{\lambda}_{i}|$ for the leading six eigenvalues are around $3\times 10^{-6}$.

First, we estimate the distance $\delta_{b}(E_{k},\widehat{E}_{k})$ between the exact space $E_{k}$ and its approximation $\widehat{E}_{k}$ by the bound (11). It is evaluated by using two-sided bounds of eigenvalues of $A$ that are obtained independently by using the same method as for the generalized eigenvalue problem in Subsection 6.1. The $\sin\theta$ theorem of Davis and Kahan davis1970rotation (i.e., Theorem 5.1) provides the estimate

[TABLE]

where $\hat{\mathbf{x}}=(\hat{u}_{n_{k}},\cdots,\hat{u}_{N_{k}})$ is the matrix with orthonormal column vectors that form $\widehat{E}_{k}$; $\hat{\Lambda}$ is the diagonal matrix with diagonal elements $\hat{\lambda}_{n_{k}},\cdots,\hat{\lambda}_{N_{k}}$; and $\delta_{\mathrm{spec}}$ is the spectral gap between the cluster of interest and the rest of the spectrum. For example, in case $k=2$, we have $\hat{\mathbf{x}}=(\hat{u}_{2},\hat{u}_{3})$, $\hat{\Lambda}=\operatorname{diag}(\hat{\lambda}_{2},\hat{\lambda}_{3})$, and $\delta_{\mathrm{spec}}=\min(\lambda_{4}-\lambda_{3},\lambda_{2}-\lambda_{1})$. Table 3 compares the computed bounds with the exact values $\delta_{b}(E_{k},\widehat{E}_{k})$ for $k=1,2,\dots,8$. The bound (11) gives better results than the Davis–Kahan’s estimate, except for the $7$th cluster $\{\lambda_{11}\}$. For the $8$th cluster, the Davis–Kahan’s estimate exceeds one; note that the value of the directed distance cannot be greater than one.

6.3 The unit square domain

Consider the Laplace eigenvalue problem with homogeneous Dirichlet boundary conditions in the unit square $\Omega=(0,1)^{2}$: find eigenvalues $\lambda_{i}\in\mathbb{R}$ and corresponding eigenfunctions $u_{i}\neq 0$ such that

[TABLE]

As we mention in Remark 1, the $a$- and $b$-norms correspond to the energy and $L^{2}(\Omega)$ norms, respectively.

The exact eigenpairs are known analytically to be

[TABLE]

These eigenvalues are either simple or multiple and we clustered them according to the multiplicity. The first four clusters are listed in Table 4. Since the exact eigenvalues are known, we use them to evaluate error bounds (10), (11), and (23). The quantity $\rho$ is chosen as $\rho=\lambda_{N+1}$. If the exact eigenvalues are not known, their two-sided bounds have to be employed.

This problem is discretized by the conforming finite element method of the first order. The finite element mesh $\mathcal{T}_{h}$ is chosen as the uniform triangulation consisting of isosceles right triangles; see the mesh with size $h=1/4$ in Figure 2.

For each cluster $K=1,2,3,4$, we compute bounds of on $\delta_{a}(E_{K},\widehat{E}_{K})$ and $\delta_{b}(E_{K},\widehat{E}_{K})$ by applying Algorithm I (Theorem 3.1 along with (23)) and Algorithm II (Theorem 5.3 along with (23)), respectively.

The convergence behavior of computed bounds from Algorithm I on a sequence of uniformly refined meshes, when applied to the exact directed distances of the four leading clusters, is shown in Figures 3 and 4. Particularly, the results confirm the expected optimal convergence rate $O(h)$ of the estimate (10), the convergence rate $O(h^{2})$ of $\delta_{b}(E_{K},\widehat{E}_{K})$, and the sub-optimal rate $O(h)$ of the bound (11).

The corresponding indices of effectivity ${\mathcal{I}(\delta)}$, for the the estimates computed by both algorithms are listed in Table 5 and Table 6. The estimate by Algorithm II provides impressively sharp bounds for the error of approximate eigenfunctions. Observed from the effectivity index, we can see that the estimate from Algorithm II on $\delta_{b}$ has an optimal convergence rate as $O(h^{2})$.

Note that the estimate computed by Algorithm I only utilizes the Rayleigh quotient of the approximate eigenfunction, while Algorithm II takes advantage of the residual error estimation along with the flux constructed by solving variational equations. The high precision of bounds on $\delta_{a}$ and $\tilde{\delta}$ based on Theorem 5.3 may be attributed to the relations (23) and (28), which come from the fundamental relation in (22). The a priori error estimates with the finite element projection can also help to obtain sharp estimates of $\delta_{a}$ and $\tilde{\delta}$ achieving effectivity indices almost one. This approach will be investigated in our subsequent paper. See also the discussion in preprint (liu-Vejchodsky-arxiv-v1, , §6).

Note that all bounds in this and in the following sections are computed in the floating-point arithmetic, and the influence of rounding errors is not taken into account. However, if needed, mathematically rigorous estimates could be obtained by employing the interval arithmetic moore2009introduction .

Remark 16

Error estimates of eigenfunctions for the Dirichlet Laplacian in the square domain are well discussed in the existing literature. For example, in CanDusMadStaVoh2017 , a computable estimate for the error in eigenfunctions is provided by solving an auxiliary problem and considering estimates of the residual, which are also used in Algorithm II. Note that paper CanDusMadStaVoh2017 measures the error of eigenfunctions with a different distance. Let us define the distance $\overline{\delta}(E,\widehat{E})$ by

[TABLE]

Note that $\overline{\delta}(E,\widehat{E})=\overline{\delta}(\widehat{E},E)$, only if $E$ and $\widehat{E}$ have the same dimension and the cluster width is zero. Paper CanDusMadStaVoh2017 utilizes this distance for the case of simple eigenvalues. It is easy to see that $\tilde{\delta}\leq\overline{\delta}$. By confirming the utilization of $\tilde{\delta}$ in Theorem 3 and Theorem 5.2, the arguments therein also hold by replacing $\tilde{\delta}$ with $\overline{\delta}$. Thus, we can apply the bound for $\tilde{\delta}$ to directly estimate $\overline{\delta}$. That is,

[TABLE]

In Table 7, the estimate of $\overline{\delta}$ through (49) and the estimate of $\delta_{b}$ from Theorem 3.1 and Theorem 5.3 is applied to the first cluster, and the effectivity indices of the obtained error estimates are compared with the ones of CanDusMadStaVoh2017 . Note that in this comparison, we follow the mesh type used in CanDusMadStaVoh2017 rather than the type in Figure 2 and the eigenvalue bounds for the clusters are different from the ones used in CanDusMadStaVoh2017 . Further note the convergence rate $O(h)$ of the overestimation factor of CanDusMadStaVoh2017 . With the analogous argument as in Remark 15, the locally reconstructed flux can also be used in (43) to obtain sharp estimates with the quadratic convergence rate of the overestimating factor.

6.4 The L-shaped domain

We consider Laplace eigenvalue problem (47) in the L-shaped domain $\Omega=(-1,1)^{2}\setminus(-1,0]^{2}$ to present the standard example with singularities of eigenfunctions and also to demonstrate the versatility of the proposed method. We solve this problem by using both the classical linear conforming finite element space and the extended space with two more basis functions

[TABLE]

Since the exact eigenvalues are not known, bounds (10), (11), and (23) are evaluated by using two-sided bounds on eigenvalues, which were computed in liu2014high and we list them in Table 8. The first four eigenvalues are simple and form trivial clusters.

The initial finite element mesh is displayed in Figure 5. First, we apply Algorithm I to the four leading eigenvalue clusters. Figure 6 shows the obtained bounds (10) and (23) on the energy distance $\delta_{a}$ and Figure 7 shows the bound (11) on the $L^{2}(\Omega)$ distance $\delta_{b}$. The results confirm that both the distance $\delta_{a}$ and its bounds converge with a slower rate for the standard finite element space and have the optimal speed of convergence $O(h)$ for the extended finite element space.

The indices of effectivity for the estimates computed by the two algorithms are listed in Table 9 and the comparison to the result of CanDusMadStaVoh2017 is shown in Table 10. Note the exact eigenfunctions are unknown and the approximate ones over quite refined mesh are used as the exact eigenfunctions in the effectivity index computation. The stable values of the effectivity indices indicate that the bound computed by Algorithm II has the optimal rate of convergence.

6.5 The dumbbell shaped domain

Finally, we consider Laplace eigenvalue problem (47) in a dumbbell shaped domain consisting of two unit squares connected by a bar of width $0.02$ and length $0.1$, see Figure 8, where also the initial mesh is depicted. This example is interesting due to singularities of eigenfunctions in re-entrant corners and especially due to tight clusters of eigenvalues.

The exact eigenpairs are not known, but the eigenvalues are expected to be close to eigenvalues for the union of two squares, i.e., two eigenvalues close to $2\pi^{2}\approx 19.739$, four eigenvalues close to $5\pi^{2}\approx 49.348$, etc. In order to compute high precision two-sided bounds for these eigenvalues, we combine the Crouzeix–Raviart nonconforming finite elements and the Lehmann–Goerisch method as proposed in Liu2015 . The resulting two-sided bounds obtained on a fine mesh and finite element spaces of the third order are presented in Table 11.

Table 11 also shows the chosen division of the first twelve eigenvalues into four clusters. Note that eigenvalues $\lambda_{3}$ and $\lambda_{4}$ are strictly separated from $\lambda_{5}$ and $\lambda_{6}$. Therefore, they could be considered as two separate clusters, but then the spectral gap between them would be small and the factor $\rho-\lambda_{n}$ in (10) and (11) would yield large overestimation. For this reason, all four eigenvalues $\lambda_{3},\dots,\lambda_{6}$ are considered in one cluster.

We apply compute the bounds on $\delta_{a}(E_{K},\widehat{E}_{K})$ and $\delta_{b}(E_{K},\widehat{E}_{K})$ for the four clusters $K=1,2,3,4$ as we did for the unit square domain. To investigate the rate of convergence of estimate computed by Algorithm I, the bounds (10), (11), and (23) are displayed in Figure 9. The first and the third cluster are very tight and we observe the first order convergence. However, the convergence curves for the second and the fourth cluster bend due to the non-negligible width of these clusters.

The comparison of effectivity indices of the estimates from Algorithm I and II are shown in Table 12. Here, the solution solved by the conforming Lagrange FEM with degree 4 is regarded as the exact solution in calculating the effectivity index. The mesh size $h$ here denotes the maximum edge length of the mesh. The presented results show that Algorithm II considerably improves the accuracy of the bounds, but has higher computational cost. The computational cost can be decreased by computing only approximate minimum in (44) using local flux reconstructions as in CanDusMadStaVoh2017 ; CanDusMadStaVoh2019 . See also study Vejchodsky2018b showing that the local flux reconstructions are considerably faster to compute than the global ones, while the two approaches provide almost the same accuracy.

7 Conclusions

For the abstractly formulated eigenvalue problem of compact operators over Hilbert spaces, we proposed two algorithms to derive guaranteed upper bounds on the directed distance between sub-spaces of exact and approximate eigenfunctions. The two presented algorithms have their own advantages for eigenvalue problems of different settings. Algorithm I only utilizes the Rayleigh quotients of the approximate eigenvectors, while Algorithm II provides sharper bounds by further considering the residual error estimation of discretized operators.

The bounds of both algorithms are independent of the discretization method and apply to arbitrary conforming approximations of eigenfunctions. Bounds on the total error of these approximations are easily computed by using solely the two-sided bounds on exact eigenvalues and the approximate eigenfunctions themselves. The derived bounds can be straightforwardly applied to, for example, the (generalized) matrix, Laplace, Steklov, and many other eigenvalue problems.

Numerical examples for the Laplace eigenvalue problem discretized by the finite element method show that the bound from Algorithm I in the $L^{2}$ norm converges with a sub-optimal rate, while the one from Algorithm II has the optimal rate. In a subsequent work we would like to employ the Aubin–Nitsche technique and the explicit a priori error estimate for the energy projection LiuOis2013 in order to derive bounds in the $L^{2}$ norm with the optimal rate of convergence, without solving the dual problem as required in Algorithm II.

Conflict of interest

The authors declare that they have no conflict of interest.

**Acknowledgement ** The authors greatly appreciate the valuable referees’ comments. Thanks to them the original manuscript considerably improved. The first author also show thanks to Dr. Yuji Nakatsukasa from Oxford University for his introduction of Davis–Kahan’s method.