Physics-Informed Gaussian Processes

• Physics-informed Gaussian processes are Bayesian models that embed physical laws into their priors, kernels, or mean functions to ensure consistency with known dynamics.
• They integrate operator constraints, simulator-informed priors, and discrepancy modeling to bridge analytical and data-driven insights across PDEs, ODEs, and control applications.
• These approaches enable efficient uncertainty quantification and scalable inference in complex systems, facilitating robust predictions in real-world scientific and engineering problems.

Searching arXiv for recent and foundational papers on physics-informed Gaussian processes. Physics-informed Gaussian processes (GPs) are Gaussian-process models in which physical knowledge is embedded into the prior, covariance construction, mean function, observation model, or posterior constraints so that inference over latent fields, dynamical trajectories, or decision boundaries reflects known mechanistic structure while retaining Bayesian uncertainty quantification. In the literature, the term spans several distinct but related constructions: operator-constrained GP regression for linear PDEs, priors derived from stochastic simulators or analytical models, ODE-constrained trajectory priors for control, multifidelity models that combine physics-based ensembles with data-driven discrepancy terms, and hybrid pipelines in which GPs denoise or regularize data used by downstream physics-informed learners (Pförtner et al., 2022, Long et al., 2022, Yang et al., 2018).

1. Scope and major formulations

The phrase “physics-informed GP” does not denote a single canonical model. In one prominent formulation, solving a linear PDE is recast as GP regression on an unknown field $u$, with the PDE operator, boundary operators, residual functionals, and measurements all treated as linear observations of a GP prior (Pförtner et al., 2022). In another, physical knowledge appears as an informative prior mean or kernel rather than as exact operator constraints: a fast accelerator simulator determines both the GP mean and a Hessian-derived multivariate RBF kernel for online optimization (Hanuka et al., 2020), aerodynamic regression models provide a mean function for estimating pitching-moment coefficients from arbitrary flight-test maneuvers (Harp et al., 2 Jan 2025), and CALPHAD or analytical materials models provide prior logits or regression means for constraint-aware alloy design (Hardcastle et al., 17 Feb 2025).

A separate branch uses stochastic physics solvers to estimate the GP prior nonparametrically. Physics-informed Kriging and Physics-informed CoKriging compute the prior mean and covariance from ensembles generated by stochastic physical models and then add a discrepancy GP to reconcile low-fidelity physics with high-fidelity observations (Yang et al., 2018). In control, the “physics-informed” component can be even stronger: all GP sample paths are constrained to satisfy a linear time-invariant ODE system by construction, so that control is performed by conditioning a trajectory prior rather than by optimizing an unconstrained surrogate (Tebbe et al., 2024).

A concise taxonomy appears below.

Formulation	How physics enters	Representative paper
Operator-constrained field inference	Linear differential and boundary operators act on a GP prior	(Pförtner et al., 2022)
Physics-informed priors from simulators	Simulator defines prior mean and/or kernel	(Hanuka et al., 2020)
Physics-informed multifidelity GP	Physics ensemble gives low-fidelity GP; discrepancy learned from data	(Yang et al., 2018)
ODE-constrained trajectory prior	GP sample paths lie in the nullspace of a differential operator	(Tebbe et al., 2024)
Residual-likelihood variational GP	Differential equations enter as Gaussian likelihood terms at collocation points	(Long et al., 2022)
Hybrid GP–PINN preprocessing	GP smooths noisy boundary or initial data before PINN training	(Bajaj et al., 2021)

This heterogeneity is significant because claims about “hard” enforcement, uncertainty calibration, or computational scaling depend strongly on which formulation is under discussion. A common misconception is that every physics-informed GP encodes the governing equations directly in the prior. The GP-PINN hybrid of noisy-data denoising is explicitly not of that type: there the GP does not encode the PDE and instead acts as a smoother and uncertainty tagger for boundary or initial data (Bajaj et al., 2021).

2. Operator-theoretic foundations

A central mathematical foundation is closure of Gaussian processes under bounded linear operators. For a prior

$u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$

and noisy operator observations

$$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$

the induced quantity $O u$ is again Gaussian, with joint covariances

$$\operatorname{Cov}\big(u,\; O u\big) = k\, O^{\!*}, \qquad \operatorname{Cov}\big(O u,\; O u\big) = O k\, O^{\!*},$$

and posterior

$$m_{\text{post}} = m + k\, O^{\!*}\big(O k\, O^{\!*} + \Sigma\big)^{-1}\big(y - O m\big),$$

$$k_{\text{post}} = k - k\, O^{\!*}\big(O k\, O^{\!*} + \Sigma\big)^{-1} O k.$$

This generalization supports point evaluations, integrals, weak residuals, collocation residuals, and boundary operators within one inference framework (Pförtner et al., 2022).

For linear PDEs on a bounded domain $\Omega \subset \mathbb{R}^d$,

$$L u = f \text{ in } \Omega, \qquad B u = c \text{ on } \partial\Omega,$$

physics enters by defining operator observations from weighted residuals and boundary conditions. With residual $R(u)=Lu-f$, one may enforce

$u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$0

where $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$1 may be weak-form integrals or strong-form point evaluations (Pförtner et al., 2022). This places collocation, finite volume, pseudospectral, Galerkin, finite element, and spectral methods inside a single probabilistic construction. Under recovery-prior conditions aligned with a given trial space, the posterior mean reproduces the classical weighted-residual solution exactly, while the posterior covariance quantifies uncertainty orthogonal to the chosen approximation space (Pförtner et al., 2022).

The same operator calculus underlies multi-output structural models. For Kirchhoff–Love plates, a GP prior on deflection $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$2 is pushed through linear differential operators to obtain a joint GP over deflection, rotations, curvatures, loads, shear forces, and bending moments, with flexural rigidity $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$3 entering the covariance analytically (Kavrakov et al., 2024). For Timoshenko beams, a latent GP on bending deflection induces cross-covariances among transverse displacement, rotation, strain, bending moment, shear force, and load, enabling Bayesian identification of $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$4 and $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$5 from heterogeneous sensors (Tondo et al., 2023).

A structurally analogous construction appears in port-Hamiltonian learning. There, a scalar GP prior is placed on the Hamiltonian $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$6, and the dynamics are induced through the port-Hamiltonian map

$u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$7

so that the learned vector field remains passive because $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$8 is skew-symmetric and $u \sim \mathrm{GP}\big(m(x),\, k(x,x')\big),$9 (Beckers et al., 2023). For robot inverse dynamics, the latent quantities are kinetic and potential energy rather than torques directly; Lagrange’s equations then define a multi-output GP over joint torques via linear differential operators on those latent energy GPs (Giacomuzzos et al., 2023).

3. Mechanisms for incorporating physics

Physics can enter a GP through several non-equivalent mechanisms. The strongest form is exact satisfaction by construction. In linear model predictive control, a GP over the stacked trajectory $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$0 is built so that every sample path satisfies the ODE $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$1. This is achieved by representing trajectories in the nullspace of a polynomial matrix differential operator using a Smith Normal Form construction; conditioning on setpoints then yields dynamically consistent control trajectories (Tebbe et al., 2024). A related MPC formulation adds truncation and Hamiltonian Monte Carlo to enforce box constraints with formal open-loop constraint-satisfaction guarantees while maintaining the ODE structure on the discretized grid (Tebbe et al., 20 Nov 2025).

A second mechanism is soft enforcement through likelihood terms on residuals. AutoIP constructs a joint GP prior over function values, derivatives, and latent source terms using kernel differentiation, then imposes the differential equation through a Gaussian residual likelihood at collocation points. Because nonlinear residuals are generally non-Gaussian under a GP prior, inference is performed with whitening and stochastic variational optimization rather than closed-form conditioning (Long et al., 2022). Physics-Informed Variational State-Space GPs place the temporal kernel in state-space form, use collocation and boundary likelihoods, and employ natural-gradient variational inference with Gauss–Newton curvature to obtain linear-in-time computation for both linear and nonlinear spatio-temporal constraints (Hamelijnck et al., 2024).

A third mechanism is prior-mean injection. In flight-test analysis, the Morelli aerodynamic regression model is used as the GP mean function for pitching moment coefficient $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$2, so the GP learns only the discrepancy between an analytical baseline and T-38 data (Harp et al., 2 Jan 2025). In alloy design, prior class probabilities derived from CALPHAD, VEC heuristics, or analytical strength models are mapped into latent logits or regression means, anchoring decision boundaries before data correction (Hardcastle et al., 17 Feb 2025). This strategy can be computationally lightweight because it leaves the kernel mostly conventional.

A fourth mechanism is prior-kernel design from physics. In online accelerator optimization, the Hessian of the log simulated objective near the simulator optimum defines the precision matrix of a multivariate RBF kernel, including off-diagonal correlations among the 13 skew quadrupole controls; the prior mean is set to the far-field simulator offset (Hanuka et al., 2020). In GP-PHS, the Hamiltonian kernel is differentiated and mapped through the port-Hamiltonian structure to produce a matrix-valued kernel for the dynamics (Beckers et al., 2023). Inverse-dynamics identification likewise uses a “Lagrangian Inspired Polynomial” kernel reflecting polynomial structure in kinetic and potential energies (Giacomuzzos et al., 2023).

A fifth mechanism is ensemble-based prior estimation. Physics-informed CoKriging estimates the low-fidelity GP mean and covariance directly from ensembles of a stochastic physics-based model, then fits a discrepancy GP and a scaling coefficient $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$3 against high-fidelity data (Yang et al., 2018). This avoids full hyperparameter optimization for the low-fidelity process and carries linear physical constraints into the posterior up to an explicit error bound (Yang et al., 2018).

Finally, topology-optimization work uses GP priors with a shared neural-network mean function and independent kernels for state and design variables. Boundary conditions are enforced via conditioning, while objective, PDE residuals, and design constraints are enforced through a composite training loss. In the Stokes/Brinkman setting, this creates a simultaneous and meshfree optimization procedure in which state and design are updated in one loop (Yousefpour et al., 2024, Yousefpour et al., 18 Mar 2025, Sun et al., 14 Jul 2025).

4. Inference, uncertainty quantification, and computation

One of the primary motivations for physics-informed GPs is uncertainty quantification. In the operator-regression view of linear PDEs, the posterior covariance quantifies discretization error, uncertainty propagated from uncertain coefficients in $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$4, and the effect of noisy physical measurements combined with operator observations (Pförtner et al., 2022). For any linear functional $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$5, the posterior variance is given by applying $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$6 to the posterior covariance operator, so uncertainty can be assessed not only pointwise but also on integrals, boundary fluxes, or other derived quantities (Pförtner et al., 2022).

In structural mechanics, uncertainty is often propagated through explicit Bayesian parameter inference. For Kirchhoff–Love plates, flexural rigidity $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$7, kernel hyperparameters, and sensor-noise levels are sampled with MCMC after constructing a multi-output GP from the governing operators, which enables stochastic estimation of deflection, moments, and loads from heterogeneous measurements (Kavrakov et al., 2024). For Timoshenko beams, Metropolis–Hastings is used to infer $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$8, $$y = O u + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \Sigma),$$9, and noise terms, after which the posterior predictive distribution over displacements, rotations, strains, and internal forces is obtained by averaging Gaussian conditionals over sampled parameters (Tondo et al., 2023).

The computational bottleneck of exact GP inference remains central. Simulator-informed accelerator optimization uses exact GP updates with $O u$0 scaling in the number of online observations and does not employ sparse approximations (Hanuka et al., 2020). The constrained-MPC formulation similarly requires inversion of dense covariance matrices on the horizon and then HMC sampling from a truncated multivariate normal; this restores guarantees but adds a second computational layer beyond closed-form Gaussian updates (Tebbe et al., 20 Nov 2025). By contrast, AutoIP introduces inducing-point sparsification and whitening so that variational inference can be performed with mini-batching and reduced complexity (Long et al., 2022), and Physics-Informed Variational State-Space GPs obtain linear-in-time scaling by exploiting state-space representations of temporal kernels and Kalman filtering/smoothing recursions (Hamelijnck et al., 2024).

Uncertainty can also guide data acquisition and sensing. CoPhIK uses greedy variance maximization to select new observation locations under the multifidelity posterior (Yang et al., 2018). The alloy-design framework uses maximum Shannon entropy on class probabilities for active learning of phase boundaries or property-threshold decision surfaces (Hardcastle et al., 17 Feb 2025). The Timoshenko-beam model uses an entropy-based sensor placement criterion that accounts for heterogeneous sensor modalities and boundary conditions built into the physics-informed covariance (Tondo et al., 2023).

A recurring practical distinction is between exact and approximate constraint satisfaction. Conditioning on boundary observations in GP-based topology optimization enforces those sampled boundary conditions essentially exactly, because the posterior mean reproduces the conditioned values (Yousefpour et al., 2024). However, when physics is imposed via residual likelihoods or penalty terms, exact satisfaction is replaced by a tunable trade-off governed by likelihood variance, penalty coefficients, or variational approximation quality (Long et al., 2022, Hamelijnck et al., 2024).

5. Representative applications

Physics-informed GPs have been deployed across a wide range of scientific and engineering settings. In numerical analysis, they reformulate linear PDE solvers for Poisson, heat, and wave equations as operator-based GP regression, recovering classical discretizations in the posterior mean while adding structured error estimates (Pförtner et al., 2022). In structural mechanics, they have been used for stochastic inference in Kirchhoff–Love plates and Timoshenko beams from deflection, curvature, rotation, strain, and load measurements, with direct relevance to structural health monitoring (Kavrakov et al., 2024, Tondo et al., 2023).

In control and dynamical systems, ODE-constrained GPs have been used as linear model predictive controllers for tracking problems in LTI systems, with soft constraints introduced through virtual setpoints (Tebbe et al., 2024). A subsequent constrained formulation augments this with sampling from truncated Gaussians to obtain safe open-loop rollouts and formal support guarantees inside state and input bounds (Tebbe et al., 20 Nov 2025). Port-Hamiltonian learning uses Hamiltonian GPs to identify passive dynamical systems while preserving compositional interconnection properties (Beckers et al., 2023).

In robotics, GP priors on kinetic and potential energies, combined with Lagrange’s equations, have been used for black-box inverse-dynamics identification on a 7-DoF Franka Emika Panda and a 6-DoF MELFA RV4FL, where the resulting model also estimates latent energies without requiring energy labels (Giacomuzzos et al., 2023). In materials and design, physics-informed latent GP classifiers with CALPHAD- or VEC-derived priors improve low-data phase-stability prediction and active learning for alloy design constraints (Hardcastle et al., 17 Feb 2025).

In large-scale scientific facilities, simulator-informed GP optimization has been applied to 13-dimensional online tuning of SPEAR3 skew quadrupoles, with the kernel precision matrix obtained from the simulator Hessian and faster convergence than routinely used online optimizers (Hanuka et al., 2020). In remote sensing, joint GPs, latent force models, and automatic GP emulators have been used to fuse in situ and simulated data, model multi-output vegetation time series, and build compact look-up tables for costly radiative-transfer models (Camps-Valls et al., 2020).

Several recent works extend the paradigm to topology optimization. In Stokes or Brinkman flow, GP priors on state and design variables with shared neural means enable simultaneous, meshfree optimization of dissipated power or compliance, with boundary conditions enforced by conditioning and design constraints handled through training losses (Yousefpour et al., 2024, Yousefpour et al., 18 Mar 2025, Sun et al., 14 Jul 2025). These studies emphasize sharp interfaces, curriculum training, and the ability to generate super-resolution topologies. Flight-test analysis provides another mean-function-based use case: a GP with Morelli-model prior estimates pitching moment coefficient from arbitrary maneuvers and then differentiates the posterior mean to recover short-period stability derivatives and modal properties of the T-38 (Harp et al., 2 Jan 2025).

6. Misconceptions, limitations, and open directions

A persistent misconception is that “physics-informed” always means exact enforcement of the full governing equations. In practice, the literature spans exact operator priors, residual penalties, variational approximations, and informative means or kernels. The GP smoothing used to stabilize PINNs under noisy boundary data is physics-adjacent rather than a full PDE-constrained GP, because the GP does not encode the PDE itself (Bajaj et al., 2021). Conversely, linear-PDE operator regression and LODE-GPs do encode the governing equations directly in the prior or observation operator (Pförtner et al., 2022, Tebbe et al., 2024).

Another misconception is that Gaussianity automatically makes nonlinear physics easy. The operator-theoretic formulas are exact for bounded linear operators, but nonlinear PDEs generally require linearization, residual likelihoods, Laplace-type approximations, or stochastic variational inference (Long et al., 2022, Hamelijnck et al., 2024). This creates a methodological divide between exact linear theory and approximate nonlinear inference. A plausible implication is that comparisons across papers are often shaped more by the chosen approximation regime than by the GP prior alone.

Scalability remains a central limitation. Exact GP updates are cubic in the number of observations, and multi-output, derivative-rich models can become large even before hyperparameter learning or MCMC is considered (Hanuka et al., 2020, Kavrakov et al., 2024). State-space formulations, inducing points, structured kernels, and sparse GP approximations partially alleviate this, but each introduces additional design choices and potential approximation error (Long et al., 2022, Hamelijnck et al., 2024).

Model misspecification is equally consequential. The flight-test framework intentionally uses an A-7E-derived Morelli mean on T-38 data to show that the GP can correct an imperfect prior, but prediction in sparse regions still depends on that prior (Harp et al., 2 Jan 2025). CALPHAD- and VEC-informed alloy classifiers likewise rely on conservative prior probabilities to hedge against physics-model error (Hardcastle et al., 17 Feb 2025). CoPhIK formalizes this trade-off: introducing a discrepancy GP can improve accuracy when the physics model is biased, but the discrepancy may weaken strict physical constraint preservation unless its kernel is chosen compatibly with the operator (Yang et al., 2018).

Current directions suggest a gradual unification of themes that were once separate. Variational state-space GPs connect physics-informed GP inference to Kalman methods and large spatio-temporal data (Hamelijnck et al., 2024). ODE-constrained control brings GP priors into control-as-inference and safe MPC (Tebbe et al., 2024, Tebbe et al., 20 Nov 2025). Topology optimization couples GP conditioning with learned mean functions and continuation-like training (Yousefpour et al., 18 Mar 2025, Sun et al., 14 Jul 2025). This suggests that the field is moving from isolated application-specific constructions toward a broader view in which physical operators, structured priors, and scalable Bayesian computation are combined modularly rather than treated as competing paradigms.