# Relaxed Gaussian process interpolation: a goal-oriented approach to Bayesian optimization

Sébastien J Petit, Julien Bect, Emmanuel Vazquez

► To cite this version:

Sébastien J Petit, Julien Bect, Emmanuel Vazquez. Relaxed Gaussian process interpolation: a goal-oriented approach to Bayesian optimization. 2022. hal-03686949v2

HAL Id: hal-03686949

https://hal.science/hal-03686949v2

Preprint submitted on 21 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - NonCommercial - NoDerivatives 4.0
International License
# Relaxed Gaussian process interpolation: a goal-oriented approach to Bayesian optimization

Sébastien Petit, Julien Bect, Emmanuel Vazquez

► To cite this version:

Sébastien Petit, Julien Bect, Emmanuel Vazquez. Relaxed Gaussian process interpolation: a goal-oriented approach to Bayesian optimization. Journal of Machine Learning Research, 2025, 26, pp.1-70. hal-03686949v4

HAL Id: hal-03686949

https://hal.science/hal-03686949v4

Submitted on 10 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Relaxed Gaussian Process Interpolation:
a Goal-Oriented Approach to Bayesian Optimization

Sébastien J. Petit *†

Julien Bect †

Emmanuel Vazquez †

SEBASTIEN.JULIEN.PETIT@GMAIL.COM

JULIEN.BECT@CENTRALESUPELEC.FR

EMMANUEL.VAZQUEZ@CENTRALESUPELEC.FR

Abstract

This work presents a new procedure for obtaining predictive distributions in the context of Gaus-
sian process (GP) modeling, with a relaxation of the interpolation constraints outside some ranges
of interest: the mean of the predictive distributions no longer necessarily interpolates the observed
values when they are outside ranges of interest, but are simply constrained to remain outside. This
method called relaxed Gaussian process (reGP) interpolation provides better predictive distribu-
tions in ranges of interest, especially in cases where a stationarity assumption for the GP model
is not appropriate. It can be viewed as a goal-oriented method and becomes particularly interest-
ing in Bayesian optimization, for example, for the minimization of an objective function, where
good predictive distributions for low function values are important. When the expected improve-
ment criterion and reGP are used for sequentially choosing evaluation points, the convergence of
the resulting optimization algorithm is theoretically guaranteed (provided that the function to be
optimized lies in the reproducing kernel Hilbert spaces attached to the known covariance of the un-
derlying Gaussian process). Experiments indicate that using reGP instead of stationary GP models
in Bayesian optimization is beneficial.

**Keywords:** Gaussian processes; Bayesian optimization; Expected improvement; Goal-oriented modeling; Reproducing kernel Hilbert spaces

1. Introduction

1.1 Context and Motivation

Gaussian process (GP) interpolation and regression (see, e.g., Stein, 1999; Rasmussen and Williams, 2006) is a very classical method for predicting an unknown function from data. It has found applications in active learning techniques, and notably in Bayesian optimization, a popular derivative-free global optimization technique for functions whose evaluations are time-consuming.

A GP model is defined by a mean and a covariance functions, which are generally selected from
data within parametric families. The most popular models assume stationarity and rely on standard
covariance functions such as the Matérn covariance. The assumption of stationarity yields models
with relatively low-dimensional parameters. However, such a hypothesis can sometimes result in
poor models when the function to be predicted has different scales of variation or different local
regularities across the domain.

<footnote>∗ . Safran Aircraft Engines, 77550, Moissy-Cramayel, France</footnote>

<footnote>†. Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France</footnote>
# Relaxed Gaussian Process Interpolation:
a Goal-Oriented Approach to Bayesian Optimization

Sébastien J. Petit

Laboratoire National de Métrologie et d'Essais,
78197, Trappes Cedex, France*

SEBASTIEN.PETIT@LNE.FR

Julien Bect

Université Paris-Saclay,
CNRS, CentraleSupélec,
Laboratoire des signaux et systèmes,
91190, Gif-sur-Yvette, France

JULIEN.BECT@CENTRALESUPELEC.FR

Emmanuel Vazquez

Université Paris-Saclay,
CNRS, CentraleSupélec,
Laboratoire des signaux et systèmes,
91190, Gif-sur-Yvette, France

EMMANUEL.VAZQUEZ@CENTRALESUPELEC.FR

**Editor:** Ryan Adams

## Abstract

This work presents a new procedure for obtaining predictive distributions in the context of Gaussian process (GP) modeling, with a relaxation of the interpolation constraints outside ranges of interest: the mean of the predictive distribution no longer necessarily interpolates the observed values when they are outside ranges of interest, but is simply constrained to remain outside. This method called relaxed Gaussian process (reGP) interpolation provides better predictive distributions in ranges of interest, especially in cases where a stationarity assumption for the GP model is not appropriate. It can be viewed as a goal-oriented method and becomes particularly interesting in Bayesian optimization, for example, for the minimization of an objective function, where good predictive distributions for low function values are important. When the expected improvement criterion and reGP are used for sequentially choosing evaluation points, the convergence of the resulting optimization algorithm is theoretically guaranteed (provided that the function to be optimized lies in the reproducing kernel Hilbert space attached to the known covariance of the underlying Gaussian process). Experiments indicate that using reGP instead of stationary GP models in Bayesian optimization is beneficial.

**Keywords:** Gaussian processes; Bayesian optimization; Expected improvement; Goal-oriented modeling; Reproducing kernel Hilbert spaces

<footnote>* Part of this work was funded by the French Agence Nationale de la Recherche et de la Technologie (ANRT) under a CIFRE grant, and was carried out while S. J. Petit was affiliated with Safran Aircraft Engines (Moissy-Cramayel, France) and L2S.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Left: the Steep function. Right: same illustration with a restrained range on the y-axis. The variations on the left overshadow the global minimum on the right.</figcaption></figure>

This is the case for instance in the motivating example given by Gramacy and Lee (2008), or in the even simpler toy minimization problem shown in Figure 1. The objective function in this example, which we shall call the Steep function, is smooth with an obvious global minimum around the point $x = 8$. However, the variations around the minimum are overshadowed by some steep variations on the left. Figure 2 shows a stationary GP fit with $n = 8$ points, where the parameters of the covariance function have been selected using maximum likelihood. Observe that the confidence bands are too large and that the conditional mean varies too much in the neighborhood of the global minimum, consistently with the stationary GP model that reflects the prior that our function oscillates around a mean value with a constant scale of variations. In this case, even if GP interpolation is consistent (Vazquez and Bect, 2010a), stationarity seems an unsatisfactory assumption for the Steep function. One expects Bayesian optimization techniques to be somehow inefficient on this problem with such a stationary model, whose posterior distributions are too pessimistic in the region of the minimum.

Nevertheless, the Steep function has the characteristics of an easy optimization problem: it has only two local minima, with the global minimum lying in a valley of significant volume. Consequently, a Bayesian optimization technique could be competitive if it relied on a model giving good predictions in regions where the function takes low values. In this work, we propose to explore goal-oriented GP modeling, where we want predictive models in regions of interest, even if it means being less predictive elsewhere.

## 1.2 Related Works

Going beyond the stationary hypothesis has been an active direction of research. With maybe a little bit of oversimplification, one can distinguish two categories of approaches that all use station-
# 1 Introduction

## 1.1 Context and Motivation

Gaussian process (GP) interpolation and regression (see, e.g., Stein, 1999; Rasmussen and Williams, 2006) is a classical method for predicting an unknown function from data. It has found applications in active learning techniques, and notably in Bayesian optimization, a popular derivative-free global optimization technique for functions whose evaluations are time-consuming.

A GP model is defined by a mean function and a covariance function, which are generally selected from data within parametric families. The most popular models assume stationarity and rely on standard covariance functions such as the Matérn covariance. The assumption of stationarity yields models with relatively low-dimensional parameters. However, such a hypothesis can sometimes result in poor models when the function to be predicted has different scales of variation or different local regularities across the domain.

This is the case for instance in the motivating example given by Gramacy and Lee (2008), or in the even simpler toy minimization problem shown in Figure 1. The objective function in this example, which we shall call the Steep function, is smooth with an obvious global minimum around the point $x = 8$. However, the variations around the minimum are overshadowed by some steep variations on the left. Figure 2 shows a stationary GP fit with $n = 8$ points, where the parameters of the covariance function have been selected using maximum likelihood. Observe that the confidence bands are too large and that the conditional mean varies too much in the neighborhood of the global minimum, consistently with the stationary GP model that reflects the prior that our function oscillates around a mean value with a constant scale of variations. In this case, even if GP interpolation is consistent (Vazquez and Bect, 2010a), stationarity seems an unsatisfactory assumption for the Steep function. One expects Bayesian optimization techniques to be somehow inefficient on this problem with such a stationary model, whose posterior distributions are too pessimistic in the region of the minimum.

Nevertheless, the Steep function has the characteristics of an easy optimization problem: it has only two local minima, with the global minimum lying in a valley of significant volume. Consequently, a Bayesian optimization technique could be competitive if it relied on a model giving good predictions in regions where the function takes low values. In this work, we propose to explore goal-oriented GP modeling, where we want predictive models in regions of interest, even if it means being less predictive elsewhere.

## 1.2 Related Works

### 1.2.1 LOCAL MODELS

Going beyond the stationary hypothesis has been an active direction of research. Local models are one popular solution, which still uses stationary Gaussian processes as a core building block.

A first class of local models is obtained by considering partitions of the input domain with different GP models on each subset. Partitions can be built by splitting the domain along the coordinate axes. This is the case of the treed Gaussian process models proposed by Gramacy and Lee (2008), which combines a fully Bayesian framework and the use of RJ-MCMC techniques for the inference, or the trust-region method by Eriksson et al. (2019). Park and Apley (2018) also propose partition-based local models built by splitting the domain along principal component directions.
<figure><img src="image_3.png" /><figcaption>Figure 2: Left: GP fit on the Steep function. Right: same illustration with a restrained range on the y-axis. The squares represent the data. The red line represents the posterior mean $\mu_n$ given by the model and the gray envelopes represent the associated uncertainties.</figcaption></figure>

ary Gaussian processes as a core building block: local models and transformation/composition of models.

### 1.2.1 LOCAL MODELS

A first class of local models is obtained by considering partitions of the input domain with different GP models on each subset. Partitions can be built by splitting the domain along the coordinate axes. This is the case of the treed Gaussian process models proposed by Gramacy and Lee (2008), which combines a fully Bayesian framework and the use of RJ-MCMC techniques for the inference, or the trust-region method by Eriksson et al. (2019). Park and Apley (2018) also propose partition-based local models built by splitting the domain along principal component directions. In such techniques, there are parameters (often many of them) related to, e.g., the way the partitions evolve with the data, the size of the partitions, or how local Gaussian processes interact with each other.

A second class of local models is obtained by spatially weighting one or several GP models. Many schemes have been proposed, including methods based on partition of unity (Nott and Dunsmuir, 2002), weightings of covariance functions (Pronzato and Rendas, 2017; Rivoirard and Romary, 2011), and convolution techniques (see, e.g., Higdon, 1998; Gibbs, 1998; Higdon, 2002; Ver Hoef et al., 2004; Stein, 2005). Let us also mention data-driven aggregation techniques: composite Gaussian process models (Ba and Joseph, 2012), and mixture of experts techniques (see, e.g., Tresp, 2001; Rasmussen and Ghahramani, 2002; Meeds and Osindero, 2006; Yuan and Neubauer, 2009; Yang and Ma, 2011; Yuksel et al., 2012). In the latter framework, the weights are called gating functions and the estimation of the parameters and the inference are usually performed using EM, MCMC, or variational techniques. Weighting methods generally have parameters specifying weighting functions, with an increased need to watch for overfitting phenomena.
<figure><img src="image_2.png" /><figcaption>Figure 1: Left: the Steep function. Right: same illustration with a restrained range on the y-axis. The variations on the left tend to overshadow the global minimum on the right.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Left: GP fit on the Steep function. Right: same illustration with a restrained range on the y-axis. The squares represent the data. The red line represents the posterior mean $\mu_n$ given by the model and the gray envelopes represent the associated uncertainties.</figcaption></figure>
1.2.2 TRANSFORMATION AND COMPOSITION OF MODELS

A first technique for composition of models consists in using a parametric transformation of a GP (Rychlik et al., 1997; Snelson et al., 2004).

Another route is to transform the input domain, using for instance a parametric density (Xiong et al., 2007), or other parametric transformations involving possible dimension reduction (Marmin et al., 2018). Bodin et al. (2020) proposed a framework that uses additional input variables, serving as nuisance parameters, to smooth out some badly behaved data. The practitioner has to specify a prior over the variance of the nuisance parameter and inference is based on MCMC.

Lázaro-Gredilla (2012) takes the step of choosing a GP prior on the output transform and resorts to variational inference techniques for inference. This type of idea can be viewed as an ancestor of deep Gaussian processes (see, e.g., Damianou and Lawrence, 2013; Dunlop et al., 2018; Hebbal et al., 2021; Jakkala, 2021; Bachoc and Lagnoux, 2021), which stack layers of linear combinations of GPs. The practitioner has to specify a network structure among other parameters and resort to variational inference.

Recently, Picheny et al. (2019) proposed another approach where prediction is made only from pairwise comparisons between data points, relying on the variational framework of ordinal GP regression proposed by Chu and Ghahramani (2005) for the inference.

**1.3 Contributions and Outline**

The brief review of the literature above reveals three types of shortcomings in methods that depart
from the stationarity hypothesis: 1) they rely on advanced techniques for deriving predictive distri-
butions; or 2) they require the practitioner to choose in advance some key parameters; or 3) they
increase the number of parameters with an increased risk of overfitting.

This article suggests a method for building models targeting regions of interest specified through function values. The main objective is to obtain global models that exhibit good predictive distributions on a range of interest. In the case of a minimization problem, the range of interest would be the values below a threshold. Outside the range of interest, we accept that the model can be less predictive by relaxing the interpolation constraints. Such a model is presented in Figure 3: compared to the situation in Figure 2, the model is more predictive in the region where the Steep function takes low values, with expected benefits for the efficiency of Bayesian optimization.

This article provides three main contributions. First, we propose a class of goal-oriented GP-based models called *relaxed Gaussian processes* (reGP). Second, we give theoretical and empirical results justifying the method and its use for Bayesian optimization. Finally, to assess the predictivity of reGP, we adopt the formalism of scoring rules (Gneiting and Raftery, 2007) and propose the use of a goal-oriented scoring rule that we call *truncated continuous ranked probability score* (tCRPS), which is designed to assess the predictivity of a model in a range of interest.

The organization is as follows. Section 2 briefly recalls the formalism of Gaussian processes and Bayesian optimization (BO). Section 3 presents reGP and its theoretical properties. The tCRPS and its use for selecting regions of interest are then presented in Section 4. Section 5 presents a reGP-based Bayesian optimization algorithm called EGO-R, together with the convergence analysis
In such techniques, there are parameters related to, e.g., the way the partitions evolve with the data, the size of the partitions, or how local Gaussian processes interact with each other.

A second class of local models is obtained by spatially weighting one or several GP models. Many schemes have been proposed, including methods based on partition of unity (Nott and Dunsmuir, 2002), weightings of covariance functions (Pronzato and Rendas, 2017; Rivoirard and Romary, 2011), and convolution techniques (see, e.g., Higdon, 1998; Gibbs, 1998; Higdon, 2002; Ver Hoef et al., 2004; Stein, 2005). Let us also mention data-driven aggregation techniques: composite Gaussian process models (Ba and Joseph, 2012), and mixture of experts techniques (see, e.g., Tresp, 2001; Rasmussen and Ghahramani, 2002; Meeds and Osindero, 2006; Yuan and Neubauer, 2009; Yang and Ma, 2011; Yuksel et al., 2012). In the latter framework, the weights are called gating functions and the estimation of the parameters and the inference are usually performed using EM, MCMC, or variational techniques. Weighting methods generally have parameters specifying weighting functions, with an increased need to watch for overfitting phenomena.

### 1.2.2 TRANSFORMATION AND COMPOSITION OF MODELS

Transformation and composition of models is another popular alternative to stationary GP modeling, which also uses stationary Gaussian processes as a core building block.

A first technique for composition of models consists in using a parametric transformation of a GP (Rychlik et al., 1997; Snelson et al., 2004).

Another route is to transform the input domain, using for instance a parametric density (Xiong et al., 2007), or other parametric transformations involving possible dimension reduction (Marmin et al., 2018). Bodin et al. (2020) proposed a framework that uses additional input variables, serving as nuisance parameters, to smooth out some badly behaved data. The practitioner has to specify a prior over the variance of the nuisance parameter and inference is based on MCMC.

Lázaro-Gredilla (2012) takes the step of choosing a GP prior on the output transform and sorts to variational inference techniques. This type of idea can be viewed as an ancestor of deep Gaussian processes (see, e.g., Damianou and Lawrence, 2013; Dunlop et al., 2018; Hebbal et al., 2021; Jakkala, 2021; Bachoc and Lagnoux, 2025), which stack layers of linear combinations of GPs. The practitioner has to specify a network structure among other parameters and resort to variational inference.

Picheny et al. (2019) proposed another approach where predictions are made only from pairwise comparisons between data points, relying on the variational framework of ordinal GP regression proposed by Chu and Ghahramani (2005) for the inference.

### 1.2.3 GOAL-ORIENTED APPROXIMATE INFERENCE TECHNIQUES

Another related area of research is goal-oriented approximate inference, where the task to be performed with the predictive model is taken into account in the approximations that are used in the inference (a.k.a. training) step.

For instance, in the context of scalable Bayesian optimization, McIntire et al. (2016) and Moss et al. (2023) have proposed goal-oriented criteria to choose inducing points. Another example is Yang et al. (2021), who propose a goal-oriented criterion to select the frequencies of a sparse spectrum approximation. Loss-calibrated approximate Bayesian computation (Lacoste-Julien et al., 2011; Morais and Pillow, 2022), which takes into account the decision-making task when approximating the posterior distribution, also falls in this category.
<figure><img src="image_4.png" /><figcaption>Figure 3: Left: prediction of the Steep function with the proposed methodology (black line: relaxation threshold $t$; blue points: relaxed observations). Right: $\mu_n$ versus $f$ (with more observations for illustration purposes). The model interpolates the data below $t$. The blue points are relaxed observations.</figcaption></figure>

of this algorithm and a numerical benchmark. Finally, Section 6 presents our conclusions and perspectives for future work.

# 2. Background and Notations

## 2.1 Gaussian Process Modeling

Consider a real-valued function $f: \mathbb{X} \to \mathbb{R}$, where $\mathbb{X} \subseteq \mathbb{R}^d$, and suppose we want to infer $f$ at a given $x \in \mathbb{X}$ from evaluations of $f$ on a finite set of points $\underline{x}_n = (x_1, \dots, x_n) \in \mathbb{X}^n$, $n \ge 1$. A standard Bayesian approach to this problem consists in using a GP model $\xi \sim \text{GP}(\mu, k)$ as a prior about $f$, where $\mu: \mathbb{X} \to \mathbb{R}$ is a mean function and $k: \mathbb{X} \times \mathbb{X} \to \mathbb{R}$ is a covariance function, which is supposed to be strictly positive-definite in this article.

The posterior distribution of $\xi$ given $\underline{Z}_n = (\xi(x_1), \dots, \xi(x_n))^\mathrm{T}$ is still a Gaussian process, whose mean and covariance functions are given by the standard kriging equations (Matheron, 1971). More precisely:

$$ \xi | \underline{Z}_n \sim \text{GP}(\mu_n, k_n), \qquad (1) $$

with

$$ \mu_n(x) = \mu(x) + k(x, \underline{x}_n) K_n^{-1} (\underline{Z}_n - \mu(\underline{x}_n)) \qquad (2) $$

and

$$ k_n(x, y) = k(x, y) - k(x, \underline{x}_n) K_n^{-1} k(y, \underline{x}_n)^\mathrm{T}, \qquad (3) $$

and where $\mu(\underline{x}_n) = (\mu(x_1), \dots, \mu(x_n))^\mathrm{T}$, $k(x, \underline{x}_n) = (k(x, x_1), \dots, k(x, x_n))$, and $K_n$ is the $n \times n$ matrix with entries $k(x_i, x_j)$. We shall also use the notation $\sigma_n^2(x) = k_n(x, x)$ for the posterior variance, a.k.a. the kriging variance, a.k.a. the squared power function, so that $\xi(x) | \underline{Z}_n \sim \mathcal{N}(\mu_n(x), \sigma_n^2(x))$.
These works focus on fine-tuning approximations of classical GP models, by taking into account
the task to be performed. This is in contrast with the method introduced in this article, which aims
to *replace* the classical GP model with a new, goal-oriented GP-based predictive distribution, which
no longer stems from conditioning a GP prior. (The development of a massively scalable version
of the proposed method is beyond the scope of this article, but could be based on approximate
goal-directed inference techniques.)

**1.3 Contributions and Outline**

The main contribution of this article is a method called *relaxed Gaussian processes* (reGP) for
building goal-oriented GP-based models targeting regions of interest specified through function
values. The objective is to obtain global models that exhibit good predictive distributions on a range
of interest—in the case of a minimization problem, the range of interest would be the values below
a threshold—while possibly being less predictive outside the range of interest. This is achieved
by relaxing interpolation constraints outside this range. Such a model is presented in Figure 3:
compared to the situation in Figure 2, the model is more predictive in the region where the Steep
function takes low values, with expected benefits for the efficiency of Bayesian optimization.

<figure><img src="image_4.png" /><figcaption>Figure 3: Left: prediction of the Steep function with the proposed methodology (black line: relaxation threshold *t*; blue points: relaxed observations). Right: *μ*<sub>*n*</sub> versus *f* (with more observations for illustration purposes). The model interpolates the data below *t*. The blue points are relaxed observations.</figcaption></figure>

This article provides two other main contributions. On the one hand, we give theoretical and
empirical results justifying the method and its use for Bayesian optimization. On the other hand,
to assess the predictivity of reGP models, we adopt the formalism of scoring rules (Gneiting and
Raftery, 2007) and propose the use of a goal-oriented scoring rule that we call *truncated continuous
ranked probability score* (tCRPS), which is designed to assess the predictivity of a model in a range
of interest. The tCRPS is used to construct a leave-one-out (LOO) goodness-of-fit criterion to auto-
The functions $\mu$ and $k$ control the posterior distribution (1) and must be chosen carefully. The standard practice is to select them from data within a parametric family $\{ (\mu_{\theta}, k_{\theta}), \theta \in \Theta \}$. A common approach is to suppose stationarity for the GP, which means choosing a constant mean function $\mu \equiv c \in \mathbb{R}$ and a stationary covariance function $k(x, y) = \sigma^2 r(x-y)$, where $r: \mathbb{R}^d \to \mathbb{R}$ is a stationary correlation function.

A correlation function often recommended in the literature (Stein, 1999) is the (geometrically anisotropic) Matérn correlation function

$$r(h) = \frac{2^{1-\nu}}{\Gamma(\nu)} (\sqrt{2\nu}\|h\|_{\rho})^{\nu} \mathcal{K}_{\nu} (\sqrt{2\nu}\|h\|_{\rho}), \quad \|h\|_{\rho}^2 = \sum_{j=1}^{d} \frac{h_{[j]}^2}{\rho_j^2}, \qquad (4)$$

for $h = (h_{[1]}, \dots, h_{[d]}) \in \mathbb{R}^d$, and where $\Gamma$ is the Gamma function and $\mathcal{K}_{\nu}$ is the modified Bessel function of the second kind. The parameters to be selected in this case are $(\sigma^2, \rho_1, \dots, \rho_d, \nu) \in (0, \infty)^{d+2}$ with $\sigma^2$ the process variance, $\rho_i$ the range parameter along the $i$-th dimension, and $\nu$ a regularity parameter controlling the smoothness of the process. Two other standard covariance functions can be recovered for specific values of $\nu$: the exponential covariance function for $\nu = 1/2$ and the squared-exponential covariance function for $\nu \to \infty$.

A variety of techniques for selecting the parameter $\theta$ have been proposed in the literature, but we can safely say that maximum likelihood estimation is the most popular and can be recommended in the case of interpolation (Petit et al., 2021). It simply consists in minimizing the negative log-likelihood

$$\mathcal{L}(\theta; \underline{Z}_n) = -\log(p(\underline{Z}_n | \theta)) \propto \log(\det(K_n)) + (\underline{Z}_n - \mu(\underline{x}_n))^\mathrm{T} K_n^{-1} (\underline{Z}_n - \mu(\underline{x}_n)) + \text{constant}, \quad (5)$$

where $p$ stands for the probability density of $\underline{Z}_n$. Other methods for selecting the parameters include the restricted maximum likelihood method and leave-one-out strategies (see, e.g., Stein, 1999; Rasmussen and Williams, 2006).

## 2.2 Bayesian Optimization

The framework of GPs is well suited to the problem of sequential design of experiments, or active learning. In particular, for minimizing a real-valued function $f$ defined on a compact domain $X$, the Bayesian approach consists in choosing sequentially evaluation points $X_1, X_2, \dots \in X$ using a GP model $\xi$ for $f$, which makes it to possible to build a sampling criterion that represents an expected information gain on the minimum of $f$ when an evaluation is made at a new point. One of the most popular sampling criterion (also called acquisition function) is the Expected Improvement (EI) (Mockus et al., 1978; Jones et al., 1998), which can be expressed as

$$\rho_n(x) = \mathrm{E}((m_n - \xi(x))_+ | \underline{Z}_n), \qquad (6)$$

where $m_n = \min(\xi(x_1), \dots, \xi(x_n))$. The EI criterion corresponds to the expectation of the excursion of $\xi$ below the minimum given $n$ observations, and can be written in closed form:
matically select the range of function values (outside the range of interest) where the interpolation constraints are relaxed.

An important feature of the proposed approach is that, with respect to classical GP modeling, it only requires the choice of one additional key parameter—namely, the range of interest where good predictions are expected. We propose several heuristics to set this range of interest in Bayesian optimization and level-set estimation problems.

The organization of this article is as follows. Section 2 briefly recalls the formalism of Gaussian processes and Bayesian optimization. Section 3 presents reGP and its theoretical properties. The tCRPS and its use for selecting the relaxation range are then presented in Section 4. Section 5 presents a reGP-based Bayesian optimization algorithm called EGO-R, together with the convergence analysis of this algorithm and a numerical benchmark. Finally, Section 6 presents our conclusions and perspectives for future work.

An open source implementation of the reGP method and the numerical experiments is available online at https://github.com/relaxedGP/regp_paper_experiments.

# 2 Background and Notations

## 2.1 Gaussian Process Modeling

Consider a real-valued function $f: \mathcal{X} \to \mathbb{R}$, where $\mathcal{X} \subseteq \mathbb{R}^d$, and suppose we want to infer $f$ at a given $x \in \mathcal{X}$ from evaluations of $f$ on a finite set of points $\underline{x}_n = (x_1, \dots, x_n) \in \mathcal{X}^n$, $n \ge 1$. A standard Bayesian approach to this problem consists in using a GP model $\xi \sim \text{GP}(\mu, k)$ as a prior about $f$, where $\mu: \mathcal{X} \to \mathbb{R}$ is a mean function and $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is a covariance function, which is supposed to be strictly positive-definite in this article.

The posterior distribution of $\xi$ given $\underline{Z}_n = (\xi(x_1), \dots, \xi(x_n))^\top$ is still a Gaussian process, whose mean and covariance functions are given by the standard kriging equations (Matheron, 1971). More precisely:

$$ \xi | \underline{Z}_n \sim \text{GP}(\mu_n, k_n), \qquad (1) $$

with

$$ \mu_n(x) = \mu(x) + k(x, \underline{x}_n) K_n^{-1} (\underline{Z}_n - \mu(\underline{x}_n)) \qquad (2) $$

and

$$ k_n(x, y) = k(x, y) - k(x, \underline{x}_n) K_n^{-1} k(y, \underline{x}_n)^\top, $$

and where $\mu(\underline{x}_n) = (\mu(x_1), \dots, \mu(x_n))^\top$, $k(x, \underline{x}_n) = (k(x, x_1), \dots, k(x, x_n))$, and $K_n$ is the $n \times n$ matrix with entries $k(x_i, x_j)$. We shall also use the notation $\sigma_n^2(x) = k_n(x, x)$ for the posterior variance, a.k.a. the kriging variance, a.k.a. the squared power function, so that $\xi(x) | \underline{Z}_n \sim \mathcal{N}(\mu_n(x), \sigma_n^2(x))$.

The functions $\mu$ and $k$ control the posterior distribution (1) and must be chosen carefully. The standard practice is to select them from data within a parametric family $\{ (\mu_\theta, k_\theta), \theta \in \Theta \}$. A common approach is to suppose stationarity for the GP, which means choosing a constant mean function $\mu \equiv c \in \mathbb{R}$ and a stationary covariance function $k(x, y) = \tau^2 r(x-y)$, where $r: \mathbb{R}^d \to \mathbb{R}$ is a stationary correlation function.

A correlation function often recommended in the literature (Stein, 1999) is the (geometrically anisotropic) Matérn correlation function

$$ r(h) = \frac{2^{1-\nu}}{\Gamma(\nu)} (\sqrt{2\nu} \|h\|_{\rho})^{\nu} \mathcal{H}_{\nu} (\sqrt{2\nu} \|h\|_{\rho}), \quad \|h\|_{\rho}^{2} = \sum_{j=1}^{d} \frac{h_{[j]}^{2}}{\rho_{j}^{2}}, \qquad (3) $$
**Proposition 1** (Jones et al., 1998; Vazquez and Bect, 2010b) The EI criterion may be written as
$$ \rho_n(x) = \gamma(m_n - \mu_n(x), \sigma_n^2(x)) \text{ with} $$

$$ \gamma: (z, s) \in \mathbb{R} \times \mathbb{R}_+ \mapsto \begin{cases} \sqrt{s}\phi\left(\frac{z}{\sqrt{s}}\right) + z\Phi\left(\frac{z}{\sqrt{s}}\right) & \text{if } s > 0, \\ \max(z, 0) & \text{if } s = 0, \end{cases} $$

where $\phi$ and $\Phi$ stand for the probability density and cumulative distribution functions of the standard Gaussian distribution. Moreover, the function $\gamma$ is continuous, verifies $\gamma(z, s) > 0$ if $s > 0$ and is non-decreasing with respect to $z$ and $s$ on $\mathbb{R} \times \mathbb{R}_+$.

When the EI criterion is used for optimization, that is, when the sequence of evaluation points $(X_n)_{n>0}$ of $f$ is chosen using the rule

$$ X_{n+1} = \arg\max_{x \in \mathcal{X}} \rho_n(x), $$

the resulting algorithm is generally called the Efficient Global Optimization (EGO) algorithm, as proposed by Jones et al. (1998). The EGO algorithm has known convergence properties (Vazquez and Bect, 2010b; Bull, 2011).

A variety of other sampling criteria for the minimization problem can be found in the literature (see, e.g., Frazier et al., 2008; Villemonteix et al., 2009; Srinivas et al., 2010; Vazquez and Bect, 2014), but we shall focus on the EI algorithm in this article.

## 2.3 Reproducing Kernel Hilbert Spaces

Reproducing kernel Hilbert spaces (RKHS, see e.g., Aronszajn, 1950; Berlinet and Thomas-Agnan, 2004) are Hilbert spaces of functions commonly used in the field of approximation theory (see, e.g., Wahba, 1990; Wendland, 2004). A Hilbert space $\mathcal{H}(X)$ of real-valued functions on $X$ with an inner product $(\cdot, \cdot)_{\mathcal{H}(X)}$ is called an RKHS if it has a reproducing kernel, that is, a function $k: X \times X \to \mathbb{R}$ such that $k(x, \cdot) \in \mathcal{H}(X)$, and

$$ (f, k(x, \cdot))_{\mathcal{H}(X)} = f(x) \qquad (7) $$

(the reproducing property), for all $x \in X$ and $f \in \mathcal{H}(X)$. Furthermore, given a (strictly) positive definite covariance function $k$, there exists a unique RKHS admitting $k$ as reproducing kernel.

Given locations $x_n = (x_1, \dots, x_n) \in X^n$, and corresponding values $z_n \in \mathbb{R}^n$, suppose we want to find a function $h \in \mathcal{H}(X)$ such that $h(x_n) = (h(x_1), \dots, h(x_n))^\top = z_n$. Then, the minimum-norm interpolation solution is given by the following proposition.

**Proposition 2** (Parzen, 1959; Matheron, 1981) The problem

$$ \left\{ \begin{array}{ll} \text{minimize} & \|h\|_{\mathcal{H}(X)} \\ \text{subject to} & h \in \mathcal{H}(X) \\ & h(x_n) = z_n \end{array} \right. \qquad (8) $$

has a unique solution given by $s_{z_n} = k(\cdot, x_n)K_n^{-1}z_n$.
for $h = (h_{[1]}, \dots, h_{[d]}) \in \mathbb{R}^d$, and where $\Gamma$ is the Gamma function and $\mathcal{K}_\nu$ is the modified Bessel function of the second kind. The covariance parameters to be selected in this case are $(\tau^2, \rho_1, \dots, \rho_d, \nu) \in (0, \infty)^{d+2}$ with $\tau^2$ the process variance, $\rho_i$ the range parameter along the $i$-th dimension, and $\nu$ a regularity parameter controlling the smoothness of the process. Two other standard covariance functions can be recovered for specific values of $\nu$: the exponential covariance function for $\nu = 1/2$ and the squared-exponential covariance function for $\nu \to \infty$.

A variety of techniques for selecting the parameter $\theta$ have been proposed in the literature, but we can safely say that maximum likelihood estimation is the most popular and can be recommended in the case of interpolation (Petit et al., 2023). It simply consists in minimizing the negative log-likelihood

$$ \mathcal{L}(\theta; \underline{Z}_n) = -\log(p(\underline{Z}_n | \theta)) \propto \log(\det(K_n)) + (\underline{Z}_n - \mu(\underline{x}_n))^\mathrm{T} K_n^{-1} (\underline{Z}_n - \mu(\underline{x}_n)) + \text{constant}, \quad (4) $$

where $p$ stands for the probability density of $\underline{Z}_n$. Other methods for selecting the parameters include the restricted maximum likelihood method and leave-one-out strategies (see, e.g., Stein, 1999; Rasmussen and Williams, 2006).

## 2.2 Bayesian Optimization

The framework of GPs is well suited to the problem of sequential design of experiments, or active learning. In particular, for minimizing a real-valued function $f$ defined on a compact domain $\mathcal{X}$, the Bayesian approach consists in sequentially choosing evaluation points $x_1, x_2, \dots \in \mathcal{X}$ using a GP model $\xi$ for $f$, which makes it to possible to build a sampling criterion that represents an expected information gain on the minimum of $f$ when an evaluation is made at a new point. One of the most popular sampling criteria (also called acquisition function) is the *Expected Improvement* (EI) (Mockus et al., 1978; Jones et al., 1998), which can be expressed as

$$ \rho_n(x) = \mathbb{E}((m_n - \xi(x))_+ | \underline{Z}_n), \quad (5) $$

where $m_n = \min(\xi(x_1), \dots, \xi(x_n))$. The EI criterion corresponds to the expectation of the excursion of $\xi$ below the minimum given $n$ observations, and can be written in closed form:

**Proposition 1** (Jones et al., 1998; Vazquez and Bect, 2010b) The EI criterion may be written as
$$ \rho_n(x) = \gamma(m_n - \mu_n(x), \sigma_n^2(x)) \text{ with} $$

$$ \gamma: (z, s) \in \mathbb{R} \times \mathbb{R}_+ \mapsto \begin{cases} \sqrt{s}\phi\left(\frac{z}{\sqrt{s}}\right) + z\Phi\left(\frac{z}{\sqrt{s}}\right) & \text{if } s > 0, \\ \max(z, 0) & \text{if } s = 0, \end{cases} $$

where $\phi$ and $\Phi$ stand for the probability density and cumulative distribution functions of the standard Gaussian distribution. Moreover, the function $\gamma$ is continuous, satisfies $\gamma(z, s) > 0$ if $s > 0$ and is non-decreasing with respect to $z$ and $s$ on $\mathbb{R} \times \mathbb{R}_+$.

When the EI criterion is used for optimization, that is, when the sequence of evaluation points $(x_n)_{n>0}$ of $f$ is chosen using the rule

$$ x_{n+1} = \arg\max_{x \in \mathcal{X}} \rho_n(x), $$
Observe that the solution $s_{\underline{z}_n}$ is equal to the posterior mean (2) when $\mu = 0$.

Moreover, for any $f \in \mathcal{H}(\mathcal{X})$ and $x \in \mathcal{X}$, the reproducing property (7) yields the upper bound

$$|f(x) - s_{\underline{z}_n}(x)| \le \sigma_n(x) \|f\|_{\mathcal{H}(\mathcal{X})}, \quad (9)$$

with $\sigma_n(x) = \sqrt{k_n(x, x)}$. Note that $\sigma_n(x)$ is the worst-case error at $x$ for the interpolation of functions in the unit ball of $\mathcal{H}(\mathcal{X})$.

# 3. Relaxed Gaussian Process Interpolation

## 3.1 Relaxed Interpolation

The example in the introduction (see Figures 1–3) suggests that, in order to gain accuracy over a range of values of interest, it can be beneficial to relax interpolation constraints outside this range. More precisely, the probabilistic model in Figure 3 interpolates data lying below a selected threshold $t$, and when data are above $t$, the model only keeps the information that the data exceeds $t$.

In the following, we consider the general setting where relaxation is carried out on a set of the form $R = \cup_{j=1}^J R_j$, where $R_1, \dots, R_J \subset \mathbb{R}$ are disjoint closed intervals with non-zero lengths. (The set $R = [t, +\infty)$ was used in the example of Figure 3).

As above, we shall write $x_n = (x_1, \dots, x_n) \in \mathcal{X}^n$ for a sequence of locations with corresponding function values $\underline{z}_n = (z_1, \dots, z_n)^\top \in \mathbb{R}^n$. Then, we introduce the set $C_{R,n} = C_1 \times \dots \times C_n \subset \mathbb{R}^n$ of relaxed constraints, where

$$\left\{ \begin{array}{ll} C_i = R_j & \text{if } z_i \in R_j \text{ for some } j, \\ C_i = \{z_i\} & \text{otherwise.} \end{array} \right. \qquad (10)$$

Let also $\mathcal{H}_{R,n} = \{h \in \mathcal{H}(\mathcal{X}) \mid h(x_n) \in C_{R,n}\}$ be the set of relaxed-interpolating functions. The following proposition gives the definition of the minimum-norm relaxed predictor.

**Proposition 3** *The problem*

$$\left\{ \begin{array}{ll} \text{minimize} & \|h\|_{\mathcal{H}(\mathcal{X})} \\ \text{subject to} & h \in \mathcal{H}_{R,n} \end{array} \right. \qquad (11)$$

has a unique solution given by $s_{\underline{z}_n}^*$, where $\underline{z}_n^*$ is the unique solution of the quadratic problem

$$\arg\min_{\underline{z} \in C_{R,n}} \underline{z}^\top K_n^{-1} \underline{z}. \qquad (12)$$

## 3.2 Relaxed Gaussian Process Interpolation

The main advantage of Gaussian processes is the possibility to obtain not only point predictions but also predictive distributions. However, Proposition 3 only defines a function approximation. We now turn relaxed interpolation into a probabilistic model providing predictive distributions whose mean is not constrained to interpolate data on a given range $R$. The following proposition makes a step in this direction.
the resulting algorithm is generally called the Efficient Global Optimization (EGO) algorithm, as proposed by Jones et al. (1998). The EGO algorithm has known convergence properties (Vazquez and Bect, 2010b; Bull, 2011).

A variety of other sampling criteria for the minimization problem can be found in the literature (see, e.g., Frazier et al., 2008; Villemonteix et al., 2009; Srinivas et al., 2010; Vazquez and Bect, 2014), but we shall focus on the EI criterion in this article. (Appendix E also considers the Upper Confidence Bound (UCB) criterion by Srinivas et al., 2010, and Appendix F deals with the task of estimating excursion sets using the straddle heuristic by Bryan et al., 2005.)

## 2.3 Reproducing Kernel Hilbert Spaces

Reproducing kernel Hilbert spaces (RKHS, see, e.g., Aronszajn, 1950; Berlinet and Thomas-Agnan, 2004) are Hilbert spaces of functions commonly used in the field of approximation theory (see, e.g., Wahba, 1990; Wendland, 2004). A Hilbert space $\mathcal{H}(X)$ of real-valued functions on $X$ with an inner product $(\cdot, \cdot)_{\mathcal{H}(X)}$ is called an RKHS if it has a reproducing kernel, that is, a function $k: X \times X \to \mathbb{R}$ such that $k(x, \cdot) \in \mathcal{H}(X)$, and

$$ (f, k(x, \cdot))_{\mathcal{H}(X)} = f(x) \qquad (6) $$

(the reproducing property), for all $x \in X$ and $f \in \mathcal{H}(X)$. Furthermore, given a (strictly) positive definite covariance function $k$, there exists a unique RKHS admitting $k$ as reproducing kernel.

Given locations $\underline{x}_n = (x_1, \dots, x_n) \in X^n$, and corresponding values $\underline{z}_n \in \mathbb{R}^n$, suppose we want to find a function $h \in \mathcal{H}(X)$ with minimal norm, such that $h(\underline{x}_n) = (h(x_1), \dots, h(x_n))^\mathrm{T} = \underline{z}_n$. Then the solution is given by the following classical result, which can be derived from the generalized representer theorem (Schölkopf et al., 2001) and appears in early work on optimal interpolation in RKHS (Kimeldorf and Wahba, 1970).

**Proposition 2 (Minimum-norm interpolant) The problem**

$$ \left\{ \begin{array}{ll} \text{minimize} & \|h\|_{\mathcal{H}(X)} \\ \text{subject to} & h \in \mathcal{H}(X) \\ & h(\underline{x}_n) = \underline{z}_n \end{array} \right. \qquad (7) $$

has a unique solution given by $s_{\underline{z}_n} = k(\cdot, \underline{x}_n)K_n^{-1}\underline{z}_n$.

Observe that the solution $s_{\underline{z}_n}$ is equal to the posterior mean (2) when $\mu = 0$.

Moreover, for any $f \in \mathcal{H}(X)$ and $x \in X$, the reproducing property (6) yields the upper bound

$$ |f(x) - s_{\underline{z}_n}(x)| \le \sigma_n(x) \|f\|_{\mathcal{H}(X)}, \qquad (8) $$

with $\sigma_n(x) = \sqrt{k_n(x, x)}$ and $\underline{z}_n = f(\underline{x}_n)$ (see, e.g., Scheurer et al., 2013, p. 13). Note that $\sigma_n(x)$ is the worst-case error at $x$ for the interpolation of functions in the unit ball of $\mathcal{H}(X)$.

# 3 Relaxed Gaussian Process Interpolation

## 3.1 Relaxed Interpolation

The example in the introduction (see Figures 1–3) suggests that, in order to gain accuracy over a range of values of interest, it can be beneficial to relax interpolation constraints outside this range.
**Proposition 4** Let $\xi \sim \text{GP}(0, k)$, $x_n = (x_1, \dots, x_n) \in \mathbb{X}^n$, $z_n \in \mathbb{R}^n$ and $x'_m = (x'_1, \dots, x'_m) \in \mathbb{X}^m$ be a set of locations of interest where predictions should be made. Write $\underline{Z}_n = (\xi(x_1), \dots, \xi(x_n))^{\top}$ and $\underline{Z}'_m = (\xi(x'_1), \dots, \xi(x'_m))^{\top}$. Then the mode of the probability density function

$$p(\underline{Z}'_m, \underline{Z}_n | \underline{Z}_n \in C_{R,n}) \qquad (13)$$

is given by $(s_{\underline{Z}_n}^*(x'_m), \underline{Z}_n^*)$.

In other words, the relaxed interpolation solution of Proposition 3 corresponds to the maximum a posteriori (MAP) estimate under the predictive model (13). Conditional distributions with respect to events of the type $\underline{Z}_n \in C_{R,n}$ have been used in Bayesian statistics for dealing with outliers and model misspecifications (see, e.g., Lewis et al., 2021, and references therein). This type of conditional distributions is also encountered for constrained GPs (see, e.g., Da Veiga and Marrel, 2012; Maatouk and Bay, 2017; López-Lopera et al., 2018), when constraints come from expert knowledge.

However, the predictive distribution (13) is non-Gaussian since the support of $\underline{Z}_n$ is truncated. In particular, no closed-form expression is available for any of its moments, and sampling requires advanced techniques (e.g., variational, MCMC). Motivated by this observation, we propose instead to build a goal-oriented probabilistic model using the following definition.

**Definition 5 (Relaxed-GP predictive distribution; fixed $\mu$ and $k$)** Given $x_n \in \mathbb{X}^n$, $z_n \in \mathbb{R}^n$, and a relaxation set $R$ (finite union of closed intervals), the relaxed-GP (reGP) predictive distribution with fixed mean function $\mu$ and covariance function $k$ is defined as the (Gaussian) conditional distribution of $\xi \sim \text{GP}(\mu, k)$ given $\underline{Z}_n = \underline{z}_n^*$, where $\underline{z}_n^*$ is given by

$$\underline{z}_n^* = \arg\min_{\underline{z} \in C_{R,n}} (\underline{z} - \mu(\underline{x}_n))^{\top} K_n^{-1} (\underline{z} - \mu(\underline{x}_n)), \qquad (14)$$

with $C_{R,n}$ defined by (10).

Observe that (14) reduces to (12) when $\mu = 0$. Consequently, the mean of the distribution is the predictor $s_{\underline{Z}_n}^*$ from Proposition 3 in this particular case, and is equal to $\mu + s_{\underline{Z}_n}^*$ in general. Moreover, the reGP predictive distribution can be seen as an approximation of (13), where $p(\underline{Z}_n | \underline{Z}_n \in C_{R,n})$ has been replaced by its mode. As discussed earlier, the main advantage of the reGP predictive distribution compared to (13) is its reasonable computational burden since it is a GP. Therefore, it makes it possible to use adaptive strategies for the choice of $R$, as in Section 4. Moreover, it also has appealing theoretical approximation properties, as discussed in Section 3.3.

As discussed in Section 2.1, the standard practice is to select the mean and the covariance functions within a parametric family $\{(\mu_\theta, k_\theta), \theta \in \Theta\}$. In this case, we propose to perform the parameter selection and the relaxation jointly. This is formalized by the following definition of relaxed Gaussian process interpolation.

**Definition 6 (Relaxed-GP predictive distribution; estimated parameters)** Given $x_n \in \mathbb{X}^n$, $z_n \in \mathbb{R}^n$, a relaxation set $R$ (finite union of closed intervals), and parametric families $(\mu_\theta)$ and $(k_\theta)$ as in Section 2.1, the relaxed-GP (reGP) predictive distribution with estimated parameters is the
More precisely, the probabilistic model in Figure 3 interpolates data lying below a selected thresh-
old $t$, and when data are above $t$, the model only keeps the information that the data exceeds $t$.

In the following, we consider the general setting where relaxation is carried out on a set of
function values of the form $R = \bigcup_{j=1}^{J} R_j$, where $R_1, \dots, R_J \subset \mathbb{R}$ are disjoint closed intervals with
non-zero lengths. (The set $R = [t, +\infty)$ was used in the example of Figure 3.)

As above, we shall write $x_n = (x_1, \dots, x_n) \in \mathbb{X}^n$ for a sequence of locations with corresponding
function values $z_n = (z_1, \dots, z_n)^\top \in \mathbb{R}^n$. Then, we introduce the set $C_{R,n} = C_1 \times \dots \times C_n \subset \mathbb{R}^n$ of
relaxed constraints, where

$$
\left\{
\begin{array}{ll}
C_i = R_j & \text{if } z_i \in R_j \text{ for some } j, \\
C_i = \{z_i\} & \text{otherwise.}
\end{array}
\right.
\qquad (9)
$$

Let also

$$
\mathcal{H}_{R,n} = \{h \in \mathcal{H}(\mathcal{X}) \mid h(\underline{x}_n) \in C_{R,n}\}
$$

be the set of relaxed-interpolating functions. The following proposition gives the definition of the
minimum-norm relaxed predictor.

**Proposition 3** *The problem*

$$
\left\{
\begin{array}{ll}
\text{minimize} & \|h\|_{\mathcal{H}(\mathcal{X})} \\
\text{subject to} & h \in \mathcal{H}_{R,n}
\end{array}
\right.
\qquad (10)
$$

has a unique solution given by $s_{\underline{z}_n}^*$, where $\underline{z}_n^*$ is the unique solution of the quadratic problem

$$
\arg\min_{\underline{z} \in C_{R,n}} \underline{z}^{\top} K_n^{-1} \underline{z}. \qquad (11)
$$

An extension of this relaxation scheme to the case of noisy observations, where the function no
longer exactly interpolates the data, is given in Proposition 41.

**3.2 Relaxed Gaussian Process Interpolation**

The main advantage of Gaussian processes is the possibility to obtain not only point predictions
but also predictive distributions. However, Proposition 3 only defines a function approximation.
We now turn relaxed interpolation into a probabilistic model providing predictive distributions
whose mean is not constrained to interpolate data on a given relaxation range $R$. The following
proposition makes a step in this direction.

**Proposition 4** Let $\xi \sim \text{GP}(0, k)$, $\underline{x}_n = (x_1, \dots, x_n) \in \mathbb{X}^n$, $\underline{z}_n \in \mathbb{R}^n$ and $\underline{x}'_m = (x'_1, \dots, x'_m) \in \mathbb{X}^m$ be a set of locations of interest where predictions should be made. Write $\underline{Z}_n = (\xi(x_1), \dots, \xi(x_n))^\top$ and $\underline{Z}'_m = (\xi(x'_1), \dots, \xi(x'_m))^\top$. Then the mode of the probability density function

$$
p(\underline{Z}'_m, \underline{Z}_n | \underline{Z}_n \in C_{R,n}) \qquad (12)
$$

is given by $(s_{\underline{z}_n}^*(\underline{x}'_m), \underline{z}_n^*)$.

In other words, the relaxed interpolation solution of Proposition 3 corresponds to the maximum
a posteriori (MAP) estimate under the predictive model (12). Conditioning on events of the form
$\underline{Z}_n \in C_{R,n}$ has been used in Bayesian statistics as a principled way to handle outliers and model
<figure><img src="image_5.png" /><figcaption>Figure 4: An example of reGP predictive distribution with $R = (-\infty, -1] \cup [1, +\infty)$ on a function $f$ represented in dashed black lines. The solid black lines represent the relaxation thresholds. The problem (15) was solved only in $\underline{z}$ as the parameters of the (constant) mean and ($v = 5/2$ Matérn) covariance functions were held fixed for illustration purposes.</figcaption></figure>

(Gaussian) conditional distribution of $\xi \sim \text{GP}(\mu_{\theta}, k_{\theta})$ given $\underline{Z}_n = \underline{z}_n^*$, where $\underline{z}_n^*$ and $\theta = \hat{\theta}_n$ are obtained jointly by minimizing the negative log-likelihood:

$$ (\hat{\theta}_n, \underline{z}_n^*) = \underset{\theta \in \Theta, \underline{z} \in C_{R,n}}{\operatorname{arg\,min}} \mathcal{L}(\theta; \underline{z}), \qquad (15) $$

with $C_{R,n}$ defined by (10).

**Remark 7 (On minimizing (15) jointly)** Let $\underline{Z}_{n,1}$ be the values within the range $R$, and $\underline{Z}_{n,0}$ the values in $R^c = \mathbb{R} \setminus R$ that are not relaxed. The negative log-likelihood can be written as

$$ \mathcal{L}(\theta; \underline{Z}_n) = -\ln(p(\underline{Z}_{n,0} | \theta)) - \ln(p(\underline{Z}_{n,1} | \theta, \underline{Z}_{n,0})), \qquad (16) $$

where the first term is a goodness-of-fit criterion based on the values in $R^c$, and where the second term can mainly be viewed as an imputation term, which “reshapes” the values in $R$ with the information from $\underline{Z}_{n,0}$. (Note also that $\theta$ appears in the second term. When this term is minimized with respect to $\underline{Z}_{n,1}$, it becomes a parameter selection term that promotes the $\theta$s compatible with the excursions in $C_{R,n}$.)

For illustration, we provide an example of a reGP predictive distribution in Figure 4, with an union of two intervals for the relaxation set $R$.
misspecification; see, e.g., Lewis et al. (2021) and references therein. Such conditional distributions also appear in other settings: in constrained Gaussian process modeling, where constraints encode expert knowledge (Da Veiga and Marrel, 2012; Maatouk and Bay, 2017; López-Lopera et al., 2018), and in classification, ordinal regression, and preference learning problems, where observations impose inequality constraints (Benavoli et al., 2021).

However, the predictive distribution (12) is non-Gaussian since the support of $\underline{Z}_n$ is truncated. In particular, no closed-form expression is available for any of its moments, and sampling requires advanced techniques (e.g., variational, MCMC; see Section 6 for a discussion). Motivated by this observation, we propose instead to build a goal-oriented probabilistic model using the following definition.

**Definition 5 (Relaxed-GP predictive distribution; fixed $\mu$ and $k$)** Given $\underline{x}_n \in \mathbb{X}^n$, $\underline{z}_n \in \mathbb{R}^n$, and a relaxation range $R$ (finite union of closed intervals), the relaxed-GP (reGP) predictive distribution with fixed mean function $\mu$ and covariance function $k$ is defined as the (Gaussian) conditional distribution of $\xi \sim \text{GP}(\mu, k)$ given $\underline{Z}_n = \underline{z}_n^*$, where

$$ \underline{z}_n^* = \arg\min_{\underline{z} \in C_{R,n}} (\underline{z} - \mu(\underline{x}_n))^\mathrm{T} K_n^{-1} (\underline{z} - \mu(\underline{x}_n)), \quad (13) $$

with $C_{R,n}$ defined by (9).

Observe that (13) reduces to (11) when $\mu = 0$. Consequently, the mean of the distribution is the predictor $s_{\underline{z}_n^*}$ from Proposition 3 in this particular case. Moreover, the reGP predictive distribution can be seen as an approximation of (12), where $p(\underline{Z}_n | \underline{Z}_n \in C_{R,n})$ has been replaced by its mode. As discussed earlier, the main advantage of the reGP predictive distribution compared to (12) is its reasonable computational cost since it is a GP. Therefore, it makes it possible to use adaptive strategies for the choice of $R$, as demonstrated in Section 4. Moreover, it also has appealing theoretical approximation properties, as discussed in Section 3.3.

As discussed in Section 2.1, the standard practice is to select the mean and the covariance functions within a parametric family $\{ (\mu_\theta, k_\theta), \theta \in \Theta \}$. Leveraging the connection between (13) and (4), we propose to perform the parameter selection and the relaxation jointly. This is formalized by the following definition of relaxed Gaussian process interpolation.

**Definition 6 (Relaxed-GP predictive distribution; estimated parameters)** Given $\underline{x}_n \in \mathbb{X}^n$, $\underline{z}_n \in \mathbb{R}^n$, a relaxation range $R$ (finite union of closed intervals), and parametric families $(\mu_\theta)$ and $(k_\theta)$ as in Section 2.1, the relaxed-GP (reGP) predictive distribution with estimated parameters is the (Gaussian) conditional distribution of $\xi \sim \text{GP}(\mu_\theta, k_\theta)$ given $\underline{Z}_n = \underline{z}_n^*$, where $\underline{z}_n^*$ and $\theta = \hat{\theta}_n$ are obtained jointly by minimizing the negative log-likelihood:

$$ (\hat{\theta}_n, \underline{z}_n^*) = \arg\min_{\theta \in \Theta, \underline{z} \in C_{R,n}} \mathcal{L}(\theta; \underline{z}), \quad (14) $$

with $C_{R,n}$ defined by (9).

An extension of this definition to the case of noisy observations is presented in Appendix G.

The mean function and the covariance function usually vary continuously with $\theta$ on an open subset of $\mathbb{R}^q$. In this case, it can be shown that the mapping $\theta \mapsto \min_{\underline{z} \in C_{R,n}} \mathcal{L}(\theta; \underline{z})$ is continuous (using, e.g., Theorem 2.2 from Best and Ding, 1995). Consequently, the minimum of (14) is reached when $\theta$ is restricted to a compact subset $\Theta \subset \mathbb{R}^q$.
**Remark 8 (Numerical details)** Minimizing (15) with respect to $z$ falls under the scope of quadratic programming (see, e.g., Nocedal and Wright, 2006) and could be solved efficiently using dedicated algorithms. This suggests that specific algorithms could be developed for the problem. In this work, we simply use a standard L-BFGS-B solver (Byrd et al., 1995) using the gradient of (15).

## 3.3 Convergence Analysis of reGP

In this section, we provide several theoretical results concerning the convergence of the method proposed. This section can be skipped on first reading.

### 3.3.1 KNOWN CONVERGENCE RESULTS ABOUT INTERPOLATION IN RKHS

Recall that the fractional-order Sobolev space $W_2^\beta(\mathbb{R}^d)$, with regularity $\beta \ge 0$, is the space of functions on $\mathbb{R}^d$ defined by

$$W_2^\beta(\mathbb{R}^d) = \{h \in L^2(\mathbb{R}^d), \|h\|_{W_2^\beta(\mathbb{R}^d)} = \|\hat{h}(1 + \|\cdot\|^2)^{\beta/2}\|_{L^2(\mathbb{R}^d)} < +\infty\},$$

where $\hat{h} \in L^2(\mathbb{R}^d)$ is the Fourier transform of $h \in L^2(\mathbb{R}^d)$.

For a given $X \subset \mathbb{R}^d$, define the Sobolev spaces $W_2^\beta(X) = \{h|_X, h \in W_2^\beta(\mathbb{R}^d)\}$ endowed with the norm

$$\|h\|_{W_2^\beta(X)} = \inf_{g \in W_2^\beta(\mathbb{R}^d), g|_X = h} \|g\|_{W_2^\beta(\mathbb{R}^d)}. \quad (17)$$

The following assumption about $X$ will sometimes be used in this section.

**Assumption 1** The domain is non-empty, compact, connected, has locally Lipschitz boundary (see, e.g., Adams and Fournier, 2003, Section 4.9), and is equal to the closure of its interior.

Assumption 1 ensures that the previous definition coincides with other common definitions, and makes it possible to use well-known results from the field of scattered data approximation, by preventing the existence of cusps. Many common domains—such as hyperrectangles or balls, for instance—verify Assumption 1.

A strictly positive-definite reproducing kernel $k: X \times X \to \mathbb{R}$ is said to have regularity $\alpha > 0$ if the associated RKHS $\mathcal{H}(X)$ coincides with $W_2^{\alpha+d/2}(X)$ as a function space, with equivalent norms. As such, the Matérn stationary kernels (4) have correlation functions $r$ whose Fourier transform verifies (see, e.g., Wendland, 2004, Theorem 6.13)

$$C_1 (1 + \|\cdot\|^2)^{-\nu-d/2} \le \hat{r} \le C_2 (1 + \|\cdot\|^2)^{-\nu-d/2}$$

for some $C_2 \ge C_1 > 0$, and have therefore Sobolev regularity $\alpha = \nu$ on $\mathbb{R}^d$ (see, e.g., Wendland, 2004, Corollary 10.13) and consequently also on $X$, using (17) and Lemma 24. Other examples are given by Wendland (2004), for instance.

We now recall a classical convergence result about interpolation in RKHS with evaluation points in a bounded domain. Consider a kernel $k: X \times X \to \mathbb{R}$, and let $(x_n)_{n \ge 1} \in X^\mathbb{N}$ be a sequence of distinct points. The following property (a minor reformulation of Theorem 4.1 of Arcangéli et al.
<figure><img src="image_5.png" /><figcaption>Figure 4: An example of reGP predictive distribution with $R = (-\infty, -1] \cup [1, +\infty)$ on a function $f$ represented in dashed black lines. The solid black lines represent the relaxation thresholds. The problem (14) was solved only in $z$ as the parameters of the (constant) mean and ($v = 5/2$ Matérn) covariance functions were held fixed for illustration purposes.</figcaption></figure>

**Remark 7 (On minimizing (14) jointly)** Let $Z_{n,1}$ be the values within the relaxation range $R$, and $Z_{n,0}$ the values in $R^c = \mathbb{R} \setminus R$ that are not relaxed. The negative log-likelihood can be written as

$$ \mathcal{L}(\theta; Z_n) = -\ln(p(Z_{n,0} | \theta)) - \ln(p(Z_{n,1} | \theta, Z_{n,0})), $$

where the first term is a goodness-of-fit criterion based on the values in $R^c$, and where the second term can mainly be viewed as an imputation term, which “reshapes” the values in $R$ with the information from $Z_{n,0}$. (Note also that $\theta$ appears in the second term. When this term is minimized with respect to $Z_{n,1}$, it becomes a parameter selection term that promotes the $\theta$s compatible with the excursions in $C_{R,n}$.)

For illustration, we provide an example of a reGP predictive distribution in Figure 4, with an union of two intervals for the relaxation range $R$.

**Remark 8 (Numerical details)** Minimizing (14) with respect to $z$ falls under the scope of quadratic programming (see, e.g., Nocedal and Wright, 2006) and could be solved efficiently using dedicated algorithms. This suggests that specific algorithms could be developed for the problem. In this work, we simply use a standard SLSQP solver (see Schittkowski, 1982) using the gradient of (14). The SciPy implementation was used with default parameters. Good performances were obtained by the following parameter initialization procedure: first, initialize the GP parameter $\theta$ using common rules while excluding the most extreme observations, and then optimize the relaxed observations $z$ using quadratic programming.
(2007)) gives error bounds that depend on the Sobolev regularity of $k$ and the so-called fill distance of $x_n \in \mathbb{X}^n$, defined by

$$h_n = \sup_{x \in \mathbb{X}} \min_{1 \le i \le n} \|x - x_i\|. \qquad (18)$$

**Proposition 9** Let $k$ be a reproducing kernel with regularity $\alpha > 0$. If $\mathbb{X}$ verifies Assumption 1, then

$$\sup_{x \in \mathbb{X}} \sigma_n(x) \lesssim h_n^{\alpha}, \quad n \ge 1, \qquad (19)$$

where $\lesssim$ denotes inequality up to a constant, that does not depend on $(x_n)_{n \ge 1}$.

Using (9) and Proposition 9, this yields the following uniform bound.

**Corollary 10** Let $k$ be a reproducing kernel with regularity $\alpha > 0$, $\mathcal{H}(\mathbb{X})$ the RKHS generated by $k$, and let $f \in \mathcal{H}(\mathbb{X})$. As above, let $s_{z_n}$ be the solution of (8) for $z_n = (f(x_1), \dots, f(x_n))^\top$, $n \ge 1$. If $\mathbb{X}$ verifies Assumption 1, then

$$\|f - s_{z_n}\|_{L^\infty(\mathbb{X})} \lesssim h_n^{\alpha} \|f\|_{\mathcal{H}(\mathbb{X})}. \qquad (20)$$

### 3.3.2 CONVERGENCE RESULTS FOR REGP

Let $k: \mathbb{X} \times \mathbb{X} \to \mathbb{R}$ be a continuous strictly positive-definite reproducing kernel. In this section, we consider the zero-mean reGP predictive distribution obtained from $\xi \sim \text{GP}(0, k)$, with relaxed interpolation constraints on a union $R = \cup_{j=1}^q R_j$ of disjoints closed intervals $R_j$ with non-zero length. Let $\mathcal{H}(\mathbb{X})$ be the RKHS attached to $k$, $f \in \mathcal{H}(\mathbb{X})$, and consider a sequence $(x_n)_{n \ge 1} \in \mathbb{X}^\mathbb{N}$ of distinct points. Furthermore, define (the) regions $\mathbb{X}_j = \{x \in \mathbb{X}, f(x) \in R_j\}$ for $1 \le j \le q$ and $\mathbb{X}_0 = \mathbb{X} \setminus \cup_{j \ge 1} \mathbb{X}_j$. We give results about the limit of the sequence of reGP predictive distributions that suggest an improved fit in $\mathbb{X}_0$.

Let $s_{R,n} = s_{z_n}^*$ be the relaxed predictor from Proposition 3 based on $(x_1, \dots, x_n)$ and $(f(x_1), \dots, f(x_n))^\top$, $n \ge 1$. The following proposition establishes the limit behavior of $(s_{R,n})_{n \ge 1}$.

**Proposition 11** Let $\mathbb{U} \subset \mathbb{X}$ and let $\mathcal{H}_{R,\mathbb{U}}$ denote the set of functions $h \in \mathcal{H}(\mathbb{X})$ such that, for all $x \in \mathbb{U}$,

$$\left\{ \begin{array}{ll} h(x) \in R_j & \text{if } f(x) \in R_j \text{ for some } j, \\ h(x) = f(x) & \text{otherwise.} \end{array} \right. \qquad (21)$$

Then the problem

$$\min_{h \in \mathcal{H}_{R,\mathbb{U}}} \|h\|_{\mathcal{H}(\mathbb{X})} \qquad (22)$$

has a unique solution denoted by $s_{R,\mathbb{U}}$. Moreover,

$$s_{R,n} \xrightarrow{\mathcal{H}(\mathbb{X})} s_{R,\mathbb{U}}, \qquad (23)$$

with $\mathbb{U}$ the closure of $\{x_n\}$.
### 3.3 Convergence Analysis of reGP

In this section, we provide theoretical results concerning the convergence of the method proposed. This section can be skipped on first reading.

#### 3.3.1 KNOWN CONVERGENCE RESULTS ABOUT INTERPOLATION IN RKHS

Recall that the fractional-order Sobolev space $W_2^β(\mathbb{R}^d)$, with regularity $β ≥ 0$, is the space of functions on $\mathbb{R}^d$ defined by

$$W_2^β(\mathbb{R}^d) = \{h \in L^2(\mathbb{R}^d), \|h\|_{W_2^β(\mathbb{R}^d)}^2 = \int_{\mathbb{R}^d} (1 + \|\omega\|^2)^β |\hat{h}(\omega)|^2 d\omega < +\infty\},$$

where $\hat{h} \in L^2(\mathbb{R}^d)$ is the Fourier transform of $h \in L^2(\mathbb{R}^d)$. For $β > d/2$, the usual identification with a set of continuous functions given by the Sobolev embedding theorem will be used in this article.

For a given $X \subset \mathbb{R}^d$, define the Sobolev spaces $W_2^β(X) = \{h|_X, h \in W_2^β(\mathbb{R}^d)\}$ endowed with the norm

$$\|h\|_{W_2^β(X)} = \inf_{g \in W_2^β(\mathbb{R}^d), g|_X = h} \|g\|_{W_2^β(\mathbb{R}^d)}. \quad (15)$$

The following assumption about $X$ will sometimes be used in this section.

**Assumption 9** The domain is non-empty, compact, connected, has locally Lipschitz boundary (see, e.g., Adams and Fournier, 2003, Section 4.9), and is equal to the closure of its interior.

Assumption 9 ensures that the previous definition coincides with other common definitions, and makes it possible to use well-known results from the field of scattered data approximation, by preventing the existence of cusps. Many common domains—such as hyperrectangles or balls, for instance—satisfy Assumption 9.

A strictly positive-definite reproducing kernel $k: X \times X \to \mathbb{R}$ is said to have regularity $α > 0$ if the associated RKHS $\mathcal{H}(X)$ coincides with $W_2^{α+d/2}(X)$ as a function space, with equivalent norms. As such, the Matérn stationary kernels (3) have correlation functions $r$ whose Fourier transform satisfies (see, e.g., Wendland, 2004, Theorem 6.13)

$$C_1 (1 + \|\cdot\|^2)^{-\nu-d/2} \le \hat{r} \le C_2 (1 + \|\cdot\|^2)^{-\nu-d/2}$$

for some $C_2 \ge C_1 > 0$, and have therefore Sobolev regularity $α = ν$ on $\mathbb{R}^d$ (see, e.g., Wendland, 2004, Corollary 10.13) and consequently also on $X$, using (15) and Lemma 31. Other examples are given by Wendland (2004), for instance.

We now recall a classical convergence result about interpolation in RKHS with evaluation points in a bounded domain. Consider a kernel $k: X \times X \to \mathbb{R}$, and let $(x_n)_{n \ge 1} \in \mathbb{X}^N$ be a sequence of distinct points. The following property (a minor reformulation of Theorem 4.1 of Arcangéli et al., 2007) gives error bounds that depend on the Sobolev regularity of $k$ and the so-called fill distance of $x_n \in \mathbb{X}^n$, defined by

$$h_n = \sup_{x \in \mathbb{X}} \min_{1 \le i \le n} \|x - x_i\|.$$
In particular, when {$x_n$} is dense in $\mathbb{X}$, then $\mathbb{U} = \mathbb{X}$ and $(s_{R,n})_{n \ge 1}$ converges to $s_{R,\mathbb{X}}$, which is the minimal-norm element of the set $\mathcal{H}_{R,\mathbb{X}}$.

The next proposition tells us that the interpolation error on $\mathbb{X}_0$ can be bounded by a term that depends on the norm of $s_{R,\mathbb{X}}$.

**Proposition 12** For any $x \in \mathbb{X}_0$ and $n \ge 1$,

$$|f(x) - s_{R,n}(x)| \le 2\sigma_{n,0}(x)\|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})}, \quad (24)$$

where $\sigma_{n,0}$ is the power function obtained using only points in $\mathbb{X}_0$ for predictions.

This yields the following error bounds when the design is dense.

**Proposition 13** Suppose that {$x_n$} is dense and that $k$ has regularity $\alpha > 0$. Let $B \subset \mathbb{X}_0$ verify Assumption 1. Then, for all $n \ge 1$,

$$\|f - s_{R,n}\|_{L^{\infty}(B)} \lesssim h_n^{\alpha} \|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})}. \quad (25)$$

Let $d(y, A)$ be the distance of $y \in \mathbb{R}$ to $A \subset \mathbb{R}$. For $j \ge 1$, $x \in \mathbb{X}_j$, and for all $n \ge 1$:

$$d(s_{R,n}(x), R_j) \lesssim h_n^{\alpha} \|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})} \quad \text{if } \alpha < 1, \quad (26)$$

$$d(s_{R,n}(x), R_j) \lesssim \sqrt{(|\ln(h_n)| + 1)} h_n \|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})} \quad \text{if } \alpha = 1, \quad (27)$$

and

$$d(s_{R,n}(x), R_j) \lesssim h_n \|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})} \quad \text{if } \alpha > 1, \quad (28)$$

where $\lesssim$ denotes inequality up to a constant, that does not depend on $f$, $n$ or $(x_n)$.

Finally, we investigate the following question: how large can be the norm of $f$ compared to the approximation $\|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})}$?

**Proposition 14** Suppose that $k$ has regularity $\alpha > 0$ and that there exists some $j \ge 1$ such that $\mathbb{X}_j$ has a non-empty interior. We have

$$\sup_{h \in \mathcal{H}_{R,\mathbb{X}}} \|h\|_{\mathcal{H}(\mathbb{X})} = +\infty, \quad (29)$$

with $\mathcal{H}_{R,\mathbb{X}}$ given by (21) for $f \in \mathcal{H}(\mathbb{X})$.

This result shows that the norm reduction obtained by approximating $f$ with relaxed interpolation constraints can therefore be arbitrarily high in the finite-smoothness case. A stronger version of Proposition 14 for the special case where $R = [t, +\infty)$ can be derived, and shows that

$$\sup_{h \in \mathcal{H}_{R,\mathbb{X}}} \|h\|_{L^{\infty}(\mathbb{X})} = +\infty.$$

Overall, no matter the element of $\mathcal{H}_{R,\mathbb{X}}$ at hand, reGP converges to a function $s_{R,\mathbb{X}}$ which: coincides with $f$ on $\mathbb{X}_0$, verifies $f(x) \in R_j \Leftrightarrow s_{R,\mathbb{X}}(x) \in R_j$ for all $x \in \mathbb{X}$, and is “nicer” than $f$ in the sense of $\|\cdot\|_{\mathcal{H}(\mathbb{X})}$. Furthermore, reGP yields error bounds carrying the norm of $s_{R,\mathbb{X}}$, which can be arbitrarily smaller than the norm of $f$ in the case of a finite-smoothness covariance function.
**Proposition 10** Let $k$ be a reproducing kernel with regularity $\alpha > 0$. If $\mathcal{X}$ satisfies Assumption 9, then

$$ \sup_{x \in \mathcal{X}} \sigma_n(x) \lesssim h_n^{\alpha}, \quad n \ge 1, \qquad (16) $$

where $\lesssim$ denotes inequality up to a constant, that does not depend on $(x_n)_{n \ge 1}$.

Using (8) and Proposition 10, this yields the following uniform bound.

**Corollary 11** Let $k$ be a reproducing kernel with regularity $\alpha > 0$, $\mathcal{H}(\mathcal{X})$ the RKHS generated by $k$, and let $f \in \mathcal{H}(\mathcal{X})$. As above, let $s_{\underline{z}_n}$ be the solution of (7) for $\underline{z}_n = (f(x_1), \dots, f(x_n))^{\top}$, $n \ge 1$. If $\mathcal{X}$ satisfies Assumption 9, then

$$ \|f - s_{\underline{z}_n}\|_{L^{\infty}(\mathcal{X})} \lesssim h_n^{\alpha} \|f\|_{\mathcal{H}(\mathcal{X})}. \qquad (17) $$

### 3.3.2 CONVERGENCE OF RELAXED INTERPOLATION WITH FIXED COVARIANCE FUNCTION

Let $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ be a continuous strictly positive-definite reproducing kernel. In this section, we consider the zero-mean reGP predictive distribution obtained from $\xi \sim \text{GP}(0, k)$, with relaxed interpolation constraints on a union $R = \cup_{j=1}^q R_j$ of disjoints closed intervals $R_j$ with non-zero length. Let $\mathcal{H}(\mathcal{X})$ be the RKHS attached to $k$, $f \in \mathcal{H}(\mathcal{X})$, and consider a sequence $(x_n)_{n \ge 1} \in \mathcal{X}^\mathbb{N}$ of distinct points. Furthermore, define the regions $\mathcal{X}_j = \{x \in \mathcal{X}, f(x) \in R_j\}$ for $1 \le j \le q$ and $\mathcal{X}_0 = \mathcal{X} \setminus \cup_{j \ge 1} \mathcal{X}_j$. We give results about the limit of the sequence of reGP predictive distributions that indicate an improved fit in $\mathcal{X}_0$.

Let $s_{R,n} = s_{\underline{z}_n}^*$ be the relaxed predictor from Proposition 3 based on $(x_1, \dots, x_n)$ and $(f(x_1), \dots, f(x_n))^{\top}$, $n \ge 1$. The following proposition establishes the limit behavior of $(s_{R,n})_{n \ge 1}$.

**Proposition 12** Let $\mathcal{U} \subset \mathcal{X}$ and let $\mathcal{H}_{\mathcal{R},\mathcal{U}}$ denote the set of functions $h \in \mathcal{H}(\mathcal{X})$ such that, for all $x \in \mathcal{U}$,

$$ \left\{ \begin{array}{ll} h(x) \in R_j & \text{if } f(x) \in R_j \text{ for some } j, \\ h(x) = f(x) & \text{otherwise.} \end{array} \right. \qquad (18) $$

Then the problem

$$ \min_{h \in \mathcal{H}_{\mathcal{R},\mathcal{U}}} \|h\|_{\mathcal{H}(\mathcal{X})} \qquad (19) $$

has a unique solution denoted by $s_{R,\mathcal{U}}$. Moreover, when $\mathcal{U}$ is the closure of $\{x_n\}$,

$$ s_{R,n} \xrightarrow{\mathcal{H}(\mathcal{X})} s_{R,\mathcal{U}}. $$

In particular, when $\{x_n\}$ is dense in $\mathcal{X}$, then $\mathcal{U} = \mathcal{X}$ and $(s_{R,n})_{n \ge 1}$ converges to $s_{R,\mathcal{X}}$, which is the minimal-norm element of the set $\mathcal{H}_{\mathcal{R},\mathcal{X}}$.

The next proposition tells us that the interpolation error on $\mathcal{X}_0$ can be bounded by a term that depends on the norm of $s_{R,\mathcal{X}}$.

**Proposition 13** For any $x \in \mathcal{X}_0$ and $n \ge 1$,

$$ |f(x) - s_{R,n}(x)| \le 2\sigma_{n,0}(x) \|s_{R,\mathcal{X}}\|_{\mathcal{H}(\mathcal{X})}, \qquad (20) $$

where $\sigma_{n,0}$ is the power function obtained using only points in $\mathcal{X}_0$ for predictions and with the convention that $\sigma_{n,0}^2(x) = k(x, x)$ if there are no points.
**Remark 15** Note that $\sigma_n \le \sigma_{n,0}$ due to the projection residuals interpretation. Empirical and theoretical results about the screening effect (see, e.g., Stein, 2011; Bao et al., 2020), suggests that $\sigma_n \approx \sigma_{n,0}$, if $k$ has smoothness $\alpha > 0$. In this case, observe that—no matter the element of $\mathcal{H}_{R,X}$ at hand—the bound (24) is larger by only a small factor compared to (9) with $f = s_{R,X}$. (However, to the best of our knowledge, no result exists concerning the screening effect for arbitrary designs.)

**Remark 16** The equality (29) does not hold in general for infinitely smooth covariance functions. For instance, Steinwart et al. (2006, Corollary 3.9) show that $\mathcal{H}_{R,X} = \{f\}$ if the interior of $X_0$ is not empty and $k$ is the squared-exponential covariance function (i.e. (4), with $v \to \infty$).

# 4. Choice of the Relaxation Set

## 4.1 Towards Goal-Oriented Cross-Validation

The framework of reGP makes it possible to predict a function $f$ from point evaluations of $f$. Suppose we are specifically interested in obtaining good predictive distributions in a range $Q \subset \mathbb{R}$ of function values, and accept degraded predictions outside this range. To achieve this goal, the idea of reGP is to relax interpolation constraints. Naturally, it makes sense to relax interpolation constraints outside the range $Q$ but it could happen that relaxing interpolation constraints does not improve predictive distributions on $Q$. Therefore, the question arises as to how to automatically select a range $R$ in $\mathbb{R} \setminus Q$, on which interpolation constraints should be relaxed.

In the following, we put $R^{(0)} = \mathbb{R} \setminus Q$, and we view the relaxation set $R$ as a parameter of the reGP model, which has to be chosen in $R^{(0)}$ along with the parameters $\theta$ of the underlying GP $\xi$. A first idea for the selection of $R$ is to rely on the standard leave-one-out cross-validation approach to select the parameters of a GP (Dubrule, 1983; Rasmussen and Williams, 2006; Zhang and Wang, 2010). Using the formalism of scoring rules (see, e.g., Gneiting and Raftery, 2007; Petit et al., 2021), selecting parameters by a leave-one-out approach amounts to minimizing a selection criterion written as

$$J_n(R) = \frac{1}{n} \sum_{i=1}^{n} S(P_{R,n,-i}, f(x_i)), \qquad (30)$$

where $P_{R,n,-i}$ is reGP predictive distribution with data $z_{n,-i} = (z_1, \dots, z_{i-1}, z_{i+1}, \dots, z_n)$ and relaxation set $R$. The function $S$ in (30) is a scoring rule, that is, a function $S: \mathcal{P} \times \mathbb{R} \to \mathbb{R} \cup \{-\infty, +\infty\}$, acting on a class $\mathcal{P}$ of probability distributions on $\mathbb{R}$, such that $S(P, z)$ assigns a loss for choosing a predictive distribution $P \in \mathcal{P}$, while observing $z \in \mathbb{R}$. Scoring rules make it possible to quantify the quality of probabilistic predictions.

Since the user is not specifically interested in good predictive distributions in $R^{(0)}$, validating the model on $R^{(0)}$ should not be a primary focus. However, simply restricting the sum (30) by removing indices $i$ such that $f(x_i) \in R^{(0)}$ would make it impossible to assess if the model is good at predicting that $f(x) \in R^{(0)}$ for a given $x \in \mathbb{X}$. For instance, in the case of minimization, with $Q = (-\infty, t^{(0)})$ and $R^{(0)} = [t^{(0)}, \infty)$, it is important to identify the regions corresponding to $f$ being above $t^{(0)}$, even if we are not interested in accurate predictions above $t^{(0)}$, because we expect that an optimization algorithm should avoid the exploration of these regions.

In the next section, we propose instead to keep the whole leave-one-out sum (30), but to choose a scoring rule $S$ that serves our goal-oriented approach.
This yields the following error bounds when the design is dense.

**Proposition 14** Suppose that $\mathcal{X}$ is compact, $\{x_n\}$ is dense, and that $k$ has regularity $\alpha > 0$. Let $B \subset \mathcal{X}_0$ satisfy Assumption 9. Then, for all $n \ge 1$,

$$ \|f - s_{R,n}\|_{L^{\infty}(B)} \lesssim h_n^{\alpha} \|s_{R,\mathcal{X}}\|_{\mathcal{H}(\mathcal{X})}. \quad (21) $$

Let $d(y, A)$ be the distance of $y \in \mathbb{R}$ to $A \subset \mathbb{R}$. For $j \ge 1$, $x \in \mathcal{X}_j$, and for all $n \ge 1$:

$$ d(s_{R,n}(x), R_j) \lesssim h_n^{\alpha} \|s_{R,\mathcal{X}}\|_{\mathcal{H}(\mathcal{X})} \quad \text{if } \alpha < 1, $$

$$ d(s_{R,n}(x), R_j) \lesssim (|\ln(h_n)| + 1)^{1/2} h_n \|s_{R,\mathcal{X}}\|_{\mathcal{H}(\mathcal{X})} \quad \text{if } \alpha = 1, $$

and

$$ d(s_{R,n}(x), R_j) \lesssim h_n \|s_{R,\mathcal{X}}\|_{\mathcal{H}(\mathcal{X})} \quad \text{if } \alpha > 1, $$

where $\lesssim$ denotes inequality up to a constant, that does not depend on $f$, $n$, $x$ or $(x_n)$.

Finally, we investigate the following question: how large can the norm of $f$ be compared to that of the approximation $\|s_{R,\mathcal{X}}\|_{\mathcal{H}(\mathcal{X})}$?

**Proposition 15** Suppose that $k$ has regularity $\alpha > 0$ and that there exists some $j \ge 1$ such that $\mathcal{X}_j$ has a non-empty interior. We have

$$ \sup_{h \in \mathcal{H}_{R, \mathcal{X}}} \|h\|_{\mathcal{H}(\mathcal{X})} = +\infty, \quad (22) $$

with $\mathcal{H}_{R,\mathcal{X}}$ given by (18) for $f \in \mathcal{H}(\mathcal{X})$.

This result shows that the norm reduction obtained by approximating $f$ with relaxed interpolation constraints can therefore be arbitrarily high in the finite-smoothness case. A stronger version of Proposition 15 for the special case where $R = [t, +\infty)$ can be derived, and shows that

$$ \sup_{h \in \mathcal{H}_{R, \mathcal{X}}} \|h\|_{L^{\infty}(\mathcal{X})} = +\infty. $$

Overall, no matter the element of $\mathcal{H}_{R,\mathcal{X}}$ at hand, reGP converges to a function $s_{R,\mathcal{X}}$ which: coincides with $f$ on $\mathcal{X}_0$, satisfies $f(x) \in R_j \Leftrightarrow s_{R,\mathcal{X}}(x) \in R_j$ for all $x \in \mathcal{X}$, and is “nicer” than $f$ in the sense of $\|\cdot\|_{\mathcal{H}(\mathcal{X})}$. Furthermore, reGP yields error bounds carrying the norm of $s_{R,\mathcal{X}}$, which can be arbitrarily smaller than the norm of $f$ in the case of a finite-smoothness covariance function.

**Remark 16** Note that $\sigma_n \le \sigma_{n,0}$ due to the projection residuals interpretation. Empirical and theoretical results about the screening effect (see, e.g., Stein, 2011; Bao et al., 2020), suggests that $\sigma_n \simeq \sigma_{n,0}$, if $k$ has smoothness $\alpha > 0$. In this case, observe that—no matter the element of $\mathcal{H}_{R,\mathcal{X}}$ at hand—the bound (20) is larger by only a small factor compared to (8) with $f = s_{R,\mathcal{X}}$. (However, to the best of our knowledge, no result exists concerning the screening effect for arbitrary designs.)

**Remark 17** (22) does not hold in general for infinitely smooth covariance functions. For instance, Steinwart et al. (2006, Corollary 3.9) show that $\mathcal{H}_{R,\mathcal{X}} = \{f\}$ if the interior of $\mathcal{X}_0$ is not empty and $k$ is the squared-exponential covariance function, i.e., (3), with $\nu \to \infty$.

Appendix B offers insights into the convergence of relaxed interpolation for functions that lie outside the RKHS attached to the covariance. Additionally, we present results on the behavior of the approximation when the covariance parameters are not fixed.
## 4.2 Truncated Continuous Ranked Probability Score

An appealing class of scoring rules for goal-oriented predictive distributions is the class of weighted scoring rules for binary predictors (Gneiting and Raftery, 2007; Matheson and Winkler, 1976), which may be written as

$$S(P, z) = \int_{-\infty}^{+\infty} s(F_P(u), 1_{z \le u}) \mu(du), \qquad (31)$$

where $s: [0, 1] \times \{0, 1\} \to \mathbb{R} \cup \{-\infty, +\infty\}$ is a scoring rule for binary predictors, and $\mu$ is a Borel measure on $\mathbb{R}$. A well-known instance of (31) is the continuous ranked probability score (Gneiting et al., 2005) written as

$$S^{\text{CRPS}}(P, z) = \int_{-\infty}^{+\infty} (F_P(u) - 1_{z \le u})^2 du,$$

which is obtained by choosing the Brier score for $s$ and the Lebesgue measure for $\mu$.

For the case where we are specifically interested in obtaining good predictive distributions in a range of interest $Q \subset \mathbb{R}$, we propose to use the following scoring rule, which we call *truncated continuous ranked probability score (tCRPS)*:

$$S_{Q}^{\text{tCRPS}}(P, z) = \int_{Q} (F_{P}(u) - 1_{z \le u})^{2} du. \qquad (32)$$

This scoring rule, proposed by Lerch and Thorarinsdottir (2013) in a different context, reduces to $S^{\text{CRPS}}$ when $Q = \mathbb{R}$. It can be seen as a special case of the weighted CRPS (Matheson and Winkler, 1976; Gneiting and Raftery, 2007; Gneiting and Ranjan, 2011), in which the indicator function $1_Q$ plays the role of the weight function—in other words, the measure $\mu$ in (31) has density $1_Q$ with respect to Lebesgue's measure.

Consider for instance the case $Q = (-\infty, t^{(0)})$:

$$S_{Q}^{\text{tCRPS}}(P, z) = \int_{-\infty}^{t^{(0)}} (F_{P}(u) - 1_{z \le u})^{2} du. \qquad (33)$$

The upper endpoint $t^{(0)}$ of the range will be referred to as the validation threshold. Note that, in this case, $S_{Q}^{\text{tCRPS}}(P, z)$ does not depend on the specific value of $z$ when $z$ is above the validation threshold. This scoring rule is thus well suited to the problem of measuring the performance of a predictive distribution in such a way as to fully assess the goodness-of-fit of the distribution when the true value is below a threshold, and only ask that the support of the predictive distribution is concentrated above the threshold when the true value is above the threshold.

We provide in Appendix A some properties of the scoring rule (32) and closed-form expressions for the case where $Q$ is an interval (or a finite union of intervals) and $P$ is Gaussian. To the best of our knowledge, these expressions are new.

## 4.3 Choosing the Relaxation Set using the tCRPS Scoring Rule

Given a range of interest $Q$, the tCRPS scoring rule makes it possible to derive a goal-oriented leave-one-out selection criterion, that we shall call the LOO-tCRPS criterion:

$$J_n(R) = \frac{1}{n} \sum_{i=1}^{n} S_{Q}^{\text{tCRPS}}(P_{R,n,-i}, f(x_i)). \qquad (34)$$
# 4 Choice of the Relaxation Range

## 4.1 Towards Goal-Oriented Cross-Validation

The framework of reGP makes it possible to predict a function $f$ from point evaluations of $f$. Suppose we are specifically interested in obtaining good predictive distributions in a range $Q \subset \mathbb{R}$ of function values, and accept degraded predictions outside this range. To achieve this goal, the idea of reGP is to relax interpolation constraints. Naturally, it makes sense to relax interpolation constraints outside the range of interest $Q$, but it can happen that relaxing all the interpolation constraints outside $Q$ does not improve predictive distributions on $Q$. Therefore, the question arises as to how to automatically select a range $R$ in $\mathbb{R} \setminus Q$, on which interpolation constraints should be relaxed.

In the following, we put $R^{(0)} = \mathbb{R} \setminus Q$, and we view the relaxation range $R$ as a parameter of the reGP model, which has to be chosen in $R^{(0)}$ along with the parameters $\theta$ of the underlying GP $\xi$. A first idea for the selection of $R$ is to rely on the standard leave-one-out cross-validation approach to select the parameters of a GP (Dubrule, 1983; Rasmussen and Williams, 2006; Zhang and Wang, 2010). Using the formalism of scoring rules (see, e.g., Gneiting and Raftery, 2007; Petit et al., 2023), selecting parameters by a leave-one-out approach amounts to minimizing a selection criterion written as

$$J_n(R) = \frac{1}{n} \sum_{i=1}^{n} S(P_{R,n,-i}, f(x_i)), \qquad (23)$$

where $P_{R,n,-i}$ is the reGP predictive distribution with data $z_{n,-i} = (z_1, \dots, z_{i-1}, z_{i+1}, \dots, z_n)$ and relaxation range $R$. The function $S$ in (23) is a scoring rule, that is, a function $S: \mathcal{P} \times \mathbb{R} \to \mathbb{R} \cup \{-\infty, +\infty\}$, acting on a class $\mathcal{P}$ of probability distributions on $\mathbb{R}$, such that $S(P, z)$ assigns a loss for choosing a predictive distribution $P \in \mathcal{P}$, while observing $z \in \mathbb{R}$. Scoring rules make it possible to quantify the quality of probabilistic predictions.

Since the user is not specifically interested in good predictive distributions in $R^{(0)}$, validating the model on $R^{(0)}$ should not be a primary focus. However, simply restricting the sum (23) by removing indices $i$ such that $f(x_i) \in R^{(0)}$ would make it impossible to assess if the model is good at predicting that $f(x) \in R^{(0)}$ for a given $x \in \mathbb{X}$. For instance, in the case of minimization, with $Q = (-\infty, t^{(0)})$ and $R^{(0)} = [t^{(0)}, +\infty)$, it is important to identify the regions corresponding to $f$ being above $t^{(0)}$, even if we are not interested in accurate predictions above $t^{(0)}$, because we expect that an optimization algorithm should avoid the exploration of these regions.

In the next section, we propose instead to keep the whole leave-one-out sum (23), but to choose a scoring rule $S$ that serves our goal-oriented approach.

## 4.2 Truncated Continuous Ranked Probability Score

An appealing class of scoring rules for goal-oriented predictive distributions is the class of weighted scoring rules for binary predictors (Gneiting and Raftery, 2007; Matheson and Winkler, 1976), which may be written as

$$S(P, z) = \int_{-\infty}^{+\infty} s(F_P(u), \mathbb{1}_{z \le u}) \mu(du), \qquad (24)$$

where $s: [0, 1] \times \{0, 1\} \to \mathbb{R} \cup \{-\infty, +\infty\}$ is a scoring rule for binary predictors, $F_P$ is the cumulative distribution function of $P$, and $\mu$ is a Borel measure on $\mathbb{R}$. A well-known instance of (24) is the
<figure><img src="image_6.png" /><figcaption>Figure 5: Illustration of the choice of a relaxation set. The range of interest $\mathcal{Q}$ is determined by the threshold $t^{(0)}$. The relaxation set $\mathcal{R}$ corresponding to the region above $t$ has been obtained by the procedure described in Section 4.3.</figcaption></figure>

Using (34), we suggest the following procedure to select a reGP model. First, choose a sequence of nested candidate relaxation sets $R^{(0)} \supset R^{(1)} \supset \dots \supset R^{(G)} = \emptyset$. The next step is the computation of $J_n(R^{(g)})$, $g = 0, \dots, G$, which involves the predictive distributions $P_{R^{(g)}, n, -i}$.

In principle, (15) should be solved again each time a data point $(x_i, z_i)$ is removed, to obtain a pair $(\hat{\theta}_{n,-i}^{(g)}, \underline{z}_{n,-i}^{(g)})$ and then the corresponding reGP distribution $P_{R^{(g)}, n, -i}$. To alleviate computational cost, a simple idea is to rely on the fast leave-one-out formulas (Dubrule, 1983) for Gaussian processes: for each set $R^{(g)}$, solve (15) to obtain $\hat{\theta}_n^{(g)}$ and $\underline{z}_n^{(g)} = (z_1^{(g)}, \dots, z_n^{(g)})^\top$, and then compute the conditional distributions $\xi(x_i) | \{\xi(x_j) = z_j^{(g)}, j \neq i\}$, where $\xi \sim \text{GP}(\mu, k)$, and where $\mu$ and $k$ have parameter $\hat{\theta}_n^{(g)}$, using the fast leave-one-out formulas. By doing so, we neglect the difference between $\hat{\theta}_{n,-i}^{(g)}$ and $\hat{\theta}_n^{(g)}$ and the difference between $\underline{z}_{n,-i}^{(g)}$ and the vector $(z_1^{(g)}, \dots, z_{i-1}^{(g)}, z_{i+1}^{(g)}, \dots, z_n^{(g)})^\top$.

The procedure ends by choosing the relaxation set $R^{(g)}$ that achieves the best LOO-tCRPS value.

Figure 5 illustrates the selection of the relaxation set used in Figure 3.
continuous ranked probability score (Gneiting et al., 2005) written as

$$S^{\text{CRPS}}(P, z) = \int_{-\infty}^{+\infty} (F_P(u) - 1_{z \le u})^2 du,$$

which is obtained by choosing the Brier score for s and the Lebesgue measure for $\mu$.

For the case where we are specifically interested in obtaining good predictive distributions in a range of interest $Q \subset \mathbb{R}$, we propose to use the following scoring rule, which we call *truncated continuous ranked probability score (tCRPS):*

$$S_{Q}^{\text{tCRPS}}(P, z) = \int_{Q} (F_{P}(u) - 1_{z \le u})^{2} du. \qquad (25)$$

This scoring rule, proposed by Lerch and Thorarinsdottir (2013) in a different context, reduces to $S^{\text{CRPS}}$ when $Q = \mathbb{R}$. It can be seen as a special case of the weighted CRPS (Matheson and Winkler, 1976; Gneiting and Raftery, 2007; Gneiting and Ranjan, 2011), in which the indicator function $1_Q$ plays the role of the weight function—in other words, the measure $\mu$ in (24) has density $1_Q$ with respect to Lebesgue’s measure.

Consider for instance the case $Q = (-\infty, t^{(0)})$:

$$S_{Q}^{\text{tCRPS}}(P, z) = \int_{-\infty}^{t^{(0)}} (F_{P}(u) - 1_{z \le u})^{2} du.$$

The upper endpoint $t^{(0)}$ of the range will be referred to as the validation threshold. Note that, in this case, $S_{Q}^{\text{tCRPS}}(P, z)$ does not depend on the specific value of $z$ when $z$ is above the validation threshold. This scoring rule is thus well suited to the problem of measuring the performance of a predictive distribution in such a way as to fully assess the goodness-of-fit of the distribution when the true value is below a threshold, and only ask that the support of the predictive distribution is concentrated above the threshold when the true value is above the threshold.

We provide in Appendix A some properties of the scoring rule (25) and closed-form expressions for the case where $Q$ is an interval (or a finite union of intervals) and $P$ is Gaussian. To the best of our knowledge, these expressions are new.

## 4.3 Choosing the Relaxation Range using the tCRPS Scoring Rule

Given a range of interest $Q$, the tCRPS scoring rule makes it possible to derive a goal-oriented leave-one-out selection criterion for the relaxation range $R$, which we call the LOO-tCRPS criterion:

$$J_n(R) = \frac{1}{n} \sum_{i=1}^{n} S_{Q}^{\text{tCRPS}}(P_{R,n,-i}, f(x_i)). \qquad (26)$$

Using (26), we suggest the following procedure to select a reGP model. First, choose a sequence of nested candidate relaxation ranges $R^{(0)} \supset R^{(1)} \supset \dots \supset R^{(G-1)} = \emptyset$. The next step is the computation of $J_n(R^{(g)})$, $g = 0, \dots, G-1$, which involves the predictive distributions $P_{R^{(g)},n,-i}$.

In principle, (14) should be solved again each time a data point $(x_i, z_i)$ is removed, to obtain a pair $(\hat{\theta}_{n,-i}^{(g)}, \underline{z}_{n,-i}^{(g)})$ and then the corresponding reGP distribution $P_{R^{(g)},n,-i}$. To alleviate computational cost, a simple idea is to rely on the fast leave-one-out formulas (Dubrule, 1983) for Gaussian processes: for each set $R^{(g)}$, solve (14) to obtain $\hat{\theta}_n^{(g)}$ and $\underline{z}_n^{(g)} = (z_1^{(g)}, \dots, z_n^{(g)})^\mathrm{T}$, and then compute the
<figure><img src="image_7.png" /><figcaption>Figure 6: Left: Histogram of the values of the function $c_6$ from the G10 problem. Right: Same illustration but for the function $c_6^\alpha$, with $\alpha = 7$. The histograms are obtained from the values of the functions on a space-filling design of size $n = 100$. On the left, the values are very separated and concentrated on two modes, yielding a function close to a piecewise constant function. After transformation, the phenomenon is mitigated.</figcaption></figure>

## 4.4 An Example for the Estimation of an Excursion Set

We illustrate the method on the problem of estimating an excursion set {$x \in X, f(x) \le 0$}. We consider the G10 optimization problem used by Regis (2014), and focus on the constraint $c_6 \le 0$. Finding solutions satisfying the $c_6 \le 0$ constraint using a GP model is difficult, probably because the values of $c_6$ are very bi-modal, as illustrated in Figure 6. However, Feliot et al. (2017) found that the difficulty could be overcome by performing an ad-hoc monotonic transformation $z \mapsto z^\alpha$, with $\alpha = 7$, on the constraint.

The estimation of an excursion set {$f \le 0$} involves capturing precisely the behavior of $f$ around zero. Thus, we define a range of interest $Q = (-t^{(0)}, t^{(0)})$ centered on zero, with $t^{(0)}$ sufficiently small (note that there may be no data in $Q$). Then, we consider relaxation range candidates $R^{(g)} = (-\infty, -t^{(g)}) \cup [t^{(g)}, +\infty)$ with a sequence of thresholds $t^{(0)} < \dots < t^{(G)} = +\infty$, and we select $t^{(g)}$ by minimizing the LOO-tCRPS as described in the previous section.

In the case of the $c_6$ constraint and a small value for $t^{(0)}$ such that there is no data in $Q$, the results are presented in Figure 7. The LOO-tCRPS chooses $t^{(g)} = t^{(0)}$, so the reGP predictive distributions use only the information of being above or below zero. Moreover, observe that the corresponding transformation after relaxation bears resemblance to the transformation $z \mapsto z^\alpha$ proposed by Feliot et al. (2017). If we apply the reGP framework on the transformed function $c_6^\alpha$ (details omitted for brevity), we find that the LOO-tCRPS chooses a large $t^{(g)}$ such that the interpolation constraints are relaxed for a few observations only.

## 5. Application to Bayesian Optimization
<figure><img src="image_6.png" /><figcaption>Figure 5: Illustration of the choice of a relaxation range. The range of interest $\mathcal{Q}$ is determined by the threshold $t^{(0)}$. The relaxation range $\mathcal{R}$ corresponding to the region above $t$ has been obtained by the procedure described in Section 4.3.</figcaption></figure>

conditional distributions $\xi(x_i) | \{\xi(x_j) = z_j^{(g)}, j \neq i\}$, where $\xi \sim \text{GP}(\mu, k)$, and where $\mu$ and $k$ have parameter $\hat{\theta}_n^{(g)}$, using the fast leave-one-out formulas. By doing so, we neglect the difference between $\hat{\theta}_{n,-i}^{(g)}$ and $\hat{\theta}_n^{(g)}$ and the difference between $z_{n,-i}^{(g)}$ and the vector $(z_1^{(g)}, \dots, z_{i-1}^{(g)}, z_{i+1}^{(g)}, \dots, z_n^{(g)})^T$. The procedure ends by choosing the relaxation range $\mathcal{R}^{(g)}$ that achieves the best LOO-tCRPS value.

Figure 5 illustrates the selection of the relaxation range used in Figure 3. Algorithm 1 summarizes the general reGP procedure using LOO-tCRPS for selecting relaxation range in a given list of candidates for relaxation range. Section 5.3.1 and Appendix F give specific implementations for Bayesian optimization and the problem of estimation of an excursion set.

**Algorithm 1** reGP with automatic selection of the relaxation range.

**Input:** Data $(\underline{x}_n, \underline{z}_n)$; a range of interest $\mathcal{Q}$; and a list $\mathbb{R} \setminus \mathcal{Q} = \mathcal{R}^{(0)} \supset \dots \supset \mathcal{R}^{(G-1)} = \emptyset$ of relaxation range candidates.

**for** $g = 0$ **to** $G-1$ **do**

Obtain $\hat{\theta}_n^{(g)}$ and $\underline{z}_n^{(g)}$ by solving (14) with $\mathcal{R}^{(g)}$

Compute $J_n(\mathcal{R}^{(g)})$ with $\mathcal{Q}, \hat{\theta}_n^{(g)}$, and $\underline{z}_n^{(g)}$ using (26)

**end for**

**Output:** The pair $(\hat{\theta}_n^{(g)}, \underline{z}_n^{(g)})$ that minimizes (26).
<figure><img src="image_8.png" /><figcaption>Figure 7: A reGP fit of $c_6$, where the relaxation thresholds have been selected by LOO-tCRPS. The observations $\tilde{z}_n$ are shown on the x-axis, whereas the “relaxed” observations $z_n^*$ are represented on the y-axis. Moreover, the green lines represent the value zero, and the brown lines represent $\pm t^{(g)}$, with $t^{(g)}$ chosen to be $t^{(0)}$ by the LOO-tCRPS. Finally, the blue line shows a best fit by $z \mapsto z^7$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Left: Histogram of the values of the function $c_6$ from the G10 problem. Right: Same illustration but for the function $c_6^\alpha$, with $\alpha = 7$. The histograms are obtained from the values of the functions on a space-filling design of size $n = 100$. On the left, the values are very separated and concentrated on two modes, yielding a function close to a piecewise constant function. After transformation, the phenomenon is mitigated.</figcaption></figure>

4.4 An Example for the Estimation of an Excursion Set

We illustrate the method on the problem of estimating an excursion set {$x \in \mathcal{X}, f(x) \le 0$}. We consider the G10 optimization problem used by Regis (2014), and focus on the constraint $c_6 \le 0$. Finding solutions satisfying the $c_6 \le 0$ constraint using a GP model is difficult, probably because the values of $c_6$ are very bi-modal, as illustrated in Figure 6. However, Feliot et al. (2017) found that the difficulty could be overcome by performing an ad-hoc monotonic transformation $z \mapsto z^\alpha$, with $\alpha = 7$, on the constraint.

The estimation of an excursion set {$f \le 0$} involves capturing precisely the behavior of $f$ around zero. Thus, we define a range of interest $Q = (-t^{(0)}, t^{(0)})$ centered on zero, with $t^{(0)}$ sufficiently small (note that there may be no data in $Q$). Then, we consider relaxation range candidates $R^{(g)} = (-\infty, -t^{(g)}) \cup [t^{(g)}, +\infty)$ with a sequence of thresholds $t^{(0)} < \dots < t^{(G-1)} = +\infty$, and we select $t^{(g)}$ by minimizing the LOO-tCRPS as described in the previous section.

Appendix F presents numerical experiments on active learning of the set {$x \in \mathcal{X}, c_6(x) \le 0$} with reGP. Figure 7 illustrates the relaxation obtained using a value of $t^{(0)}$ such that $Q$ contains 25% of the observations of a design of size $n = 300$ in the experiments presented in Appendix F. Here, the LOO-tCRPS chooses $t^{(g)} = t^{(0)}$. Observe that the transformation after relaxation resembles the transformation $z \mapsto z^\alpha$ proposed by Feliot et al. (2017). Applying the reGP framework on the transformed function $c_6^\alpha$ (details omitted for brevity), we find that the LOO-tCRPS chooses a large $t^{(g)}$ such that the interpolation constraints are relaxed for only a few observations.
## 5.1 Efficient Global Optimization with Relaxation

The first motivation for introducing reGP models is Bayesian optimization, where obtaining good predictive distributions over ranges corresponding to optimal values is a key issue. In this article, we focus more specifically on the minimization problem

$$ \min_{x \in \mathcal{X}} f(x), \qquad (35) $$

where f is a real-valued function defined on a compact set $\mathcal{X} \subset \mathbb{R}^d$, but the methodology can be generalized to constrained and/or multi-objective formulations.

Given f, our objective is to construct a sequence of evaluation points $X_1, X_2, \dots \in \mathcal{X}$ by choosing each point $X_{n+1}$ as the maximizer of the expected improvement criterion (6) computed with respect to the reGP predictive distribution, with a relaxation set $R_n = [t_n, +\infty)$. More precisely, the sequence $(X_n)$ is constructed sequentially using the rule

$$ X_{n+1} = \arg\max_{x \in \mathcal{X}} E_n((m_n - \xi(x))_+), \qquad (36) $$

where $m_n = f(X_1) \wedge \dots \wedge f(X_n)$, and $E_n$ is the expectation under the reGP predictive distribution with relaxation set $R_n$ and data $z_n = (f(X_1), \dots, f(X_n))^T$.

As in Section 4.3, the relaxation threshold $t_n$ at iteration $n$ is chosen using the LOO-tCRPS criterion (34) among candidate values

$$ t_n^{(0)} < t_n^{(1)} < \dots < t_n^{(G)}, \qquad (37) $$

where $t_n^{(0)}$ is the validation threshold, which delimits the range of interest $Q_n = (-\infty, t_n^{(0)})$ used at iteration $n$. In the following, the optimization method just described will be called efficient global optimization with relaxation (EGO-R), in reference to the EGO name proposed by Jones et al. (1998).

Implementation specifics are detailed in Section 5.3. In the next section, we show that using the EI criterion with a reGP model yields a convergent algorithm.

## 5.2 Convergence of EGO-R with Fixed Parameters and Varying Threshold

In this section, we extend the result of Vazquez and Bect (2010b) and show the convergence of the EGO-R algorithm, in the case where the predictive distributions derive from a zero-mean Gaussian process with fixed covariance function.

We suppose that $\mathcal{X}$ is a compact domain and that $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is continuous, strictly positive-definite, and has the NEB (no-empty ball) property (Vazquez and Bect, 2010b), which says that the posterior variance cannot go to zero at a given point if there is no evaluation points in a ball centered on this point. In other words, the NEB property requires that the posterior variance $\sigma_n^2(x)$ at $x \in \mathcal{X}$ remains bounded away from zero for any $x$ not in the closure of the sequence of points $(X_n)$ evaluated by the optimization algorithm. A stationary covariance function with smoothness $\alpha > 0$ verifies the NEB property (Vazquez and Bect, 2010b), whereas the squared-exponential covariance function does not (Vazquez and Bect, 2010a).
<figure><img src="image_8.png" /><figcaption>Figure 7: A reGP fit of $c_6$, where the relaxation thresholds have been selected by LOO-tCRPS. The observations $\hat{z}_n$ are shown on the x-axis, whereas the “relaxed” observations $z_n^*$ are represented on the y-axis. The black squares represent interpolated observations and the blue points show relaxed observations. Moreover, the green lines represent the value zero, and the black lines represent $\pm t^{(g)}$, with $t^{(g)}$ chosen to be $t^{(0)}$ by the LOO-tCRPS.</figcaption></figure>
**Proposition 17** Let $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ be a continuous strictly positive-definite covariance function that verifies the NEB property, $\mathcal{H}(\mathcal{X})$ the corresponding RKHS and $f \in \mathcal{H}(\mathcal{X})$. Let $n_0 > 0$. Let $(X_n)_{n \ge 1}$ be a sequence in $\mathcal{X}$ such that, for each $n \ge n_0$, $X_{n+1}$ is obtained by (36) with $t_n > m_n$. Then the sequence $(X_n)_{n \ge 1}$ is dense in $\mathcal{X}$.

Proposition 17 implies the convergence of EGO-R with a fixed threshold $t > \min_{i \le n_0} f(X_i)$. In this case, the theoretical insights from Section 3.3 suggest a faster convergence might be achieved due to the improved error estimates (24) and (25) in a neighborhood of the global minimum.

The convergence of EGO-R also holds in the case of a varying relaxation set $R_n = [t_n, +\infty)$, with $t_n > m_n$, and in particular when $t_n$ is selected at each step using the LOO-tCRPS criterion (34) with a validation threshold $t_n^{(0)} > m_n$. In this case, the norm term in (24) gets smaller if $(t_n)_{n \ge 1}$ is decreasing.

## 5.3 Optimization Benchmark

In this section, we run numerical experiments to demonstrate the interest of using EGO-R instead of EGO for minimization problems.

### 5.3.1 METHODOLOGY

In practice, we must choose the sequence of thresholds (37). The validation threshold $t_n^{(0)}$ should be set above $m_n$ to ensure there is enough data to carry out the validation. We propose two different heuristics: a) a *constant heuristic*, where $t_n^{(0)}$ is kept constant through the iterations and set to an empirical quantile of an initial data set constructed before EGO-R is run, and b) a *concentration heuristic*, where $t_n^{(0)}$ corresponds to an empirical quantile of $z_n$.

In the case of the constant heuristic, we set $t_n^{(0)}$ to the $\alpha$-quantile of the function values on an initial design, which is typically built to fill $\mathcal{X}$ as evenly as possible with, e.g., maximin Latin hypercube sampling (McKay et al., 2000). The numerical experiments were conducted with $\alpha = 0.25$ in this article.

In the case of the concentration heuristic, we consider an $\alpha$-quantile of the values of $f$ at the points visited by the algorithm (again with $\alpha = 0.25$). As the optimization algorithm makes progress, the evaluations will likely concentrate around the global minimum. Thus, $t_n^{(0)}$ will get closer to the minimum value, and the range $Q_n = (-\infty, t_n^{(0)})$ of validation values will get smaller. Besides, since we expect better predictive distributions in this range, a better convergence may be obtained.

Both heuristics can be justified by the idealized setting of the convergence result from the previous section. Proposing alternative adaptive strategies to the concentration heuristic, or more generally conducting a theoretical study on the performance of such adaptive strategies, is out of the scope of this article.

For a given $t_n^{(0)}$, the candidate relaxation thresholds $t_n^{(g)}$, $g = 1, \dots, G$, are chosen so that $t_n^{(g)} - m_n$ ranges logarithmically from $t_n^{(0)} - m_n$ to $\max f(X_i) - m_n$ (with $G = 10$, in the experiments below).

To assess the performances of EGO-R with the two heuristics for choosing $t_n^{(0)}$, we compare them to the standard EGO algorithm. For all three algorithms, we use a first initial design of size $n_0 = 3d$, and we consider GPs with constant mean and a Matérn covariance function with regularity
# 5 Application to Bayesian Optimization

## 5.1 Efficient Global Optimization with Relaxation

The first motivation for introducing reGP models is Bayesian optimization, where obtaining good predictive distributions over ranges corresponding to optimal values is a key issue. In this article, we focus more specifically on the minimization problem

$$ \min_{x \in \mathcal{X}} f(x), $$

where $f$ is a real-valued function defined on a compact set $\mathcal{X} \subset \mathbb{R}^d$, but the methodology can be generalized to constrained and/or multi-objective formulations.

Given $f$, our objective is to construct a sequence of evaluation points $x_1, x_2, \dots \in \mathcal{X}$ by choosing each point $x_{n+1}$ as the maximizer of the expected improvement criterion (5) computed with respect to the reGP predictive distribution, with a relaxation range $R_n = [t_n, +\infty)$. More precisely, the sequence $(x_n)$ is constructed sequentially using the rule

$$ x_{n+1} = \arg\max_{x \in \mathcal{X}} E_n((m_n - \xi(x))_+), \qquad (27) $$

where $m_n = f(x_1) \wedge \dots \wedge f(x_n)$, and $E_n$ is the expectation under the reGP predictive distribution with relaxation range $R_n$ and data $\underline{z}_n = (f(x_1), \dots, f(x_n))^{\top}$.

As in Section 4.3, the relaxation threshold $t_n$ at iteration $n$ is chosen using the LOO-tCRPS criterion (26) among candidate values

$$ t_n^{(0)} < t_n^{(1)} < \dots < t_n^{(G-1)}, \qquad (28) $$

where $t_n^{(0)}$ is the validation threshold, which delimits the range of interest $Q_n = (-\infty, t_n^{(0)})$ used at iteration $n$. In the following, the optimization method just described will be called efficient global optimization with relaxation (EGO-R), in reference to the EGO name proposed by Jones et al. (1998).

Implementation specifics are detailed in Section 5.3. In the next section, we show that using the EI criterion with a reGP model yields a convergent algorithm.

## 5.2 Convergence of EGO-R with Fixed Parameters and Varying Threshold

In this section, we extend the result of Vazquez and Bect (2010b) and show the convergence of the EGO-R algorithm, in the case where the predictive distributions derive from a zero-mean Gaussian process with fixed covariance function.

We suppose that $\mathcal{X} \subset \mathbb{R}^d$ is a compact domain and that $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ is continuous, strictly positive-definite, and has the NEB (no-empty ball) property (Vazquez and Bect, 2010b), which says that the posterior variance cannot go to zero at a given point if there is no evaluation points in a ball centered on this point. In other words, the NEB property requires that the posterior variance $\sigma_n^2(x)$ at $x \in \mathcal{X}$ remains bounded away from zero for any $x$ not in the closure of the sequence of points $(x_n)$ evaluated by the optimization algorithm. A stationary covariance function with smoothness $\alpha > 0$ satisfies the NEB property (Vazquez and Bect, 2010b), whereas the squared-exponential covariance function does not (Vazquez and Bect, 2010a).
<table><caption>Table 1: Optimization benchmark. The acronym G-P stands for Goldstein-Price.</caption><thead><tr><td>Problem</td><td><i>d</i></td></tr></thead><tbody><tr><td>Branin</td><td>2</td></tr><tr><td>Six-hump Camel</td><td>2</td></tr><tr><td>Three-hump Camel</td><td>2</td></tr><tr><td>Hartman</td><td>3, 6</td></tr><tr><td>Ackley</td><td>4, 6, 10</td></tr><tr><td>Rosenbrock</td><td>4, 6, 10</td></tr><tr><td>Shekel</td><td>5, 7, 10</td></tr><tr><td>Goldstein-Price</td><td>2</td></tr><tr><td>Log-Goldstein-Price</td><td>2</td></tr><tr><td>Cross-in-Tray</td><td>2</td></tr><tr><td>Beale</td><td>2</td></tr><tr><td>Dixon-Price</td><td>4, 6, 10</td></tr><tr><td>Perm</td><td>4, 6, 10</td></tr><tr><td>Michalewicz</td><td>4, 6, 10</td></tr><tr><td>Zakharov</td><td>4, 6, 10</td></tr></tbody></table>

$v = 5/2$. The maximization of the sampling criteria (6) and (36) is carried out using a sequential Monte Carlo approach (Benassi et al., 2012; Feliot et al., 2017).

For reference, we also run the Dual simulated Annealing algorithm (inspired by Xiang et al. (1997)) from SciPy (Virtanen et al., 2020), with the default settings and with a random initialization.

The optimization algorithms are tested against a benchmark of test functions from Surjanovic and Bingham (2013) summarized in Table 1, with $n_{rep} = 30$ (random) repetitions, and a budget of $n_{tot} = 300$ evaluations for each repetition. This benchmark is partly inspired by Jones et al. (1998) and Merrill et al. (2021). In particular, we also use a log-version of the Goldstein-Price function as Jones et al. (1998).

To evaluate the algorithms we use, for each test function, several targets defined as spatial quantiles of the function and estimated with a subset simulation algorithm (see, e.g., Au and Beck, 2001). Then, the performances of the algorithms are evaluated using the fractions of runs that manage to reach the targets and the average numbers of evaluations to reach the targets (with unsuccessful runs counted as $n_{tot}$).

### 5.3.2 FINDINGS

The full set of results is provided in Appendix C. In Figure 8, we present a representative subset of these results.

First, observe in Figure 8 that the EGO-R methods can be considerably helpful and can outperform EGO largely on functions that are difficult to model with stationary GPs, such as Goldstein-Price, Perm (10), and Beale.
**Proposition 18** Let $k: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ be a continuous strictly positive-definite covariance function that satisfies the NEB property, $\mathcal{H}(\mathcal{X})$ the corresponding RKHS and $f \in \mathcal{H}(\mathcal{X})$. Let $n_0 > 0$. Let $(x_n)_{n \ge 1}$ be a sequence in $\mathcal{X}$ such that, for each $n \ge n_0$, $x_{n+1}$ is obtained by (27) with $t_n > m_n$. Then the sequence $(x_n)_{n \ge 1}$ is dense in $\mathcal{X}$.

Proposition 18 implies the convergence of EGO-R with a fixed threshold $t > \min_{i \le n_0} f(x_i)$. In this case, the theoretical insights from Section 3.3 suggest a faster convergence might be achieved due to the improved error estimates (20) and (21) in a neighborhood of the global minimum.

The convergence of EGO-R also holds in the case of a varying relaxation range $R_n = [t_n, +\infty)$, with $t_n > m_n$, and in particular when $t_n$ is selected at each step using the LOO-tCRPS criterion (26) with a validation threshold $t_n^{(0)} > m_n$. In this case, the norm term in (20) gets smaller if $(t_n)_{n \ge 1}$ is decreasing.

## 5.3 Optimization Benchmark

In this section, we run numerical experiments to demonstrate the interest of using EGO-R instead of EGO for minimization problems. Appendix E presents a similar study based on the UCB sampling criterion (Cox and John, 1992; Srinivas et al., 2010). Appendix F considers the task of estimating an excursion set.

### 5.3.1 METHODOLOGY

In practice, we must choose the sequence of thresholds (28). The validation threshold $t_n^{(0)}$ should be set above $m_n$ to ensure there is enough data to carry out the validation. We propose three different heuristics to this end: a) a *constant heuristic*, where $t_n^{(0)}$ remains constant throughout the iterations and is set to an empirical quantile of an initial data set constructed prior to running EGO-R, b) a *concentration heuristic*, where $t_n^{(0)}$ corresponds to an empirical quantile of $z_n$, and c) a *spatial heuristic*, where $t_n^{(0)}$ corresponds to a spatial quantile estimate of the function $f$. By spatial quantile, we mean a quantile of $f(U)$ with $U$ uniformly distributed over the domain $\mathcal{X}$.

In the case of the constant heuristic, we set $t_n^{(0)}$ to the $\alpha$-quantile of the function values on an initial design, which is typically built to cover $\mathcal{X}$ as uniformly as possible using, for example, maximin Latin hypercube sampling (McKay et al., 2000). In this article, the numerical experiments were conducted with $\alpha = 0.25$.

For the spatial heuristic, the spatial quantile is estimated by building a one-nearest-neighbor regression model trained on the data $(x_n, z_n)$, and then taking an $\alpha$-quantile of its predictions on a uniform sample on $\mathcal{X}$. The numerical value $\alpha = 0.25$ is also used here. This can be viewed as a refinement of the constant heuristic.

For the concentration heuristic, we consider the $\alpha$-quantile of the values of $f$ at the points visited by the algorithm (again with $\alpha = 0.25$). As the optimization algorithm progresses, the evaluations are expected to cluster around the global minimum. Thus, $t_n^{(0)}$ will approach the minimum value, and the validation range $Q_n = (-\infty, t_n^{(0)})$ will shrink. As a result, since better predictive distributions are expected in this range, a better convergence may be achieved.

All three heuristics are supported by the idealized convergence result from the previous section. Proposing alternative heuristics, or conducting a theoretical comparison of their performances, is out of the scope of this article.
Observe also that the EGO-R methods have about the same performance as EGO on functions that are easy to model with stationary GPs. This is the case of the Log-Goldstein-Price and the Branin functions, for which the LOO-tCRPS criterion for choosing the relaxation set detects that the larger values help predict near the minimum, and that no relaxation is needed as a result.

Furthermore, it is instructive to compare the performances of the EGO-R algorithms on the Goldstein-Price function on the one hand, with the performance of the EGO algorithm on Log-Goldstein-Price function on the other hand. Using reGP modeling enables to perform as well as with the logarithmic transform, but in an automatic way. This is also illustrated by Figure 9, which shows that the (non-parametric) transform learned by reGP resembles a logarithmic transform.

Finally, observe that the constant heuristic performs as well as EGO on Ackley (10), whereas the concentration heuristic lags behind. A closer look at the results for this function shows that the concentration heuristic get sometimes stuck in a local minimum. We explain this by the fact that the reGP model with the concentration heuristic can become very predictive in a small region around the local minimum, and underestimate the function variations elsewhere (the variance of the predictive distributions above $t_n^{(0)}$ are too small, and the optimization algorithm does not sufficiently explore unknown regions). To this regard, the constant heuristic is probably more conservative. Overall, taking the results from Appendix C into account, the concentration heuristic appears to be more (resp. less) efficient than the constant heuristic when there are few (resp. many) local minima.

# 6. Conclusion

This article presents a new technique called reGP to build predictive distributions for a function observed on a sequence of points. This technique can be applied when a user wants good predictive distributions in a range of function values, for example below a given threshold, and accepts degraded predictions outside this range. The technique relies on Gaussian process interpolation, and operates by relaxing the interpolation constraints outside the range of interest. This goal-oriented technique is kept simple and cheap: there are no additional parameters to set compared to the standard Gaussian process framework. The user only needs to specify a range of function values where good predictions should be obtained. The relaxation range can be selected automatically, using a scoring rule adapted to reGP models.

Such goal-oriented models can then be used in Bayesian sequential search algorithms. Here we are interested in the problem of mono-objective optimization and we propose to study the EI / EGO algorithm with such models. In a first step, we guarantee the convergence of the reGP-based algorithm on the RKHS attached to the underlying GP covariance. Then, we provide a benchmark that shows very clear benefits of using reGP models for the optimization of various functions.

A key element of the reGP approach is the definition of a range of interest $Q$, for instance a range of the form $(-\infty, t^{(0)})$ in a minimization problem. In some use cases the range will be provided by the user, but in others it is desirable to set it automatically. Two simple heuristics have been proposed in Section 5.3 to achieve this goal in our optimization benchmark, and it has been observed that the choice of heuristic has an impact on the exploratory behaviour of the resulting Bayesian optimization algorithm. Finding better heuristics, studying their properties, and assessing their impact in Bayesian optimization applications, is an important direction for future research.
<table><caption>Table 1: Optimization benchmark.</caption><thead><tr><td>Problem</td><td>d</td></tr></thead><tbody><tr><td>Branin</td><td>2</td></tr><tr><td>Six-hump Camel</td><td>2</td></tr><tr><td>Three-hump Camel</td><td>2</td></tr><tr><td>Hartman</td><td>3, 6</td></tr><tr><td>Ackley</td><td>4, 6, 10</td></tr><tr><td>Rosenbrock</td><td>4, 6, 10</td></tr><tr><td>Shekel</td><td>5, 7, 10</td></tr><tr><td>Goldstein-Price</td><td>2</td></tr><tr><td>Log-Goldstein-Price</td><td>2</td></tr><tr><td>Cross-in-Tray</td><td>2</td></tr><tr><td>Beale</td><td>2</td></tr><tr><td>Dixon-Price</td><td>4, 6, 10</td></tr><tr><td>Perm</td><td>4, 6, 10</td></tr><tr><td>Michalewicz</td><td>4, 6, 10</td></tr><tr><td>Zakharov</td><td>4, 6, 10</td></tr></tbody></table>

For a given $t_n^{(0)}$, the candidate relaxation thresholds $t_n^{(g)}$, $g = 0, \dots, G-1$, are chosen such that $t_n^{(g)} - m_n$ ranges logarithmically from $t_n^{(0)} - m_n$ to $\max f(x_i) - m_n$ (with $G = 10$ in the experiments below). We adopt the convention that $t_n^{(G-1)} > \max f(x_i)$, so that the G-th model is a standard GP.

To assess the performances of EGO-R with these three heuristics, we compare them to the standard EGO algorithm. For all four algorithms, we use an initial design of size $n_0 = 10d$, and we consider GPs with a constant mean function and a Matérn covariance function with regularity $\nu = 5/2$. The maximization of the sampling criteria (5) and (27) is performed using a sequential Monte Carlo approach (Benassi et al., 2012; Feliot et al., 2017).

The optimization algorithms are tested against a benchmark of test functions from Surjanovic and Bingham (2013) summarized in Table 1, with $n_{rep} = 100$ (random) repetitions, and a budget of $n_{tot} = 300$ evaluations for each repetition. This benchmark is partly inspired by Jones et al. (1998) and Merrill et al. (2021). In particular, we also use a log-version of the Goldstein-Price function as Jones et al. (1998).

To evaluate the algorithms, we construct for each test function a list of target values, defined as spatial quantiles of the function (estimated with a subset simulation algorithm; see, e.g., Au and Beck, 2001). The probability levels associated with this list of spatial quantiles define a non-linear scale that we use instead of raw function values. We plot in this scale the 10%/50%/90% quantiles of the best function value found, as a function of $n$. We also represent, for each probability level, the fraction of runs that reach the corresponding target.

### 5.3.2 FINDINGS

The full set of results is provided in Appendix D. In Figure 8, we present a representative subset of these results.

First, observe in Figure 8 that the EGO-R methods outperform EGO significantly on functions that are difficult to model with stationary GPs, such as Goldstein-Price, Perm (10), and Beale.
<figure><img src="image_9.png" /><figcaption>Figure 8: For each case, the top plot is the average number of iterations to reach the spatial quantile, and the bottom plot is the proportion of successful runs. Both are represented versus the level of the spatial quantile. The gray dotted line stands for the Dual simulated Annealing algorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the “Constant” and “Concentration” heuristics. The results are shown until the success rates of the three previous methods fall below the 66% threshold materialized by the black dashed line.</figcaption></figure>
In these cases, relaxing the constraints for the highest observations proves beneficial. Notice also that, among the three EGO-R variants, the one based on the concentration heuristic converges faster in such situations. (Approximately half of the cases of our benchmark fall into this category; see Appendix D for details.)

In the other cases (Log-Goldstein-Price, Ackley (4) and Ackley (10) on Figure 8) the EGO-R methods perform similarly to EGO. This corresponds to functions where classical stationary GP modeling already provides good predictions. In these cases, the LOO-tCRPS criterion for selecting the relaxation range detects that the largest values also contribute to predicting near the minimum, resulting in little to no relaxation being necessary. The concentration heuristic sometimes provides a small advantage on such functions (see Log-Goldstein-Price and Ackley (4)), but occasionally degrades performances compared to EGO and the other EGO-R variants, as seen with the Ackley (10) function. A closer examination shows that the concentration heuristic sometimes gets temporarily trapped in a local minimum. We explain this by the fact that the reGP model with the concentration heuristic can become overly predictive in a small region around the local minimum, but underestimate the function variations elsewhere (the variance of the predictive distributions above $t_n^{(0)}$ is too low, and the optimization algorithm does not explore unknown regions sufficiently). In this sense, the constant and spatial heuristics can be considered more conservative than the concentration heuristic.

Finally, it is useful to compare the performance of the EGO-R algorithms on the Goldstein-Price function with the performance of the EGO algorithm on the Log-Goldstein-Price function. Using reGP modeling allows for performance comparable to that achieved with a logarithmic transformation, but in an automatic manner. This is illustrated in Figure 9, where the (non-parametric) transform learned by reGP resembles a logarithmic transformation.

### 5.3.3 RUNNING TIME OF REGP

Training a reGP model takes longer than training a standard GP for several reasons. First, the optimization problem in (14) is an extension of (4) and involves a higher dimension of the search space, as the number of relaxed observations increases, which likely increases the time required to solve it. Second, the automatic selection procedure described in Section 4 involves testing $G$ relaxation ranges, making reGP at least $G$ times more costly to compute.

For $G = 10$, Figure 10 provides a quantitative assessment of the increase in computational time as a function of the number of evaluations $n$, across four representative test cases. With automatic selection of the relaxation range, training a reGP model is, roughly speaking, between $G = 10$ and 1000 times more expensive than training a standard GP. This trend holds for the other test functions in Table 1 (figures omitted for brevity). Among the heuristics, the concentration heuristic is the most computationally expensive version of reGP in this benchmark. This is because the number of observations outside the range of interest—and thus the maximum number of relaxed observations—grows linearly with $n$. In contrast, for the spatial and constant heuristics, the computational overhead remains around a factor of 10, which suggests that the primary cost increase comes from having to test $G$ relaxation ranges.

Note that the procedure described in Section 4.3 involves an exhaustive search over the set of relaxation range candidates. Assessing the performance of reGP with a smaller set of candidates, or using adaptive strategies such as dichotomies, could help reduce computation time. Note also that a factor $G$ in computation time can be saved if the user already knows how to choose a relaxation
<figure><img src="image_10.png" /><figcaption>Figure 9: A reGP fit to the Goldstein-Price function with $n = 30$ points, with $R = [10^3, +\infty)$. The observations $z_n$ are shown on the x-axis, whereas the relaxed observations $z_n^*$ are represented on the y-axis.</figcaption></figure>

More generally, the goal-oriented approach proposed in this article is not limited to single-objective (Bayesian) optimization. The example of Section 4.4 shows that it is also readily applicable, for instance, to level set estimation problems, for which a number of GP-based sequential design—aka active learning—strategies have been proposed in the literature (see, e.g., Chevalier et al., 2014; Bogunovic et al., 2016, and references therein). Other extensions are possible but will require more work. Noisy observations are one such example: considering that the interpolation constraints are already relaxed by the presence of noise, how should we transpose the goal-oriented modeling approach to this setting? Constrained and/or multi-objective optimization is another interesting but challenging direction for future research: in this case the function of interest in multivariate—one objective and several constraints, or several objectives—which requires significant adaptations to proposed methodology.
Goldstein-Price

<figure><img src="image_9.png" /></figure>

Log-Goldstein-Price

<figure><img src="image_10.png" /></figure>

Perm (10)

<figure><img src="image_11.png" /></figure>

Ackley (4)

<figure><img src="image_12.png" /></figure>

Beale

<figure><img src="image_13.png" /></figure>

Ackley (10)

<figure><img src="image_14.png" /><figcaption>Figure 8: For each case, the left plot shows the evolution of quantiles of the best function value against *n*, and the right plot represents on the x-axis the fraction of runs that reach each target value. Both plots use the probability level of the target as a scale for the y-axis. Left plot: solid lines represent the median; shaded areas are delimitated by the 10% and 90% quantiles. Both plots: EGO (red), EGO-R + constant heuristic (blue), EGO-R + concentration heuristic (green), EGO-R + spatial heuristic (black).</figcaption></figure>
## Appendix A. Properties of the Truncated CRPS

We shall now write (32) more explicitly for the case where the range of interest is an interval $Q = (a, b)$, $-\infty \le a < b \le +\infty$, and provide closed-form expressions for the case where, in addition, the predictive distribution $P$ is Gaussian.

**Remark 18** The value of the tCRPS for an interval $Q = (a, b)$ remains unchanged if the interval is closed at one or both of its endpoints.

**Remark 19** The value of the tCRPS for a finite (or countable) union of disjoint intervals follows readily from its values on intervals, since $Q \mapsto S_Q^{\text{tCRPS}}(P, z)$ is $\sigma$-additive.

We shall start by defining a quantity that shares similarities with (6).

**Definition 20** $\text{EI}_q^\uparrow(P, z) = E((N_1 \vee \cdots \vee N_q - z)_+)$ with $N_j \stackrel{\text{iid}}{\sim} P$.

The following expressions hold for a general predictive distribution $P$.

**Proposition 21** Suppose that $P$ has a first order moment.

* Let $a, b \in \mathbb{R}$ with $a \le b$. Then,

$$ S_{a,b}^{\text{tCRPS}}(P, z) = (b \wedge z - a)_+ + \text{EI}_2^\uparrow(P, b) - \text{EI}_2^\uparrow(P, a) - 2 \mathbb{1}_{z \le b} (\text{EI}_1^\uparrow(P, b) - \text{EI}_1^\uparrow(P, a \vee z)). \quad (38) $$

* Let $b \in \mathbb{R}$ and $N_1, N_2 \stackrel{\text{iid}}{\sim} P$. Then,

$$ S_{-\infty, b}^{\text{tCRPS}}(P, z) = b \wedge z + \text{EI}_2^\uparrow(P, b) - E(N_1 \vee N_2) - 2 \mathbb{1}_{z \le b} (\text{EI}_1^\uparrow(P, b) - \text{EI}_1^\uparrow(P, z)). \quad (39) $$

* Finally, if $a \in \mathbb{R}$, then

$$ S_{a, +\infty}^{\text{tCRPS}}(P, z) = S_{-\infty, -a}(P, z), \quad (40) $$

where $P$ is the distribution of $-U$ if $U$ is $P$-distributed.

Now, leveraging well-known analytic expressions (see, e.g., Nadarajah and Kotz, 2008; Chevalier and Ginsbourger, 2013), we have the following closed-form expressions in the Gaussian case.

**Proposition 22** (Nadarajah and Kotz, 2008; Chevalier and Ginsbourger, 2013)

Suppose that $P = \mathcal{N}(\mu, \sigma^2)$ and let $\phi$ and $\Phi$ denote respectively the pdf and the cdf of the standard Gaussian distribution. Then

$$ \text{EI}_1^\uparrow(P, z) = \sigma h_1\left(\frac{\mu-z}{\sigma}\right), \quad \text{with} \qquad (41) $$
<figure><img src="image_15.png" /><figcaption>Figure 9: A reGP fit to the Goldstein-Price function with $n = 30$ points, with $R = [10^3, +\infty)$. The observations $z_n$ are shown on the x-axis, whereas the relaxed observations $z_n^*$ are represented on the y-axis.</figcaption></figure>

range for his task. Another approach, left for future work, would be to solve (14) with dedicated algorithms (see Remark 8).

Finally, it should be stressed that the computation time is less of a concern for active learning tasks, since evaluating the objective function—which can take several hours—is usually the time bottleneck in this setting.

# 6 Conclusion

This article presents a new technique called reGP to build predictive distributions for a function observed on a sequence of points. This technique can be applied when a user wants good predictive distributions in a range of function values, for example below a given threshold, and accepts degraded predictions outside this range. The technique relies on Gaussian process interpolation, and operates by relaxing interpolation constraints outside the range of interest. The relaxation range can be selected automatically, using a scoring rule adapted to reGP models. With respect to classical GP modeling, this goal-oriented technique simply involves the choice of one additional key parameter, which is the range of interest—i.e., the user only needs to specify a range of function values where good predictions should be obtained.

Such goal-oriented models can then be used in Bayesian sequential search algorithms. Here we are more specifically interested in the problem of mono-objective optimization and we propose to study the EI / EGO algorithm with such models. We provide theoretical guarantees for the convergence of the reGP-based algorithm on the RKHS attached to the underlying GP covariance,
• for $q \ge 2$, we have $EI_q^\dagger(P, z) = \sigma h_q(\frac{\mu-z}{\sigma})$, where

$$h_q(t) = qt \Phi_q(\delta_q^t; 0, D_q D_q^T) + q \Phi(-t)^{q-1} \phi(t) + \frac{q(q-1)}{2\sqrt{\pi}} \Phi_{q-1}(\delta_{q-1}^t; 0, \frac{1}{2} B_q), \quad (42)$$

where $\Phi_q(\cdot; m, \Sigma)$ is the cdf of the multivariate $N(m, \Sigma)$ distribution,

$$B_q = 2 \operatorname{diag}(0, \mathbb{1}_{q-2}^T) + \mathbb{1}_{q-1} \mathbb{1}_{q-1}^T,$$

$D_q$ is the matrix representing the linear map

$$\mathbb{R}^q \to \mathbb{R}^q, \quad (y_1, \dots, y_q)^T \mapsto (-y_1, y_2 - y_1, y_3 - y_1, \dots, y_q - y_1)^T,$$

and $\delta_q^t = (t, 0_{q-1}^T)^T$,

• finally for $q = 2$ we have

$$E(N_1 \vee N_2) = \mu + \frac{\sigma}{\sqrt{\pi}}. \quad (43)$$

For a scoring rule $S: \mathcal{P} \times \mathbb{R} \to \mathbb{R}$ and $P_1, P_2 \in \mathcal{P}$ such that $y \in \mathbb{R} \mapsto S(P_1, y)$ is $P_2$-integrable, write $S(P_1, P_2) = E_{U \sim P_2}(S(P_1, U))$. The propriety of scoring rules is an important notion that formalizes “well-calibration” in the sense that a generating distribution must be identified to be optimal on average.

**Definition 23** (see, e.g. Gneiting and Raftery, 2007) A scoring rule $S: \mathcal{P} \times \mathbb{R} \to \mathbb{R}$ is said to be (strictly) proper with respect to $\mathcal{P}$ if, for all $P_1, P_2 \in \mathcal{P}$, the mapping $y \in \mathbb{R} \mapsto S(P_1, y)$ is $P_2$-integrable and the mapping $P_1 \in \mathcal{P} \mapsto S(P_1, P_2)$ admits $P_2$ as a (unique) minimizer.

A strictly proper scoring rule $S$ on a class $\mathcal{P}$ induces a divergence

$$(P_1, P_2) \mapsto S(P_1, P_2) - S(P_2, P_2),$$

which is non-negative on $\mathcal{P} \times \mathcal{P}$, and vanishes if and only if $P_1 = P_2$. In the case of the truncated CRPS, simple calculations lead to (Matheson and Winkler, 1976):

$$S_Q^{\text{CRPS}}(P_1, P_2) - S_Q^{\text{CRPS}}(P_2, P_2) = \int_Q (F_{P_1}(u) - F_{P_2}(u))^2 du.$$

It follows that $S_Q^{\text{CRPS}}$ is proper for any measurable $Q \subset \mathbb{R}$, and is strictly proper with respect to the class of non-degenerate Gaussian measures on $\mathbb{R}$ as soon as $Q$ has non-empty interior.
<figure><img src="image_16.png" /><figcaption>Figure 10: Ratio between the time needed to train a reGP model (including the automatic selection of the relaxation range) and that required to train a GP during Bayesian optimization with the EI criterion, across four representative test cases, as a function of *n*. The curves show median ratios, while the shaded areas represent intervals bounded by the 10% and 90% quantiles. Blue, green, and black represent reGP using the constant, concentration, and spatial heuristics, respectively.</figcaption></figure>
## Appendix B. Proofs

**Lemma 24** (Aronszajn, 1950, Section 1.5) Let $k: X \times X \to \mathbb{R}$ be a positive-definite covariance function, $U \subset X$, and $\mathcal{H}(U)$ be the RKHS attached to the restriction of $k$ to $U \times U$. The RKHS $\mathcal{H}(U)$ is the space of restrictions of functions from $\mathcal{H}(X)$ and the norm of $g \in \mathcal{H}(U)$ is given by

$$ \inf_{h \in \mathcal{H}(X), h|_{U}=g} \|h\|_{\mathcal{H}(X)}. \qquad (44) $$

**Proof of Proposition 3.** First the existence and the uniqueness of the solution are given by the first statement of Proposition 11 (with $\mathcal{H}_{R,n} = \mathcal{H}_{R, \{x_1, \dots, x_n\}}$).

Furthermore let $\underline{z} \in \mathbb{R}$ and write $\underline{\alpha} = K_n^{-1} \underline{z}$, the reproducing property (7) gives

$$ \|s_{\underline{z}}\|^2_{\mathcal{H}(X)} = \underline{\alpha}^{\mathrm{T}} K_n \underline{\alpha} = \underline{z}^{\mathrm{T}} K_n^{-1} \underline{z}, \qquad (45) $$

and therefore

$$ \min_{h \in \mathcal{H}_{R,n}} \|h\|^2_{\mathcal{H}(X)} = \inf_{\underline{z} \in C_{R,n}} \min_{h \in \mathcal{H}(X), h(x_n)=\underline{z}} \|h\|^2_{\mathcal{H}(X)} = \inf_{\underline{z} \in C_{R,n}} \underline{z}^{\mathrm{T}} K_n^{-1} \underline{z}, $$

where the last infimum is uniquely reached by the evaluation of the solution on $x_n$.

**Proof of Proposition 4.** Write $K_{m,n}$ for the covariance matrix of the random vector $(\underline{z}'^\mathrm{T}, \underline{z}_n^\mathrm{T})^\mathrm{T}$. Using the equalities (5) and (45), and a slight abuse of notation by dropping irrelevant constants with respect to $\underline{z}'$ and $\underline{z} \in C_{R,n}$, we have

$$ -2 \ln \left( p \left( \underline{z}', \underline{z} \mid \underline{Z}_n \in C_{R,n} \right) \right) = \left( \underline{z}'^\mathrm{T}, \underline{z}^\mathrm{T} \right) K_{m,n}^{-1} \left( \underline{z}'^\mathrm{T}, \underline{z}^\mathrm{T} \right)^\mathrm{T} = \min_{h \in \mathcal{H}(X), h(x_n)=\underline{z}, h(x_m)=\underline{z}'} \|h\|^2_{\mathcal{H}(X)}. $$

This gives

$$ \inf_{\underline{z}' \in \mathbb{R}^m, \underline{z} \in C_{R,n}} -2 \ln \left( p \left( \underline{z}', \underline{z} \mid \underline{Z}_n \in C_{R,n} \right) \right) = \min_{h \in \mathcal{H}(X), h(x_n) \in C_{R,n}} \|h\|^2_{\mathcal{H}(X)}, $$

which is reached by taking $\underline{z} = \underline{z}_n^*$ and $\underline{z}' = (s_{\underline{z}_n^*}^{\mathrm{T}}(x_1'), \dots, s_{\underline{z}_n^*}^{\mathrm{T}}(x_m'))^\mathrm{T}$.

**Proof of Proposition 9.** First, one has

$$ \sup_{x \in X} \sigma_n(x) = \sup_{x \in X} \sup_{\|f\|_{\mathcal{H}(X)}=1} |f(x) - s_{\underline{z}_n}(x)| = \sup_{\|f\|_{\mathcal{H}(X)}=1} \|f - s_{\underline{z}_n}\|_{L^\infty(X)}. $$

Now, let $f \in \mathcal{H}(X)$ such that $\|f\|_{\mathcal{H}(X)} = 1$, and $X^o$ be the interior of $X$. The boundary of $X^o$ is the one of $X$ under Assumption 1, and the Sobolev space $W_2^{\alpha+d/2}(X^o)$ defined by (17) is norm-equivalent to the Sobolev–Slobodeckij space (see, e.g., Di Nezza et al. (2012, Proposition 3.4) for a statement on $\mathbb{R}^d$ and Grisvard (1985, Theorem 1.4.3.1) for the existence of an extension operator).

Then, one can apply (Arcangéli et al., 2007, Theorem 4.1) to $f - s_{\underline{z}_n}$ restricted to $X^o$ to show that, for $h_n$ lower than some $h_0$ (not depending on $f$ or $(x_n)_{n \ge 1}$), we have:

$$ \|f - s_{\underline{z}_n}\|_{L^\infty(X)} = \|f - s_{\underline{z}_n}\|_{L^\infty(X^o)} \lesssim h_n^\alpha \|f - s_{\underline{z}_n}\|_{W_2^\alpha(X^o)} \lesssim h_n^\alpha \|f - s_{\underline{z}_n}\|_{\mathcal{H}(X)} \lesssim h_n^\alpha \qquad (46) $$
and a numerical benchmark that shows very clear benefits of using reGP models for the optimization
of various functions.

A key element of the reGP approach is the definition of the range of interest, which is for
instance of the form $(-\infty, t^{(0)})$ in a minimization problem. In some use cases the range will be
provided by the user, but in others it is desirable to set it automatically. Three simple heuristics
have been proposed to achieve this goal in our optimization benchmark, and it has been observed
that the choice of heuristic has an impact on the exploratory behaviour of the resulting Bayesian
optimization algorithm. Finding better heuristics, studying their properties, and assessing their
impact in Bayesian optimization applications, is an important direction for future research.

More generally, the goal-oriented approach proposed in this article is not limited to single-
objective (Bayesian) optimization. The example from Section 4.4 and the benchmark results of
Appendix F show that it is also readily applicable, for instance, to level set estimation problems.
Appendix G provides a proof of concept for extending the approach to noisy observations. A thor-
ough investigation of the best modeling choices and theoretical guarantees is left for future work.
Other extensions are possible but will require more work. Constrained and/or multi-objective opti-
mization is another interesting but challenging direction for future research: in this case the func-
tion of interest is multivariate (one objective and several constraints, or several objectives), which
requires significant adaptations to the proposed methodology.

Another possible direction for future research would be to infer the unknown function using the
full predictive distribution (12), of which reGP is a mode approximation. This distribution defines
a Skew GP, for which various inference methods have been proposed in the literature, using for
instance quasi Monte-Carlo sampling, linear elliptical slice sampling or Gibbs sampling (see, e.g.,
Da Veiga and Marrel, 2012; Botev, 2017; Benavoli et al., 2021; Gessner et al., 2020; Takeno et al.,
2023). Adapting these ideas to the reGP setting may be relevant in a context where uncertainty
quantification is critical.

## Acknowledgments and Disclosure of Funding

The authors thank the action editor and three anonymous reviewers for their insightful comments, which significantly improved both the content and the presentation of the article. This work was partially supported by the French Agence Nationale de la Recherche et de la Technologie (ANRT) through a CIFRE grant.
by continuity of $f - s_{z_n}$, since $\|\cdot\|_{W_2^{\alpha+d/2}(X^o)} \le \|\cdot\|_{W_2^{\alpha+d/2}(X)}$ due to the definition (17), $W_2^{\alpha+d/2}(X)$ being norm equivalent to $\mathcal{H}(X)$, and because of the projection interpretation of $s_{z_n}$ (see, e.g., Wendland, 2004, Theorem 13.1). Finally, one can get rid of the condition $h_n \le h_0$ for simplicity by increasing the constant eventually, since $\sup_{x \in X} \sigma_n(x) \le \sup_{x \in X} \sqrt{k(x, x)} < +\infty$ by compacity. ■

**Proof of Proposition 11.** First observe that $\mathcal{H}_{R,U}$ is not empty since it contains $f$. Furthermore, it is easy to verify that $\mathcal{H}_{R,U}$ is convex and that it is closed since pointwise evaluation functionals are continuous on an RKHS. The problem is then the one of projecting the null function on a convex closed subset; hence the existence and the uniqueness.

Then, the function $s_{R,n}$ is the projection of the null function on the closed convex set $\mathcal{H}_{R,n}$. Moreover, the sequence $(\mathcal{H}_{R,n})_{n \ge 1}$ is non-increasing, so the sequence $(s_{R,n})_{n \ge 1}$ converges in $\mathcal{H}(X)$ to the projection of the null function on $\bigcap_{n \ge 1} \mathcal{H}_{R,n}$ (see, e.g., Brezis, 2011, Exercise 5.5), i.e. the solution of (22), with $U = \{x_n\}$. But this last solution is also the solution on the closure since it verifies the constraints by continuity. ■

**Proof of Proposition 12.** Define $x_n^0$ and $z_n^0$ to be data points within $X_0$, and $s_{x_n^0, z_n^0}$ the associated (interpolation) predictor, i.e. the solution of (8). Observing that $s_{x_n^0, z_n^0}$ interpolates $s_{R,n}$, we have:

$$ \begin{aligned} |f(x) - s_{R,n}(x)| &\le |f(x) - s_{x_n^0, z_n^0}(x)| + |s_{x_n^0, z_n^0}(x) - s_{R,n}(x)| \\ &\le \sigma_{n,0}(x) \|f\|_{\mathcal{H}(X_0)} + \sigma_{n,0}(x) \|s_{R,n}\|_{\mathcal{H}(X_0)} \\ &\le 2\sigma_{n,0}(x) \|s_{R,n}\|_{\mathcal{H}(X)}, \end{aligned} $$

since $f$ coincides with $s_{R,X}$ on $X_0$, $s_{R,X} \in \mathcal{H}_{R,n}$, and by Lemma 24. ■

**Lemma 25** If $k$ has smoothness $\alpha > 0$, then there exists $h_0 > 0$ depending only on $\alpha$ such that, for all $x, y \in X$ verifying $\|x - y\| \le h_0$, we have

$$ \|\delta_y - \delta_x\|_{\mathcal{H}^*(X)} \lesssim \|x - y\|^\alpha, \text{ for } \alpha < 1, \quad (47) $$

$$ \|\delta_y - \delta_x\|_{\mathcal{H}^*(X)} \lesssim \sqrt{|\ln(\|x - y\|)|} \|x - y\|, \text{ for } \alpha = 1, \quad (48) $$

and

$$ \|\delta_y - \delta_x\|_{\mathcal{H}^*(X)} \lesssim \|x - y\|, \text{ for } \alpha > 1. \quad (49) $$

**Proof.** Since equivalent norms give equivalent operator norms on the topological dual of a normed space, it suffices to show the result for a unit-variance isotropic Matérn covariance function (4) of regularity $\alpha$.

In this case, we have

$$ \|\delta_y - \delta_x\|^2_{\mathcal{H}^*(X)} = k(x, x) + k(y, y) - 2k(x, y) = 2(1 - r_\alpha(\|x - y\|)), $$

with $r_\alpha$ the corresponding isotropic correlation function. Standard results on principal irregular terms (see, e.g., Stein, 1999, Chapter 2.7) give the results. ■
**Appendix A. Properties of the Truncated CRPS**

We shall now write (25) more explicitly for the case where the range of interest is an interval $Q = (a, b)$, $-\infty \le a < b \le +\infty$, and provide closed-form expressions for the case where, in addition, the predictive distribution $P$ is Gaussian.

**Remark 19** The value of the tCRPS for an interval $Q = (a, b)$ remains unchanged if the interval is closed at one or both of its endpoints.

**Remark 20** The value of the tCRPS for a finite (or countable) union of disjoint intervals follows readily from its values on intervals, since $Q \mapsto S_Q^{\text{tCRPS}}(P, z)$ is $\sigma$-additive.

We shall start by defining a quantity that shares similarities with (5).

**Definition 21** $\text{EI}_q^{\uparrow}(P, z) = E((N_1 \vee \cdots \vee N_q - z)_+)$ with $N_j \stackrel{\text{iid}}{\sim} P$.

The following expressions hold for a general predictive distribution $P$.

**Proposition 22** Suppose that $P$ has a first order moment.

* Let $a, b \in \mathbb{R}$ with $a \le b$. Then,

$$S_{a,b}^{\text{tCRPS}}(P, z) = (b \wedge z - a)_+ + \text{EI}_2^{\uparrow}(P, b) - \text{EI}_2^{\uparrow}(P, a) - 2 \mathbb{1}_{z \le b} (\text{EI}_1^{\uparrow}(P, b) - \text{EI}_1^{\uparrow}(P, a \vee z)).$$

* Let $b \in \mathbb{R}$ and $N_1, N_2 \stackrel{\text{iid}}{\sim} P$. Then,

$$S_{-\infty, b}^{\text{tCRPS}}(P, z) = b \wedge z + \text{EI}_2^{\uparrow}(P, b) - E(N_1 \vee N_2) - 2 \mathbb{1}_{z \le b} (\text{EI}_1^{\uparrow}(P, b) - \text{EI}_1^{\uparrow}(P, z)).$$

* Finally, if $a \in \mathbb{R}$, then

$$S_{a, +\infty}^{\text{tCRPS}}(P, z) = S_{-\infty, -a}^{\text{tCRPS}}(P, -z),$$

where $P$ is the distribution of $-U$ if $U$ is $P$-distributed.

Now, leveraging well-known analytic expressions (see, e.g., Nadarajah and Kotz, 2008; Chevalier and Ginsbourger, 2013), we have the following closed-form expressions in the Gaussian case.

**Proposition 23** (Nadarajah and Kotz, 2008; Chevalier and Ginsbourger, 2013)

Suppose that $P = \mathcal{N}(\mu, \sigma^2)$ and let $\phi$ and $\Phi$ denote respectively the pdf and the cdf of the standard Gaussian distribution. Then

* $\text{EI}_1^{\uparrow}(P, z) = \sigma h_1(\frac{\mu-z}{\sigma})$, with $h_1(t) = t\Phi(t) + \phi(t)$,
**Lemma 26** Let $B \subset \mathbb{X}$ verifying Assumption 1 and $h_{n,B}$ be the fill distance of $X_{n,B} = \{x_1, \cdots, x_n\} \cap B$ within $B$, with the convention $h_{n,B} = \text{diam}(B)$ if $X_{n,B}$ is empty. Then, $h_{n,B} \lesssim h_n$.

**Proof.** The idea of the proof is given by Wendland (2004, Lemma 11.31), but it is interlinked with a much more sophisticated construction, so we provide a specific version here for completeness. Adams and Fournier (2003, Section 4.11) shows that $B$ verifies a cone condition with radius $\rho > 0$ and angle $\phi \in (0, \pi/2)$. If $X_{n,B}$ is not empty, then the compacity of $B$ ensures the existence of an $x \in B$ such that $d(x, X_{n,B}) = h_{n,B}$. (If $X_{n,B}$ is empty, then the rest of the proof is also valid taking an arbitrary $x \in B$.)

A cone $C$ originating from $x$ with angle $\phi$ and radius $\delta = \min(h_{n,B}, \rho)$ is contained in $B$ and its interior do not contains observations. Furthermore, Wendland (2004, Lemma 3.7) shows that there exists a $y \in C$ such that the open ball $B(y, \delta \sin(\phi)(1+\sin(\phi))^{-1})$ is subset of $C$, and therefore contains no observations as well. Thus, $h_n \ge \delta \sin(\phi)(1+\sin(\phi))^{-1}$. Now, if $h_{n,B} \le \rho$, then the desired result follows. If not, the result holds as well since $h_{n,B} \le \text{diam}(B)$. ■

**Proof of Proposition 13.** For the first assertion, let $\sigma_{n,B}$ be the power function using only the observations within $B$. Using Proposition 12, the inequality $\sigma_{n,0} \le \sigma_{n,B}$ given by the projection residuals interpretation, and applying Proposition 9 to $B$ yields a bound depending on the fill distance $h_{n,B}$ of $\{x_1, \cdots, x_n\} \cap B$ within $B$. Finally, Lemma 26 allows us to conclude.

Regarding second assertion, $f$ is continuous so the sets $X_j$ are compact for $j \ge 1$. In addition, they are disjoint so

$$ \delta = \min_{1 \le j < p} \inf_{x \in X_j, y \in X_p} \|x-y\| > 0. $$

Suppose now that $h_n < \delta$ and let $j \ge 1$, $x \in X_j$ and $1 \le i \le n$ the index of the closet $x_i$ to $x$. By definition, $\|x-x_i\| \le h_n$ and therefore $x \in X_0 \cup X_j$. Let $\mathcal{H}^*(X)$ be the (topological) dual of $\mathcal{H}(X)$ and $\delta_y: h \in \mathcal{H}(X) \mapsto h(y)$, which lies in $\mathcal{H}^*(X)$ for all $y \in X$. Then using the reproducing property (7), we have

$$ |s_{R,n}(x_i) - s_{R,n}(x)| \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(X)} \|s_{R,n}\|_{\mathcal{H}(X)} \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(X)} \|s_{R,X}\|_{\mathcal{H}(X)}, $$

and therefore

$$ d(s_{R,n}(x), R_j) \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(X)} \|s_{R,X}\|_{\mathcal{H}(X)} + d(s_{R,n}(x_i), R_j). $$

Now, if $x_i \in X_j$, then $d(s_{R,n}(x_i), R_j) = 0$. Otherwise, $x_i \in X_0$ necessarily and then, using the fact that $s_{R,X}(x) \in R_j$, we have:

$$ d(s_{R,n}(x_i), R_j) \le |s_{R,n}(x_i) - s_{R,X}(x)| = |s_{R,X}(x_i) - s_{R,X}(x)| \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(X)} \|s_{R,X}\|_{\mathcal{H}(X)}. $$

So one can use Lemma 25 along with the previous statements to conclude if $h_n < \min(\delta, h_0)$.

Finally, treating the case $h_n \ge \min(\delta, h_0)$ is straightforward using the fact that $\sup_{x \in X} \sqrt{k(x, x)}$ is finite thanks to the compacity of $X$ and $d(s_{R,n}(x), R_j) \le |s_{R,n}(x) - s_{R,X}(x)|$ for $j \ge 1$ and $x \in X_j$. ■

**Lemma 27** If $g, h \in W_2^\gamma(X)$ for $\gamma > d/2$, then $gh \in W_2^\gamma(X)$.
• for $q \ge 2$, we have $\mathrm{EI}_q^\dagger(P, z) = \sigma h_q(\frac{\mu-z}{\sigma})$, where

$$h_q(t) = qt \Phi_q(\delta_q^t; 0, D_q D_q^T) + q \Phi(-t)^{q-1} \phi(t) + \frac{q(q-1)}{2\sqrt{\pi}} \Phi_{q-1}(\delta_{q-1}^t; 0, \frac{1}{2} B_q),$$

where $\Phi_q(\cdot; m, \Sigma)$ is the cdf of the multivariate $N(m, \Sigma)$ distribution,

$$B_q = 2 \mathrm{diag}(0, \mathbb{1}_{q-2}^T) + \mathbb{1}_{q-1} \mathbb{1}_{q-1}^T,$$

$D_q$ is the matrix representing the linear map

$$\mathbb{R}^q \to \mathbb{R}^q, \quad (y_1, \dots, y_q)^T \mapsto (-y_1, y_2 - y_1, y_3 - y_1, \dots, y_q - y_1)^T,$$

and $\delta_q^t = (t, 0_{q-1}^T)^T$,

• finally for $q = 2$ we have

$$E(N_1 \vee N_2) = \mu + \frac{\sigma}{\sqrt{\pi}}.$$

The propriety of scoring rules is an important notion that formalizes “well-calibration” in the sense that a generating distribution must be identified to be optimal on average.

**Definition 24** (see, e.g., Gneiting and Raftery, 2007) A scoring rule $S: \mathcal{P} \times \mathbb{R} \to \mathbb{R}$ is said to be (strictly) proper with respect to $\mathcal{P}$ if, for all $P_1, P_2 \in \mathcal{P}$, the mapping $y \in \mathbb{R} \mapsto S(P_1, y)$ is $P_2$-quasi-integrable and the mapping

$$P_1 \in \mathcal{P} \mapsto S(P_1, P_2) := E_{U \sim P_2} (S(P_1, U))$$

admits $P_2$ as a (unique) minimizer.

In the case of the truncated CRPS, simple calculations lead to (Matheson and Winkler, 1976):

$$S_Q^{\mathrm{CRPS}}(P_1, P_2) = S_Q^{\mathrm{CRPS}}(P_2, P_2) + \int_Q (F_{P_1}(u) - F_{P_2}(u))^2 du.$$

It follows that $S_Q^{\mathrm{CRPS}}$ is proper with respect to the class of all Borel probability measures on $\mathbb{R}$ for any measurable $Q \subset \mathbb{R}$, and is strictly proper with respect to the class of non-degenerate Gaussian measures on $\mathbb{R}$ as soon as $Q$ has non-empty interior.
**Proof.** By the definition (17) of $W_2^\gamma(\mathbb{X})$, the functions $g$ and $h$ can be extended into functions on $\mathbb{R}^d$, having their product in $W_2^\gamma(\mathbb{R}^d)$ (Strichartz, 1967, Theorem 2.1). Taking the restrictions shows the desired result. ■

**Proof of Proposition 14.** We use a bump function argument. Let $B(x_0, r) \subset \mathbb{X}_j$ (with $r > 0$) be an open ball. There exists a $C^\infty$ function $\phi: \mathbb{R}^d \to \mathbb{R}$ such that

$$ \begin{cases} 0 \le \phi \le 1, \\ \phi(x) = 1 & \text{only if } x = 0, \\ \phi(x) = 0 & \text{if } x \in \mathbb{X} \setminus B(0, r). \end{cases} $$

Let $c \in R_j \setminus \{f(x_0)\}$, $\phi_n = \phi(n(\cdot - x_0))$ as a function on $\mathbb{X}$, and $f_n = (1 - \phi_n)f + c\phi_n$, for $n \ge 1$. We have $\phi_n \in W_2^{\alpha+d/2}(\mathbb{R}^d)$ as a function on $\mathbb{R}^d$, so it belongs to $W_2^{\alpha+d/2}(\mathbb{X})$ as a function on $\mathbb{X}$, and Lemma 27 ensures that $f_n \in \mathcal{H}(\mathbb{X})$. Moreover, it is easy to check that $f_n \in \mathcal{H}_{R, \mathbb{X}}$. Observe that the sequence $(f_n)_{n \ge 1}$ converges pointwise to a discontinuous function that lies thus outside $\mathcal{H}(\mathbb{X})$.

Suppose now that $\|f_n\|_{\mathcal{H}(\mathbb{X})} \to +\infty$ then one can extract a bounded subsequence of norms and a classical weak compacity argument would yield a weakly convergent subsequence, which is impossible since the pointwise limit is not in $\mathcal{H}(\mathbb{X})$. ■

**Lemma 28** Use the notations of Proposition 17 and let $(y_n)_{n \ge 1}$ be a sequence in $\mathbb{X}$. Assume that the sequence $(y_n)_{n \ge 1}$ is convergent, denote by $y^*$ its limit and assume that $y^*$ is an adherent point of the set $\{x_n\}$. Let $t_\infty = \liminf t_n$, then

* $s_n(y_n) \to f(y^*)$ if $f(y^*) < t_\infty$,

* $\liminf s_n(y_n) \ge t_\infty$, otherwise.

In particular, we have

$$ \liminf s_n(y_n) \ge \min(f(y^*), t_\infty). \qquad (50) $$

**Proof.** Suppose that $y^* \notin \{x_n\}$. Then, let $(x_{\phi(n)})_{n \ge 1}$ be a subsequence converging to $y^*$ and let $\psi(n) = \max\{\phi(k), \phi(k) \le n\}$. We proceed as in Proposition 13 to have:

$$ |s_n(x_{\psi(n)}) - s_n(y_n)| \le \|\delta_{x_{\psi(n)}} - \delta_{y_n}\|_{\mathcal{H}^*(\mathbb{X})} \|s_n\|_{\mathcal{H}(\mathbb{X})} \to 0, \qquad (51) $$

thanks to Lemma 25 and the inequality $\|s_n\|_{\mathcal{H}(\mathbb{X})} \le \|f\|_{\mathcal{H}(\mathbb{X})}$.

Finally, have $s_n(x_{\psi(n)}) \ge \min(f(x_{\psi(n)}), t_n)$ by construction, so $\liminf s_n(x_{\psi(n)}) \ge \min(f(y^*), t_\infty)$, which gives the second assertion thanks to (51). Observe that $f(x_{\psi(n)}) < t_n$ ultimately if $f(y^*) < t_\infty$ for the first assertion.

If $y^* \in \{x_n\}$, then the result also follows similarly. ■

**Lemma 29** Using the notations of Proposition 17 and writing $v_n = \sup_{x \in \mathbb{X}} \rho_{n,t_n}(x)$, we have $\liminf v_n = 0$.
# Appendix B. Theoretical Insights on reGP

This appendix provides further theoretical study of relaxed interpolation. We start by recalling a well-known approximation result for the interpolation of classes of functions outside the RKHS. Then, we give a corresponding adaptation for relaxed interpolation. Finally, Appendix B.3 presents results on approximation and uncertainty quantification for reGP with estimated parameters. The notations are that of Section 3.3.

## B.1 Escape Theorem for Gaussian Process Interpolation

In GP interpolation, error bounds are sometimes available for interpolation of functions not belonging to the RKHS. Let $q_n = \min_{i \neq j} \|x_i - x_j\|/2$ be the separation distance and $\rho_n = h_n/q_n$ be the mesh ratio. A sequence $(x_n)_{n \ge 1}$ is quasi-uniform if there exists a positive constant $C_{qu}$ such that $q_n \le h_n \le C_{qu}q_n$. The following proposition is a reminder of the well-known “escape” theorem, which provides error bounds for standard Gaussian process interpolation of classes of functions outside the RKHS attached to a kernel with regularity $\alpha > 0$ (see Theorem 4.2 of Narcowich et al., 2006, and Theorem 4.2 of Karvonen et al., 2020, and Theorem 1 of Wynne et al., 2021, for recent extensions).

**Proposition 25** Let $0 < \beta < \alpha < +\infty$ be such that $f \in W_2^{\beta+d/2}(\mathbb{X})$ and assume that $\mathbb{X}$ satisfies Assumption 9. If $k$ has regularity $\alpha$, then

$$ \|f - s_{z_n}\|_{L^\infty(\mathbb{X})} \lesssim h_n^\beta (1 + \rho_n^{\alpha-\beta}) \|f\|_{W_2^{\beta+d/2}(\mathbb{X})}, $$

where $\lesssim$ denotes inequality up to a constant, that does not depend on $f$, $n$ or $(x_n)$. In particular, if $(x_n)_{n \ge 1}$ is quasi-uniform, we have

$$ \|f - s_{z_n}\|_{L^\infty(\mathbb{X})} \lesssim n^{-\beta/d} \|f\|_{W_2^{\beta+d/2}(\mathbb{X})}. $$

## B.2 Convergence of Relaxed Interpolation of Functions Outside the RKHS

Section 3.3.2 only considers functions belonging to the RKHS. There are other situations in which reGP is expected to offer advantages. For the remainder of this appendix, assume $f: \mathbb{X} \to \mathbb{R}$ is an arbitrary function. For a relaxation range $R = \bigcup_{j=1}^q R_j$, the subsets $X_j \subset \mathbb{X}$ are defined as in Section 3.3.2.

**Definition 26** Let $\mathcal{F}$ be a set of functions $\mathbb{X} \to \mathbb{R}$. We say that $f$ is an $R$-element of $\mathcal{F}$ if there exists a function $g \in \mathcal{F}$ such that $g$ and $f$ coincide on $X_0$ and $f(x) \in R_j \Leftrightarrow g(x) \in R_j$, for all $j \ge 1$ and $x \in \mathbb{X}$.

Naturally, if $f \in \mathcal{F}$, then $f$ is an $R$-element of $\mathcal{F}$. Note also that $f$ can be an $R$-element of a space of continuous functions and be discontinuous in the interior of an $X_j$, for $j \ge 1$. When $\mathcal{F} = \mathcal{H}(\mathbb{X})$, this condition is equivalent to requiring that the set $\mathcal{H}_{R, \mathbb{X}}$, defined by (18), is non-empty.

Suppose $f \notin \mathcal{H}(\mathbb{X})$ but $f$ is an $R$-element of $\mathcal{H}(\mathbb{X})$. In this case, the error bounds (8) and (17) are invalid, whereas the conclusions of Propositions 12–14 still hold. (Indeed, observe that the reGP predictor of $f$ is that of a function $g$ given by Definition 26. The propositions apply to $g$.)

Moreover, even if $f$ is not an $R$-element of $\mathcal{H}(\mathbb{X})$, error bounds can sometimes be stated on $X_0$. Recall that the space $W_2^{\beta+d/2}(\mathbb{X})$ is an RKHS with a continuous kernel for $\beta > 0$. If $f$ is
**Proof.** This is an adaptation of Lemma 12 from (Vazquez and Bect, 2010b).

Let $y^*$ be a cluster point of $(x_n)_{n \ge 1}$ and let $(x_{\phi(n)})_{n \ge 1}$ be a subsequence converging to $y^*$. We are going to prove that $v_{\phi(n)-1} \to 0$. Define

$$z_{\phi(n)-1} = m_{\phi(n)-1} - s_{\phi(n)-1}(x_{\phi(n)}).$$

Then, Lemma 28 gives<fnref n="1" /> $\liminf_{n \to \infty} s_{\phi(n)-1}(x_{\phi(n)}) \ge \min(f(y^*), t_\infty)$. Moreover we have

$$m_{\phi(n)-1} \le \min(f(x_{\phi(n-1)}), t_{\phi(n-1)}),$$

since $\phi(n-1) \le \phi(n) - 1$, so $\lim_{n \to \infty} m_{\phi(n)-1} \le \min(f(y^*), t_\infty)$ holds because $(m_{\phi(n)-1})_{n \ge 1}$ is non-increasing. The previous arguments show that $\limsup_{n \to \infty} z_{\phi(n)-1} \le 0$.

Moreover, one can use (Vazquez and Bect, 2010b, Proposition 10, (i) ⇒ (ii)) similarly to show that $\sigma_{\phi(n)-1}^2(x_{\phi(n)}) \to 0$ and therefore

$$v_{\phi(n)-1} = \gamma(z_{\phi(n)-1}, \sigma_{\phi(n)-1}^2(x_{\phi(n)})) \le \gamma\left(\sup_{k \ge n} z_{\phi(k)-1}, \sigma_{\phi(n)-1}^2(x_{\phi(n)})\right) \to 0,$$

since $\gamma$ is non-decreasing with respect to its first argument and continuous.

<figure><img src="image_11.png" /></figure>

**Proof of Proposition 17.** This is an adaptation of Theorem 6 from (Vazquez and Bect, 2010b). For $f \in \mathcal{H}(X)$, write $(x_n)_{n \ge 1}$ for the corresponding sequence $(X_n)_{n \ge 1}$ generated by EGO-R. Moreover, write $s_n = s_{R_n,n}$ for the reGP predictor at the step $n$ to avoid cumbersome notations. Then, for $x \in X$, write $\rho_{n,t_n}(x) = \gamma(m_n - s_n(x), \sigma_n^2(x))$ for the expected improvement under the reGP predictive distribution, with $\gamma$ the function defined in Proposition 1.

Suppose that there exists some $x_0 \in X$ such that $\sigma_n^2(x_0) \ge C_1 > 0$. The sequence $(m_n)_{n \ge 1}$ converges and the reproducing property (7) yields

$$|s_n(x_0)| \le \sqrt{k(x_0, x_0)} \|s_n\|_{\mathcal{H}(X)} \le \sqrt{k(x_0, x_0)} \|f\|_{\mathcal{H}(X)},$$

so the sequence $(|m_n - s_n(x_0)|)_{n \ge 1}$ is bounded by, say $C_2$. We have then

$$v_n = \sup_{x \in X} \rho_{n,t_n}(x) \ge \gamma(m_n - s_n(x_0), \sigma_n^2(x_0)) \ge \gamma(-C_2, C_1) > 0$$

by Proposition 1. But this yields a contradiction with Lemma 29, so the decreasing sequence $(\sigma_n^2)_{n \ge 1}$ converges pointwise on $X$ to zero. Proposition 10 from Vazquez and Bect (2010b) then implies that every $x \in X$ is adherent to $\{x_n\}$. ■

**Lemma 30** Assume that $b$ is finite, and that either $a$ is finite too or $\int_{-\infty}^0 F_P(u) du = \int_{-\infty}^0 |u| P(du)$ is finite. Then

$$S_{a,b}^{\text{CRPS}}(P, z) = (b - a \vee z)_+ + R_2(a, b) - 2R_1(a \vee z, b), \quad (52)$$

where

$$R_q(a, b) = \int_{-\infty}^{+\infty} 1_{a \le u \le b} F_P(u)^q du. \quad (53)$$

<footnote n="1">1. Lemma 28 yields $\liminf_{n \to \infty} s_{n-1}(x_{\psi(n)}) \ge \min(f(y^*), t_\infty)$ with $\psi(n) = \max\{\phi(k), \phi(k) \le n\}$, and the claim follows by extracting a $\phi$-subsequence.</footnote>
an R-element of $W_2^{\beta+d/2}(\mathcal{X})$, then Proposition 12 ensures the existence of a minimum $W_2^{\beta+d/2}(\mathcal{X})$-norm element $s_{R,\mathcal{X}}^{(\beta)}$ of the set of functions $g \in W_2^{\beta+d/2}(\mathcal{X})$ coinciding with $f$ on $\mathcal{X}_0$ and such that $f(x) \in R_j \Leftrightarrow g(x) \in R_j$, for all $j \ge 1$ and $x \in \mathcal{X}$. The following proposition is an adaptation of Proposition 25 for bounding the error of the reGP predictor on $\mathcal{X}_0$. It relies on membership in the sense of Definition 26.

**Proposition 27** Suppose that $\mathcal{X}$ is bounded and measurable. Let $0 < \beta < \alpha < +\infty$ be such that $f$ is an R-element of $W_2^{\beta+d/2}(\mathcal{X})$ and $B \subset \mathcal{X}_0$ satisfy Assumption 9. If $k$ has regularity $\alpha$, then it holds that:

$$ \|f - s_{R,n}\|_{L^{\infty}(B)} \lesssim h_n^{\beta}(1 + \rho_n^{\alpha-\beta}) \|s_{R,\mathcal{X}}^{(\beta)}\|_{W_2^{\beta+d/2}(\mathcal{X})} \quad (29) $$

where $\lesssim$ denotes inequality up to a constant, that does not depend on $f$, $n$ or $(x_n)$. In particular, if $(x_n)_{n \ge 1}$ is quasi-uniform, then it holds that:

$$ \|f - s_{R,n}\|_{L^{\infty}(B)} \lesssim n^{-\beta/d} \|s_{R,\mathcal{X}}^{(\beta)}\|_{W_2^{\beta+d/2}(\mathcal{X})} \quad (30) $$

Contrasting Proposition 27 with Proposition 25 suggests an improved fit with reGP in $\mathcal{X}_0$ for two reasons. First, Proposition 27 applies to a larger class of functions since $f$ can be an R-element of $W_2^{\beta+d/2}(\mathcal{X})$ without lying in $W_2^{\beta+d/2}(\mathcal{X})$. Furthermore, Proposition 15 shows that $\|s_{R,\mathcal{X}}^{(\beta)}\|_{W_2^{\beta+d/2}(\mathcal{X})}$ can be arbitrarily smaller than $\|f\|_{W_2^{\beta+d/2}(\mathcal{X})}$ when $f \in W_2^{\beta+d/2}(\mathcal{X})$.

## B.3 reGP with Estimated Parameters

The reGP predictive distribution is obtained by conditioning the GP on the relaxed observations and the model parameters given by (14). The results from Section 3.3.2 do not take parameter selection into account. This appendix presents adaptations of recent theoretical results about parameter selection of Matérn covariance functions for Gaussian process interpolation of deterministic functions.

We still assume a zero mean GP and use the notations from Section 3.3.2, but in this appendix, we consider the RKHS $\mathcal{H}(\mathcal{X})$ with reproducing kernel $r$, not $k$, and the posterior variance will be $\tau^2 \sigma_n^2(\cdot)$, not $\sigma_n^2(\cdot)$ (recall from Section 2.1 that the covariance function is $k(\cdot, \cdot) = \tau^2 r(\cdot - \cdot)$).

First, we consider the maximum likelihood estimation of the variance parameter $\tau^2$. Suppose that $r$ is a Matérn correlation function with known regularity parameter $\nu$ and range parameters $(\rho_1, \dots, \rho_d)$. Using the standard Gaussian process model, and writing $z_n = (f(x_1), \dots, f(x_n))^T$, it is well-known that the maximum likelihood estimate of $\tau^2$ is given by:

$$ \hat{\tau}_n^2 = \frac{\|s_{z_n}\|_{\mathcal{H}(\mathcal{X})}^2}{n}. $$

The estimated posterior variance at $x \in \mathcal{X}$ is $\hat{\tau}_n^2 \sigma_n^2(x)$. Suppose that $(x_n)_{n \ge 1}$ is dense in $\mathcal{X}$. It holds that $\|s_{z_n}\|_{\mathcal{H}(\mathcal{X})} \to \|f\|_{\mathcal{H}(\mathcal{X})}$ if $f \in \mathcal{H}(\mathcal{X})$ (see, e.g., Iske, 2018, Theorem 8.37). Therefore, if $f \in \mathcal{H}(\mathcal{X})$ and $x \in \mathcal{X}$, then (8) implies that:

$$ \frac{|f(x) - s_{z_n}(x)|}{\hat{\tau}_n \sigma_n(x)} \le C \sqrt{n} $$

for some $C > 0$ and $n$ large enough and using the convention $0/0 = 1$. In other words, the GP can be overconfident by a factor of magnitude at most $\sqrt{n}$ (Karvonen et al., 2020; Karvonen, 2022).
**Proof.**

$$
\begin{align*}
S_{a,b}^{\text{tCRPS}}(P, z) &= \int_{-\infty}^{+\infty} \mathbb{1}_{a \le u \le b} (\mathbb{1}_{z \le u} - F_P(u))^2 \, du \\
&= \int_{-\infty}^{+\infty} \mathbb{1}_{a \le u \le b} (\mathbb{1}_{z \le u} + F_P(u)^2 - 2 \mathbb{1}_{z \le u} F_P(u)) \, du \\
&= \int_{-\infty}^{+\infty} (\mathbb{1}_{a \lor z \le u \le b} + \mathbb{1}_{a \le u \le b} F_P(u)^2 - 2 \mathbb{1}_{a \lor z \le u \le b} F_P(u)) \, du \\
&= (b - a \lor z)_+ + R_2(a,b) - 2R_1(a \lor z, b).
\end{align*}
$$

**Lemma 31** Let $a,b \in \mathbb{R}$ with $a \le b$. Let $q \ge 1$. Then

$$
R_q(a,b) = b - a + \mathrm{EI}_q^{\uparrow}(P,b) - \mathrm{EI}_q^{\uparrow}(P,a). \quad (54)
$$

**Proof.**

$$
\begin{align*}
R_q(a,b) &= \int \mathbb{1}_{a \le u \le b} F_P(u)^q \, du \\
&= \int \mathbb{1}_{a \le u \le b} \prod_{j=1}^q E(\mathbb{1}_{N_j \le u}) \, du \quad \text{with } N_j \stackrel{\text{iid}}{\sim} P \\
&= E\left(\int \mathbb{1}_{a \lor N_1 \lor \cdots \lor N_q \le u \le b} \, du\right) \\
&= E\left((b - a \lor N_1 \lor \cdots \lor N_q)_+\right) \\
&= b - a + \mathrm{EI}_q^{\uparrow}(P,b) - \mathrm{EI}_q^{\uparrow}(P,a).
\end{align*}
$$

**Proof of Proposition 21.** The first result is given by Lemma 30 and Lemma 31.

Then using the dominated convergence theorem, it is easy to see that, when $a \to -\infty$

$$
\mathrm{EI}_q^\uparrow(P,a) = \mathrm{E}(N_1 \vee \cdots \vee N_q) - a + o(1), \quad (55)
$$

and therefore

$$
R_q(-\infty, b) = \lim_{a \to -\infty} R_q(a, b) = b + \mathrm{EI}_q^{\uparrow}(P, b) - \mathrm{E}(N_1 \vee \cdots \vee N_q), \quad (56)
$$

which gives the second statement.

Finally, a change of variable gives

$$
\int_a^{+\infty} (F_P(u) - \mathbb{1}_{z \le u})^2 du = \int_{-\infty}^{-a} (F_P(u) - P(U = -u) - \mathbb{1}_{-z < u})^2 du,
$$

and the last statement follows by observing that a probability measure admits at most a countable
number of atoms.
$\blacksquare$
If $f \notin \mathcal{H}(X)$, then (8) is not valid and one can typically expect larger errors. Karvonen et al. (2020) study how $\|s_{z_n}\|_{\mathcal{H}(X)}$ diverges to compensate in this case.

Returning to relaxed Gaussian process interpolation with relaxation range $R$, notice that the likelihood can be optimized with respect to the relaxed observation independently of $\tau^2$. Then, given the optimal relaxed observations $z_n^*$, the reGP variance estimate writes

$$ \hat{\tau}_{n,R}^2 = \frac{\|s_{R,n}\|_{\mathcal{H}(X)}^2}{n}, $$

where $s_{R,n} = s_{z_n}^*$ is the relaxed predictor. Using (7) and (10), it is easy to show that $\hat{\tau}_{n,R}^2 \le \hat{\tau}_n^2$. The following proposition shows that the situations where $\hat{\tau}_{n,R}^2 = \mathcal{O}(n^{-1})$ are characterized by membership of $f$ in $\mathcal{H}(X)$ in the sense of Definition 26.

**Proposition 28** Suppose that $f$ is an $R$-element of the space of continuous functions from $X$ to $\mathbb{R}$ and that $(x_n)_{n \ge 1}$ is dense in $X$. Using the notations of Proposition 12, it holds that $\|s_{R,n}\|_{\mathcal{H}(X)} \to \|s_{R,X}\|_{\mathcal{H}(X)}$ if $f$ is an $R$-element of $\mathcal{H}(X)$ and $\|s_{R,n}\|_{\mathcal{H}(X)} \to +\infty$ otherwise.

Coordinating Proposition 28 with (20) then shows that, for $x \in X_0$, the reGP predictor can be over-confident by a factor of magnitude at most

$$ \sqrt{n} \frac{\sigma_{n,0}(x)}{\sigma_n(x)}, $$

when $f$ is an $R$-element of $\mathcal{H}(X)$. As mentioned in Remark 16, empirical and theoretical evidence about the screening effect suggests that the ratio $\sigma_{n,0}(x)/\sigma_n(x)$ is asymptotically close to one, for finite-smoothness covariance functions. (Alternatively, matching bounds on power functions, such as Equation 16 and Theorem 4.4 from Karvonen et al., 2020, can be used to prove the existence of a constant $C$ such that $\sigma_{n,0}(x) \le C\sigma_n(x)$, for $n$ large enough, in the case of a quasi-uniform sequence of points.) In this case, reGP can be overconfident by a factor of magnitude at most $\sqrt{n}$ as for vanilla GPs.

If $f$ is not an $R$-element of $\mathcal{H}(X)$, then (20) does not hold anymore and Proposition 28 shows that $\|s_{R,n}\|_{\mathcal{H}(X)}$ diverges to compensate. The following proposition is an adaptation of Proposition 4.5 from Karvonen et al. (2020), giving an upper bound of $\hat{\tau}_{n,R}^2$ for quasi-uniform sequences.

**Proposition 29** Suppose that $X$ is bounded and measurable, that $(x_n)_{n \ge 1}$ is quasi-uniform, and that there exists a $\beta > 0$ such that $f$ is an $R$-element of $W_2^{\beta+d/2}(X)$. Then, with a Matérn covariance function with regularity $\nu > \beta$, it holds that:

$$ \hat{\tau}_{n,R}^2 \lesssim n^{2(\nu-\beta)/d-1} \|s_{R,X}^{(\beta)}\|_{W_2^{\beta+d/2}(X)}^2, \quad (31) $$

with $s_{R,X}^{(\beta)}$ defined in the paragraph preceding Proposition 27.

Following Wang and Jing (2022)—but using a slightly different convention—we define the smoothness $\beta(f) = \sup\{\beta > 0 : f \in W_2^{\beta+d/2}(X)\}$ of a function $f$. Given a relaxation range $R$, we use the definition:

$$ \beta_R(f) = \sup \left\{ \beta > 0 : f \text{ is an } R\text{-element of } W_2^{\beta+d/2}(X) \right\}, \quad (32) $$

taking the relaxation into account. It holds that $\beta(f) = \beta_0(f) \le \beta_R(f)$.
**Appendix C. Optimization Benchmark Results**

The results are provided in Figure 10, Figure 11, Figure 12, and Figure 13, for the other test func-
tions from Table 1, using the same format as in Figure 8. Observe that the two heuristics for reGP
yield (sometimes very) substantial improvements on Zakharov (4), Zakharov (10), Three-hump
Camel, Perm (4), Perm (6), Rosenbrock (4), and Rosenbrock (6). However, only the “Concentra-
tion” heuristic yields clear benefits for Zakharov (6), Dixon-Price (4), Dixon-Price (6), and Rosen-
brock (10). Furthermore, the “Concentration” heuristic underperforms slightly on the (multimodal)
Shekel problems. Finally, the EGO and EGO-R algorithms yield indistinguishable results on the
remainings cases.
Let $v > \beta_R(f) > 0$ and assume the supremum (32) defining $\beta_R(f)$ is reached for simplicity. The rate in the upper bound (31) is optimized by taking $\beta = \beta_R(f)$. Under more restrictive assumptions on the spectrum and the support of $f$, Karvonen et al. (2020) give matching lower bounds for traditional Gaussian process interpolation and discuss uncertainty quantification for functions outside the RKHS. Obtaining a similar result for reGP is more involved and thus left for future work. Considering (30) and Theorem 4.4 from Karvonen et al. (2020), it can be seen that a class of functions such that the bound (31) with $\beta = \beta_R(f)$ is sharp would also essentially be such that reGP can be overconfident by a factor of magnitude at most $\sqrt{n}$ on $\mathbb{X}_0$ for quasi-uniform observations.

Finally, we consider the selection of the smoothness parameter $v$ given fixed prior variance and range parameters. For Gaussian process interpolation, Karvonen (2023) gives asymptotic bounds for maximum likelihood estimates of $v$ with quasi-uniform observations from a deterministic function. The Matérn covariance function depends continuously on $v$. Therefore, for $0 < v_{\min} < v_{\max} < +\infty$ and $\Theta = [v_{\min}, v_{\max}]$, the minimum of (14) is reached on $\Theta \times C_{R,n}$ according to the paragraph following Definition 6. For each $n$, let $(\hat{v}_n^R, \underline{z}_n^*)$ be such a minimizer, with ties broken arbitrarily. The following asymptotic lower bound is an adaptation of the first statement of Theorem 3.12 from Karvonen (2023).

**Proposition 30** Let $0 < v_{\min} < v_{\max} < +\infty$ and $\Theta = [v_{\min}, v_{\max}]$ and assume $\beta_R(f) + d/2 \in \Theta$, with $\beta_R(f)$ given by (32). If $\mathbb{X}$ satisfies Assumption 9 and the sequence $(x_n)_{n \ge 1}$ is quasi-uniform, it holds that $\liminf \hat{v}_n^R \ge \beta_R(f) + d/2$.

A smoothness parameter estimate based on reGP is asymptotically strictly larger than $\beta_R(f)$. Consequently, the function $f$ is not an $R$-element of the RKHS. This echoes the results obtained for Gaussian process interpolation in the previous reference. These results must be seen considering that parameters are selected using the likelihood of a model whose sample paths are almost surely not contained in the RKHS (see, e.g., Lukić and Beder, 2001). Petit (2025, Section 6) studies asymptotically the effect of a joint estimation of $v$ and $\tau^2$ in a simplified Gaussian process interpolation framework.

Assume the supremum (32) is reached for simplicity and $0 < \beta_R(f) < +\infty$. With a Matérn covariance function with regularity $v > 0$, quasi-uniform observations, and $B \subset \mathbb{X}_0$ satisfying Assumption 9, the upper bounds (21) and (30) yield:

$$ \left\{ \begin{array}{ll} \|f - s_{R,n}\|_{L^\infty(B)} \lesssim n^{-v/d} & \text{if } v \le \beta_R(f), \\ \|f - s_{R,n}\|_{L^\infty(B)} \lesssim n^{-\beta_R(f)/d} & \text{otherwise.} \end{array} \right. \qquad (33) $$

As with standard Gaussian processes, this suggests that, for quasi-uniform observations, the most accurate reGP models on $\mathbb{X}_0$ are obtained by taking $v \ge \beta_R(f)$. Proposition 30 shows that reGP selects such models automatically.

The previous remark also applies with $R = \emptyset$, recovering a conclusion of Karvonen (2023, Section 4.1) for Gaussian process interpolation. Consequently, the inequality $\beta(f) = \beta_0(f) \le \beta_R(f)$ and (33) suggest an improved fit with reGP compared to a traditional Gaussian process when both models are used with estimated smoothness parameters.

Finally, note that the quasi-uniformity assumption is rather strong and prevents applications to active learning tasks (see Wynne et al., 2021, Section 5.2, for a possible application to a Bayesian optimization policy designed to ensure quasi-uniformity).
Zakharov (4)

<figure><img src="image_12.png" /></figure>

Zakharov (10)

<figure><img src="image_13.png" /></figure>

Perm (4)

<figure><img src="image_14.png" /></figure>

Zakharov (6)

<figure><img src="image_15.png" /></figure>

Three-hump Camel

<figure><img src="image_16.png" /></figure>

Perm (6)

<figure><img src="image_17.png" /><figcaption>Figure 10: Same as Figure 8. The gray dashed line stands for the Dual simulated Annealing algorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the “Constant” and “Concentration” heuristics.</figcaption></figure>
## Appendix C. Proofs

**Lemma 31** (Aronszajn, 1950, Section 1.5) Let $k: \mathbb{X} \times \mathbb{X} \to \mathbb{R}$ be a positive-definite covariance function, $\mathbb{U} \subset \mathbb{X}$, and $\mathcal{H}(\mathbb{U})$ be the RKHS attached to the restriction of $k$ to $\mathbb{U} \times \mathbb{U}$. The RKHS $\mathcal{H}(\mathbb{U})$ is the space of restrictions of functions from $\mathcal{H}(\mathbb{X})$ and the norm of $g \in \mathcal{H}(\mathbb{U})$ is given by

$$ \inf_{h \in \mathcal{H}(\mathbb{X}), h|_{\mathbb{U}}=g} \|h\|_{\mathcal{H}(\mathbb{X})}. $$

**Proof of Proposition 3.** First the existence and the uniqueness of the solution are given by the first statement of Proposition 12 (with $\mathcal{H}_{R,n} = \mathcal{H}_{R, \{x_1, \dots, x_n\}}$).

Furthermore let $\underline{z} \in \mathbb{R}^n$ and write $\underline{\alpha} = K_n^{-1}\underline{z}$, the reproducing property (6) gives

$$ \|s_{\underline{z}}\|_{\mathcal{H}(\mathbb{X})}^2 = \underline{\alpha}^{\mathrm{T}} K_n \underline{\alpha} = \underline{z}^{\mathrm{T}} K_n^{-1} \underline{z}, \quad (34) $$

and therefore

$$ \min_{h \in \mathcal{H}_{R,n}} \|h\|_{\mathcal{H}(\mathbb{X})}^2 = \inf_{\underline{z} \in C_{R,n}} \min_{h \in \mathcal{H}(\mathbb{X}), h(x_n)=\underline{z}} \|h\|_{\mathcal{H}(\mathbb{X})}^2 = \inf_{\underline{z} \in C_{R,n}} \underline{z}^{\mathrm{T}} K_n^{-1} \underline{z}, $$

where the last infimum is uniquely reached by the evaluation of the solution on $x_n$.

**Proof of Proposition 4.** Write $K_{m,n}$ for the covariance matrix of the random vector $(\underline{z}'_m^{\mathrm{T}}, \underline{z}'_n^{\mathrm{T}})^{\mathrm{T}}$. Using the equalities (4) and (34), and a slight abuse of notation by dropping irrelevant constants with respect to $\underline{z}'$ and $\underline{z} \in C_{R,n}$, we have

$$ -2 \ln \left( p \left( \underline{z}', \underline{z} \mid \underline{Z}_n \in C_{R,n} \right) \right) = \left( \underline{z}'^{\mathrm{T}}, \underline{z}^{\mathrm{T}} \right) K_{m,n}^{-1} \left( \underline{z}'^{\mathrm{T}}, \underline{z}^{\mathrm{T}} \right)^{\mathrm{T}} = \min_{h \in \mathcal{H}(\mathbb{X}), h(x_n)=\underline{z}, h(x_m)=\underline{z}'} \|h\|_{\mathcal{H}(\mathbb{X})}^2. $$

This gives

$$ \inf_{\underline{z}' \in \mathbb{R}^m, \underline{z} \in C_{R,n}} -2 \ln \left( p \left( \underline{z}', \underline{z} \mid \underline{Z}_n \in C_{R,n} \right) \right) = \min_{h \in \mathcal{H}(\mathbb{X}), h(x_n) \in C_{R,n}} \|h\|_{\mathcal{H}(\mathbb{X})}^2, $$

which is reached by taking $\underline{z} = \underline{z}_n^*$ and $\underline{z}' = (s_{\underline{z}_n^*}^{\mathrm{T}}(x_1'), \dots, s_{\underline{z}_n^*}^{\mathrm{T}}(x_m'))^{\mathrm{T}}$.

**Proof of Proposition 10.** First, one has

$$ \sup_{x \in \mathbb{X}} \sigma_n(x) = \sup_{x \in \mathbb{X}} \sup_{\|f\|_{\mathcal{H}(\mathbb{X})}=1} |f(x) - s_{\underline{z}_n}(x)| = \sup_{\|f\|_{\mathcal{H}(\mathbb{X})}=1} \|f - s_{\underline{z}_n}\|_{L^\infty(\mathbb{X})}. $$

Now, let $f \in \mathcal{H}(\mathbb{X})$ such that $\|f\|_{\mathcal{H}(\mathbb{X})} = 1$, and $\mathbb{X}^o$ be the interior of $\mathbb{X}$. The boundary of $\mathbb{X}^o$ is the one of $\mathbb{X}$ under Assumption 9, and the Sobolev space $W_2^{\alpha+d/2}(\mathbb{X}^o)$ defined by (15) is norm-equivalent to the Sobolev–Slobodeckij space (see, e.g., Proposition 3.4 of Di Nezza et al., 2012, for a statement on $\mathbb{R}^d$ and Theorem 1.4.3.1 of Grisvard, 1985, for the existence of an extension operator).

Then, one can apply Theorem 4.1 from Arcangéli et al. (2007) to $f - s_{\underline{z}_n}$ restricted to $\mathbb{X}^o$ to show that, for $h_n$ lower than some positive $h_0$ (not depending on $f$ or $(x_n)_{n \ge 1}$), we have:

$$ \|f - s_{\underline{z}_n}\|_{L^\infty(\mathbb{X})} = \|f - s_{\underline{z}_n}\|_{L^\infty(\mathbb{X}^o)} \lesssim h_n^\alpha \|f - s_{\underline{z}_n}\|_{W_2^{\alpha+d/2}(\mathbb{X}^o)} \lesssim h_n^\alpha \|f - s_{\underline{z}_n}\|_{\mathcal{H}(\mathbb{X})} \lesssim h_n^\alpha \quad (35) $$
Six-hump Camel

<figure><img src="image_18.png" /></figure>

Dixon-Price (6)

<figure><img src="image_19.png" /></figure>

Rosenbrock (6)

<figure><img src="image_20.png" /></figure>

Cross-in-Tray

<figure><img src="image_21.png" /></figure>

Rosenbrock (4)

<figure><img src="image_22.png" /></figure>

Rosenbrock (10)

<figure><img src="image_23.png" /><figcaption>Figure 11: Same as Figure 8. The gray dashed line stands for the Dual simulated Annealing algorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the “Constant” and “Concentration” heuristics.</figcaption></figure>
by continuity of $f - s_{z_n}$, since $\|\cdot\|_{W_2^{\alpha+d/2}(X^0)} \le \|\cdot\|_{W_2^{\alpha+d/2}(X)}$ due to the definition (15), $W_2^{\alpha+d/2}(X)$ being norm equivalent to $\mathcal{H}(X)$, and because of the projection interpretation of $s_{z_n}$ (see, e.g., Wendland, 2004, Theorem 13.1). Finally, one can get rid of the condition $h_n \le h_0$ for simplicity by possibly increasing the constant, since $\sup_{x \in X} \sigma_n(x) \le \sup_{x \in X} \sqrt{k(x, x)} < +\infty$ by compacity. ■

**Proof of Proposition 12.** First observe that $\mathcal{H}_{R,U}$ is not empty since it contains $f$. Furthermore, it is easy to verify that $\mathcal{H}_{R,U}$ is convex and that it is closed since pointwise evaluation functionals are continuous on an RKHS. The problem is then the one of projecting the null function on a convex closed subset; hence the existence and the uniqueness.

Then, the function $s_{R,n}$ is the projection of the null function on the closed convex set $\mathcal{H}_{R,n}$. Moreover, the sequence $(\mathcal{H}_{R,n})_{n \ge 1}$ is non-increasing, so the sequence $(s_{R,n})_{n \ge 1}$ converges in $\mathcal{H}(X)$ to the projection of the null function on $\cap_{n \ge 1} \mathcal{H}_{R,n}$ (see, e.g., Brezis, 2011, Exercice 5.5), i.e., the solution of (19), with $U = \{x_n\}$. But this last solution is also the solution on the closure since it satisfies the constraints by continuity. ■

**Proof of Proposition 13.** Define $x_n^0$ and $z_n^0$ to be data points within $X_0$, and $s_{x_n^0, z_n^0}$ the associated (interpolation) predictor, i.e., the solution of (7). Take the null function if there are no observations within $B$. Observing that $s_{x_n^0, z_n^0}$ interpolates $s_{R,n}$, we have:

$$ \begin{aligned} |f(x) - s_{R,n}(x)| & \le |f(x) - s_{x_n^0, z_n^0}(x)| + |s_{x_n^0, z_n^0}(x) - s_{R,n}(x)| \\ & \le \sigma_{n,0}(x) \|f\|_{\mathcal{H}(X_0)} + \sigma_{n,0}(x) \|s_{R,n}\|_{\mathcal{H}(X_0)} \\ & \le 2\sigma_{n,0}(x) \|s_{R,n}\|_{\mathcal{H}(X)}, \end{aligned} $$

since $f$ coincides with $s_{R,X}$ on $X_0$, $s_{R,X} \in \mathcal{H}_{R,n}$, and by Lemma 31. ■

**Lemma 32** Use the notation (36). If $k$ has smoothness $\alpha > 0$, then there exists $h_0 > 0$ depending only on $\alpha$ such that, for all $x, y \in X$ satisfying $\|x - y\| \le h_0$, we have

$$ \begin{aligned} \|\delta_y - \delta_x\|_{\mathcal{H}^*(X)} & \lesssim \|x - y\|^\alpha, \quad \text{for } \alpha < 1, \\ \|\delta_y - \delta_x\|_{\mathcal{H}^*(X)} & \lesssim \sqrt{|\ln(\|x - y\|)|} \|x - y\|, \quad \text{for } \alpha = 1, \end{aligned} $$

and

$$ \|\delta_y - \delta_x\|_{\mathcal{H}^*(X)} \lesssim \|x - y\|, \quad \text{for } \alpha > 1. $$

**Proof.** Since equivalent norms give equivalent operator norms on the topological dual of a normed space, it suffices to show the result for a unit-variance isotropic Matérn covariance function (3) of regularity $\alpha$.

In this case, we have

$$ \|\delta_y - \delta_x\|^2_{\mathcal{H}^*(X)} = k(x, x) + k(y, y) - 2k(x, y) = 2(1 - r_\alpha(\|x - y\|)), $$

with $r_\alpha$ the corresponding isotropic correlation function. Standard results on principal irregular terms (see, e.g., Stein, 1999, Chapter 2.7) give the results. ■
<figure><img src="image_24.png" /><figcaption>Figure 12: Same as Figure 8. The gray dashed line stands for the Dual simulated Annealing algorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the “Constant” and “Concentration” heuristics.</figcaption></figure>
**Lemma 33** Let $B \subset \mathbb{X}$ satisfying Assumption 9 and $h_{n,B}$ be the fill distance of $X_{n,B} = \{x_1, \dots, x_n\} \cap B$ within $B$, with the convention $h_{n,B} = \text{diam}(B)$ if $X_{n,B}$ is empty. Then, $h_{n,B} \lesssim h_n$.

**Proof.** The idea of the proof is given by Wendland (2004, Lemma 11.31), but it is interlinked with a much more sophisticated construction, so we provide a specific version here for completeness. Adams and Fournier (2003, Section 4.11) show that $B$ satisfies a cone condition with radius $\rho > 0$ and angle $\phi \in (0, \pi/2)$. If $X_{n,B}$ is not empty, then the compacity of $B$ ensures the existence of an $x \in B$ such that $d(x, X_{n,B}) = h_{n,B}$. (If $X_{n,B}$ is empty, then the rest of the proof is also valid taking an arbitrary $x \in B$.)

A cone $C$ originating from $x$ with angle $\phi$ and radius $\delta = \min(h_{n,B}, \rho)$ is contained in $B$ and its interior do not contains observations. Furthermore, Wendland (2004, Lemma 3.7) shows that there exists a $y \in C$ such that the open ball $B(y, \delta \sin(\phi)(1+\sin(\phi))^{-1})$ is subset of $C$, and therefore contains no observations as well. Thus, $h_n \ge \delta \sin(\phi)(1+\sin(\phi))^{-1}$. Now, if $h_{n,B} \le \rho$, then the desired result follows. If not, the result holds as well since $h_{n,B} \le \text{diam}(B)$. ■

**Proof of Proposition 14.** For the first assertion, suppose first that $B$ contains observations and let $\sigma_{n,B}$ be the power function using only the observations within $B$. It can be shown that the restriction of the kernel to $B$ also has regularity $\alpha$ using (15) and Lemma 31. Using Proposition 13, the inequality $\sigma_{n,0} \le \sigma_{n,B}$ given by the projection residuals interpretation, and applying Proposition 10 to $B$ yields a bound depending on the fill distance $h_{n,B}$ of $\{x_1, \dots, x_n\} \cap B$ within $B$. Finally, Lemma 33 allows us to conclude. Dealing with the case where $B$ contains no observations using (6) makes no difficulty.

Regarding the second assertion, $f$ is continuous so the sets $X_j$ are compact for $j \ge 1$. In addition, they are disjoint so

$$ \delta = \min_{1 \le j < p} \inf_{x \in X_j, y \in X_p} \|x - y\| > 0. $$

Suppose now that $h_n < \delta$ and let $j \ge 1$, $x \in X_j$ and $1 \le i \le n$ the index of the closest $x_i$ to $x$. By definition, $\|x - x_i\| \le h_n$ and therefore $x_i \in X_0 \cup X_j$. Let $\mathcal{H}^*(\mathbb{X})$ be the (topological) dual of $\mathcal{H}(\mathbb{X})$ and

$$ \delta_y: h \in \mathcal{H}(\mathbb{X}) \mapsto h(y), \qquad (36) $$

which lies in $\mathcal{H}^*(\mathbb{X})$ for all $y \in \mathbb{X}$. Then using the reproducing property (6), we have

$$ |s_{R,n}(x_i) - s_{R,n}(x)| \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(\mathbb{X})} \|s_{R,n}\|_{\mathcal{H}(\mathbb{X})} \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(\mathbb{X})} \|s_{R,n}\|_{\mathcal{H}(\mathbb{X})}, $$

and therefore

$$ d(s_{R,n}(x), R_j) \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(\mathbb{X})} \|s_{R,n}\|_{\mathcal{H}(\mathbb{X})} + d(s_{R,n}(x_i), R_j). $$

Now, if $x_i \in X_j$, then $d(s_{R,n}(x_i), R_j) = 0$. Otherwise $x_i \in X_0$ necessarily and, then, using the fact that $s_{R,\mathbb{X}}(x) \in R_j$, we have:

$$ d(s_{R,n}(x_i), R_j) \le |s_{R,n}(x_i) - s_{R,\mathbb{X}}(x)| = |s_{R,\mathbb{X}}(x_i) - s_{R,\mathbb{X}}(x)| \le \|\delta_{x_i} - \delta_x\|_{\mathcal{H}^*(\mathbb{X})} \|s_{R,\mathbb{X}}\|_{\mathcal{H}(\mathbb{X})}. $$

So one can use Lemma 32 along with the previous statements to conclude if $h_n < \min(\delta, h_0)$. (It may be necessary to reduce $\min(\delta, h_0)$ to use the monotonicity of $x \mapsto x\sqrt{-\ln(x)}$, for $x > 0$ close to zero.)
<figure><img src="image_25.png" /><figcaption>Figure 13: Same as Figure 8. The gray dashed line stands for the Dual simulated Annealing algorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the “Constant” and “Concentration” heuristics.</figcaption></figure>
Finally, treating the case $h_n \ge \min(\delta, h_0)$ is straightforward using the fact that $\sup_{x \in X} \sqrt{k(x, x)}$ is finite thanks to the compacity of $X$ and $d(s_{R,n}(x), R_j) \le |s_{R,n}(x) - s_{R,X}(x)|$ for $j \ge 1$ and $x \in X_j$. ■

**Lemma 34** If $g, h \in W_2^\gamma(X)$ for $\gamma > d/2$, then $gh \in W_2^\gamma(X)$.

**Proof.** By the definition (15) of $W_2^\gamma(X)$, the functions $g$ and $h$ can be extended into functions on $\mathbb{R}^d$, having their product in $W_2^\gamma(\mathbb{R}^d)$ (Strichartz, 1967, Theorem 2.1). Taking the restrictions shows the desired result. ■

**Proof of Proposition 15.** We use a bump function argument. Let $B(x_0, r) \subset X_j$ (with $r > 0$) be an open ball. There exists a $C^\infty$ function $\phi: \mathbb{R}^d \to \mathbb{R}$ such that

$$ \begin{cases} 0 \le \phi \le 1, \\ \phi(x) = 1 & \text{only if } x = 0, \\ \phi(x) = 0 & \text{if } x \in X \setminus B(0, r). \end{cases} $$

Let $c \in R_j \setminus \{f(x_0)\}$, $\phi_n = \phi(n \cdot x_0)$ as a function on $X$, and $f_n = (1 - \phi_n)f + c\phi_n$, for $n \ge 1$. We have $\phi_n \in W_2^{\alpha+d/2}(\mathbb{R}^d)$ as a function on $\mathbb{R}^d$, so it belongs to $W_2^{\alpha+d/2}(X)$ as a function on $X$, and Lemma 34 ensures that $f_n \in \mathcal{H}(X)$. Moreover, it is easy to check that $f_n \in \mathcal{H}_{R,X}$. Observe that the sequence $(f_n)_{n \ge 1}$ converges pointwise to a discontinuous function that lies thus outside $\mathcal{H}(X)$.

Suppose now that $\|f_n\|_{\mathcal{H}(X)} \to +\infty$. Then, one can extract a bounded subsequence of norms and a classical weak compacity argument would yield a weakly convergent subsequence, which is impossible since the pointwise limit is not in $\mathcal{H}(X)$. ■

**Lemma 35** Use the notations from the proof of Proposition 18 and let $(y_n)_{n \ge 1}$ be a sequence in $X$. Assume that the sequence $(y_n)_{n \ge 1}$ is convergent, denote by $y^*$ its limit and assume that $y^*$ is an adherent point of the set $\{x_n\}$. Let $t_\infty = \liminf_{n} t_n$, then

* $s_n(y_n) \to f(y^*)$ if $f(y^*) < t_\infty$,

* $\liminf s_n(y_n) \ge t_\infty$, otherwise.

In particular, we have

$$ \liminf s_n(y_n) \ge \min(f(y^*), t_\infty). $$

**Proof.** Suppose that $y^* \notin \{x_n\}$. Then, let $(x_{\phi(n)})_{n \ge 1}$ be a subsequence converging to $y^*$ and let $\psi(n) = \max\{\phi(k), \phi(k) \le n\}$. We proceed as in Proposition 14 to have:

$$ |s_n(x_{\psi(n)}) - s_n(y_n)| \le \|\delta_{x_{\psi(n)}} - \delta_{y_n}\|_{\mathcal{H}^*(X)} \|s_n\|_{\mathcal{H}(X)} \to 0, \quad (37) $$

thanks to the continuity of $k$ and the inequality $\|s_n\|_{\mathcal{H}(X)} \le \|f\|_{\mathcal{H}(X)}$.

Finally, we have $s_n(x_{\psi(n)}) \ge \min(f(x_{\psi(n)}), t_n)$ by construction and therefore $\liminf s_n(x_{\psi(n)}) \ge \min(f(y^*), t_\infty)$, which gives the second assertion thanks to (37). Observe that $f(x_{\psi(n)}) < t_n$ ultimately if $f(y^*) < t_\infty$ for the first assertion.

If $y^* \in \{x_n\}$, then the result follows similarly. ■
References

R. A. Adams and J. J. F. Fournier. Sobolev spaces. Elsevier, 2003.

R. Arcangéli, M. C. López de Silanes, and J. J. Torrens. An extension of a bound for functions in Sobolev spaces, with applications to (m, s)-spline interpolation and smoothing. *Numerische Mathematik*, 107(2):181–211, 2007.

N. Aronszajn. Theory of reproducing kernels. *Transactions of the American mathematical society*, 68(3):337–404, 1950.

S.-K. Au and J. L. Beck. Estimation of small failure probabilities in high dimensions by subset simulation. *Probabilistic engineering mechanics*, 16(4):263–277, 2001.

S. Ba and V. R. Joseph. Composite Gaussian process models for emulating expensive functions. *The Annals of Applied Statistics*, 6(4):1838–1860, 2012.

F. Bachoc and A. Lagnoux. Posterior contraction rates for constrained deep Gaussian processes in density estimation and classication. *arXiv preprint arXiv:2112.07280*, 2021.

J. Y. Bao, F. Ye, and Y. Yang. Screening effect in isotropic Gaussian processes. *Acta Mathematica Sinica*, 36(5), 2020.

R. Benassi, J. Bect, and E. Vazquez. Bayesian optimization using sequential Monte Carlo. In *Learning and Intelligent Optimization – 6th International Conference, LION 6*, pages 339–342, Paris, France, January 16-20 2012. Springer.

A. Berlinet and C. Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2004.

E. Bodin, M. Kaiser, I. Kazlauskaitė, Z. Dai, N. Campbell, and C. H. Ek. Modulating surrogates for Bayesian optimization. In *International Conference on Machine Learning*, pages 970–979. PMLR, 2020.

I. Bogunovic, J. Scarlett, A. Krause, and V. Cevher. Truncated variance reduction: A unified approach to Bayesian optimization and level-set estimation. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, editors, *Advances in neural information processing systems*, volume 29, 2016.

H. Brezis. Functional analysis, Sobolev spaces and partial differential equations, volume 2. Springer, 2011.

A. D. Bull. Convergence rates of efficient global optimization algorithms. *Journal of Machine Learning Research*, 12(10), 2011.

R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained optimization. *SIAM Journal on scientific computing*, 16(5):1190–1208, 1995.