Gaussian Process Regression Emulator

• Gaussian Process Regression Emulators are nonparametric surrogate models that use Gaussian processes to interpolate, extrapolate, and quantify uncertainty in expensive simulations.
• They integrate Bayesian inference, active learning, and nonstationary adaptations such as latent augmentation and composite strategies to enhance predictive accuracy.
• Scalable variants including Sparse, Deep, and local approximate GP models enable effective handling of large datasets and high-dimensional inputs in computational science.

A Gaussian Process Regression Emulator is a nonparametric, probabilistic surrogate model that uses Gaussian processes (GPs) to interpolate, extrapolate, and quantify uncertainty in the predictions of expensive deterministic or stochastic simulators. GPs define a distribution over functions and, after conditioning on observed simulator runs, yield closed-form predictive means and variances at new inputs. The GP emulator framework naturally supports Bayesian inference, active experimental design, calibration, and error quantification, making it a central tool in modern computational science, engineering, and statistics.

1. Mathematical Foundations and Standard GP Emulation

Given a function $f:\mathcal{X}\subset\mathbb{R}^p\to\mathbb{R}$ approximated via computationally intensive simulations, a GP emulator begins by specifying

$$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$

where $m$ is a mean function (typically $m(x)=h(x)^\top\beta$ for a basis $h(x)$ and coefficients $\beta$) and $C$ is a covariance kernel, usually taken as stationary (e.g., squared-exponential, Matérn).

For a set of $n$ training points $X = [x_1,\ldots,x_n]$ and outputs $f = (f_1,\ldots,f_n)^\top$, the GP prior yields the joint Gaussian

$$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$0

where $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$1 is the $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$2 design matrix with rows $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$3, and $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$4 is the kernel matrix parameterized by hyperparameters $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$5.

The posterior predictive mean and variance at a new $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$6 are

$$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$7

where $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$8 is the vector $$f(\cdot)\sim \mathcal{GP}\left(m(\cdot), \, C(\cdot, \cdot)\right)$$9 and $m$0 is the generalized least squares estimator. Marginalizing $m$1 and $m$2 (with conjugate priors) produces a Student-$m$3 posterior for predictions (Montagna et al., 2013, Schaechtle et al., 2015).

Once trained, the emulator offers efficient, uncertainty-aware predictions at any untried input and can serve as a surrogate for expensive simulation calls in downstream analyses such as optimization or calibration (Stetzler et al., 2022).

2. Nonstationary and Composite GP Emulators

Standard GP emulators with stationary covariance kernels perform poorly when the response surface is nonstationary—exhibiting sharp transitions, discontinuities, or highly localized features. Nonstationary extensions are critical for such regimes.

Nested/Latent-Augmented GP (Nonstationary GP):

Montagna & Tokdar (Montagna et al., 2013) propose augmenting the input space with a learned latent coordinate $m$4, where $m$5 is itself drawn from a second GP:

• The main GP for $m$6 is then indexed by $m$7,
• The covariance incorporates both the original and latent dimensions.

Formally, letting $m$8 denote the main kernel and $m$9 that of the latent process, joint inference is performed over $m(x)=h(x)^\top\beta$0: $m(x)=h(x)^\top\beta$1 Inference is conducted via particle learning, allowing for fully Bayesian sequential updating, closed-form marginalization over linear trend and scale, and marginal Student-$m(x)=h(x)^\top\beta$2 or Gaussian predictions at new points. The latent augmentation stretches the input space near sharp features, maintaining high predictive variance and driving design points toward complex regimes. Empirically, this approach reduces RMSE by 25–40% and delivers sharper error bands than stationary GPs or treed GPs (Montagna et al., 2013).

Composite Gaussian Process (CGP):

Ba & Joseph (Ba et al., 2013) present a GP model decomposed into an additive sum of two independent GPs:

• A global smooth trend $m(x)=h(x)^\top\beta$3
• A local detail process $m(x)=h(x)^\top\beta$4

Optionally, the local variance can vary spatially: $m(x)=h(x)^\top\beta$5. This allows the emulator to adapt both mean and local signal variance, improving stability and accuracy for sparse or heteroscedastic data. CGPs empirically reduce RMSPE by 20–40% and yield prediction intervals with up to 50% shorter interval-score compared to stationary GPs (Ba et al., 2013).

3. Local, High-dimensional, and Functional Emulation

Local Approximate GP (laGP):

To address computational bottlenecks for $m(x)=h(x)^\top\beta$6, local GPs fit a distinct, small-dataset GP tailored to each prediction point. The local subdesign is selected greedily to maximally reduce predictive variance at the point of interest. Sung, Gramacy & Haaland (Sung et al., 2016) introduce algorithmic accelerations such as:

• Maximum-distance screening (using a KD-tree to ignore distant candidates).
• Feature-approximation (random Fourier or Nyström features with LSH queries to reduce high-dimensional searches).

Applied to datasets of size $m(x)=h(x)^\top\beta$7, these techniques yield speedups of 10–20× over exhaustive search with negligible variance increase (<2%) (Sung et al., 2016).

High-dimensional Structure:

• Clustered Active Subspace Local GP: For problems where different low-dimensional structures dominate in different parameter regimes, CAS-LGP first clusters the input space based on gradient information and then fits a local, low-dimensional GP emulator within each cluster. This approach combines local projection-based dimension reduction with adaptive clustering for highly complex surfaces (Xiong et al., 2020).
• ANOVA–GP Decomposition: Decomposing the simulator response into low-dimensional ANOVA terms, then performing PCA-and-GP emulation on each, yields efficient surrogates for high-dimensional (dozens to hundreds of inputs) models, achieving order-of-magnitude error improvements over global PCA+GP surrogates (Chen et al., 2019).
• Mesh-clustered GP for PDEs: For parameterized PDEs solved on spatial meshes, node-wise GP priors can be clustered (Dirichlet process) to share hyperparameters among similar regions, giving interpretable, accurate, and computationally feasible surrogates (Sung et al., 2023).

Functional and Tensorial Outputs:

• Tensor-variate GP Regression: For matrix or higher-order array outputs, a TVGP model with Kronecker-structured covariance exploits output index correlations. Outer Product Emulator (OPE) exploits output similarity, while Parallel Partial Emulator (PPE) fits independent GPs to each output. OPE yields lower predictive RMSE when outputs are strongly correlated, whereas PPE is computationally scalable and handles heteroskedasticity (Semochkina et al., 14 Feb 2025).

Dynamic Simulator Emulation:

• One-step-ahead dynamic GPs model the system's discrete flow map, propagating uncertainty through time. Recent advances provide analytic moment-propagation schemes (linked GP) for these models, which are orders of magnitude faster than Monte Carlo alternatives and preserve error bars even in chaotic regimes (Heo, 26 Mar 2025, Mohammadi et al., 2018).

4. Large-scale, Deep, and Generalized GP Emulators

Sparse and Deep GP Methods:

To scale GP emulators to $m(x)=h(x)^\top\beta$8–$m(x)=h(x)^\top\beta$9 with moderate loss of accuracy:

• Sparse Variational GP (SVGP): Inducing-point models approximate the full GP posterior by conditioning on a small set of pseudo-inputs, yielding computational cost $h(x)$0 with $h(x)$1. Variational inference maximizes a tractable lower bound (ELBO) (Stetzler et al., 2022).
• Stochastic Variational GP: Minibatch-based SVGP reduces per-iteration cost to $h(x)$2 for batch-size $h(x)$3, supporting large $h(x)$4 (Stetzler et al., 2022).
• Deep Kernel GP (DKLGP): Pass inputs through a trained neural-network $h(x)$5, then apply a stationary GP kernel in latent space, enabling the GP to address complex nonstationary structure (Stetzler et al., 2022). For functional data and high-dimensional outputs, scalable local-approximate GP emulators with low-rank decompositions (basis SVD or PCA followed by local GPs for latent weights) deliver rapid prediction and calibration in applications with $h(x)$6 (Hutchings et al., 2024).

Generalized Deep GP (GDGP) Framework:

For heterogeneous data types and responses—including Poisson, negative-binomial, categorical, or heteroscedastic Gaussian outputs—a generalized DGP can propagate uncertainty through multiple GP layers and link functions. Scalable inference is achieved by the Vecchia approximation (sparse precision matrices via low-order conditioning) and efficient sampling strategies. GDGP unifies probabilistic surrogate modeling across a diverse range of simulator outputs (Ming et al., 25 Mar 2026).

5. Uncertainty Quantification, Calibration, and Active Learning

Predictive Uncertainty:

GP emulators provide full posterior predictive distributions, yielding local uncertainty quantification and error bars. When hyperparameters are estimated by empirical Bayes (marginal likelihood maximization), convergence rates are provably unaffected under mild regularity conditions (Teckentrup, 2019).

Bayesian Calibration:

GP surrogates enable Bayesian parameter inference by replacing the expensive simulator in the likelihood with the emulator's predictive mean and variance, reducing the computational burden in MCMC and enabling practical uncertainty propagation. Modular calibration schemes (separating surrogate training from parameter inference) are efficient and maintain statistical validity (Stetzler et al., 2022, Hutchings et al., 2024).

Active Learning and Experimental Design:

Acquisition criteria such as maximum predictive variance (Active Learning MacKay) or expected improvement are used for sequential design. The nonstationary GP maintains high variance near regions with rapid surface changes, thus adaptively refining emulator fidelity in the most informative regions. Multi-output extensions (e.g., AMOGAPE) incorporate gradient and joint-output metrics to concentrate design points in regions of high response curvature and low data density (Svendsen et al., 2019).

6. Robustness, Implementation, and Best Practices

Parameter Estimation:

Robust parameter estimation is critical for stable GP emulation. Marginal posterior mode estimators with objective or jointly-robust priors are designed to prevent collapsed or degenerate kernels. Identification of "inert" (inactive) inputs can be performed at no extra computational cost via normalized inverse-range parameters at the mode (Gu et al., 2018).

Probabilistic Programming and Automation:

Embedding GP regression as a statistical memoizer in probabilistic programming languages allows the dynamic, online construction of surrogates that adapt as new data is observed, supports fully hierarchical hyperparameter inference, and enables Bayesian structure learning over composite kernels (Schaechtle et al., 2015).

Implementation Tooling:

Numerous R and Python packages implement the above methodologies. Tools such as RobustGaSP, FlaGP, dynemu, and dgpsi provide robust fitting, support multi-output and functional data, modular calibration, and efficient high-dimensional inference (Gu et al., 2018, Hutchings et al., 2024, Heo, 26 Mar 2025, Ming et al., 25 Mar 2026).

Model/Method	Complexity/Scalability	Regime
Full GP	$h(x)$7 time, $h(x)$8 mem	$h(x)$9
SVGP/Sparse GP	$\beta$0	$\beta$1 large, $\beta$2
Local GP (laGP)	$\beta$3 per test pt	Massive $\beta$4, local approx
Functional PCA+GP	$\beta$5 ($\beta$6 basis), local GP	Large outputs, $\beta$7 large
GDGP/Vecchia	$\beta$8, $\beta$9 neighbors	Large $C$0, non-Gauss resp
OPE (Kronecker)	$C$1	Matrix/array outputs

7. Empirical Performance and Theoretical Guarantees

Nonstationary and composite GPs consistently demonstrate significant reductions in predictive RMSE in scenarios with localized features or nonhomogeneous variance (Montagna et al., 2013, Ba et al., 2013). Local approximation and inducing-point models yield practical emulators for medium and large $C$2, trading off small increases in RMSE for the addressability of large data volumes (Stetzler et al., 2022, Hutchings et al., 2024). Convergence theory indicates that empirical-Bayes learning of hyperparameters does not degrade the asymptotic accuracy of GP emulators for function approximation or downstream posterior inference (Teckentrup, 2019).

The cumulative impact is that GP regression emulators—by leveraging model-based uncertainty quantification, scalable implementations, and extensible modeling (nonstationary, high-dimensional, multi-output)—define the computational and statistical foundation for modern surrogate modeling in simulation-based science and engineering.