Amandine Marrel

Bertrand Iooss

Probabilistic surrogate modeling by Gaussian process: A review on recent insights in estimation and validation

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

In the framework of risk assessment, computer codes are increasingly used to understand, model and predict physical phenomena. As these codes can be very time-consuming to run, which severely limit the number of possible simulations, a widely accepted approach consists in approximating the CPU-time expensive computer model by a so-called "surrogate model". In this context, the Gaussian Process regression (also called kriging) is one of the most popular technique. It offers the advantage of providing a predictive distribution for all new evaluation points. An uncertainty associated with any quantity of interest (e.g. a probability of failure in reliability studies) to be estimated can thus be deduced and adaptive strategies for choosing new points to run with respect to this quantity can be developed. This paper focuses on the estimation of the Gaussian process covariance parameters by reviewing recent works on the analysis of the advantages and disadvantages of usual estimation methods, the most relevant validation criteria (for detecting poor estimation) and recent robust and corrective methods.

Introduction

In the framework of risk assessment, computer codes (or numerical simulators) are developed and increasingly used to understand, model and predict physical, engineering or biological phenomena [START_REF] Sokolowski | Modeling and Simulation Fundamentals: Theoretical Underpinnings and Practical Domains[END_REF]. They usually take a large number of input parameters driving the phenomenon of interest or related to its physical and numerical modeling. However, the available information about some of these parameters is often limited or uncertain. The uncertainties come mainly from the lack of knowledge about the underlying physics, the characterization of the input parameters of the model (e.g. due to the lack of experimental data) or to the choice of scenario parameters. Therefore, it is essential to take the uncertainties tainting the results of computer simulations into account in order to perform "Uncertainty Quantification" [START_REF] Smith | Uncertainty quantification[END_REF].

A probabilistic framework where the input uncertainties are modeled by fully or partially known probability distributions, based on available data, expert opinions or bibliographic databases is usually considered [START_REF]Uncertainty in industrial practice[END_REF]. The uncertainty quantification process therefore relies on Monte Carlo techniques: a sample of code simulations is performed, where the inputs are drawn according to their probabilistic distributions. Estimators of the target statistical quantities, also called statistical quantities of interest (e.g. the variance, a probability of exceeding a threshold, or some quantiles) are then computed from the sample of code outputs. Depending on the nature of the quantities to be estimated and the expected confidence in the estimators, a very large number of simulations of the code can be necessary: from a few hundred to several tens of thousands for example. This number can also depends on the dimension of the inputs, when performing sensitivity analysis for instance [START_REF] Veiga | Basics and Trends in Sensitivity Analysis[END_REF]. There are numerous examples of sensitivity analysis and uncertainty quantification based on probabilistic approaches, particularly in the case of risk assessment using environmental models (see, e.g., the motivating examples in the textbooks [START_REF]Uncertainty in industrial practice[END_REF][START_REF] Smith | Uncertainty quantification[END_REF][START_REF] Veiga | Basics and Trends in Sensitivity Analysis[END_REF]).

In this context, one key issue is that the numerical model under study can be very timeconsuming 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). These metamodels can be based on polynomials, splines, random forests, neural networks, etc. [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF][START_REF] Afshari | Machine learning-based methods in structural reliability: A review[END_REF], in fact on any machine learning techniques [START_REF] Hastie | The elements of statistical learning[END_REF]. Built from a set of computer code simulations, they must be as representative as possible of the code outputs in the domain of variation of the uncertain parameters while having good prediction capabilities. Nowadays, metamodels are extensively used in several engineering fields to solve industrial issues as it provides a multi-purpose tool [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF]: once fitted, the metamodel can be used, possibly in conjunction with the costly computer code, to perform sensitivity analysis, as well as uncertainty propagation, optimization, or calibration studies. Such techniques have been extensively developed for instance in nuclear engineering (see, e.g., [START_REF] Wu | Kriging-based inverse uncertainty quantification of nuclear fuel performance code BISON fission gas release model using time series measurement data[END_REF][START_REF] Marrel | The ICSCREAM methodology: Identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics[END_REF]). However, to be confident with this approximation-based approach in support of the different uncertainty quantification tasks, it is crucial to develop accurate and reliable metamodels to approximate the computer model.

Among the metamodels classically used for numerical experiments, the Gaussian Process (GP) regression, also called kriging model, is a popular tool for non-parametric function estimation. Historically introduced in the context of geostatistics for spatial interpolation (see, e.g., [START_REF] Chilès | Geostatistics: Modeling spatial uncertainty[END_REF]), GP regression has been extended to interpolation of numerical simulation outputs (see, e.g., [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]) and machine learning approximation (see, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]). Its intuitive idea is to start from a prior over random functions (a GP is characterized by its mean and covariance functions), then the GP regression yields a posterior over functions given the observed data. On one hand, this makes it a very flexible non-parametric regression tool, suitable for modeling of numerical simulators and whose effectiveness has been illustrated in many applications [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF]. On the other hand, it also offers a probabilistic framework: the GP metamodel yields a predictive distribution for the code output at each prediction point, with a simple analytic formulation. From this, a prediction but also an uncertainty via prediction intervals can be analytically derived.

Figure 1 illustrates this principle of using a (probabilistic) metamodel, by the way of the GP, to emulate computer code from a set of code simulations. The main idea is that, in addition to a predicted value for the quantity of interest, the stochastic information provided by the GP metamodel allows to obtain confidence bands on this quantity. In many engineering problems, and especially those related to safety, it is essential to have access to prediction intervals in addition to the predictions of machine learning algorithms, in order to understand and assess the degree of confidence in the algorithm's prediction. This is part of the more general framework called "Trustworthy Machine Learning" [START_REF] Varshney | Trustworthy Machine Learning[END_REF]. This major advantage, intrinsic to GP, is therefore highly appealing for risk and safety assessment applications. Its value has been illustrated, in simulationbased nuclear engineering and safety analyses, for estimating probabilities of failure [START_REF] Marrel | Sensitivity analysis of safety factor predictions for nuclear component behaviour under accidental conditions[END_REF], sensitivity indices [START_REF] Gratiet | A Bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF], excursion sets [START_REF] Chevalier | Fast krigingbased stepwise uncertainty reduction with application to the identification of an excursion set[END_REF], functional risk curves [START_REF] Iooss | Uncertainty and sensitivity analysis of functional risk curves based on Gaussian processes[END_REF], maximal values [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF] and high-order quantiles [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF].

Moreover, the GP-associated prediction intervals allow to develop so-called adaptive (also called "active learning" or "goal-oriented") strategies: the idea is to sequentially find, from the current set of simulations, a new set of points to run in the input space in order to most efficiently estimate a statistical quantity of interest (see, e.g., [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF]). This approach has been shown to drastically increase convergence with respect to standard Monte Carlo or quasi-Monte Carlo algorithms (see 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] for overviews and benchmarks of GP-based adaptive algorithms). Among the numerous algorithms that have been recently developed, one can distinguish those based on pointwise criteria, which use uniquely the conditional mean and variance of the GP at a given point (see, e.g., [START_REF] Echard | AK-MCS: An active learning reliability method combining kriging and Monte Carlo simulation[END_REF]) and those based on integral criteria, which integrate functions of the conditional mean and variance of the GP over the whole input domain (see, e.g., [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]). Among the different topics that use such algorithms, all recently published in the Reliability Engineering and System Safety journal, one can cite prediction-forecast [START_REF] Betancourt | Gaussian process metamodeling of functional-input code for coastal flood hazard assessment[END_REF], calibration [START_REF] Perrin | Adaptive calibration of a computer code with time-series output[END_REF], structural reliability [START_REF] Huang | A new active learning kriging metamodel for structural system reliability analysis with multiple failure modes[END_REF][START_REF] Afshari | Machine learning-based methods in structural reliability: A review[END_REF], reliability-based design optimization [START_REF] Ma | Reliability-based design optimization using adaptive kriging-a single-loop strategy and a double-loop one[END_REF] and robust optimization [START_REF] Ribaud | Robust optimization: A kriging-based multi-objective optimization approach[END_REF]. The main aim of this paper is to focus strong attention on a key point, often omitted in the above mentioned studies, in the practical building and use of GP. Indeed, when fitting a GP metamodel on a given dataset, the specification of the covariance structure of the process is particularly sensitive. The power of the GP lies to a large extent in the wise choice of a covariance model and a good estimate of its parameters (often called the hyperparameters). A poorly specified covariance affects the predictivity of the GP, but above all, and to a greater extent, the correctness of its confidence intervals. Although the confidence intervals' issue is at the core of this work, the validation metrics that will be presented and discussed attempt to validate GP predictivity and confidence intervals.

Concerning the choice of mean and covariance functions, they are usually chosen within parametric families (for instance, class of Matérn covariance functions). The estimation of GP then consists only in the estimation of hyperparameters. For this, different methods exist based on either likelihood maximization, cross-validation technique or Bayesian approach. From a theoretical point of view and under the hypothesis of a well-specified covariance model, some results exist concerning the consistency of the different types of estimators [START_REF] Zhang | Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics[END_REF][START_REF] Bachoc | Asymptotic analysis of maximum likelihood estimation of covariance parameters for gaussian processes: An introduction with proofs[END_REF][START_REF] Bachoc | Cross-validation estimation of covariance parameters under fixed-domain asymptotics[END_REF]. From a practical standpoint, some studies [START_REF] Sundararajan | Predictive approaches for choosing hyperparameters in Gaussian processes[END_REF][START_REF] Martin | Use of kriging models to approximate deterministic computer models[END_REF][START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF] have proposed comparisons on analytical functions, but no consensus really emerges: the estimation of hyperparameters is often unstable regardless of the method. Moreover, most of the time, the authors only focus on the accuracy of GP predictor to assess the impact of hyperparameter estimates. The reliability of prediction intervals is rarely considered (except in Petit et al. [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF]), while it is often affected by misspecification of GP covariance and poor estimation of hyperparameters (as highlighted by Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF]). To mitigate this problem, other authors such as Acharki et al. [START_REF] Acharki | Robust prediction interval estimation for Gaussian processes by cross-validation method[END_REF] propose to correct the hyperparameter estimates to obtain more robust and accurate prediction intervals (that are called reliable prediction intervals in the following). Whatever the estimation method or correction considered, this calls for validation indicators to control the performance and robustness of GP regression.

In this scope, this paper reviews recent works dealing with the difficulties inherent in estimating GP hyperparameters. The objectives are to analyze the advantages and disadvantages of estimation methods, to list and propose relevant validation criteria (to detect poor hyperparameter estimation), and to study some recently proposed robust and corrective methods. The rest of the document is organized as follows. Reminders on GP regression, parameterization and estimation of parameters are given in Section 2. Section 3 reviews the different hyperparameter estimation algorithms. Section 4 lists important criteria for GP validation. Section 5 then proposes a review of very recent papers dealing with the robust estimation of GP hyperparameters and allows to explain the orientation chosen for our research work, w.r.t. the application context. In particular, validation criteria to control the performance and robustness of GP regression are detailed. The last section gives some conclusions and prospects of this work. From this extensive review, a companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A new estimation algorithm for more robust prediction[END_REF] proposes a new algorithm that solves some of the identified drawbacks of the previous ones.

Reminders on Gaussian process regression

Throughout the rest of this paper, the numerical model (computer code or simulator) is represented by the following input-output relationship:

M : X -→ Y X -→ Y = M(X) (1)
where the uncertain output variable Y and the d input parameters X = (X 1 , . . . , X d ) ⊤ belong to some measurable spaces respectively denoted by Y and X ⊂ R d . As part of the probabilistic approach, the inputs are considered as random variables with probability distributions denoted by P X on X [START_REF] Smith | Uncertainty quantification[END_REF]. It is therefore assumed that we have a n-size 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 (also called the "experimental design"), and Y s = {y 1 , . . . , y n } denotes the corresponding outputs observations with y i = M(x (i) ). (X s , Y s ) constitutes the learning sample.

GP metamodel conditioned by the learning sample

In the GP regression framework [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF], the data are modeled as discrete observations of a GP sample path. The prior knowledge on observations is modeled by a GP completely specified by its mean function m

(x) = E[Y (x)] and its covariance function k(x, x) = COV(Y (x), Y (x)) = E[(Y (x) -m(x))(Y (x) -m(x))]. k(•,
•) is also called the covariance kernel and is assumed to be a positive definite kernel. The predictive GP distribution is therefore naturally given by the GP conditioned by the known observations Y s , denoted [Y (x)|Y(X s ) = Y s ]. Its distribution can be obtained analytically from the following joint distribution:

Y (x) Y(X s ) ∼ N m(x) m(X s ) , k(x, x) k(x, X s ) T k(x, X s ) K , (2)
where:

• Y(X s ) = (Y (x (i) )) 1≤i≤n ∈ R n is the vector of output value at sample locations,

• m(X s ) = (m(x (i) )) 1≤i≤n ∈ R n is the vector of mean function evaluated at sample locations,

• K = (k(x (i) , x (j) )) 1≤i,j≤n ∈ R n×n is the covariance matrix at sample locations,

• k(x, X s ) = (k(x, x (i) )) 1≤i≤n ∈ R n is the covariance vector between x and sample locations.

By applying the conditioning theorem of Gaussian vectors to the joint distribution, the conditional field Y (x)|Y(X s ) = Y s is still a GP whose mean is given by:

ŷ (x) = E[Y (x)|Y(X s ) = Y s ] = m(x) + k(x, X s ) T K -1 Y s -m(X s ) , (3)
and its covariance function:

ĉ(x, x) = COV[Y (x), Y (x)|Y(X s ) = Y s ] = k(x, x) -k(x, X s ) T K -1 k(x, X s ). (4)
In the following, the conditioning notation [•|Y(X s ) = Y s ] will be reduced to [•|Y s ] for the sake of brevity. Therefore, the predictive distribution for a new (unobserved) point x is the Gaussian distribution N (ŷ(x), ŝ2 (x)) where ŝ2 (x) = ĉ(x, x).

The conditional expectation ŷ(x) is used as the predictor of the GP regression model and its mean-square error is given by the conditional variance ŝ2 (x), while the Gaussian predictive distribution can be used to build predictive intervals of any level α ∈]0, 1[. More generally, conditional simulations (i.e. simulation of conditional GP trajectories) can be used to estimate, with a confidence interval, any statistical quantity of interest derived from the output (probability of exceeding a threshold, quantiles, etc.). The predictive distribution of some quantities of interest can also be defined analytically: this is obviously the case for a vector of prediction points, but also for derivatives or excursion sets. This possibility offered by GP regression is of particular interest in uncertainty quantification studies or for developing optimization strategies.

Covariance function, hyperparameters and nugget effect

The prior knowledge in GP regression consists in specifying the mean m(x) and the covariance function k(x, x) which is certainly the most important ingredient of a GP regression as it describes the dependence structure and controls the smoothness of the approximation.

Usual covariance functions and consideration on a priori choice

In the GP regression of computer experiments, the most popular choice is undoubtedly the class of stationary ν-Matérn functions defined in one dimension (x ∈ R) by:

k σ,ν,θ (x, x) = σ 2 2 1-ν Γ(ν) √ 2νh θ ν K ν √ 2νh θ , ( 5
)
where h = |x -x|, σ 2 and θ ∈ R + are respectively the variance parameter and the correlation hyperparameter (also called correlation length or length-scale). K ν is a modified Bessel function of second kind with parameter ν ∈ R + , and Γ is the Euler Gamma function. The parameter ν controls the smoothness of the GP: Y (x) is at least k-time mean-square differentiable if and only if ν > k. ν = 1/2 corresponds to the exponential covariance with continuous but not differentiable GP trajectories, while the limiting case ν → ∞ yields the Gaussian covariance function with infinitely differentiable trajectories. Between these two extreme cases, two popular ν-Matérn covariances in the statistical learning community are the ones for ν = 3/2 and ν = 5/2, which respectively yield to GP trajectories once and twice differentiable (see Table 1). Hence, choosing a correlation function most often consists in assuming a prior regularity for the model to be emulated since it directly defines the space of the possible trajectories in which the "real" function (or model) is supposed to belong. The success of the GP metamodel is conditioned to an adapted covariance model [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF]. To ensure a relevant choice, one solution might be to consider the family of ν-Matérn functions and to integrate the regularity parameter ν in the set of hyperparameters to be estimated from the dataset. But, as pointed out among others by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF], this is not a very relevant idea in practice notably for the emulation of computer experiments in high dimension. First, most often a space-filling design [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF] is used to generate X s in order to have simulation points with good input space covering. Although this type of design optimizes in practice the predictivity of the GP metamodel, the absence of points' aggregates prevents from properly characterizing and therefore estimating the ν parameter. Indeed, the covariance functions mainly differ at the origin (i.e. for a distance between points tending towards 0) and in particular by the decay rate at this point. Furthermore, it turns out that jointly estimate ν with (σ 2 , θ) may cause in practice computational and inferential difficulties in the estimation processes. It is therefore preferable (and commonly adopted) to estimate (σ 2 , θ), conditionally to a specified value of ν.

v = 1 2 v = 3 2 v = 5 2 v = +∞ Usual name exponential 3/2-Matérn 5/2-Matérn Gaussian k σ,ν,θ (x, x) σ 2 e -h Differentiability C 0 C 1 C 2 C ∞ of GP trajectories
Hence, as suggested by Petit [START_REF] Petit | Improved Gaussian process modeling: Application to Bayesian optimization[END_REF] and in the direct line of Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF], an interesting compromise (that we also recommend from our experience) is to consider a finite collection of covariance functions (those of Table 1), then estimate the hyperparameters (σ 2 , θ) for each of them, and finally use a validation criterion (different from criterion used for the estimation) to select the best covariance. It is also reasonable to consider the 5/2-Matérn covariance as the default practical choice because, as highlighted by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF], it has very interesting behavior of the 5/2-Matérn covariance w.r.t. to the distance between two input points. On the one hand, when distance tends towards zero, 5/2-Matérn covariance behaves like Gaussian covariance, maintaining the smoothness for nearby inputs while ensuring better conditioning number of the covariance matrix. On the other hand, when distance tends to infinity, it behaves like the exponential covariance, preventing from decreasing quickly with distance (as does the Gaussian correlation). This can be useful for sparse data (as is often the case in numerical simulator emulation) or for non-influential inputs for which it is logical that the correlation is quasi-constant with distance.

Finally, let us mention another approach that goes beyond the usual covariance kernels and consists in building a mixed covariance function in a single GP model. A new kernel function is constructed by mixing several usual kernels (such as those in Table 1) through operators such as addition and multiplication [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. However, such an approach complicates the hyperparameter estimation step as this might increase the total number of hyperparameters and add some weighting parameters which must also be estimated. It might be impractical for high-dimensional test cases and in case of insufficient number of data. Problems of hyperparameter identifiability might also be increased in practice. Nevertheless, a tractable process for learning composite kernels has recently been proposed in Palar et al. [START_REF] Palar | Gaussian process surrogate model with composite kernel learning for engineering design[END_REF]. The method is currently reserved for combinations of kernels sharing the same hyperparameter vector. However, from an applications point of view, it would be more valuable with different hyperparameters for each of the kernels in order to capture phenomena occurring at different scales. Unfortunately, this extension is not recommended by Palar et al. [START_REF] Palar | Gaussian process surrogate model with composite kernel learning for engineering design[END_REF], at least not with their algorithm, as it presents poorer results. More generally, even though this composite kernel approach doesn't seem to be widely used in practical studies, it clearly stands as a future avenue of research and improvement for GP metamodeling.

Extension to multivariate case

In order to extend to multi-dimensional inputs x ∈ R d , a widely used approach consists in considering a tensorized covariance defined as a product of univariate covariances:

k σ,ν,θ (x, x) = σ 2 d i=1 k 1,ν,θ i (x i -xi ). ( 6
)
The d 1-D covariance functions can be of different natures (with different smoothness parameters ν i for instance). 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.

Additional variance modeled by nugget effect

An additional nugget effect can also be considered in the covariance: it assumes an additive white noise effect, whose variance denoted σ 2 ϵ constitutes the nugget parameter. Most often, σ 2 ϵ is assumed to be constant, independent from the inputs (homoscedastic hypothesis). The covariance matrix then becomes K ′ = K + σ 2 ϵ I n where I n is the identity matrix. From a purely parametric point of view, the variance of the nugget effect is often considered and parameterized relatively to the variance of the GP with λ = σϵ σ 2 ∈ R + . Even for noiseless data as in the case of a deterministic simulator, the nugget effect is often used in GP metamodeling because its practical interest is twofold: both to relax the interpolation property of the GP regression and to improve the conditioning number of the covariance matrix (also referred to as GP regularization). Conceptually, it means that the model function M (numerical simulator in computer experiments) is supposed to be a slightly noised version of a smoother and deterministic simulator. This regularity aspect can also be considered in presence of sparse problems, characterized by a weak density of observations in the input parameter space.

Considerations around the GP trend

Finally, let us say now a brief word about the GP mean (or trend) m(x). A constant m(x) = β 0 or a one-degree polynomial trend m(x) = β 0 + i β i x i is usually considered in practice. But any linear regression model on a set of known basis functions could be used instead. For simplicity, it is assumed, in the rest of the section, that the prior mean is a constant and more exactly equals to zero (assuming that data are centered, for instance). This assumption is only made to simplify some equations (e.g. reminders on maximum likelihood estimation).

Positioning of the GP regression model w.r.t. other metamodels

As mentioned in the introduction, the GP metamodel offers many advantages over other metamodels with regard to the applications considered in our work. On this point, let's recall that we're interested in the emulation of computer experiments from learning bases of a few dozen to a thousand data, with a few dozen explanatory variables as input. Moreover, in the framework of risk and safety assessment applications, it is crucial to have an uncertainty associated with the metamodel predictions in order to take into account the metamodeling error in risk calculations. This is also necessary in adaptive strategies, where the metamodel is used to enrich the learning sample. This is of course what is offered by the GP metamodel which, in addition to its high flexibility and generalization capability, encloses an intrinsic measure of uncertainty.

There are, of course, other metamodels commonly used for emulating computational code output (see, e.g., [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF][START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Smith | Uncertainty quantification[END_REF]). The literature is abundant on all of them and a comparative analysis is beyond the scope of this paper. For a more complete analysis, we can refer to Teixeira et al. [START_REF] Teixeira | Adaptive approaches in metamodel-based reliability analysis: A review[END_REF] and Kianifar and Campean [START_REF] Kianifar | Performance evaluation of metamodelling methods for engineering problems: towards a practitioner guide[END_REF], which also list other previous papers proposing comparative studies. Essentially, we can point out some key elements:

1. Few of metamodels, as the GP does, naturally provide a predictive output distribution, without using bootstrap or bagging techniques which can be costly depending on the type of metamodel to be estimated, and unreliable if the sample size is small; 2. Polynomial regression-based metamodels have the advantage of simplicity but cannot represent complex relationships, while regression based on splines, more flexible, is limited to low dimensions;

3. Radial basis functions-based regression presents obvious connections with GP regression [START_REF] Kanagawa | Gaussian processes and kernel methods: A review on connections and equivalences[END_REF] but does not naturally enclose prediction uncertainty like the GP;

4. This drawback is the same for polynomial chaos expansion metamodels (ans their sparse versions) which have nevertheless many other useful properties for computer experiments applications;

5. Random forests metamodels are particularly robust to the curse of dimensionality, but produce predictors that are often irregular for small learning sample sizes;

6. Support vector regressors can handle non-linearity and high-dimensional inputs thanks to the kernel trick, but require a particularly delicate tuning step of several hyperparameters (kernel parameter, regularization parameter, margin of tolerance);

7. Neural networks are capable of representing very complex phenomena, but at the cost of a large number of data. Even if regularization techniques (e.g. early stopping, regularization, dropout) are used to avoid overfitting, most of application experiences with limited sample sizes (100 to 1000) and several tens of inputs show that the predictivity of a neural network is unsatisfactory.

Furthermore, to obtain a predictive law at the output of a neural network, more complex approaches such as Bayesian neural networks (BNN) have to be considered. BNN has been studied for computer code emulation in Staber and Da Veiga [START_REF] Staber | Benchmarking bayesian neural networks and evaluation metrics for regression tasks[END_REF]. This study concludes that fine-tuning hyperparameters to produce reliable confidence intervals and predictivity remains a challenge, and that only few approximation method for BNN achieve good accuracy (namely stochastic gradient Monte Carlo Markov Chain and deep ensembles). Moreover, the reliability of BNN for higherdimensional problem remains to be studied. Finally, let us mention the promising approach of conformal predictions, that provide prediction intervals with strong theoretical guarantees. They should be developed in the context of computer code emulation (see, e.g., a first step towards this objective in Jaber et al. [START_REF] Jaber | Conformal approach to Gaussian process surrogate evaluation with coverage guarantees[END_REF]).

Estimation of Gaussian process hyperparameters

The type of covariance function is generally fixed among the usual choices of Table 1. It then remains to estimate the covariance parameters (σ 2 , θ) (and eventually λ if a nugget effect is considered). The main estimation procedures are based either on minimization of the squared prediction error calculated by cross-validation (CV), or on maximization of likelihood (denoted MLE for maximum likelihood estimation). Note that Petit et al. [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF] propose a review and comparison of a larger panel of criteria to be optimized to estimate hyperparameters. In particular, they suggest a generalization of the likelihood criterion called the Fasshauer's Hölderized likelihood and which is based on the orthogonal decomposition of the covariance matrix. In parallel to MLE and CV, a Bayesian estimation approach also exists: a prior distribution is assumed for the hyperparameters and combined with MLE to obtain a posterior distribution of the hyperparameters which is then propagated into the GP predictive distribution. The following subsections detail these different approaches and summarize recent works on the theoretical analysis and empirical comparison of these methods.

Cross-validation-based approach

A first estimation approach relies on the mean squared error (MSE) in prediction computed by cross-validation (CV) [START_REF] Hastie | The elements of statistical learning[END_REF]. More precisely, in the case of leave-one-out (LOO) method, the GP hyperparameters are computed by minimizing the LOO-MSE:

LOO-MSE σ 2 , θ := 1 n n i=1 (ŷ -i -y i ) 2 ,
where ŷ-i denotes the GP predictor (mean of predictive distribution) in x (i) when (x (i) , y i ) is removed from the set of observations (this comes down to consider the GP conditioned by Y s,-i ).

Using the CV formulas of Dubrule [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF], the predictive mean ŷ-i and variance ŝ2 -i are given by ŷ-i -

y i = (K σ 2 ,θ y) i (K σ 2 ,θ ) i,i , (7) and ŝ2
-i = 1 (K σ 2 ,θ ) i,i . ( 8
)
The LOO-based estimators of GP hyperparameters are thus given by σ2 M SE , θMSE = arg min

σ 2 ,θ y ⊤ K σ 2 ,θ Diag(K σ 2 ,θ ) -2 K σ 2 ,θ y. ( 9
)
In practice, there are no closed-form expressions for σ2 M SE , θMSE , and the quantity to be optimized is not convex and may have several local optima. The optimization has to be done numerically. Note that under the hypothesis of a well-specified covariance model, some theoretical results exist concerning the consistency of LOO (and more generally CV) estimators [START_REF] Bachoc | Cross-validation estimation of covariance parameters under fixed-domain asymptotics[END_REF].

Maximum likelihood-based approach

The most widely used approach is the MLE which consists in identifying the values of θ which minimizes the negative log-likelihood of the dataset:

ℓ(Y s ) = 1 2 n log(2π) + log K σ 2 ,θ + Y s T K -1 σ 2 ,θ Y s
with |A| denoting the determinant of matrix A. Provided that θ is known and writing K σ 2 ,θ = σ 2 R θ , the MLE estimator of the variance parameter is given by:

σ 2 M LE = 1 n Y s T R -1 θ Y s . ( 10
)
Plugging back σ 2 M LE into ℓ(Y s ) to get a concentrated (or profile) log likelihood involving just θ, MLE results in the following minimization problem for θ:

θ M LE = arg min θ σ 2 M LE R θ 1 n . ( 11
)
Even if calculating the derivative of the above expression is analytic, solving it is not, and no closed form solution can be obtained. As for LOO-MSE approach, numerical methods are thus required to estimate θ M LE . The interested reader can refer to Santner et al. [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] for the MLE equations in the case where a regression model is considered for the trend m(x). In addition, Zhang [START_REF] Zhang | Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics[END_REF] and Bachoc [START_REF] Bachoc | Asymptotic analysis of maximum likelihood estimation of covariance parameters for gaussian processes: An introduction with proofs[END_REF] provide theoretical results on the consistency of MLE (still under the hypothesis of a well-specified covariance model).

As an alternative approach for MLE, Li and Sudjianto [START_REF] Li | Analysis of computer experiments using penalized likelihood in Gaussian kriging models[END_REF] propose to optimize a penalized likelihood function: a penalty term (function of the hyperparameters and depending on a regularization parameter) is added to the standard likelihood. This amounts to maximizing the likelihood function subject to a constraint on the hyperparameters. This can also be interpreted as a Bayesian approach where the penalty term corresponds to a prior on the hyperparameters. Zhang et al. [START_REF] Zhang | A penalized blind likelihood Kriging method for surrogate modeling[END_REF] have recently extended this approach by simultaneously improve the trend and covariance parameter estimation via regularization techniques. However, the problem of such an approach is its cost, as it requires an additional cross-validation loop to estimate the penalization parameters. This is why it does not seem to be very used in practice.

Maximum likelihood estimation in Gaussian process framework may be ill-posed

In their recent papers, 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] discuss how the MLE is wellor conversely ill-posed and how this notion can be defined in the framework of GP hyperparameter estimation, and in the non-asymptotic setting. Note that the scalar notation will be used in what follows for θ without loss of generality for an extension to the case of a vector θ.

Infinite θ and lack of continuity with respect to the training dataset

Under the assumption of a well-specified prior (mean and covariance), the predictive performance of the GP is well understood in an asymptotic setting [START_REF] Bachoc | Cross-validation estimation of covariance parameters under fixed-domain asymptotics[END_REF][START_REF] Karvonen | Maximum likelihood estimation and uncertainty quantification for Gaussian process approximation of deterministic functions[END_REF]. But, there is not much in the literature about the non-asymptotic setting when the hyperparameters θ are estimated from the learning sample, and in a deterministic interpolation framework (interpolating GP to emulate a deterministic function).

First, Karvonen and Oates [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF] defines the MLE ill-posedness as the occurrence of

θ M LE = ∞, which yields k(x, x) = 1 ∀ (x, x) ∈ R × R and estimated correlation matrix R θ M LE = 1 n 1 T n .
The GP predictive distribution is degenerated for each point of prediction (predictive variance is zero and all the probability is assigned to a single value). Having infinite precision from a finite dataset is undesirable and the GP metamodel loses its interest as a tool for uncertainty quantification. It is therefore important to have validation criteria that take into account the whole predictive distribution and thus enable this situation to be detected, such as those presented in section 4.3.

From this definition of ill-posedness, Karvonen and Oates [39, Theorem 2.3] demonstrate that if the data are m-constant, i.e. shifted from the mean function m(x) by a constant c ∈ R:

y i = m x (i) + c for i = 1, . . . , n, with n ≥ 2,
and if the covariance function has a polynomial decaying Fourier transform (like the Matérn functions with smoothness ν > 0), then θ M LE = ∞. On the contrary, if the data are not m-constant, then θ M LE < ∞ so that the predictive distributions are non-degenerate.

Then, the authors propose to consider the classical definition of well-posedness of Hadamard for an inference or estimation problem. More precisely, it is well-posed if (i) a solution exists, (ii) the solution is unique, and (iii) the solution depends continuously on the data. If these conditions are not met, the problem is ill-posed. The (iii) condition about the sensitivity of the GP estimate to the training dataset is particularly relevant for sensitive applications such as nuclear safety applications. From this definition, and under the same assumptions of regularity of the covariance function (polynomial decaying Fourier transform), Karvonen and Oates [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF] demonstrate that if θ = ∞, the GP regression problem is ill-posed by violation of (iii), in the sense that the resulting predictive distributions are not locally Lipschitz in the data w.r.t. the Hellinger distance, which means that predictive inference can be sensitive to small perturbations of the dataset.

In order to find solutions to this problem, the authors examine several alternative modelling and estimation methods and show that:

• the LOO-CV estimator of θ shares the same undesirable property of MLE when the data are m-constant, which is not surprising given the close connection between MLE and CV (Fong and Holmes [START_REF] Fong | On the marginal likelihood and cross-validation[END_REF]);

• the addition of a parametric prior mean function m(x), also estimated from the data (i.e., as in universal kriging), does not help;

• simultaneous MLE of σ and θ does not prevent ill-posedness (as usually done in GP regression, see Eqs. [START_REF] Marrel | The ICSCREAM methodology: Identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics[END_REF][START_REF] Chilès | Geostatistics: Modeling spatial uncertainty[END_REF]);

• the addition of a nugget effect does not prevent from having θ M LE = ∞ when the data are m-constant but prevents the GP predictive distribution from being degenerate. Moreover, when the data are not m-constant but close to it, θ M LE → ∞ and the condition number of the correlation matrix R θ M LE increases with a rate related to the smoothness of the covariance. This results in a numerical issue in the MLE process (likelihood for large value of θ cannot be computed). The introduction of a λ > 0 allows to mitigate this issue (by upper bounding the condition number). The price to pay is therefore to relax the interpolation constraint.

In conclusion, the practical recommendations to be retained from Karvonen and Oates [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF] are:

▶ the addition of a nugget effect λ > 0, also estimated by MLE, is recommended as a regularization parameter, taking care to limit its value by a reasonable upper bound. Another solution can be to bound the parameter θ and find the MLE estimate in (0, θ max ] (constrained MLE). However, it implies the arbitrary choice of a θ max which turn out to be the estimated value of θ M LE if the data are m-constant.

▶ Finally, reasoning in a rather rudimentary way, this would argue for not assuming a too complex model for the GP trend m(x) so that the observed data are not m(x)-constant. In practice, m(x) is assumed to be a constant or a one-degree polynomial (universal kriging), with parameters also estimated by MLE. Another practical recommendation could be to check if the data are m-constant (or close to be) in order to detect problematic datasets which cause GP regression to be ill-posed. However, it is reasonable to expect that in our industrial applications, the probability of having strictly m(x)-constant data (without any deviation from m(x) or noise) is relatively low due to the complexity of the considered models (or codes), the large input dimension and the sparsity of the data (low sample size).

Lack of robustness

Previously to Karvonen and Oates [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF], Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] had already proposed to define the illposed MLE problem, which the authors prefer to refer to as a lack of robustness. Their starting point is that the likelihood (Eq. ( 11)) is sometimes very flat in the tails. This is clearly illustrated on our analytical example of Figure 4 -case (1), as in Karvonen and Oates [START_REF] Karvonen | Maximum likelihood estimation in Gaussian process regression is ill-posed[END_REF]: θ = ∞ and consequently R θ = 1 n 1 T n . In multidimensional case (i.e. d > 1), this case corresponds to

θ i = ∞ for all 1 ≤ i ≤ d; -case (2): θ = 0 and consequently R θ = I n . For d > 1, this is encountered when ∃i, 1 ≤ i ≤ d, for which θ i = 0.
Case (2) is not more desirable than case (1) because it means that R θ M LE is near I n and the GP predictor is an impulse function interpolating the observations, while following the GP mean m(x) elsewhere. The authors state and numerically show that even if "such degeneracies are somewhat unusual in one-dimension, they are not particularly unusual with higher dimensional inputs" (large dimension d).

The authors also mention that MLE instability can often be "overcome by adding a nugget effect, but studies have found that the features of the emulator can significantly change when a nugget is added [START_REF] Berger | Objective Bayesian analysis of spatially correlated data[END_REF]." So, to circumvent the problem of robustness, Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] prefer to focus on Bayesian approaches and demonstrate that certain prior and parameterizations for the GP parameters result in a more robust estimation than others (see Section 5.2 dedicated to the so-called RobustGaSP Bayesian approach).

Bayesian approach

A third approach is to consider a full-Bayesian approach where a prior is assumed on the GP hyperparameters. The marginal posterior distribution is then inferred by Bayes' rule from marginal likelihood of data, and with regard to the prior. The resulting posterior uncertainty is then integrated in the GP predictive distribution. More precisely, assuming a prior on the hyperparameters

(σ 2 , θ) ∼ π (σ 2 , θ), their posterior distribution writes π (σ 2 , θ | Y s ) ∝ π (Y s | σ 2 , θ) π (σ 2 , θ).
The probability density function of the GP predictive distribution of Y (x) is then given by:

p (y(x)|Y s ) = p y(x) | Y s , σ 2 , θ π σ 2 , θ | Y s dσ 2 dθ. ( 12
)
Full Bayesian approach thus allows to take into account the uncertainty on the estimation of the GP hyperparameters and to propagate it in the GP predictive law. It has been illustrated in dimension two by Wieskotten et al. [START_REF] Wieskotten | A comparison between Bayesian and ordinary kriging based on validation criteria: application to radiological characterisation[END_REF] showing that it is relevant and can outperform the ordinary GP in terms of both predictivity and accuracy of predictive intervals. This is especially true when the sample size is small, the benefit decreasing as the size increases, as one might expect.

However, the tractability of the full Bayesian approach in higher dimension remains a major obstacle to its use. Indeed, in practice, the computation of π (σ 2 , θ | Y s ) and p (y(x)|Y s ) requires the use of Markov chain Monte Carlo (MCMC) methods like Metropolis-Hastings algorithm [START_REF] Robert | Monte Carlo statistical methods[END_REF] or Hamiltonian Methods [START_REF] Neal | Priors for infinite networks[END_REF]. The calculation cost of the predictive distribution with MCMC techniques becomes expensive in large dimension. To circumvent this limitation, some plug-in approaches can be considered as in Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF]: p (y(x)|Y s ) is computed with the GP hyperparameters fixed at the maximum a posteriori probability (MAP) estimate (that equals the mode of the posterior distribution π (σ 2 , θ | Y s )). In this case, the Bayesian framework is only used to compute this posterior distribution and is then discarded to calculate the predictive law p (y(x)|Y s ). Basically, this means replacing π (σ 2 , θ | Y s ) with a Dirac distribution whose mass is concentrated on the MAP value. This somewhat brutal solution facilitates the intensive use of the predictive distribution. Only the problem of the estimation of π (σ 2 , θ | Y s ) remains. Finally, the Bayesian approach (whether full or plug-in) also requires to define the prior distribution for the hyperparameters, choice which can be of prime importance as detailed in the work of Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] (See Section 5.2).

Discussion on the relative practical performance of the different approaches

Discussions on the choice of a method between MLE and CV (or LOO) methods are not new, but recent work and especially intensive benchmarks are shedding new light. Let us try to make a brief synthesis.

As shown by Bachoc [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF], the MLE method is optimal when the covariance function is wellspecified, i.e. when the "true" covariance function belongs to the assumed parametric set of covariance functions. If this is not the case (misspecification case), there is no more guarantee that the MLE method would perform correctly and optimally. Bachoc [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF] illustrated that MLE may not be very robust in this case, especially if the number of data is small, while the CV-based approach performs better. On the other side, the availability of gradients in the MLE case (without significant additional computational cost) is an advantage in the implementation of numerical optimization algorithms required for hyperparameter estimation. However, this advantage must be qualified by the recent works on the computational complexity of cross-validation schemes and more precisely the fast computation of gradients of LOO criteria [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF][START_REF] Petit | Improved Gaussian process modeling: Application to Bayesian optimization[END_REF]. Still in Petit [START_REF] Petit | Improved Gaussian process modeling: Application to Bayesian optimization[END_REF], an intensive benchmark on analytic functions of different dimensions shows that MLE is often preferable to its competitors (not only in well specified cases but also in case of overestimated regularity), and that the choice of regularity (ν in Matérn class) might be often more important than the estimation of GP hyperparameters. To conclude the comparison between MLE and LOO, let us mention the remark of Zhang and Wang [START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF] concerning the flatness of both MLE or LOO-based criteria around their optimal value. The authors argue that the flatness of LOO-based criteria is less damaging since it indicates that the predictive distribution is less sensitive to the hyperparameter value in the flatness region.

Faced with this lack of consensus between MLE and LOO, the full-Bayesian approach could appear as a relevant solution as it may yield more robust predictions (see, e.g., [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF]). But, this approach strongly depends on the prior distribution of the hyperparameters, as highlighted by Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF], and has a much higher computational cost, especially if the number of hyperparameters is high. Besides, to the best of our knowledge, there are no applications using this approach to emulate numerical simulators in large dimension (d ≥ 10 for instance). Only the RobustGaSP Bayesian method proposed by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] and detailed in Section 5.2 with its specific priors and approximations, could overcome these limitations.

However, another solution would also be relevant: both the MLE and CV criteria would be integrated into the estimation of the hyperparameters. A multi-objective procedure could be developed where the MLE criterion would remain the main reference objective and another LOObased criterion could be considered as a complementary criterion. This is the purpose of the work proposed in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A new estimation algorithm for more robust prediction[END_REF] of the present article. More generally, we also think that a nugget effect has to be considered and estimated jointly with the other hyperparameters as it allows to enrich the family of covariance functions (this having not been considered in the benchmark of Petit et al. [START_REF] Petit | Parameter selection in Gaussian process interpolation: An empirical study of selection criteria[END_REF]). It also facilitates the MLE by regularizing the likelihood function, improving the conditioning of the correlation matrix and numerical convergence of algorithms. However, this nugget effect can be double-edged into a Bayesian approach as it may increase identifiability problems. Its use should be restricted to simple MLE.

As a synthesis of this section, including more recent methods presented in Section 5, Table 2 provides a summary of the key elements relating to GP estimation methods, to provide a clearer synthetic overview encompassing advantages and practical refinements, disadvantages and limitations, and other relevant details.

Quantitative criteria for Gaussian process validation

Once the GP metamodel has been estimated, its predictive capabilities need to be checked to ensure confidence in its use (as a substitute for the simulator). Thus, validation criteria must be defined to assess the accuracy of the GP predictor, but also of its prediction variance, its covariance and more generally of the whole GP conditional distribution. Checking the reliability of the entire predictive law thus enables simultaneous validation of both the chosen covariance model and the estimated hyperparameters. These criteria can also be used in practice to choose between several covariance models. For example, Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF] have illustrated the significant impact of an unsuitable choice of covariance. Hence, different quantitative criteria have been proposed (see, e.g., [START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF][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]) and are listed in the following. A new criterion, the IAEα criterion, is also proposed. Note that we are focusing on the validation criteria that can be applied to a real test case, in the sense that they can be estimated on a single learning sample using cross-validation or LOO techniques. They are therefore formulated in the following in their LOO version. Of course, they can be defined (and computed) in a similar way on a test sample, different and independent from the learning sample. Similar expressions can also be obtained with K-fold cross-validation [START_REF] Hastie | The elements of statistical learning[END_REF].

Criteria to assess the accuracy of the Gaussian process predictor ŷ (x)

Classically, the root mean squared error (RMSE) writes

RMSE = 1 n n i=1 y x (i) -ŷ-i x (i) 2 1/2 (13)
and its counterpart expressed in terms of the proportion of the variance explained, namely the predictivity coefficient Q 2 (see, e.g., [START_REF] Fekhari | Model predictivity assessment: incremental test-set selection and accuracy evaluation[END_REF]):

Q 2 = 1 - RM SE 2 1 n n i=1 y i -1 n n i=1 y i 2 . ( 14
)
The closer to one the Q 2 , the better the accuracy of the metamodel predictor. On the contrary, a zero Q 2 indicates very poor predictive abilities, i.e. equivalent to what would be obtained with the empirical mean of the observations. Note that both RMSE and Q 2 correspond to averaged indicators and should be complemented with a plot of observed data versus predicted values for more detailed analysis. Find the posterior distribution of (𝝈 𝟐 , 𝜽) inferred by Bayes' rule from marginal likelihood of data and prior distribution of (𝝈 𝟐 , 𝜽) Use Maximum a posteriori (MAP) estimate Robustness increased and flatness reduced thanks to the prior distribution. Some theoretical guarantees when using reference prior and for specific covariance model In practice, outperforms MLE when 𝑛 is small (𝑛 < 50) Quantification of the residual uncertainty on (𝜎 2 , 𝜽) estimates (information given by the posterior distribution of (𝜎 2 , 𝜽) itself) Plug-in approach with MAP estimates of (𝜎 2 , 𝜽) can be used to simplify the computation of a GP predictive distribution No closed-form expression for the posterior distribution, MCMC methods required Very high computation cost of MCMC (in addition to the calculation cost of each MLE evaluation). Untractable when the number of hyperparameters 𝑑 is large (𝑑 ≥ 10) Problems of mixing the Markov chain (results in correlated marginal distributions) when 𝑑 is large Require a relevant/expert choice of prior distribution RobustGaSP Bayesian method Section 5.2.2 Bayesian approach with specific prior distribution of (𝝈 𝟐 , 𝜽) and approximations Efficient and tractable variant of Bayesian approach, based on robust prior, plug-in approximation with MAP estimates of (𝜎 2 , 𝜽) Robustness demonstrated for some (usual) covariance models Computationally tractable even when 𝑑 is large with the approximation of robust prior proposed by the author Available R code library with C++ blocks for higher performance Complexity of the method, difficult to implement from scratch High computation cost of reference prior and MAP estimates of (𝜎 2 , 𝜽) Still few applications on industrial cases Choice of certain parameters a priori (for reference prior in particular). Dependence of results on these choices?

Robust Prediction Intervals

Estimation (RPIE) Section 5.1

Correct initial (𝝈 𝟐 , 𝜽) estimates to ensure a good LOO-coverage probability of a given level 𝜶

The only existing method proposing a correction of (𝜎 2 , 𝜃) to directly control the reliability of prediction intervals Improvement based on prior estimation of another MLE or CV method (small risk of performance degradation w.r.t. them) Strong dependence to the initial values of (𝜎 2 , 𝜽) estimates (obtained from MLE or CV method, e.g.), ratio of components of 𝜽 preserved Procedure to be performed for each given level 𝛼: a different set of corrected values of (𝜎 2 , 𝜽) for each 𝛼, no single corrected GP metamodel Numerical benchmark provided by the author is too poor to conclude that the method is efficient Table 2: Synthesis of main estimation methods and recent corrective approaches.

Criterion to evaluate if the conditional Gaussian process variance is of the right order of magnitude

Other important indicators propose to deal with the model variance [START_REF] Bastos | Diagnostics for Gaussian process emulators[END_REF]. We focus here on the predictive variance adequacy (PVA) factor (see, e.g., [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF][START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF]):

PVA = log 1 n n i=1 (y i -ŷ-i ) 2 ŝ2 -i . ( 15
)
In order to get reliable prediction intervals from the model, the prediction variances should be of the same order of the prediction errors so that the PVA should be close to zero. In summary, the smaller the PVA, the more reliable the prediction intervals. On the contrary, too low prediction variances w.r.t. the prediction errors (i.e. an "overconfident" predictive model) or too large prediction variances ("underconfident" or too uncertain predictive model) yield poor PVA. A more detailed analysis and interpretation of its values is available in Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF].

Criteria to assess the accuracy of the whole Gaussian process predictive distribution

The logarithmic score [START_REF] Bernardo | Expected Information as Expected Utility[END_REF] is defined as the negative logarithm of the predictive density evaluated on the observations:

LogS = n 2 log(2π) + n i=1 log ŝ-i + 1 2 (y i -ŷ-i ) 2 ŝ2 -i . ( 16
)
Some similarities exist between LogS and PVA, both depending on the standardized residuals. But Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF] point out that PVA is preferable for GP validation since it will mitigate the effect of extreme values and it will similarly penalize models with too large or too small predictive variances. In contrast, the weighting of the two terms in LogS will less penalize the too large predictive variances. Then, we have all the class of criterion based on the reliability of predictive intervals. As recalled by Zhang and Wang [START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF], from the Brier score defined for the predictive cumulative distribution

F -i by BS(y) = 1 n n i=1 F -i (y) -1 Y x (i) ≤ y 2 ,
with 1{A} the indicator function of A, the continuous ranked probability score (CRPS, [START_REF] Gneiting | Probabilistic forecasts, calibration and sharpness[END_REF]) is defined as the integration of BS:

CRPS = 1 n n i=1 ∞ -∞ F -i (y) -1 Y x (i) ≤ y 2 dy = ∞ -∞ BS(y)dy. ( 17
)
For a GP metamodel, it can be demonstrated that:

CRPS = 1 n n i=1 ŝ-i (y i -ŷ-i ) ŝ-i 2Φ (y i -ŷ-i ) ŝ-i -1 + 2ϕ (y i -ŷ-i ) ŝ-i - 1 √ π ( 18
)
with ϕ and Φ respectively denoting the probability and cumulative functions of the standard Gaussian distribution. CRPS is more robust than LogS, but it will tend to favour models with small predictive variance, subject to similar calibration performance. This is not desirable for GP validation since the objective is to have the most reliable predictive distribution, more than the most concentrated, especially in a safety study framework. Moreover, as illustrated by Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF] on their application, the CRPS-based criterion does not allow to identify an inaccurate covariance model when this mismodeling only affects the predictive variance (and not the predictor), unlike the PVA or other criteria based on the predictive interval and presented immediately afterwards.

By focusing on the validation of GP prediction intervals (PI), the level α ∈]0, 1[ of any PI can be compared to the proportion of observations that actually lie within this interval. This proportion also called empirical coverage function [START_REF] Schruben | A coverage function for interval estimators of simulation response[END_REF] is defined as:

∆(α) = 1 n n i=1 1 y i ∈ P I α,-i x (i) , ( 19
)
where P I α,-i is the α-level prediction interval for the point x (i) built from the Gaussian distribution N (ŷ -i , ŝ2 -i ). ∆(α) corresponds to the empirical LOO estimator of what is called the prediction interval coverage probability in the deep learning literature.

From this, a graphical tool referred to as α-PI plot can be built by plotting ∆(α) against α (see, e.g., [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF][START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF]). By definition, the more the points should be located around the y = x line, the more reliable the GP predictive intervals are. In order to have a quantitative indicator summarizing the quality of the α-PI plot, we naturally propose to consider the following IAEα criterion:

IAEα = 1 0 ∆(α) -α dα. ( 20
)
This criterion denoted IAEα for integrated absolute error on α corresponds to the area between the alpha plot and the reference line. IAEα lies in [0, 1] and the closer to zero the IAEα, better the PI in average. Note that this criterion is very close to the so-called "mean squared error α" of Wieskotten et al. [START_REF] Wieskotten | A comparison between Bayesian and ordinary kriging based on validation criteria: application to radiological characterisation[END_REF], defined with L 2 norm instead of L 1 norm. As in Wieskotten et al. [START_REF] Wieskotten | A comparison between Bayesian and ordinary kriging based on validation criteria: application to radiological characterisation[END_REF],

IAEα will be computed in practice with a regular discretization of α over ]0, 1[. The L 1 norm is preferred here to give a homogeneous weight whatever α, avoiding to give too much weight to the strongest deviations and having a direct interpretation w.r.t. the α-PI plot. More generally, validation metrics based on coverage probability of the predictive distribution, such as ∆(α) and IAEα, go beyond simple validation of the predictive variance. They are therefore better adapted to assessing the validity of confidence intervals, regardless of the nature of the predictive laws, whether Gaussian or non-Gaussian.

Remark 1. Practical recommendation.

The last group of criteria composed of ∆(α) and those derived from it, namely the α-PI plot and IAEα, are perfectly adapted to the control of the reliability of prediction intervals. However, they should not be used alone but in addition to a prior control of predictivity with Q 2 (or RMSE). To better understand this recommendation, let us take for example the extreme case with the following metamodel: a constant predictor corresponding to the empirical mean of the data, a constant prediction variance, equal to the empirical variance of the data, and a Gaussian predictive distribution. Finally, let us assume that the sample of observed data follows a distribution with a Gaussian shape. We would then obtain for the predictive law a Q 2 equal to 0 but a ∆(α) very close to the right level for each α. Consequently, ∆(α) is only of interest and should only be used if and only if the predictivity of the metamodel has already been checked and controlled (via RMSE or Q 2 for instance).

As a synthesis of this section, Table A.3 in Appendix A, complementary of the one previously proposed by Demay et al. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF], gives a synthetic overview of the various validation criteria listed above.

Illustration of validation criteria

To illustrate the interest of considering several criteria in the validation process, let us consider the example of the "re-scaled" Branin function [START_REF] Picheny | A benchmark of kriging-based infill criteria for noisy optimization[END_REF]. It is defined for two independent uniform inputs X 1 and X 2 on [0, 1] by:

M Branin (X 1 , X 2 ) = 1 51.95 a(V 2 -bV 2 1 + cV 1 -r) 2 + s(1 -t) cos(V 1 ) + s -44.81 , ( 21
) with V 1 = 15X 1 -5, V 2 = 15X 2 , a = 1, b = 5.1 (4π 2 )
, c = 5 π , r = 6, s = 10 and t = 1 8π . This function is illustrated in Figure 2. Consider a GP metamodel with a 2-D tensorized Gaussian covariance and conditioned by a fixed random learning sample of n = 30 points (red points on Figure 2). A random test sample of 1000 points is used to compute the validation criteria. Considering three different sets of hyperparameters [θ 1 , θ 2 ], Figure 3 illustrates that they can lead to metamodels with a very similar Q 2 (Q 2 ∼ 0.9) but significantly different IAEα (namely 0.2 and 0.05), and, on the contrary, a similar IAEα ( IAEα ∼ 0.05) and very different Q 2 (namely 0.9 and 0.62). Moreover, considering the two first sets of hyperparameters, it also illustrates the resulting major difference in the level of the GP prediction intervals compared to the desired level. An increase of 0.15 in IAEα results here in much more conservative and less realistic intervals. For instance, prediction intervals of level α = 0.7 actually include on average nearly 95% of the data. Evaluating only Q 2 (or RMSE), as is often done, does not make it possible to dissociate the two metamodels, whereas one of the two offers a much better assessment of prediction uncertainty.

This simple example illustrates that the criteria should be used in a complementary way. It is also a good opportunity to illustrate the difficulty of finding an optimum for MLE estimation, as shown by Figure 4, in the middle column. It clearly shows the flatness of the objective function derived from the MLE (Eq. ( 11)) in the most interesting regions (in blue). Remember that the aim is to minimize this function. A red square indicates the global minimum of the function on [0, 2] 2 and a green square shows the MLE estimate found by a standard optimization algorithm, probably stuck in a local minimum. A logarithmic scale (base 10) is also proposed in the bottom middle graph to highlight the level lines in the low-value areas. However, not having found the global optimum of the log-likelihood is not detrimental to the GP model's performance here, since its predictivity is almost perfect (Q 2 ∼ 0.99), as well as its confidence intervals (see bottom and top left graphs). Note that the results obtained with another estimation method, the Bayesian RobustGaSP approach, are also shown in Figure 4: they are discussed in Section 5.2.2 dedicated to this approach. By way of illustration, over the range of variation [0, 2] 2 , 22% (resp. 15%) of θ values lead to a poor IAEα ≥ 0.2 with a Q 2 ≥ 0.8 (resp. Q 2 ≥ 0.9). At least a Q 2 ≥ 0.95 is required to reduce to 5% this number of poor IAEα cases.

More generally, for a more in-depth study of the links and connections between estimation and validation criteria, we refer the interested reader to our companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A new estimation algorithm for more robust prediction[END_REF] in which other analytical examples are proposed, particularly in higher dimensions (from 3 to dimension 20). Furthermore, for various testcases, the boxplots of Q 2 , PVA and IAEα values are plotted according to the size n of the training sample, for various covariance functions and for different hyperparameter optimization methods. This notably allows to appreciate the range of variation of IAEα w.r.t. the one of Q 2 , for values found as optimal for the MLE (values found by the optimization algorithms).

Recent developments for more reliable Gaussian process predictions

In the case of emulation of deterministic functions from a numerical simulator, there is no guarantee that the function to be emulated is part of the set of possible trajectories generated by the assumed GP model (case of model misspecifications). As a result, there is no guarantee either that the MLE (or any other estimation approach) will work, or that the estimated prediction variance will control the metamodel error. More generally, the GP predictive distribution may not accurately cover unobserved data. Based on this statement, Acharki et al. [START_REF] Acharki | Robust prediction interval estimation for Gaussian processes by cross-validation method[END_REF] were interested in correcting the GP estimated hyperparameters to adjust the GP prediction intervals and ensure better coverage probabilities of the GP predictive distribution. In the following, we outline this approach, which is further detailed in Appendix B, before explaining the limits of the method in our application context.

A corrective approach to directly modify the bounds of Gaussian process prediction intervals

First of all, Acharki et al. [START_REF] Acharki | Robust prediction interval estimation for Gaussian processes by cross-validation method[END_REF] assume that he has a first set of estimated GP hyperparameters obtained by MLE or CV method and denoted ( σ 2 0 , θ 0 ). Then, for a given level α of prediction intervals, they consider what they refer to as "Leave-One-Out Coverage Probability" and which corresponds to the ∆(α) given by Eq. [START_REF] Demay | Model selection for Gaussian process regression: an application with highlights on the model variance validation[END_REF]. ∆(α) can be rewritten as a function of the quantiles of the GP predictive distribution, quantiles of level (1 -α)/2 and (1 + α)/2. So, to ensure a good value of ∆(α) (i.e. close to α), the authors propose to find two new sets of hyperparameters, respectively denoted σ2 , θ and (σ 2 , θ) which guarantee that the GP quantiles respectively of level (1 -α)/2 and (1 + α)/2 will yield ∆(α) = α. In addition to satisfying this condition, these two sets of hyperparameters are found to be as close as possible to the initial values ( σ 2 0 , θ 0 ), in the sense of a similarity measure (see Appendix B for more details). Hence, the initial set ( σ 2 0 , θ 0 ) is only used to build the GP predictor. σ2 , θ and (σ 2 , θ) (which can be viewed as corrections of the initial set) then yields two other GP metamodels which are used to compute the two bounds (i.e. the lower and upper GP quantiles) corresponding to the α-prediction interval. The resulting method is called Robust Prediction Intervals Estimation (RPIE) by the authors.

The numerical tests proposed by the authors show that when the GP metamodel is well specified (good coverage probability of intervals obtained with initial MLE or CV-estimates θ 0 ), the RPIE method does not bring any added value. In the opposite case, the RPIE method is relevant and corrects efficiently the prediction intervals.

To the best of our knowledge, the RPIE method is the first to focus on a correction of the hyperparameters to control the quality of prediction intervals but, it seems perfectible on several points: ▶ First of all, the new sets of hyperparameter strongly depend on the set of initial values (estimated by MLE or CV), since the searched solutions are expressed as an isotropic shift of θ 0 . In multidimensional case, the ratio between the different θ 0,i i=1...d will thus be preserved, even if it was initially badly estimated. The method thus assumes that the initial estimation has been carried out correctly, under the assumption of a well-specified covariance.

▶ Secondly, the procedure proposed must be performed for the two bounds of any desired interval of level α, yielding two sets of corrected hyperparameters. RPIE does not look for a single hyperparameter correction that would provide a single GP that would both adjust the data well and provide reliable prediction intervals. Not having a single GP metamodel is not satisfactory in the perspective of using the GP metamodel as a multi-objective tool, i.e. to predict a quantile or a probability [START_REF] Marrel | The ICSCREAM methodology: Identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics[END_REF], estimate sensitivity indices [START_REF] Gratiet | A Bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF], perform an optimization [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences[END_REF], or more generally to implement a SUR (Stepwise Uncertainty Reduction) approach [START_REF] Chevalier | Fast krigingbased stepwise uncertainty reduction with application to the identification of an excursion set[END_REF]. These studies cannot be implemented directly from Acharki et al. [START_REF] Acharki | Robust prediction interval estimation for Gaussian processes by cross-validation method[END_REF]'s approach. We argue that it is simpler and preferable to have a single GP and to address the misspecification problem apart from the hyperparameter estimation, by testing and comparing different covariances.

▶ More fundamentally, it is not a correction of the GP predictive law but only of its quantiles (high and low) for a given level of prediction interval α It seems more interesting to us to correct the predictive law to simultaneously control ∆(α) whatever α.

▶ Finally, the numerical benchmark carried out by the authors is not exhaustive enough. Only 3 numerical examples are considered, all in dimension d = 10 and with a relatively large sample size n = 600. Moreover, only one Monte Carlo sample is drawn. More thorough tests appear necessary to evaluate the robustness of RPIE method to a poor estimation of the initial hyperparameters, to sampling variability, and to a smaller sample size (especially for the correction of high or very low quantiles).

A more efficient Bayesian approach

As mentioned in Section 3.4, the Bayesian approach used to estimate the GP hyperparameters relies on two key ingredients: first, the choice of a prior distribution π (σ 2 , θ) and second, the estimation of the posterior distribution and its propagation in the GP metamodel (see Eq.( 12)).

Discussion on the choice of prior and focus on reference and Jeffreys priors

As highlighted by Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF], prior knowledge about the GP hyperparameters is often lacking. It then seems natural to use non-informative priors but these may fail to lead to a proper posterior (i.e., a distribution that integrates to a finite mass), this condition being necessary in our context of quantifying the uncertainty of GP parameters. Note that, for the 1-D or isotropic cases (i.e. θ ∈ R), proving that the reference posterior is proper amounts to finding appropriate upper bounds on the tail rates of π (θ | Y s ) π(θ) as θ → 0 and as θ → +∞ (see Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF] for details).

With regard to this issue and still for the isotropic case, Berger et al. [START_REF] Berger | Objective Bayesian analysis of spatially correlated data[END_REF] first showed that among several prior distributions the reference prior of Bernardo [START_REF] Bernardo | Reference analysis[END_REF] is the most satisfying default choice. Reference priors [START_REF] Berger | The formal definition of reference priors[END_REF] aim to formalize the notion of "uninformative prior" and are defined so as to maximize a measure of distance or divergence between the posterior and prior, as data are observed (this choice allows the data to have maximum effect on the posterior estimates). Berger et al. [START_REF] Berger | Objective Bayesian analysis of spatially correlated data[END_REF] demonstrated that reference priors yield proper posterior for isotropic rough correlations that include the exponential correlation and the set of the Matérn family with smoothness parameter ν ⩾ 1. Their demonstration notably relies on the fact that the correlation kernel cannot then be twice continuously differentiable at 0. Very recently, Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF] succeeds in extending this to a large class of smooth kernels, which includes the Gaussian and the Matérn family with smoothness parameter ν < 1 [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF]Theorem 4.4]. But this extension remains in the only case where the dimension of θ is equal to one (isotropic covariance). In a nutshell, the key behind Berger et al. [START_REF] Berger | Objective Bayesian analysis of spatially correlated data[END_REF] and Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF]'s demonstrations is that the reference prior should "compensate for the marginal likelihood" so that the integrated likelihood (i.e. posterior distribution) has the right decay rates on the distribution tails.

The extension of previous results to the anisotropic case (θ ∈ R d with d > 1) is obtained by defining anisotropic correlations as products of one-dimensional rough correlations and a possible additional nugget effect [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF]. The demonstration relies on the use of Jeffreys prior (obtained as the square root of the determinant of the Fisher information matrix) which is here a reference prior since the authors consider a separable product correlation function. Note therefore that this extension is not valid for anisotropic geometric covariances. Moreover, the one-dimensional correlation functions considered to build the tensorized covariance are again assumed to be rough. And, unfortunately, the proof used in Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF] to deal with smoother kernels cannot easily be adjusted to the much more complex prior considered by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF].

In summary, and to the best of our knowledge, there are only two configurations in which obtaining a proper posterior law π (θ | Y s ) from a reference prior π(θ) could be established:

▶ isotropic covariance (θ ∈ R) with the works of Berger et al. [START_REF] Berger | Objective Bayesian analysis of spatially correlated data[END_REF] for rough covariances (e.g. exponential, spherical or Matérn with ν < 1) and Muré [START_REF] Muré | Propriety of the reference posterior distribution in Gaussian process modeling[END_REF] for smooth covariances (e.g. Gaussian or Matérn with ν ⩾ 1), ▶ anisotropic covariance defined from a product of one-dimensional correlation functions [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF].

Focus on RobustGaSP method

Despite the fact that the demonstration of obtaining a proper posterior is not theoretically established for product of smooth covariances, the most relevant works in this context seem to be those of Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] which is referred to as the RobustGaSP method. So, to build their approach, the authors rely on the three following ingredients.

• A robust prior π R (β, σ 2 , θ) is assigned to the set of GP parameters (details in Appendix C). From this, a marginal likelihood can be deduced as a function of θ alone. The marginal posterior distribution of θ is then obtained.

• Sampling this marginal posterior distribution calls for the use of a Metropolis algorithm (MCMC sampling). Because of the cost of each likelihood evaluation (in O (n 3 )) and the associated computational error which can be very large especially when the correlation matrix is close to the matrix 1 n 1 T n , Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] do not advocate this method. Instead, they recommends to estimate (and consider) only the MAP (i.e. the mode) of the posterior distribution π (θ | Y s ). The Bayesian approach then stops at this step: the MAP is directly used to compute the GP predictive distribution in a so-called "plug-in" approach. This avoids a too much computationally prohibitive MCMC algorithm to estimate the posterior distribution of Eq. ( 12). In a nutshell, the reference prior distribution acts only as a penalizing factor on the likelihood and the Bayesian framework only ensures a robust estimation of the hyperparameters, hence the name Robust GaSP method (GaSP is the acronym for Gaussian Stochastic Process). Note that an alternative solution (not detailed here) is proposed by Muré [START_REF] Muré | Optimal compromise between incompatible conditional probability distributions, with application to objective Bayesian kriging[END_REF] to make the full-Bayesian procedure tractable: univariate conditional Jeffreys-rule posterior distributions and pseudo-Gibbs sampler are notably used.

• Finally, Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] consider some specific choices of reparameterization. Even if MLE is invariant under a reparameterization (injective transformation of GP hyperparameters), this is not the case of the MAP of posterior distribution because of the presence of the Jacobian for the prior. As a result, the authors consider other common ways of parameterizing the hyperparameters θ, in particular the inverse and log-inverse transformations.

From these elements, Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] establish two theorems related to the robustness of the estimated MAP, one theorem with nugget effect and another without. More precisely, under the assumption of the reference prior defined by the equation (C.2) with a = 1 and using the standard parametrization or the log-inverse reparametrisation, the authors demonstrate that the estimation of the MAP of the posterior distribution (Eq. (C.3)) is robust for the tensorized form of Matern, spherical and exponential correlation functions [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF]Theorem 3.1]. Recall that the robustness defined by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] refers to the two extreme cases described in Section 3.3 and leading to R θ = 1 n 1 T n and R θ = I n . However, a serious drawback remains in RobustGaSP method: the computational cost required to compute the reference prior and, even more, to compute the mode of the posterior distribution π (θ | Y s ). Even if in fine the MAP estimate is used instead of the posterior sampling, the use of mode search algorithms, such as the quasi-Newton optimization method, typically relies on the information of the derivatives. Computing the derivative requires more evaluations of the likelihood and therefore more inversions of the correlation matrix. The total cost of the procedure then becomes prohibitive. To overcome this limitation, Gu [START_REF] Gu | Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection[END_REF] introduces an approximation of the reference prior π R (β 0 , σ 2 , θ), which he calls the jointly robust prior. The author demonstrates that this prior retains the robustness property while being computationally simpler than the reference prior for the purpose of hyperparameter estimation. As an illustration, the RobustGaSP method is applied to the analytical example of the Branin function proposed in Section 4.4, with the approximated reference prior and the default parameter values of Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF]'s algorithm (package from [START_REF] Gu | RobustGaSP: Robust Gaussian Stochastic Process Emulation[END_REF]). The estimated hyperparameters are indicated on the middle plots of Figure 4 (by a magenta square) and the performances of the GP metamodel obtained with these hyperparameters are illustrated by the plots on the right. Note that we observe the poor quality of the prediction intervals (bottom right plot). Indeed, even in this case of a quasi-perfect PG predictor (Q 2 ∼ 0.99), this issue can arise due to an improper estimation of σ (directly connected to a poor estimation of θ).

The approach proposed by Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] combined with the approximation of the reference prior π R (β 0 ) of Gu [START_REF] Gu | Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection[END_REF] appears to be one of the existing approaches most in line with the idea of having a more robust estimation of hyperparameters (rather than adopting an intractable full-Bayesian approach)(and indeed tractable compared to a full Bayesian approach). Moreover, the proposed approximation of the reference prior allows to consider the application of the method even in the case of a large number of input variables. Finally, its availability in a dedicated software package [START_REF] Gu | RobustGaSP: Robust Gaussian Stochastic Process Emulation[END_REF] makes it the most interesting existing method. For all these reasons, RobustGaSP method will be compared with the new algorithm proposed in the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A new estimation algorithm for more robust prediction[END_REF].

The alternative corrective approaches detailed in this section are included in Table 2, which enables them to be clearly positioned in relation to the standard approaches.

Conclusion

The value of GP regression for emulating costly computational codes in the context of uncertainty management is well established, and explains why it is now widely used. Having a probabilistic metamodel, in the sense that it provides a predictive distribution 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 GP-based approaches for active learning, robust optimization, reliability assessment, etc. In this context, it is essential to guarantee confidence in the GP predictive distribution, and not just in its mean value. This confidence requires, on the one hand, a reliable estimation of the GP metamodel and in particular of its hyperparameters, and, on the other hand, a rigorous validation of the entire GP predictive distribution.

The present paper has reviewed recent works dealing with the estimation of GP hyperparameters, from theoretical and practical points of view. It appears that the usual methods sometimes lead to poor-quality and not very robust estimates. MLE, the most widely used method, often leads to ill-posed problems. Although it leads in practice to good metamodel predictivity, the associated uncertainties and prediction intervals can be of poor quality. It is therefore essential to have validation indicators to detect this unreliability of the predictive distribution. Typically, it is insufficient to check only the GP's predictive capabilities: the accuracy of the entire GP predictive distribution needs to be assessed. To this end, we have reviewed the most relevant indicators and have proposed some derivatives. Emphasis has thus been put on GP validation that requires careful and informed consideration.

Concerning the estimation process, recent alternatives to standard estimation approaches have been explained. In particular, Bayesian approaches are theoretically very attractive, offering a kind of regularization of likelihood. However their cost in terms of complexity and required expertise, particularly in the definition of so-called robust priors and its tractability in large dimension (large number of inputs), refrain their use. Others approaches rely on ad-hoc corrections of the quantiles of the GP predictive distribution to ensure reliable prediction intervals for a given level, but these approaches do not necessarily seem relevant to our application context of multi-objective use of the metamodel. As an illustration, the RobustGaSP method is applied to the analytical example of the Branin function proposed in Section 4.4, with the approximated reference prior and the default parameter values of Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF]'s algorithm (package from [START_REF] Gu | RobustGaSP: Robust Gaussian Stochastic Process Emulation[END_REF]). The estimated hyperparameters are indicated on the middle plots of Figure 4 (by a magenta square) and the performances of the GP metamodel obtained with these hyperparameters are illustrated by the plots on the right. In the companion paper [START_REF] Marrel | Probabilistic surrogate modeling by Gaussian process: A new estimation algorithm for more robust prediction[END_REF], a new technique of using MLE for estimation is proposed, in particular by considering other criteria in the estimation procedure (criteria that until now have been reserved for the validation procedure). It also includes an intensive benchmark to test this new multi-objective optimization algorithm and compares its results with those obtained with other more standard estimation algorithms. The method is then applied on a real test case modeling an aquatic ecosystem, and used for environmental assessment.

More generally, beyond a more reliable estimation of hyperparameters, several challenges remain to enable an even wider use of the GP metamodel. The first is its estimation for very high dimensional problems (for example beyond 30 to 50 input variables), without using prior reduction techniques (e.g. preliminary screening of input variables). To achieve it and bypass the problem of MLE in such large dimension, Appriou et al. [START_REF] Appriou | Combination of optimization-free kriging models for highdimensional problems[END_REF] propose to consider a weighted combination of several GP regression sub-models with fixed hyperparameters. It will be interesting to further evaluate the performance of this approach and compare it to a standard approach (combined with a preliminary screening) on a benchmark of various test functions.

Extending the GP to inputs of a more complex nature is also an important prospect. For example, in solid and fluid mechanics, many studies involve varying mesh of a 2D or 3D shape. Handling such complex and very high-dimensional inputs in GP regression calls for the definition of suitable covariance kernels. Ongoing work is focusing in particular on the use of graph kernels, trying on the one hand to reduce the complexity of the usual graph kernels, which is extremely penalizing in relation to the number of nodes, and on the other hand to adapt them to handle inputs with different edges and even a different number of nodes (between two model evaluations or simulations.

Finally, it often arises that simulated scenarios considered in industrial safety assessment are affected by physical threshold phenomena, leading to the occurrence of phenomenological bifurcations or at least strong irregularities in the simulation outputs. Extensions of the GP, such as non-stationary GP [START_REF] Sauer | Non-stationary Gaussian process surrogates[END_REF] or Treed GP [START_REF] Gramacy | Bayesian treed Gaussian process models with an application to computer modeling[END_REF], exist to deal with this type of non-stationarity, but most have been implemented and applied in small dimensions (one or two input variables). The challenge would be to extend these methods to a larger dimension, more representative of real cases.

The optimization problem is finally reformulated as a one-dimensional problem 1 with d 2 ((σ 2 , θ) , ( σ 2 0 , θ 0 )) = W 2 2 (N (m, K), N (m 0 , K 0 )). The purpose is to find the two solutions ζ * and ζ * of P ζ , for a = (1 + α)/2 and a = (1 -α)/2 respectively. These two solutions yields two GP metamodels which are used to build the two bounds of the predictive interval whose coverage probability is demonstrated to be optimal. The resulting method is called Robust Prediction Intervals Estimation (RPIE) by the authors. The numerical tests propsed by the authors show that when the GP metamodel is well specified (good coverage probability of intervals obtained with initial MLE or CV-estimates θ 0 ), the RPIE method does not bring any added value. In the opposite case, the RPIE method is relevant and corrects efficiently the prediction intervals.

Appendix C. Details on prior and marginal posterior for the RobustGaSP method

To simplify, we suppose a constant GP mean m(x) = β 0 . The formulas generalized to the case of a q-dimensional vector of basis functions with parameter β ∈ R q+1 for m(x) are available in Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF].

To define their "robust" prior, Gu et al. [START_REF] Gu | Robust Gaussian stochastic process emulation[END_REF] first assign the objective prior for the regression and variance parameters:

π β 0 , σ 2 ∝ 1 (σ 2 ) a , (C.1)
with a > 0. a = 1 corresponds to the standard reference prior. Then, the authors consider the Jeffrey's prior which is a reference prior for their parametric model (separable product of 1-D correlation functions as in Eq. ( 6)):

π R β 0 , σ 2 , θ ∝ |I(θ)| 1/2 (σ 2 ) a , (C.2)
where I(•) ∈ R d×d is the expected Fisher information matrix as below, ,

I(θ) =          n -
where

W l = δR θ δθ l R -1 θ P R θ with P R θ = I n -1 n (1 n ) T R -1 θ 1 n -1 (1 n ) T R -1 θ ,

Figure 1 :

 1 Figure 1: Metamodel and probabilistic metamodel principles, illustrated with GP mean (GP metamodel predictor) and GP 90%-prediction intervals (quantiles of GP predictive distribution).

  in Section 4.4. Then, Gu et al. [40, Definition 3.1 and Lemma 3.2] propose to define the lack of robustness of GP hyperparameters by the occurrence of two cases:

1 Find

 1 the value of hyperparameters that minimize the MSE of GP predictions computed by CV (w.r.t. the observations) Use of fast CV Dubrule formulas to reduce the computation cost from 𝑂(𝐾𝑛 3 ) to 𝑂(𝑛 3 ) Theoretical results under the hypothesis of well-specified GP model. In practice, more robust (than MLE) when the Gaussian process model is misspecified No closed-form expression for the optimal values, numerical optimization required Few work on efficient computation of gradients Local minima More sensitive to the smoothness of function than MLE Ill-posed in non-asymptotic setting (finite dataset), as MLE under the GP model assumption Available gradients Theoretical results on the consistency of MLE under the hypothesis of well-specified GP model In practice, outperforms CV when the model is well-specified Close connection between MLE and CV Ills-posedness mitigatable by nugget effect and hyperparameter bounding No closed-form expression for the optimal values, numerical optimization required Computation cost of MLE and its gradients w.r.t. (𝜎 2 , 𝜽) in 𝑂(𝑛 3 + 𝑑𝑛 2 ) for a tensorized covariance with 𝑑 hyperparameters Local maxima and flatness in tails Ill-posed in non-asymptotic setting (finite dataset), as CV Bayesian Section 3.4

Figure 2 :

 2 Figure 2: Re-scaled Branin function on [0, 1] 2 (surface plot on the left and contour plot with isolines on the right).

Figure 3 :

 3 Figure 3: MBranin Function -Impact of hyperparameters values on graphical and quantitative validation criteria for a GP with a tensorized Gaussian covariance and conditioned by a fixed random sample of n = 30 simulations. Each column corresponds to a specific set of hyperparameters (specified in the colored frame at the top). The first line plots the predicted values of the metamodel against the true values of MBranin (on a test sample of 1000 points). The second line shows the α-plot of prediction intervals (also computed on the test sample).

Figure 4 :

 4 Figure 4: MBranin Function -Validation criteria obtained with MLE (left column) and RobustGaSP (right column) estimation with the same GP parameterization and training sample than results of Figure 3. The objective function given by Eq. (11) derived from likelihood objective function is plotted in the middle (surface plot above and contour plot with isolines below). Global optimal values, and MLE and RobustGaSP estimates, are plotted with squared.

  P ζ : P ζ : argmin ζ∈(0,+∞) L(ζ) := d 2 σ 2 opt (ζ), ζ θ 0 , σ 2 0 , θ 0 (B.3)

θ 1 n - 1 2 (S 2 )-( n- 1 2

 121 for 1 ≤ l ≤ d. The proofs of Eq. (C.2) and of the formula of I(•) are given in Paulo [74] (precisely in Proposition 2.1 and Appendix A.0.2).The prior on mean and variance parameter (Eq. (C.1)) allows to marginalize out these parameters in the likelihood function to obtain the marginal likelihood according to θ: L (Y s | θ). The marginal posterior of θ is then obtained by:π (θ | Y s ) ∝ L (Y s | θ) |I(θ)| 1/2 (C.3)whereL (Y s | θ) ∝ |R θ | -1 2 (1 n ) T R -1 +a-1) and S 2 = (Y s ) T R -1 θ P R θ Y s .1 Under some hypothesis on the GP trend, the authors demonstrate for all λ > 0,H ζ (λ) := σ 2 ∈ [0, +∞), ψ (δ) a σ 2 ,ζθ 0 = a is a non-empty and compact subset of R + . Assuming additional assumption of regularity of H ζ (λ), they deduce that it provides the continuity of σ 2 opt on (0, +∞) and that P ζ admits at least one global minimizer ζ * in (0, +∞).

Table 1 :

 1 Usual covariance functions and properties of associated GP trajectories.

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.

Appendix A. Synthesis on main validation criteria

Criterion

Prediction variances are of the same order of the observed prediction errors (for new location points). For example, a PVA≈0.2 corresponds to prediction variances around 20% too large or too small w.r.t the prediction errors.

High value

Unreliable prediction variances (and so confidence intervals): "overconfident" or "underconfident" metamodel, to be distinguished thanks to 𝛼-CI plot. For example, a PVA≈0.7 (resp. 1.2) corresponds to a prediction variances in average 2 (resp. 3) times too large or too small w.r.t the prediction errors.

Logarithmic score

LogS ≥ 0 (To be minimized)

Sharpness and variance of GP predictive distribution

Value close to 0 Observed values for unobserved location points have a high density w.r.t. GP predictive distribution

High value

No quantitative or absolute interpretation, only relative between several metamodels for example (and on the same test basis). Too small predictive variances (overconfident metamodel) are more penalized than too large ones.

Continuous ranked probability score

CRPS ≥ 0 (To be minimized)

Reliability of GP predictive intervals

Value close to 0 Observed values for unobserved location points are closed to the median (here mean for GP) of predictive distribution and sharpness of GP predictive distribution. So, in a nutshell, predictive and confident metamodel.

High value

No quantitative or absolute interpretation, only relative between several metamodels for example (and on the same test basis). Tends to favor metamodels with small predictive variance, subject to similar predictive performance.

Integrated Absolute Error

of 𝜶-PI plot IAE𝛼 ∈ [0 , 0.5] (To be minimized)

Reliability of GP predictive intervals

Value close to 0 Only if Q² is also high, reliable predicted confidence intervals Value close to 0.5 Unreliable predicted confidence intervals ("underconfident" or "overconfident" model, invalid Gaussianity hypothesis). The origin of the lack of reliability must be found in the cross-interpretation of Q², PVA and 𝛼-PI plot.

Table A.3: Summary of the main validation criteria with the purpose of the validation and details of interpretation.

Appendix B. Details on RPIE method

From ∆(α) given by Eq. ( 19), the authors introduce what they call the "quasi-Gaussian proportion" ψ a to describe how close the a-quantile q a of the standardized predictive distribution is to the level a (ideally, it should correspond to a). More precisely, ψ a results from the rewriting of ∆(α) for GP as

, where q a denotes the a-quantile of the standard normal distribution and h is the Heaviside step function

Considering a nominal level of quantile a, ψ a is then defined as a map from [0, +∞) × (0, +∞) d to [0, 1]:

.

The objective is then to find the two sets of hyperparameters (σ 2 , θ) such as to obtain the ideal values of ψ a for the two bounds of the α predictive interval. These two pairs denoted σ2 , θ and (σ 2 , θ) are respectively defined by

These GP parameters allow to get the optimal value ∆(α) = α. The authors modify the function ψ a into a new ψ (δ) a (with δ > 0 to obtain one formulation for the optimization problem, whatever the value of α, and define the set of solution A a,δ w.r.t. ∆(α):

Finally, as a correction of the initial estimated hyperparameters ( σ 2 0 , θ 0 ), the authors propose to find the hyperparameters in A a,δ which are the closest to ( σ 2 0 , θ 0 ) in the sense of a continuous similarity measure d sim between the multivariate Gaussian distributions N (m(X s ), K) generated with the two sets of parameters. This results in the following optimization problem:

This problem should then be solved for a = (1 + α)/2 and a = (1 -α)/2 to obtain estimates of σ2 , θ and (σ 2 , θ), respectively. However, two problems arise. First, the resolution of the problem (B.1) may be too costly and heavy to solve, especially as dimension d increases. Secondly, depending on the metric chosen, there is no guarantee that the barycenters of the two prediction intervals (generated from the two GP predictive distributions) are close. Similarly, there is no control over the sharpness of the prediction intervals obtained with the optimal solutions of Problem (B.1). This could result in a solution with a good coverage function but a poor Q 2 , cf. Remark 1. To mitigate the second drawback, the authors first recommend the use of second Wasserstein distance W 2 for d sim . Then, to address the two drawbacks, the authors propose a relaxed problem denoted P ζ where the optimal θ is defined as shifted values of θ 0 : ζ θ 0 with ζ > 0. Moreover, considering that σ 2 should be as small as possible to reduce the uncertainty of the GP predictions, the author define the optimal value of σ