On the convergence of iterated penalty methods for structure-preserving discretizations of saddle point problems

Problem Statement
This paper addresses the numerical solution of linear generalized saddle point systems of the form:
$$\begin{aligned} a(u,v) + (Dv, p)_Q &= F(v) \quad \forall v \in V, \\ (Du, q)_Q &= G(q) \quad \forall q \in Q, \end{aligned}$$
where $V$ and $Q$ are Hilbert spaces, $a(\cdot,\cdot)$ is a continuous bilinear form, and $D: V \to Q$ is a surjective bounded linear operator. The authors focus specifically on structure-preserving discretizations, where the discrete spaces $V_h \subset V$ and $Q_h \subset Q$ satisfy the condition $D V_h = Q_h$. While such schemes offer advantages like automatic satisfaction of the discrete inf-sup condition and improved error estimates, constructing an explicit basis for the image space $Q_h = D V_h$ is often computationally prohibitive or an open problem.

The paper investigates the iterated penalty method (also known as the Uzawa iteration on an augmented Lagrangian formulation) as a solution strategy that avoids the need for an explicit basis of $Q_h$. The method iteratively solves a sequence of penalized problems to converge to the solution of the original saddle point system.

Methodology
The authors analyze the method through three primary lenses:

1. Stability of the Penalized System: They first establish stability estimates for a single step of the penalty method (a perturbed saddle point problem) with penalty parameter $\varepsilon$. They derive explicit stability constants dependent on the continuity of $a(\cdot,\cdot)$, the inf-sup constant of $D$, and the coercivity properties of $a(\cdot,\cdot)$ on the kernel of $D$.
2. Convergence Analysis: The iterated penalty method is formulated as an operator iteration. The authors derive error propagation operators for both the primal variable ($u$) and the dual variable ($p$). They establish geometric convergence rates under various assumptions on the bilinear form $a(\cdot,\cdot)$:
   • General Case: $a(\cdot,\cdot)$ satisfies an inf-sup condition on the kernel of $D$.
   • Symmetric Non-negative (A1): $a(\cdot,\cdot)$ is symmetric and non-negative, coercive on the kernel.
   • Coercive (A2/A3): $a(\cdot,\cdot)$ is coercive on the whole space or symmetric and coercive.
3. Implementation Strategy: A critical contribution is the development of an implementation scheme (Algorithm 1) that eliminates the need for a basis of $Q_h$. The authors show that the Lagrange multiplier updates and stopping criteria can be computed using only the primal space $V_h$ and the operator $D$, provided $G=0$ or $G$ can be approximated via a secondary penalty iteration.

Key Contributions

• Novel Convergence Estimates: The paper provides new, explicit geometric convergence rates for the iterated penalty method. Unlike prior literature which often assumed $a(\cdot,\cdot)$ was symmetric positive definite (SPD) or coercive on the kernel, this work assumes only that $a(\cdot,\cdot)$ satisfies an inf-sup condition on the kernel of $D$.
• Sharp Stability Constants: The authors derive sharper stability estimates for the penalized/perturbed saddle point problems than those currently available in standard references (e.g., Brenner and Scott, Fortin and Glowinski). These estimates are explicit in the stability parameters ($\alpha_X, \beta_X, M_a$) and remain valid as the penalty parameter approaches zero or the inf-sup constant becomes small.
• General $G \neq 0$ Implementation: The paper presents a practical algorithm for the general case where the constraint source term $G(q)$ is non-zero. This includes a procedure to approximate the Riesz map $R_Q^{-1}G$ without a basis for $Q$, making the method applicable to a broader class of problems.
• Refined Analysis for Indefinite Cases: For the case where $a(\cdot,\cdot)$ is symmetric and satisfies the inf-sup condition (A3), the authors perform a refined analysis using generalized eigenfunctions. This analysis explains the behavior of the method when the inf-sup constant $\beta_X$ is arbitrarily small, revealing that while the constraint residual may converge, the error in the primal and dual variables may remain large if $\beta_X$ is unknown.

Numerical Results
The theory is validated through three model applications:

1. Hodge Decompositions: Solving saddle point systems arising from discrete de Rham complexes (Nédélec and Raviart–Thomas elements) to compute Hodge decompositions. The method successfully computed orthogonal components with iteration counts consistent with the theoretical geometric rate.
2. Fourth-Order Wave Equation: A $C^1$-spline discretization of a fourth-order problem. The study demonstrated that as the frequency $\omega$ approaches an eigenvalue (causing the inf-sup constant to degrade), the iteration count increases as predicted by the theory.
3. Incompressible Flow (Scott–Vogelius): The method was applied to the Stokes equations using Scott–Vogelius elements on meshes where the inf-sup constant $\beta_X$ could be made arbitrarily small by varying a geometric parameter $\delta$. The results confirmed that while the method converges geometrically, the stopping criterion based on the constraint residual can be satisfied prematurely when $\beta_X$ is very small, leading to large errors in the solution variables.

Significance and Claims
The authors claim that their work provides a more robust theoretical foundation for the iterated penalty method in structure-preserving discretizations. Specifically:

• The convergence results are explicit in the standard stability parameters, allowing practitioners to predict convergence behavior based on problem data.
• The assumption that $a(\cdot,\cdot)$ only needs to satisfy an inf-sup condition on the kernel (rather than being SPD or coercive) generalizes previous results found in works by Fortin & Glowinski, Brenner & Scott, and Awanou & Lai.
• The derived stability constants for the penalized system are sharper than existing literature, particularly in the limits of small penalty parameters or small inf-sup constants.
• The proposed implementation scheme removes the bottleneck of constructing bases for the image space $Q_h$, enabling the use of high-order structure-preserving methods (like $C^1$ splines or Scott–Vogelius elements) in practical software without explicit basis generation.

The paper concludes by noting that while the method is effective, caution is required when $\beta_X$ is small; the stopping criterion based on the constraint residual may be satisfied while the actual solution error remains significant.