PDE-Constrained Variational Models Overview

• PDE-constrained variational models are formulations that integrate a PDE directly into the variational principle to enforce state, control, or parameter constraints.
• They employ methodologies such as reduced-space optimization, total variation and L1-regularization, and physics-informed penalties to manage inverse and control problems.
• These models have broad applications in optimal control, image registration, and quantum simulation, often utilizing adjoint techniques and accelerated numerical solvers.

PDE-constrained variational models are formulations in which a state, control, parameter, latent field, or constitutive law is determined by minimizing a functional while respecting a partial differential equation either exactly as a constraint, implicitly through a reduced solution operator, or softly through residual or energy penalties. Across the literature represented here, the term covers classical optimal control and inverse problems, state-regularized nonsmooth PDEs, weak-form finite-element formulations, total-variation and $L^1$-regularized models, accelerated optimization flows derived from variational principles, and recent physics-informed machine-learning and variational-quantum constructions (Funke et al., 2013, Cristinelli et al., 2023, Caflisch et al., 2013, Poudel et al., 1 Feb 2026, Surana et al., 2024). The unifying feature is that the governing PDE is not merely background structure: it is built directly into the variational principle and thereby determines admissibility, regularity, optimization geometry, and computational complexity.

1. Classical formulation and reduced-space structure

A standard PDE-constrained optimization problem is written as

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$

where $u$ is the state variable, $m$ is the control or parameter variable, $J(u,m)\in \mathbb R$ is the scalar objective functional, and $F(u,m)=0$ is the PDE constraint (Funke et al., 2013). When the PDE has a unique solution operator $u(m)$, the formulation reduces to the parameter-only problem

$\hat J(m)\equiv J(u(m),m),$

with optimization carried out over $m$ alone (Funke et al., 2013). This reduced viewpoint is also central in total-variation-regularized control,

$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$

where the PDE enters through the control-to-state map $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$0 and the smooth part $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$1 is composed with a nonsmooth regularizer (Cristinelli et al., 2023).

A canonical example is distributed elliptic control: $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$2 This displays the standard ingredients of the field: misfit, regularization, a PDE state law, and admissibility constraints on the control (Funke et al., 2013). A different but related reduced structure appears in transport-based registration and mass-preserving transport, where the state satisfies an advection or continuity equation and the optimization variable is a smooth stationary velocity field $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$3 regularized in a Sobolev norm (He et al., 9 Oct 2025).

In state-dependent inverse problems, the unknown is not a field over physical space but a constitutive relation indexed by the state itself. The abstract formulation is

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$4

with the specific specialization

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$5

so that the optimization variable lives on a state interval rather than on $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$6 or $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$7 (Bukshtynov, 16 Jan 2026). This reorganizes identifiability, gradient computation, and regularization around the attained range of the state.

2. Variational origins of the governing PDE

Many PDE-constrained variational models arise because the PDE is itself the Euler–Lagrange equation or gradient flow of an energy. In “PDEs with Compressed Solutions,” the starting point is the convex but nonsmooth functional

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$8

with $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$9 symmetric positive definite and $u$0 (Caflisch et al., 2013). Defining

$u$1

minimization yields the elliptic variational inclusion

$u$2

while the $u$3-gradient flow gives the parabolic differential inclusion

$u$4

In this template, the state itself is regularized by an $u$5 term, and the PDE is the first-order condition of a nonsmooth variational principle rather than an externally imposed state law (Caflisch et al., 2013).

A second variational mechanism is acceleration in function space. For energies of the form

$u$6

standard $u$7-gradient descent produces

$u$8

whereas PDE acceleration derives from the action

$u$9

and yields the second-order flow

$m$0

For constant coefficients this reduces to

$m$1

turning diffusion-type optimization PDEs into damped wave equations (Benyamin et al., 2018).

Weak-form variational formulations are equally central. For Poisson control, the strong PDE

$m$2

is expressed weakly as

$m$3

and the resulting UFL representation becomes both the modeling language and the differentiation substrate in automated adjoint frameworks (Funke et al., 2013). A more probabilistic variant appears in “Variational Autoencoding of PDE Inverse Problems,” where the weak residual is embedded into a variational inference objective through the Gaussian density

$m$4

so that the PDE acts as a weak-form decoder regularizer rather than a hard solve (Tait et al., 2020).

These constructions show that “PDE-constrained” does not refer to a single formalism. It includes exact constraints, reduced maps, weak residuals, and gradient-flow dynamics, provided the variational model is organized around PDE structure.

3. Regularization, sparsity, total variation, and state-space constraints

Regularization is not ancillary in this literature; it often defines the qualitative behavior of admissible solutions. In $m$5-regularized state models, the subgradient relation

$m$6

and the pointwise form

$m$7

create exact zero regions, compact support, and finite speed of propagation (Caflisch et al., 2013). In one dimension with $m$8, the elliptic and parabolic inclusions become

$m$9

and support estimates such as

$J(u,m)\in \mathbb R$0

quantify the compression effect (Caflisch et al., 2013). The same paper rewrites the divisible sandpile problem by replacing the positivity constraint with an $J(u,m)\in \mathbb R$1 term,

$J(u,m)\in \mathbb R$2

whose Euler–Lagrange equation is

$J(u,m)\in \mathbb R$3

This replaces a hard state constraint by a nonsmooth variational regularizer (Caflisch et al., 2013).

Total variation regularization occupies a parallel but distinct role. For piecewise constant controls on a triangulation $J(u,m)\in \mathbb R$4,

$J(u,m)\in \mathbb R$5

and the extremal points of the corresponding zero-mean TV unit ball are normalized characteristic functions of simple sets (Cristinelli et al., 2023). This convex-geometric characterization explains why TV-regularized controls tend to be piecewise constant with simple jump sets, and it enables atom-based algorithms driven by extreme points (Cristinelli et al., 2023). The same work proves that $J(u,m)\in \mathbb R$6 discretizations do not in general approximate isotropic TV; instead they converge to an anisotropic functional

$J(u,m)\in \mathbb R$7

with anisotropy induced by the fine-scale mesh geometry (Cristinelli et al., 2023).

State-dependent inverse problems introduce another regularization layer. The canonical reduced problem is

$J(u,m)\in \mathbb R$8

where the constitutive law $J(u,m)\in \mathbb R$9 is identifiable only on the attained state interval

$F(u,m)=0$0

Outside $F(u,m)=0$1, recovery is regularization- or model-class-driven (Bukshtynov, 16 Jan 2026). The survey emphasizes both $F(u,m)=0$2 and $F(u,m)=0$3 regularization,

$F(u,m)=0$4

as well as Sobolev gradients obtained from

$F(u,m)=0$5

on the state interval $F(u,m)=0$6 (Bukshtynov, 16 Jan 2026). A plausible implication is that regularization in PDE-constrained variational models is often inseparable from identifiability itself: it does not merely stabilize optimization, but selects among structurally indistinguishable solutions.

4. Optimality systems, adjoints, and computational frameworks

The standard first-order optimality conditions are encoded by the Lagrangian

$F(u,m)=0$7

with stationarity conditions

$F(u,m)=0$8

(Funke et al., 2013). In automated finite-element frameworks, the forward model is recorded as a tape of equation solves, and the adjoint is derived from that discrete representation rather than from a manually derived continuous adjoint. The resulting workflow repeatedly solves the forward PDE, evaluates the reduced functional, computes $F(u,m)=0$9 via the adjoint, and updates $u(m)$0 using a numerical optimizer (Funke et al., 2013).

Transport-dominated models make this reduced-space structure especially explicit. For image registration by stationary advection, the reduced gradient is

$u(m)$1

where $u(m)$2 solves the state equation

$u(m)$3

and $u(m)$4 solves the backward adjoint

$u(m)$5

(He et al., 9 Oct 2025). The solver contribution in that work is GA-NGMRES, which accelerates first-order fixed-point iterations by alternating nonlinear GMRES steps with plain fixed-point updates and yields runtimes up to $u(m)$6 faster than Newton–Krylov baselines in the reported experiments (He et al., 9 Oct 2025).

For diffeomorphic registration in LDDMM, the control may instead be the initial velocity $u(m)$7, with the whole trajectory generated by the EPDiff equation

$u(m)$8

The objective

$u(m)$9

is then constrained simultaneously by EPDiff and image or deformation transport equations (Hernandez, 2018). Gauss–Newton–Krylov optimization acts in the reduced variable $\hat J(m)\equiv J(u(m),m),$0, while band-limited vector fields reduce memory and runtime enough to make the geodesic PDE-constrained formulation practical (Hernandez, 2018).

The following table summarizes representative computational styles documented in the literature.

Model class	Constraint handling	Characteristic computation
Classical optimal control (Funke et al., 2013)	Hard PDE constraint, reduced functional	Automated discrete adjoints in UFL/FEniCS
TV-regularized control (Cristinelli et al., 2023)	Hard PDE constraint through $\hat J(m)\equiv J(u(m),m),$1, nonsmooth TV objective	FC-GCG with graph-cut atom oracle
$\hat J(m)\equiv J(u(m),m),$2-regularized state models (Caflisch et al., 2013)	PDE as Euler–Lagrange inclusion or gradient flow	Proximal shrinkage, Douglas–Rachford, FFT-based solves
Transport-constrained registration (He et al., 9 Oct 2025)	Hard transport PDE, reduced-space optimization	Preconditioned gradient descent accelerated by GA-NGMRES
Momentum-constrained LDDMM (Hernandez, 2018)	Hard EPDiff plus transport constraints	Inexact Gauss–Newton–Krylov in band-limited space
Weak-form latent inference (Tait et al., 2020)	Soft weak residual in probabilistic decoder	FEM assembly without repeated global solves

A recurring theme is that the adjoint remains the principal tool whenever exact PDE constraints are retained, whereas proximal splitting, graph cuts, or residual backpropagation dominate when the PDE is absorbed into a nonsmooth penalty or soft residual term.

5. Learning-based and quantum extensions

Recent work extends PDE-constrained variational modeling into scientific machine learning by treating network outputs as continuous fields and imposing physics through differentiable residuals or energies. In microscopy segmentation, the field

$\hat J(m)\equiv J(u(m),m),$3

is optimized with the composite objective

$\hat J(m)\equiv J(u(m),m),$4

where the reaction–diffusion residual is

$\hat J(m)\equiv J(u(m),m),$5

and the phase-field energy is

$\hat J(m)\equiv J(u(m),m),$6

(Poudel et al., 1 Feb 2026). The PDE is not solved exactly inside the training loop; it is weakly imposed as a soft residual and energy penalty on the network output. The reported gains are especially strong in the low-sample regime, including OOD Dice improvement from $\hat J(m)\equiv J(u(m),m),$7 to $\hat J(m)\equiv J(u(m),m),$8 at 10% training data (Poudel et al., 1 Feb 2026).

A related quantum line of work treats measured observables of a variational quantum circuit as a discretized field and penalizes PDE residuals at collocation points. The general physics loss is

$\hat J(m)\equiv J(u(m),m),$9

with measured outputs

$m$0

For finite-difference discretizations of a $m$1-th order PDE, each residual depends on at most $m$2 neighboring outputs, so the total loss is a sum of local cost terms (Hewage et al., 10 Apr 2026). That locality is the basis for the paper’s claim that PDE-constrained losses mitigate barren plateaus by inheriting the polynomial trainability of local cost functions and by inducing “constraint-induced landscape narrowing” (Hewage et al., 10 Apr 2026).

A second quantum direction addresses hardware noise in variational PDE solvers. There the soft constrained objective

$m$3

is studied under depolarizing, amplitude damping, and bit-flip noise. The reported result is that zero-noise extrapolation reduces absolute error by $m$4–$m$5 at $m$6, while constrained circuits retain $m$7–$m$8 higher fidelity than unconstrained ones at $m$9 (Hewage et al., 11 Apr 2026). This suggests that physics-based constraints may function as robustness mechanisms, not only as inductive biases.

Hybrid quantum-classical PINN-style hydrological models extend the same soft-constrained pattern to flood prediction. There the loss

$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$0

combines focal classification loss with a Saint-Venant continuity residual

$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$1

and Manning consistency penalty

$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$2

The reported model converges in about $\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$3 fewer training epochs and uses about $\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$4 fewer trainable parameters than the classical PINN baseline (Hewage et al., 10 Apr 2026).

6. Structural limits, computability, and open directions

A major theme across the literature is that PDE-constrained variational structure helps explain not only algorithms but also the limits of what can be identified or computed. In state-dependent inverse problems, the gradient with respect to a constitutive law is concentrated on state level sets, for example

$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$5

so only the interval of state values actually attained by the dynamics is identifiable (Bukshtynov, 16 Jan 2026). If a perturbation $\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$6 is supported outside the attained range, then the first variation vanishes (Bukshtynov, 16 Jan 2026). This is a structural limitation, not a numerical artifact.

Supremal formulations expose a different structural regime. For Navier–Stokes-constrained data assimilation, the finite-$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$7 objective

$\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$8

converges, as $\min_{u\in V} J(u):=F(Ku)+TV(u,\Omega),$9, to the worst-case functional

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$00

The resulting $$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$01 optimality system involves measure coefficients rather than a classical Euler–Lagrange PDE and concentrates on the sets where the maximal mismatch is attained (Clark et al., 2021). This shows that PDE-constrained variational models can move from integral control to uniform control, but at the price of markedly more singular optimality systems.

The most abstract structural analysis in the present corpus concerns algorithmic complexity. There the PDE is rewritten as minimization of a least-squares residual

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$02

and one studies computability of the minimizer through finite-dimensional approximation and discrete gradient flow (Cardona et al., 24 Oct 2025). Under convexity, coercivity, QGC, and Lipschitz gradient assumptions, the error decomposes as

$$\begin{aligned} \min_{u,m}\;& J(u,m) \ \text{subject to }& F(u,m)=0, \ & h(m)=0, \ & g(m)\le 0, \end{aligned}$$03

with the optimization error decaying exponentially in the iteration index under the stated assumptions (Cardona et al., 24 Oct 2025). When the PDE preserves analyticity, the solution is reported to be polynomial-time computable in the Sobolev norm, whereas loss of analyticity produces super-polynomial complexity blowup even when the input data are polynomial-time computable (Cardona et al., 24 Oct 2025). A plausible implication is that variational formulation alone does not guarantee tractability; tractability also depends on the regularity transmitted by the PDE solution operator.

Several open directions recur across the papers. The survey on state-dependent inverse problems identifies nonconvexity, state-space coverage, model-class dependence, scalable adjoint computation, sparse and indirect data, and uncertainty quantification as persistent challenges (Bukshtynov, 16 Jan 2026). The neural segmentation paper notes that its PDE priors are morphology priors rather than imaging-physics laws and that exact constraint satisfaction is not guaranteed under soft residual training (Poudel et al., 1 Feb 2026). The quantum papers likewise emphasize small-scale simulated validation, soft rather than hard enforcement, and unresolved scaling questions (Hewage et al., 10 Apr 2026, Hewage et al., 11 Apr 2026, Hewage et al., 10 Apr 2026). Taken together, these works support a broad but technically precise conclusion: PDE-constrained variational models constitute a family of formulations rather than a single method, and their behavior is determined jointly by constraint enforcement mode, regularization geometry, adjoint or proximal structure, discretization, and the analytical properties of the underlying PDE.