Boundary-Respecting GP Priors

• Boundary-respecting GP priors are Gaussian processes whose sample paths are constructed to exactly satisfy specified boundary conditions such as Dirichlet or Neumann.
• They enhance Bayesian surrogate modeling for PDEs by incorporating physical constraints, resulting in improved convergence rates and robust uncertainty quantification.
• Construction techniques include harmonic eigenfunctions, constrained random fields, and algebraic frameworks, enabling efficient applications in spatial statistics and physics-informed emulation.

Boundary-respecting Gaussian process (GP) priors are a class of probabilistic models constructed so that all sample paths satisfy specified boundary conditions—such as Dirichlet, Neumann, or more general linear boundary operators—on the domain boundary. These priors are central to enforcing physical or structural properties in Bayesian surrogate modeling for differential equations, scientific emulation, statistical learning on spatial domains, and physics-informed machine learning. Their defining property is exact satisfaction of boundary constraints, either by analytic construction of the kernel or by algebraic or operator-theoretic transformations of standard GP priors. Development in this area encompasses kernel engineering, spectral approaches, operator-theoretic algebra, and connections to classical numerical analysis.

1. Mathematical Formulation and Key Principles

A boundary-respecting GP prior is defined as a Gaussian process $f: \overline{\Omega} \to \mathbb{R}$ (or $\mathbb{R}^p$), for a domain $\Omega \subset \mathbb{R}^d$ with boundary $\partial \Omega$, such that sample paths satisfy a system of PDEs and/or boundary conditions: $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$ where $L(\partial)$ is a (typically linear, constant-coefficient) differential operator, and $B(\partial)$ is a (possibly vector-valued) boundary operator, such as identity (Dirichlet), normal derivative (Neumann), or more general (Robin, mixed).

Construction strategies differ among the literature, but the principal goal is that the space of realizations of the GP prior lies within the solution set

$$\mathcal{H} = \{u \in C^\infty(\overline{\Omega}): L(\partial)u = 0,\, B(\partial)u = 0 \}.$$

Exact enforcement of boundary conditions is achieved through basis modification, kernel engineering, feature-space constrained representations, or operator conditioning (Huang et al., 2024, Solin et al., 2019, Ma et al., 28 Nov 2025, Lange-Hegermann, 2020, Ding et al., 2019).

2. Construction Techniques and Frameworks

Several rigorous frameworks for constructing boundary-respecting GP priors are established in the literature:

A. Boundary Ehrenpreis-Palamodov GP (B-EPGP) (Huang et al., 2024)

• Relies on the Ehrenpreis-Palamodov fundamental principle: for linear constant-coefficient PDEs, exponential-polynomial functions $e^{x\cdot z}$, $z$ in the characteristic variety $\mathbb{R}^p$0, span the solution space.
• The prior takes the form $\mathbb{R}^p$1 with basis elements $\mathbb{R}^p$2 crafted to satisfy both $\mathbb{R}^p$3 in the domain and $\mathbb{R}^p$4 on the boundary.
• The kernel is $\mathbb{R}^p$5 for a diagonal weight matrix $\mathbb{R}^p$6, or in the infinite-rank limit by spectral integration over $\mathbb{R}^p$7.
• Boundary satisfaction is algebraic: for each frequency, coefficients $\mathbb{R}^p$8 are found such that $\mathbb{R}^p$9, where $\Omega \subset \mathbb{R}^d$0 is a finite fiber over the boundary variety.

B. Harmonic Feature (Spectral) GPs (Solin et al., 2019)

• Constructs GPs using Laplace eigenfunctions $\Omega \subset \mathbb{R}^d$1 of $\Omega \subset \mathbb{R}^d$2 on $\Omega \subset \mathbb{R}^d$3 under specified boundary conditions.
• The kernel is approximated as $\Omega \subset \mathbb{R}^d$4, where $\Omega \subset \mathbb{R}^d$5 are eigenvalues of $\Omega \subset \mathbb{R}^d$6 and $\Omega \subset \mathbb{R}^d$7 is the spectral density of the target stationary kernel.
• This ensures exact vanishing (Dirichlet) or zero normal derivative (Neumann) on $\Omega \subset \mathbb{R}^d$8.

C. Constrained Gaussian Random Fields (cGRFs) (Ma et al., 28 Nov 2025)

• Enforces continuous linear boundary restrictions $\Omega \subset \mathbb{R}^d$9 for $\partial \Omega$0 by constructing a corrected field

$\partial \Omega$1

where $\partial \Omega$2 is the unconstrained field, $\partial \Omega$3 a boundary projection, $\partial \Omega$4 vanishes off $\partial \Omega$5.

• The resulting process has closed-form mean and covariance and satisfies all boundary operators exactly.

D. Algebraic/Gröbner Frameworks (Lange-Hegermann, 2020)

• For multi-output systems, solution spaces are parametrized via chained nullspace constructions in appropriate operator algebras—first, the nullspace of the PDE (using Weyl algebra and Gröbner bases), and second, nullspace of the boundary module. Their intersection via pullback gives a generator set for the solution space.
• GPs are built on this parametrization: $\partial \Omega$6, where $\partial \Omega$7 is the basis matrix.

E. BdryGP Model (Ding et al., 2019)

• For Dirichlet boundaries on a hypercube, constructs a mean interpolant matching boundary values and a tensor-product kernel $\partial \Omega$8 that vanishes whenever any argument is on a known boundary.
• The kernel is built from 1D building blocks, explicit in terms of $\partial \Omega$9 and $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$0 functions, ensuring exact satisfaction of Dirichlet boundaries.

3. Exactness, Theoretical Guarantees, and Convergence

Boundary-respecting priors possess several rigorously established properties:

• Exactness: All sample paths are in the intersection of the null spaces of the governing operator and the boundary operator: $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$1. No realization violates the enforced boundary conditions, in contrast to standard GPs with unconstrained kernels (Huang et al., 2024, Ma et al., 28 Nov 2025, Lange-Hegermann, 2020).
• Posterior Existence and Uniqueness: For GPs with finite-dimensional (e.g., harmonic or exponential) bases, posteriors given pointwise Gaussian observations remain Gaussian with closed-form mean and covariance (Huang et al., 2024).
• Density: The constructed basis elements (e.g., boundary-respecting exponentials, separated Laplacian eigenfunctions) are dense in the space of solutions with given boundary conditions. By reflection and symmetrization, classical Fourier-type bases are recovered (Huang et al., 2024).
• Improved Convergence: Incorporation of boundary information leads to accelerated convergence rates in function approximation and regression tasks. For the BdryGP model on $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$2 with full Dirichlet information and Smolyak sparse grids, $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$3-error decays as $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$4 and sup-norm as $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$5, an exponential improvement over standard rates (Ding et al., 2019).
• Uncertainty Quantification: Uncertainty collapses to zero at constrained boundaries, yielding physically interpretable credible intervals and risk assessments (Ma et al., 28 Nov 2025).

4. Implementation and Computational Considerations

Implementation varies with basis and kernel type:

• Basis Computation: Boundary-respecting exponentials and Laplace/Helmholtz eigenfunctions require symbolic construction or numerical eigenvalue decompositions (e.g., ARPACK). Gröbner basis methods are data-independent but can be computationally intensive in high dimensions or with many operators (Huang et al., 2024, Lange-Hegermann, 2020, Solin et al., 2019).
• Kernel Assembly: Finite-rank harmonic or exponential kernels can be evaluated with $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$6 cost ($$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$7 data points, $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$8 basis elements). Infinite-dimensional (integral) forms are typically truncated for computation.
• Boundary Projections: Construction of cGRF priors requires continuous projection maps from the domain to the boundary; guaranteed for convex domains but complicated for nonconvex geometries (Ma et al., 28 Nov 2025).
• Feature Space Efficiency: Low-rank approximations using harmonic features or sparse grid designs enable scalable inference, accommodating both Gaussian and non-Gaussian likelihoods via variational optimization (Solin et al., 2019, Ding et al., 2019).
• Regularity and Kernel Choice: For derivative boundary conditions, base kernels must admit derivatives of order $$L(\partial)f(x) = 0, \ x \in \Omega; \quad B(\partial)f(x) = 0, \ x \in \partial \Omega,$$9 (for $L(\partial)$0-th order derivatives to be enforced). This restricts the choice to sufficiently smooth RKHSs.

5. Practical Applications and Empirical Performance

Boundary-respecting GP priors are deployed in:

• Physics-Informed Surrogates: Solution of linear PDEs with constant coefficients and linear boundary conditions, often for forward or inverse modeling in computational science (Huang et al., 2024, Ma et al., 28 Nov 2025).
• Scientific Emulation: High-fidelity emulators for codes that simulate physical systems, with boundary data known from governing physics (Ding et al., 2019).
• Probabilistic Numerics: Bayesian solvers for PDEs—quantifying discretization error, ensuring uncertainty contracts near boundaries (Ma et al., 28 Nov 2025).
• Learning Dynamical Systems: Physics-informed data-driven discovery (e.g., of nonlinear evolution equations) with reduced false discoveries and improved parameter recovery by strict enforcement of boundary information (Ma et al., 28 Nov 2025).
• Spatial Statistics: Regression and classification on spatial domains with complex or irregular boundaries, employing Laplace eigenfunctions tailored to domain topology (Solin et al., 2019).

Empirically, these priors deliver uniformly lower errors compared to unconstrained alternatives and exhibit robustness to increasing dimensionality. For example, B-EPGP achieves $L(\partial)$1-absolute errors $L(\partial)$2 in canonical PDE setups, outperforming neural operators by orders of magnitude (Huang et al., 2024). BdryGP outpaces standard GP predictors by a factor of two in slope on log–error vs. log–$L(\partial)$3 plots and damps the curse of dimensionality to a polylogarithmic penalty when boundary information is complete (Ding et al., 2019).

6. Limitations and Open Challenges

• Geometry Restrictions: Many frameworks depend on convexity or specific boundary geometry (flat, polygonal) to guarantee continuous projections or terminate reflection constructions. Nonconvex domains or curved boundaries necessitate hybrid or numerical treatment (Ma et al., 28 Nov 2025, Huang et al., 2024).
• Computational Cost: Computing eigenbases, symbolic nullspaces, or large numbers of derivative kernel terms can be burdensome in high dimension ($L(\partial)$4). Sparse grid, low-rank, or multigrid approximations partially mitigate this cost (Ding et al., 2019, Solin et al., 2019).
• Regularity Requirements: Contraints on derivative orders mandate kernels with correspondingly smooth sample paths (e.g., Matérn kernels with high $L(\partial)$5). For rough processes or limited smoothness, boundary enforcement can fail.
• Domain Generality: Arbitrarily complex or disconnected domains challenge current methods; general solutions exist only for classes where projections and weight functions can be globally defined and regular (Solin et al., 2019, Ma et al., 28 Nov 2025).
• Non-Uniqueness: Multiple representations of the constrained prior (choice of projections, weights, or basis orderings) yield the same law, but implementation efficiency and numerical stability can vary substantially.

7. Theoretical and Computational Connections

Boundary-respecting GP priors reveal deep connections between probabilistic inference and classical numerical analysis:

• Finite Element Links: The BdryGP model with Brownian-type kernels, for sparse or full grids, reproduces piecewise-linear finite element interpolants; hence, GP regression with such priors can be interpreted as probabilistic analogues of FEM (Ding et al., 2019).
• Spectral Duality: Spectral constructions using domain-respecting eigenfunctions (harmonic, Laplacian) generalize random Fourier features to arbitrary domains and information-enforcing settings (Solin et al., 2019).
• Operator Algebra: Gröbner-based parametrization of solution spaces in boundary-respecting frameworks provides an algebraic underpinning for symbolic enforcement (Lange-Hegermann, 2020).

These developments position boundary-respecting GP priors at the interface of functional analysis, Bayesian inference, computational mathematics, and scientific machine learning. The field is characterized by ongoing refinement in extending to nonlinear operators, more general domains, lower regularity, and efficient numerical realization.