Physics-Constrained Gaussian Process Regression

Updated 11 April 2026

Physics-constrained Gaussian Process Regression is a method that integrates physical laws, such as PDEs and boundary conditions, into GP models for precise and data-efficient inference.
It employs eigenfunction-based kernel construction and spectral expansion to ensure that all samples rigorously satisfy the imposed physical and boundary constraints.
The approach supports multifidelity and hybrid deep-kernel designs, offering scalable solutions with robust uncertainty quantification in high-dimensional or scarce data scenarios.

Physics-constrained Gaussian Process Regression (GPR) refers to a class of Gaussian process modeling techniques in which physical laws—typically in the form of partial differential equations (PDEs), boundary conditions, or more generally linear operator constraints—are exactly or probabilistically enforced on the GP prior, the kernel, the mean, or the structure of inference. These methods allow the solution or surrogate to be informed not only by data but also by mechanistic knowledge, thereby enabling data-efficient, stable, and physically consistent regression for high-dimensional or scarce-data problems commonly encountered in scientific and engineering domains.

1. Foundations: GP Priors under Differential and Boundary Constraints

Physics-constrained GPR formalizes the inference of latent functions $u(x)$ that are known, a priori, to satisfy a (typically linear) boundary-value problem

$L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$

where $L$ is a linear second-order differential operator (e.g., elliptic, parabolic, hyperbolic), and $B$ denotes homogeneous boundary conditions such as Dirichlet ( $u=0$ ), Neumann ( $\partial u/\partial n=0$ ), or mixed types (Gulian et al., 2020). The essence of the physics-constrained paradigm is to choose a GP prior on $u$ such that every prior sample $u\sim GP(m_u, K_u)$ satisfies—either exactly or to a quantifiable error—the required operator and boundary constraints (Pförtner et al., 2022, Henderson et al., 2021).

For linear operator constraints, necessary and sufficient conditions for a stochastic process to satisfy $T[u]=0$ in the sense of distributions is that both the mean $m_u$ and each kernel section $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 0 are annihilated by $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 1; i.e., $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 2 and $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 3 in $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 4 (Henderson et al., 2021).

2. Physics-Constrained Kernel Construction

A canonical approach to building physics-constrained kernels is to expand $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 5 in a truncated set of eigenfunctions $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 6 of the homogeneous BVP,

$L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 7

and define

$L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 8

with $L u(x) = f(x), \quad x\in\Omega,\qquad B u(x) = 0, \quad x\in\partial\Omega$ 9 determined from the spectral density of a parent kernel or learned from data. By restricting to $L$ 0, all samples and the posterior mean exactly satisfy $L$ 1 (Gulian et al., 2020). The kernel for $L$ 2 is then $L$ 3. Cross-covariances $L$ 4 are similarly constructed.

Alternatively, for vector-valued PDEs (e.g., incompressible fluid flow), divergence-free and boundary-admissible covariance structures are constructed by taking vector differential operators of a scalar kernel and applying continuous boundary projections via the Karhunen–Loève (Mercer) expansion restricted to the domain boundary (Padilla-Segarra et al., 23 Jul 2025). Recent frameworks also enable arbitrary bounded linear operator constraints, allowing the use of observation functionals corresponding to point evaluation, domain integrals, or weighted residuals (Pförtner et al., 2022).

3. Posterior Inference and Reduced-Rank Computation

Given observations $L$ 5 potentially of both $L$ 6 and $L$ 7 at locations $L$ 8, the joint GP for $L$ 9 has a block covariance assembled by the above operator actions. The posterior mean and covariance for any new $B$ 0 are then computed via standard GP conditioning: $B$ 1

$B$ 2

with $B$ 3 denoting the augmented covariance with noise terms (Gulian et al., 2020).

For spectral expansions, computation can be reduced to $B$ 4 (assembly), $B$ 5 (inversion), and $B$ 6 (per-prediction) via reduced-rank matrix identities (Woodbury formula). This provides significant acceleration when $B$ 7 (number of modes) is much less than $B$ 8 (number of observations) (Gulian et al., 2020).

4. Comparison to Unconstrained or Softly Constrained Approaches

Standard "physics-informed" GPR constructs the joint kernel by acting differential operators on stationary kernels (e.g., squared-exponential) and typically enforces PDE or boundary constraints via penalization, "soft" virtual points, or by adding data at the boundary (Gulian et al., 2020, Cross et al., 2023). Such approaches cannot ensure zero posterior variance at $B$ 9, nor exact satisfaction of boundary conditions—posterior means may drift with sparse data, and credible intervals inflate near the boundaries.

By leveraging eigenfunction-based kernels or direct projection in function space, physics-constrained GPR ensures that both the posterior mean and samples lie precisely within the constraint manifold. This yields sharper extrapolation, precise uncertainty quantification, and increased stability, especially in the sparse-data regime (Gulian et al., 2020, Pförtner et al., 2022).

5. Integration with Multifidelity and Data-Driven Enhancements

Modern frameworks extend physics-constrained GPR with multifidelity and discrepancy models. A low-fidelity GP is constructed via physics-informed Kriging from stochastic realizations of imperfect simulations, enforcing the operator constraint up to a quantifiable error. A data-driven GP for the discrepancy is then co-kriged with the physics GP, with the joint model trained via marginal likelihood. Error bounds for constraint preservation, as well as active learning strategies targeting uncertainty maxima, are naturally accommodated (Yang et al., 2018, Yang et al., 2018).

Hybrid schemes also combine deep-kernel learning with physics regularization, using neural-network-based feature maps where the loss function includes both the standard GP likelihood and a physics-induced penalty (e.g., Boltzmann–Gibbs weight for the PDE residual). This approach enables scalable inference in high-dimensional or complex-geometry settings where traditional spectral methods may be prohibitive (Chang et al., 2022).

6. Applications and Extensions

Physics-constrained GPR has been deployed for:

Inverse problems and state estimation in PDEs, such as inferring source terms, initial/boundary conditions, or unobserved states (Gulian et al., 2020, Ma et al., 2020, Henderson et al., 2023);
Surrogate modeling in engineering, combining physical structure with data-driven flexibility, e.g., in structural mechanics and fluid flow (Cross et al., 2023, Padilla-Segarra et al., 23 Jul 2025);
Multifidelity model hierarchies, where heterogeneous codes or simulation levels are integrated (Yang et al., 2018, Yang et al., 2018);
Quantification and propagation of uncertainty due to input coefficients or boundary data in PDE-constrained modeling (Pförtner et al., 2022, Ma et al., 2020);
Problems requiring exactly divergence-free or boundary-constrained vector fields, as in incompressible aerodynamics (Padilla-Segarra et al., 23 Jul 2025);
Specialized enforcement of physics-inspired inequality, monotonicity, or derivative constraints in material models or constitutive laws (Li et al., 6 Jan 2026, Tran et al., 2022, Kochan et al., 2024).

7. Computational Considerations and Limitations

Exact enforcement of physics and boundary constraints comes at the price of computing eigenfunctions on arbitrary domains, which is tractable in regular geometries or via numerical eigensolvers. For large-scale or high-dimensional domains, scalable variants using inducing points, deep-kernel surrogates, or stochastic approximations are under development (Chang et al., 2022). The choice of basis size $u=0$ 0 sets a tradeoff between expressivity and computational cost. In problems where the PDE coefficients themselves are uncertain, the joint GP is extended to include these as random fields, enabling full Bayesian propagation of parametric uncertainty to the solution (Pförtner et al., 2022).

Physics-constrained GPR is not immune to model-form misspecification; if the physical operator or boundary description omits crucial mechanisms, the data-misfit will localize near these regions, and the model may be over-constrained or under-adaptive. In hybrid designs, identifiability between the physics-informed and data-driven kernel components may require careful regularization or cross-validation (Cross et al., 2023).

Physics-constrained Gaussian Process Regression constitutes a unifying probabilistic framework that rigorously interpolates between classic numerical PDE discretization, statistical machine learning, and inverse modeling. By exactly encoding physics into the model, it substantially reduces data requirements, provides robust extrapolation guarantees, and equips practitioners with a principled mechanism for uncertainty quantification and error control across a wide range of scientific and engineering inference tasks (Gulian et al., 2020, Pförtner et al., 2022, Padilla-Segarra et al., 23 Jul 2025, Li et al., 6 Jan 2026).