Tanguy Appriou
email: tanguy.appriou@emse.fr

Didier Rullière

David Gaudrie

High-dimensional Bayesian Optimization with a Combination of Kriging models

Keywords: Bayesian Optimization, Kriging, Gaussian Process Regression, High Dimension, Maximum Likelihood Estimation, Model Aggregation

In Bayesian optimization (BO), Kriging surrogate models reducing the number of function evaluations to reach the optimum. This method has been successfully applied to many real-world applications in low dimensions (less than 30 design parameters). However, in high dimension, building an accurate Kriging model is difficult, especially when the number of samples is limited as is the case when dealing with numerical simulators. This is due to the inner optimization of the Kriging length-scale hyperparameters which can lead to inaccurate models and impacts the performances of the optimization. In this paper, we introduce a new method for high-dimensional BO which bypasses the length-scales optimization by combining sub-models with random length-scales, and whose expression, obtained in closed-form, avoids any inner optimization. We also describe how to sample suitable length-scales for the sub-models using an entropy-based criterion, in order to avoid degenerated sub-models having either too large or too small length-scales. Finally, the variance of the combination being not directly available, we present a method to compute the prediction variance for any weighting method. We apply our combined Kriging model to high-dimensional BO for analytical test functions and for the design of an electric machine. We show that our method builds more accurate surrogate models than ordinary Kriging when the number of samples is small. This results in faster convergence for BO using the combination.

Introduction

In engineering design optimization, surrogate models (also called metamodels) are widely used to emulate black-box functions we want to optimize [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] because such functions are often obtained via a computationally expensive computer simulation (e.g. computational fluid dynamics or finite-elements solvers). This makes the optimization prohibitively expensive to perform with usual optimization methods, such as evolutionary algorithmsor gradient-based approaches, due to the large number of function evaluations required. Global Bayesian optimization (BO) strategies hinging on surrogate models were developed to solve this type of optimization problems. In particular, Efficient Global Optimization (EGO) introduced by [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] has been applied to a wide variety of design optimization problems (e.g. [START_REF] Candelieri | Bayesian optimization of pump operations in water distribution systems[END_REF][START_REF] Meliani | Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design[END_REF][START_REF] Picheny | A bayesian optimization approach to find nash equilibria[END_REF]. This algorithm uses a non-parametric class of statistical models called Kriging models [START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF] as a surrogate. Kriging models have proven to be effective in modelling black-box functions, not only for the approximation of numerical experiments [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF], but also in other fields such as machine learning [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], where it is known as Gaussian Process regression.

Although EGO has been applied successfully to a number of low-dimensional applications (dimension d < 20), one of its main drawback is that Kriging scales poorly to highdimensional cases. Originating from geostatistics [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF][START_REF] Matheron | Principles of geostatistics[END_REF] where it was used as an interpolator for spatial fields, Kriging was initially devised for problems with dimension limited to either d = 2 or 3. However in real-world design optimization, engineering designs may commonly be parametrized by 50 or more parameters [START_REF] Shan | Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions[END_REF][START_REF] Gaudrie | Modeling and optimization with gaussian processes in reduced eigenbases[END_REF]. In higher dimensions, Kriging suffers from the curse of dimensionality [START_REF] Bellman | Dynamic programming[END_REF] and building an accurate surrogate model is met with various setbacks. This, in turn, is problematic for design optimization since the convergence speed of EGO is related to the accuracy of the surrogate. Thus, using accurate models is instrumental in order to reduce the number of function evaluations to reach the optimum. One of the main challenges in Kriging is the inner optimization of the covariance length-scale hyperparameters. The latter regulate the decay of the correlation between observations when their distance increases. Usually Kriging models implement anisotropy by considering one length-scale per dimension. Estimating these hyperparameters correctly is essential to obtain a model with a good accuracy and they are typically determined by Maximum Likelihood Estimation (MLE), an extremely popular method for fitting Kriging models whose theoretical properties have been studied extensively (e.g. [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF][START_REF] Zhang | Inconsistent estimation and asymptotically equal interpolations in modelbased geostatistics[END_REF][START_REF] Van Der Vaart | Information rates of nonparametric gaussian process methods[END_REF][START_REF] Kaufman | The role of the range parameter for estimation and prediction in geostatistics[END_REF][START_REF] Karvonen | Maximum likelihood estimation and uncertainty quantification for gaussian process approximation of deterministic functions[END_REF]. However, most of these results are based on asymptotic considerations on the number of samples, and fewer works study practical cases where the number of samples is limited. For example, [START_REF] Karvonen | Maximum likelihood estimation in gaussian process regression is ill-posed[END_REF] showed that, in certain conditions, MLE is ill-posed in the sense that it leads to an infinite estimate of the length-scales. Practical MLE of the Kriging hyperparameters is difficult, especially for high-dimensional problems, in particular because the size of the search space grows exponentially with the dimension. Since the optimization is typically solved using gradient-based methods (e.g BFGS) with multi-start or evolutionary algorithms [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF], there is no better solution than increasing the number of iterations which might result in an increased computational effort as the cost of the likelihood and its gradient scale with O(n 3 ). However, most of the times in design optimization, the number of samples is very limited and the cost of the length-scale optimization remains negligible compared to the cost of obtaining the samples. In these cases where the number of samples is small, [START_REF] Ginsbourger | A note on the choice and the estimation of kriging models for the analysis of deterministic computer experiments[END_REF] and [START_REF] Mohammed | Over-fitting in model selection with gaussian process regression[END_REF] showed empirically that MLE can lead to very dispersed results, and [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF] that it can fail to recover the true hyperparameters.

While not as frequently used as MLE, other methods exist to determine the length-scales. For example, [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of gaussian processes with model misspecification[END_REF] suggested that cross-validation is more adapted than MLE when the model is ill-specified. [START_REF] Li | Analysis of computer experiments using penalized likelihood in gaussian kriging models[END_REF] and [START_REF] Yi | Penalized gaussian process regression and classification for high-dimensional nonlinear data[END_REF] proposed penalized versions of the likelihood to reduce the variance of the length-scale estimator, and [START_REF] Gu | Robust gaussian stochastic process emulation[END_REF] used maximum a posteriori estimation with a reference prior to obtain a robust estimator avoiding the setting of lower and upper bounds which can also be a tricky point of the length-scale optimization. However, all these length-scale estimation procedures also face difficulties in high-dimension.

High-dimensional Kriging and Bayesian optimization has recently gained attention, see Binois and Wycoff (2022) for a review. One classical approach is to reduce the problem's dimension by embedding the design space into a lower dimension space (e.g. [START_REF] Constantine | Active subspaces: Emerging ideas for dimension reduction in parameter studies[END_REF][START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by partial least squares dimension reduction[END_REF], and by building the Kriging model in this low-dimension space. Other methods consider simplifying hypotheses such as additive models [START_REF] Durrande | Additive Covariance kernels for high-dimensional Gaussian Process modeling[END_REF], where the function is assumed to decompose into a sum of one-dimensional components, enabling a sequential optimization of the length-scales. However, all these additional assumptions (low dimension representation or additive structure) are not necessarily satisfied in practice and may not generalize to any design engineering problem. [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF] introduced another method for high-dimensional Kriging based on a combination of Kriging sub-models, each one having fixed random length-scales, thus bypassing the difficult hyperparameter optimization, with no additional assumption on the black-box. However in this paper, the focus was on comparing different weighting methods for combining the sub-models, and on studying the influence of the number of sub-models.

In this paper, we apply this method to high-dimensional Bayesian optimization. Notation and main concepts of Kriging and Bayesian optimization are introduced in Section 2. Then, in Section 3, we present a motivating example illustrating how classical MLE fails to produce a good model when only few observations are available. In Section 4, we detail our new model for high-dimensional Bayesian optimization: in Section 4.1, we present a new entropy-based criterion to sample suitable length-scales for the sub-models to avoid degenerated cases with either too large or too small length-scales; in Section 4.2 we describe the weighting method to combine the sub-models; and in Section 4.3 and 4.4, we explain how we compute the prediction variance of the combination, required for Bayesian optimization, and which is not readily available as the correlation between the sub-models is unknown. Finally, Section 5 discusses the performances of our method for the optimization of two high-dimensional analytical test functions, and for one real-world application.

2 Kriging surrogate models and Bayesian optimization

Ordinary Kriging

Originally, Ordinary Kriging (OK) was developed for interpolation of spatial data by [START_REF] Matheron | Principles of geostatistics[END_REF] who named the method after the South African mining engineer D.G. Krige. The method was later applied to the approximation of computer experiments [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF]. This method is also known as Gaussian Process Regression [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. This section briefly covers the basics of the method and introduces the notation used throughout this paper. We denote by y : x ∈ X ⊂ R d → y(x) ∈ R the d-dimensional black-box function that we want to approximate by a surrogate model. We have n training points X = (x 1 , . . . , x n ) ⊤ where the function y is known, and we denote Y = (y(x 1 ), . . . , y(x n )) ⊤ the function values at these locations. In Kriging, y is assumed to be the realization of a Gaussian process (GP) on X : Y (.) ∼ GP µ(.), σ 2 k θ (., .) . µ(.) is the mean function of the GP which we consider constant in this paper: µ(.) = µ. This is a standard choice for the approximation of computer codes where no we have prior knowledge about the true function [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF][START_REF] Ginsbourger | A note on the choice and the estimation of kriging models for the analysis of deterministic computer experiments[END_REF] known as ordinary Kriging. k θ : X × X → [-1, 1] is the positive definite correlation function indexed by the hyperparameters θ ∈ R d , called the correlation length-scales vector (also range or scale parameters), with one length-scale per dimension of the input space. Finally, σ 2 ∈ R + is another hyperparameter calibrating the amplitude of the variance, and σ 2 k(., .) is called the covariance function or kernel. A stationary GP with a Matérn-class covariance function is often recommended [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF][START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. One example of which, used throughout this paper, is the radial Matérn 5/2 correlation defined as:

k θ (x, x ′ ) := 1 + √ 5 x -x ′ θ + 5 3 x -x ′ θ 2 exp - √ 5 x -x ′ θ , (1)
where x-x ′ θ is the scaled distance between two points x, x ′ ∈ X using component-wise

division: x-x ′ θ 2 := d ℓ=1 x (ℓ) -x ′(ℓ) θ (ℓ)

2

. This is a typical choice for design optimization [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF] when there is no prior information on the function. This is because GP trajectories with this correlation are twice differentiable [START_REF] Abrahamsen | A review of gaussian random fields and correlation functions[END_REF], which may be better suited than the Gaussian kernel (infinitely differentiable trajectories) or than the exponential correlation (trajectories not differentiable). Other covariance functions can be used if prior knowledge about the unknown function is available (e.g. cylindrical kernels), and, notice that even when the covariance is misspecified, a proper estimation of the hyperparameters can still yield a model with good predictive capacities [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of gaussian processes with model misspecification[END_REF].

The ordinary Kriging predictor is a linear combination of the observations which is obtained by conditioning the Gaussian process Y over D = (X, Y ):

ŷ(x) := E(Y (x)|D) = µ + k θ (x, X)k θ (X, X) -1 (Y -µ), (2)
where k(x, X) is the vector of correlations between the prediction point x and the sample points X, and k(X, X) is the n × n matrix of correlations between the components of X.

Note that this predictor does not depend on σ 2 . One of the main advantage of Kriging is that it not only provides a prediction, but also the prediction variance, interpreted as a model uncertainty. This is particularly important in Bayesian optimization where, alongside with the prediction, this model uncertainty helps exploring under-visited parts of the design space. The prediction variance is:

ŝ2 (x) := Var(Y (x)|D) = σ 2 k θ (x, x) -k θ (x, X)k θ (X, X) -1 k θ (X, x) . (3)
In the following, we will sometimes denote the correlation matrix as K θ := k θ (X, X).

Hyperparameter estimation

The estimation of the covariance hyperparameters of a Kriging model drastically affects its precision. As illustrated in Figure 1, an appropriate choice of length-scales yields much better accuracy than arbitrary length-scales. As often in parametric statistic models, the usual approach to estimate the hyperparameters is to use maximum likelihood estimation (MLE) [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], which consists in maximizing the marginal likelihood of the model:

L(σ, θ) := 1 (2π) n/2 det(σ 2 K θ ) 1/2 exp - 1 2σ 2 (Y -µ) ⊤ K -1 θ (Y -µ) . (4)
This is equivalent to minimizing -log(L(σ, θ)). Given a fixed θ, the maximum likelihood estimator for µ and σ 2 are:

μ = 1 ⊤ K -1 θ Y 1 ⊤ K -1 θ 1
, and σ2

M LE = (Y -μ) ⊤ K θ -1 (Y -μ) n , (5)
The length-scales θ are estimated by minimizing the concentrated log-likelihood obtained by injecting (5) into (4):

θMLE = arg min θ n 2 log σ2 M LE + 1 2 log (det(K θ )) . (6)
In this expression, the length-scales are involved in both terms σ2 M LE and det(K θ ) through the correlation matrix. The inner optimization ( 6) is solved numerically, typically using a gradient-based method (e.g BFGS) with multi-start as the gradient of the likelihood with respect to the length-scales can be in closed-form, or using evolutionary algorithms [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF].

Bayesian optimization

As presented in the previous section, Kriging gives a way to build a surrogate model to approximate a black-box function based on some observations. Bayesian Optimization (BO) aims at finding the global optimum of this function x * = arg min x y(x) in as few evaluations of the black-box as possible. One common Bayesian optimization framework is the Efficient Global Optimization (EGO) algorithm introduced by [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], summarized in Algorithm 1.

Algorithm 1 Efficient Global Optimization (EGO) algorithm

Create an initial design plan:

X = (x 1 , . . . , x n ) ⊤ .
Compute the associated values: Y = (y(x 1 ), . . . , y(x n )) ⊤ . Fit the Kriging model to the data D = (X, Y ): find μ, σ2 , and θ. repeat

x n+1 = arg max x EI(x). X ← X ∪ x n+1 . Y ← Y ∪ y(x n+1 ).
Re-estimate the hyperparameters μ, σ2 , and θ. until convergence or budget exhaust

The EGO algorithm begins by evaluating a small set of designs x i called the initial design of experiments (DoE), usually space-filling, and sequentially adds new observations by maximizing a so-called acquisition criterion which quantifies the worth of any unevaluated

x. The most popular one is the Expected Improvement (EI), which is the expectation of the improvement I(x) = max(0, y min -Y (x)) over y min = min(Y ), the best value observed so far. For x ∈ X , the EI can be easily computed as:

EI(x) := E [I(x)] = (y min -ŷ(x))Φ y min -ŷ(x) ŝ(x) + ŝ(x)ϕ y min -ŷ(x) ŝ(x) , (7)
where which Φ(.) and ϕ(.) are respectively the cumulative distribution function and the density of a standard normal distribution. Note that the EI is computed using both the mean predictor ŷ and the prediction variance ŝ. It allows for a trade-off between exploitation and exploration by selecting new points near the current optimum (where ŷ is small) and far from any observation (where ŝ is large).

However, in high-dimensional problems, EGO suffers from some drawbacks related to the geometry of high-dimensional spaces. With a quick calculation, one can see that the volume of a thin layer of thickness ε in the border of the hypercube [0, 1] d tends to 1 when d grows to infinity. In other words, in high-dimensional spaces, almost all the volume is located on the boundary of the hypercube. Added the fact that the EI is a highly multi-modal function, its optimization can also be difficult and often leads to new points located on one face of the hypercube, resulting in a optimizer too explorative. To avoid this behavior, it is common practice in high-dimensional Bayesian optimization to use trust regions in order to reduce the size of the search space (see e.g. [START_REF] Eriksson | Scalable global optimization via local bayesian optimization[END_REF][START_REF] Diouane | Trego: a trust-region framework for efficient global optimization[END_REF][START_REF] Binois | A survey on high-dimensional gaussian process modeling with application to bayesian optimization[END_REF]. In these methods, the surrogate model and the optimization are restricted to a local neighborhood of the current best solution whose size increases if a new best solution is discovered, or decreases otherwise.

An example of Kriging failing in high dimension

In this section, we illustrate the inability of the classical Kriging approach to correctly estimate the length-scale hyperparameters for high-dimensional problems when not enough observations are available. While having too few observations can occur in any dimension, it is especially prevalent in higher dimensions where the thumb rule of 10d observations to build the model [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] is often not affordable. In the context of Bayesian optimization, it is also common practice to start with a small initial DoE (about 2-3d) [START_REF] Gaudrie | High-Dimensional Bayesian Multi-Objective Optimization[END_REF][START_REF] Garnett | Bayesian optimization[END_REF] and to allocate most of the budget for the acquisition points.

For this example, we consider two test functions for which we fit a Kriging model. First, the sphere test function defined as:

f sphere (x 1 , . . . , x d ) := d ℓ=1 (x ℓ -0.5) 2 , 0 ≤ x ℓ ≤ 1, ℓ = 1, . . . , d. ( 8
)
This function is simply a parabola centered in the middle of the hypercube [0, 1] d , and it is thus reasonable to expect a Kriging model to accurately approximate this simple smooth and convex function. Second, we consider trajectory samples of an isotropic GP:

Y iso (.) ∼ GP (0, k iso (., .)) , (9)
where k iso is a Matérn 5/2 covariance with isotropic length-scale θ true = 3 and amplitude σ 2 true = 1. These GP trajectories correspond to functions where the Kriging prior assumption is exactly satisfied. For the numerical experiment, we consider a dimension d = 50, and we build an anisotropic Ordinary Kriging model using a varying number of samples n ∈ [100, 1000]. The d hyperparameters are optimized using the DiceKriging R package [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF] by MLE using a maximum of 500 L-BFGS-B iterations and 3 restarts. We compare the estimated hyperparameters with "real" ones, the latter being obtained after fitting a "reference GP" to 5000 samples of the Sphere function. For the GP trajectories, these reference hyperparameters are simply those used for sampling the trajectory, θ true = 3. To measure the global accuracy of the models, we compute the Q 2 coefficient based on n test = 10000 random test points

x (t) 1 , . . . , x (t)
ntest ∈ [0, 1] d :

Q 2 := 1 - ntest k=1 ŷ(x (t) k ) -f (x (t) k ) 2 ntest k=1 f (x (t) k ) -1 ntest ntest l=1 f (x (t) l ) 2 , ( 10
)
in which f is either of the test functions.

The results for 10 different initial DoEs are shown in Figure 2. Figures 2a and2c show that the Ordinary Kriging fails in approximating the simple sphere function in the small data regime as highlighted by the very poor global accuracy, the Q 2 being even negative for less than 200 observations. This is especially problematic in early stages of Bayesian optimization since first surrogates with a poor global accuracy will negatively impact the discovery of areas of interest and decrease the convergence speed of the optimization.

The good performance of the reference hyperparameters (blue boxplots), able to achieve a very good accuracy with only 100 observations, shows that the issue does not lie in the modelling capability of Kriging surrogates but rather in the hyperparameter values. Moreover, the log-likelihood plots in Figure 2b and 2d inform us that the issue does not reside in the convergence of the 50 dimensional hyperparameter optimization as the loglikelihood of the estimated hyperparameters is superior to the log-likelihood of the reference hyperparameters as expected. It rather seems that with few points, maximum likelihood is no longer a relevant metric to select good hyperparameters, since high likelihoods do not correspond to accurate models.

Finally, although the experiment was performed for high-dimensional test functions (d = 50), the issue is not directly the dimension itself, but rather the low number of observations relative to the dimension. In fact, we retrieve in Figure 2 the usual empirical rule of 10d observations to obtain a precision close to that of the model with reference hyperparameters. As such, similar results can be obtained in lower dimension by scaling down the number of samples. However, while it seems reasonable to obtain more than 10d samples to build the model in lower dimension, it is more difficult in higher dimensions with limited observations. Therefore, the described lack-of-accuracy issue is more prevalent for high-dimensional problems.

To avoid this issue, a usual solution is to reduce the dimension of the problem. While this facilitates the hyperparameter optimization, it raises the question of the accurate representation of the true function in a lower dimension space and of the information loss incurred. In the next section, an alternative approach based on a combination of Kriging models with random length-scales is presented. This method keeps the information about correlations among all variables and generalizes to any problem since it does not require the low-dimensional representation hypothesis.

Sphere function

(a) Q 2 (b) Log-likelihood GP trajectories (c) Q 2 (d) Log-likelihood

Combination of Kriging models with random length-scales

To address the issues raised in the previous section, [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF] introduced a method to bypass the hyperparameter optimization by using a combination of Kriging models with random length-scales. The idea was to replace the difficult and sometimes costly length-scale estimation by the easier optimization of the weights in the combination. This results in an easier-to-build model, and which is more accurate in such cases where the length-scales are wrongly estimated by maximum likelihood estimation. In the aforementioned paper, the focus was on the weighting methods and on studying the impact of the number of sub-models. In the present paper, we develop a more effective methodology to sample the random sub-models length-scales, and the combination is further developed to accommodate to Bayesian optimization by providing a method to obtain the prediction variance.

The combined model M tot proposed in [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF] writes as:

M tot (x) := p i=1 w i (x)M i (x), (11)
where w i , i = 1, . . . , p, are the weights of the p sub-models M i . In the following sub-sections, we first discuss the choice of the sub-models, and in particular how to sample their random length-scales. Then, the second sub-section describes how the sub-models are weighted, and finally, the last two sub-sections detail how the prediction variance of the combined model is obtained, by modeling the correlation between sub-models and discussing the estimation of the variance amplitude hyperparameter.

Choice of the sub-models

The sub-models we consider are ordinary Kriging models with length-scales θ i chosen randomly, hence not as the result of a maximum likelihood estimation:

M i (x) := E(Y θ i (x)|D i ) = µ i + k θ i (x, X i )k θ i (X i , X i ) -1 (Y i -µ i ),
In the following, we will assume that the training data set of each sub-model is the entire data set:

D i = (X i , Y i ) = (X, Y ) though each sub-model could consider a specific subset of (X, Y ).
The sampling of the length-scales θ i is therefore critical as it is the only source of variability between the sub-models.

The sampling of length-scales should enable a wide diversity of sub-model behaviors so that the combination can select the appropriate ones through their weights w i . To this aim, the length-scales must be dissimilar enough, but at the same time, they need to lie in an appropriate range to avoid degenerated sub-models occurring when the lengthscales are either too large or too small compared to observed distances between samples.

In particular, small length-scales will result in all cross-correlations being close to 0. In this case, as illustrated in Figure 3a, the Kriging prediction equals to its mean function µ(.) almost everywhere with spikes to interpolate the observations. In the contrary, too large length-scales will result in all correlation being close to 1 resulting in a correlation matrix which can be ill-conditioned and inverting it -when possible -can cause numerical instabilities. This issue can be resolved by adding a small nugget term in the diagonal of K θ at the cost of loosing the interpolating property as illustrated in Figure 3b. The sampling procedure must not produce such non-informative sub-models. Note that this corresponds to only accepting robust sub-models as defined in [START_REF] Gu | Robust gaussian stochastic process emulation[END_REF]. In the original random length-scales combination paper [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF], to avoid the degenerated cases, the length-scales were sampled uniformly in a bounded interval computed according to their influence on the correlations. However, this method tended to produce too many lengthscales on the smaller side of the acceptable range which tends to worsen the sub-models. Instead of trying to design a more suited distribution to sample from, we propose an other non-parametric approach to better sample the length-scales based on the entropy of the correlation.

First, we will derive the entropy for a Gaussian correlation and detail the sampling strategy. Then we will extend the method to any other correlation. Let us consider that design points are distributed as a random vector X = (X (1) , . . . , X (d) ), with i.i.d components with common variance σ 2 X and kurtosis κ X . Let D 2 be the random squared distance between two independent points X and X ′ . For a large enough dimension d:

D 2 := d ℓ=1 X ℓ -X ′ ℓ 2 ∼ N 2dσ 2 X , 2dσ 4 X (κ X + 1) . ( 12
)
For a Gaussian correlation function, the random correlation between these two points for a given length-scale θ is then:

R θ := exp - D 2 2θ 2 ∼ log N - σ 2 X θ 2 d, σ 4 X 2θ 4 (κ X + 1)d . (13)
Finally, the entropy of the correlation isthe entropy of the log-normal distribution in (13):

H(R θ ) := E[-log f R θ (R θ )] = - σ 2 X θ 2 d + 1 2 log 2π σ 4 X 2θ 4 (κ X + 1)d + 1 2 , ( 14
)
where f R θ is the density of R θ . Figure 4 shows the entropy as a function of the length-scale θ for a dimension d = 50. The entropy decreases rapidly for length-scales below a critical value, while its drop is less pronounced for large length-scales. It is therefore especially important to avoid too small length-scales. Additionally, we can remark on the length-scale θ * which maximizes the entropy:

θ * := arg max θ H(R θ ) = σ X √ d.
This is the length-scale which maximizes the potential to capture the information of the training set, prior to any observation. We also note that this length-scale is such that the square of the length-scale is of the same order as the squared-distances between samples: [START_REF] Obrezanova | Gaussian processes: a method for automatic qsar modeling of adme properties[END_REF] used this definition as an initial value prior to any observation.

2θ * 2 = E[D 2 ].
The entropy obtained in equation ( 14) measures the variability of the correlations. The two degenerated cases described above (correlation is always 0 for small length-scales or always 1 for large length-scales) correspond to cases where there is no variability and where the entropy is small: For a general non-Gaussian correlation function, e.g. a Matérn one, the entropy cannot be obtained analytically. In this case, we use instead a non-parametric estimate of the entropy [START_REF] Ahmad | A nonparametric estimation of the entropy for absolutely continuous distributions (corresp.)[END_REF][START_REF] Beirlant | Nonparametric entropy estimation: An overview[END_REF] given by:

lim θ→0 H(R θ ) = -∞, and lim θ→∞ H(R θ ) = -∞.
Ĥ(R θ ) = - 1 n corr ncorr k=1 log fR θ k (k) θ , (15)
where n corr = n(n-1) 2 is the number of observed correlations k

(k)
θ , k = 1, . . . , n corr between distinct samples, and fR θ is a kernel density estimate of correlations density f R θ .

To obtain varied length-scales in an appropriate range for the sub-models, we sample length-scales corresponding to high entropies:

f (θ (ℓ) ) ∝ exp (H(R θ (ℓ) )) , ℓ = 1, . . . , d, (16)
where f (θ (ℓ) ) is the density from which θ (ℓ) is sampled. Note that each of the d lengthscales of the vector θ = (θ (1) , . . . , θ (d) ) are sampled independently. Note also that since the differential entropy can be negative (contrarily to the entropy for discrete random variables), a positive transformation of the entropy (exponential) is used in equation ( 16).

Weighting scheme

To obtain the weights of the sub-models, Appriou et al. ( 2023) compared different methods coming from the model combination literature (see e.g. [START_REF] Ginsbourger | Discrete mixtures of kernels for kriging-based optimization[END_REF][START_REF] Cao | Generalized product of experts for automatic and principled fusion of gaussian process predictions[END_REF][START_REF] Deisenroth | Distributed gaussian processes[END_REF][START_REF] Rullière | Nested kriging predictions for datasets with a large number of observations[END_REF]. In this paper, we follow the definition proposed in [START_REF] Acar | Ensemble of metamodels with optimized weight factors[END_REF] and Viana et al. (2009) which aim at minimizing the global mean-square error (MSE) of the combination given by:

E (M tot (X) -y(X)) 2 . ( 17
)
As the MSE is a global measure, we keep the weights constant throughout the design space: w(x) = w. Since we only have access to a finite number of observations, a discrete approximation of the MSE is obtained using leave-one-out cross-validation (LOOCV):

e LOOCV (M tot ) := 1 n n k=1 p i=1 w i M i -k (x k ) -y(x k ) 2 = w ⊤ Cw, (18)
where M i -k is the ith Kriging sub-model built by removing the kth sample x k . The components of the matrix C ∈ R p×p are c ij = 1 n e ⊤ i e j with e i the LOOCV vector for the ith sub-model: e i = (e (1) i , . . . , e (n) i ). For Kriging models, these residuals can be computed easily without the need to build n distinct models using the conditional Gaussian and block-matrix inversion formulas [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF][START_REF] Ginsbourger | Fast calculation of gaussian process multiple-fold cross-validation residuals and their covariances[END_REF]):

e (k) i = [K -1 θ i (Y -µ i )] k [K -1 θ i ] k,k . ( 19
)
The weights w are obtained by minimizing (18):

w * = arg min w w ⊤ Cw, subject to p i=1 w i = 1.
Using a Lagrange multiplier and setting the derivatives to zero:

w * = C -1 1 1 ⊤ C -1 1 . ( 20
)
One issue pointed out in Viana et al. (2009) and [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF] is that w can be negative or greater than one, which can lead to over-fitting, especially when the number of sub-models p is large. To solve this issue, instead of combining all sub-models at once, we combine them two-by-two, following a binary tree structure illustrated in Figure 5. With this definition, we can enforce w ∈ [0, 1] by imposing w = 0 when w < 0 and w = 1 when w > 1, while keeping the sum of weights equal to one. In practice, it is also possible to implement a threshold to induce more sparsity in the weights. In addition, this two-bytwo combination scheme will also prove itself convenient in the following subsection for computing the variance of the combination.

M tot M 12 M 1 w 1 M 2 w 2 w 12 M 34 M 3 w 3 M 4 w 4 w 34 . . . . . . . . .

Variance of the combination

As noted in section 2.3, the prediction error is essential in Bayesian optimization as it is used inside most acquisition functions. For Kriging models, it is obtained naturally in analytical form through the variance of the prediction as in equation ( 3). Each individual model in the combination is a Kriging model with Y i ∼ GP µ i , σ 2 k θ i (., .) and has a prediction variance:

ŝ2 i (x) = E (M i (x) -Y i (x)) 2 = σ 2 k θ i (x, x) -k θ i (x, X)k θ i (X, X) -1 k θ i (X, x) .
However, since the correlation structure between the sub-models is unspecified, we cannot access directly to the variance of the combination. In the literature, several works aim at estimating prediction errors for surrogate models where the latter is not readily available. For example, [START_REF] Viana | Efficient global optimization algorithm assisted by multiple surrogate techniques[END_REF] and Viana and Haftka (2009) import the uncertainty estimate of a Kriging sub-model to sub-models not equipped with a measure of uncertainty (e.g. RBFs) to perform Bayesian optimization. Den [START_REF] Hertog | The correct kriging variance estimated by bootstrapping[END_REF] use parametric bootstrapping instead of the classical formula to estimate the variance of a Kriging model.

Bootstrapping is also used in [START_REF] Kleijnen | Expected improvement in efficient global optimization through bootstrapped kriging[END_REF] and [START_REF] Kleijnen | Simulation-optimization via kriging and bootstrapping: a survey[END_REF] to perform EGO using this corrected error estimate. In the field of neural networks, various methods were also developed to estimate the uncertainty of the prediction (see for instance [START_REF] Papadopoulos | Confidence estimation methods for neural networks: A practical comparison[END_REF][START_REF] Khosravi | Lower upper bound estimation method for construction of neural network-based prediction intervals[END_REF][START_REF] Pearce | High-quality prediction intervals for deep learning: A distribution-free, ensembled approach[END_REF][START_REF] Abdar | A review of uncertainty quantification in deep learning: Techniques, applications and challenges[END_REF]. Conformal prediction [START_REF] Lei | Distributionfree predictive inference for regression[END_REF][START_REF] Romano | Conformalized quantile regression[END_REF] is an other empirical approach to build prediction interval without making distributional assumptions. [START_REF] Berk | Valid post-selection inference[END_REF] and [START_REF] Bachoc | Uniformly valid confidence intervals post-model-selection[END_REF] propose a procedure to construct valid confidence intervals post model selection for linear models. [START_REF] Acharki | Robust prediction interval estimation for gaussian processes by cross-validation method[END_REF] use cross-validation to calibrate prediction intervals for GPs in the case of model misspecification. Another related class of method are the mixture models (see [START_REF] Yuksel | Twenty years of mixture of experts[END_REF], for a review) which is a Bayesian approach to combining models in which the posterior predictive distribution of the combination is given as a mixture of posterior predictive distributions of several sub-models. In the context of Kriging, Pronzato and Rendas (2017) builds a fully-Bayesian mixture of submodels with different covariances to obtain a non-stationary model, and [START_REF] Ginsbourger | Discrete mixtures of kernels for kriging-based optimization[END_REF] uses a mixture of Kriging models in the context of Bayesian optimization. While the mixture models give the same mean prediction as the linear combination of models used in this paper, their variance differs. The variance of a mixture of Kriging models can be obtained analytically and the mixed expected improvement can be expressed as the convex combination of the individual expected improvements. However, the variance of a mixture is larger than that of a linear combination of models. In addition, the weights of the mixture can be interpreted as the probability to randomly select one model among all others, even in high-dimension where no individual sub-model is expected to be clearly better than all others. It is preferable to have a combination as it outperforms the best submodel as shown in [START_REF] Appriou | Combination of optimization-free kriging models for high-dimensional problems[END_REF]. In this section, we introduce a method to obtain the variance of the linear combination of Kriging models introduced in the previous section. This approach is based on a LOOCV strategy and can be seen as a natural extension of the linear combination of models hypothesis.

The variance of the combination depends on the global covariance structure between submodels. We make the hypothesis that the underlying Gaussian process Y tot is also a linear combination of independent Gaussian processes Y i , i = 1, . . . , p with the same fixed lengthscales as those of the M i :

Y tot = σ 2 tot p i=1 α i Y i , with Y i ∼ GP (µ tot , k θ i (., .)) , p i=1 α i = 1. ( 21
)
Here, σ 2 tot is the amplitude of the variance whose estimation is discussed in the next subsection. Without loss of generality, we consider centered GPs (µ tot = 0) as we are only interested in the associated covariance:

Cov(Y tot (x), Y tot (x ′ ) := σ 2 tot k tot (x, x ′ ) = σ 2 tot p i=1 α 2 i k θ i (x, x ′ ). ( 22
)
Note that the weights α i in the combination of GPs ( 21) are different from the weights w i of the combination of models in (11). Similarly to the combination of models in Figure 5, we also assume that the GPs are combined two-by-two as illustrated in Figure 6. This enables the derivation of an analytical expression for the weights α using a LOOCV procedure.

Notice the link between this approach and a mixture of models. Indeed, if α 2 = w, the MSE of a mixture equals the MSE of the combination:

E (M mix (x) -Y mix (x)) 2 = E (M tot (x) -Y tot (x)) 2 , when k tot (x, x ′ ) = p i=1 w i k θ i (x, x ′ ).
(23) Details on the proof of this proposition are given in Appendix A.1. To obtain the weights α, we aim at minimizing the discrepancy between the combination of sub-models that we obtained, and the new hypothesis we introduced about the combination of GPs. This discrepancy can be measured as the expected MSE of the combined model with respect to Y tot . For a combination of two models, with Y tot = αY 1 + (1 -α)Y 2 , the discrepancy is minimal for:

M tot M 12 M 1 w 1 M 2 w 2 w 12 M 34 M 3 w 3 M 4 w 4 w 34 . . . . . . . . . Y tot Y 12 Y 1 α 1 Y 2 α 2 α 12 Y 34 Y 3 α 3 Y 4 α 4 α 34 . . . . . . . . .
α * = arg min α E E (wM 1 (x) + (1 -w)M 2 (x) -αY 1 (x) -(1 -α)Y 2 (x)) 2 |Y 1 , Y 2 . (24)
Finally, by approximating the global MSE using the LOOCV error as in equation ( 18), we obtain an analytical expression for α * :

α * = a 1 (w) a 1 (w) + a 2 (w) , ( 25
)
where a 1 (w) and a 2 (w) depend on the weight w of the combination of models and on the LOOCV vector of the sub-models, and are given by:

a 1 (w) = w 2 E(e LOOCV (M 1 )|Y 2 ) + (1 -w 2 )E(e LOOCV (M 2 )|Y 2 ), (26)
a 2 (w) = (1 -w 2 )E(e LOOCV (M 1 )|Y 1 ) + (1 -(1 -w 2 ))E(e LOOCV (M 2 )|Y 1 ). ( 27
)
Here, e LOOCV (M i )|Y j , i, j ∈ {1, 2}, is the LOOCV error of the model M i when the global GP is Y j (.) ∼ GP(0, k θ j (., .)), which can be obtained with the Kriging LOO formula ( 19):

e LOOCV (M i )|Y j = n k=1 [K -1 θ i ] k, [K -1 θ i ] k,k Y j (X) 2 , with Y j (X) ∼ N (0, k θ j (X, X)). ( 28
)
Details of the proof are given in Appendix A.2.

Finally, the variance of the combination is the standard Kriging variance using the global covariance in ( 22) with the weights α computed with (25):

ŝ2 tot (x) = Var(Y tot (x)|D) = σ 2 tot k tot (x, x) -k tot (x, X)k tot (X, X) -1 k tot (X, x) . ( 29
)
Notice that if the combined GP hypothesis in ( 21) was used from the start to obtain the mean prediction, we would have recovered an additive Kriging model [START_REF] Durrande | Additive Covariance kernels for high-dimensional Gaussian Process modeling[END_REF]:

M (x) := E(Y tot (x)|D) = µ tot + k tot (x, X)k tot (X, X) -1 (Y -µ tot ). ( 30
)
However, the issue in ( 30) is that it involves k tot (X, X) -1 , the inverse of a sum of correlation matrices. As no direct formula exists for the inverse of a sum of matrices, the optimization of the weights for this alternative model (using a LOOCV procedure for instance) would involve a large number of matrix inversions and an inner optimization loop similarly to the classical hyperparameter optimization in ordinary Kriging. This is precisely what we aim at avoiding with the proposed method, in which both the weights of the combination in ( 20) and the weights for the variance in ( 25) are obtained analytically.

Amplitude of the variance

The final step to completely obtain the variance of the combination is to set the hyperparameter σ 2 tot in (21) which is used to calibrate the amplitude of the variance. In leave-oneout strategies, this hyperparameter is typically set by observing that the normalized LOO residuals should be normally distributed if the model is well-specified. Based on this, the usual method is to set the empirical variance of the normalized LOO residuals to 1 (see for instance Cressie, 1993, p. 102;[START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of gaussian processes with model misspecification[END_REF]:

σ2 LOO = 1 n n k=1 e (k) LOO s (k) LOO 2 , (31)
where, for k = 1, . . . , n, e

LOO are the LOO residuals of the model that can be obtained using the LOO formula (19) for the combination of Kriging models. s

(k)
LOO are the LOO standard deviations, obtained similarly:

e (k) LOO := y(x k )-M -k (x k ) = p i=1 w i [K -1 θ i (Y -µ i )] k [K -1 θ i ] k,k
, and s

(k) LOO := ŝ-k (x k ) = 1 [K -1 tot ] k,k
.

However, when the model is ill-defined (which is often the case in practice), the LOO residuals are not normally distributed and in particular, the presence of outliers can lead to an over-estimation the variance amplitude. Thus, we propose a similar approach more robust to outliers by fitting the amplitude hyperparameter using the empirical interquartile distance instead of the empirical variance:

σtot = IQ e LOO s LOO IQ norm , (32)
where IQ(.) designates the empirical interquartile distance, and IQ norm is the interquartile distance of a standard normal distribution.

Numerical results

In this section, we investigate the performances of the combination of Kriging models introduced in the previous section, and we compare it to ordinary Kriging. The focus is on high-dimensional problems where the number of samples used to build the model is limited, since we showed in section 3 that ordinary Kriging can fail to provide accurate models when the length-scales hyperparameters are estimated using MLE. As in section 3, we consider two analytical test functions in dimension 50 which are the sphere function (equation ( 8)) and trajectory samples from an isotropic GP (equation ( 9)). The sphere function is an example of a simple convex function which should be simple to model and to optimize. The GP trajectories are more complex multimodal functions, more difficult to optimize, and closer to practical usecases. In addition to these analytical functions, to validate the methods on a more realistic problem, they are also compared on a real-world application: the design of an electrical machine. The shape of the machine is parameterized by d = 37 design variables representing the position and size of air holes and magnets as well as the radius of the machine. The layout of the machine is illustrated in Figure 7. The performance of a machine is assessed by two objective functions which are its consumption and its cost. In addition, a valid machine must satisfy 10 constraints related to the dynamic of the car (e.g. maximum speed, acceleration, ...), to the dynamic of the machine (e.g. oscillation amplitudes, ...), and to the dimensioning of the reducer.

First, we will compare the global accuracy of the different models on these functions for a fixed number of observations. We will also assess the quality of the prediction intervals constructed. These two aspects are important in Bayesian optimization as an accurate surrogate model with correct error estimation should result in faster convergence, thus requiring less function evaluations to reach the optimum. Finally, since we are rather interested in faster convergence towards optimal designs than in obtaining a better global accuracy, we compare the performance of Bayesian optimization (EGO procedure) using the combination of Kriging models to that of Bayesian optimization equipped with standard ordinary Kriging.

Model precision and confidence intervals assessment

The global accuracy and confidence interval precision of the different methods are evaluated for the three test functions described above. The reference method is ordinary Kriging using a Matérn 5/2 covariance with hyperparameters estimated by MLE. The optimization is performed using the package DiceKriging in the R language [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF] with 300 maximum iterations of L-BFGS-B, and the length-scales lower and upper bounds are fixed to 0.1 and 20 respectively. The second method is the combination introduced in this paper using p = 16 sub-models (5 levels in the tree structure) with random lengthscales sampled following the procedure in section 4.1, with weights given by equation ( 20), and with variance obtained using equations ( 25), ( 22), ( 29) and ( 32). Finally, for the GP trajectories where the true length-scale θ true = 3 is known, we also compare both methods to a Kriging model with length-scales θ true as a reference. Each model is built using the same space-filling set of training points x 1 , . . . , x n train ∈ [0, 1] d obtained with Latin Hypercube Sampling (LHS). For the sphere function, the number of training points is n train = 250, and for the GP trajectories and the real-world application which are more complex we have n train = 500. The accuracy of each model is evaluated by computing the Q 2 (10) on a set of n test = 5000 test points

x (t) 1 , . . . , x (t)
ntest ∈ [0, 1] d sampled uniformly on the design space. To assess the precision of the confidence intervals, the coverage probabilities, corresponding to the proportion of test data that lies in the α% confidence interval, are computed on this same test set. For different levels of confidence α, we plot the coverage probabilities which should be equal to their theoretical level if the confidence intervals are well-calibrated. Each experiment is repeated for 10 different random seeds and theresults are presented in boxplots. For the electrical machine, as there are 12 functions to model (2 objectives and 10 constraints), the boxplot shows the performance averaged over the 12 functions.

The results for the three test functions are given in Figure 8. In each case, the combination of Kriging models is globally more accurate than ordinary Kriging with hyperparameters estimated by MLE as measured by the Q 2 indicator. In Figure 8d, the combination even manages to achieve a precision similar to that of the reference model (with the true lengthscales θ true ). The poor performances of ordinary Kriging can be explained by what was observed in Figure 2: MLE-optimal hyperparameters do not necessarily correspond to good Q 2 for high-dimensional problems with a limited number of observations. This is verified explicitly for the GP trajectories test case where the MLE-estimated hyperparameters differ from their true values resulting in a worse accuracy. In addition, unlike in Section 3 where the parameters of the hyperparameters optimizer were set to ensure convergence with 500 iterations and several restarts, here we take a lighter though more common setting in BO, with 300 iterations and no restart. This can lead to additional variability in the results for ordinary Kriging as the difficult length-scale optimization may not have fully converged in every cases. Still, for the GP trajectories, the log-likelihood at the end of the inner optimization is superior to the log-likelihood of the reference hyperparameters, even if the maximum likelihood optimization may not have fully converged. Overall, this shows that estimating the hyperparameters by MLE may not guarantee a good model.

Regarding the prediction intervals, we see in Figures 8b and8e that the ordinary Kriging is overconfident in its predictions as all probability levels CPs are lower than their theoretical values. This phenomenon was expected as the Kriging variance (3) usually underestimates the true prediction variance since it is obtained using an estimated covariance whose uncertainty on the estimation is not considered (Cressie, 1993, p. 127;[START_REF] Hertog | The correct kriging variance estimated by bootstrapping[END_REF]. In contrast, the prediction intervals for the combination in Figures 8c and8f are well calibrated as the CPs are aligned with the line y = x. For the electrical machine in Figure 8h and 8i, we observe a different behavior: the CPs are higher than their theoretical values at low probability levels and lower at high probability levels. One possible explanation is that the true prediction distribution is better approximated by a heavier-tailed distribution rather than by a Gaussian one. Still, we see that the confidence intervals for the combination are better calibrated than the original ones.

Bayesian optimization

The previous results show that the combination of Kriging models produces more accurate surrogates than ordinary Kriging with few points in high-dimension. In this section, we investigate whether more accurate models translate into faster convergence in BO. We perform EGO for the same three test functions: the sphere function given in equation 8, GP trajectories in equation 9, and the real-world application of an electrical machine. Note that in this paper, we only consider single-objective problems. As such, for the optimization of the electrical machine, we will only maximize on the first constraint which gives the maximum speed of the vehicle, as high values of this constraint are well correlated with good performances of the machine in general. However, as for ordinary Kriging, the combination can also be employed for constrained multi-objective Bayesian optimization using adapted acquisition criterion such as the Expected Hypervolume Improvement (EHVI) (see for instance [START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF]. For the sphere function and the GP trajectories, we consider 3 cases with varying dimension d = 15, 30, and 50. Each time, an initial design plan with n init = 2d is built by LHS, then n iter iterations of EGO are performed to add new samples. For the dimensions d = 15 and 30, we take n iter = 10d and for d = 50, we take n iter = 8d. At each iteration, the ordinary Kriging model and the combination are built using the same settings as described in the Section 5.1. The acquisition criterion is the Expected Improvement given in equation ( 7), and the EI maximization is performed with the R package DiceOptim [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]. To improve the performances, especially for high dimensions, we also use the TREGO [START_REF] Diouane | Trego: a trust-region framework for efficient global optimization[END_REF] algorithm to implement trust regions. Note that while the global minimum of the sphere function is known (min x f sphere (x) = 0), the global optimum for the GP trajectories and for the electrical machine are unknown, hence we only compare the performances of the two methods. The optimization is repeated for 10 different random seeds, with different initial design of experiments each time.

Figure 9 shows the evolution of the current best solution during the EGO iterations for the two methods. For the sphere function (Figure 9a, 9b and 9c) both methods are able to approach the global optimum of the function. However, for all dimensions, EGO with the combination converges faster than with ordinary Kriging, especially at the start of the optimization. The solution found after n iter iterations with ordinary Kriging is obtained within 2 to 3 times less iterations with the combination because the combination is globally more accurate than ordinary Kriging when very few samples are present, as is the case at the beginning of the optimization. This is confirmed by Figure 10 giving the global accuracy of both models during the EGO optimization for the sphere function with d = 30. After only a few iterations, the combination is able to achieve a good global accuracy in contrast Sphere function to ordinary Kriging where the Q 2 is still negative. Because it is more globally accurate at the start of the optimization, the combination of Kriging models is able to identify the interesting areas of the design space to explore, near the true global optimum of the function, more rapidly. Similarly, the convergence is improved with the combination on the GP trajectories with d = 15 (Figure 9d). For the cases with d = 30 and 50 in Figures 9e and9f, in addition to a faster convergence, we also observe that the final solution is significantly better than the one obtained with ordinary Kriging because, since the GP trajectories are multimodal functions, the less globally accurate ordinary Kriging models get trapped by a local optimum and miss the better one found by the combination. The same behavior is observed on the real-world application in Figure 9g (maximization problem in this case), where the EGO with the combination discovers a better solution faster.

Sphere function

Conclusion

In this paper, we have proposed a method to build a combination of Kriging sub-models for Bayesian optimization. We have presented a way to sample the sub-models suitably to avoid degenerate cases. A simple method based on LOOCV provides weights for the sub-models in closed-form, avoiding any costly inner optimization. Finally, we presented an approach to compute the prediction variance of the model, which is not available directly, by introducing a global covariance whose weights are computed analytically from the weights of the sub-models using an LOOCV strategy. An approach based on the empirical interquartile distance of the normalized LOO residuals, found to be more robust to outliers than the commonly employed empirical variance, was used to provide the amplitude hyperparameter of the variance.

The proposed method avoids the cumbersome optimization of the Kriging length-scales hyperparameters, which often fails to provide correct values in high-dimensional cases or when the number of training samples is limited as was demonstrated in a first numerical experiment. This is of interest for Bayesian optimization as such low samples regimes are common in design optimization, especially at the start of the optimization. The numerical results we obtained show that the combination of Kriging models enables both the construction of more accurate models and the faster convergence towards optimal solutions in Bayesian optimization, especially in the start when the number of samples is very limited. 2 to 3 times less expensive function evaluations are required with our method and in multimodal cases, better solutions never found by the original method are obtained. Contrarily to other high-dimensional BO methods, such as dimension reduction or additive Kriging, no additional hypothesis on the underlying function, such as a low dimension representation or an additive structure, is required. Thus, the proposed method is easily generalizable to any design engineering problem.

Several aspects still need to be explored in further research. First, in this paper the proposed method was compared with a vanilla EGO using an ordinary Kriging model. It would be interesting to compare it with other existing high-dimensional BO approaches (e.g. [START_REF] Eriksson | Scalable global optimization via local bayesian optimization[END_REF][START_REF] Diouane | Trego: a trust-region framework for efficient global optimization[END_REF][START_REF] Bouhlel | Efficient global optimization for high-dimensional constrained problems by using the kriging models combined with the partial least squares method[END_REF][START_REF] Antonov | High dimensional bayesian optimization with kernel principal component analysis[END_REF][START_REF] Binois | On the choice of the low-dimensional domain for global optimization via random embeddings[END_REF] on a more complete benchmark. Furthermore, we have only considered single objective optimization here. Extending the method to multi-objective and/or constrained problems is a perspective that would allow it to tackle a wider range of realistic problems such as the electrical machine design optimization.

A More details on the variance of the combination

A.1 Relation between mixture and combination

In this appendix, we will detail the proof of the statement in equation 23 that the MSE of a mixture model is equal to the MSE of the linear combination of the same models when the weights α 2 of the global covariance in equation ( 21) are equal to the weights w of the combination. First, for the combination of Kriging models, we have:

M tot (x) := p i=1
w i M i (x), with M i (x) = k θ i (x, X)k θ i (X, X) -1 Y, i = 1, . . . , p.

Here, without loss of generality, we suppose the models centered. Using the definition for the global GP in equation ( 21), we obtain the following global covariance:

Cov(Y (x), Y (x ′ )) = σ 2 tot p i=1 α 2 i k θ i (x, x ′ ). ( 33
)
Then the MSE of the combination is:

E (M tot (x) -Y (x)) 2 = Var (Y (x)) + Var (M tot (x)) -2Cov (M tot (x), Y (x)) = σ 2 tot p i=1 α 2 i k θ i (x, x) + w ⊤ K M (x)w -2w ⊤ k M (x) , (34) with:
(K M (x)) i,j = k θ i (x, X)k θ i (X, X) -1 p ℓ=1 α 2 ℓ k θ ℓ (X, X)k θ j (X, X) -1 k θ j (X, x), i, j = 1, . . . , p, (k M (x)) i = k θ i (x, X)k θ i (X, X) -1 p ℓ=1 α 2 ℓ k θ ℓ (X, x), , i = 1, . . . , p.

For the mixture now, we have a mixture of posterior GPs: Re-injecting these two expressions into (35):

E (M mix (x) -Y I (x)) 2 = σ 2 tot p i=1 w i k θ i (x, x) + w ⊤ K M (x)w -2w ⊤ k M (x) , (36) with:
(K M (x)) i,j = k θ i (x, X)k θ i (X, X) -1 p ℓ=1 w ℓ k θ ℓ (X, X)k θ j (X, X) -1 k θ j (X, x), i, j = 1, . . . , p,

(k M (x)) i = k θ i (x, X)k θ i (X, X) -1 p ℓ=1
w ℓ k θ ℓ (X, x), , i = 1, . . . , p.

We have equality of the MSE of the combination in (34) and of the mixture model in (36) if α 2 = w, hence the result announced.

A.2 Coefficients for the variance of the combination

In this appendix, we give the proof for the analytical expression of the optimal weights α of the global covariance given in equations ( 25), ( 26) and ( 27).

As defined in equation ( 24), the optimal weights are obtained by minimizing the expected MSE of the combined model with respect to Y tot . For a combination of two models, with Y tot = αY 1 + (1 -α)Y 2 , we recall the expression:

α * = arg min α E E (wM 1 (x) + (1 -w)M 2 (x) -αY 1 (x) -(1 -α)Y 2 (x)) 2 |Y 1 , Y 2 .
Similarly to the derivation of the combination weights w in equations ( 17) and ( 18), the global MSE is approximated using the LOOCV error:

α * = arg min α E [e LOOCV (wM 1 + (1 -w)M 2 )|Y 1 , Y 2 ] , (37)
where the LOOCV error is obtained using the Kriging LOO formula ( 19):

e LOOCV (wM 1 + (1 -w)M 2 ) = n k=1 wM 1 -k (x) + (1 -w)M 2 -k (x) -αY 1 (x k ) -(1 -α)Y 2 (x k ) 2 = n k=1 (W k Y tot (X)) 2 ,
where:

W k := w [K -1 θ 1 ] k, [K -1 θ 1 ] k,k + (1 -w) [K -1 θ 2 ] k, [K -1 θ 2 ] k,k . (38)

Figure 1 :

 1 Figure 1: Example of an ordinary Kriging model in 1D. Left: model with an arbitrary length-scale θ = 0.03, right: model with length-scale obtained by MLE.

Figure 2 :

 2 Figure 2: Q 2 of Kriging models (the higher the better) and corresponding log-likelihood for different number of training points over 10 independent runs. The top two figures give the results for the Sphere test function and the bottom two figures the results for GP trajectories. On the left are the accuracies of the models measured by the Q 2 , on the right are the corresponding log-likelihood. The results in red correspond to the Kriging model with hyperparameters estimated by maximum likelihood, those in blue to the model with reference hyperparameters.

Figure 3 :

 3 Figure 3: Examples of the two degenerated cases in 1D. Left: ordinary Kriging with a small lengthscale θ = 10 -2 , right ordinary Kriging with a large length-scale θ = 3 and a small nugget term ε 2 = 10 -8 .

Figure 4 :

 4 Figure 4: Entropy of a Gaussian correlation as a function of the length-scale θ. In this example, d = 50 and X are uniformly distributed (σ 2 X = 1 12 and κ X = 9 5 ).

Figure 5 :

 5 Figure 5: Illustration of the sub-models combination following a binary tree structure (3 levels for 4 sub-models here).

Figure 6 :

 6 Figure 6: Illustration of the two-by-two GPs combination following a binary tree structure (3 levels for 4 individuals GPs here) with weights α (right) associated to the two-by-two combination of sub-models with weights w (left).

Figure 7 :

 7 Figure 7: Layout of the electrical machine. The 37 design parameters are the size and position of the air holes (in white) and of the magnets (in black), as well as the radius of the machine.

Figure 8 :

 8 Figure 8: Left column: accuracy of the models measured by the Q 2 . Central and right column: precision of the confidence intervals for the ordinary Kriging (center) and the combination (right), the closest to the y = x line the better. First row: sphere function, middle row: GP trajectories, bottom row: electrical machine. Red boxes stand for Kriging models with hyperparameters obtained by MLE and blue boxes for the combination of Kriging models.

  (a) d = 15 (b) d = 30 (c) d = 50 GP trajectories (d) d = 15 (e) d = 30 (f) d = 50 Electrical machine (g) d = 37

Figure 9 :

 9 Figure 9: Evolution of the current best solution during the EGO iterations. The x-axis corresponds to the number of function evaluations. Top row: sphere function, middle row: GP trajectories, bottom row: electrical machine. Red: EGO using an ordinary Kriging model, blue: EGO using the combination of Kriging models.

Figure 10 :

 10 Figure 10: Q 2 of the ordinary Kriging model (red) and the combination of Kriging models (blue) during the optimization (sphere function with d = 30, Figure 9b).

Y

  I (x)|Y I (X), with density p i=1w i p(Y i (x)|Y i (X)),where I is the indicator such that P (I = i) = w i , i = 1, . . . , p. The corresponding mixture model is equal to the combination:M mix (x) = E (Y I (x)|Y I (X)) = p i=1 w i M i (x).The MSE of the mixture model is obtained by:E (M mix (x) -Y I (x)) 2 = Var(Y I (x)) + Var(M mix (x)) -2Cov(M mix (x), Y I (x)). (35)Using the lax of total variance:Var(Y I (x)) = E(Var(Y I (x)|I)) + Var(E(Y I (x)|I)) = σ 2 tot p i=1 w i k θ i (x, x) + 0,and: Cov(M mix (x), Y I (x)) = E(Cov(M mix (x), Y I (x)|I)) + Cov(E(M mix (x)|I), E(Y I (x)|I)) j Cov(M i (x), Y j (x)).

Acknowledgments

This research was conducted with the support of the consortium in Applied Mathematics CIROQUO (https://doi.org/10.5281/zenodo.6581217), gathering partners in technological and academia in the development of advanced methods for Computer Experiments. This research was partly funded by a CIFRE grant (convention #2021/1284) established between the ANRT and Stellantis for the doctoral work of Tanguy Appriou.

Let K tot := α 2 K θ 1 + (1 -α) 2 K θ 2 be the global correlation matrix. Then, Y tot (X) ∼ N (0, K tot ), and W k Y tot (X) ∼ N (0, W k K tot W ⊤ k ). Thus, the expected LOOCV error is:

Taking the derivative with respect to α:

Plugging back the expression for W k (38) in (39):

, and:

.

We can also rewrite these results using the expected LOOCV error seeing that:

Finally, by setting the partial derivative equal to 0 in (39), and using the expression in (40), we obtain the result given in equations ( 24) for the optimal weights α * solution of (37).


Apurva V Oza
email: apurva.v.oza@jpl.caltech.edu

Francois Leblanc

Robert E Johnson

DUSK/DAWN ATMOSPHERIC ASYMMETRIES ON TIDALLY-LOCKED SATELLITES II: THERMAL TIDES AND OUTGASSING AT THE GALILEAN SATELLITES

Satellite Tides Background

Thermal tides on a tidally-locked satellite have been studied as gas flows on a rotating surface since the '60s [START_REF] Hodges | Lateral transport in planetary exospheres[END_REF]. The authors derived two relations nT 5/2 and nT = constant, by allowing non-uniform gas concentrations (number density (n)gradients)) to drive lateral flows at a given temperature, T. In effect, when tenuous gas in an exosphere is coupled to a rotating surface, the semidiurnal thermal tide (∝ rotation rate Ω, and temperature gradient dT/dϕ) pulls the gas away from the expected thermal flux maxima F at noon (see Figure 1 sketch, motivated by tidal density perturbations at close-in exoplanets [START_REF] Arras | Thermal Tides in Fluid Extrasolar Planets[END_REF]). The expressions studied for lateral transport in planetary exospheres are unable to reproduce the location ϕ and timing t, of the peak O 2 gas column densities N, simulated on the tidally-locked satellites Europa (JII) [START_REF] Apurva | Dusk over dawn O 2 asymmetry in Europa's nearsurface atmosphere[END_REF] and Ganymede (JIII) [START_REF] Leblanc | On the orbital variability of Ganymede's atmosphere[END_REF], observed so far only in oxygen emission at JII [START_REF] Roth | Europa's far ultraviolet oxygen aurora from a compre-hensive set of HST observations[END_REF].

Atmospheric Evolution and Escape Model Here, we build on the rotating 1-D mass conservation model (nommé dishoom (desorbing interiors via satellite heating to observe outgassing model)) in Paper I (Oza, Johnson, and Leblanc [START_REF] Apurva | Dusk/dawn atmospheric asymmetries on tidally-locked satellites: O 2 at Europa[END_REF]) where we showed that the density peaks consistently at dusk only if a thermal source is used to source the oxygen aurorae observed by the Hubble Space Telescope (HST). Since orbital longitude ϕ probes the axis of time (and the associated surface heating dT/dϕ) our analytic model is fundamentally tidal in nature. Thermal tides can therefore be useful in describing the exospheric accumulation of O 2 /H 2 O and volatiles generally given our recent understanding on the thermal nature of the O 2 population (binding energy, U s = 0.14 eV (Johnson et al. [START_REF] Johnson | The Origin and Fate of O 2 in Europa's Ice: An Atmospheric Perspective[END_REF]) ). Furthermore, recent HST evidence suggests that at the sunlit trailing hemispheres of JII and JIII (Roth et al. [START_REF] Roth | A sublimated water atmosphere on Ganymede detected from Hubble Space Telescope observations[END_REF]) H 2 O may generate a locally-collisional atmosphere, whose density remains to be accurately constrained by future ground, space, or insitu spacecraft. Although, it is generally agreed these O Implications Studying atmospheric evolution on surface-bounded atmospheres is valuable in that properties of the icy surface Johnson et al. [START_REF] Johnson | The Origin and Fate of O 2 in Europa's Ice: An Atmospheric Perspective[END_REF], and its interior Birds-eye view illustrating the two tidal components acting on a tidally-locked satellite atmosphere/exosphere. Adapted from Arras and Socrates [START_REF] Arras | Thermal Tides in Fluid Extrasolar Planets[END_REF].

Hesse et al. [START_REF] Hesse | Transport of surface oxidants into internal oceans by brine migration through ice shells[END_REF] can be revealed. These constraints may also be able to inform the formation of primordial icy bodies (e.g comets) in the protosolar nebula (Oza and Johnson [START_REF] Oza | A common origin of oxygen at Jupiter's icy moons and comets: Modeling thermal outgassing at Europa[END_REF]).

Thermal Outgassing of O 2 /H 2 O at Europa, Ganymede, and Callisto Unlike JI, where the freezing point of SO 2 frost (201 K) poses no challenge to our understanding of the Ionian surface-atmosphere boundary layer, the trapped O 2 observed at JII, JIII, and JIV continues to be puzzling (Spencer, Calvin, and Person [START_REF] Spencer | CCD Spectra of the Galilean Satellites: Molecular Oxygen on Ganymede[END_REF]; Spencer & Calvin 2002) as the trapped O 2 in amorphous or crystalline ice grains must thermally outgas since P vap,O2 ≫ P JII-JIV . Figure 2 provides a model considering the diurnal tide acting on the surface ice, from a range of regolith temperatures representative of the Galilean satellite surfaces. If the diurnal tide is able to sufficiently anneal and release trapped volatiles from inclusions/bubbles Johnson and Jesser [START_REF] Johnson | O 2 /O 3 Microatmospheres in the Surface of Ganymede[END_REF], this model is a reasonable feedback mechanism for the icy Galilean satellite atmospheres, providing a direct parallel to volcanic SO 2 frost on JI. This continues the idea that O 2 is indeed accessible to the atmosphere as a surface frost in quasi-vapor pressure equilibrium (Paper I).

Summary If a resonance exists between the atmospheric lifetime and rotation rate a dusk-over-dawn at- mospheric asymmetry appears on tidally-locked satellites as shown for ultraviolet HST observations in Paper I [START_REF] Apurva | Dusk/dawn atmospheric asymmetries on tidally-locked satellites: O 2 at Europa[END_REF]. Evidence of thermal outgassing of trapped volatiles may be present in spectra of the newly launched JWST, equipped with the mid-infrared detector MIRI [START_REF] Lagage | Characterisation of exoplanet atmosphere with the JWST and then ARIEL[END_REF]. Future observations may reveal thermal tide signatures, in the form of phase-curve variability as also studied for close-in, asynchronous exoplanets [START_REF] Arras | Thermal Tides in Fluid Extrasolar Planets[END_REF])

  2 /H 2 O atmospheres are more tenuous than SO 2 (≪ 10 17 cm -2 ) at Io (JI) the simulated near-surface atmospheres at JII & JIII indicate that JI, JII and JIII are all indeed asymmetric towards dusk ((Oza et al. [3]; Leblanc et al. [4]; Walker et al. [9]; Lellouch et al. [10]).

Figure 1 :

 1 Figure 1:Birds-eye view illustrating the two tidal components acting on a tidally-locked satellite atmosphere/exosphere. Adapted from Arras and Socrates[START_REF] Arras | Thermal Tides in Fluid Extrasolar Planets[END_REF].

Figure 2 :

 2 Figure 2: Surface-Atmosphere Exchange system modeled by dishoom where the trapped column is estimated following Johnson & Jesser 1997 and Johnson et al. 2019. The sublimation and ice crystallization temperatures are illustrated for O 2 sublimation (54.36 K) and H 2 O ice crystallization (≈ 150 K) [15]. Thermal outgassing/sublimation of water is observed in the lab to be significant ≫ 150 K [16]. The outgassed column densities are normalized to rough constraints by HST oxygen aurorae observations (Hall et al. 1998).

Acknowledgments:

Part of this work was conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract