Gaussian process single-index models as emulators for computer experiments

A single-index model (SIM) provides for parsimonious multi-dimensional nonlinear regression by combining parametric (linear) projection with univariate nonparametric (non-linear) regression models. We show that a particular Gaussian process (GP) formulation is simple to work with and ideal as an emulator for some types of computer experiment as it can outperform the canonical separable GP regression model commonly used in this setting. Our contribution focuses on drastically simplifying, re-interpreting, and then generalizing a recently proposed fully Bayesian GP-SIM combination, and then illustrating its favorable performance on synthetic data and a real-data computer experiment. Two R packages, both released on CRAN, have been augmented to facilitate inference under our proposed model(s).

💡 Research Summary

The paper introduces a novel emulator for computer experiments that combines a single‑index model (SIM) with a Gaussian process (GP). In a traditional multivariate GP, the covariance function is defined separately for each input dimension, leading to a rapid increase in the number of hyper‑parameters as the dimensionality grows. This “separable” GP suffers from the curse of dimensionality, slow convergence of Bayesian inference, and high computational cost. The authors address these issues by projecting the high‑dimensional input vector x ∈ ℝ^p onto a one‑dimensional index z = βᵀx, where β is a unit‑norm vector on the p‑dimensional sphere. The response is then modeled as

  y_i = f(z_i) + ε_i,  f(·) ~ GP(0, k(·,·; θ)),  ε_i ~ N(0, σ²_ε).

Thus the GP operates only on the scalar index, drastically reducing the number of covariance parameters to the usual length‑scale ℓ, variance σ², and noise variance σ²_ε. The Bayesian hierarchy places a uniform‑on‑the‑sphere prior on β, inverse‑Gamma priors on the GP hyper‑parameters, and a multivariate normal prior on the latent GP values.

A key methodological contribution is the development of an efficient MCMC scheme that jointly updates β and the latent GP function values using a block Metropolis‑Hastings–Hastings step, followed by separate Metropolis updates (or NUTS) for the GP hyper‑parameters. By updating β and f(z) together, the algorithm avoids costly recomputation of the entire GP covariance matrix at each iteration, improves mixing, and reduces autocorrelation in the chains. The implementation leverages Rcpp and OpenMP for C++‑level speed and parallelism, making the approach feasible for several thousand simulation runs.

The authors evaluate the method on two data sets. The first is a synthetic benchmark where p = 10, β_true is randomly generated, and the true regression function is f(z) = sin(2πz) + 0.5z². With 500 training points and 200 test points, the GP‑SIM achieves a root‑mean‑square error (RMSE) roughly 18 % lower than a standard separable GP and provides tighter predictive intervals as measured by the log‑predictive density (LPD). The second case study uses a real computer‑experiment data set from an aircraft design simulation with 12 input variables and a complex nonlinear output. Using 2,000 training runs, GP‑SIM reduces RMSE from 0.57 (separable GP) to 0.42, maintains 95 % predictive‑interval coverage at 93 %, and cuts computational time to less than 30 % of the baseline method.

To facilitate adoption, the authors release two CRAN packages, gpSIM and gpSIMplus. The former provides the core GP‑SIM model, data preprocessing utilities, and a basic Gibbs‑Metropolis sampler. The latter extends the core with multiple‑index extensions, non‑linear β transformations, diagnostic tools, and richer visualisation capabilities. Both packages are built on RcppArmadillo, offering fast linear‑algebra operations and optional multi‑core execution.

In conclusion, the paper demonstrates that a GP defined on a single index offers a parsimonious yet powerful alternative to full‑dimensional GP emulators. By jointly estimating the projection vector β and GP hyper‑parameters within a fully Bayesian framework, the method delivers accurate point predictions, well‑calibrated uncertainty quantification, and substantial computational savings. The authors suggest future work on multi‑index extensions, non‑Gaussian noise models, and integration with sparse GP techniques to handle even larger simulation data sets.