Efficient constrained Gaussian process approximation using elliptical slice sampling (ESS):

This repository contains R code and functions for implementing part of the numerical results from the paper 'Efficient Constrained Gaussian Process Approximation Using Elliptical Slice Sampling' (Maatouk et al., 2025b). The proposed \textsf{LS-ESS} approach developed in this repository can incorporate multiple shape constraints such as monotonicity, boundedness, and convexity, in the entire domain and estimate the hyperparameters. It handles one- and multi-dimensional input settings. It is important to note that the proposed \textsf{LS-ESS} algorithm has been tested in cases of extreme constraints, a task that becomes impossible for the HMC sampler. This repository contains functions for implementing the finite-dimensional Gaussian Process (GP) approximation (Maatouk and Bay, 2017; Maatouk et al., 2024b; Maatouk et al. 2025b) using different prior samplers such as Fast.LS (Maatouk et al. 2025a), LS.KLE (Maatouk et al. 2024a), and FFT-WC (Murray et al. 2010), as well as different posterior samplers such as Hamiltonian Monte Carlo (HMC) (Pakman and Paninski, 2014) and Botev's sampler (Botev, 2017). It is important to note that the proposed approach can handle higher-dimensional input spaces (d>=2) with large datasets (see the image Monotonicity 7D). It is also important to note that the posterior distribution of the proposed finite-dimensional GP approximation has been sampled using different Gibbs sampling approaches.

⚠️ Under Development: This repository is a work in progress and will continue to be updated and improved over time.

Description of the associated R files:

1. 'all_base_functions.R' file contains all the required base functions to generate the efficient constrained GP approximation using elliptical slice sampling, including functions like Fast.LS (draws a very large Gaussian vector prior based on Maatouk et al. (2025a)), LS.KLE (draws a very large Gaussian vector prior based on Maatouk et al. (2024a)), samp.WC (draws a sample based on the Fast Fourier Transform (FFT) sampling scheme developed in Wood and Chan (1994)), ESS (draws a sample using elliptical slice sampler by Murray et al. (2010)).

One can find detailed description of each functions in the file.

2. 'all_models.R' file contains all models implementing efficient constrained Gaussian Processes for shape-restricted function estimation, such as:
   • linCGP.ESS — the proposed, highly efficient approach based on elliptical slice sampling,
   • linCGP.HMC — an approach based on the Hamiltonian Monte Carlo (HMC) sampler, and
   • linCGP.WC.ESS — an approach based on elliptical slice sampler and the Fast Fourier Transform.

These models can incorporate multiple shape constraints and estimate hyperparameters. The proposed \textsf{LS-ESS} function can utilize various efficient prior samplers, including Fast.LS, LS.KLE, and Cholesky. This function has been tested under extreme constraint scenarios where the HMC sampler fails to perform reliably. The file also includes the model 'Bound.CGP.general.R', which uses Botev's sampler to generate samples from the posterior truncated multivariate normal (MVN) distribution (see Botev, 2017, below).

3. 'Botev_boundedness.R' file contains code for boundedness constraints using the minimax tilting accept-reject method proposed by Botev (2017).

4. 'bound_extreme.R' file contains code for implementing boundedness constraints in an extreme case (see Maatouk et al. (2025b)). Through this numerical example, we investigate the performance and flexibility of the proposed approach, referred to as \textsf{LS-ESS}. It is worth noting that the HMC sampler fails to generate sample paths from the posterior distribution in this extreme case.

5. 'bound_syn.R' file contains code for implementing boundedness constraints using different posterior and prior samples. Through this numerical example, we compare the proposed approach, \textsff{LS-ESS}, with the HMC sampler and the FFT-WC in terms of prediction accurary and computational run-time.

6. 'mon_syn.R' file contains code for implementing monotonicity constraints using different posterior and prior samples. Through this numerical example, we compare the proposed approach, \textsff{LS-ESS}, with the HMC sampler and the FFT-WC in terms of prediction accurary and computational run-time. It is important to note that the proposed approach reduces the effect of the mass-shifting phenomenon well described in Zhou et al. (2024).

7. 'mon_bData.R' and 'mon_bData_N.R' files contain code for implementing monotonicity constraints with large datasets using the proposed \textsf{LS-ESS} approach. Through these two numerical examples, we investigate the performance of the proposed approach in terms of computational running time and prediction accuracy when the number of training samples n is large.

8. 'multi_monot_bd_syn.R' file contains code for implementing multiple shape-constraints (monotonicity and boundedness). Through this numerical example, we compare the proposed \textsf{LS-ESS} approach with the HMC sampler.

   For more details on the codes or the functions, refer to the associated R files.

Usage:

setwd("path to all_base_functions.R and all_models.R")
source('all_base_functions.R')
source('all_models.R')

## true monotone function 
f <- function(x) {
  3/(1 + exp(-10 * x + 2.1))
}
## synthetic training data
set.seed(12345)
n <- 100 # nb of training samples
sigN <- 0.5 # noise sd
xtr <- runif(n = n, min = 0, max = 1) # covariates
ytr <- sapply(X = xtr, FUN = f) + rnorm(n = n, mean = 0, sd = sigN)
N1 <- 9 # size of the 1st subdomain
M <- 3 # nb of subdomains
N <- N1 * M # nb of knot points
eta <- 50 # smooth approximate parameter
nsim <- 5000 # nb of mcmc retained samples
brn <- 1000 # burn-in
thin <- 1
nu <- 1.5 # smoothness parameter for Matérn cov function
l <- round(x = l_est(nu = nu, range = c(0, 1), val = 0.05), digits = 2) # length-scale parameter
tol <- 1e-12
## Illustration LS-ESS
post_mon_LS <- linCGP.ESS(y = ytr, x = xtr, N1 = N1, M = M, nu = nu, l = l, eta_fix = eta,
                          nsim = nsim, burn.in = brn, thin = thin,
                          sigsq.in = 1, tausq.in = 1, constrType = 'increasing',
                          prior = 'Fast.LS', return.plot = T, tol = tol, sseed = 12345)

mtext(text = 'LS-ESS with monotonicity', side = 3, line = 0.8, cex = 0.8)
mtext(text = paste("Running Time (s) =", round(x = post_mon_LS$time[3], digits = 2), "and dimension N =", N), side = 3, line = 0.1, cex = 0.8)
x <- seq(from = 0, to = 1, length = 100)
lines(x = x, y = f(x), type = 'l', lwd = 2)

## end

Illustrative examples:

Multiple constraints: GP approximation using the proposed LS-ESS approach under only monotonicity constraints (left) and under both monotonicity and boundedness constraints (middle). The HMC approach is employed in the right panel under both constraints. The black dashed horizontal line represents the upper bound constraint.

Monotonicity 7D: From left to right: the mean a posteriori (mAP) estimate based on 6,000 MCMC iterations, the Maximum a posteriori (MAP) estimate based on the posterior mode, and the target function. All these functions are evaluated at (x1,x2,0,...,0) (top) and at (x1,0,x3,0,...,0) (bottom), where (x1,x2,x3) in [0,1]^3.

Botev Boundedness: Bounded function estimation using the Botev's sampler from the \texttt{R} package \textsf{TruncatedNormal} with $N=27$ basis function (i.e. $N=27$ knots). The computation running time (in seconds) for generating 6,000 bounded sample paths.

Boundedness in an extreme case: Bounded function estimation in an extreme case using the proposed \textsf{LS-ESS} approach.

Monotonicity with large datasets: Monotone function estimation using the proposed \textsf{LS-ESS} approach with n = 1,000 (left) and n = 5,000 (right) training samples. The computation time (in seconds) for 6,000 MCMC iterations is shown in each panel, with the first 1,000 iterations discarded as burn-in.

Run-time large datasets: The runtime in seconds as a function of the number of samples (up to n=100,000) for generating 6,000 MCMC monotone and bounded surfaces using the proposed generalization, where the first 1,000 are discarded as burn-in.

Prediction accuracy for monotonicity with large datasets: The average mean squared prediction error (MSPE) over ten replicates as a function of the number of knots N for the two proposed estimators, MAP and mAP, in monotone function estimation. The number of training samples is fixed at n=1,000 (left) and n=5,000 (right).

Real world data application: age and income: Estimation accuracy of the two competing methods applied on the age-income data. The WAIC values corresponding to the two competing methods are displayed in the main of each panel based on 6,000 MCMC posterior sample paths.

Mass-shifting phenomenon (bivariate truncated normal distribution): Histogram of the first variable, from a zero-mean bivariate Gaussian vector, restricted to the positive orthant, for different values of the correlation parameter ρ. Through this image, we highlight the mass-shifting phenomenon observed when generating a truncated multivariate normal (MVN) distribution. This phenomenon manifests in the marginal densities of a dependent truncated MVN (tMVN) distribution. It is important to note that the effect becomes more pronounced as the correlation between components increases.

Approx Indicator fct: Smooth approximations of the indicator function on [−2,1] for different values of η.

Authors:

Hassan Maatouk (LAMPS, UPVD), Didier Rullière (Mines de St-Étienne), and Xavier Bay (Mines de St-Étienne)

Maintainer:

Hassan Maatouk, hassan.maatouk@univ-perp.fr

References:

Botev, Z. (2017). "The normal law under linear restrictions: simulation and estimation via minimax tilting". Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(1): 125–148

Cao, J. and Katzfuss, M. (2023). "Linear-cost Vecchia approximation of multivariate normal probabilities". arXiv preprint arXiv:2311.09426

Maatouk, H. and Bay, X. (2017). "Gaussian process emulators for computer experiments with inequality constraints". Mathematical Geosciences, 49(5): 557-582 doi

Maatouk, H. and Rullière, D., and Bay, X. (2024a). "Sampling large hyperplane‐truncated multivariate normal distributions". Computational Statistics, 39: 1779–1806 doi

Maatouk, H. and Rullière, D., and Bay, X. (2024b). "Bayesian analysis of constrained Gaussian processes". Bayesian Analysis, doi: 10.1214/24-BA1429 pdf Supplementary Material

Maatouk, H. and Rullière, D., and Bay, X. (2025a). "Large-scale constrained Gaussian processes for shape-restricted function estimation". Statistics and Computing, 35(7) doi

Maatouk, H. and Rullière, D. and Bay, X. (2025b). "Efficient constrained Gaussian processes using elliptical slice sampling" under review, preprint

Murray I., Adams R., and MacKay D. (2010). "Elliptical slice sampling". In: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, pp 541–548

Pakman, A. and Paninski, L. (2014). "Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians". Journal of Computational and Graphical Statistics, 23(2): 518–542

Ray, P., Pati, D., and Bhattacharya, A. (2020). "Efficient Bayesian shape-restricted function estimation with constrained Gaussian process priors". Statistics and Computing, 30(4): 839–853

Zhou, S., Ray, P., Pati, D., and Bhattacharya, A. (2024). "A mass-shifting phenomenon of truncated multivariate normal priors". Journal of the American Statistical Association, 119(545): 582–596