Boundary-Constrained Gaussian Processes

• Boundary-Constrained Gaussian Processes are Gaussian process models that embed exact boundary and shape constraints to reflect underlying physical or geometric properties.
• They utilize methods like eigenfunction expansion, SPDE formulations, and algebraic parametrization to ensure all sample paths adhere to prescribed constraints.
• Applications include surrogate modeling, uncertainty quantification, and inverse problems, leading to improved convergence and predictive performance over standard GPs.

Boundary-Constrained Gaussian Processes (BCGPs) are Gaussian process (GP) models designed so that all sample paths exactly—or in a specified sense—satisfy boundary conditions and/or shape constraints (such as Dirichlet, Neumann, Robin, positivity, monotonicity, or convexity) over the input domain. BCGPs combine probabilistic modeling with rigorous incorporation of prior or physical boundary knowledge, making them fundamentally different from unconstrained GPs and providing increased reliability, predictive accuracy, and deeper connections to the underlying physics or structure of the modeled phenomenon.

1. Core Methodologies for Enforcing Boundary Constraints

Numerous BCGP frameworks have been developed, all sharing the aim of embedding information about boundaries into the GP prior or posterior. The principal methods are:

• Covariance Modification via Eigenfunction Expansion: By expanding the covariance operator in the basis of Laplace (or more general differential operator) eigenfunctions satisfying prescribed boundary conditions, the resulting GP inherits the desired behavior globally. For a domain $\Omega$ and eigenpairs $(\lambda_j, \phi_j)$ of $-\Delta$ with, e.g., Dirichlet ($\phi_j|_{\partial\Omega}=0$) or Neumann ($$\frac{\partial\phi_j}{\partial n}|_{\partial\Omega}=0$$) conditions, the covariance is truncated as

$$k(x,x') \approx \sum_{j=1}^m s(\lambda_j) \phi_j(x)\phi_j(x').$$

This structure ensures the covariances and hence all sample paths respect the boundary constraints (Solin et al., 2019, Swiler et al., 2020, Riutort-Mayol et al., 2020, Jones et al., 2022).

• SPDE and Path Integral Formulations: For Matérn-type GPs, the covariance is interpreted as the Green’s function of the associated stochastic partial differential equation (SPDE) equipped with boundary conditions. The BdryMatérn GP, for example, modifies the SPDE to incorporate Dirichlet, Neumann, or Robin boundaries and uses the Feynman–Kac formula to express the covariance as an expected value over Brownian paths conditioned to satisfy the boundary (Ding et al., 12 Jul 2025).
• Spectral Projection on Constrained Subspaces: Physics-informed BCGPs for, e.g., incompressible flows on domains with immersed profiles, use spectral projection to “zero-out” components of a base kernel on the boundary, ensuring that function values there are constant (homogeneous Dirichlet) or match prescribed behaviors (Padilla-Segarra et al., 23 Jul 2025).
• Algorithmic Parametrization via Computer Algebra: For linear PDEs and boundary value problems with constant coefficients, BCGPs may be constructed by computing a parametrization of the solution set of the PDE and the admissible functions for the boundary constraints using Gröbner bases; the intersection yields a pushforward map from a white noise GP to the set of physically admissible functions (Lange-Hegermann, 2020, Huang et al., 25 Nov 2024).
• Finite-Dimensional Approximation and Convex Restriction: For inequality constraints (boundedness, monotonicity, convexity), the GP is approximated by a finite linear basis expansion wherein the infinite collection of constraints in function space is reduced to finite (often convex) constraints on the expansion coefficients. Posterior inference employs (truncated) multivariate Gaussian sampling under these restrictions (Maatouk, 2017, Agrell, 2019, Zhang et al., 2018).

2. Classes of Boundary and Shape Constraints

BCGPs have been developed for several boundary and shape constraint classes:

Type	Imposition Mechanism	Typical Domains
Dirichlet	Function vanishes (or equals prescribed value)	Regular and irregular, all dim
Neumann	Normal derivative vanishes	All domains, needs normal vectors
Robin	Linear combination of value and derivative	All domains, physics-motivated
Periodic	Covariance function wrapped with period $N$	Tori, circles, spatial lattices
Positivity, boundedness	Truncation/projection or convex constraints	Any
Monotonicity	Convex constraints on derivatives coefficients	1D, sometimes multidimensional
PDE-constrained	Parametrization or spectral kernel construction	Physics (heat, wave, div-free fields)

The adoption of one or more of these constraints is problem-specific and typically tied to known physical, geometric, or statistical properties.

3. Mathematical Formulation and Covariance Kernel Construction

Several canonical constructions appear across the literature:

• Path-integral Matérn kernels (BdryMatérn GP): The covariance is given, for Dirichlet boundaries, via

$$k_{\nu+2,\mathcal{B}}(x,x') = k_{\nu+2}(x,x') - \mathbb{E}[k_{\nu+2}(x, B_\tau) e^{-\kappa^2\tau}] - \mathbb{E}[k_{\nu+2}(B_\gamma, x') e^{-\kappa^2\gamma}] + \mathbb{E}[k_{\nu+2}(B_\tau, B_\gamma) e^{-\kappa^2(\tau+\gamma)}],$$

where $B_t$ are Brownian motions and $\tau$, $\gamma$ their hitting times of $\partial\mathcal{X}$ (Ding et al., 12 Jul 2025).

• Spectral kernel expansion with constrained basis: The expansion

$$k(x, x') = \sum_{j=1}^m s(\lambda_j)\phi_j(x)\phi_j(x')$$

uses eigenfunctions $\phi_j$ that satisfy the boundary conditions, with $s$ the spectral density of the base kernel. The boundary is thus enforced by the vanishing of $\phi_j$ or their derivatives on $\partial\Omega$ (Solin et al., 2019, Swiler et al., 2020, Gulian et al., 2020).

• Physics-informed tensorized kernels: For high-dimensional cubes with Robin boundaries, the full kernel is built as a tensor product of 1D constrained kernels:

$$k^{\rm (T)}(x, x') = \prod_{l=1}^d k_{\nu,\mathcal{B},c_l}(x_l, x_l')$$

for each $x \in [0,1]^d$ (Ding et al., 12 Jul 2025).

• Projection in RKHS: For boundary data on an uncountable set $T_0$, the conditional mean and covariance are

$$\mu_0(s) = \mu(s) + \langle k_s|_{T_0}, (g-\mu)|_{T_0}\rangle_{\mathcal{H}(T_0)},\qquad k_0(s,s') = k(s,s') - \langle k_s|_{T_0}, k_{s'}|_{T_0}\rangle_{\mathcal{H}(T_0)},$$

where $g$ is the imposed constraint and $\mathcal{H}(T_0)$ the RKHS restricted to $T_0$ (Brown et al., 2022).

4. Implementation and Approximation Strategies

The high computational cost of enforcing constraints on irregular domains or in high dimensions has motivated sophisticated approximation schemes:

• Finite Element Modeling (FEM): The true kernel is projected onto a finite basis (e.g., B-splines on a mesh). Coupled Monte Carlo (MC) estimators are employed to estimate kernel entries via simulations of stochastic processes (e.g., Brownian motion hitting times) so that the positive-definiteness and boundary properties are preserved (Ding et al., 12 Jul 2025). FEM-based kernel approximation yields error bounds that decay at rates competitive with problem smoothness and basis density.
• Low-Rank and Reduced-Rank Formulations: Direct expansion into $m \ll n$ features allows GP prediction at $O(n m^2)$ and hyperparameter learning at $O(m^3)$, much faster than $O(n^3)$ for full GPs (Riutort-Mayol et al., 2020, Solin et al., 2019).
• Tensor Construction in High Dimensions: For product domains, tensorizing the one-dimensional boundary-constrained kernel reduces the complexity of constructing global kernels while guaranteeing per-coordinate boundary satisfaction (Ding et al., 12 Jul 2025).
• Sampling from Truncated Gaussians: When bounding constraints reduce to convex restrictions on expansion coefficients (e.g., for bounded or monotonic GPs), efficient sampling and MAP estimation rely on numerical algorithms for high-dimensional truncated Gaussians (Maatouk, 2017, Agrell, 2019).

5. Theoretical Guarantees: Smoothness and Convergence

BCGPs provide rigorous control over smoothness and convergence properties. For the BdryMatérn GP, for instance, it is proven that:

• Smoothness Control: With smoothness parameter $\nu$, sample paths are almost surely $\lceil \nu \rceil-1$ times differentiable, matching the unconstrained Matérn case (Ding et al., 12 Jul 2025). This is established using spectral decomposition and elliptic PDE (Sobolev embedding) results.
• Approximation Error Analysis: The coupled MC–FEM approach comes with explicit $L^2$-error bounds, which guarantee that the numerical approximation of the kernel—and thus the posterior—can be made arbitrarily close to the (intractable) true boundary-constrained process as the basis is refined (Ding et al., 12 Jul 2025).
• Improved Convergence Rates: When boundary information is included—either as known values or structural constraints—the resulting GP approximation can achieve deterministic $L^p$ error rates as fast as $o(n^{-1})$ and uniform errors $O(n^{-1}[{\log n}]^{2d-1})$, which are substantially improved over unconstrained GPs, particularly in high dimensions (Ding et al., 2019). Resistance to the curse of dimensionality is an important consequence.

6. Representative Applications

BCGPs are deployed in a wide spectrum of scientific and engineering domains where boundary constraints are physically or empirically mandated:

• Surrogate Modeling for Computer Experiments: Many computer codes, especially those simulating PDE-governed systems (e.g., fluid flow, elasticity), return outputs with known behavior on the boundary. BCGPs enable robust emulation with limited internal data by leveraging this knowledge (Ding et al., 2019, Ding et al., 12 Jul 2025).
• Physics-Informed Inverse Problems: Reconstruction of incompressible 2D velocity fields around aerodynamic profiles is achieved by constructing GPs over stream functions, ensuring both divergence-free velocity and slip conditions at solid boundaries (Padilla-Segarra et al., 23 Jul 2025).
• Uncertainty Quantification in Structural Health Monitoring: GPs constrained with boundary-informed covariance kernels provide accurate and stable prediction of acoustic emission time-of-arrival maps on complex plate geometries, improving error and variance without requiring boundary observations (Jones et al., 2022).
• Constrained Density Estimation and Regression: Projection and truncation schemes ensure credible intervals and posterior samples respect nonnegativity, boundedness, or monotonicity, with closed-form posterior mean and variance expressions for efficient computation (Zhang et al., 2018, Agrell, 2019).
• Probabilistic Solvers for PDEs: B-EPGPs use spectral or algebraic basis constructions so that every sample is a solution to the specified PDE and boundary conditions, providing accuracy and physical correctness unattainable via neural operators or unconstrained surrogates, especially on complex domains or with variable coefficients (Huang et al., 25 Nov 2024).

7. Numerical Case Studies and Empirical Results

Empirical studies highlight key outcomes:

• On irregular 2D domains (T-shaped, ring, disk, or holed rectangles), the BdryMatérn GP produces significantly lower log-MSE and visually superior interpolants at boundaries versus traditional Matérn GPs (Ding et al., 12 Jul 2025).
• In $d=30$ dimensions, tensorized BdryMatérn GPs outperform high-dimensional product GPs on structured test functions, with gains in predictive error and variance quantification.
• For PDE-constrained problems, B-EPGP shows 1–2 orders of magnitude improvement in both absolute and relative errors compared to leading neural operator models, while retaining exact satisfaction of energy conservation and boundary invariants (Huang et al., 25 Nov 2024).
• Reconstructions of flow fields around airfoils achieve error reductions up to $10^{-4}$ for the slip condition accuracy, effective even without boundary observations (Padilla-Segarra et al., 23 Jul 2025).
• Monte Carlo/FEM estimator preserves the positive-definiteness of the approximate kernel almost surely, and error bounds scale nearly optimally with the number of basis points (Ding et al., 12 Jul 2025).

Conclusion

Boundary-Constrained Gaussian Processes provide a comprehensive probabilistic modeling toolkit for domains where physical, geometric, or shape constraints must be incorporated at the model level. Recent advances span from spectral and SPDE-based kernel engineering to computer algebraic frameworks and efficient FEM-based numerical implementations. BCGPs rigorously satisfy boundary conditions, offer strong smoothness and convergence guarantees, and dramatically enhance predictive performance and credible uncertainty quantification over unconstrained or naively regularized GPs, especially on irregular domains, in high dimensions, or where data is scarce away from boundaries. These models have set a new standard for surrogate modeling in scientific and engineering contexts where boundary and shape constraints are intrinsic to the problem.