---
title: PDE-Constrained Variational Models Overview
url: https://www.emergentmind.com/topics/pde-constrained-variational-models
type: topic
---

# PDE-Constrained Variational Models Overview

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 [1302.3894] [2310.19777] [1311.5850] [2602.01069] [2405.16651]. 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 [1302.3894]. 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 [1302.3894]. 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 \(K\) and the smooth part \(F(Ku)\) is composed with a nonsmooth regularizer [2310.19777].

A canonical example is distributed elliptic control:
\[
\begin{aligned}
\min_{u,m}\;& \frac{1}{2}\|u-u_d\|_{L^2(\Omega)}^2 + \frac{\alpha}{2}\|m\|_{L^2(\Omega)}^2 \\
\text{subject to }& -\nabla^2 u = m \quad \text{in }\Omega, \\
&u=0 \quad \text{on }\partial\Omega, \\
&a\le m\le b \quad \text{in }\Omega.
\end{aligned}
\]
This displays the standard ingredients of the field: misfit, regularization, a PDE state law, and admissibility constraints on the control [1302.3894]. 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 \(v\) regularized in a Sobolev norm [2510.08782].

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
\[
\min_{p(\cdot)\in P} J(u,p) \quad\text{subject to}\quad \mathcal F(u,p)=0,
\]
with the specific specialization
\[
p=p(\zeta(u)),
\]
so that the optimization variable lives on a state interval rather than on \(\Omega\) or \(\Omega\times(0,T)\) [2601.10920]. 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
\[
E(u)=\int \frac12 (\nabla u)\cdot M(\nabla u)-u f+\gamma |u|\,dx,
\]
with \(M(x)\) symmetric positive definite and \(\gamma\ge 0\) [1311.5850]. Defining
\[
Au=-\nabla\cdot(M\nabla u),
\]
minimization yields the elliptic variational inclusion
\[
Au=f-\gamma p(u), \qquad p(u)\in \partial \|u\|_{L^1},
\]
while the \(L^2\)-gradient flow gives the parabolic differential inclusion
\[
u_t+Au=f-\gamma p(u),\qquad u(x,0)=g(x).
\]
In this template, the state itself is regularized by an \(L^1\) term, and the PDE is the first-order condition of a nonsmooth variational principle rather than an externally imposed state law [1311.5850].

A second variational mechanism is acceleration in function space. For energies of the form
\[
E[u] := \int_{\Omega} \Phi(x,\nabla u) + \Psi(x,u)\,dx,
\]
standard \(L^2\)-gradient descent produces
\[
u_t + L_z(x,u,\nabla u) - \operatorname{div}\big(\nabla_p L(x,u,\nabla u)\big)=0,
\]
whereas PDE acceleration derives from the action
\[
J[u] = \int_{t_0}^{t_1} k(t)\left( \frac12 \int_\Omega \rho(x)\,u_t^2\,dx - b(t)E[u] \right)\,dt
\]
and yields the second-order flow
\[
u_{tt}+a(t)u_t=-b(t)\rho(x)^{-1}\nabla E[u].
\]
For constant coefficients this reduces to
\[
u_{tt}+au_t=-\nabla E[u],
\]
turning diffusion-type optimization PDEs into damped wave equations [1810.00410].

Weak-form variational formulations are equally central. For Poisson control, the strong PDE
\[
-\nabla^2 u = m
\]
is expressed weakly as
\[
\langle \nabla u, \nabla v\rangle_\Omega = \langle m,v\rangle_\Omega \qquad \forall v\in V,
\]
and the resulting UFL representation becomes both the modeling language and the differentiation substrate in automated adjoint frameworks [1302.3894]. 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
\[
p_\epsilon(u\mid z)\propto \exp\left\{ -\frac{1}{2\epsilon^2} (\mathbf{L}[z]\mathbf{u}-\mathbf{f})^\top \mathbf{A}^{-1} (\mathbf{L}[z]\mathbf{u}-\mathbf{f}) \right\},
\]
so that the PDE acts as a weak-form decoder regularizer rather than a hard solve [2006.15641].

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 \(L^1\)-regularized state models, the subgradient relation
\[
\|v\|_1\ge \|u\|_1+\langle v-u,p(u)\rangle
\]
and the pointwise form
\[
p(u)= \begin{cases} \operatorname{sign}(u), & \text{if } |u|>0,\\[4pt] q,\quad |q|\le 1, & \text{if } u=0 \end{cases}
\]
create exact zero regions, compact support, and finite speed of propagation [1311.5850]. In one dimension with \(M=1\), the elliptic and parabolic inclusions become
\[
u_{xx}=-f+\gamma p(u), \qquad
u_t-u_{xx}=f-\gamma p(u),
\]
and support estimates such as
\[
|\operatorname{supp}(u)|\le \gamma^{-1}\int_{\operatorname{supp}(u)} |f|\,dx
\]
quantify the compression effect [1311.5850]. The same paper rewrites the divisible sandpile problem by replacing the positivity constraint with an \(L^1\) term,
\[
\min_u \int \frac12 |\nabla u|^2+|u|-uf\,dx,
\]
whose Euler–Lagrange equation is
\[
\Delta u = p(u)-f.
\]
This replaces a hard state constraint by a nonsmooth variational regularizer [1311.5850].

Total variation regularization occupies a parallel but distinct role. For piecewise constant controls on a triangulation \(T\),
\[
TV(u,\Omega)=\frac12\sum_{i,j}|u(T_i)-u(T_j)|\,\mathcal H^{d-1}(\partial T_i\cap\partial T_j),
\]
and the extremal points of the corresponding zero-mean TV unit ball are normalized characteristic functions of simple sets [2310.19777]. 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 [2310.19777]. The same work proves that \(P_0\) discretizations do not in general approximate isotropic TV; instead they converge to an anisotropic functional
\[
TV_\varphi(u) = \int_\Omega \varphi\!\left(\frac{Du}{|Du|}\right)\,d|Du|,
\]
with anisotropy induced by the fine-scale mesh geometry [2310.19777].

State-dependent inverse problems introduce another regularization layer. The canonical reduced problem is
\[
\hat p=\operatorname*{argmin}_{p\in P} \left\{ \frac12\|O(u(p))-y_{\mathrm{obs}}\|_Y^2 + R(p) \right\},
\]
where the constitutive law \(p(\zeta)\) is identifiable only on the attained state interval
\[
I=[\zeta_\alpha,\zeta_\beta].
\]
Outside \(I\), recovery is regularization- or model-class-driven [2601.10920]. The survey emphasizes both \(L_2\) and \(H^1\) regularization,
\[
R(p)=\frac{\lambda_1}{2}\|p-\bar p\|_{L_2(D)}^2,
\qquad
R(p)=\frac{\lambda_2}{2}\|p-\bar p\|_{\dot H^1(D)}^2,
\]
as well as Sobolev gradients obtained from
\[
\nabla_p^{H^1}J - \ell^2 \partial_s^2 \nabla_p^{H^1}J = \nabla_p^{L_2}J
\]
on the state interval \(D\) [2601.10920]. 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
\[
\mathcal L(u,m,\lambda) \equiv J(u,m) + \lambda^T F(u,m),
\]
with stationarity conditions
\[
\frac{\partial J}{\partial u} + \lambda^T \frac{\partial F}{\partial u}=0, \qquad
\frac{\partial J}{\partial m} + \lambda^T \frac{\partial F}{\partial m}=0, \qquad
F(u,m)=0
\]
[1302.3894]. 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 \(\nabla \hat J(m)\) via the adjoint, and updates \(m\) using a numerical optimizer [1302.3894].

Transport-dominated models make this reduced-space structure especially explicit. For image registration by stationary advection, the reduced gradient is
\[
g(v)=\alpha \mathcal L^\ast \mathcal L\, v + \int_0^1 \lambda \nabla m\,dt,
\]
where \(m\) solves the state equation
\[
\partial_t m + v\cdot \nabla m = 0, \qquad m(\cdot,0)=m_0,
\]
and \(\lambda\) solves the backward adjoint
\[
-\partial_t \lambda - v\cdot \nabla \lambda = 0, \qquad \lambda(\cdot,1)=m(\cdot,1)-m_1
\]
[2510.08782]. 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 \(5\times\) faster than Newton–Krylov baselines in the reported experiments [2510.08782].

For diffeomorphic registration in LDDMM, the control may instead be the initial velocity \(v_0\), with the whole trajectory generated by the EPDiff equation
\[
\partial_t v = -ad_v^\dagger v.
\]
The objective
\[
E(v_0) = \frac{1}{2} \langle L v_0, v_0 \rangle_{L^2} + \frac{1}{\sigma^2}\|m(1)-I_1\|_{L^2}^2
\]
is then constrained simultaneously by EPDiff and image or deformation transport equations [1807.11560]. Gauss–Newton–Krylov optimization acts in the reduced variable \(v_0\), while band-limited vector fields reduce memory and runtime enough to make the geodesic PDE-constrained formulation practical [1807.11560].

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

| Model class | Constraint handling | Characteristic computation |
|---|---|---|
| Classical optimal control [1302.3894] | Hard PDE constraint, reduced functional | Automated discrete adjoints in UFL/FEniCS |
| TV-regularized control [2310.19777] | Hard PDE constraint through \(K\), nonsmooth TV objective | FC-GCG with graph-cut atom oracle |
| \(L^1\)-regularized state models [1311.5850] | PDE as Euler–Lagrange inclusion or gradient flow | Proximal shrinkage, Douglas–Rachford, FFT-based solves |
| Transport-constrained registration [2510.08782] | Hard transport PDE, reduced-space optimization | Preconditioned gradient descent accelerated by GA-NGMRES |
| Momentum-constrained LDDMM [1807.11560] | Hard EPDiff plus transport constraints | Inexact Gauss–Newton–Krylov in band-limited space |
| Weak-form latent inference [2006.15641] | 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
\[
u_\theta = \mathcal{N}_\theta(I), \qquad u(x)\in(0,1),
\]
is optimized with the composite objective
\[
\mathcal{L} = \mathcal{L}_{\text{Dice}} + \mathcal{L}_{\text{BCE}} + \lambda_{\text{RD}} \mathcal{L}_{\text{RD}} + \lambda_{\text{PF}} \mathcal{L}_{\text{PF}},
\]
where the reaction–diffusion residual is
\[
r_{\text{RD}}(u)= D \nabla^2 u + u(1-u)(u-a),
\qquad
\mathcal{L}_{\text{RD}} = \int_\Omega r_{\text{RD}}(u)^2\,dx,
\]
and the phase-field energy is
\[
E_{\text{PF}}(u)= \int_\Omega \left( \frac{\varepsilon}{2} |\nabla u|^2 + \frac{1}{\varepsilon} u^2(1-u)^2 \right) dx
\]
[2602.01069]. 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 \(0.3602\) to \(0.7681\) at 10% training data [2602.01069].

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
\[
\mathcal{L}_{\text{PDE}} = \frac{1}{N_c}\sum_{j=1}^{N_c} R_j(\mathbf f)^2,
\]
with measured outputs
\[
f_k(\boldsymbol{\phi}) = \braket{\psi(\boldsymbol{\phi})|\sigma_Z^{(k)}|\psi(\boldsymbol{\phi})}.
\]
For finite-difference discretizations of a \(p\)-th order PDE, each residual depends on at most \(2p+1\) neighboring outputs, so the total loss is a sum of local cost terms [2604.09957]. 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” [2604.09957].

A second quantum direction addresses hardware noise in variational PDE solvers. There the soft constrained objective
\[
\mathcal{L}(\boldsymbol{\theta}) = \mathcal{L}_{\mathrm{data}}(\boldsymbol{\theta}) + \lambda\,\mathcal{L}_{\mathrm{phys}}(\boldsymbol{\theta}),
\qquad
\mathcal{L}_{\mathrm{phys}}(\boldsymbol{\theta}) = \frac{1}{N_c}\sum_{i=1}^{N_c} \left| \mathcal{N}[u_{\boldsymbol{\theta}}(x_i,t_i)] \right|^2
\]
is studied under depolarizing, amplitude damping, and bit-flip noise. The reported result is that zero-noise extrapolation reduces absolute error by \(82\)–\(96\%\) at \(p=0.001\), while constrained circuits retain \(25\)–\(47\%\) higher fidelity than unconstrained ones at \(p=0.01\) [2604.10099]. 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
\[
\mathcal{L}_{\text{HQC}} = \mathcal{L}_{\text{data}} + \lambda_{\text{SV}}\mathcal{L}_{\text{SV}} + \lambda_{\text{M}}\mathcal{L}_{\text{Manning}}
\]
combines focal classification loss with a Saint-Venant continuity residual
\[
\mathcal{L}_{\text{SV}} = \frac{1}{N_c}\sum_{j=1}^{N_c} \left\| \frac{\partial \hat{A}_j}{\partial t} + \frac{\partial \hat{Q}_j}{\partial x} - \hat{q}_{l,j} \right\|^2
\]
and Manning consistency penalty
\[
\mathcal{L}_{\text{Manning}} = \frac{1}{N_c}\sum_{j=1}^{N_c} \left\| \hat{Q}_j - \frac{1}{n_j}\hat{A}_j\hat{R}_{h,j}^{2/3}\hat{S}_{f,j}^{1/2} \right\|^2.
\]
The reported model converges in about \(3\times\) fewer training epochs and uses about \(44\%\) fewer trainable parameters than the classical PINN baseline [2604.09374].

## 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
\[
\nabla_p J(\zeta_0) = \int_{\Omega\times(0,T)} \delta\bigl(\zeta(x,t)-\zeta_0\bigr)\, \nabla u^*(x,t)\cdot \nabla u(x,t)\,dx\,dt,
\]
so only the interval of state values actually attained by the dynamics is identifiable [2601.10920]. If a perturbation \(\delta p\) is supported outside the attained range, then the first variation vanishes [2601.10920]. This is a structural limitation, not a numerical artifact.

Supremal formulations expose a different structural regime. For Navier–Stokes-constrained data assimilation, the finite-\(p\) objective
\[
E_p(u,p,y) := (1-\lambda)\,\|K(\cdot,\cdot,u,\nabla u,p)\|_{\dot L^p(\Omega_T)} + \lambda\,\|y\|_{\dot L^p(\Omega_T)}
\]
converges, as \(p\to\infty\), to the worst-case functional
\[
E_\infty(u,p,y) := (1-\lambda)\, \|K(\cdot,\cdot,u,\nabla u,p)\|_{L^\infty(\Omega_T)} + \lambda\, \|y\|_{L^\infty(\Omega_T)}.
\]
The resulting \(L^\infty\) optimality system involves measure coefficients rather than a classical Euler–Lagrange PDE and concentrates on the sets where the maximal mismatch is attained [2105.06547]. 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
\[
\mathcal{L}[u] := \|\mathcal{N}[u]-f\|_{H^a(\Omega)}^2 + \|u-g\|_{H^b(\partial\Omega)}^2,
\]
and one studies computability of the minimizer through finite-dimensional approximation and discrete gradient flow [2510.21290]. Under convexity, coercivity, QGC, and Lipschitz gradient assumptions, the error decomposes as
\[
\|u^*-\hat u_{\theta^j}\|_{H(\Omega)}^2 \lesssim \epsilon_{\mathrm{app}} + \epsilon_{\mathrm{int}} + \epsilon_{\mathrm{opt}},
\]
with the optimization error decaying exponentially in the iteration index under the stated assumptions [2510.21290]. 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 [2510.21290]. 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 [2601.10920]. 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 [2602.01069]. The quantum papers likewise emphasize small-scale simulated validation, soft rather than hard enforcement, and unresolved scaling questions [2604.09957] [2604.10099] [2604.09374]. 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.

Source: https://www.emergentmind.com/topics/pde-constrained-variational-models