🔢 mathematics

Hierarchical proximal Galerkin: a fast $hp$-FEM solver for variational problems with pointwise inequality constraints

This paper introduces a fast, mesh-independent $hp$-FEM solver for variational problems with pointwise inequality constraints by combining the proximal Galerkin algorithm with hierarchical $p$ -finite element bases and block preconditioning, achieving $hp$-robust convergence and superior wall-clock performance compared to low-order methods.

Original authors: Ioannis P. A. Papadopoulos

Published 2026-08-07

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Ioannis P. A. Papadopoulos

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are an engineer trying to design the perfect bridge, a new drug delivery system, or a more efficient airplane wing. To do this, you use a powerful mathematical tool called a "solver." Think of a solver as a super-smart calculator that tries to find the best possible shape or path for your design by minimizing energy, much like how a ball naturally rolls to the lowest point in a valley. However, real life is messy. You can't just let the ball roll anywhere; it might hit a rock (an obstacle) or be forced to stay within a specific speed limit (a gradient constraint). These rules turn a simple math problem into a complex puzzle where the solution has to obey strict, point-by-point rules.

For decades, scientists have solved these puzzles using a method called "Finite Element Method" (FEM). The most common version of this method is like using a low-resolution digital photo: it breaks the problem into tiny, simple squares (like pixels) and uses simple, straight-line math to guess the answer. It's fast and reliable, but it's not very detailed. To get a sharper picture, you usually have to add millions of these tiny squares, which makes the computer work incredibly hard and take a long time. The big question in this field has been: Can we use "high-resolution" math—using fancy, curvy, high-order polynomials instead of simple straight lines—to get a much better answer without the computer crashing? The answer has traditionally been "no," because high-resolution math usually creates a massive, messy tangle of numbers that is too hard to untangle quickly.

This paper introduces a new, clever way to solve these constrained puzzles called the "Hierarchical Proximal Galerkin" (hpG) algorithm. The author, who led the work, has built a solver that acts like a master chef who can cook a gourmet, high-resolution meal (using very high-order math) just as fast as someone making a simple sandwich (using low-order math). They achieved this by combining a recent mathematical trick called the "proximal point method" with a special way of organizing the math that keeps the numbers from getting messy.

Here is the magic: The author discovered that by using a specific type of high-order math basis (think of it as a set of building blocks that are very good at hiding complexity), they could keep the computer's workload surprisingly light. They found that their new solver doesn't just work; it is often faster than the old, low-resolution methods. In their tests, they solved problems with polynomial degrees as high as 82 (which is like having 82 layers of detail in a single block) and found that they could reach the same level of accuracy as low-resolution methods in a fraction of the time. In some cases, their high-resolution solver was up to 100 times faster per step than the traditional low-resolution solver.

The paper explicitly argues against the long-held belief that high-order methods are too slow or too complex for problems with strict "obstacle" or "gradient" rules. They show that while low-resolution methods are simple, they are often inefficient because they need millions of tiny blocks to get a good answer. Their new method proves that you can get a much more accurate answer with far fewer blocks, and the computer actually finishes the job quicker. They tested this on various scenarios, including a 3D obstacle problem and a "thermoforming" problem (which models how a heated metal sheet is pressed into a mold), and the results held up. The author observed that the number of steps their solver needed to find the answer didn't get worse as they made the math more detailed or the grid finer, which is a rare and valuable property.

In short, this paper presents a fast, high-order solver that breaks the usual trade-off between speed and accuracy. It suggests that for complex problems with strict rules, using a "high-definition" mathematical approach isn't just possible—it's the smartest, fastest way to go. The author provides concrete evidence, including timing data, showing that their method outperforms traditional low-order solvers, offering a new, powerful tool for engineers and scientists who need precise solutions to difficult physical problems.

Technical Summary: Hierarchical Proximal Galerkin (hpG)

Problem Statement
The paper addresses the numerical solution of variational problems involving pointwise inequality constraints, specifically the obstacle problem ( $u \leq \phi$ ) and the generalized elastic-plastic torsion problem (gradient-type constraints, $|\nabla u| \leq \phi$ ). These problems arise in diverse fields such as optimal control, topology optimization, and elasticity theory. A significant challenge in solving these problems is the limited regularity of the solution ( $u \in H^s$ for $s < 5/2$ ), even when the domain and data are smooth. This regularity cap has traditionally favored low-order finite element methods (FEM), as high-order methods were perceived to suffer from increased computational costs, loss of sparsity, and mesh-dependent convergence rates. Furthermore, existing high-order approaches for these constraints often rely on direct solvers or exhibit mesh-dependent iteration counts, limiting their scalability.

Methodology
The author proposes the Hierarchical Proximal Galerkin (hpG) algorithm, a high-order solver that leverages the Latent Variable Proximal Point (LVPP) method introduced by Keith and Surowiec. The core methodology involves:

Discretization Strategy: The LVPP saddle-point systems are discretized using a hierarchical $p$ -FEM basis.
- The primal variable $u$ is discretized in a continuous $H^1$ -conforming space using hierarchical shape functions (weighted Jacobi polynomials and hat functions).
- The latent variable $\psi$ is discretized in a discontinuous $L^2$ -conforming space using scaled and translated Legendre polynomials.
- This specific pairing (continuous $p$ for $u$ , discontinuous $p-2$ or $p-1$ for $\psi$ ) is chosen to satisfy inf-sup stability while preserving sparsity.
Newton Linearization: The resulting nonlinear systems are linearized via Newton's method, yielding symmetric saddle-point linear systems. Crucially, the choice of discontinuous basis for $\psi$ ensures that the linearization of the nonlinear term results in a block-diagonal matrix structure, preventing the stiffness and mass matrices from becoming dense as the polynomial degree $p$ increases.
Preconditioning and Solvers:
- The linear systems are solved using a block-preconditioned Krylov subspace method (FGMRES for the outer loop, GMRES for inner Schur complement solves).
- A block-diagonal preconditioner is constructed for the Schur complement. For obstacle problems, this involves approximating the Schur complement using a spectral Galerkin basis to maintain sparsity. For gradient-type constraints, a simpler approximation is used.
- The solver exploits fast transforms (via Discrete Cosine Transform) for the assembly and application of matrices, particularly for the discontinuous basis, achieving quasi-optimal complexity.
- The approach supports $hp$-adaptivity, utilizing a posteriori error estimators to guide mesh refinement and polynomial degree increases.

Key Contributions
The paper makes the following specific contributions:

First Preconditioners for High-Order Constraints: The author develops the first preconditioners for Newton systems arising in proximal Galerkin methods for obstacle and gradient-type constraints. These preconditioners are shown to be exceptionally effective.
Sparsity Preservation: By utilizing a specific hierarchical basis pairing, the method retains sparsity in the discretized Newton systems even as $p \to \infty$ , overcoming a major barrier to high-order methods for these problems.
$hp$-Robustness: The algorithm exhibits mesh ( $h$ ) and degree ( $p$ ) independence in the number of nonlinear (Newton) iterations required for convergence. The inner Krylov iterations show only mild polylogarithmic growth with respect to $h \to 0$ and $p \to \infty$ .
Competitive Performance: The paper provides the first reported competitive wall-clock timing comparisons showing that high-order discretizations can be significantly faster than low-order counterparts for these specific problems, particularly when the problem exhibits a tensor-product structure.
High-Degree Capability: The solver successfully handles discretizations with very high polynomial degrees (up to $p=82$ on each element, resulting in total degrees up to 164 in 2D), which were previously unattainable for such constrained problems using iterative solvers.

Results
The author validates the hpG algorithm through five numerical examples:

1D Oscillatory Obstacle Problem: Demonstrates that $hp$-adaptive strategies achieve convergence rates of $O((h/p)^{10})$ , significantly outperforming uniform refinement. The solver requires a bounded number of Newton iterations regardless of refinement.
2D Obstacle Problem: Shows that high-order discretizations achieve smaller errors per degree of freedom (dof) and per unit of wall-clock time compared to low-order ( $p=1$ ) methods. The hpG solver was observed to be up to 24 times faster per linear solve than a low-order primal-dual active set (PDAS) strategy to reach the same error tolerance.
2D Gradient-Type Constraint: Confirms $hp$-robustness (74–78 Newton iterations) and superconvergence ( $O(h^p)$ ) when the mesh aligns with the discontinuity of the constraint. High-order methods ( $p > 1$ ) were found to be roughly 100 times faster than low-order methods for reaching a specific error.
Thermoforming Quasi-Variational Inequality (QVI): Applies the solver to a complex obstacle-type QVI with a fixed-point iteration. The method successfully handles the highest-order discretization of an elliptic obstacle-type QVI reported in the literature to date ( $p=82$ ).
3D Obstacle Problem: Extends the method to three dimensions, demonstrating that while Cholesky factorization becomes infeasible for large systems, the proposed iterative solver with AMG preconditioning remains effective.

Significance and Claims
The paper claims that the hpG algorithm fundamentally challenges the prevailing notion that low-order discretizations are the only viable option for pointwise-constrained problems due to solution regularity limits. The author asserts that:

High-order discretizations can induce errors orders of magnitude smaller per degree of freedom, particularly in the $H^1$ -norm.
Contrary to expectations, high-order methods can be computationally cheaper (in terms of wall-clock time) than low-order methods for these problems when combined with the specific sparsity-preserving discretization and preconditioning strategies proposed.
The method provides a robust, mesh-independent framework for solving variational inequalities with pointwise constraints, enabling the use of high-order accuracy in applications like topology optimization and elasticity where such constraints are prevalent.

The author remains modest regarding generalizations, noting that the current implementation relies on tensor-product (Cartesian) meshes and that extensions to non-Cartesian cells or $p$ -multigrid strategies are open areas for future work.

More like this