A Survey of Constrained Gaussian Process Regression: Approaches and Implementation Challenges

Laura Swiler; Mamikon Gulian; Ari Frankel; Cosmin Safta; John Jakeman

Source-linked AI summary

A Survey of Constrained Gaussian Process Regression: Approaches and Implementation Challenges

Laura Swiler, Mamikon Gulian, Ari Frankel, Cosmin Safta, John Jakeman

arXiv:2006.09319v3Jun 16, 2020cs.LGmath.STstat.ML

TL;DR

Scientific applications often have limited or expensive data, motivating constrained GPR methods that incorporate physical and other prior information. The survey classifies approaches for bound, shape, PDE, and boundary constraints, compares their implementations, and discusses computational challenges, including multidimensional scaling and sampling limitations.

Problem
Scientific machine learning often must build models from limited or expensive data while respecting physical principles and prior information.
Method
The survey organizes constrained GPR approaches by constraint type and enforcement strategy, covering bounds, monotonicity, convexity, linear PDEs, and boundary conditions.
Results
The survey provides an overview of major constrained-GPR categories, implementation details, numerical examples, and computational challenges.
Takeaways & Limitations
Constraints can supplement limited data and help models remain consistent with prior physical knowledge, including in extrapolatory regions.
Takeaways & Limitations
The survey is not exhaustive, and multidimensional spline inference can scale as O(M^(3d)) with M knots per dimension.

Abstract

from arXiv · show

Gaussian process regression is a popular Bayesian framework for surrogate modeling of expensive data sources. As part of a broader effort in scientific machine learning, many recent works have incorporated physical constraints or other a priori information within Gaussian process regression to supplement limited data and regularize the behavior of the model. We provide an overview and survey of several classes of Gaussian process constraints, including positivity or bound constraints, monotonicity and convexity constraints, differential equation constraints provided by linear PDEs, and boundary condition constraints. We compare the strategies behind each approach as well as the differences in implementation, concluding with a discussion of the computational challenges introduced by constraints.

1 Introduction

The survey examines constrained Gaussian process regression as a way to incorporate physical principles and other prior information when scientific data are limited or expensive. It classifies constraint types and enforcement strategies, while emphasizing implementation guidance and computational challenges.

Scientific machine learning often lacks abundant training data because computational and physical experiments are expensive.
Physical principles and a priori information can supplement limited data and regularize model behavior.
Gaussian process regression is useful for expensive data sources because accurate surrogates can often be built from tens to hundreds of training pairs.
The survey covers bound, monotonicity, convexity, linear PDE, and boundary-condition constraints in constrained GPR.
Its taxonomy is designed to help practitioners identify relevant constraints and implementation approaches, with cross-references for overlapping methodologies.
The survey is selective rather than exhaustive because the constrained-GP literature is extensive and rapidly expanding.

2 Gaussian Process Regression

Gaussian process regression models distributions over functions using mean and covariance functions, then updates them with observed data through Bayesian inference. Standard inference is analytically convenient under Gaussian likelihoods but faces cubic scaling, conditioning, and nonconvex hyperparameter optimization.

A Gaussian process defines finite collections of function values with joint multivariate Gaussian distributions determined by mean and covariance functions.
The covariance kernel must be symmetric and positive semidefinite, and its choice can substantially affect predictions and enable constraint-specific constructions.
GPR combines a GP prior with an observation likelihood and uses Bayes’ rule to obtain posterior predictions at new inputs.
For Gaussian likelihoods, the posterior is Gaussian, so exact computation and sampling are straightforward; non-Gaussian likelihoods generally remove this convenience.
Hyperparameters govern the covariance kernel and noise model, and training commonly estimates them by optimizing the marginal likelihood.
Inference requires inversion of an N×N matrix with O(N^3) scaling, while conditioning and nonconvex hyperparameter optimization can make computation costly and unstable.

3 Strategies for Constraints

Constrained GP methods must approximate or encode conditions that may apply globally across a continuous domain and across function realizations. The survey organizes strategies by modifying outputs, likelihoods, kernels, function representations, or data.

Global constraints create an infinite set of conditions over a continuous domain and an infinite-dimensional function space, requiring numerical approximations.
Some methods enforce constraints at finite virtual points, while others transform GP outputs to guarantee admissible predictions.
Constraint-specific covariance kernels can encode physical relationships through joint GPs and four-block covariance structures.
Spline or other function-space representations translate global constraints into finite coefficient constraints, while alternative representations can make every sampled function admissible.
Data can impose constraints implicitly when posterior predictions inherit exact or approximate properties satisfied by the supplied data.
Some constrained methods produce non-Gaussian posteriors, so the term constrained Gaussian process is technically a misnomer in those cases.

4 Bound Constraints

Bound constraints require GP predictions to remain within specified limits over a region, including nonnegativity as a special case. The survey presents three approaches for enforcing these bounds.

Bound constraints require predictions to satisfy a ≤ f(x) ≤ b over a region of interest.
Nonnegativity is the special case f ≥ 0, obtained with a = 0 and b = ∞.
The survey motivates bounds with chemical concentration regression, where predicted values should remain between 0 and 1.
Three approaches for enforcing bound constraints are presented.

4.1 Transformed Output and Likelihood

Bound constraints can be imposed by transforming GP outputs or replacing the Gaussian likelihood with a bounded non-Gaussian likelihood. Warping preserves standard GP treatment, whereas likelihood formulations require approximate inference because their posteriors are not analytically tractable.

Warping functions: Warping transforms bounded observations into unbounded values before fitting a traditional Gaussian process.The probit function maps values in [0,1] to the real line.
Warping functions: The probit warping function is the inverse cumulative distribution function of a standard normal random variable.It is commonly used for observations bounded in [0,1].
Likelihood formulations: Bounded likelihood formulations replace the Gaussian likelihood with a finite-support distribution so GP observations satisfy the constraints.Examples include truncated Gaussian and beta likelihoods.
Likelihood formulations: Truncated Gaussian and beta likelihoods produce posteriors that are not analytically tractable, requiring approximate inference.Laplace approximation and expectation propagation approximate the posterior with a multivariate Gaussian.

4.2 Discrete Constraints using Truncated Gaussian Distributions

Discrete constraint methods approximate global GP constraints by conditioning predictions at finitely many auxiliary points on a truncated multivariate normal distribution. They support location-dependent bounds and other constraints, but require costly sampling and leave the necessary number of auxiliary points uncertain.

Discrete constraints: Global continuous-domain constraints are approximated by imposing bounds at a finite set of auxiliary or virtual points.The constrained values satisfy a ≤ f(x_i) ≤ b at the selected points.
Discrete constraints: Conditioning the GP values at auxiliary points on the bounds changes their distribution from multivariate normal to truncated multivariate normal.The construction first conditions on observed data and then applies the inequality constraints.
Moment computation: In special one-dimensional cases, truncated-normal moments can be derived analytically using standardized bounds and normal density and distribution functions.The general case requires estimating moments empirically.
Computational challenges: Sampling and computing moments of the truncated multivariate normal is computationally demanding.The method described uses an expensive rejection-sampling procedure based on a modified Gibbs sampler.
Scope and limitations: The approach permits location-dependent bounds and can also enforce monotonicity and convexity when the auxiliary-point set is chosen adequately.However, the number of points needed for a prescribed global confidence level is unclear.

4.3 Constrained maximum likelihood optimization to enforce nonnegativity constraints

Constrained maximum-likelihood optimization selects GP hyperparameters subject to probabilistic nonnegativity and data-fit constraints. This preserves the Gaussian posterior form, while increasing optimization cost and reducing the need for more expensive inference methods.

Constrained optimization: Hyperparameters are chosen by constraining log-marginal-likelihood optimization so predictions satisfy nonnegativity with small violation probability at finite constraint points.The approach uses posterior mean and standard deviation constraints.
Probabilistic constraints: The method sets ϵ = 2.3%, making the posterior mean minus two standard deviations nonnegative at constraint points.This follows from Φ^-1(ϵ) = −2.
Probabilistic constraints: The optimization also constrains predictions at observed inputs to remain within ν of the observations.The reported formulation uses ν > 0 as a data-fit tolerance.
Results: Using ν = 0.03, a nonlinear interior-point solver enforced nonnegativity with high probability and significantly reduced posterior variance.These outcomes were reported for the constrained optimization problem.
Computational implications: Although constrained optimization is more expensive than unconstrained marginal-likelihood optimization, it preserves the Gaussian posterior form and avoids methods such as MCMC.The examples considered nonnegative constraints, though the framework can in principle address other constraints.

4.4 Splines

Spline-based constrained GPs reduce continuous-domain constraint enforcement to finite-dimensional constraints on spline coefficients. Piecewise linear interpolation can guarantee global bounds or monotonicity, but higher-dimensional inference suffers severe computational scaling and difficult truncated-normal sampling.

Spline formulation: Maatouk and Bay place a multivariate Gaussian prior on spline functions and impose constraints on their coefficients.The spline representation reduces the continuous GP to a finite-dimensional problem involving knot values.
Global constraint enforcement: Because piecewise linear splines form a partition of unity, bounds or monotonicity imposed on knot values hold between neighboring knots and therefore globally.This converts global continuous-domain constraints into finite-dimensional coefficient constraints.
Spline formulation: Piecewise linear splines form a C0 interpolant of knot values, with coefficients serving as the spline representation.The knots are located at x_i = i/M on [0,1].
Interpolation: Spline coefficients satisfying the interpolation system reproduce the observed data, provided the spline basis is sufficiently fine.The construction requires enough coefficients for each basis function to span no more than two data points.
Inference: The constrained knot values are inferred from a GP prior and observations, then conditioned on inequality constraints through a truncated normal distribution.Sampling from this distribution creates an additional computational challenge.
Sampling: Simple rejection sampling becomes inefficient at larger scale because rejection probability increases exponentially with the number of splines.Imprecise mode evaluation can further reduce acceptance.
Multidimensional setting: In d dimensions with M knots per dimension, the number of coefficients scales as M^d and inference cost as O(M^3d).The construction may also require M ≥ N, while rejection and Gibbs sampling efficiency deteriorate as problem size grows.

5 Monotonicity Constraints

Monotonicity constraints are imposed by requiring selected derivatives of a Gaussian process to be nonnegative, using either constrained likelihoods, truncated derivative distributions, or globally constrained splines. These methods introduce implementation and computational challenges, including virtual-point selection, non-Gaussian inference, and higher-dimensional convexity limitations.

5.1 Constrained Likelihood with Derivative Information: A probit likelihood converts derivative values into probabilities that the derivative is positive, allowing soft monotonicity constraints rather than a hard step function.The reviewed formulation uses ν = 10^-6, making the likelihood nearly step-like while remaining more tolerant of error.
5.1 Constrained Likelihood with Derivative Information: Monotonicity can be enforced by modeling GP derivatives jointly with function values and constraining derivative observations to be positive at selected virtual points.The expanded covariance includes function-function, derivative-derivative, and function-derivative blocks.
5.1 Constrained Likelihood with Derivative Information: The resulting posterior is analytically intractable because the derivative likelihood is non-Gaussian, so inference uses expectation propagation or MCMC sampling.The survey’s example applies MCMC to a challenging non-monotonic dataset subject to a monotonic GP requirement.
5.1 Constrained Likelihood with Derivative Information: Virtual operating points determine where derivative constraints are imposed, but grid-based placement suffers from the curse of dimensionality and may require adaptive selection.A suggested strategy adds points where the current GP most likely has negative derivatives.
5.2 Truncated Derivative Distributions and 5.3 Monotonic splines: Alternative methods constrain derivative distributions at finite virtual points through truncated densities or extend spline constraints globally using finite-difference matrices.Spline-based monotonicity uses neighboring knot differences to impose the constraint over the represented domain.
5.4 Convexity Constraints: Convexity becomes more difficult beyond one dimension because it requires a positive-semidefinite Hessian rather than a simple linear derivative inequality.Sylvester’s criterion yields equivalent inequalities involving Hessian minors.

6 Differential Equation Constraints

Differential-equation constraints enter GPR either through joint modeling of solution and forcing data or through kernels whose structure enforces the constraint globally. These approaches improve constrained prediction but introduce implementation and computational trade-offs.

6.1 Block Covariance Kernel: Joint GPR uses observations of both u and the forcing f through a four-block covariance matrix containing solution, forcing, and cross-covariance terms.Hyperparameters can be trained by assembling the block matrix and minimizing the negative log-marginal likelihood, with separate noise terms when measurements are noisy.
Operator-based formulation: Linear differential operators transform a GP for u into a GP for Lu, with covariance L_xL_x′k(x,x′) when the required regularity conditions hold.This provides the basis for modeling PDE-constrained solutions and forcing functions.
6.1 Block Covariance Kernel: Additional forcing observations can improve prediction accuracy and reduce variance for u, especially where direct solution data are sparse.In the example, ten u observations on [0.2, 0.8] plus ten PDE-data observations across [0, 1] improve reconstruction outside the observed subinterval.
Implementation challenges: Closed-form operator-applied kernels simplify implementation, whereas numerical differentiation can create covariance ill-conditioning, prediction artifacts, or failure as datasets grow.Diagonal regularization can reduce ill-conditioning but may reduce regression accuracy.
6.2 Transformed Covariance Kernel: Transformed covariance kernels enforce differential constraints globally through kernel structure without auxiliary constraint data, but require solving an operator equation tailored to the constraint.The method avoids the four-block covariance burden but has more limited applicability and specialized implementation.
6.2 Transformed Covariance Kernel: Solving the operator equation can be difficult and case-dependent, with integral-operator constructions making kernel evaluation harder than block-covariance computation.Specialized basis kernels can improve reconstructions for PDEs such as Laplace, heat, and Helmholtz equations, but require analytical construction and are nonstationary.
6.3 Empirical Mean and Covariance: Empirical mean-and-covariance approaches preserve deterministic linear-operator constraints only within a certain error and bypass pointwise kernel prediction.They therefore cannot interpolate the field at arbitrary points between available stochastic realizations.

7 Boundary Condition Constraints

Boundary information can be incorporated intrinsically into Gaussian processes, including through Laplacian spectral expansions that enforce zero Dirichlet conditions. The approach requires computing the domain’s spectrum and has extensions and limitations tied to boundary type and kernel isotropy.

Boundary values can be incorporated intrinsically into Gaussian processes rather than treated only as scattered data, supporting accuracy and well-posed modeling.
7.1 Spectral Expansion Approach: Spectral expansion methods construct covariance kernels from Laplacian eigenfunctions to enforce zero Dirichlet boundary values.
7.1 Spectral Expansion Approach: The resulting covariance is zero when either input lies on the boundary, so predictions satisfy zero boundary conditions when the GP mean does.
7.1 Spectral Expansion Approach: Implementing the method requires computing the Dirichlet eigenvalues and eigenfunctions, analytically on basic domains or numerically on general domains.
7.1 Spectral Expansion Approach: The spectral formulation supports efficient GP regression and hyperparameter optimization through covariance diagonalization and derivatives of the spectral density.
7.3 Extensions: The surveyed boundary method is primarily developed for periodic or zero Dirichlet conditions and isotropic kernels, while nonzero Dirichlet, Neumann, and nonisotropic extensions are suggested.

8 Computational Considerations

Constrained Gaussian processes are more expensive than unconstrained ones because constraints enlarge the data, posterior, and numerical integration problems. The survey reviews truncated-normal methods, constrained hyperparameter estimation, and scalable linear-algebra strategies, while noting important open questions.

Constrained GP construction is typically more expensive because derivative, bound, and virtual constraint data enlarge the computational problem.
8.1 The Truncated Multivariate Normal Distribution: The truncated normal conditions a Gaussian random variable on rectangular bounds, but normalization and sampling become difficult for general covariance matrices and dimensions.
8.1 The Truncated Multivariate Normal Distribution: Simple Monte Carlo methods scale poorly in high dimensions, while MCMC for spline knots may converge slowly as the number of splines grows.
8.2 Hyperparameter Estimation: Constrained maximum-likelihood estimation adds Gaussian orthant probabilities to the objective, and their numerical approximation can affect optimization accuracy.
8.2 Hyperparameter Estimation: For limited data, cMLE was found to estimate hyperparameters and confidence intervals more accurately than MLE, whereas constraint effects diminish under dense fixed-domain data.
8.3 Reducing Computational Cost: Standard GP inference scales as O(N^3), and virtual constraint points can further enlarge covariance matrices, motivating subset, low-rank, hierarchical, and numerical-linear-algebra approaches.
8.3 Reducing Computational Cost: For tensor grids, Kronecker methods reduce storage from O(N^2) to O(dN^(2/d)) and inversion complexity from O(N^3) to O(dN^((d+1)/d)).

9 Conclusion

The survey organizes constrained Gaussian processes by constraint type and enforcement strategy, covering bounds, shape constraints, linear differential equations, and boundary values. It concludes that these methods broaden physically consistent modeling but introduce substantial computational challenges, and that the non-exhaustive literature continues to expand.

The survey covers positivity and bounds, monotonicity and convexity, linear differential equations, and boundary-value constraints.
Constraints can be enforced implicitly through data, by constructing a tailored sample space or covariance kernel, or by modifying GP outputs or likelihoods.
Constraint enforcement may be global, imposed at finite virtual points, or approximate, providing key distinctions among surveyed approaches.
Constrained GPs introduce challenges including constrained covariance construction, truncated-normal sampling, non-Gaussian posteriors, expanded virtual-data matrices, and eigenvalue computation.
The survey is non-exhaustive and aims to catalog common approaches, guide method selection, and identify computational challenges and possible responses.