Technical Summary: Optimal Complexity of Adaptive FEM for Symmetric Second-Order Linear Elliptic PDEs Driven by Non-Residual Estimators

1. Problem Statement

The paper addresses the adaptive finite element method (AFEM) for symmetric second-order linear elliptic partial differential equations (PDEs) of the form:
$$-\text{div}(A \nabla u^\star) + c u^\star = f - \text{div}(\mathbf{f}) \quad \text{in } \Omega, \quad u^\star = 0 \text{ on } \partial\Omega$$
where $\Omega \subset \mathbb{R}^d$ is a bounded polyhedral Lipschitz domain. The primary challenge addressed is the simultaneous control of two sources of error:

1. Discretization Error: Arising from the finite element approximation on a mesh.
2. Algebraic Error: Arising from the inexact solution of the resulting linear systems via iterative solvers.

Unlike previous works that often assume exact solutions to the discrete systems or rely exclusively on residual-based error estimators, this work focuses on non-residual error estimators (such as averaging-based ZZ-type estimators and equilibrated flux estimators) while integrating an inexact iterative algebraic solver into the adaptive loop. The goal is to prove that the adaptive algorithm achieves unconditional full R-linear convergence and optimal complexity with respect to the total computational cost.

2. Methodology and Framework

2.1. The Adaptive Algorithm

The authors propose an adaptive algorithm (Algorithm A) that loops through four modules: SOLVE, ESTIMATE, MARK, and REFINE.

• SOLVE & ESTIMATE: These modules are intertwined. An iterative algebraic solver (contractive map $\Psi_\ell$) is applied until a stopping criterion based on a computable algebraic error estimator $\zeta_\ell$ is met relative to the discretization error estimator $\mu_\ell$. Specifically, the solver stops when $\zeta_\ell(u^k_\ell) \leq \lambda \mu_\ell(u^k_\ell)$.
• MARK: A set of elements is marked using the Dörfler marking criterion based on the non-residual estimator $\mu_\ell$.
• REFINE: The marked elements are refined using Newest-Vertex Bisection (NVB).

2.2. Key Assumptions

The analysis relies on abstract properties of the error estimator $\mu_\ell$ and the algebraic solver:

• Contractive Solver: The iterative solver satisfies a contraction property $|||u^\star_\ell - \Psi_\ell(v_\ell)||| \leq q_{ctr} |||u^\star_\ell - v_\ell|||$ with $0 < q_{ctr} < 1$, independent of the mesh size.
• Non-Residual Estimator Properties: The estimator $\mu_\ell$ must satisfy:
  1. Local Equivalence: $\mu_\ell$ is locally equivalent to the standard residual-based estimator $\eta_\ell$ for the exact discrete solution $u^\star_\ell$. Specifically, $\eta_\ell(U_\ell; u^\star_\ell) \lesssim \mu_\ell(T^m_\ell[U_\ell]; u^\star_\ell)$ and vice versa, where $T^m_\ell$ denotes a patch of level $m$.
  2. Weak Stability: $\mu_\ell$ satisfies a stability condition relating the estimator at two different discrete functions, involving a patch level $r$.
• Mesh Refinement: The use of NVB ensures standard properties like shape regularity, overlay estimates, and mesh-closure estimates.

2.3. The Quasi-Error

The central quantity analyzed is the quasi-error $M^k_\ell$, defined as the sum of the algebraic error and the discretization error estimator:
$$M^k_\ell := |||u^\star_\ell - u^k_\ell||| + \mu_\ell(u^\star_\ell)$$
Note that $u^\star_\ell$ (the exact FE solution) is never computed; the term $\mu_\ell(u^\star_\ell)$ is a theoretical construct used for analysis, while the computable counterpart $\zeta_\ell(u^k_\ell) + \mu_\ell(u^k_\ell)$ is shown to be equivalent.

3. Key Contributions and Results

3.1. Unconditional Full R-Linear Convergence

The primary theoretical result (Theorem 8) establishes that the quasi-error $M^k_\ell$ converges unconditionally and R-linearly. That is, there exist constants $C_{lin} > 0$ and $0 < q_{lin} < 1$ such that for any two indices $(\ell, k)$ and $(\ell', k')$ in the adaptive history:
$$M^k_\ell \leq C_{lin} q_{lin}^{|(\ell, k)| - |(\ell', k')|} M^{k'}_{\ell'}$$
Crucially, this convergence holds for any choice of adaptivity parameters ($\theta, C_{mark}, \lambda$). This removes the need for "sufficiently small" parameters to guarantee convergence, a common restriction in prior literature involving inexact solvers.

3.2. Optimal Complexity

The paper proves that the decay rate of the quasi-error with respect to the total computational cost (measured by the cumulative number of degrees of freedom and solver steps) is optimal.

• Theorem 15: If the adaptivity parameters $\theta$ and $\lambda$ are chosen sufficiently small, the algorithm achieves the optimal convergence rate. Specifically, the decay rate of the quasi-error matches the best possible approximation rate in the nonlinear approximation class $\mathcal{A}^s$.
• The result implies that the algorithm does not waste computational resources on unnecessary solver iterations or mesh refinements, provided the parameters are tuned correctly.

3.3. Application to Specific Estimators

The abstract framework is applied to two specific classes of non-residual estimators, proving they satisfy the required local equivalence and stability assumptions:

1. ZZ-Type Averaging Estimators: Based on the seminal work of Zienkiewicz and Zhu. The paper proves local equivalence to the residual estimator for arbitrary polynomial degrees $p \geq 1$ (Theorem 16).
2. Equilibrated Flux Estimators: Based on local flux reconstruction (e.g., Raviart-Thomas elements). The paper establishes local equivalence and weak stability for these estimators (Theorem 22), noting their $p$-robustness.

3.4. Numerical Experiments

Section 6 presents 2D numerical experiments on an L-shaped domain (a problem with a singularity). The experiments compare:

• Standard residual-based estimators.
• ZZ-type estimators.
• Equilibrated flux estimators.

The results confirm:

• All three estimators generate comparable meshes with refinement concentrated at the singularity.
• Both non-residual estimators achieve optimal convergence rates with respect to the number of degrees of freedom and cumulative runtime.
• The equilibrated flux estimator demonstrates superior efficiency indices (close to 1) and $p$-robustness, though it may require more solver iterations due to tighter stopping criteria.

4. Significance and Relation to Literature

The authors position their work as a unification and extension of existing literature:

• vs. [KS11, CN12]: Unlike these works which focus on optimal rates with exact solvers, this paper includes inexact solvers and focuses on optimal complexity (cost vs. error). Furthermore, this work avoids the restrictive assumptions of [CN12] (sufficiently fine initial mesh, interior-node property) and [KS11] (lowest-order FEM, refining neighbors).
• vs. [CKNS08, BM09, CFPP14]: These works typically require exact FE solutions or rely on perturbation arguments that only guarantee convergence for small parameters. This paper provides unconditional convergence for any parameter choice.
• vs. [BFM+25]: While [BFM+25] extends optimal complexity to inexact solvers, it is restricted to residual-based estimators. This paper is the first to extend these results to non-residual estimators (ZZ and equilibrated flux), which are widely used in practice but analytically more challenging due to their lack of direct residual structure.

Core Innovation: The paper overcomes the difficulty that local equivalence between non-residual and residual estimators typically holds only for the exact discrete solution (which is never computed). By employing a subtle modification of the analysis in [BFM+25] and leveraging the weak stability of the non-residual estimator, the authors bridge the gap between the computed inexact solution and the theoretical exact discrete solution, thereby proving unconditional convergence and optimal complexity.

5. Conclusion

This work provides a rigorous mathematical foundation for using non-residual error estimators in adaptive finite element methods with inexact solvers. It demonstrates that, under general assumptions, these methods are not only unconditionally convergent but also achieve optimal computational complexity. This validates the practical use of popular estimators like ZZ-type and equilibrated flux estimators in adaptive algorithms where solving linear systems exactly is computationally prohibitive.