Rigorous certification of auxiliary-subspace a posteriori error estimates for p/hp finite element approximations of elliptic eigenvalue clusters

High-order (p and hp) finite element discretizations of self-adjoint elliptic eigenvalue problems admit computable a posteriori estimates of collective cluster errors — the Hausdorff distance between true and computed eigenvalue clusters and the gap between the corresponding invariant subspaces — via approximate error functions computed in an auxiliary, hierarchically enriched subspace (Giani–Grubišić–Hakula–Ovall, 2021). Two pillars of this estimator are unproven: the reliability constant in the two-sided source-error bound depends on the polynomial degree p in an unquantified way, and the computable Gram-matrix surrogate is used in place of the uncomputable Gram matrix on purely empirical evidence. The problem is to give the estimator a rigorous foundation — p-independent (or sharply p-explicit) reliability constants, computable bounds on the surrogate defect, oscillation-inclusive certification for general coefficients, and cluster-size-robust alternatives to pessimistic trace bounds — with constants whose dependencies on mesh, polynomial degree, cluster size, and coefficient regularity are explicit.

Background

Consider a self-adjoint, positive definite elliptic eigenvalue problem on a bounded domain $\Omega$, written in operator form $A\psi=\lambda\psi$ with $Av=-\nabla\cdot(\mathbf{A}\nabla v)+bv$ and boundary conditions — homogeneous Dirichlet, Neumann, or mixed Dirichlet/Neumann — incorporated in the trial space $\mathcal H\subset H^1(\Omega)$; $A$ is defined through the energy bilinear form $B(v,w)=(\mathbf{A}\nabla v,\nabla w)+(bv,w)$, and $|\!\|\!v\|\!\|=B(v,v)^{1/2}$ denotes the energy norm. In the $p$- and $hp$-versions of the finite element method ($p$/hp-FEM), accuracy is improved by raising the local polynomial degree $p$ of the element functions, on a fixed or graded mesh; the discretization theory of the source is developed for $\Omega\subset\mathbb R^2$ with Lipschitz boundary, partitioned into shape-regular (curvilinear) triangles and quadrilaterals, while the abstract estimator theory underlying it is not dimension-dependent. Hierarchical bases, in which the degree-$m$ basis is nested inside the degree-$(m+1)$ basis, are the standard tool, and high-order discretizations can deliver exponential convergence for eigenvalue computations on suitably graded meshes.

When several eigenvalues form a tight cluster $\Lambda=\{\lambda_1,\dots,\lambda_r\}$ (by symmetry, near-multiplicity, or isospectral domain constructions), errors of individual eigenpairs are ill-conditioned and the natural error measures are collective: the Hausdorff distance $\operatorname{dist}(\Lambda,\hat\Lambda)$ between the true and computed clusters $\hat\Lambda=\{\hat\mu_1,\dots,\hat\mu_r\}$, and the subspace gap between the invariant subspace $E=E(\Lambda)$ and its approximation $\hat E$, i.e. the operator norm of the difference of the orthogonal projectors onto the two subspaces. Grubišić and Ovall (2009) reduced such estimates to source-problem residuals: eigenpair errors are controlled by $\|\hat\varphi-A^{-1}(\hat\mu\hat\varphi)\|$ times a cluster separation factor $C(\hat\mu,\Lambda)$ determined by the distance of $\hat\mu$ from the rest of the spectrum.

Giani, Grubišić, Hakula and Ovall (2021), extending Hakula, Neilan and Ovall (2017), made these estimates computable for $p$/hp discretizations by the auxiliary subspace technique. Given data $f=\hat\mu_j\hat\varphi_j$ built from a computed eigenpair, let $u(f)=A^{-1}f$ solve the associated source problem and $\hat u(f)\in V$ its finite element approximation in the eigenvalue space $V$. One computes an approximate error function $\varepsilon(f)$ in an auxiliary space $W$ on the same mesh (a hierarchical complement): where $V$ uses element functions of degree $m$ on an element, $W$ uses element functions of degree $m+2$ on that element, and where $V$ uses edge functions of degree $m$ on an edge, $W$ uses edge functions of degree $m+1$ on that edge; $\varepsilon(f)$ is determined by the computable projection equation $B(\varepsilon(f),v)=(f,v)-B(\hat u(f),v)$ for all $v\in W$. The central two-sided bound (Theorem 3.1 of the source) is $\|\varepsilon(f)\|\ \le\ \|u(f)-\hat u(f)\|\ \le\ c\left(\|\varepsilon(f)\|+\operatorname{osc}(R,r,\mathcal T)\right),$ where all norms are energy norms, $c$ provably depends only on the shape-regularity of the mesh $\mathcal T$ (a uniform bound on the aspect ratios of its elements) and the polynomial degree $p$ — while the numerical evidence that $c$ does not grow with $p$ has never been turned into a theorem — and $\operatorname{osc}$ is the residual oscillation, measuring how far the element and edge residuals are from polynomials of degree $p-1$. From the functions $\varepsilon_j=\varepsilon(\hat\mu_j\hat\varphi_j)$ one builds the computable matrix $\tilde H$ with entries $\tilde H_{ij}=B(\varepsilon_j,\varepsilon_i)$, a surrogate for the Gram matrix $H_{ij}=B((I-S)\hat\varphi_j,(I-S)\hat\varphi_i)$ of the true eigenvector errors ($S$ the orthogonal projector onto $E$); together with $G_{ij}=B(\hat\varphi_j,\hat\varphi_i)$ it enters the exact identity $\sup_{v\in\hat E}\|(I-S)v\|^2/\|v\|^2=\lambda_{\max}(G^{-1}H)$ for the squared subspace gap. Trace-type estimates such as $\sum_j(\hat\mu_j-\mu_j)\le c^2C^2\sum_j\|\varepsilon_j\|^2$ are rigorous modulo the constant $c$ (for piecewise constant coefficients, where the oscillation vanishes); the sharper estimates $\operatorname{dist}(\Lambda,\hat\Lambda)\approx\lambda_{\max}(\tilde H)$ and $\operatorname{gap}(E,\hat E)\approx\lambda_{\max}(G^{-1}\tilde H)$ rest only on the heuristic that the eigenvalues of $\tilde H$ are close to those of the uncomputable $H$. The paper proves a Bauer–Fike bound (the classical matrix-perturbation theorem bounding how far eigenvalues of a matrix can move under an additive perturbation) converting $\|H-\tilde H\|$ into a spectral distance between the pairs $(H,G)$ and $(\tilde H,G)$, but gives no bound on $\|H-\tilde H\|$ itself. Extensive 2D experiments (slit disk, mixed-boundary and isospectral "bridge" domains — domain pairs constructed so that their spectra coincide exactly) show effectivities — ratios of estimate to true error — near 1.

What is missing is certification. Every one of the limitations itemized by the authors — the unproven $p$-independence of $c$, the unbounded surrogate defect $\|H-\tilde H\|$ (the sharper estimates rest only on empirical evidence), the lack of oscillation control for general coefficients (the certified estimates assume piecewise constant $\mathbf A$ and $b=0$), and the pessimism of trace bounds when the cluster size $r$ is large — is a gap in rigor of the same estimator rather than a separate application. The paper's experiments are all 2D and use fixed, non-adaptive discretizations, but adaptive refinement driven by these indicators is not proposed as an open problem in the source and is not part of the target stated here.

Problem Statement

Give the auxiliary-subspace error estimator of Giani–Grubišić–Hakula–Ovall (2021) for $p$/hp finite element discretizations of self-adjoint elliptic eigenvalue clusters a complete rigorous foundation. Concretely, in the estimator setting of the source — a bounded domain $\Omega\subset\mathbb R^2$ with Lipschitz boundary, shape-regular meshes $\mathcal T$ of (curvilinear) triangles and quadrilaterals, homogeneous Dirichlet, Neumann, or mixed Dirichlet/Neumann boundary conditions as incorporated in $\mathcal H\subset H^1(\Omega)$, and the hierarchical auxiliary space $W$ (element-function degrees exceeding those of the eigenvalue space $V$ by two, edge-function degrees by one; the abstract theory behind the estimator is dimension-independent):

1. ($p$-robustness) Prove that the constant $c$ in $\|\varepsilon(f)\|\le\|u(f)-\hat u(f)\|\le c(\|\varepsilon(f)\|+\operatorname{osc})$ can be chosen independent of the polynomial degree $p$, depending only on mesh shape-regularity — or determine the sharp rate at which it must grow with $p$.
2. (Surrogate defect) Establish a computable bound on $\|H-\tilde H\|$ that implies the guarantee that $\lambda_{\max}(\tilde H)$ and $\lambda_{\max}(G^{-1}\tilde H)$ approximate $\lambda_{\max}(H)$ and $\lambda_{\max}(G^{-1}H)$ with quantified error, in terms of computable quantities and the mesh/polynomial parameters.
3. (Oscillation) Extend the certified eigenvalue and subspace-gap estimates to general (non-piecewise-constant) coefficients by incorporating the residual oscillation term, with controlled effect on reliability and effectivity.
4. (Large clusters) Derive rigorously justified computable bounds for the cluster error measures whose quality does not degrade pessimistically as the cluster size $r$ grows.

Components 1–4 are the certification of the estimator's existing guarantees. The problem as posed is complete only when the estimator's cluster error bounds hold as theorems — with constants whose dependencies on mesh, polynomial degree, cluster size, and coefficient regularity are explicit.

References

1. Stefano Giani, Luka Grubišić, Harri Hakula, Jeffrey S. Ovall, A Posteriori Error Estimates for Elliptic Eigenvalue Problems Using Auxiliary Subspace Techniques, Journal of Scientific Computing 88(3), article 55, 2021. DOI: 10.1007/s10915-021-01572-2, https://doi.org/10.1007/s10915-021-01572-2 (preprint: arXiv:2009.06677, https://arxiv.org/abs/2009.06677)
2. Harri Hakula, Michael Neilan, Jeffrey S. Ovall, A Posteriori Estimates Using Auxiliary Subspace Techniques, Journal of Scientific Computing 72(1):97–127, 2017. DOI: 10.1007/s10915-016-0352-0, https://doi.org/10.1007/s10915-016-0352-0
3. Luka Grubišić, Jeffrey S. Ovall, On estimators for eigenvalue/eigenvector approximations, Mathematics of Computation 78(266):739–770, 2009. DOI: 10.1090/s0025-5718-08-02181-9, https://doi.org/10.1090/s0025-5718-08-02181-9
4. Eric Cancès, Geneviève Dusson, Yvon Maday, Benjamin Stamm, Martin Vohralík, Guaranteed a posteriori bounds for eigenvalues and eigenvectors: multiplicities and clusters, Mathematics of Computation 89(326):2563–2611, 2020. DOI: 10.1090/mcom/3549, https://doi.org/10.1090/mcom/3549
5. Claudio Canuto, Adaptive hp-FEM for eigenvalue computations, Calcolo 56(4), article 39, 2019. DOI: 10.1007/s10092-019-0335-2, https://doi.org/10.1007/s10092-019-0335-2
6. Stefano Giani, Jeffrey Ovall, Gabriel Pinochet-Soto, Adaptive Refinement for Eigenvalue Problems Based on an Associated Source Problem, Journal of Scientific Computing 105(2), article 32, 2025. DOI: 10.1007/s10915-025-03065-y, https://doi.org/10.1007/s10915-025-03065-y (preprint: arXiv:2501.05311, https://arxiv.org/abs/2501.05311)
7. Jay Gopalakrishnan, Gabriel Pinochet-Soto, Reliable Eigenspace Error Estimation Using Source Error Estimators, Computational Methods in Applied Mathematics 26(2):185–209, 2026. DOI: 10.1515/cmam-2025-0161, https://doi.org/10.1515/cmam-2025-0161 (preprint: arXiv:2601.08051, https://arxiv.org/abs/2601.08051)
8. Harri Hakula, Harmonic extension elements: Eigenproblems and error estimation, Computers & Mathematics with Applications 163:27–41, 2024. DOI: 10.1016/j.camwa.2024.03.012, https://doi.org/10.1016/j.camwa.2024.03.012
9. Harri Hakula, On high-order finite element solution of eigenvalue problems on isospectral surfaces, Computers & Mathematics with Applications 168:22–32, 2024. DOI: 10.1016/j.camwa.2024.05.022, https://doi.org/10.1016/j.camwa.2024.05.022
10. Jiachuan Zhang, Ran Zhang, Xiaoshen Wang, A Posteriori Estimates of Taylor-Hood Element for Stokes Problem Using Auxiliary Subspace Techniques, Journal of Scientific Computing 93(1), article 16, 2022. DOI: 10.1007/s10915-022-01924-6, https://doi.org/10.1007/s10915-022-01924-6
11. Théophile Chaumont-Frelet, Alexandre Ern, Martin Vohralík, On the derivation of guaranteed and p-robust a posteriori error estimates for the Helmholtz equation, Numerische Mathematik 148(3):525–573, 2021. DOI: 10.1007/s00211-021-01192-w, https://doi.org/10.1007/s00211-021-01192-w
12. Geneviève Dusson, Yvon Maday, An overview of a posteriori error estimation and post-processing methods for nonlinear eigenvalue problems, Journal of Computational Physics 491:112352, 2023. DOI: 10.1016/j.jcp.2023.112352, https://doi.org/10.1016/j.jcp.2023.112352 (preprint: arXiv:2303.01273, https://arxiv.org/abs/2303.01273)
13. Luka Grubišić, Harri Hakula, High order approximations of the operator Lyapunov equation have low rank, BIT Numerical Mathematics 62(4):1433–1459, 2022. DOI: 10.1007/s10543-022-00917-z, https://doi.org/10.1007/s10543-022-00917-z
14. Yonah Conjungo Taumhas, Geneviève Dusson, Virginie Ehrlacher, Tony Lelièvre, François Madiot, Reduced basis method for non-symmetric eigenvalue problems: application to the multigroup neutron diffusion equations, 2023. arXiv:2307.05978, https://arxiv.org/abs/2307.05978