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

Motivation

• Many PDE solutions must satisfy physical or economic inequalities that equality conditioning cannot enforce, leaving posterior mass on inadmissible states.Examples include nonnegative concentrations, monotone transport solutions, and arbitrage-free option surfaces.Source: §1 · p. 1

Research gap

• Existing physics-informed GP solvers and constrained GPs handle PDE equalities and inequalities separately, leaving their combination in one posterior unaddressed.The passage contrasts operator-conditioned PDE solvers with inequality-constrained data-fit priors.Source: §1 · p. 2

Contribution

• The proposed two-stage construction first conditions on PDE, boundary, and data equalities, then applies linear inequality information at virtual locations.The resulting posterior is non-Gaussian, with finite-dimensional marginal means and covariances computed through truncated-Gaussian moments.Source: §1 · p. 2

Contributions

• Sequential and joint conditioning yield the same Borel probability measure on the function space, while the solver retains operator cross-covariances.The paper also evaluates changes in posterior-mean accuracy and constraint diagnostics across observation regimes.Source: §1 · p. 2

Probabilistic numerics and GP PDE solvers

• Physics-informed GP solvers condition on linear PDE operators to obtain tractable solution posteriors; Pförtner et al. (2022) formalize the equality stage used here.Their operator framework generalizes weighted residual methods and unifies collocation, finite-element, and spectral solvers.Source: §2 · p. 2

Encoding linear constraints in the GP prior

• Equality-constrained GP priors can satisfy linear systems by construction, but these approaches target equalities rather than inequalities and build priors for particular operators.The proposed solver retains a generic LinPDE-GP recipe and adds inequalities.Source: §2 · p. 2

Inequality and shape-constrained GPs

• Virtual-location methods impose GP inequalities using truncated-Gaussian moments, sampling, or finite-basis approximations; Agrell (2019) extends the approach to arbitrary linear operators.The paper adopts Agrell’s operator-based construction, while noting that the family scales poorly with dimension.Source: §2 · p. 2

Physics-informed neural networks

• PINNs impose governing equations through soft residual penalties, whereas the proposed method conditions on an inequality event at noisy virtual points.The cited discussion characterizes PINN uncertainty as approximate rather than a conditional law.Source: §2 · p. 3

Arbitrage-free option-price surfaces

• Existing GP option-surface approaches combine quote fitting or calibration with shape constraints, while the proposed equality stage conditions on the Dupire operator itself.The surface’s admissibility conditions include price bounds, maturity monotonicity, and strike convexity.Source: §2 · p. 3

Positioning

• The paper links operator-level PDE conditioning with linear inequalities, a combination its authors say prior work had not put in one posterior.Proposition 1 makes the LinPDE-GP posterior a valid prior for Agrell’s update and establishes equivalence with joint conditioning.Source: §2 · p. 3

Equality information

• The prior is a Gaussian process whose sample paths lie in a separable reproducing kernel Banach space, making point evaluation a bounded linear functional.The path map is assumed to be a function-space-valued Gaussian random element.Source: §3 · p. 3
• A bounded linear operator collects finite-dimensional equality information, and independent Gaussian noise permits noise-free or noisy coordinates.Conditioning on the realized equalities leaves a Gaussian process under the stated assumptions.Source: §3 · p. 3

Finite-dimensional inequality information

• The inequality operator stacks bounded linear maps across constraint families and virtual locations into a finite vector of scalar quantities.The vector dimension counts the pooled virtual variables; their bounds form a componentwise box.Source: §3 · p. 3
• Independent positive-definite virtual noise makes the box restriction soft at selected locations rather than enforcing inequalities throughout the domain.The noisy virtual vector is defined by applying the inequality operator to the solution and adding virtual noise.Source: §3 · p. 4

Target law and diagnostic

• Conditioning on the equality event and the virtual box event gives a generally non-Gaussian target law; P(C | E) is the Gaussian box probability before the inequality update.The diagnostic depends on the operator, bounds, virtual noise, locations, and conditional mean and covariance, so comparisons are most defensible within a fixed design.Source: §3 · p. 4

Theoretical Framework

• Joint conditioning on the linear equalities and inequality event gives the same Borel probability measure on U as sequentially updating the LinPDE-GP posterior.The agreement holds for every finite-dimensional linear image and is lifted to the full path-space law.Source: §4 · pp. 4–5
• The constrained finite-dimensional marginals are Gaussian conditionals mixed over a truncated Gaussian virtual-observation variable.This compound-Gaussian representation is the basis for computing posterior summaries.Source: §4 · pp. 4–5
• The proof fixes a regular conditional Gaussian law on U for each observed value y, avoiding almost-sure qualifications on y.The construction uses a version of the equality-conditioned law as a measure on U.Source: §4 · p. 4
• Agreement of every finite-dimensional linear image extends to all Borel sets of U because the cylinder sets generate B(U).The proof uses separability of U and the Dynkin π–λ theorem.Source: §4 · p. 5
• Strictly positive virtual-noise covariance and componentwise nonempty bounds ensure the inequality event has positive probability under the equality posterior.The proposition assumes Σv ≻ 0 and ai < bi componentwise.Source: §4 · p. 5

5 Proposed Method

• The solver applies an equality update followed by an inequality update at virtual points, with posterior summaries computed from truncated-Gaussian moments.The constrained law is non-Gaussian and its finite-dimensional marginals use Agrell’s compound-Gaussian construction.Source: §5 · p. 6
• The constraints use bounded identity or first-order partial-derivative operators; second-order constraints are left to future work.Nonnegativity uses one family, Poisson monotonicity uses two families, and option-surface conditions use three families on separate grids.Source: §5 · p. 6
• Operator cross-covariances are available in closed form, and multiple constraint families produce a block-structured virtual-observation covariance.The off-diagonal blocks are cross-family covariances.Source: §5 · p. 6
• The predictive covariance conditions on the virtual constraint vector through a cross-covariance correction.The covariance expression is given for the constrained prediction at u* conditional on Z.Source: §5 · p. 6
• Truncated-Gaussian moments are evaluated by randomized quasi-Monte Carlo, so the posterior summaries are subject to numerical-integration error.The reported posterior summaries include the mean and a 1.96-standard-deviation band; sample paths use an exact minimax-tilting sampler.Source: §5 · p. 6
• The stage-one mean is handled by shifting the inequality bounds at the virtual points.The implementation gives this as an equivalent transformation.Source: §5 · p. 6

6.1 Two-dimensional Poisson equation: Monotonicity improves reconstruction with incomplete data

• With incomplete boundary data, monotonicity reduced full-grid Poisson RMSE from 0.33 to 0.12 (63%) and average posterior band width from 0.93 to 0.52.Absolute error at the unobserved corner (1, 1) fell from 1.07 to 0.04, but residual uncertainty remained largest there and overshoot reached 0.24 along unobserved edges. Monotonicity alone had RMSE 0.70, showing that PDE and boundary observations were needed to anchor magnitude.Source: Table 1 · p. 8, §6.1 · pp. 8–9
• For Poisson data with 16 noisy observations, monotonicity reduced mean RMSE from 0.095 to 0.071 and average band width from 0.39 to 0.27.At σobs = 0.1, RMSE reduction was 30% with 10 observations and 13% with 40; at nobs = 16, it ranged from 6% at σobs = 0.02 to 26% at σobs = 0.2. Exact observations yielded only a 0.5% reduction, with P(C | E) = 0.998 before the update.Source: §6.1 · p. 10

6.2 Arbitrage-constrained option-price surfaces

• For the synthetic option surface at 16 quotes and σobs = 0.01, inequalities reduced mean held-out RMSE from 0.047 to 0.019 (59%).Across realizations, upper spread-bound violations fell from 21% to 15%, calendar violations from 15% to 13%, and un imposed butterfly-convexity violations from 19% to 17%; at 10 versus 60 quotes, RMSE reductions were 56% and 42%.Source: §6.2 · p. 11, Figure 5 · p. 12
• At 20 retained S&P 500 quotes, inequalities reduced held-out RMSE from 0.032 to 0.025 (22%) and calendar violations from 10% to 1%.Across retained-quote counts, mean held-out RMSE fell; the reduction rose from 9% with 300 quotes to 23% with 12. With all 3658 quotes, PDE conditioning reduced RMSE from 0.075 to 0.027, and inequalities left it at 0.027.Source: §6.2 · pp. 11–12, Figure 6 · p. 13

6.3 Nonnegative reconstruction of an air-quality field

• Across five synthetic air-quality noise realizations, nonnegativity reduced the Q0 sub-zero-tail proxy from 0.160 to 0.107 (33%), while mean held-out RMSE increased from 0.076 to 0.088.Average band width narrowed from 0.35 to 0.30; the virtual constraint event probability before updating was P(C | E) = 0.09.Source: §6.3 · p. 13
• In the synthetic air-quality noise sweep, increasing σobs from 0.05 to 0.3 raised Q0 reduction from 23% to 56%, while held-out RMSE reduction shifted from 7% to −43%.At σobs = 0.1, reducing monitors from 28 to 8 raised Q0 reduction from 28% to 41%, while RMSE reduction fell from 32% to −19%; stronger tail-proxy reductions could accompany poorer mean predictions.Source: §6.3 · p. 14
• The real-data analysis used four daytime OpenAQ PM2.5 and NO2 snapshots from Maryland and New York, with 90 evaluations across seeded splits and sampled monitor counts.Each snapshot used its own ERA5-derived wind direction and a source-free transport model with S = 0; monitor counts were sampled before stations were withheld for scoring.Source: §6.3 · p. 14
• At fixed nobs = 14, tail-proxy gains grew with observation noise while RMSE reduction became negative.The curves compare gains over five noise realizations against the PDE-conditioned prediction.Source: Figure 8 · p. 15
• Tail-proxy reductions rose from 67% at 10 monitors to (74 ± 3)% at 6 monitors across four snapshots.Mean P(C | E) decreased from 0.55 to 0.47 over the same range; clipping negative equality-posterior means reduced Q0 much less than adding nonnegativity.Source: §6.3 · p. 14, Figure 9 · p. 15
• Across four snapshots, thinning from ten to six monitors changed held-out RMSE reduction from approximately −1.5% to 2.6%.Nonnegativity consistently reduced the moment-based tail diagnostic, but its effect on held-out mean accuracy depended on the sampled network.Source: §6.3 · p. 14, Figure 9 · p. 15

Discussion and limitations

• Constraint benefits depend on data quality: sparse or noisy observations can benefit, while abundant noise-free Poisson and dense option data gain little.The comparisons support treating inequalities as an additional information source whose value depends on the number and quality of observations.Source: §7 · pp. 15–16
• Accuracy, constraint violations, and dispersion can diverge; nonnegativity sometimes lowered Q0 while increasing held-out RMSE, and narrower bands do not guarantee accuracy.The reported synthetic option sweeps also include settings where RMSE decreased while derivative-violation rates increased.Source: §7 · p. 16
• P(C|E) measures pre-update agreement with the specified bounds, but alone neither identifies why the event is unlikely nor predicts accuracy gains.Its value depends on the equality posterior, constraint operators, bounds, virtual locations, and virtual noise.Source: §7 · p. 16
• Proposition 1 makes the equality posterior a valid prior for the inequality update and establishes equivalence between sequential and joint conditioning.The construction retains the chosen GP prior and Gaussian equality update while inequalities update posterior means and covariances.Source: §7 · p. 16
• Constraints apply to a finite vector of noisy operator evaluations, so they are soft at selected locations and feasibility elsewhere is assessed on the evaluation grid.Virtual-point placement and noise are part of the model specification; truncated-Gaussian moments are evaluated numerically, with computation becoming expensive as virtual dimension grows.Source: §7 · p. 16

Conclusion

• Sequential Gaussian equality conditioning followed by a finite virtual constraint event yields the same law as joint conditioning.The result is stated for Borel probability measures on the solution space.Source: §8 · p. 16
• The constrained solver retains operator cross-covariances and computes posterior moments through truncated-Gaussian calculations.The resulting law is non-Gaussian.Source: §8 · p. 16
• Monotonicity, option-price bounds, and nonnegativity add information beyond PDE observations, but effects on error, constraint diagnostics, and dispersion differ across settings.P(C | E) measures agreement with the specified virtual constraints before the update; predictive benefits are evaluated against reference solutions or held-out data.Source: §8 · pp. 16–17
• With soft constraints at finitely many virtual locations, physical and economic admissibility is incorporated directly into probabilistic PDE inference.The experiments include Poisson, option-price, and air-quality settings.Source: §8 · p. 17

Measure-Theoretic Background

• A regular conditional distribution is a measurable map from observed values to Borel probability measures, and the proof fixes an explicit Gaussian version for every prescribed value.This gives Proposition 1 a definite meaning even at observations with zero probability under the observation law.Source: A Measure-Theoretic Background · p. 17
• Image measures preserve equality in law: random variables with the same law have images under a measurable map with the same law.This is the consequence used in the proof's image-measure steps.Source: A Measure-Theoretic Background · p. 17
• Cylinder sets record finitely many bounded linear measurements, and the π–λ theorem extends agreement on them to the generated Borel σ-algebra.For separable U, the passage states that these cylinder sets generate B(U).Source: A Measure-Theoretic Background · p. 17
• The equality Gram matrix may be singular; the pseudoinverse construction still makes the residual cross-covariance vanish.The proof uses the fact that directions in the Gram matrix's kernel have zero variance.Source: A Measure-Theoretic Background · p. 17

• The air-quality operator is assembled from weighted Laplacian, directional-derivative, and identity terms, yielding closed-form squared-exponential cross-covariances that agreed with nested automatic differentiation within 10^-6.A single merged operator instead falls back to bare first-order partials, which the kernel does not dispatch in closed form.Source: B Implementation Details · p. 17
• Option constraints use price bounds, two-sided call-spread bounds, and calendar monotonicity on a 3×3 virtual grid; butterfly convexity is only a finite-difference diagnostic.Spread, calendar, and butterfly violation rates use valid finite-difference locations, with τ = 10^-4; the spread diagnostic counts only upper-bound violations.Source: B Implementation Details · pp. 17–18
• At nobs = 10 and σobs = 0.01, adding inequalities reduced mean held-out RMSE by 56%, while calendar violations rose from 6.0% to 12.9% and butterfly violations from 9.9% to 16.7%.These changes were averaged over the same 5 noise realizations; the grid-based diagnostics measure changes beyond the finite virtual constraint event.Source: B Implementation Details · p. 18
• The inequality update centers the process at the physics-informed mean and shifts bounds by its linear image, using virtual-observation noise of 10^-4.When shared derivative blocks make the virtual covariance numerically singular, the noise is increased geometrically until factorization succeeds.Source: B Implementation Details · p. 18
• The synthetic air-quality reference field uses a second-order finite-difference solve, whose manufactured-solution test observed convergence rate approximately 2.Source: B Implementation Details · p. 18
• Constrained posterior bands use the mean plus or minus 1.96 marginal standard deviations as a Gaussian-equivalent dispersion summary, not a 95% credible interval.The constrained standard deviation comes from the truncated-Gaussian covariance; randomized quasi-Monte Carlo evaluates the moments, leaving numerical-integration error despite a fixed seed.Source: B Implementation Details · p. 18
• The truncated-Gaussian moment evaluation dominates the inequality update, whose cost depends on virtual constraint dimension d rather than prediction-point count and grows at least quadratically in d.The experiments used d = 32 for Poisson, 27 for options, and 9 for air quality.Source: B Implementation Details · p. 18

• The real-data experiments draw on a single-day ^SPX listing snapshot from Yahoo Finance and pollutant measurements from the OpenAQ v3 API.Source: C Data and Preprocessing · p. 18
• The transport winds are sampled hourly at 10 m from the Open-Meteo historical archive, which serves ECMWF ERA5 reanalysis.Source: C Data and Preprocessing · p. 18
• Snapshot files and fetch and preprocessing scripts are released in the accompanying code repository.Source: C Data and Preprocessing · p. 18

• The ^SPX snapshot contains 18,352 raw quotes across 49 expiries; filtering leaves 11,382 quotes over 30 expiries after bid, spread, and maturity criteria.Quotes with nonpositive bids or relative bid–ask spreads above 25% are dropped, and maturities are restricted to 0.05–1 years.Source: §C.1 · p. 18
• The ^SPX chain is a single-day Yahoo Finance snapshot, and both calls and puts are retained initially to recover forwards and discount rates from paired quotes.Source: §C.1 · p. 18
• For each expiry, put–call-parity regressions on paired call and put strikes estimate rates and forwards; across-expiry medians give spot ≈7.4 × 10^3 index points and rate ≈7%.Prices and strikes are normalized to c = Ccall/S0 and x = ln(K/S0).Source: §C.1 · p. 19
• The constant-coefficient Dupire operator uses average at-the-money implied volatility, while quote noise is set to the median normalized bid–ask half-spread, about 2 × 10^-4.At-the-money quotes have |x| < 0.03; each quote's noise estimate is half its normalized bid–ask spread.Source: §C.1 · p. 19
• The held-out extrapolation block contains 971 retained calls, while 3,658 bulk quotes form the training pool and the score target stays fixed across arms and subsamples.The held-out block uses log-moneyness |x| > ln 1.15 or maturity T > 0.8 years.Source: §C.1 · p. 19
• The sparse-observation sweep retains 12, 20, 35, 60, 120, or 300 quotes, using 20 independently seeded subsets per count and a shared noise level within each subset.The dense control uses all 3,658 bulk quotes at their original mid prices and chain-derived noise.Source: §C.1 · p. 19

• Each snapshot averages hourly OpenAQ readings at each station over a daytime window of up to six hours, matching the steady-state field reconstruction.Station coordinates are projected to a local east–north frame and rescaled to the unit square.Source: §C.2 · p. 19
• The four pooled snapshots differ in day, region, and pollutant; two further live pulls were excluded because they had fewer than six monitors.A London sample with real station coordinates but synthetic concentrations is provided for offline reproducibility and excluded from the pool.Source: §C.2 · p. 19
• Each snapshot's advection direction comes from the vector-mean 10 m wind during local daytime hours, while wind magnitude is normalized away.Source: §C.2 · p. 19
• Four snapshots yield 90 snapshot–count–split evaluations, with sampled monitor counts from 6 to 22 where station availability permits and held-out monitors scored within each split.Five seeded subsets are drawn at each retained count; nmon includes both training and held-out monitors.Source: §C.2 · p. 19
• At each monitor count, pooled ratio-of-means reductions summarize held-out RMSE and the sub-zero-tail diagnostic Q0; for the non-Gaussian posterior, Q0 is a Gaussian moment-based proxy.Bootstrap standard errors use 2,000 resamples of snapshot–split rows, preserving paired equality and constrained results.Source: §C.2 · p. 20
• The all-snapshot curve pools four snapshots at monitor counts 6, 8, and 10, while a Maryland PM2.5 curve pools two snapshots across all six sampled counts.The fixed compositions separate changes with monitor count from changes in contributing snapshots.Source: §C.2 · p. 20