Bayesian Deep Gaussian Processes using MCMC

Performs Bayesian posterior inference for deep Gaussian processes following Sauer, Gramacy, and Higdon (2023, ). See Sauer (2023, < http://hdl.handle.net/10919/114845>) for comprehensive methodological details and < https://bitbucket.org/gramacylab/deepgp-ex/> for a variety of coding examples. Models are trained through MCMC including elliptical slice sampling of latent Gaussian layers and Metropolis-Hastings sampling of kernel hyperparameters. Gradient-enhancement and gradient predictions are offered following Booth (2026, ). Vecchia approximation for faster computation is implemented following Sauer, Cooper, and Gramacy (2023, ). Optional monotonic warpings are implemented following Barnett et al. (2025, ). Downstream tasks include sequential design through active learning Cohn/integrated mean squared error (ALC/IMSE; Sauer, Gramacy, and Higdon, 2023), optimization through expected improvement (EI; Gramacy, Sauer, and Wycoff, 2022, ), and contour location through entropy (Booth, Renganathan, and Gramacy, 2025, ). Models extend up to three layers deep; a one layer model is equivalent to typical Gaussian process regression. Incorporates OpenMP and SNOW parallelization and utilizes C/C++ under the hood.

deepgp Package

Maintainer: Annie S. Booth ([email protected])

Performs Bayesian posterior inference for deep Gaussian processes following Sauer, Gramacy, and Higdon (2023). See Sauer (2023) for comprehensive methodological details and https://bitbucket.org/gramacylab/deepgp-ex/ for a variety of coding examples. Models are trained through MCMC including elliptical slice sampling of latent Gaussian layers and Metropolis-Hastings sampling of kernel hyperparameters. Gradient-enhancement and gradient predictions are offered following Booth (2026). Vecchia approximation for faster computation is implemented following Sauer, Cooper, and Gramacy (2023). Optional monotonic warpings are implemented following Barnett et al. (2025). Downstream tasks include sequential design through active learning Cohn/integrated mean squared error (ALC/IMSE; Sauer, Gramacy, and Higdon, 2023), optimization through expected improvement (EI; Gramacy, Sauer, and Wycoff, 2022), and contour location through entropy (Booth, Renganathan, and Gramacy, 2025). Models extend up to three layers deep; a one layer model is equivalent to typical Gaussian process regression.
Incorporates OpenMP and SNOW parallelization and utilizes C/C++ under the hood.

Run help("deepgp-package") or help(package = "deepgp") for more information.

References

Sauer, A. (2023). Deep Gaussian process surrogates for computer experiments. Ph.D. Dissertation, Department of Statistics, Virginia Polytechnic Institute and State University. http://hdl.handle.net/10919/114845

Booth, A. S. (2026). Gradient-enhancement and gradient predictions for deep Gaussian process modeling of expensive computer experiments. arXiv:2512.18066

Sauer, A., Gramacy, R.B., & Higdon, D. (2023). Active learning for deep Gaussian process surrogates. Technometrics, 65, 4-18. arXiv:2012.08015

Sauer, A., Cooper, A., & Gramacy, R. B. (2023). Vecchia-approximated deep Gaussian processes for computer experiments. Journal of Computational and Graphical Statistics, 32(3), 824-837. arXiv:2204.02904

Gramacy, R. B., Sauer, A. & Wycoff, N. (2022). Triangulation candidates for Bayesian optimization. Advances in Neural Information Processing Systems (NeurIPS), 35, 35933-35945. arXiv:2112.07457

Booth, A., Renganathan, S. A. & Gramacy, R. B. (2025). Contour location for reliability in airfoil simulation experiments using deep Gaussian processes. Annals of Applied Statistics, 19(1), 191-211. arXiv:2308.04420

Barnett, S., Beesley, L. J., Booth, A. S., Gramacy, R. B., & Osthus D. (2025). Monotonic warpings for additive and deep Gaussian processes. Statistics and Computing, 35(3), 65. arXiv:2408.01540

Version History

What's new in version 1.2.3?

• Minor bug fix when the package is installed without OpenMP

What's new in version 1.2.2?

• Gradient features (enhancement and predictions) now support cov = "matern" with v = 2.5
• Bug fix for predict with one-layer GP when return_all = TRUE and cores = 1. Thanks Steven Barnett!

What's new in version 1.2.1?

• Replaced span with arma::span throughout C++ code to avoid namespace conflict when compiled with g++15 in C++20 mode.
• Reduced the computation time required to compile the vignette by lowering the number of MCMC samples used in each illustrative example.

What's new in version 1.2.0?

• Gradients! All gradient implementations come with optional Vecchia approximation.
  • Gradient-enhancement is enabled through the dydx argument in fit_one_layer and fit_two_layer (gradient-enhancement is not yet offered for three layer models). Gradient-enhancement requires the exp2 kernel.
  • Setting grad = TRUE in the predict or post_sample functions will return posterior predictions/samples of the gradient (one and two layer only). Again, this requires the exp2 kernel.
• A new post_sample function offers joint posterior draws with optional Vecchia approximation (compared to the predict functions that return summarized posterior moments).
• A new to_vec function will convert a non-Vecchia fit to a Vecchia approximated fit. This is helpful when training data sizes are relatively small (not requiring Vecchia) but testing data sizes are large (requiring Vecchia).
• Additional user controls.
  • The user may now specify the scale (tau2_w and/or tau2_z) value to use on latent layers using the settings argument in fit_two_layer and fit_three_layer.
  • The user may now specify the number of cores to use for Vecchia OpenMP parallelization using the optional cores argument in any of the fit functions. If not specified, this defaults to min(4, maxcores - 1).
• Monotonic warpings are no longer scaled to [0, 1], and they now use an isotropic lengthscale on the outer layer.
• Internal code improvements (external interface is unchanged):
  • Vecchia approximation objects and implementation have been streamlined.
  • The ordering of covariance parameters as arguments to internal functions has been unified. They are now always presented in the order: tau2, theta, g, v.
  • • If true_g is specified, the returned fit object no longer stores an entire vector of g values. The specified g value is only stored once.
• User-specified hyperparameter priors are now set using a list with the parameter name first, e.g., settings = list(theta = list(alpha = 1, beta = 1), g = list(alpha = 1, beta = 1)).
• New warning message if cores used for SNOW parallelization is greater than the number of provided nmcmc iterations. If so, cores is overwritten with cores = nmcmc.
• Minor bug fix for use of rand_mvn_vec function when pmx = TRUE and the scale on the latent layer was not equal to 1.
• Bug fix for storage of tau2 values within Gibbs sampling. Thanks Parul Patil!

What's new in version 1.1.3?

• Option to force monotonic warpings in the two-layer DGP with the argument monowarp = TRUE to fit_two_layer. Monotonic warpings trigger separable lengthscales on the outer layer.
• Updated default prior values on lengthscales and nugget for noisy settings (when true_g = NULL)
• Minor bug fix in Gibbs updating of separable lengthscale sampling in fit_one_layer
• Some improvements to default plotting
• Updated package examples and vignette

What's new in version 1.1.2?

• Option for user to specify ordering for Vecchia approximation (through ord argument in fit functions)
• lite = TRUE predictions have been sped up
  • bypassing the cov(t(mu_t)) computation altogether (this is only necessary for lite = FALSE)
  • removing d_new calculations
  • using diag_quad_mat Cpp function more often
• Expected improvement is now available for Vecchia-approximated fits
• Internally, predict functions have been consolidated (removing hundreds of lines of redundant code)
• Removed internal clean_prediction function as it was no longer needed
• Minor bug fixes
  • Fixed error in fit_one_layer with vecchia = TRUE and sep = TRUE caused by the arma::mat covmat initialization in the vecchia.cpp file
  • Fixed error in predict.dgp2 with return_all = TRUE (replaced out with object - thanks Steven Barnett!)
  • Fixed storage of ll in continue functions (thanks Sebastien Coube!)

What's new in version 1.1.1?

• Entropy calculations for contour locating sequential designs are offered through the specification of an entropy_limit in any of the predict functions.
• In posterior predictions, there is now an option to return point-wise mean and variance estimates for all MCMC samples through the specification of return_all = TRUE.
• To save on memory and storage, predict functions no longer return s2_smooth or Sigma_smooth. If desired, these quantities may be calculated by subtracting tau2 * g from the diagonal.
• The vecchia = TRUE option may now utilize either the Matern (cov = "matern") or squared exponential kernel (cov = "exp2"").
• Performance improvements for cores = 1 in predict, ALC, and IMSE functions (helps to avoid a SNOW conflict when running multiple instances on the same machine).
• Fit functions now return the outer log likelihood value along with MCMC samples.
  Used in trace plots to assess burn-in.
• In fit_two_layer, the intermediate latent layer may now have either a prior mean of zero (default) or a prior mean equal to x (pmx = TRUE). If pmx is set to a constant, this will be the scale parameter on the inner Gaussian layer.

What's new in version 1.1.0?

• Package vignette
• Option to specify sep = TRUE in fit_one_layer to fit a GP with separable/anisotropic lengthscales
• Default cores in predict are now 1 (this avoids a conflict when running multiple sessions simultaneously on a single machine)

What's new in version 1.0.1?

• Minor bug fixes/improvements
• New warning message when OpenMP parallelization is not utilized for the Vecchia approximation. This happens when the package is downloaded from CRAN on a Mac. To set up OpenMP, download package source and compile with gcc/g++ instead of clang.

What's new in version 1.0.0?

• Models may now leverage the Vecchia approximation (through the specification of vecchia = TRUE in fit functions) for faster computation. The speed of this implementation relies on OpenMP parallelization (make sure the -fopenmp flag is present with package installation).
• SNOW parallelization now uses less memory/storage
• tau2 is now calculated at the time of MCMC, not at the time of prediction.
  This avoids some extra calculations.

What's new in version 0.3.0?

• The Matern kernel is now the default covariance. The smoothness parameter is user-adjustable but must be either v = 0.5, v = 1.5, or v = 2.5 (default). The squared exponential kernel is still required for use with ALC and IMSE (set cov = "exp2" in fit functions).
• Expected improvement (EI) may now be computed concurrently with predictions. Set EI = TRUE inside predict calls. EI calculations are nugget-free and are for minimizing the response (negate y if maximization is desired).
• To save memory, hidden layer mappings used in predictions are no longer stored and returned by default. To store them, set store_latent = TRUE inside predict.