Amandine Marrel

Bertrand Iooss

Probabilistic surrogate modeling by Gaussian process: A new estimation algorithm for more robust prediction

Keywords: Computer experiments, Gaussian process regression, Machine learning, Optimization, Uncertainty, Validation criteria

In reliability engineering studies, computer codes are increasingly used to model physical phenomena which, in many cases, can be very time-consuming to run. A widely accepted approach consists in approximating the CPU-time expensive computer model by a surrogate model. One of the most popular surrogate model is the Gaussian Process regression, as it provides, additionally to a prediction at an unobserved point, an uncertainty around this prediction (a predictive distribution). However, in practice, the quality of this metamodel depends on several choices, as the estimation and validation algorithms. The present work aims at proposing a new algorithm, based on constrained optimization multi-objective techniques, to estimate the Gaussian process hyperparameters in order to ensure robust and accurate (i.e. reliable) predictive distribution of the Gaussian process. An intensive numerical benchmark on various analytical functions, with different input dimensions and learning sample sizes, shows its good performance in comparison with standard estimation algorithms. The new algorithm is also applied to a real test case modeling an aquatic ecosystem. It is compared with a recent robust and sophisticated Bayesian method; it proves to be as efficient while being less sensitive to the specification of the Gaussian process model.

Introduction

In several engineering fields, numerical models are used to simulate physical phenomena of industrial or environmental processes in order to improve its understandings, optimize some performances or perform a risk analysis related to safety criteria. In this context, one key issue is that the numerical model under study can be very time-consuming to run, which can drastically limit the number of possible simulations. To solve this cost issue, a widely accepted approach consists in approximating the CPU-time expensive computer model by a CPU-time inexpensive mathematical function called "surrogate model" (or "metamodel", term that is used in the following). Fit from a set of inputs and outputs of computer code simulations, these metamodels can come from any machine learning technique, from the simplest as the polynomial regression, to more complex as the random forests and neural networks, see the reviews of Villa-Vialaneix et al. [START_REF] Villa-Vialaneix | A comparison of eight metamodeling techniques for the simulation of N2O fluxes and N leaching from corn crops[END_REF] in the environmental domain and Afshari et al. [START_REF] Afshari | Machine learning-based methods in structural reliability: A review[END_REF] in structural reliability.

In a first paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF], companion of the present one, the value of Gaussian process regression (GPR) metamodel for emulating costly computational codes has been emphasized. Thanks to its great flexibility, this non-parametric regression tool has proved highly effective in modelling numerical simulators, in a wide range of application fields [START_REF]Springer Handbook on Uncertainty Quantification[END_REF]. Moreover, it yields an analytical predictive distribution for the code output at each prediction point. Having a such probabilistic metamodel, in the sense that it provides a predictive law for each new evaluation point, is of great added value, particularly for safety, reliability or risk assessment studies. It also enables the deployment of sophisticated Gaussian process-based approaches for active learning, robust optimization, reliability assessment, etc, as reviewed for example in Fuhg et al. [START_REF] Fuhg | State-of-the-art and comparative review of adaptive sampling methods for kriging[END_REF], Moustapha et al. [START_REF] Moustapha | Active learning for structural reliability: Survey, general framework and benchmark[END_REF]. In this context, it is essential to guarantee confidence in the Gaussian process (GP) metamodel predictive law, and not just in its mean value (i.e. the GP predictor). In many engineering problems, having access to prediction intervals in addition to predictions is essential. This major advantage of GP has been shown to be particularly useful for various objectives such as prediction-forecast [START_REF] Betancourt | Gaussian process metamodeling of functional-input code for coastal flood hazard assessment[END_REF], sensitivity analysis [START_REF] Marrel | Calculations of the Sobol indices for the Gaussian process metamodel[END_REF], calibration [START_REF] Perrin | Adaptive calibration of a computer code with time-series output[END_REF], probabilities of failure [START_REF] Huang | A new active learning kriging metamodel for structural system reliability analysis with multiple failure modes[END_REF], reliabilitybased design optimization [START_REF] Ma | Reliability-based design optimization using adaptive kriging-a single-loop strategy and a double-loop one[END_REF], robust optimization [START_REF] Ribaud | Robust optimization: A kriging-based multi-objective optimization approach[END_REF] and functional risk curves assessment [START_REF] Iooss | Uncertainty and sensitivity analysis of functional risk curves based on Gaussian processes[END_REF]. The main idea is that, in addition to the mean quantity of interest, the stochastic information brought by the GP metamodel allows to obtain confidence bands on the quantity of interest.

This confidence first calls for a robust estimation of the GP metamodel and, more precisely, of its parameters, referred to as hyperparameters. Secondly, a rigorous validation of the entire GP predictive distribution is required, as extensively outlined by Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF] and Petit et al. [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF]. With regard to the first requirement, the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF] has reviewed recent works dealing with the estimation of GP hyperparameters, from a theoretical and empirical point of view. It appears that the usual methods, based mainly on likelihood maximization or on minimizing an error estimated by cross-validation, sometimes lead to poor-quality and provide non-robust estimates. In such cases, the validation step (the second requirement) must detect any inaccuracy in the GP predictive distribution. For this, it must be based on different criteria [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF][START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF] evaluating the GP's predictive capabilities, as well as the correctness of the associated prediction intervals. Furthermore, particular care and informed consideration must also be taken when jointly analyzing these criteria.

To ensure a more robust estimation of hyperparameters and then a more robust and accurate (denoted "reliable" in the following) prediction distribution, some recent alternatives to standard estimation approaches (which are based on maximizing the likelihood or minimizing the crossvalidation prediction error) have been proposed. Among them, Bayesian approaches are theoretically very attractive, as they offer a kind of likelihood regularization. However, their deployment in practice is limited by their complexity, not only in terms of computational cost especially in large dimension (large number of inputs), but also in terms of the expertise required to define so-called robust priors and analyze the resulting posteriors. Other alternative approaches based on ad-hoc corrections of the quantiles of the GP predictive distribution have also been proposed, as in Acharki et al. [START_REF] Acharki | Robust prediction interval estimation for Gaussian processes by cross-validation method[END_REF], in order to ensure reliable prediction intervals for a given level. However, these approaches do not necessarily seem relevant for obtaining a metamodel for multi-objective use.

The presented work aims at proposing a new algorithm to estimate the GP hyperparameters in order to ensure robust and reliable GP predictive distribution. This algorithm is based on a preliminary thorough analysis of estimation and validation criteria, both through the recent literature on the subject (reviewed in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]) and through an empirical exploration of their links on a large benchmark of analytic functions. On this basis, the new algorithm proposes to jointly maximize the likelihood of observations and the accuracy of GP prediction intervals, under the constraint of not degrading the GP predictivity.

The rest of the document is organized as follows. First, a brief reminder of the formalism of GP regression, estimation methods and validation criteria are first proposed. The interested reader is invited to refer to the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF] for more details on these subjects. Then, an empirical study of the connections and links between estimation and validation criteria is presented in Section 3. From this, the new estimation algorithm is proposed and detailed in Section 4. An intensive numerical benchmark is performed in Section 5 on various analytical functions from the literature on emulation and prediction of computer experiments. This benchmark makes it possible to evaluate, for different input dimensions and learning sample sizes, the additional value brought by the new algorithm (compared to standard algorithms). Finally, in Section 6, the algorithm is applied to a real test case modeling an aquatic ecosystem and more precisely a prey-predator chain. The last section gives some conclusions and prospects of this work.

Reminders on Gaussian process regression, estimation and validation

All the notations introduced in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF] are retained, and those necessary for the present paper are recalled in what follows. The numerical model (computer code or simulator) is represented by the following input-output relationship:

M : X -→ Y X -→ Y = M(X) (1)
where X = (X 1 , . . . , X d ) ⊤ are the d uncertain input parameters belonging to some measurable space X ⊂ R d . In the context of given-data and data-driven GP metamodeling, we only have a n-size learning sample of inputs and associated outputs denoted by (X s , Y(X s )) where X s = {x (1) , . . . , x (n) }, with

x (i) = x (i) 1 , . . . , x (i)
d , denotes the matrix of n-size sample locations, and Y s = {y 1 , . . . , y n } denotes the corresponding outputs observations with y i = M(x (i) ).

Gaussian process metamodel and associated parametric choices

In the GPR, Y is considered as a realization of a Gaussian stochastic process:

Y (•) ∼ N (m(•), k(•, •)) , ( 2
)
where m(•) is the mean function and k(•, •) is the covariance kernel function. The predictive GP distribution is then the GP conditioned by the learning sample (X s , Y s )

Y (•) | Y s ∼ N (ŷ(•), ĉ(•, •)) , ( 3
) with ŷ (x) = E[Y (x)|Y s ] = m(x) + k(x, X s ) T K -1 Y s -m(X s ) , (4) and ĉ
(x, x) = COV[Y (x), Y (x)|Y s ] = k(x, x) -k(x, X s ) T K -1 k(x, X s ), (5)
where

m(X s ) = (m(x (i) )) 1≤i≤n ∈ R n , K = k(x (i) , x (j) )) 1≤i,j≤n ∈ R n×n and k(x, X s ) = (k(x, x (i) )) 1≤i≤n ∈ R n .
The efficiency of GP metamodeling strongly depends on the specifications of its regression and covariance functions. In practice, this consists of making parametric choices for m(•) and k(x, x) among a panel of usual functions. First, in the absence of prior knowledge, a constant m(x) = β 0 or a one-degree polynomial trend m(x) = β 0 + i β i x i is usually considered.

Secondly, concerning k(x, x), as the isotropy hypothesis is often too restrictive to emulate complicated models with inputs of very different kinds and influence, an anisotropic covariance is generally assumed. More precisely, the following tensorized covariance is considered:

k σ,θ (x, x) = σ 2 d i=1 k θ i (x i -xi ). ( 6
)
where σ 2 is the variance parameter and θ i ∈ R + correlation hyperparameter (also called correlation length or length-scale) associated to the i th input. As discussed in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF], the d 1-D covariance functions can be of different natures. But in practice, given the large number of inputs and without any prior knowledge, the usual practice is to use the same function for all variables.

A very popular choice in machine learning, and in particular for the emulation of computer code, is the Matérn function defined in one dimension by:

k ν,θ (x, x) = 2 1-ν Γ(ν) √ 2ν|x -x| θ ν K ν √ 2ν|x -x| θ , ( 7
)
where K ν is a modified Bessel function of second kind with parameter ν ∈ R + , and Γ is the Euler Gamma function. Parameter ν controls the smoothness of the GP (cf. Table 1 of the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]). Usual choices in machine learning are ν = 1/2 (exponential covariance), ν = 3/2 and ν = 5/2 (referred to as 3/2-Matérn and 5/2-Matérn covariances), and the Gaussian covariance which can be viewed as the the limiting case of Matérn function when ν → ∞. However, for the emulation of calculation codes, a Gaussian covariance (with infinitely differentiable GP trajectories) is not recommended, as few physical phenomena are infinitely smooth. At the opposite extreme, an exponential covariance with non differentiable trajectories will be probably too rough, since the outputs of numerical simulators are often the results of solving partial differential equations.

A reasonable choice is rather to consider intermediate choices such as 3/2-Matérn or 5/2-Matérn covariances, as highlighted by Gramacy [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF] or Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] (See also [3, Section 2.2.1] for further details). Finally, an additional nugget effect can be considered. It consists in assuming an additive white noise effect of variance σ 2 ϵ which results in an additional term σ 2 ϵ δ xx in the covariance function. In practice, it is characterized by the parameter λ = σϵ σ 2 ∈ R + . Remark 1. The kriging and GPR frameworks are often considered equivalent, but there are a few differences beyond the vocabulary, which lie mainly in inference for unknown quantities (to be predicted) and in practical estimation. First, they are built-upon different points of view: kriging focuses on estimating the best linear unbiased estimator, while GPR adopts a Bayesian viewpoint where a prior GP is updated by conditioning with the observed data. The formulas obtained in both cases for the predictive mean and variance match (due to Gaussian conditioning properties): for a GPR with prior constant mean function (resp. polynomial function), we obtain ordinary (resp. universal) kriging formulas. The main difference thus lies in the additional assumption of Gaussian distribution for the output in GPR, which allows to build confidence intervals. Secondly, the approaches differ in their practical implementation: kriging, originally developed for the spatial framework, relies on the analysis and adjustment of the empirical variogram to select the covariance kernel and estimate its hyperparameters, whereas GPR most often uses the likelihood of the data, as detailed in the next section. As a result, kriging focuses on problems with low input dimensions (d ≤ 3), whereas GPR can handle much higher-dimensional problems and offer much greater flexibility in the choice of covariance models (e.g. anisotropic choice even for d > 3). For more details on differences between kriging and GPR frameworks, we refer the interested reader to the recent paper of Christianson et al. [START_REF] Christianson | Traditional kriging versus modern gaussian processes for large-scale mining data[END_REF].

Main estimation methods

Assuming the aforementioned parametric choices, estimating the GP from the learning sample therefore boils down to estimating its parameters (β, σ 2 , θ) and eventually λ if a nugget effect is considered. As described in the companion paper [3, Section 3], three main estimation procedures are used: minimization of the squared prediction error calculated by cross-validation (CV), maximization of likelihood (denoted MLE for maximum likelihood estimation), and Bayesian approach.

For the first method, and considering the leave-one-out (LOO) procedure (specific case of CV), the GP parameters are computed by minimizing the LOO-MSE:

LOO-MSE β, σ 2 , θ := 1 n n i=1 (y i -ŷ-i ) 2 ,
where ŷ-i denotes the mean of GP predictive distribution in x (i) when (x (i) , y i ) is removed from the learning sample (this comes down to consider the GP conditioned by Y s,-i ). Note that the calculation of the LOO predictive distribution is greatly facilitated by the use of CV Dubrule [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF]'s formulas, which provide analytical formulation for LOO mean ŷ-i and LOO variance ŝ2

-i = ĉ-i (x (i) , x (i) ).
The MLE method consists in identifying the values of the parameters that minimize the negative log-likelihood (NLL) of the learning sample. Provided that θ is know, analytical solutions are available for (β, σ 2 ) and MLE boils down to minimize the profile NLL involving just θ. But, as for LOO-MSE, there is no closed-form expressions for the optimal values of parameters θ and numerical methods are thus required to estimate them.

The last approach, the Bayesian one is close to the MLE as it considers the likelihood of the learning sample but it also incorporates a prior distribution for the GP parameters. Bayesian estimation thus relies on maximizing the marginal posterior distribution of the parameters with regard to this prior distribution.

As discussed in the companion paper [3, Section 3.5], even if no consensus really emerges from theoretical analysis and empirical comparison of estimation methods, the recent empirical results of Petit [START_REF] Petit | Improved Gaussian process modeling: Application to Bayesian optimization[END_REF] tend to argue that MLE is often preferable to its competitors. But MLE (as LOO-MSE) can be theoretically an ill-posed problem [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF][START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] and flatness of MLE (or LOO-based criteria) around the optimal hyperparameter value can lead to poor performance of optimization algorithms. For more details on the MLE ill-posedness problem and its lack of robustness, we refer the interested reader to the dedicated section in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]Section 3.3], where recent papers of Karvonen and Oates [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF] and Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] are reviewed. An illustration of the flatness of the objective function derived from the MLE is also provided on a very simple analytical example (see [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]Figure 4]). The full-Bayesian approach under the assumption of suitable priors, offer better properties but suffer from poor tractability. However, the recent RobustGaSP Bayesian method proposed by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] (see also [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]Section 5.2]) with its specific robust priors and approximations, could overcome these limitations. Finally, considering a nugget effect jointly estimated with the other hyperparameters can facilitates the estimation by MLE by regularizing the likelihood function, improving the conditioning of the correlation matrix and numerical convergence of algorithms. It is of much less interest in the Bayesian approach, where the prior already plays a regularizing role and the additional estimation of the nugget can increase the identifiability problems.

Beyond the problem of obtaining a robust estimate of the hyperparameters θ, it also appears that the choice of regularity of covariance function (ν in Matérn class) might be a key element to ensure GP's predictive capabilities.

Validation criteria

Whatever the parametric choices and the estimation method chosen, it is essential to have criteria for assessing the reliability of the obtained GP predictive distribution. These validation criteria can be used both to select the most appropriate covariance function and to check that the estimated hyperparameters lead to a reliable predictive distribution. They can be estimated from a test sample (different from the learning sample), or, as is often the case in our small-data application context (i.e. limited number of simulations), by cross-validation (or LOO, as a special case of cross validation). Hence, as described in Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF], different validation criteria have been proposed to assess the whole GP conditional distribution:

• predictivity coefficient Q 2 to assess the accuracy of the GP predictor ŷ (x):

Q 2 = 1 - 1 n n i=1 (y i -ŷ-i ) 2 1 n n i=1 y i -1 n n i=1 y i 2 = 1 - LOO-MSE 1 n n i=1 y i -1 n n i=1 y i 2 .
The closer to one the Q 2 , the better the accuracy of the metamodel predictor.

• predictive variance adequacy (PVA) factor introduced by Bachoc [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF], in order to check if the conditional GP variance is of the right order of magnitude: PVA = log

1 n n i=1 (y i -ŷ-i ) 2 /ŝ 2 -i .
The smaller the PVA, the more reliable the prediction intervals.

• IAEα criterion to assess the reliability of GP prediction intervals: IAEα = 1 0 ∆(α) -α dα where ∆(α) is the empirical coverage function [START_REF] Schruben | A coverage function for interval estimators of simulation response[END_REF]. IAEα quantitatively summarizes the α-PI plot of Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF]. It lies in [0, 1] and the closer to zero the IAEα, better the PI in average.

For a synthetic overview of the various validation criteria, we refer the interested reader to the companion paper [3, Appendix A]. To tackle the problems of estimating hyperparameters, and given the importance of having a reliable GP predictive distribution in the context of uncertainty quantification, we are convinced of the benefits of considering these criteria, rather reserved for validation, directly in the estimation process, in addition to NLL. And to find out how best to use and combine them, the links between them are firstly investigated, according to the values of the hyperparameters.

Empirical study of connections between likelihood and validation criteria

To design the new estimation algorithm, a preliminary study considers the links between some of the estimation and validation criteria, according to the values of the hyperparameters. Based on the review of the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF] about recent works on GP estimation, the MLE is considered as the main estimation criterion. The interest is then to assess how this criterion can really control the quality of the predictive distribution (under the assumption that the optimization of this criterion has converged well). For this, the three following validation criteria (to be minimized) are considered: the MSE (∝ 1 -Q 2 ), the PVA and the IAEα.

A wide range of analytical functions, commonly used in metamodeling benchmarks, are then considered. For each function, random learning samples of different size n are simulated to build the conditional GP (Eqs 4,5). A GP with constant mean is considered in all that follows: β := β 0 ∈ R. The domains of variation for the hyperparameters θ are discretized finely (on a grid for example if d = 2). For each set of θ values (i.e. point of the grid if d = 2), the corresponding MLE estimates of (β, σ) are computed (from analytical MLE formulas of (β, σ) as functions of θ), and the associated GP conditioned by the learning sample is built. Therefore, for each set of θ values, Different criteria are then computed:

• on the one hand, the negative log-likelihood (NLL) of the learning sample (criterion that one would seek to minimize in practice in the case of MLE),

• on the other hand, three validation criteria of the conditional GP computed on a large test sample of 1000 points: (1 -Q 2 ), PVA and IAEα. Ideally, these three criteria should be as low as possible to ensure GP accuracy and predictivity (as detailed in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]).

In summary, for all the aforementioned criteria, the objective is to have the minimum value. Different sample sizes n are considered, well chosen w.r.t. the analytical functions tested in order to encompass cases of underlearning (n too small w.r.t. the complexity of the function) and overlearning (n much larger than necessary). For each test case, all the procedure is repeated 100 times for different i.i.d. random learning samples.

Illustration on an infinitely differentiable toy function

This section provides some results representative of the general trends observed over the 100 repetitions. First, the "re-scaled" Branin function [START_REF] Picheny | A benchmark of kriging-based infill criteria for noisy optimization[END_REF], illustrated in the companion paper [3, Section 4.2 and Figure 2], is considered. The results obtained for the four quantities of interest aforementioned are given for the exponential, 3/2-Matérn, 5/2-Matérn and Gaussian covariances, by Figure 1 for n = 30. Similar plots are given for n = 50 by Figure A.12 in Appendix A. Note that NLL is plotted in logarithmic scale (hence the legend logNLL). First of all, the NLL and especially its optimal values are closely correlated to those of the predictivity coefficient Q 2 , whatever the covariance considered. Minimizing the NLL for the estimation of the GP hyperparameters is therefore a good way to control the quality of the metamodel predictor and to ensure a good predictivity. On the other hand, the optimal points of the NLL do not correspond to the optimal points of the two other criteria controlling the quality of the whole predictive distribution, namely PVA and IAEα. Their behaviors w.r.t. the hyperparameters is relatively similar and more irregular than NLL and Q 2 , with notably different zones with local minima. Moreover, the optimal points of PVA and IAEα do not correspond to those of NLL and Q 2 . It would therefore be unwise to consider only PVA or IAEα in the optimization of hyperparameters, without first ensuring that the NLL is close to its optimal values. This example is representative of the close links often observed between NLL and Q 2 , and between PVA and IAEα. It also illustrates that these same two groups can be antagonistic, particularly in the case of a rough covariance (exponential for instance), while for smoother covariances, some interesting compromises between the two groups can be found. For instance, for Gaussian covariance, without degrading Q 2 too much (by staying on the same level line for it), the second group of criteria can be improved. This concordance between the criteria, which is observed all the more when the covariance is regular, makes sense here, given the high regularity of the Branin function. In the case of a well-specified covariance, the NLL optimization could thus be relevantly completed by an additional optimization of PVA or IAEα (under the constraint, for example, of not degrading the Q 2 ).

Illustration on a non-differentiable function

This section considers a more irregular function (non-differentiable and non-monotonic relationships), namely the G-Sobol function, defined in dimension d for uniform independent inputs on [0, 1] by:

M Sobol (X) = d k=1 g k (X k ) where g k (X k ) = |4X k -2| + a k 1 + a k and a k ≥ 0. ( 8
)
The same procedure as for Branin function is applied, with d = 2, a k = k for k = 1, 2 and n = 50. An illustration representative of the results obtained is given by Figure 2. Similar results are found, with agreement between NLL and Q 2 and between PVA and IAEα, but here, there is more compatibility between the two groups in particular for covariance 3/2-Matérn. For more regular covariance, it's more difficult to optimize both groups simultaneously. However, for Gaussian covariance, the optimum area of PVA and IAEα corresponds to a local optimum zone of Q 2 and NLL, whose value is close to that of the global optimum. On this example, it could be interesting to shift slightly from the optimal point of the NLL towards more optimal values for the IAEα under the constraint of remaining in a region of Q 2 values close to the optimal Q 2 .

Concluding remarks

We have carried out several other simulations on different analytical functions and for different dimensions, which are not presented here for the sake of brevity. In summary, similar observations often emerge with:

• a close behavior between NLL and Q 2 which pleads in favor of keeping NLL as the main estimation method to control the predictivity of the metamodel, which is consistent with the results obtained by Petit et al. [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF];

• PVA and IAEα validation criteria often have a similar and more irregular behavior w.r.t. the hyperparameter values, and according to the covariance structure. They often present areas with local minima: some of them can be very compatible with rather optimal values of NLL and Q 2 , while other minima correspond to very degraded values of NLL and Q 2 . It would therefore be unwise to optimize only PVA or IAEα, without considering NLL or Q 2 ;

• In many examples, it appears possible to find in the neighborhood of the optimal NLL point, a better point w.r.t. PVA and IAEα. But this second optimization must be done while controlling the possible degradation of Q 2 value, which is the first guarantee of the good prediction capacity of the metamodel.

In addition, we also studied the shape of the Pareto fronts of the different pairs of criteria, considering their values computed either on a test basis (to have the "true" value), or by cross validation (LOO estimators) to be more representative of the results of a multi-objective optimization that would be carried out on the learning sample. In conclusion, it first emerges that the priority must be given to the minimization of the NLL to have accurate GP predictions. In addition, Q 2 and IAEα appear as complementary criteria but without redundancy in the information they convey: the Q 2 validating the predictive mean of the GP and the IAEα validating the quality of the prediction intervals if the mean of the interval has been previously validated by the Q 2 .

The IAEα taken alone does not constitute a validation of the predictive distribution, it is only under condition of a good Q 2 . All these considerations have guided the choice of the constrained multi-objective algorithm proposed in the next section.

New algorithm based on MLE and IAEα under constraint on Q 2

Description of the algorithm

In line with the reparametrization considered by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] and based on our expertise acquired on the optimization of GP hyperparameters, the inverse reparametrization is used in the rest of the document, namely θi = 1 θ i for i = 1, . . . , d. Even if the MLE is invariant under any injective reparameterization, the inverse reparametrization impacts the Jacobian of the likelihood function and we empirically observe a better behaviour of the gradient-based algorithms usually used for MLE (such as Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm). However, the question of the practical impact of certain standard reparametrizations (e.g. inverse or log) would benefit from further investigation. In the same vein, the set of input variables are systematically re-scaled on [0, 1] to allow an homogeneous interpretation of bounds, initial or estimated values of the hyperparameters. Hence, from the previous analysis, our new algorithm for optimizing the GP hyperparameters is summarized by Figure 3 and based on the following three main steps:

▶ Step 1: initial optimization based on NLL. The aim is to have an initial set of hyperparameters provided by usual estimation methods, which will be a starting point for the multi-objective algorithm, and, above all, a reference value for GP predictivity Q 2 . So, a first set of GP hyperparameters denoted θinit M LE is obtained by minimizing the NLL of the learning sample. For this, a multistart procedure is used combined with a BFGS algorithm, on the bounded domain of hyperparameters D θ. The LOO estimator of Q 2 corresponding to θinit

M LE
is computed and denoted Q 2,init LOO . Note that the BFGS algorithm is one of the most widely used quasi-Newton method for solving unconstrained nonlinear optimization problems.

▶ Step 2: constrained multi-objective optimization based on NLL and LOO-IAEα.

A second optimization is carried out by considering now two objectives: the NLL and IAEα computed by LOO on the learning sample (using Dubrule's formula [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF]). The aim is to find a new set of hyperparameters that would improve the accuracy of GP prediction intervals but without significantly degrading its predictivity. To control this last point, a constraint on Q 2 (computed by LOO) is added into the optimization problem: the respect of a minimum threshold value of Q 2 denoted c Q 2 , is considered. Remember that Q 2 can be interpreted as the proportion of the output variance explained by the metamodel, it is therefore reasonable to assume that a reduction of 5% of the variance explained, so a reduction of Q 2,init LOO of 0.05, is a relevant limit. This limit allows to obtain a metamodel with similar predictive power to those of Step 1, while allowing the exploration of θ values enabling a potential and substantial improvement in the quality of prediction intervals. In practice, the simplest way is to express this limit as c

Q 2 = Q 2,init
LOO -γ with a threshold parameter γ ∈ [0, 1], γ = 0.05 being the recommended default value. Note that it could also be expressed relatively to

Q 2,init LOO by c Q 2 = γ Q 2,init LOO , with γ ∈ [0, 1]
, and default values γ = 0.95 or γ = 0.9 for instance. The constraint is then expressed for any new candidate

θnew by Q 2,new ≥ c Q 2 with Q 2,new
the Q 2 associated to the GP of hyperparameters θnew . To perform this optimization, a constrained multi-objective evolutionary algorithm, namely the constrained NSGA-II (Nondominated Sorting Genetic Algorithm), is used. Proposed by Deb et al. [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-ii[END_REF], it allows to reduce the computational complexity, and alleviate the non-elitism approach (of classic genetic algorithm) by using a modified mating and survival selection. Moreover, to ensure a search in the neighborhood of the point estimated by NLL at Step 1, the initial population is built upon θinit M LE . More precisely, half of the individuals are selected by randomly drawing a relative perturbation of ±10% of the values of θinit M LE . The other half is chosen via a space-filling design on D θ, namely a maximin Latin Hypercube Sampling (LHS), for more exploratory purposes [START_REF] Morris | Exploratory designs for computationnal experiments[END_REF][START_REF] Fang | Design and Modeling for Computer Experiments[END_REF].

▶ Step 3: Selection of the best individual in the Pareto front. Iterative evolution of a population of candidate solutions in NSGA-II algorithm results in a Pareto front. This is a set of so-called non-dominated solutions, i.e. offering the best possible trade-offs between the two objectives (here NLL and LOO-IAEα). The aim of step 3 is then to choose the best solution from this set, in relation to our problem. For this, a classification algorithm is first used to bring out homogeneous groups of solutions: a k-means algorithm is carried out where the optimal number of clusters is determined by the Elbow method [START_REF] Thorndike | Who belongs in the family?[END_REF][START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF]. A fundamental step for any unsupervised algorithm is to determine the optimal number of clusters into which the data may be clustered. The Elbow method is one of the most popular heuristic to determine this optimal value k to ensure a trade-off between the complexity of classification model (more clusters to capture more variance) and avoiding over-fitting the data (fewer clusters). To achieve it, the Elbow method computes the within-cluster sum of squares (WCSS, also called intra-cluster variance) as a function of k. The optimal value of k is then identified as the "elbow" point where WCSS starts to level off, i.e. where adding more clusters does not significantly reduce WCSS. Other methods such as silhouette method or Gap statistic-based approach [START_REF] Maechler | Cluster: cluster analysis basics and extensions[END_REF] could also be used to determine the optimal value of k. They could be substituted for the Elbow method without any further changes to our algorithm.

Then, several "best" candidates could be chosen: a prudent approach (w.r.t. predictivity) is first to select the cluster (denoted C1) with minimal value of the 1 st objective NLL and then, to select the best individual in C1 according to the 2 nd objective IAEα. This choice, which we recommend, results in the solution denoted θBestC1 N SGA in Figure 3. Other choices could be made, such as selecting the C1 cluster's centroid (which is more conservative in terms of the potential improvement in IAEα), or on the contrary, selecting the best individual for IAEα among all the points of the Pareto Front. This latter choice is more penalizing for Q 2 degradation.

In all the numerical experiments presented in the following, the algorithmic parameters of the different steps are set as follows:

• Step 1: D θ = [0.01; 10] d and a BFGS algorithm is used with a multistart approach of 10 points chosen by maximin LHS on D θ. An additional central initial point is considered with all the hyperparameters set at θ = 2.

• Step 2: for the NSGA-II algorithm the crossover and mutation percentages are both set at 50%, the mutation rate and step size are respectively µ NSGA = 0.02 and σ NSGA = 0.1∆ θ with ∆ θ the width of the definition interval of θ (same values are considered for all the ( θi ) 1≤i≤d ). Moreover, the total size of the initial population of NSGA-II algorithm is 80 (40 from random perturbation of θinit M LE an 40 from maximin LHS design on D θ). The maximum number of iterations is set at 50. The rest of NSGA-II parameters are set at their default values (see

Table B.2 in Appendix B). Concerning the constraint, it is here defined as

c Q 2 = Q 2,init
LOO -0.05. This parameter is fixed just for the automation of the procedure. In practice, this choice is really to be made by the user w.r.t. the initial value Q 2,init LOO and the acceptable degradation of Q 2 .

• Step 3: usual default values of the elbow algorithm are used. More precisely, the maximum number of clusters to try is fixed at ⌈ √ m⌉ with m the number of individuals in the Pareto front. The optimal number of clusters is selected so that a fraction equal to 0.95% of the variance is explained and the k-means procedure is repeated three times, taking the best result at the end. For all the numerical tests in the following, all the parameters of the optimization algorithm are fixed, independently of the model M, the dimension d and the sample size n. Thus, the results obtained will be less dependent on these choices. It aims to illustrate that a use with standard values is quite possible, without real optimization of the algorithm parameters. It is obvious that in practice a more expert choice is be possible and will improve the results.

Analysis of the algorithm computational cost

It should be noted that the proposed algorithm has a higher computational cost than standard MLE optimization approaches such as a simple BGFS or a multistart-BFGS. More precisely, each likelihood evaluation cost is in O (n 3 ). If we add the gradients with respect to θ ∈ R d , the cost is in O (n 3 + dn 2 ) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. The cost of fast formulas (Dubrule's formulas [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF]) for LOO cross-validation indicators Q 2 and IAEα is in O (n 3 ) [START_REF] Petit | Improved Gaussian process modeling: Application to Bayesian optimization[END_REF]. As a result, the complexities (w.r.t. n) of one iteration and one candidate point are therefore similar for both algorithms.

However, the NSGA-II evolutionary algorithm will perform many more iterations and evaluate many more candidates in practice, due to the nature of genetic algorithms. Moreover, it doesn't use the gradients of objectives and constraints to converge faster (unlike BFGS). The total computational cost of each algorithm obviously depends on the user's choice of optimization parameters (kind and value of stopping criterion, number of re-starts for multistart BFGS, size of initial population, mutation percentage, rate or step size for NSGA-II). For information purposes only, with the optimization parameters previously stated and fixed for the whole benchmark, this translates into a factor ranging from 10 to 20, for learning sample sizes of a few hundreds and less than ten of inputs. For instance, on a standard machine (Processor Intel Xeon W-10855M CPU 2.80GHz and RAM 32,0 Go), the CPU cost of the estimation process between multistart BFGS and our algorithm increases from 4s to 45s (resp. from 50s to 800s) for d = 5 inputs and learning sample of size n = 100 (resp. n = 500).

This additional cost is judged quite reasonable in our application context (learning samples of a few hundred simulations and a dozen inputs, typically), and especially with regard to the CPU-cost of each model evaluation. Remember that in the engineering applications considered and targeted by this work, each evaluation of the simulator (or model) can range from a few tens of minutes to several hours, or even several days of calculation for some simulators (CFD calculations for instance). The benefit-cost ratio of the new procedure is therefore favorable in view of the stakes involved and for the applications addressed. It becomes less so for large sample sizes (typically n > 1000). But all this could be considerably mitigated by parallelizing and distributing tasks over several processors, as genetic algorithms are particularly well-suited to this process.

Numerical benchmark on analytical functions

This section focuses on different analytical functions usually used in the literature on emulation of computer experiments, from dimension d = 3 to d = 20. All the functions have been redefined for uniform independent inputs on [0, 1] and the learning sample sizes n have been chosen to obtain from medium to high GP predictivities: Q 2 around 0.7 and Q 2 ≥ 0.9, respectively. Hence, a large number of functions have been tested, and a representative sample is given below for:

• Friedman function in dimension d = 5 proposed by Friedman [START_REF] Friedman | Multivariate Adaptive Regression Splines[END_REF]:

M Fried-d5 (X) = 10 sin(πX 1 X 2 ) + 20 (X 3 -0.5) 2 + 10X 4 + 5X 5 ; (9)
• Dette & Pepelyshev function in dimension d = 8 defined by Dette and Pepelyshev [START_REF] Dette | Generalized latin hypercube design for computer experiments[END_REF]:

M Dette-d8 (X) = 4 X 1 -2 + 8X 2 -8X 2 2 2 +(3 -4X 2 ) 2 +16 X 3 + 1 (2X 3 -1) 2 + 8 i=4 i ln   1 + i j=3 X j   ; (10)
• a function in dimension d = 20 proposed by Marrel et al. [START_REF] Marrel | The ICSCREAM methodology: Identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics[END_REF] and inspired from the Friedman function:

M M arrel-d20 (X) = a 1 sin 6π X 5/2
1 (X 2 -0.5) +a 2 (X 3 -0.5) 2 +a 3 X 4 +a 4 X 5 +r X6,...,X 15 [START_REF] Ma | Reliability-based design optimization using adaptive kriging-a single-loop strategy and a double-loop one[END_REF] where r X 6 ,...,X 15 = We also consider the G-Sobol (introduced in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF]), Ishigami and Becker's functions for dimension d = 8, d = 3 and d = 9, respectively. Details on the two last functions and associated results are given in Appendix C, for the sake of brevity.

Results without additional nugget effect

No nugget effect is considered at first and a constant mean m(x) = β 0 is assumed for all the GP. The learning samples are randomly generated according to space-filling LHS and, for each configuration, all the procedure is repeated 100 times by generating independent LHS designs. Hence, for each test case, three methods for estimating the GP hyperparameters are compared:

• a simple BLFG algorithm, as a very rudimentary estimate (not to be used in practice).

• A "multistart" BFGS (denoted "Multi-BFGS "), as a commonly used reference method offering in general a good compromise between simplicity (theoretical and in its practical use), efficiency and execution time.

• Our constrained multi-objective algorithm denoted "C-NSGA-II-BestC1 ": the best individual (in the Pareto front) according to the IAEα is chosen in the cluster denoted C 1 with the minimal mean value of the NLL (see the description of Step 3 in Section 4).

For each configuration, validation criteria Q 2 , PVA and IAEα, are computed on a random test sample of 10 4 simulations (independent from the learning sample). The results obtained are given in boxplot form by Figures 4 to 6, for the different functions ordered by increasing input dimension. For each function, different covariance functions are presented by order of increasing regularity, to facilitate interpretation. Thus, several observations can be made from these initial tests.

▶ Concerning the GP predictivity (Q 2 ), for all functions and whatever the covariance considered, the Multi-BFGS and C-NSGA-II-BestC1 algorithms are largely better in terms of predictivity than a simple BFGS which presents significantly lower and more variable Q 2 (especially when the sample size n is small). A bad estimation of the hyperparameters can thus be very penalizing for the quality of the metamodel. Then, C-NSGA-II-BestC1 does not yield a degradation of Q 2 w.r.t. Multi-BFGS. This is obviously explained by the candidate selected in the Pareto front in Step 3 of the algorithm (see 4): choosing the best candidate in C 1 being the most conservative and safe from the point of view of Q 2 . However, even with a threshold of 0.05 (on the Q 2 computed by LOO on the learning sample), the observed loss of Q 2 is often much lower or even almost zero. Note that other possibilities for the choice of best candidate in Pareto front have also been tested: the centroid of cluster C 1 gives equivalent results (sometimes a little worse) while the choice of the optimal point in terms of IAEα on the whole Pareto Front causes too much degradation of the Q 2 on average. ▶ Regarding now the accuracy of the GP predictive distribution, C-NSGA-II-BestC1 makes it possible to improve PVA and IAEα criteria, almost all the time and sometimes in a very significant way. The improvement brought by the algorithm is even more interesting in the following cases:

-If the Q 2 is very good and/or n is very large compared to the complexity of the function.

The NLL function is probably very flat in a rather large area of the hyperparameters. As a result, its optimization is not discriminative enough, many values of the hyperparameters leading to interesting NLL (and Q 2 ) values. Taking into account the IAEα criterion allows to focus on more interesting hyperparameters to obtain more reliable prediction intervals. The example of M Fried-d5 with 3/2-Matérn covariance (Figure 5 plot(a)) is a good illustration.

-If the model is misspecified, i.e. if the covariance function does not match the regularity of the function. This is what is observed for example for M Fried-d5 with 3/2-Matérn covariance (Figure 5

Results with an additional estimated nugget effect

We now consider an additional homoscedastic nugget effect in the covariance in order to evaluate the robustness of the method with this additional parameter to estimate. This nugget effect is included in the proposed algorithm in the same way as the other hyperparameters: the nugget parameter λ is also estimated in the multi-objective optimization. Its variation domain is D λ = [10 -3 ; 0.5] with an initial value set at λ = 0.3. A few results are given by Figures 7 to 9. As in the case without the nugget effect, a general improvement brought by the C-NSGA-II-BestC1 algorithm on PVA and IAEα criteria is observed, with no degradation of Q 2 . In large dimension, the contribution of the constrained multi-objective algorithm is all the more significant on the reliability of the GP predictive law, provided that the metamodel is already sufficiently predictive (e.g. Q 2 > 0.8), as it can be observed for M Marrel-d20 from n = 200 (see Figure 9). Note that on this last function, we observe a re-increase of PVA and IAEα indicators when the model becomes highly predictive (Q 2 > 0.95 for n = 400), this increase being mitigated by the use of C-NSGA-II-BestC1. In this case, the GP metamodel predicts the function almost perfectly but the width of the prediction intervals does not decrease quickly enough with n and they become too conservative.

Finally, similar benchmarks have also been performed considering simple Monte Carlo designs for the learning samples (instead of space-filling LHS designs). An extract of the results obtained is given in Appendix D. Overall, the results are close to those obtained with LHS designs, but more variable, due to the more variable sampling and, above all, less efficient for metamodeling purpose. Similar conclusions are drawn about the relative performance of the methods: the constrained multi-objective approach again gives better results (larger Q 2 and more reliable prediction intervals) and is significantly more robust to the sampling variance.

Application to an aquatic ecosystem model

In this section, a test case modeling a prey-predator chain in an aquatic ecosystem is studied. This so-called model MELODY (for MESocosm structure and functioning for representing LOtic DYnamic ecosystems) simulates the functioning of aquatic mesocosms and the impact of toxic substances on the dynamics of populations. Inside the model, two compartments linked within a prey-predator chain are considered: the Periphyton and the Grazers. The Periphyton-Grazers sub-model is representative of processes involved in dynamics of primary producers and primary consumers, i.e. photosynthesis, excretion, respiration, egestion, mortality, sloughing and predation. More details are available in Ciric et al. [START_REF] Ciric | Use of sensitivity analysis to discriminate non-influential and influential parameters within an aquatic ecosystem model[END_REF] and Iooss et al. [START_REF] Iooss | Some new insights in derivative-based global sensitivity measures[END_REF]. In this test case, a total number of d = 20 uncertain and independent input parameters are considered. These parameters characterize, among others, the photosynthesis, consumption, respiration, mortality or excretion rates of both populations. In absence of any expert opinion or relevant information on their uncertainty, these parameters are assumed to follow uniform distributions whose variation intervals are given in Ciric et al. [START_REF] Ciric | Use of sensitivity analysis to discriminate non-influential and influential parameters within an aquatic ecosystem model[END_REF]. As output from the model, we focus on the Grazers biomass at a given reference time (day 60 of simulations), denoted Y G .

A sample of n = 100 simulations of the model MELODY is available, drawn from a LHS with low discrepancy [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF]. Previous sensitivity analysis studies [START_REF] Roustant | Parseval inequalities and lower bounds for variance-based sensitivity indices[END_REF] have revealed the presence of strong non-linear and interaction effects for Y G . Considering the inputs' high dimensionality, these outputs are consequently very complex to emulate, making this test case relevant to evaluate the contribution of the multi-objective algorithm.

Description of the dataset and lognormal-kriging approach

The highly dispersed nature of the data requires to apply a preliminar logarithmic transformation to represent them. The histograms of the obtained values, denoted Z G = log(Y G ), is given by Figure 10. A kernel density estimation plot is also added to provide a graphical illustration of the probability density function. This logarithmic transformation will also be considered for the metamodeling, considering a lognormal-kriging approach to return to the initial space of the output values [START_REF] Journel | The lognormal approach to predicting local distributions of selective mining unit grades[END_REF][START_REF] Cressie | Lognormal kriging: Bias adjustment and kriging variances[END_REF]. Note that a more general Box-Cox transformation could also be considered. Hence, a predictive Gaussian distribution N (ẑ G (x), ŝ2

z G (x)) is obtained for Z G (see Eqs ( 4) and ( 5)). Lognormal-kriging backtransformations [START_REF] Cressie | Lognormal kriging: Bias adjustment and kriging variances[END_REF] are then used to obtain the predictor and associated variance prediction for the original data Y G :

ŷG (x) = e (ẑG(x)+0.5ŝ 2 z G (x)) (12) ŝ2 Y (x) = e ŝ2 z G (x) -1 e (2ẑG(x)+ŝ 2 z G (x)) . ( 13
)

Results of Gaussian process metamodeling

The lognormal-kriging approach is applied to the dataset: a GP is fitted on Z G sample with either the usual multistart-BFGS or our new C-NSGA-II multi-objective algorithm (C-NSGA-II-BestC1 algorithm). Different covariance functions are considered, with or without an estimated nugget effect, and a constant mean is considered in all cases. Validation criteria are computed by LOO on the learning sample, for both Z G and Y G . It should be noted that the computation of criteria by cross-validation for small data results in a double loop cross-validation structure for the C-NSGA-II algorithm, since it already includes a cross-validation loop. However this internal loop remains largely tractable for n = 100 and thanks to Dubrule's formulas [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF]. The obtained results, with the most interesting highlighted, are given in Table 1 (columns "Multi-BFGS" and "New C-NSGA-II Algorithm" of Tables 1(a) and 1(b)).

First, for the transformed variable Z G , the predictivity is similar for the different metamodels with nearly 90% of output variance explained. Results differ when we go back to the variables Y G . This is explained by the expression of the log-kriging predictor (Eq. ( 12)) which depends on the prediction variance ŝ2

Z . The multi-objective algorithm, which has more reliable prediction variances for Z, has consequently a much better predictivity for Y , with Q 2 on average increased by 0.04. Log-kriging is a further illustration of the need for reliable predictive variances in GP regression. Moreover, C-NSGA-II algorithm also yields more reliable prediction intervals, with a IAEα significantly reduced, whatever the covariance considered and whether or not there is a nugget effect. Note that for the IAEα criterion, the results are identical in both spaces, the empirical coverage function in GP framework being invariant by strictly monotonic transformation of the output variable. Note that IAEα is invariant by logarithmic transformation and more generally by strictly monotonic transformations since they preserve (for increasing transformations) or reverse (for decreasing transformations) the order of data. For example, in the case of the logarithmic transformation, the quantiles of the predictive law for Y are equal to the exponential of the quantiles of the GP predictive law obtained for Z G = log(Y G ). The values of ∆(α) computed with the set of z G or y G = e z G are therefore equal.

Some diagnostic plots are also given by Figure 11 with 5/2-Matérn covariance, to illustrate more clearly the improvement brought by the C-NSGA-II algorithm on the accuracy of prediction intervals. As an indication of the degree of non-linearity of the output, it may be mentioned that a linear regression with an ElasticNet-type penalization [START_REF] Hastie | The elements of statistical learning[END_REF]

Comparison with the RobustGaSP simplified Bayesian approach

To the best of our knowledge, and as detailed in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF], the RobustGaSP approach proposed by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] is the most interesting existing Bayesian approach to perform robust hyperparameter estimation, while being tractable in large dimension. To achieve that, RobustGaSP relies on three key ingredients. First, specific priors for the hyperparameters ensures theoretical guarantees on the estimation of hyperparameter posterior distribution. Second, some approximations simplify the computation of GP posterior distribution and avoid costly Markov So, too carry out the comparison with our approach, the R package RobustGaSP [START_REF] Gu | RobustGaSP: Robust Gaussian Stochastic Process Emulation[END_REF] is used. The jointly robust prior is considered to efficiently approximate the reference prior (argument prior choice set at ref approx in the rgasp function). In addition, we have also combined our multi-objective algorithm with RobustGaSP. To do so, the new algorithm is adapted by, replacing the NLL in the multi-objective procedure, by the marginal posterior of θ proposed by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] with their jointly robust priors.The obtained results for RobustGaSP approach and C-NSGA-II combined with RobustGaSP are given by Table 1 in columns 'RobustGaSP" and "New C-NSGA-II with Robust NLL", respectively.

If a nugget effect is considered and estimated, the C-NSGA-II algorithm gives significantly better results than RobustGaSP, for both Q 2 and IAEα criteria. Best results are obtained with 5/2-Matérn covariance. If no nugget effect is considered, RobustGaSP performs much better for the smoother covariances, namely 5/2-Matérn and Gaussian ones. Compared with C-NSGA-II algorithm, it yields an equivalent Q 2 for 5/2-Matérn covariance but much better Q 2 for Gaussian covariance. But in any case, it performs less well in terms of prediction interval reliability with a IAEα two to three times higher. RobustGasp results are therefore very sensitive to GP specifications (covariance function and nugget), whereas our C-NSGA-II algorithm is more robust to these choices. If we finally look at the combination of both approaches in the so-called "New C-NSGA-II with Robust NLL" algorithm, its results are promising: it seems to be a good option since the results of C-NSGA-II algorithm are systematically improved. Moreover, it leads to the best overall result with 5/2-Matérn covariance and nugget effect (Q 2 = 0.88 and IAEα = 0.02).

Conclusion

In the framework of emulation of numerical simulators with GP regression, this work has introduced a new algorithm for the estimation of GP covariance parameters, referred to as GP hyperparameters. Since the use of the entire GP predictive distribution is one of the major assets of this metamodel, a robust estimation of the hyperparameters is a cornerstone of successful metamodeling. For this purpose, and convinced of the value of going beyond the usual maximum likelihood approach, a new constrained multi-objective algorithm has been proposed. Based on a thorough analysis of the links between different estimation and validation criteria, this algorithm consists in jointly maximizing the likelihood of the observations as well as the empirical coverage function of GP prediction intervals computed by a leave-one-out (LOO) procedure, under the constraint of not degrading the GP predictivity. Particular care has been taken to detail and justify the parametric choices of the algorithm to facilitate its implementation and favor reproducible results.

A large benchmark on analytical functions of variable input dimension (d = 1 to 20) has been performed, considering different designs of experiments and different covariance models. The new algorithm was compared to usual simple and multi-start quasi-Newton algorithms (which focuses only on the likelihood optimization). The constrained multi-objective algorithm provides better results in terms of predictivity and reliability of prediction intervals, and is much more robust to the sampling variance. The improvement brought by the algorithm is all the more interesting when the covariance model is misspecified, when the number of GP hyperparameters is large or even when the size of the data is large.

The application relevance of this algorithm has also been shown on a real test case modeling an aquatic ecosystem (prey-predator chain). The GP metamodeling is used to predict the biomass of the two species at a given time. Once again, the multi-objective algorithm performs better than standard algorithms. Furthermore, the log-kriging approach notably illustrates the need for well-estimated and reliable prediction variances in GP regression. Finally, these good results have been confirmed by a comparison with the RobustGaSP method [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF], which is to our knowledge the only efficient and tractable Bayesian method. The multi-objective algorithm achieves results at least as good as those of RobustGaSP, while being less sensitive to GP specifications (covariance function and nugget). The combination of both approaches in a modified version of our algorithm yields promising results, taking the best of both methods.

In conclusion, as widely illustrated, the proposed new algorithm enables a more robust estimation of GP hyperparameters, by penalizing, to a certain extent, the hyperparameter values that would not enable reliable prediction intervals. This "penalty" is conveyed by the training data itself, unlike a Bayesian approach where the penalty is conveyed by the prior law, the choice of which can have a major impact on the results and requires a certain amount of expertise. The approach is also simpler and more tractable in practice than a Bayesian approach.

Only some parameters, intrinsic to the evolutionary optimization algorithm, need to be chosen by the user. Although their impact is fairly understandable, and default choices have been proposed, their influence should be further studied and quantified. The approach also requires more CPU time than standard algorithms due to the higher number of likelihood function evaluations, but this cost remains completely tractable on datasets of a few hundred to a thousand observations. Beyond that, the cost-benefit ratio is no longer attractive, especially since, on a "real" dataset, rigorous and fair validation of the final metamodel calls for a nested cross-validation (double algorithmic loop), as validation criteria are used in the optimization algorithm. This means relaunching the optimization algorithm as many times as there are folds in the external cross-validation loop. This problem can easily be solved by parallelizing the procedure. As another perspective, it might also be worthwhile accelerating the procedure by focusing multi-objective optimization only on the main influential hyperparameters or by using the gradients of LOO-based criteria (taking inspiration from the work of Petit et al. [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF]). Finally, as discussed in the conclusion of the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation[END_REF], other research directions for GP concern its estimation for very high dimensional problems (for example beyond 30 to 50 input variables), its extension to more complex inputs (graph structure for instance), and the learning of outputs with highly irregular or even chaotic behavior (due for example to physical threshold phenomena and phenomenological bifurcations).

Beyond these algorithmic refinements, the main objective of this work was to demonstrate the value of going beyond standard estimation approaches by combining several criteria, and of taking particular care to control its quality in fine. GP metamodel is a powerful tool but its implementation requires a certain expertise which can be a limiting factor in practice for its automation and its use on complex industrial cases. While this paper attempts to provide guidance, the fact remains that GP validation also requires careful and informed consideration of the topic, to ensure confidence in its use for predictive purposes.

Appendix C.3. Results for Becker's function in dimension d = 20

We also considered the Becker's function in dimension d = 20 proposed by Becker [START_REF] Becker | Metafunctions for benchmarking in sensitivity analysis[END_REF]. This function relies on a list of univariate basis functions that may possibly appear, representing response features observed in physical models. This list of basis functions is the following:

f 1 (x) = x (linear) f 2 (x) = x 2
(quadratic) f 3 (x) = x 3 (cubic) f 4 (x) = (e x -1) /(e -1) (exponential) f 5 (x) = 1 2 sin(2πx) + 1 2 (periodic)

f 6 (x) = 1 if x > 1 2
and 0 otherwise (discontinuity) f 7 (x) = 0 (no effect)

f 8 (x) = 4 x - 1 2 2
(quadratic, non-monotonic) f 9 (x) = (10 -1/1.1) -1 (x + 0.1) -1 -0.1 (inverse with small shift)

where each basis function has been scaled so that inputs in [0, 1] also map to outputs in [0, 1] (see Figure of Becker [44]). Given these inputs, the function M Becker is built as a sum of the main effects and interactions:

M Becker (x, u, V, w, Θ) = d i=1 a i f u i (x i ) + d 2 i=1 b i f u V i,1 x V i,1 f u V i,2 x V i,2 + d 3 i=1 c i f u W i,1 x W i,1 f u W i,2 x W i,2 f u W i,3 x W i,3 ,
where Θ = {a, b, c} are weighting coefficients applied to the main effect, second and third-order interaction terms respectively. One basis function is therefore specified for each input variable as its main effect. The interaction terms are generated as mixtures of these main effect terms. Note that this restriction on the construction of the interaction terms is done to avoid additional parameters, but it could be set otherwise.

To cover a large spectrum of models encountered in physical applications, random functions are then generated from the definition of M Becker by sampling the parameters of {u, V, W, Θ} from an appropriate distribution. More precisely, as suggested by Becker [START_REF] Becker | Metafunctions for benchmarking in sensitivity analysis[END_REF], u and V are sampled independently from uniform discrete distributions on 1, d . Each parameter of Θ is sampled independently from a mixture of two zero-mean Gaussian distributions: one with a low variance equals to 0.5, and another with a high variance equals to 5. The mixture parameters are set at 0.7 and 0.3, respectively for the two distributions. Note that for our numerical benchmark, we only considered functions with interaction effects of second order at most (i.e., we have imposed that

Figure 1 :

 1 Figure 1: MBranin Function -Comparison of NLL (computed on the learning sample) and validation criteria (computed on a test sample), for a GP built on a Monte Carlo learning sample of size n = 30. The optimal value for each quantity is indicated by a red square.

Figure 2 :

 2 Figure 2: M Sobol Function -Comparison of NLL (computed on the learning sample) and validation criteria (computed on a test sample), for a GP built on a Monte Carlo learning sample of size n = 50. The optimal value for each quantity is indicated by a red square.

Figure 3 :

 3 Figure 3: Flowchart of multi-objective constrained algorithm.

a 5 √(

 5 i=6...15 i 2 ) i=6...15 √ 12i(X i -0.5). More details and illustration of main effects of M M arrel-d20 are given in Appendix C.1.

Figure 4 :

 4 Figure 4: M Dette-d8 Function -Evolution of validation criteria, according to sample size n, for different hyperparameter estimation methods (GP with different covariances and all without nugget effect).

Figure 5 :Figure 6 :

 56 Figure 5: M Fried-d5 Function -Evolution of validation criteria, according to sample size n, for different hyperparameter estimation methods (GP with different covariances and all without nugget effect).(a) 3/2-Matérn covariance

-

  plot (a)) or for the infinitely differentiable M Ishig-d3 function with a 3/2-Matérn covariance whose trajectories are only once differentiable (Figure C.14 plot (a), in Appendix C.2). In this particular case, the C-NSGA-II-BestC1 algorithm allows to divide by two the IAEα for n = 150. When the dimension of the inputs d is large (and consequently the number of hyperparameters in the case of a tensorized anisotropic stationary covariance). The C-NSGA-II-BestC1 algorithm provides significantly better and more robust results, less sensitive to sampling. This is what is observed for M Marrel-d20 in dimension d = 20 and M Becker-d9 in dimension d = 9 (see Appendix C.3).

Figure 7 :

 7 Figure 7: M Fried-d5 Function -Evolution of validation criteria, according to sample size n, for different hyperparameter estimation methods (GP with different covariances and all with an estimated nugget effect).

Figure 8 :

 8 Figure 8: M Dette-d8 Function -Evolution of validation criteria, according to sample size n, for different hyperparameter estimation methods (GP with different covariances and all with an estimated nugget effect).

Figure 9 :

 9 Figure 9: M Marrel-d20 Function -Evolution of validation criteria, according to sample size n, for different hyperparameter estimation methods (GP with different covariances and all with an estimated nugget effect).

Figure 10 :

 10 Figure 10: MELODY test case -Histogram of the output Y G after log-transformation, for the n = 100-size dataset.

  leads to Q 2 = 0.57 for Z G and to zero Q 2 if directly fitted on Y G . (a) With nugget parameter (estimated) (b) Without nugget parameter

Table 1 :

 1 MELODY Output Y G -Validation criteria computed by cross validation for GP with different covariances and estimation methods, without and with (estimated) nugget effect. The significantly better (resp. best) results are framed in black (resp. red).

Figure 11 :

 11 Figure 11: MELODY Output Y G -Diagnostic plots for GP metamodeling with 5/2-Matérn covariance and nugget effect: LOO predicted values versus exact values on the left (with log-scale) and α-plot on the right.

Figure C. 13 :

 13 Figure C.13: M Marrel-d20 -Illustration of the one-dimensional mean effects for the first five inputs, figure extracted from Marrel et al. [35].

Figure C. 14 :

 14 Figure C.14: M Ishig-d3 Function -Evolution of validation criteria, according to sample size n, for different hyperparameter estimation methods (GP with different covariances and without nugget effect).

Acknowledgments

This work was funded by the French ANR project SAMOURAI (ANR-20-CE46-0013). We warmly thank the associate editor and three reviewers whose remarks have greatly helped to improve the paper. We are also grateful to Yann Richet for fruitful discussions and help with implementation of RobustGaSP approach.

As in Section 3, a further illustration of the behavior of NLL and validation criteria for the "re-scaled" Branin function is given below by Figure A.12, for a dataset of size n = 50.

Appendix B. Parameters of NSGA-II algorithm

The default (and recommended) values used in all the paper for the NSGA-II algorithm are specified in

Appendix C. Benchmark of Section 5: details on test functions and additional results

Appendix C.1. Details on a modified version of Friedman function

The analytical model proposed by Marrel et al. [START_REF] Marrel | The ICSCREAM methodology: Identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics[END_REF] and inspired from the Friedman function [START_REF] Friedman | Multivariate Adaptive Regression Splines[END_REF] is defined in dimension d = 20 by:

where r X 6 ,...,X 15 = a 5 ( i=6...15 i 2 ) i=6... [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF] √ 12i(X i -0.5) and X = (X 1 , . . . , X 20 ) are independent and uniform random variables on [0, 1]. The model depends only on the 15 first inputs. The first term represents a strong and non monotonic interaction between the two first inputs. The second term is a quadratic function of X 3 while the other ones are linear. The parameters for tuning the influence of the different inputs are chosen as follows: a = (5, 20, 8, 5, 1.5) ⊤ . Under this parametrization, X 2 explains alone around 10% of the output variance, X 1 has no individual effect but its interaction effect with X 2 is strong (around 30% of the output variance). X 3 , X 4 and X 5 only have individual effects (no interaction) and explain, respectively, around 11%, 28% and 10% of the output variance. The effects of the ten remaining inputs (X 6 to X 15 ) represent around 11.5% of the output variance. The one-dimensional mean effects of the five first (and main influential) inputs are plotted on Figure (C.13).

Appendix C.2. Results for Ishighami function in dimension d = 3

Following a protocol strictly similar to that described in Section 5, tests were carried out on other analytic functions, among them the Ishighami and Becker's functions. The Ishighami function introduced by Ishigami and Homma [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF] is defined in dimension d = 3 by: