Sébastien J Petit
email: sebastien.julien.petit@gmail.com

Julien Bect
email: julien.bect@centralesupelec.fr

Emmanuel Vazquez
email: emmanuel.vazquez@centralesupelec.fr

Relaxed Gaussian Process Interpolation: a Goal-Oriented Approach to Bayesian Optimization

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

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 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 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 spaces 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.

Introduction 1.1 Context and Motivation

Gaussian process (GP) interpolation and regression (see, e.g., [START_REF] Stein | Interpolation of Spatial Data: Some Theory for kriging[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]) 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. The variations on the left overshadow the global minimum on the right. This is the case for instance in the motivating example given by [START_REF] Gramacy | Bayesian treed Gaussian process models with an application to computer modeling[END_REF], 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.

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- given by the model and the gray envelopes represent the associated uncertainties.

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

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 [START_REF] Gramacy | Bayesian treed Gaussian process models with an application to computer modeling[END_REF], which combines a fully Bayesian framework and the use of RJ-MCMC techniques for the inference, or the trust-region method by [START_REF] Eriksson | Scalable global optimization via local Bayesian optimization[END_REF]. [START_REF] Park | Patchwork kriging for large-scale Gaussian process regression[END_REF] 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 [START_REF] Nott | Estimation of nonstationary spatial covariance structure[END_REF], weightings of covariance functions [START_REF] Pronzato | Bayesian local kriging[END_REF][START_REF] Rivoirard | Continuity for kriging with moving neighborhood[END_REF], and convolution techniques (see, e.g., [START_REF] Higdon | A process-convolution approach to modelling temperatures in the north atlantic ocean[END_REF][START_REF] Gibbs | Bayesian Gaussian processes for regression and classification[END_REF][START_REF] Higdon | Space and space-time modeling using process convolutions[END_REF][START_REF] Ver Hoef | Flexible spatial models for kriging and cokriging using moving averages and the fast fourier transform (fft)[END_REF][START_REF] Stein | Nonstationary spatial covariance functions[END_REF]. Let us also mention data-driven aggregation techniques: composite Gaussian process models [START_REF] Ba | Composite Gaussian process models for emulating expensive functions[END_REF], and mixture of experts techniques (see, e.g., [START_REF] Tresp | Mixtures of Gaussian processes[END_REF][START_REF] Rasmussen | Infinite mixtures of Gaussian process experts[END_REF][START_REF] Meeds | An alternative infinite mixture of Gaussian process experts[END_REF][START_REF] Yuan | Variational mixture of Gaussian process experts[END_REF][START_REF] Yang | An efficient em approach to parameter learning of the mixture of Gaussian processes[END_REF][START_REF] Yuksel | Twenty years of mixture of experts[END_REF]. 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.

TRANSFORMATION AND COMPOSITION OF MODELS

A first technique for composition of models consists in using a parametric transformation of a GP [START_REF] Rychlik | Modelling and statistical analysis of ocean-wave data using transformed Gaussian processes[END_REF][START_REF] Snelson | Warped Gaussian processes[END_REF].

Another route is to transform the input domain, using for instance a parametric density [START_REF] Xiong | A non-stationary covariance-based kriging method for metamodelling in engineering design[END_REF], or other parametric transformations involving possible dimension reduction [START_REF] Marmin | Warped Gaussian processes and derivativebased sequential designs for functions with heterogeneous variations[END_REF]. [START_REF] Bodin | Modulating surrogates for Bayesian optimization[END_REF] 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., [START_REF] Damianou | Deep Gaussian processes[END_REF][START_REF] Dunlop | How deep are deep Gaussian processes[END_REF][START_REF] Hebbal | Bayesian optimization using deep Gaussian processes with applications to aerospace system design[END_REF][START_REF] Jakkala | Deep Gaussian processes: A survey[END_REF][START_REF] Bachoc | Posterior contraction rates for constrained deep Gaussian processes in density estimation and classication[END_REF], 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, [START_REF] Picheny | Ordinal Bayesian optimisation[END_REF] 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 [START_REF] Chu | Gaussian processes for ordinal regression[END_REF] for the inference.

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 distributions; 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 GPbased 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 [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] 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 of this algorithm and a numerical benchmark. Finally, Section 6 presents our conclusions and perspectives for future work.

Background and Notations

Gaussian Process Modeling

Consider a real-valued function f : X → R, where X ⊆ R d , and suppose we want to infer f at a given x ∈ X from evaluations of f on a finite set of points x n = (x 1 , . . . , x n ) ∈ X n , n ≥ 1. A standard Bayesian approach to this problem consists in using a GP model ξ ∼ GP (µ, k) as a prior about f , where µ : X → R is a mean function and k : X × X → R is a covariance function, which is supposed to be strictly positive-definite in this article.

The posterior distribution of ξ given

Z n = (ξ (x 1 ), • • • , ξ (x n ))
T is still a Gaussian process, whose mean and covariance functions are given by the standard kriging equations [START_REF] Matheron | The theory of regionalized variables and its applications[END_REF]. More precisely:

ξ | Z n ∼ GP (µ n , k n ) , (1)
with µ n (x) = µ(x) + k (x, x n ) K -1 n (Z n -µ(x n )) (2) and k n (x, y) = k(x, y) -k (x, x n ) K -1 n k (y, x n ) T , (3)
and where µ(

x n ) = (µ(x 1 ), . . . , µ(x n )) T , k (x, x n ) = (k(x, x 1 ), . . . , k(x, x n ))
, and K n is the n×n matrix with entries k(x i , x j ). We shall also use the notation σ 2 n (x) = k n (x, x) for the posterior variance, a.k.a. the kriging variance, a.k.a. the squared power function, so that ξ

(x) | Z n ∼ N µ n (x), σ 2 n (x) .
The functions µ 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 {(µ θ , k θ ) , θ ∈ Θ}. A common approach is to suppose stationarity for the GP, which means choosing a constant mean function µ ≡ c ∈ R and a stationary covariance function k(x, y) = σ 2 r(xy), where r : R d → R is a stationary correlation function.

A correlation function often recommended in the literature [START_REF] Stein | Interpolation of Spatial Data: Some Theory for kriging[END_REF] is the (geometrically anisotropic) Matérn correlation function

r(h) = 2 1-ν Γ(ν) √ 2ν h ρ ν K ν √ 2ν h ρ , h 2 ρ = d ∑ j=1 h 2 [ j] ρ 2 j , (4)
for h = (h [1] , . . . , h [d] ) ∈ R d
, and where Γ is the Gamma function and K ν is the modified Bessel function of the second kind. The parameters to be selected in this case are (σ 2 , ρ 1 , . . . , ρ d , ν) ∈ (0, ∞) d+2 with σ 2 the process variance, ρ i the range parameter along the i-th dimension, and ν a regularity parameter controlling the smoothness of the process. Two other standard covariance functions can be recovered for specific values of ν: the exponential covariance function for ν = 1/2 and the squared-exponential covariance function for ν → ∞.

A variety of techniques for selecting the parameter θ 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 [START_REF] Petit | Model parameters in Gaussian process interpolation: an empirical study of selection criteria[END_REF]. It simply consists in minimizing the negative loglikelihood

L (θ ; Z n ) = -log (p (Z n | θ )) ∝ log (det (K n )) + (Z n -µ(x n )) T K -1 n (Z n -µ(x n )) + constant, ( 5
)
where p stands for the probability density of Z n . Other methods for selecting the parameters include the restricted maximum likelihood method and leave-one-out strategies (see, e.g., [START_REF] Stein | Interpolation of Spatial Data: Some Theory for kriging[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

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 , . . . ∈ X using a GP model ξ 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) [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], which can be expressed as

ρ n (x) = E (m n -ξ (x)) + | Z n , (6)
where

m n = min(ξ (x 1 ), • • • , ξ (x n )).
The EI criterion corresponds to the expectation of the excursion of ξ below the minimum given n observations, and can be written in closed form:

Proposition 1 [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]Vazquez and Bect, 2010b) The EI criterion may be written as

ρ n (x) = γ m n -µ n (x), σ 2 n (x) with γ : (z, s) ∈ R × R + → √ sφ z √ s + zΦ z √ s if s > 0, max(z, 0) if s = 0,
where φ and Φ stand for the probability density and cumulative distribution functions of the standard Gaussian distribution. Moreover, the function γ is continuous, verifies γ(z, s) > 0 if s > 0 and is non-decreasing with respect to z and s on R × 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∈X ρ n (x) ,
the resulting algorithm is generally called the Efficient Global Optimization (EGO) algorithm, as proposed by [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. The EGO algorithm has known convergence properties (Vazquez and Bect, 2010b;[START_REF] Bull | Convergence rates of efficient global optimization algorithms[END_REF].

A variety of other sampling criteria for the minimization problem can be found in the literature (see, e.g., [START_REF] Frazier | A knowledge-gradient policy for sequential information collection[END_REF][START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF][START_REF] Srinivas | Gaussian process optimization in the bandit setting: No regret and experimental design[END_REF][START_REF] Vazquez | A new integral loss function for Bayesian optimization[END_REF], but we shall focus on the EI algorithm in this article.

Reproducing Kernel Hilbert Spaces

Reproducing kernel Hilbert spaces (RKHS, see e.g., [START_REF] Aronszajn | Theory of reproducing kernels[END_REF][START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF] are Hilbert spaces of functions commonly used in the field of approximation theory (see, e.g., [START_REF] Wahba | Spline Models for Observational Data[END_REF][START_REF] Wendland | Scattered data approximation[END_REF]. A Hilbert space H (X) of real-valued functions on X with an inner product (• , •) H (X) is called an RKHS if it has a reproducing kernel, that is, a function k : X × X → R such that k(x, •) ∈ H (X), and

( f , k(x, •)) H (X) = f (x) (7)
(the reproducing property), for all x ∈ X and f ∈ 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 , • • • , x n ) ∈ X n
, and corresponding values z n ∈ R n , suppose we want to find a function h ∈ H (X) such that h(x n ) = (h(x 1 ), . . . , h(x n )) T = z n . Then, the minimum-norm interpolation solution is given by the following proposition.

Proposition 2 [START_REF] Parzen | Statistical inference on time series by Hilbert space methods[END_REF][START_REF] Matheron | Splines and kriging; their formal equivalence[END_REF] The problem

     minimize h H (X) subject to h ∈ H (X) h(x n ) = z n (8)
has a unique solution given by s

z n = k(•, x n )K -1 n z n .
Observe that the solution s z n is equal to the posterior mean (2) when µ = 0. Moreover, for any f ∈ H (X) and x ∈ X, the reproducing property (7) yields the upper bound

f (x) -s z n (x) ≤ σ n (x) f H (X) , (9)
with σ n (x) = k n (x, x). Note that σ n (x) is the worst-case error at x for the interpolation of functions in the unit ball of H (X).

Relaxed Gaussian Process Interpolation

Relaxed Interpolation

The example in the introduction (see Figures 123) 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 = J j=1 R j , where R 1 , . . . , R J ⊂ R are disjoint closed intervals with non-zero lengths. (The set R = [t, +∞) was used in the example of Figure 3).

As above, we shall write

x n = (x 1 , • • • , x n ) ∈ X n for a sequence of locations with corresponding function values z n = (z 1 , • • • , z n ) T ∈ R n . Then, we introduce the set C R, n = C 1 × • • • × C n ⊂ R n of relaxed constraints, where C i = R j if z i ∈ R j for some j, C i = {z i } otherwise. ( 10
) Let also H R, n = {h ∈ H (X) | h(x n ) ∈ 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 minimize h H (X) subject to h ∈ H R, n (11)
has a unique solution given by s z n , where z n is the unique solution of the quadratic problem

arg min z∈C R, n z T K -1 n z. ( 12
)

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.

Proposition 4 Let ξ ∼ GP(0, k), x n = (x 1 , • • • , x n ) ∈ X n , z n ∈ R n and x m = (x 1 , • • • , x m ) ∈ X m
be a set of locations of interest where predictions should be made. Write Z n = (ξ (x 1 ), . . . , ξ (x n )) T and Z m = (ξ (x 1 ), . . . , ξ (x m )) T . Then the mode of the probability density function

p Z m , Z n | Z n ∈ C R, n (13)
is given by s z n (x m ), 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 Z n ∈ 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., [START_REF] Veiga | Gaussian process modeling with inequality constraints[END_REF][START_REF] Maatouk | Gaussian process emulators for computer experiments with inequality constraints[END_REF][START_REF] López-Lopera | Finite-dimensional Gaussian approximation with linear inequality constraints[END_REF], when constraints come from expert knowledge.

However, the predictive distribution ( 13) is non-Gaussian since the support of 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 µ and k) Given x n ∈ X n , z n ∈ R n , and a relaxation set R (finite union of closed intervals), the relaxed-GP (reGP) predictive distribution with fixed mean function µ and covariance function k is defined as the (Gaussian) conditional distribution of ξ ∼ GP(µ, k) given Z n = z n , where z n is given by

z n = arg min z∈C R, n (z -µ(x n )) T K -1 n (z -µ(x n )) , (14)
with C R, n defined by (10).

Observe that (14) reduces to (12) when µ = 0. Consequently, the mean of the distribution is the predictor s z n from Proposition 3 in this particular case, and is equal to µ + s z n in general. Moreover, the reGP predictive distribution can be seen as an approximation of (13), where p (Z n | Z n ∈ 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 {(µ θ , k θ ) , θ ∈ Θ}. 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 ∈ X n , z n ∈ R n , a relaxation set R (finite union of closed intervals), and parametric families (µ θ ) and (k θ ) as in Section 2.1, the relaxed-GP (reGP) predictive distribution with estimated parameters is the (Gaussian) conditional distribution of ξ ∼ GP(µ θ , k θ ) given Z n = z n , where z n and θ = θ n are obtained jointly by minimizing the negative log-likelihood:

θ n , z n = arg min θ ∈Θ, z∈C R, n L (θ ; z) , (15)
with C R, n defined by (10).

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

L (θ ; Z n ) = -ln p Z n, 0 | θ -ln p Z n, 1 | θ , Z n, 0 , (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 Z n, 0 . (Note also that θ 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 θ 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.

Remark 8 (Numerical details) Minimizing (15) with respect to z falls under the scope of quadratic programming (see, e.g., [START_REF] Nocedal | Quadratic programming. Numerical optimization[END_REF] 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 [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] using the gradient of (15).

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.

KNOWN CONVERGENCE RESULTS ABOUT INTERPOLATION IN RKHS

Recall that the fractional-order Sobolev space W β 2 (R d ), with regularity β ≥ 0, is the space of functions on R d defined by

W β 2 (R d ) = h ∈ L 2 (R d ), h W β 2 (R d ) = h 1 + • 2 β /2 L 2 (R d ) < +∞ , where h ∈ L 2 (R d ) is the Fourier transform of h ∈ L 2 (R d ). For a given X ⊂ R d , define the Sobolev spaces W β 2 (X) = h |X , h ∈ W β 2 (R d ) endowed with the norm h W β 2 (X) = inf g∈W β 2 (R d ), g |X =h g W β 2 (R d ) . (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., [START_REF] Adams | Sobolev spaces[END_REF], Section 4.9), and is equal to the closure of its interior.

Assumption 1 ensures that the previous definition coincides with other commons 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 × X → R is said to have regularity α > 0 if the associated RKHS H (X) coincides with W α+d/2 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 + • 2 -ν-d/2 ≤ r ≤ C 2 1 + • 2 -ν-d/2
for some C 2 ≥ C 1 > 0, and have therefore Sobolev regularity α = ν on 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 [START_REF] Wendland | Scattered data approximation[END_REF], 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 × X → R, and let (x n ) n≥1 ∈ X N be a sequence of distinct points. The following property (a minor reformulation of Theorem 4.1 of [START_REF] Arcangéli | An extension of a bound for functions in Sobolev spaces, with applications to (m, s)-spline interpolation and smoothing[END_REF] gives error bounds that depend on the Sobolev regularity of k and the so-called fill distance of x n ∈ X n , defined by

h n = sup x∈X min 1≤i≤n x -x i . ( 18
)
Proposition 9 Let k be a reproducing kernel with regularity α > 0. If X verifies Assumption 1, then

sup x∈X σ n (x) h α n , n ≥ 1 , ( 19
)
where denotes inequality up to a constant, that does not depend on (x n ) n≥1 .

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

Corollary 10 Let k be a reproducing kernel with regularity α > 0, H (X) the RKHS generated by k, and let f ∈ H (X). As above, let s z n be the solution of (8

) for z n = ( f (x 1 ), • • • , f (x n )) T , n ≥ 1. If X verifies Assumption 1, then f -s z n L ∞ (X) h α n f H (X) . (20)

CONVERGENCE RESULTS FOR REGP

Let k : X × X → R be a continuous strictly positive-definite reproducing kernel. In this section, we consider the zero-mean reGP predictive distribution obtained from ξ ∼ GP(0, k), with relaxed interpolation constraints on a union R = q j=1 R j of disjoints closed intervals R j with non-zero length. Let H (X) be the RKHS attached to k, f ∈ H (X), and consider a sequence (x n ) n≥1 ∈ X N of distinct points. Furthermore, define (the) regions X j = {x ∈ X, f (x) ∈ R j } for 1 ≤ j ≤ q and X 0 = X \ j≥1 X j . We give results about the limit of the sequence of reGP predictive distributions that suggest an improved fit in X 0 .

Let s R, n = s z n be the relaxed predictor from Proposition 3 based on

(x 1 , • • • , x n ) and ( f (x 1 ), • • • , f (x n )) T , n ≥ 1. The following proposition establishes the limit behavior of (s R, n ) n≥1 . Proposition 11 Let U ⊂ X and let H R, U denote the set of functions h ∈ H (X) such that, for all x ∈ U, h(x) ∈ R j if f (x) ∈ R j for some j, h(x) = f (x) otherwise. ( 21
)
Then the problem

min h∈H R, U h H (X) (22)
has a unique solution denoted by s R,U . Moreover, s R,n

H (X) ----→ s R,U , (23)
with U the closure of {x n }.

In particular, when {x n } is dense in X, then U = X and (s R, n ) n≥1 converges to s R, X , which is the minimal-norm element of the set H R, X .

The next proposition tells us that the interpolation error on X 0 can be bounded by a term that depends on the norm of s R, X .

Proposition 12 For any x ∈ X 0 and n ≥ 1,

| f (x) -s R, n (x)| ≤ 2σ n, 0 (x) s R, X H (X) , ( 24
)
where σ n, 0 is the power function obtained using only points in 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 α > 0. Let B ⊂ X 0 verify Assumption 1. Then, for all n ≥ 1,

f -s R, n L ∞ (B) h α n s R, X H (X) . ( 25
)
Let d(y, A) be the distance of y ∈ R to A ⊂ R. For j ≥ 1, x ∈ X j , and for all n ≥ 1:

d(s R, n (x), R j ) h α n s R, X H (X) if α < 1, (26)
d(s R, n (x), R j ) (|ln(h n )| + 1) h n s R, X H (X) if α = 1, (27) and d
(s R, n (x), R j ) h n s R, X H (X) if α > 1 , ( 28
)
where 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, X H (X) ?

Proposition 14 Suppose that k has regularity α > 0 and that there exists some j ≥ 1 such that X j has a non-empty interior. We have

sup h∈H R, X h H (X) = +∞, (29)
with H R, X given by (21) for f ∈ H (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, +∞) can be derived, and shows that

sup h∈H R, X h L ∞ (X) = +∞ .
Overall, no matter the element of H R, X at hand, reGP converges to a function s R, X which: coincides with f on X 0 , verifies f (x) ∈ R j ⇔ s R, X (x) ∈ R j for all x ∈ X, and is "nicer" than f in the sense of • H (X) . Furthermore, reGP yields error bounds carrying the norm of s R, X , which can be arbitrarily smaller than the norm of f in the case of a finite-smoothness covariance function.

Remark 15 Note that σ n ≤ σ n, 0 due to the projection residuals interpretation. Empirical and theoretical results about the screening effect (see, e.g., [START_REF] Stein | Rietz lecture: When does the screening effect hold?[END_REF][START_REF] Bao | Screening effect in isotropic Gaussian processes[END_REF], suggests that σ n σ n, 0 , if k has smoothness α > 0. In this case, observe that-no matter the element of 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 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 ν → ∞).

Choice of the Relaxation Set

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 ⊂ 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 R \ Q, on which interpolation constraints should be relaxed.

In the following, we put R (0) = R \ 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 θ of the underlying GP ξ . 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 [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF]. Using the formalism of scoring rules (see, e.g., [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF][START_REF] Petit | Model parameters in Gaussian process interpolation: an empirical study of selection criteria[END_REF], selecting parameters by a leave-one-out approach amounts to minimizing a selection criterion written as

J n (R) = 1 n n ∑ i=1 S (P R, n, -i , f (x i )) , (30)
where P R, n, -i is reGP predictive distribution with data z n,-i = (z 1 , . . . , z i-1 , z i+1 , . . . , z n ) and relaxation set R. The function S in (30) is a scoring rule, that is, a function S :

P × R → R ∪ {-∞, +∞},
acting on a class P of probability distributions on R, such that S(P, z) assigns a loss for choosing a predictive distribution P ∈ P, while observing z ∈ 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 ) ∈ R (0) would make it impossible to assess if the model is good at predicting that f (x) ∈ R (0) for a given x ∈ X. For instance, in the case of minimization, with Q = (-∞, t (0) ) and R (0) = [t (0) , ∞), 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.

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 [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF][START_REF] Matheson | Scoring rules for continuous probability distributions[END_REF], which may be written as

S (P, z) = +∞ -∞ s(F P (u), 1 z≤u ) µ(du) , (31)
where s : [0, 1] × {0, 1} → R ∪ {-∞, +∞} is a scoring rule for binary predictors, and µ is a Borel measure on R. A well-known instance of ( 31) is the continuous ranked probability score [START_REF] Gneiting | Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation[END_REF] written as

S CRPS (P, z) = +∞ -∞ (F P (u) -1 z≤u ) 2 du ,
which is obtained by choosing the Brier score for s and the Lebesgue measure for µ.

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

S tCRPS Q (P, z) = Q (F P (u) -1 z≤u ) 2 du. ( 32
)
This scoring rule, proposed by [START_REF] Lerch | Comparison of non-homogeneous regression models for probabilistic wind speed forecasting[END_REF] in a different context, reduces to S CRPS when Q = R. It can be seen as a special case of the weighted CRPS [START_REF] Matheson | Scoring rules for continuous probability distributions[END_REF][START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF][START_REF] Gneiting | Comparing density forecasts using threshold-and quantile-weighted scoring rules[END_REF], in which the indicator function 1 Q plays the role of the weight function-in other words, the measure µ in (31) has density 1 Q with respect to Lebesgue's measure. Consider for instance the case Q = -∞, t (0) :

S tCRPS Q (P, z) = t (0) -∞ (F P (u) -1 z≤u ) 2 du. ( 33
)
The upper endpoint t (0) of the range will be referred to as the validation threshold. Note that, in this case, S tCRPS Q (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.

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: 0) . The relaxation set R corresponding to the region above t has been obtained by the procedure described in Section 4.3.

J n (R) = 1 n n ∑ i=1 S tCRPS Q (P R, n, -i , f (x i )) . (34)
Using (34), we suggest the following procedure to select a reGP model. First, choose a sequence of nested candidate relaxation sets

R (0) ⊃ R (1) ⊃ • • • ⊃ R (G) = /
0. The next step is the computation of J n (R (g) ), g = 0, . . . , 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 ( θ

(g) n,-i , z (g)
n,-i ) 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 [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF] for Gaussian processes: for each set R (g) , solve (15) to obtain θ

(g) n and z (g) n = (z (g) 1 , . . . , z (g)
n ) T , and then compute the conditional distributions ξ (x i ) | {ξ (x j ) = z (g) j , j = i}, where ξ ∼ GP(µ, k), and where µ and k have parameter θ (g) n , using the fast leave-one-out formulas. By doing so, we neglect the difference between θ (g) n,-i and θ (g) n and the difference between z (g) n,-i and the vector (z

(g) 1 , . . . , z (g) i-1 , z (g) i+1 , . . . , z (g) n ) T .
The procedure ends by choosing the relaxation set R (g) that achieves the best LOO-tCRPS value. illustration but for the function c α 6 , with α = 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.

An Example for the Estimation of an Excursion Set

We illustrate the method on the problem of estimating an excursion set {x ∈ X, f (x) ≤ 0}. We consider the G10 optimization problem used by [START_REF] Regis | Constrained optimization by radial basis function interpolation for high-dimensional expensive black-box problems with infeasible initial points[END_REF], and focus on the constraint c 6 ≤ 0. Finding solutions satisfying the c 6 ≤ 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, [START_REF] Feliot | A Bayesian approach to constrained single-and multi-objective optimization[END_REF] found that the difficulty could be overcome by performing an ad-hoc monotonic transformation z → z α , with α = 7, on the constraint.

The estimation of an excursion set { f ≤ 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) = (-∞, -t (g) ] ∪ [t (g) , +∞) with a sequence of thresholds t (0) < • • • < t (G) = +∞, 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 → z α proposed by [START_REF] Feliot | A Bayesian approach to constrained single-and multi-objective optimization[END_REF]. If we apply the reGP framework on the transformed function c α 6 (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.

Application to Bayesian Optimization

-t (g) 0 +t (g) z n -t (g) 0 +t (g) z n R R Q R R Q Figure 7
: A reGP fit of c 6 , where the relaxation thresholds have been selected by LOO-tCRPS.

The observations 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 ± t (g) , with t (g) chosen to be t (0) by the LOO-tCRPS. Finally, the blue line shows a best fit by z → z 7 .

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∈X f (x) , ( 35
)
where f is a real-valued function defined on a compact set X ⊂ 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 . . . ∈ 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 , +∞). More precisely, the sequence (X n ) is constructed sequentially using the rule

X n+1 = arg max x∈X E n (m n -ξ (x)) + , ( 36
)
where

m n = f (X 1 ) ∧ • • • ∧ 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 ), . . . , 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 (0) n < t (1) n < • • • < t (G) n , ( 37
)
where t

(0)
n is the validation threshold, which delimits the range of interest Q n = -∞,t (0) n 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 [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF].

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.

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 X is a compact domain and that k : X × X → R is continuous, strictly positivedefinite, 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 σ 2 n (x) at x ∈ 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 α > 0 verifies the NEB property (Vazquez and Bect, 2010b), whereas the squared-exponential covariance function does not (Vazquez and Bect, 2010a).

Proposition 17 Let k : X × X → R be a continuous strictly positive-definite covariance function that verifies the NEB property, H (X) the corresponding RKHS and f ∈ H (X). Let n 0 > 0. Let (X n ) n≥1 be a sequence in X such that, for each n ≥ n 0 , X n+1 is obtained by (36) with t n > m n . Then the sequence (X n ) n≥1 is dense in X.

Proposition 17 implies the convergence of EGO-R with a fixed threshold t > min i≤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 , +∞), 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 (0) n > m n . In this case, the norm term in (24) gets smaller if (t n ) n≥1 is decreasing.

Optimization Benchmark

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

METHODOLOGY

In practice, we must choose the sequence of thresholds (37). The validation threshold t (0) n 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 (0) n 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 (0) n corresponds to an empirical quantile of z n . In the case of the constant heuristic, we set t (0) n to the α-quantile of the function values on an initial design, which is typically built to fill X as evenly as possible with, e.g., maximin Latin hypercube sampling [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. The numerical experiments were conducted with α = 0.25 in this article.

In the case of the concentration heuristic, we consider an α-quantile of the values of f at the points visited by the algorithm (again with α = 0.25). As the optimization algorithm makes progress, the evaluations will likely concentrate around the global minimum. Thus, t (0) n will get closer to the minimum value, and the range

Q n = (-∞, t (0)
n ) 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 , the candidate relaxation thresholds t

(g) n , g = 1, . . . , G, are chosen so that t (g) n -m n ranges logarithmically from t (0) n -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 (0) n , 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/2.
The maximization of the sampling criteria ( 6) and ( 36) is carried out using a sequential Monte Carlo approach [START_REF] Benassi | Bayesian optimization using sequential Monte Carlo[END_REF][START_REF] Feliot | A Bayesian approach to constrained single-and multi-objective optimization[END_REF].

For reference, we also run the Dual simulated Annealing algorithm (inspired by [START_REF] Xiang | Generalized simulated annealing algorithm and its application to the Thomson model[END_REF]) from SciPy [START_REF] Virtanen | SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF], with the default settings and with a random initialization.

The optimization algorithms are tested against a benchmark of test functions from [START_REF] Surjanovic | Virtual library of simulation experiments: Test functions and datasets[END_REF] 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 [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] and [START_REF] Merrill | An empirical study of Bayesian optimization: Acquisition versus partition[END_REF]. In particular, we also use a log-version of the Goldstein-Price function as [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF].

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., [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]. 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 ).

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.

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 (0) n 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.

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 -∞,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. More generally, the goal-oriented approach proposed in this article is not limited to singleobjective (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., [START_REF] Chevalier | Fast parallel krigingbased stepwise uncertainty reduction with application to the identification of an excursion set[END_REF]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.

• for q ≥ 2, we have EI ↑ q (P, z) = σ h q µ-z σ , where h q (t) = qt Φ q (δ t q ; 0, D q D T q ) + q Φ (-t) q-1 φ (t) + q(q -1)

2 √ π Φ q-1 (δ t q-1 ; 0, 1 2 B q ), (42)
where Φ q (•; m, Σ) is the cdf of the multivariate N (m, Σ) distribution, B q = 2 diag(0, 1 T q-2 ) + 1 q-1 1 T q-1 , D q is the matrix representing the linear map R q → R q , (y 1 , . . . , y q ) T → (-y 1 , y 2y 1 , y 3y 1 , . . . , y qy 1 ) T , and δ t q = (t, 0 T q-1 ) T ,

• finally for q = 2 we have

E (N 1 ∨ N 2 ) = µ + σ √ π . ( 43
)
For a scoring rule S : P × R → R and P 1 , P 2 ∈ P such that y ∈ R → S(P 1 , y) is P 2 -integrable, write S(P 1 , P 2 ) = E U∼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. [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF] A scoring rule S : P × R → R is said to be (strictly) proper with respect to P if, for all P 1 , P 2 ∈ P, the mapping y ∈ R → S(P 1 , y) is P 2integrable and the mapping P 1 ∈ P → S(P 1 , P 2 ) admits P 2 as a (unique) minimizer.

A strictly proper scoring rule S on a class P induces a divergence (P 1 , P 2 ) → S(P 1 , P 2 ) -S(P 2 , P 2 ), which is non-negative on P × P, and vanishes if and only if P 1 = P 2 . In the case of the truncated CRPS, simple calculations lead to [START_REF] Matheson | Scoring rules for continuous probability distributions[END_REF]:

S tCRPS Q (P 1 , P 2 ) -S tCRPS Q (P 2 , P 2 ) = Q (F P 1 (u) -F P 2 (u)) 2 du.
It follows that S tCRPS Q is proper for any measurable Q ⊂ R, and is strictly proper with respect to the class of non-degenerate Gaussian measures on R as soon as Q has non-empty interior.

Appendix B. Proofs

Lemma 24 [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], Section 1.5) Let k : X × X → R be a positive-definite covariance function, U ⊂ X, and H (U) be the RKHS attached to the restriction of k to U × U. The RKHS H (U) is the space of restrictions of functions from H (X) and the norm of g ∈ H (U) is given by

inf h∈H (X), h |U =g h H (X) . (44)
Proof of Proposition 3. First the existence and the uniqueness of the solution are given by the first statement of Proposition 11 (with

H R, n = H R, {x 1 ,••• ,x n } ).
Furthermore let z ∈ R and write α = K -1 n z, the reproducing property (7) gives

s z 2 H (X) = α T K n α = z T K -1 n z, (45)
and therefore

min h∈H R, n h 2 H (X) = inf z∈C R, n min h∈H (X), h(x n )=z h 2 H (X) = inf z∈C R, n z T K -1 n 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 Z m T , Z T n T . Using the equalities ( 5) and ( 45), and a slight abuse of notation by dropping irrelevant constants with respect to z and z ∈ C R, n , we have

-2 ln p z , z | Z n ∈ C R, n = z T , z T K -1 m, n z T , z T T = min h∈H (X), h(x n )=z, h(x m )=z h 2 H (X) .

This gives inf

z ∈R m , z∈C R, n -2 ln p z , z | Z n ∈ C R, n = min h∈H (X), h(x n )∈C R, n h 2 H (X) ,
which is reached by taking z = z n and z = s z n (x 1 ), . . . , s z n (x m ) T .

Proof of Proposition 9. First, one has

sup x∈X σ n (x) = sup x∈X sup f H (X) =1 | f (x) -s z n (x)| = sup f H (X) =1 f -s z n L ∞ (X) . Now, let f ∈ H (X) such that f 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 α+d/2 2 (X o ) defined by ( 17) is normequivalent to the Sobolev-Slobodeckij space (see, e.g., Di Nezza et al. (2012, Proposition 3.4) for a statement on 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 fs 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≥1 ), we have:

f -s z n L ∞ (X) = f -s z n L ∞ (X o ) h α n f -s z n W α 2 (X o ) h α n f -s z n H (X) h α n (46) by continuity of f -s z n , since • W α+d/2 2 (X o ) ≤ • W α+d/2 2 (X) due to the definition (17), W α+d/2 2 (X)
being norm equivalent to 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 ≤ h 0 for simplicity by increasing the constant eventually, since sup x∈X σ n (x) ≤ sup x∈X k(x, x) < +∞ by compacity.

Proof of Proposition 11. First observe that H R, U is not empty since it contains f . Furthermore, it is easy to verify that 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 H R, n . Moreover, the sequence (H R, n ) n≥1 is non-increasing, so the sequence (s R, n ) n≥1 converges in H (X) to the projection of the null function on n≥1 H R, n (see, e.g., Brezis, 2011, Exercice 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 0 n and z 0 n to be data points within X 0 , and s x 0 n , z 0 n the associated (interpolation) predictor, i.e. the solution of (8). Observing that s x 0 n , z 0 n interpolates s R, n , we have:

f (x) -s R, n (x) ≤ f (x) -s x 0 n , z 0 n (x) + s x 0 n , z 0 n (x) -s R, n (x) ≤ σ n, 0 (x) f H (X 0 ) + σ n, 0 (x) s R, n H (X 0 ) ≤ 2σ n,0 (x) s R, X H (X) ,
since f coincides with s R, X on X 0 , s R, X ∈ H R, n , and by Lemma 24.

Lemma 25 If k has smoothness α > 0, then there exists h 0 > 0 depending only on α such that, for all x, y ∈ X verifying xy ≤ h 0 , we have

δ y -δ x H (X)
xy α , for α < 1, (47)

δ y -δ x H (X) |ln( x -y )| x -y , for α = 1, (48)
and

δ y -δ x H (X) x -y , for α > 1. ( 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 α.

In this case, we have

δ y -δ x 2 H (X) = k(x, x) + k(y, y) -2k(x, y) = 2 (1 -r α ( x -y )) ,
with r α the corresponding isotropic correlation function. 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 ρ > 0 and angle φ ∈ (0, π/2). If X n, B is not empty, then the compacity of B ensures the existence of an x ∈ 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 ∈ B.) A cone C originating from x with angle φ and radius δ = min(h n, B , ρ) is contained in B and its interior do not contains observations. Furthermore, Wendland (2004, Lemma 3.7) shows that there exists a y ∈ C such that the open ball B y, δ sin(φ ) (1 + sin(φ )) -1 is subset of C, and therefore contains no observations as well. Thus, h n ≥ δ sin(φ )(1 + sin(φ )) -1 . Now, if h n, B ≤ ρ, then the desired result follows. If not, the result holds as well since h n, B ≤ diam(B).

Proof of Proposition 13. For the first assertion, let σ n, B be the power function using only the observations within B. Using Proposition 12, the inequality σ n, 0 ≤ σ 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 , • • • , x n } ∩ 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 ≥ 1. In addition, they are disjoint so δ = min 1≤ j<p inf x∈X j , y∈X p

xy > 0.

Suppose now that h n < δ and let j ≥ 1, x ∈ X j and 1 ≤ i ≤ n the index of the closet x i to x. By definition, xx i ≤ h n and therefore x ∈ X 0 ∪ X j . Let H (X) be the (topological) dual of H (X) and δ y : h ∈ H (X) → h(y), which lies in H (X) for all y ∈ X. Then using the reproducing property (7), we have |s R, n (x i )s R, n (x)| ≤ δ x iδ x H (X) s R, n H (X) ≤ δ x iδ x H (X) s R, X H (X) , and therefore

d(s R, n (x), R j ) ≤ δ x i -δ x H (X) s R, X H (X) + d(s R, n (x i ), R j ).
Now, if x i ∈ X j , then d(s R, n (x i ), R j ) = 0. Otherwise, x i ∈ X 0 necessarily and then, using the fact that s R, X (x) ∈ R j , we have:

d(s R, n (x i ), R j ) ≤ |s R, n (x i ) -s R, X (x)| = |s R, X (x i ) -s R, X (x)| ≤ δ x i -δ x H (X) s R, X H (X) .
So one can use Lemma 25 along with the previous statements to conclude if h n < min (δ , h 0 ). Finally, treating the case h n ≥ min (δ , h 0 ) is straightforward using the fact that sup x∈X k (x, x) is finite thanks to the compacity of X and d(s R, n (x), R j ) ≤ |s R, n (x)-s R, X (x)| for j ≥ 1 and x ∈ X j .

Lemma 27 If g, h ∈ W Lemma 31 Let a, b ∈ R with a ≤ b. Let q ≥ 1. Then R q (a, b) = ba + EI ↑ q (P, b) -EI ↑ q (P, a).

(54)

Proof.

R q (a, b) = 1 a≤u≤b F P (u) q du = 1 a≤u≤b q ∏ j=1 E 1 N j ≤u du with N j

iid

∼ P = E 1 a∨N 1 ∨•••∨N q ≤u≤b du = E (b -a ∨ N 1 ∨ • • • ∨ N q ) +
= ba + EI ↑ q (P, b) -EI ↑ q (P, a).

.

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 → -∞

EI ↑ q (P, a) = E (N 1 ∨ • • • ∨ N q ) -a + o(1), (55)
and therefore

R q (-∞, b) = lim a→-∞ R q (a, b) = b + EI ↑ q (P, b) -E (N 1 ∨ • • • ∨ N q ) , (56)
which gives the second statement. Finally, a change of variable gives (F P (u) -P(U = -u) -1 -z<u ) 2 du, and the last statement follows by observing that a probability measure admits at most a countable number of atoms.

  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.

Figure 2 :

 2 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 µ n given by the model and the gray envelopes represent the associated uncertainties.

Figure 4 :

 4 Figure 4: An example of reGP predictive distribution with R = (-∞, -1] ∪ [1, +∞) on a function f represented in dashed black lines. The solid black lines represent the relaxation thresholds. The problem (15) was solved only in z as the parameters of the (constant) mean and (ν = 5/2 Matérn) covariance functions were held fixed for illustration purposes.

Figure 5 :

 5 Figure 5: Illustration of the choice of a relaxation set. The range of interest Q is determined by the threshold t (0) . The relaxation set R corresponding to the region above t has been obtained by the procedure described in Section 4.3.

Figure 5

 5 Figure 5 illustrates the selection of the relaxation set used in Figure 3.

Figure 6 :

 6 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 , with α = 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.

Figure 8 :Figure 9 :

 89 Figure8: 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.

  γ 2 (X) for γ > d/2, then gh ∈ W γ 2 (X).

1

  a≤u≤b 1 z≤u + F P (u) 2 -2 1 z≤u F P (u) du = +∞ -∞ 1 a∨z≤u≤b + 1 a≤u≤b F P (u) 2 -2 1 a∨z≤u≤b F P (u) du = (ba ∨ z) + + R 2 (a, b) -2R 1 (a ∨ z, b).

Figure 10 :

 10 Figure 10: Same as Figure 8. The gray dashed line stands for the Dual simulated Annealing al-gorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the "Constant" and "Concentration" heuristics.

Figure 11 :Figure 13 :

 1113 Figure 11: Same as Figure 8. The gray dashed line stands for the Dual simulated Annealing al-gorithm, the red line for the standard EGO algorithm, and the blue and green lines for EGO-R, with the "Constant" and "Concentration" heuristics.

Table 1 :

 1 Optimization benchmark. The acronym G-P stands for Goldstein-Price.

	Problem	d
	Branin	2
	Six-hump Camel	2
	Three-hump Camel	2
	Hartman	3, 6
	Ackley	4, 6, 10
	Rosenbrock	4, 6, 10
	Shekel	5, 7, 10
	Goldstein-Price	2
	Log-Goldstein-Price	2
	Cross-in-Tray	2
	Beale	2
	Dixon-Price	4, 6, 10
	Perm	4, 6, 10
	Michalewicz	4, 6, 10
	Zakharov	4, 6, 10

  Lemma 26 Let B ⊂ X verifying Assumption 1 and h n, B be the fill distance of X n, B = {x 1 , • • • , x n } ∩ B within B, with the convention h n, B = diam(B) if X n, B is empty. Then, h n, B h n .

Standard results on principal irregular terms (see, e.g.,

Stein, 1999, Chapter 2.7)

give the results.

Lemma 28 yields lim inf s n-1 (x ψ(n) ) ≥ min( f (y ), t ∞ ) with ψ(n) = max{φ (k), φ (k) ≤ n},and the claim follows by extracting a φ -subsequence.

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), -∞ ≤ a < b ≤ +∞, 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 → S tCRPS Q (P, z) is σ -additive.

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

The following expressions hold for a general predictive distribution P.

Proposition 21 Suppose that P has a first order moment.

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

Now, leveraging well-known analytic expressions (see, e.g., [START_REF] Nadarajah | Exact distribution of the max/min of two Gaussian random variables[END_REF][START_REF] Chevalier | Fast computation of the multi-points expected improvement with applications in batch selection[END_REF], we have the following closed-form expressions in the Gaussian case.

Proposition 22 [START_REF] Nadarajah | Exact distribution of the max/min of two Gaussian random variables[END_REF][START_REF] Chevalier | Fast computation of the multi-points expected improvement with applications in batch selection[END_REF] Suppose that P = N (µ, σ 2 ) and let φ and Φ denote respectively the pdf and the cdf of the standard Gaussian distribution. Then

Proof. By the definition (17) of W γ 2 (X), the functions g and h can be extended into functions on R d , having their product in W γ 2 (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) ⊂ X j (with r > 0) be an open ball. There exists a

) as a function on X, and

(X) as a function on X, and Lemma 27 ensures that f n ∈ H (X). Moreover, it is easy to check that f n ∈ H R, X . Observe that the sequence ( f n ) n≥1 converges pointwise to a discontinuous function that lies thus outside H (X).

Suppose now that f n H (X) +∞ 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 H (X).

Lemma 28 Use the notations of Proposition 17 and let (y n ) n≥1 be a sequence in X. Assume that the sequence (y n ) n≥1 is convergent, denote by y its limit and assume that y is an adherent point of the set {x n }. Let t ∞ = lim inft n , then

In particular, we have

Proof. Suppose that y / ∈ {x n }. Then, let x φ (n) n≥1 be a subsequence converging to y and let

We proceed as in Proposition 13 to have:

thanks to Lemma 25 and the inequality s n H (X) ≤ f H (X) . Finally, have s n (x ψ(n) ) ≥ min f x ψ(n) , t n by construction, so lim inf s n (x ψ(n) ) ≥ min ( f (y ) , t ∞ ), which gives the second assertion thanks to (51). Observe that f (x ψ(n) ) < t n ultimately if f (y ) < t ∞ for the first assertion.

If y ∈ {x n }, then the result also follows similarly.

Lemma 29 Using the notations of Proposition 17 and writing v n = sup x∈X ρ n,t n (x), we have lim inf v n = 0.

Proof. This is an adaptation of Lemma 12 from (Vazquez and Bect, 2010b). Let y be a cluster point of (x n ) n≥1 and let x φ (n) n≥1 be a subsequence converging to y . We are going to prove that v φ (n)-1 → 0. Define

Then, Lemma 28 gives 1 lim inf s φ (n)-1 (x φ (n) ) ≥ min( f (y ), t ∞ ). Moreover we have

Moreover, one can use (Vazquez and Bect, 2010b, Proposition 10, (i)

since γ is non-decreasing with respect to its first argument and continuous.

Proof of Proposition 17. This is an adaptation of Theorem 6 from (Vazquez and Bect, 2010b). For f ∈ H (X), write (x n ) n≥1 for the corresponding sequence (X n ) n≥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 ∈ X, write ρ n,t n (x) = γ(m ns n (x), σ 2 n (x)) for the expected improvement under the reGP predictive distribution, with γ the function defined in Proposition 1.

Suppose that there exists some x 0 ∈ X such that σ 2 n (x 0 ) ≥ C 1 > 0. The sequence (m n ) n≥1 converges and the reproducing property (7) yields

by Proposition 1. But this yields a contradiction with Lemma 29, so the decreasing sequence σ 2 n n≥1 converges pointwise on X to zero. Proposition 10 from Vazquez and Bect (2010b) then implies that every x ∈ X is adherent to {x n }.

Lemma 30 Assume that b is finite, and that either a is finite too or

where

Appendix C. Optimization Benchmark Results

The results are provided in Figure 10, Figure 11, Figure 12, and Figure 13, for the other test functions 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 "Concentration" heuristic yields clear benefits for Zakharov (6), Dixon-Price (4), Dixon-Price (6), and Rosenbrock (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.