Constrained Probabilistic Numerics: Gaussian Process Solvers for Linear PDEs under Linear Inequality Constraints

Boyuan Deng; Kshitiz Upadhyay; Michael Shields

stat.MLOct 9, 2026

Constrained Probabilistic Numerics: Gaussian Process Solvers for Linear PDEs under Linear Inequality Constraints

Authors: Boyuan Deng, Kshitiz Upadhyay, Michael Shields

Organizations: Department of Civil and Systems Engineering, Johns Hopkins University · Department of Aerospace Engineering and Mechanics, University of Minnesota

Abstract

Gaussian process solvers for linear partial differential equations (PDEs) return posteriors rather than point estimates, but they condition only on equality information: the PDE residual, boundary data, and measurements. Many solutions of interest also obey inequality constraints, including nonnegativity, monotonicity, and absence of arbitrage, which equality conditioning cannot enforce. We prove that conditioning a Gaussian process prior jointly on linear equality observations and a finite-dimensional inequality event is equivalent to conditioning sequentially in terms of Borel measures on the solution space. The physics-informed posterior is therefore a valid prior for the inequality update at virtual points. The combined law is non-Gaussian yet computable through truncated-Gaussian moments, and neither the PDE conditioning nor the inequality update requires modification. We evaluated the construction on a Poisson benchmark, option-price surfaces, and air-quality fields. Under incomplete boundary data, a monotonicity constraint lowered the error of a two-dimensional Poisson reconstruction by $63\%$ . On a real S&P 500 option chain thinned to $20$ retained quotes, no-arbitrage constraints reduced held-out extrapolation error by $22\%$ and cut the calendar-arbitrage rate from $10\%$ to $1\%$ . On four real air-quality snapshots, nonnegativity reduced a Gaussian moment-based sub-zero-tail proxy by $74\%$ at the sparsest monitor count. Changes in mean accuracy and constraint diagnostics differed across data regimes. The probability of the noisy virtual constraint event under the equality posterior measures agreement with the specified bounds before the update. The construction turns admissibility from a post-processing step into conditioning information for probabilistic PDE solvers.

Figures & tables

Y=L_{\rm eq}u+\varepsilon_{\rm eq}=y,\qquad\varepsilon_{\rm eq}\sim\mathcal{N}(0,\Sigma_{\rm eq}),

Algorithm 1 LinPDE-GP with linear inequality constraints

Figure 1: Monotonicity improves Poisson reconstruction with incomplete boundary data. (a) Diagonal posterior means for equality-only ( $u|E$ ), inequality-only ( $u|C$ ), and combined conditioning ( $u|(E,C)$ ), compared with the truth ( $u_{\star}$ ). Shading shows $m\pm 1.96s$ for the equality-only and combined posteriors. (b) Manufactured solution with boundary observations (squares), PDE collocation points (open circles), and monotonicity virtual points (filled circles). (c–e) Posterior means conditioned on the PDE and partial boundary data, monotonicity alone, and both, respectively. Panels (b–e) share a color scale, with values outside $[0,1]$ shown in the endpoint colors.

Figure 2: Signed reconstruction errors and posterior uncertainty for the Poisson problem with incomplete boundary data. Equality data $E$ comprise the PDE and partial boundary data; $C$ denotes the monotonicity constraints. (a,b) Signed errors $m-u_{\star}$ of the equality-only and combined posterior means, respectively. (c,d) Corresponding posterior half-widths $1.96s$ , where $s$ is the posterior standard deviation, showing uncertainty magnitude alone. Each pair shares a color scale. Uncertainty is largest at the unobserved corner $(1,1)$ and is reduced by conditioning on monotonicity.

	Equality only	With monotonicity	Reduction
	$u\mid E$	$u\mid(E,C)$
RMSE	0.328	0.122	$63\%$
Maximum absolute error	1.067	0.242	$77\%$
Mean full band width	0.930	0.516	$45\%$
Absolute error at $(1,1)$	1.067	0.043	$96\%$

Table 1: Poisson reconstruction errors and posterior dispersion with partial boundary data. Both predictions use the same prior, PDE collocation points, and boundary observations. Band width refers to the full width of

m\pm 1.96s

, averaged over the evaluation grid.

Figure 3: Poisson reconstruction with PDE and partial boundary data (

E

), monotonicity alone (

C

), and both

(E,C)

. Bars show full-grid RMSE, mean full width of the

m\pm 1.96s

bands, and absolute error at

(1,1)

Figure 4: Effect of observation noise magnitude

\sigma_{\mathrm{obs}}

and number of observations

n_{\mathrm{obs}}

on the benefit of adding monotonicity constraints. (a) Noise sweep at

n_{\mathrm{obs}}=16

. (b) Observation-count sweep at

\sigma_{\mathrm{obs}}=0.1

. Solid curves show full-grid RMSE against the true solution

u_{\star}

, averaged over 5 noise realizations. Ribbons show

\pm 1

standard deviation across realizations. Green shading gives the absolute RMSE reduction. Dashed curves report the relative reduction from Eq. ( 20 ) on the right-hand axes. Monotonicity gives larger relative improvements with noisy observations and at the lower observation counts.

Figure 5: Price and derivative constraints improve extrapolation of a synthetic option-price surface. Derivatives in panels (b)–(f) are computed from the posterior means by central differences. (a) Black–Scholes reference surface with training quotes, held-out quotes, and virtual points. (b,c) Log-moneyness derivatives $\partial_{x}m$ before and after adding inequality constraints to the PDE-conditioned model. Positive values violate the upper spread bound. The dashed blue contour marks the lower bound. (d) Maximum $\partial_{x}m$ and minimum $\partial_{T}m$ over maturities, plotted against log-moneyness. Adding inequality constraints reduces both the magnitude and spatial extent of violations. Red shaded regions indicate infeasible values. (e,f) Maturity derivatives $\partial_{T}m$ before and after adding inequality constraints. Negative values shown in red violate calendar monotonicity. (g) Held-out RMSE for quotes alone, quotes with the PDE, quotes with inequality constraints alone, and quotes with both. (h) Constraint-violation rates of the posterior means, including butterfly convexity, which was not imposed. Panels (b)–(h) show one noise realization. (i) Mean held-out RMSE versus quote count at $\sigma_{\mathrm{obs}}=0.01$ , with $\pm 1$ standard-deviation ribbons across 5 realizations. Green shading gives the absolute RMSE reduction; the annotation gives the relative reduction.

Figure 6: Effects of inequality constraints on an S&P 500 option chain with experimentally thinned and perturbed training quotes. (a) Mean held-out RMSE across 20 subsets at each retained-quote count. Gray shading shows the

25

th–

75

th percentiles. Adding the inequality constraints reduces mean error throughout the thinning experiment. The separate dense control (square gray point) uses all

3658

training quotes at their bid–ask-derived noise level, where the two RMSE values coincide. (b) Calendar-violation rates of the posterior mean surfaces, averaged over the same subsets, decrease after adding the inequality constraints at every tested quote count.

Figure 7: Nonnegativity reduces negative posterior means in a synthetic plume reconstruction. Results use $n_{\mathrm{obs}}=14$ , $\sigma_{\mathrm{obs}}=0.1$ , and one noise realization. (a) Finite-difference reference field with monitor, held-out, and virtual locations. (b) Posterior mean conditioned on monitors and the PDE, with negative values overlaid in red. (c) Adding nonnegativity reduces the negative region. The dashed white contour marks its boundary before the update. Table 2 reports metrics for this realization.

Conditioning	Held-out RMSE	Fraction $m<0$	Band width
Monitors	0.086	0.102	1.602
Monitors + PDE	0.142	0.120	0.355
Monitors + nonnegativity	0.314	0.078	1.282
Monitors + PDE + nonnegativity	0.075	0.084	0.305

Table 2: Synthetic air-quality reconstruction for the noise realization shown in Figure 7 . Negative-mean fractions and full

m\pm 1.96s

band widths are evaluated over the grid. Bold values mark the minimum in each column.

Figure 8: Observation noise separates improvements in the sub-zero tail diagnostic from improvements in mean accuracy. Curves show percentage gains in the sub-zero-tail proxy

Q_{0}

and held-out RMSE after adding the nonnegativity constraint, for fixed

n_{\mathrm{obs}}=14

. Gains compare means over 5 noise realizations using Eq. ( 20 ). The gains in

Q_{0}

grow with noise, whereas the RMSE reduction becomes negative, indicating reduced accuracy. The horizontal line marks no change relative to the PDE-conditioned prediction.

Figure 9: Effects of nonnegativity constraints on real OpenAQ reconstructions. (a) Gains in the Gaussian moment-based sub-zero-tail proxy

Q_{0}

. The dashed baseline clips the equality-posterior mean at zero and retains its variance, using the two snapshots that span the full monitor-count range. (b) Held-out RMSE reduction, which can be positive or negative. The dark blue curves pool all 4 snapshots over their common monitor counts. Light blue curves pool the two Maryland

\mathrm{PM}_{2.5}

snapshots over the full range. Monitor counts decrease from left to right and include stations subsequently held out for scoring. Gains compare pooled means. Shading denotes bootstrap standard errors from 2,000 resamples of snapshot–split pairs at each count.

Appendix figures & tables1 asset

Supplementary material from the paper’s appendix.

Appendix

Region	Pollutant	Date	Monitors
Maryland (state)	$\mathrm{PM}_{2.5}$	2021-07-15	$22$
Maryland (state)	$\mathrm{PM}_{2.5}$	2023-07-20	$31$
New York City	$\mathrm{PM}_{2.5}$	2023-07-20	$12$
Maryland (state)	$\mathrm{NO}_{2}$	2024-06-20	$13$

Appendix

Table 3: OpenAQ daytime snapshots in the pooled sparsity study. The monitor count is the number of live stations returned inside the bounding box after per-station aggregation.

Explore similar work

Efficient Message Passing for Partial Differential Equation Priors
Sep 26, 2026Anna Kazachkova, Leonhard Hennicke, Rainer Schlosser +1PDE Solving Physics-Informed ML
Hasso Plattner Institute University of Potsdam Potsdam, Germany
Multi-output Gaussian process prediction of physical fields under linear equality constraints
Aug 26, 2026Mahamat Hamdan Nassouradine, Clément Gauchy, Pierre-Emmanuel Angeli +1Surrogate Modeling Gaussian Process Regression
Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation, Gif-sur-Yvette, 91191, France · Univ Rennes, Ensai, CNRS, CREST - UMR 9194, Rennes, F-35000, France · Université Paris-Saclay, CEA, Service de Thermohydraulique et de Mécanique des Fluides, Gif-sur-Yvette, 91191, France
Error-Conditioned Neural Solvers
Jun 25, 2026Haina Jiang, Liam Wang, Peng-Chen Chen +4Neural PDE Solvers PDE Surrogate Modeling
University of Michigan · KAIST AI · Los Alamos National Laboratory

Sep 26, 2026, cs.LGJ/K move · Enter open · S save

Efficient Message Passing for Partial Differential Equation Priors

Anna Kazachkova, Leonhard Hennicke, Rainer Schlosser, Ralf Herbrich

Hasso Plattner Institute University of Potsdam Potsdam, Germany

PDE Solving Physics-Informed ML Bayesian Inference

Open paper AskPDF arXiv

Prior information for real-world physical quantities is most elegantly expressed via partial differential equations (PDEs). In this paper, we propose a novel way to solve PDEs using probabilistic inference on a factor graph. In general, factor graphs provide a natural way to encode prior knowledge into a model as explicit factors; here, this knowledge is provided by a governing PDE, which narrows the solution space, while observed data further shape the posterior over the parameters. The approximate parameter posterior is inferred using message passing based on moment matching, without posterior sampling or global gradient-based optimization. We demonstrate our approach on the first-order advection and the second-order semi-linear Fisher-KPP equations, where it achieves predictive accuracy comparable to a standard baseline while providing structured predictive uncertainty. Moreover, the inferred posterior marginal means and uncertainty structure match more closely those obtained using Hamiltonian Monte Carlo than the evaluated mean-field variational inference baseline, while requiring up to 10x less training time in our experiments, with inference speed comparable to variational inference.

Abstract

Explore similar work

Efficient Message Passing for Partial Differential Equation Priors

Multi-output Gaussian process prediction of physical fields under linear equality constraints

Error-Conditioned Neural Solvers