Elias Fekhari
email: elias.fekhari@edf.fr

Bertrand Iooss
email: bertrand.iooss@edf.fr

Joseph Muré
email: joseph.mure@edf.fr

Luc Pronzato
email: luc.pronzato@i3s.unice.fr

Maria-João Rendas
email: ren-das@i3s.unice.fr

Luc +33130877969

Pronzato

Model predictivity assessment: incremental test-set selection and accuracy evaluation

Keywords: Design of experiments, Discrepancy, Gaussian process, Machine learning, Metamodel, Validation

Unbiased assessment of the predictivity of models learnt by supervised machine learning (ML) methods requires knowledge of the learned function over a reserved test set (not used by the learning algorithm). The quality of the assessment depends, naturally, on the properties of the test set and on the error statistic used to estimate the prediction error. In this work we tackle both issues, proposing a new predictivity criterion that carefully weights the individual observed errors to obtain a global error estimate, and using incremental experimental design methods to "optimally" select the test points on which the criterion is computed. Several incremental constructions are studied, including greedy-packing (coffee-house design), support points and kernel herding techniques. Our results show that the incremental and weighted versions of the latter two, based on Maximum Mean Discrepancy concepts, yield superior performance. An industrial test case provided by the historical French electricity supplier (EDF) illustrates the practical relevance of the methodology, indicating that it is an efficient alternative to expensive cross-validation techniques.

Introduction

The development of tools for automatic diagnosis relying on learned models imposes strict requirements on model validation. For example, in industrial non-destructive testing (e.g. for the aeronautic or the nuclear industry), generalized automated inspection, which increases efficiency and lowers costs, must provide high performance guarantees [START_REF] Eniq | Qualification of an AI / ML NDT system -Technical basis[END_REF][START_REF] Hawkins | Guidance on the assurance of machine learning in autonomous systems (AMLAS)[END_REF]. Establishing these guarantees requires availability of a reserved test set, i.e. a data set that has not been used either to train or to select the machine learning (ML) model [START_REF] Borovicka | Selecting representative data sets[END_REF][START_REF] Xu | On splitting training and validation set: A comparative study of cross-validation, bootstrap and systematic sampling for estimating the generalization performance of supervised learning[END_REF][START_REF] Iooss | Sample selection from a given dataset to validate machine learning models[END_REF]. Using the prediction residuals on this test set, an independent evaluation of the proposed ML model can be done, enabling the estimation of relevant performance metrics, such as the mean-squared error for regression problems, or the misclassification rate for classification problems.

The same need for independent test sets arises in the area of computer experiments, where computationally expensive simulation codes are often advantageously replaced by ML models, called surrogate models (or metamodels) in this context [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF][START_REF] Fang | Design and Modeling for Computer Experiments[END_REF]. Such surrogate models can be used, for instance, to estimate the region of the input space that maps to specific values of the model outputs [START_REF] Chevalier | Fast krigingbased stepwise uncertainty reduction with application to the identification of an excursion set[END_REF] with a significantly decreased computational load when compared to direct use of the original simulation code. Validation of these surrogate models consists in estimating their predictivity, and can either rely on a suitably selected validation sample, or be done by cross-validation [START_REF] Kleijnen | A methodology for fitting and validating metamodels in simulation[END_REF][START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF][START_REF] Demay | Model selection for Gaussian Process regression: an application with highlights on the model variance validation[END_REF]]. One of the numerical studies presented in this paper will address an example of this situation of practical industrial interest in the domain of nuclear safety assessment, concerning the simulation of thermal-hydraulic phenomena inside nuclear pressurized water reactors, for which finely validated surrogate models have demonstrated their usefulness [START_REF] Lorenzo | Assessment of an isolation condenser of an integral reactor in view of uncertainties in engineering parameters[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].

In this paper, we present methods to choose a "good" test set, either within a given dataset or within the input space of the model, as recently motivated in [START_REF] Iooss | Sample selection from a given dataset to validate machine learning models[END_REF][START_REF] Joseph | SPlit: An optimal method for data splitting[END_REF]. A first choice concerns the size of the test set. No optimal choice exists, and, when only a finite dataset is available, classical ML handbooks [START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Goodfellow | Deep learning[END_REF] provide different heuristics on how to split it, e.g., 80%/20% between the training and test samples, or 50%/25%/25% between the training, validation (used for model selection) and test samples. We shall not formally address this point here (see [START_REF] Xu | On splitting training and validation set: A comparative study of cross-validation, bootstrap and systematic sampling for estimating the generalization performance of supervised learning[END_REF] for a numerical study of this issue), but in the industrial case-study mentioned above we do study the impact of the ratio between the sizes of the training and test sets on the ability of assessing the quality of the surrogate model. A second issue concerns how the test sample is picked within the input space. The simplest -and most common -way to build a test sample is to extract an independent Monte Carlo sample [START_REF] Hastie | The elements of statistical learning[END_REF]. For small test sets, these randomly chosen points may fall too close to the training points or leave large areas of the input space unsampled, and a more constructive method to select points inside the input domain is therefore preferable. Similar concerns motivate the use of space-filling designs when choosing a small set of runs for cpu-time expensive computer experiments on which a model will be identified [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF][START_REF] Pronzato | Design of computer experiments: space filling and beyond[END_REF].

When the test set must be a subset of an initial dataset, the problem amounts to selecting a certain number of points within a finite collection of points. A review of classical methods for solving this issue is given in [START_REF] Borovicka | Selecting representative data sets[END_REF]. For example, the CADEX and DUPLEX algorithms [START_REF] Kennard | Computer aided design of experiments[END_REF][START_REF] Snee | Validation of regression models: Methods and examples[END_REF] can sequentially extract points from a database to include them in a test sample, using an inter-point distance criterion. In ML, identifying within the dataset "prototypes" (set of data instances representative of the whole data set) and "criticisms" (data instances poorly represented by the prototypes) has recently been proposed to help model interpretation [START_REF] Molnar | Interpretable machine learning[END_REF]; the extraction of prototypes and criticisms relies on a Maximum Mean Discrepancy (MMD) criterion [START_REF] Smola | A Hilbert space embedding for distributions[END_REF] (see e.g. [START_REF] Pronzato | Bayesian quadrature and energy minimization for spacefilling design[END_REF], and [START_REF] Pronzato | Performance analysis of greedy algorithms for minimising a maximum mean discrepancy[END_REF] for a performance analysis of greedy algorithms for MMD minimization).

Several algorithms have also been proposed for the case where points need to be added to an already existing training sample. When the goal is to assess the quality of a model learnt using a known training set, one may be tempted to locate the test points the furthest away from the training samples, such that, in some sense, the union of the training and test sets is space-filling. As this paper shows, test sets built in this manner do enable a good assessment of the quality of models learnt with the training set if the observed residuals are appropriately weighted. Moreover, the incremental augmentation of a design can be useful when the assessed model turns out to be of poor quality, or when an additional computational budget is available after a first study [START_REF] Sheikholeslami | Progressive Latin hypercube sampling: An efficient approach for robust sampling-based analysis of environmental models[END_REF][START_REF] Shang | Fully-sequential space-filling design algorithms for computer experiments[END_REF]. Different empirical strategies have been proposed for incremental space-filling design [START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF][START_REF] Crombecq | Efficient space-filling and non-collapsing sequential design strategies for simulation-based modelling[END_REF][START_REF] Li | A novel extension algorithm for optimized Latin hypercube sampling[END_REF], which basically entail the addition of new points in the zones poorly covered by the current design. Shang and Apley [START_REF] Shang | Fully-sequential space-filling design algorithms for computer experiments[END_REF] have recently proposed an improvement of the CADEX algorithm, called the Fully-sequential space-filling (FSSF) design; see also [START_REF] Nogales Gómez | Incremental space-filling design based on coverings and spacings: improving upon low discrepancy sequences[END_REF] for an alternative version of coffee-house design enforcing boundary avoidance. Although they are developed for different purposes, nested space filling designs [START_REF] Qian | Construction of nested space-filling designs[END_REF] and sliced space filling designs [START_REF] Qian | Sliced space filling designs[END_REF] can also be used to build sequential designs.

In this work, we provide new insights into these subjects in two main directions: (i) definition of new predictivity criteria through an optimal weighting of the test points residuals, and (ii) use of test sets built by incremental space-filling algorithms, namely FSSF, support points [START_REF] Mak | Support points[END_REF] and kernel herding [START_REF] Chen | Super-samples from kernel herding[END_REF], the latter two algorithms being typically used to provide a representative sample of a desired theoretical or empirical distribution. Besides, this paper presents a numerical benchmark analysis comparing the behaviour of the three algorithms on a selected set of test cases.

This paper is organized as follows. Section 2 defines the predictivity criterion considered and proposes different methods for its estimation. Section 3 presents the three algorithms used for test-point selection: FSSF, support points and kernel herding. Our numerical results are presented in Sections 4 and 5: in Section 4 a test set is freely chosen within the entire input space, while in Section 5 an existing data set can be split into a training sample and a test set. Finally, Section 6 concludes and outlines some perspectives.

Predictivity assessment criteria for an ML model

In this section, we propose a new criterion to assess the predictive performance of a model, derived from a standard model quality metric by suitably weighting the errors observed on the test set. We denote by X ⊂ R 𝑑 the space of the input variables x = (𝑥 1 , . . . , 𝑥 𝑑 ) of the model. Then let 𝑦(x) ∈ R (resp. 𝑦(x ) ∈ R) be the observed output at point x ∈ X (resp. x ∈ X). We denote by (X 𝑚 , y 𝑚 ) the training sample, with y 𝑚 = [𝑦(x (1) ), . . . , 𝑦(x (𝑚) )] . The test sample is denoted by (X 𝑛 , y 𝑛 ) = (x (𝑚+𝑖) , 𝑦(x (𝑚+𝑖) )) 1≤𝑖 ≤𝑛 .

The predictivity coefficient

Let 𝜂 𝑚 (x) denote the prediction at point x of a model learned using (X 𝑚 , y 𝑚 ) [START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. A classical measure for assessing the predictive ability of 𝜂 𝑚 , in order to evaluate its validity, is the predictivity coefficient. Let 𝜇 denote the measure that weights how comparatively important it is to accurately predict 𝑦 over the different regions of X. For example the input could be a random vector with known distribution: in that case, this distribution would be a reasonable choice for 𝜇. The true (ideal) value of the predictivity is defined as the following normalization of the Integrated Square Error (ISE):

𝑄 2 ideal (𝜇) = 1 - ISE 𝜇 (X 𝑚 , y 𝑚 ) 𝑉 𝜇 , (1)
where

ISE 𝜇 (X 𝑚 , y 𝑚 ) = ∫ X [𝑦(x) -𝜂 𝑚 (x)] 2 d𝜇(x) , 𝑉 𝜇 = ∫ X 𝑦(x) - ∫ X 𝑦(x )𝑑𝜇(x ) 2 d𝜇(x) .
The ideal predictivity 𝑄 2 ideal (𝜇) is usually estimated by its empirical version calculated over the test sample (X 𝑛 , y 𝑛 ), see [10, p. 32]:

𝑄 2 𝑛 = 1 - x∈X 𝑛 [𝑦(x) -𝜂 𝑚 (x)] 2 x∈X 𝑛 𝑦(x) -𝑦 𝑛 2 , (2)
where 𝑦 𝑛 = (1/𝑛) 𝑛 𝑖=1 𝑦(x (𝑚+𝑖) ) denotes the empirical mean of the observations in the test sample. Note that the calculation of 𝑄 2 𝑛 only requires access to the predictor 𝜂 𝑚 (•). To compute 𝑄 2 𝑛 , one does not need to know the training set which was used to build 𝜂 𝑚 (•). 𝑄 2 𝑛 is the coefficient of determination (a standard notion in parametric regression) common in prediction studies [START_REF] Kleijnen | A methodology for fitting and validating metamodels in simulation[END_REF][START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF], often called "Nash-Sutcliffe criterion" [START_REF] Nash | River flow forecasting through conceptual models part I-A discussion of principles[END_REF]: it compares the prediction errors obtained with the model 𝜂 𝑚 with those obtained when prediction equals the empirical mean of the observations. Thus, the closer 𝑄 2 𝑛 is to one, the more accurate the surrogate model is (for the test set considered). On the contrary, 𝑄 2 𝑛 close to zero (negative values are possible too) indicates poor predictions abilities, as there is little improvement compared to prediction by the simple empirical mean of the observations. The next section shows how a suitable weighting of the residual on the training sample may be key to improving the estimation of 𝑄 2 𝑛 .

Weighting the test sample

Let 𝜉 𝑛 = 𝑛 -1 𝑛 𝑖=1 𝛿 x (𝑖) be the empirical distribution of the prediction error, with 𝛿 x the Dirac measure at x. In 𝑄 2 𝑛 , ISE 𝜇 (X 𝑚 , y 𝑚 ) is estimated by the empirical average of the squared residuals

ISE 𝜉 𝑛 (X 𝑚 , y 𝑚 ) = 1 𝑛 𝑛 ∑ 𝑖=1 𝑦(x (𝑚+𝑖) ) -𝜂 𝑚 (x (𝑚+𝑖) ) 2 .
When the points x (𝑚+𝑖) of the test set X 𝑛 are distant from the points of the training set X 𝑚 , the squared prediction errors |𝑦(x (𝑚+𝑖) ) -𝜂 𝑚 (x (𝑚+𝑖) )| 2 tend to represent the worst possible error situations, and ISE 𝜉 𝑛 (X 𝑚 , y 𝑚 ) tends to overestimate ISE 𝜇 (X 𝑚 , y 𝑚 ). In this section, we postulate a statistical model for the prediction errors in order to be able to quantify this potential bias when sampling the residual process, enabling its subsequent correction .

In [START_REF] Pronzato | Validation design I: construction of validation designs via kernel herding[END_REF], the authors propose a weighting scheme for the test set when the ML model interpolates the train set observations. They suggest several variants corresponding to different constraints on the weights (e.g., non-negativity, summing to one). In the following, we consider the unconstrained version only, which in our experience works best. Let 𝛿 𝑚 (x) = 𝑦(x) -𝜂 𝑚 (x) denote the predictor error and assume it is a realization of a Gaussian Process (GP) with zero mean and covariance kernel 𝜎 2 𝐾 |𝑚 , which we shall note 𝛿 𝑚 (x) ∼ GP(0, 𝜎 2 𝐾 |𝑚 ), with

𝜎 2 𝐾 |𝑚 (x, x ) = E{𝛿 𝑚 (x)𝛿 𝑚 (x )} = 𝜎 2 𝐾 (x, x ) -k 𝑚 (x)K -1 𝑚 k 𝑚 (x ) .
Here, k 𝑚 (x) denotes the column vector [𝐾 (x, x (1) ), . . . , 𝐾 (x, x (𝑚) )] and K 𝑚 is the 𝑚 × 𝑚 matrix whose element (𝑖, 𝑗) is given by {K 𝑚 } 𝑖, 𝑗 = 𝐾 (x (𝑖) , x ( 𝑗) ), with 𝐾 a positive definite kernel. The rationale for using this model is simple: assume first a prior GP model GP(0, 𝜎 2 𝐾) for the error process 𝛿(x); if 𝜂 𝑚 interpolates the observations y 𝑚 , the errors observed at the design points x (𝑖) equal zero, 𝑖 = 1, . . . , 𝑚, leading finally to the posterior GP(0, 𝜎 2 𝐾 |𝑚 ) for 𝛿 𝑚 (x). However, the predictor 𝜂 𝑚 is not always an interpolator, see Section 5 for an example, so we extend the approach of [START_REF] Pronzato | Validation design I: construction of validation designs via kernel herding[END_REF] to the general situation where 𝜂 𝑚 does not necessarily interpolate the training data y 𝑚 . The same prior GP(0, 𝜎 2 𝐾) for 𝛿(x) yields 𝛿 𝑚 (x) ∼ GP( 𝛿 𝑚 (x), 𝜎 2 𝐾 |𝑚 ), where

𝛿 𝑚 (x) = k 𝑚 (x)K -1 𝑚 (y 𝑚 -𝜂 𝑚 ) (3)
is the Kriging interpolator for the errors, with 𝜂 𝑚 = [𝜂 𝑚 (x (1) ), . . . , 𝜂 𝑚 (x

) = E ISE 𝜉 𝑛 (X 𝑚 , y 𝑚 ) -ISE 𝜇 (X 𝑚 , y 𝑚 ) 2 = E ∫ X 𝛿 2 𝑚 (x) d(𝜉 𝑛 -𝜇) (x) 2 = E ∫ X 2 𝛿 2 𝑚 (x)𝛿 2 𝑚 (x ) d(𝜉 𝑛 -𝜇) (x)d(𝜉 𝑛 -𝜇) (x ) .
Tonelli's theorem gives

Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚 ) = ∫ X 2 E{𝛿 2 𝑚 (x)𝛿 2 𝑚 (x )} d(𝜉 𝑛 -𝜇) (x)d(𝜉 𝑛 -𝜇) (x ) . Since E 𝑈 2 𝑉 2 = 2 (E{𝑈𝑉 }) 2 + E 𝑈 2 E 𝑉 2
for any one-dimensional normal centered random variables 𝑈 and 𝑉, when 𝜂 𝑚 (x) interpolates y 𝑚 , we obtain

Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚 ) = 𝜎 2 𝑑 2 𝐾 |𝑚 (𝜉 𝑛 , 𝜇) , (4)
where

𝑑 2 𝐾 |𝑚 (𝜉 𝑛 , 𝜇) = ∫ X 2 𝐾 |𝑚 (x, x ) d(𝜉 𝑛 -𝜇) (x)d(𝜉 𝑛 -𝜇) (x ) , 𝐾 |𝑚 (x, x ) = 2 𝐾 2 |𝑚 (x, x ) + 𝐾 |𝑚 (x, x)𝐾 |𝑚 (x , x ) , (5)
and we recognize 𝑑

w * 𝑛 = K -1 |𝑚 (X 𝑛 )p 𝐾 |𝑚 ,𝜇 (X 𝑛 ) . (6)
We shall denote by 𝜁 * 𝑛 the measure supported on X 𝑛 with the optimal weights (6) and

𝑄 2 𝑛 * = 1 - ISE 𝜁 * 𝑛 (X 𝑚 , y 𝑚 ) 1 𝑛 𝑛 𝑖=1 𝑦(x (𝑚+𝑖) ) -𝑦 𝑛 2 = 1 - 𝑛 𝑖=1 𝑤 * 𝑖 𝑦(x (𝑚+𝑖) ) -𝜂 𝑚 (x (𝑚+𝑖) ) 2 1 𝑛 𝑛 𝑖=1 𝑦(x (𝑚+𝑖) ) -𝑦 𝑛 2 , (7)
with ) ). Notice that the weights 𝑤 * 𝑖 do not depend on the variance parameter 𝜎 2 of the GP model.

𝑦 𝑛 = (1/𝑛) 𝑛 𝑖=1 𝑦(x (𝑚+𝑖
Remark 1 When X 𝑛 is constructed by kernel herding, see Section 3.3, 𝐾 can be chosen identical to the kernel used there. This will be the case in Sections 4 and 5, but it is not mandatory.

Conversely, one may think of choosing a design X 𝑛 that minimizes 𝑑 (𝜁 * 𝑛 , 𝜇), also with the objective to obtain a precise estimation of ISE 𝜇 (X 𝑚 , y 𝑚 ). This can be achieved by kernel herding using the kernel 𝐾 |𝑚 and is addressed in [START_REF] Pronzato | Validation design I: construction of validation designs via kernel herding[END_REF]. However, the numerical results presented there show that the precise choice of the test set X 𝑛 has a marginal effect compared to the effect of non-uniform weighting with w * 𝑛 , provided that X 𝑛 fills the holes left in X by the training design X 𝑚 .

⊳

Remark 2 When the observations 𝑦(x (𝑖) ), 𝑖 = 1, . . . , 𝑛, are available at the validation stage, an alternative version of 𝑄 2 𝑛 would be

𝑄 2 𝑛 = 1 - x∈X 𝑛 [𝑦(x) -𝜂 𝑚 (x)] 2 x∈X 𝑛 𝑦(x) -𝑦 𝑚 2 , (8)
where

𝑦 𝑚 = (1/𝑚) 𝑚 𝑖=1 𝑦(x 𝑖 )
, which compares the performance on the test set of two predictors 𝜂 𝑚 and 𝑦 𝑚 based on the same training set. It is then possible to also apply a weighting procedure to the denominator of 𝑄 2 𝑛 ,

𝐷 𝜉 𝑛 (X 𝑚 , y 𝑚 ) = 1 𝑛 𝑛 ∑ 𝑖=1 𝑦(x (𝑚+𝑖) ) -𝑦 𝑚 2 ,
in order to make it resemble its idealized version

𝑉 𝜇 (y 𝑚 ) = ∫ X 𝑦(x) -𝑦 𝑚 2 d𝜇(x).
The GP model is now

𝜀 𝑚 (x) = 𝑦(x) -𝑦 𝑚 ∼ GP(𝜂 𝑚 (x) -𝑦 𝑚 , 𝜎 2 𝐾 |𝑚 )
. Similar developments to those used above for the numerator of 𝑄 2 𝑛 yield

Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚 ) = E 𝐷 𝜉 𝑛 (X 𝑚 , y 𝑚 ) -𝑉 𝜇 (y 𝑚 ) 2 = E ∫ X 2 𝜀 2 𝑚 (x)𝜀 2 𝑚 (x ) d(𝜉 𝑛 -𝜇) (x)d(𝜉 𝑛 -𝜇) (x ) = 𝜎 2 𝑑 2 𝐾 |𝑚 (𝜉 𝑛 , 𝜇) ,
where

𝐾 |𝑚 (x, x ) = 𝐾 |𝑚 (x, x ) + [𝜂 𝑚 (x) -y 𝑚 ] 2 [𝜂 𝑚 (x ) -y 𝑚 ] 2 + [𝜂 𝑚 (x) -y 𝑚 ] 2 𝐾 |𝑚 (x , x ) + [𝜂 𝑚 (x ) -y 𝑚 ] 2 𝐾 |𝑚 (x, x) + 4 [𝜂 𝑚 (x) -y 𝑚 ] [𝜂 𝑚 (x ) -y 𝑚 ]𝐾 |𝑚 (x, x ) .
We can then substitute 𝐷 𝜁 𝑛 * (X 𝑚 , y 𝑚 ) for 𝐷 𝜉 𝑛 (X 𝑚 , y 𝑚 ) in ( 8), where 𝜁 𝑛 * allocates the

weights w 𝑛 * = K -1 |𝑚 (X 𝑛 )p 𝐾 |𝑚 ,𝜇 (X 𝑛 ) to the 𝑛 points in X 𝑛 , with {K |𝑚 (X 𝑛 )} 𝑖, 𝑗 = 𝐾 |𝑚 (x (𝑚+𝑖) , x (𝑚+ 𝑗) ), 𝑖, 𝑗 = 1, . . . , 𝑛. ⊳

Test-set construction

In the previous section we assumed the test set as given, and proposed a method to estimate ISE 𝜇 (X 𝑚 , y 𝑚 ) by a weighted sum of the residuals. In this section we address the choice of the test set. Below we give an overview of the three methods used in this paper, all relying on the concept of space-filling design [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF][START_REF] Pronzato | Design of computer experiments: space filling and beyond[END_REF]. While most methods for the construction of such designs choose all points simultaneously, the methods we consider are incremental, selecting one point at a time.

Our objective is to construct an ordered test set of size 𝑛, denoted by X 𝑛 = x (1) , . . . , x (𝑛) ⊂ X. When there is no restriction on the choice of X 𝑛 , the advantage of using an incremental construction is that it can be stopped once the estimation of the predictivity of an initial model, built with some given design X 𝑚 , is considered sufficiently accurate. In case the conclusion is that model predictions are not reliable enough, the full design X 𝑚+𝑛 = X 𝑚 ∪ X 𝑛 and the associated observations y 𝑚+𝑛 can be used to update the model. This updated model can then be tested at additional design points, elements of a new test set to be constructed. All methods presented in this section (except the Fully Sequential Space-Filling method) are implemented in the Python package otkerneldesign which is based on the OpenTURNS library for uncertainty quantification [START_REF] Baudin | Open TURNS: An industrial software for uncertainty quantification in simulation[END_REF].

Fully-Sequential Space-Filling design

The Fully-Sequential Space-Filling forward-reflected (FSSF-fr) algorithm [START_REF] Shang | Fully-sequential space-filling design algorithms for computer experiments[END_REF] relies on the CADEX algorithm [START_REF] Kennard | Computer aided design of experiments[END_REF] (also called the "coffee-house" method [START_REF] Müller | Collecting Spatial Data[END_REF]). It constructs a sequence of nested designs in a bounded set X by sequentially selecting a new point x as far away as possible from the x (𝑖) previously selected. New inserted points are selected within a set of candidates S which may coincide with X or be a finite subset of X (which simplifies the implementation, only this case is considered here). The improvement of FSSF-fr when compared to CADEX is that new points are selected at the same time far from the previous design points as well as far from the boundary of X.

The algorithm is as follows:

1. Choose S, a finite set of candidate points in X, with size 𝑁 𝑛 in order to allow a fairly dense covering of X. When X = [0, 1] 𝑑 , [START_REF] Shang | Fully-sequential space-filling design algorithms for computer experiments[END_REF] recommends to take S equal to the first 𝑁 = 1 000 𝑑 + 2 𝑛 points of a Sobol sequence in X. 2. Choose the first point x (1) randomly in S and define X 1 = {x (1) }. 3. At iteration 𝑖, with X 𝑖 = {x (1) , . . . , x (𝑖) }, select

x (𝑖+1) ∈ Arg max x∈S\X 𝑖 min min 𝑗 ∈ {1,...,𝑖 } x -x ( 𝑗) , √ 2 𝑑 dist(x, 𝑅(x)) , (9)
where 𝑅(x) is the symmetric of x with respect to its nearest boundary of X, and set X 𝑖+1 = X 𝑖 x (𝑖+1) . 4. Stop the algorithm when X 𝑛 has the required size.

The standard coffee-house (greedy packing) algorithm simply uses x (𝑖+1) ∈ Arg max x∈S\X 𝑖 min 𝑗 ∈ {1,...,𝑖 } xx ( 𝑗) . The role of the reflected point 𝑅(x) is to avoid selecting x (𝑖+1) too close to the boundary of X, which is a major problem with standard coffee-house, especially when X = [0, 1] 𝑑 with 𝑑 large. The factor √ 2𝑑 in [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF] proposed in [START_REF] Shang | Fully-sequential space-filling design algorithms for computer experiments[END_REF] sets a balance between distance to the design X 𝑖 and distance to the boundary of X. Another scaling factor, depending on the target design size 𝑛 is proposed in [START_REF] Nogales Gómez | Incremental space-filling design based on coverings and spacings: improving upon low discrepancy sequences[END_REF].

FSSF-fr is entirely based on geometric considerations and implicitly assumes that the selected set of points should cover X evenly. However, in the context of uncertainty quantification [START_REF] Smith | Uncertainty quantification[END_REF] it frequently happens that the distribution 𝜇 of the model inputs is not uniform. It is then desirable to select a test set representative of 𝜇. This can be achieved through the inverse probability integral transform: FSSFfr constructs X 𝑛 in the unit hypercube [0, 1] 𝑑 , and an "isoprobabilistic" transform 𝑇 : [0, 1] 𝑑 → X is then applied to the points in X 𝑖 , 𝑇 being such that, if 𝑼 is a https://pypi.org/project/otkerneldesign/ random variable uniform on [0, 1] 𝑑 , then 𝑇 (𝑼) follows the target distribution 𝜇. The transformation can be applied to each input separately when 𝜇 is the product of its marginals, a situation considered in our second test-case of Section 4, but is more complicated in other cases, see [START_REF] Lemaire | Structural reliability[END_REF]Chap. 4]. Note that FSSF-fr operates in the bounded set [0, 1] 𝑑 even if the support of 𝜇 is unbounded. The other two algorithms presented in this section are able to directly choose points representative of a given distribution 𝜇 and do not need to resort to such a transformation.

Support points

Support points [START_REF] Mak | Support points[END_REF] are such that their associated empirical distribution 𝜉 𝑛 has minimum Maximum-Mean-Discrepancy (MMD) with respect to 𝜇 for the energydistance kernel of Székely and Rizzo [START_REF] Székely | Testing for equal distributions in high dimension[END_REF][START_REF] Székely | Energy statistics: A class of statistics based on distances[END_REF],

𝐾 𝐸 (x, x ) = 1 2 ( x + x -x -x ) . (10)
The squared MMD between 𝜉 𝑛 and 𝜇 for the distance kernel equals

𝑑 2 𝐾 𝐸 (𝜉 𝑛 , 𝜇) = 2 𝑛 𝑛 ∑ 𝑖=1 E x (𝑖) -𝜁 - 1 𝑛 2 𝑛 ∑ 𝑖=1 𝑛 ∑ 𝑗=1 x (𝑖) -x ( 𝑗) -E 𝜁 -𝜁 , (11)
where 𝜁 and 𝜁 are independently distributed with 𝜇; see [START_REF] Sejdinovic | Equivalence of distance-based and RKHS-based statistics in hypothesis testing[END_REF]. A key property of the energy-distance kernel is that it is characteristic [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF]: for any two probability distributions 𝜇 and 𝜉 on X, 𝑑 2 𝐾 𝐸 (𝜇, 𝜉) equals zero if and only if 𝜇 = 𝜉, and so it defines a norm in the space of probability distributions. Compared to more heuristic methods for solving quantization problems, support points benefit from the theoretical guarantees of MMD minimization in terms of convergence of 𝜉 𝑛 to 𝜇 as 𝑛 → ∞.

As E x (𝑖) -𝜁 is not known explicitly, in practice 𝜇 is replaced by its empirical version 𝜇 𝑁 for a given large-size sample (x (𝑘) ) 𝑘=1...𝑁 . The support points X 𝑠 𝑛 are then given by X 𝑠 𝑛 ∈ Arg min

x (1) ,...,x (𝑛)

2 𝑛𝑁 𝑛 ∑ 𝑖=1 𝑁 ∑ 𝑘=1 x (𝑖) -x (𝑘) - 1 𝑛 2 𝑛 ∑ 𝑖=1 𝑛 ∑ 𝑗=1 x (𝑖) -x ( 𝑗) . (12)
The function to be minimized can be written as a difference of functions convex in X 𝑛 , which yields a difference-of-convex program. In [START_REF] Mak | Support points[END_REF], a majorization-minimization procedure, efficiently combined with resampling, is applied to the construction of large designs (up to 𝑛 = 10 4 ) in high dimensional spaces (up to 𝑑 = 500). The examples treated clearly show that support points are distributed in a way that matches 𝜇 more closely than Monte-Carlo and quasi-Monte Carlo samples [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF]. The method can be used to split a dataset into a training set and a test set [START_REF] Joseph | SPlit: An optimal method for data splitting[END_REF]: the 𝑁 points X 𝑁 in [START_REF] Demay | Model selection for Gaussian Process regression: an application with highlights on the model variance validation[END_REF] are those from the dataset, X 𝑠 𝑛 gives the test set and the other 𝑁 -𝑛 points are used for training. There is a serious additional difficulty though, as choosing X 𝑠 𝑛 among the dataset corresponds to a difficult combinatorial optimization problem. A possible solution is to perform the optimization in a continuous domain X and then choose X 𝑠 𝑛 that corresponds to the closest points in X 𝑁 (for the Euclidean distance) to the continuous solution obtained [START_REF] Joseph | SPlit: An optimal method for data splitting[END_REF].

The direct determination of support points through [START_REF] Demay | Model selection for Gaussian Process regression: an application with highlights on the model variance validation[END_REF] does not allow the construction of a nested sequence of test sets. One possibility would be to solve [START_REF] Demay | Model selection for Gaussian Process regression: an application with highlights on the model variance validation[END_REF] sequentially, one point at a time, in a continuous domain, and then select the closest point within X 𝑁 as the one to be included in the test set. We shall use a different approach here, based on the greedy minimization of the MMD [START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF] for the candidate set S = X 𝑁 : at iteration 𝑖, the algorithm chooses

x 𝑠 𝑖+1 ∈ Arg min x∈S 1 𝑁 𝑁 ∑ 𝑘=1 x -x (𝑘) - 1 𝑖 + 1 𝑖 ∑ 𝑗=1 x -x ( 𝑗) . ( 13
)
The method requires the computation of the 𝑁 (𝑁 -1)/2 distances x (𝑖)x ( 𝑗) , 𝑖, 𝑗 = 1, . . . , 𝑁, 𝑖 ≠ 𝑗, which hinders its applicability to large-scale problems (a test-case with 𝑁 = 1 000 is presented in Section 5). Note that we consider support points in the input space X only, with X ⊆ R 𝑑 , in contrast with [START_REF] Joseph | SPlit: An optimal method for data splitting[END_REF] which considers couples (x (𝑖) , 𝑦(x (𝑖) )) in R 𝑑+1 to split a given dataset into a training set and a test set.

Greedy MMD minimization can be applied to other kernels than the distance kernel [START_REF] Veiga | Basics and Trends in Sensitivity Analysis[END_REF], see [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF][START_REF] Pronzato | Performance analysis of greedy algorithms for minimising a maximum mean discrepancy[END_REF]. In the next section we consider the closely related method of Kernel Herding (KH) [START_REF] Chen | Super-samples from kernel herding[END_REF], which corresponds to a conditional-gradient descent in the space of probability measures supported on a candidate set S; see, e.g., [START_REF] Pronzato | Bayesian quadrature and energy minimization for spacefilling design[END_REF] and the references therein.

Kernel herding

Let 𝐾 be a positive definite kernel on X × X. At iteration 𝑖 of kernel herding, with 𝜉 𝑖 = (1/𝑖) 𝑖 𝑗=1 𝛿 x (𝑖) the empirical measure for X 𝑖 , the next point x 𝑖+1 minimizes the directional derivative 𝐹 𝐾 (𝜉 𝑖 , 𝜇, 𝛿 x ) of the squared MMD 𝑑 2 𝐾 (𝜉, 𝜇) at 𝜉 = 𝜉 𝑖 in the direction of the delta measure 𝛿 x , see Appendix A. Direct calculation gives

𝐹 𝐾 (𝜉 𝑖 , 𝜇, 𝛿 x ) = 2 [𝑃 𝐾 , 𝜉 (x) -𝑃 𝐾 ,𝜇 (x)] -2 ∫ X [𝑃 𝐾 , 𝜉 (x) -𝑃 𝐾 ,𝜇 (x)
] d𝜉 (x), with 𝑃 𝐾 , 𝜉 (x) (resp. 𝑃 𝐾 , 𝜇 (x)) the potential of 𝜉 (resp. 𝜇) at x, see [START_REF] Goodfellow | Deep learning[END_REF], and thus

x 𝑖+1 ∈ Arg min x∈S 𝑃 𝐾 , 𝜉 𝑖 (x) -𝑃 𝐾 ,𝜇 (x) , (14)
with S ⊆ X a given candidate set. Here, 𝑃 𝐾 , 𝜉 𝑖 (x) = (1/𝑖) 𝑖 𝑗=1 𝐾 (x, x ( 𝑗) ). When an empirical measure 𝜇 𝑁 based on a sample (x (𝑘) ) 𝑘=1...𝑁 is substituted for 𝜇, we get 𝑃 𝐾 , 𝜇 𝑁 (x) = (1/𝑁) 𝑁 𝑘=1 𝐾 (x, x (𝑘) ), which gives

x 𝑖+1 ∈ Arg min x∈S       1 𝑖 𝑖 ∑ 𝑗=1 𝐾 (x, x ( 𝑗) ) - 1 𝑁 𝑁 ∑ 𝑘=1 𝐾 (x, x (𝑘) )       .
When 𝐾 is the energy-distance kernel [START_REF] Veiga | Basics and Trends in Sensitivity Analysis[END_REF] we thus obtain [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF] with a factor 1/𝑖 instead of 1/(𝑖 + 1) in the second sum.

The candidate set S in ( 14) is arbitrary and can be chosen as in Section 3.1. A neat advantage of kernel herding over support points in that the potential 𝑃 𝐾 ,𝜇 (x) is sometimes explicitly available. When S = X 𝑁 , this avoids the need to calculate the 𝑁 (𝑁 -1)/2 distances x (𝑖)x ( 𝑗) and thus allows application to very large sample sizes. This is the case in particular when X is the cross product of onedimensional sets X

[𝑖 ] , X = X [1] × • • • × X [𝑑 ] , 𝜇
𝐾 5/2, 𝜃 (𝑥 -𝑥 ) = 1 + √ 5 𝜃 |𝑥 -𝑥 | + 5 3𝜃 2 (𝑥 -𝑥 ) 2 exp - √ 5 𝜃 |𝑥 -𝑥 | , (15)
the one-dimensional potentials are given in Appendix B for 𝜇 [𝑖 ] uniform on [0, 1] or 𝜇 [𝑖 ] the standard normal 𝒩(0, 1). When no observation is available, which is the common situation at the design stage, the correlation lengths have to be set to heuristic values. We empirically found the values of the correlation lengths to have a large influence over the design. A reasonable choice for X = [0, 1] 𝑑 is 𝜃 𝑖 = 𝑛 -1/𝑑 for all 𝑖, with 𝑛 the target number of design points; see [START_REF] Pronzato | Bayesian quadrature and energy minimization for spacefilling design[END_REF].

Numerical illustration

We apply FSSF-fr (denoted FSSF in the following), support points and kernel herding algorithms to the situation where a given initial design of size 𝑚 has to be completed by a series of additional points x (𝑚+1) , . . . , x (𝑚+𝑛) . The objective is to obtain a full design X 𝑚+𝑛 that is a good quantization of a given distribution 𝜇.

Figures 1 and2 correspond to 𝜇 uniform on [0, 1] 2 and 𝜇 the standard normal distribution 𝒩(0, I 2 ), with I 2 the 2-dimensional identity matrix, respectively. All methods are applied to the same candidate set S.

The initial designs X 𝑚 are chosen in the class of space-filling designs, well suited to initialize sequential learning strategies [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. When 𝜇 is uniform, the initial design is a maximin Latin hypercube design [START_REF] Morris | Exploratory designs for computationnal experiments[END_REF] with 𝑚 = 10 and the candidate set is given by the 𝑁 = 2 12 first points S 𝑁 of a Sobol sequence in [0, 1]. When 𝜇 is normal, the inverse probability transform method is first applied to S 𝑁 and X 𝑚 (this does not raise any difficulty here as 𝜇 is the product of its marginals). The candidate points S are marked in gray on Figures 1 and2 and the initial design is indicated by the red crosses. The index 𝑖 of each added test point x (𝑚+𝑖) is indicated (the font size decreases with 𝑖). In such a small dimension (𝑑 = 2), a visual appreciation gives the impression that the three methods have comparable performance. We can notice, however, that FSSF tends to choose points closer to the boundary of S than the other two, and that support points seem to sample more freely the holes of X 𝑚 than kernel herding, which seems to be closer to a space-filling continuation of the training set. We will come back to these designs when analysing the quality of the resulting predictivity metric estimators in the next section.

Numerical results I: construction of a training set and a test set

This section presents numerical results obtained on three different test-cases, in dimension 2 (test-cases 1 and 2) and 8 (test-case 3), for which 𝑦(x) = 𝑓 (x) with 𝑓 (x) having an easy to evaluate analytical expression, see Section 4.1. This allows a good estimation of 𝑄 2 ideal (𝜇) by 𝑄 2 𝑀 𝐶 = 𝑄 2 ideal (𝜇 𝑀 ), see [START_REF] Baudin | Open TURNS: An industrial software for uncertainty quantification in simulation[END_REF], where 𝜇 𝑀 is the empirical measure for a large Monte-Carlo sample (𝑀 = 10 6 ), that will serve as reference when assessing the performance of each of the other estimators. We consider the validation designs built by FSSF, support points and kernel herding, presented in Sections 3.1, 3.2, and 3.3, respectively, and, for each one, we compare the performances obtained for both the uniform and the weighted estimator of Section 2.2.

Test-cases

The training design X 𝑚 and the set S of potential test set points are as in Section 3.4. For test-cases 1 and 3, 𝜇 is the uniform measure on X = [0, 1] 𝑑 , with 𝑑 = 2 and 𝑑 = 8, respectively; X 𝑚 is a maximin Latin hypercube design in X, and S corresponds to the first 𝑁 points S 𝑁 of Sobol' sequence in X, complemented by the 2 𝑑 vertices. In the second test-case, 𝑑 = 2, 𝜇 is the normal distribution 𝒩(0, I 2 ), and the sets X 𝑚 and S 𝑁 must be transformed as explained in section 3.1. There are 𝑁 = 2 14 candidate points for test-cases 1 and 2 and 𝑁 = 2 15 for test-case 3 (this value is rather moderate for a problem in dimension 8, but using a larger 𝑁 yields numerical difficulties for support points; see Section 3.2).

For each test-case, a GP regression model is fitted to the 𝑚 observations using ordinary Kriging [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] (a GP model with constant mean), with an anisotropic Matérn kernel with regularity parameter 5/2: we substitute [(x-x ) D(x-x )] 1/2 for |𝑥 -𝑥 | in [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF], with D a diagonal matrix with diagonal elements 1/𝜃 2 𝑖 , and the correlation lengths 𝜃 𝑖 are estimated by maximum likelihood via a truncated Newton algorithm. All calculations were done using the Python package OpenTURNS for uncertainty quantification [START_REF] Baudin | Open TURNS: An industrial software for uncertainty quantification in simulation[END_REF]. The kernel used for kernel herding is different and corresponds to the tensor product of one-dimensional Matérn kernels [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF], so that the potentials 𝑃 𝐾 , 𝜇 (•) are known explicitly (see Appendix B); the correlations lengths are set to 𝜃 = 0.2 in test-cases 1 and 3 (𝑑 = 2) and to 𝜃 = 0.7 in test-case 3 (𝑑 = 8).

Assuming that a model is classified, in terms of the estimated value of its predictivity index 𝑄 2 as "poor fitting" if 𝑄 2 ∈ [0.6, 0.8], "reasonably good fitting", when 𝑄 2 ∈ (0.8, 0.9], and "very good fitting" if 𝑄 2 > 0.9, we selected, for each test-case three different sizes 𝑚 of the training set such that the corresponding models cover all three possible situations. For all test-cases, the impact of the size 𝑛 of the test set is studied in the range 𝑛 ∈ {4, . . . , 50}.

Test-case 1.

This test function is 𝑓 1 (x) = ℎ(2 𝑥 1 -1, 2 𝑥 2 -1), (𝑥 1 , 𝑥 2 ) ∈ X = [0, 1] 2 , with ℎ(𝑢 1 , 𝑢 2 ) = exp(𝑢 1 ) 5 - 𝑢 2 5 + 𝑢 6 2 3 + 4𝑢 4 2 -4𝑢 2 2 + 7𝑢 2 1 10 + 𝑢 4 1 + 3 4𝑢 2 1 + 4𝑢 2 2 + 1 .
Color coded 3d and contour plots of 𝑓 1 for X ∈ X are shown on the left panel of Figure 3, showing that the function is rather smooth, even if its behaviour along the boundaries of X, in particular close to the vertices, may present difficulties for some regression methods. The size of the training set for this function are: 𝑚 ∈ {5, 15, 30}. Test-case 2.

The second test function, plotted in the right panel of Figure 3 for

x ∈ [0, 1] 2 , is 𝑓 2 (x) = cos 5 + 3 2 𝑥 1 + sin 5 + 3 2 𝑥 1 + 1 100 5 + 3 2 𝑥 1 5 + 3 2 𝑥 2 .
Training set sizes for this test-case are 𝑚 ∈ {8, 15, 30}.

Test-case 3.

The third function is the so-called "gSobol" function, defined over X = [0, 1] 8 by

𝑓 3 (x) = 8 𝑖=1 |4𝑥 𝑖 -2| + 𝑎 𝑖 1 + 𝑎 𝑖 , 𝑎 𝑖 = 𝑖 2 .
This parametric function is very versatile as both the dimension of its input space and the coefficients 𝑎 𝑖 can be freely chosen. The sensitivity to input variables is determined by the 𝑎 𝑖 : the larger 𝑎 𝑖 is, the less 𝑓 is sensitive to 𝑥 𝑖 . Larger training sets are considered for this test-case: 𝑚 ∈ {15, 30, 100}.

Results and analysis

The numerical results obtained in this section are presented in Figures 4,5, and 6. Each figure corresponds to one of the test-cases and gathers three sub-figures, corresponding to test sets with sizes 𝑚 yielding poor (left), reasonably good (centre) or very good (right) fittings.

The baseline value of 𝑄 2 𝑀 𝐶 , calculated with 10 6 Monte-Carlo points, is indicated by the black diamonds (the black horizontal lines). We assume that the error of 𝑄 2 𝑀 𝐶 is much smaller than the errors of all other estimators, and compare the distinct methods through their ability to approximate 𝑄 2 𝑀 𝐶 . For each sequence of nested test-sets (𝑛 ∈ {4, . . . , 50}), the observed values of 𝑄 2 𝑛 (equation ( 2)) and 𝑄 2 𝑛 *

(equation ( 7)), are plotted as the solid and dashed lines, respectively. The figures also show the value 𝑄 2 𝐿𝑂𝑂 obtained by Leave-One-Out (LOO) cross validation, which is indicated at the left of each figure by a red diamond (values smaller than 0.25 are not shown). Note that, contrarily to the other methods considered, for LOO the test set is not disjoint from the training set, and thus the method does not satisfy the conditions set in the Introduction. As we repeat the complete model-fitting procedure for each training sample of size 𝑚 -1, including the maximum-likelihood estimation of the correlation lengths of the Matérn kernel, the closed-form expressions of [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF] cannot be used, making the computations rather intensive. As the three figures show, and as we should expect, 𝑄 2 𝐿𝑂𝑂 tends to underestimate 𝑄 2 ideal : by construction of the training set, LOO cross validation relies on model predictions at points x (𝑖) far from the other 𝑚 -1 design points used to build the model, and thus tends to systematically overestimate the prediction error at x (𝑖) . The underestimation of 𝑄 2 ideal can be particularly severe when 𝑚 is small, the training set being then necessarily sparse; see Figure 4 where 𝑄 2 𝐿𝑂𝑂 < 0.3 for 𝑚 = 5 and 15. Let us first concentrate on the non-weighted estimators (solid curves). We can see that the two MMD-based constructions, support points (in orange) and kernel herding (in blue), generally produce better validation designs than FSSF (green curves), leading to values of 𝑄 2 𝑛 that approach 𝑄 2 ideal quicker as 𝑛 increases. This is particularly noticeable for "good" and "very good" models (central and rightmost panels of all three figures). This supports the idea that test sets should complement the training set X 𝑚 by populating the holes it leaves in X while at the same time be able to mimic the target distribution 𝜇, this second objective being more difficult to achieve for FSSF than for the MMD-based constructions.

Comparison of the two MMD based estimators reveals that support points tend to under-estimate ISE, leading to an over-confident assessment of the model predictivity, while kernel herding displays the expected behaviour, with a negative bias that decreases with 𝑛. The reason for the positive bias of estimates based on support points designs is not fully understood, but may be linked to the fact that support points tend to place validation points at "mid-range" from the designs (and not at the furthest points like FSSF or kernel herding), see central and rightmost panels in Figure 1, and thus residuals at these points are themselves already better representatives of the local average errors.

We consider now the impact of the GP-based weighting of the residuals when estimating 𝑄 2 (by 𝑄 2 𝑛 * ), which is related to the relative training-set/validation-set geometry (the manner in which the two designs are entangled in ambient space). The improvement resulting of applying residual weighting is apparent on all panels of the three figures, the dashed curves lying closer to 𝑄 2 ideal than their solid counterparts; see in particular kernel herding (blue curve) in Figure 4 and FSSF (green curve) in Figure 5. Unexpectedly, the estimators based on support points seem to be rather insensitive to residual weighting, the dashed and solid orange curves being most of the time close to each other (and in any case, much closer that the green and blue ones). While the reason for this behavior deserves a deeper study, the fact that the support point designs -see Figure 1 -sample in a better manner the range of possible training-to-validation distances, being in some sense less space-filling than both FSSF and kernel herding, is again a plausible explanation for this weaker sensitivity to residual weighting.

Consider now comparison of the behaviour across test-cases. Setting aside the strikingly singular situation of test-case 2, for which kernel herding displays a pathological (bad) behaviour for the "very good" model, and all methods present an overall astonishing good behaviour, we can conclude that the details of the tested function do not seem to play an important role concerning the relative merits of the estimators and validation designs.

We finally observe how the methods behave for models of distinct quality (𝑚 leading to poor, good or very good models), comparing the three panels in each figure. On the left panels, 𝑚 is too small for the model 𝜂 𝑚 to be accurate, and all methods and test-set sizes are able to detect this. For models of practical interest (good and very good), the test sets generated with support points and kernel herding allow a reasonably accurate estimation of 𝑄 2 with a few points. Note, incidentally, that except for test-case 2 (where the interplay with a non-uniform measure 𝜇 complicates the analysis), it is in general easier to estimate the quality of the very good model (right-most panel) than that of the good model (central panel), indicating that the expected complexity (the entropy) of the residual process should be a key factor determining how large the validation set must be. In particular, it may be that larger values of 𝑚 allow for smaller values of 𝑛.

Numerical results II: splitting a dataset into a training set and a test set

In this section, we illustrate the performance of the different designs and estimators considered in this paper when applied in the context of an industrial application, to split a given dataset of size 𝑁 into training and test sets, with 𝑚 and 𝑛 points respectively, 𝑚 + 𝑛 = 𝑁. In contrast with [START_REF] Joseph | SPlit: An optimal method for data splitting[END_REF], the observations 𝑦(x (𝑖) ), 𝑖 = 1, . . . , 𝑁, are not used in the splitting mechanism, meaning that it can be performed before the observations are collected and that there cannot be any selection bias related to observations (indeed, the use of observation values in a MMD-based splitting criterion may favour the allocation of the most different observations to different sets, training versus validation).

A ML model is fitted to the training data, and the data collected on the testset are used to assess the predictivity of the model. The influence of the ratio 𝑟 𝑛 = 𝑛/𝑁 = 1 -𝑚/𝑁 on the quality assessment is investigated. We also consider Random Cross-Validation (RCV), where 𝑛 points are chosen at random among the 𝑁 points of the dataset: for each 𝑛, there are 𝑁 𝑛 possible choices, and we randomly select 𝑅 = 1 000 designs among them. We fit a model to each of the 𝑚-point complementary designs (𝑚 = 𝑁 -𝑛), which yields an empirical distribution of 𝑄 2 values for each ratio 𝑛/𝑁 considered.

Industrial test-case CATHARE

The test-case corresponds to the computer code CATHARE2 (for "Code Avancé de ThermoHydraulique pour les Accidents de Réacteurs à Eau"), which models the thermal-hydraulic behavior inside nuclear pressurized water reactors [START_REF] Geffraye | CATHARE2 V2.5_2: A single version for various applications[END_REF]. The studied scenario simulates a hypothetical large-break loss of primary coolant accident for which the output of interest is the peak cladding temperature [START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF][START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF]. The complexity of this application lies in the large run-time of the computer model (of the order of twenty minutes) and in the high dimension of the input space: the model involves 53 input parameters 𝑧 𝑖 , corresponding mostly to constants of physical laws, but also coding initial conditions, material properties and geometrical modeling. The 𝑧 𝑖 were independently sampled according to normal or log-normal distributions (see axes histograms in Figure 7 corresponding to 10 inputs). These characteristics make this test-case challenging in terms of construction of a surrogate model and validation of its predictivity.

We have access to an existing Monte Carlo sample Z 𝑁 of 𝑁 = 1 000 points in R 53 , that corresponds to 53 independent random input configurations; see [START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF] for details. The output of the CATHARE2 code at these 𝑁 points is also available. To reduce the dimensionality of this dataset, we first performed a sensitivity analysis [START_REF] Veiga | Basics and Trends in Sensitivity Analysis[END_REF] to eliminate inputs that do not impact the output significantly. This dimensionreduction step relies on the Hilbert-Schmidt Independence Criterion (HSIC), which is known as a powerful tool to perform input screening from a single sample of inputs and output values without reference to any specific ML regression model [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF][START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF]. HSIC-based statistical tests and their associated 𝑝-values are used to identify (with a 5%-threshold) inputs on which the output is significantly dependent (and therefore, also those of little influence). They were successfully applied to similar datasets from thermal-hydraulic applications in [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: Application on safety studies for nuclear reactor[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]. The screened dataset only includes 10 influential inputs, over which the candidate set X 𝑁 used for the construction of the test-set X 𝑛 (and therefore of the complementary training set X 𝑁 -𝑛 ) is defined. An input-output scatter plot is presented in Figure 7, showing that indeed the retained factors are correlated with the code output. The marginal distributions are shown as histograms along to the axes of the plots.

To include RCV in the methods to be compared, we need to be able to construct many (here, 𝑅 = 1 000) different models 𝜂 𝑚 for each considered design size 𝑚. Since Gaussian Process regression proved to be too expensive for this purpose, we settled for the comparatively cheaper Partial Least Squares (PLS) method [START_REF] Wold | PLS-regression: a basic tool of chemometrics[END_REF], which retains acceptable accuracy. For each given training set, the model obtained is a sum of monomials in the 10 input variables. Note that models constructed with different training sets may involve different monomials and have different numbers of monomial terms. Consider RCV first. For each value of 𝑟 𝑛 = 𝑛/𝑁, the empirical distribution of 𝑄 2 𝑅𝐶𝑉 obtained from 𝑅 = 10 3 random splittings of X 𝑁 into X 𝑚 ∪ X 𝑛 is summarized by a boxplot. Depending on 𝑟 𝑛 , we can roughly distinguish three behaviors. For 0.1 ≤ 𝑟 𝑛 0.3 the distribution is bi-modal, with the lower mode corresponding to unlucky test-set selections leading to poor performance evaluations. When 0.3 𝑛/𝑁 0.7, the distribution looks uni-modal, revealing a more stable performance evaluation. Note that this is (partly) in line with the recommendations discussed in section 1. For 𝑟 𝑛 0.7, the variance of the distribution increases with 𝑟 𝑛 : many unlucky training sets lead to poor models. Note that the median of the empirical distribution slowly decreases as 𝑟 𝑛 increases, which is consistent with the intuition that the model predictivity should decrease when the size of the training set decreases.

Benchmark results and analysis

For completeness, we also show by a red diamond on the left of Figure 8 the value of 𝑄 2 𝐿𝑂𝑂 computed by LOO cross-validation. In principle, being computed using the entire dataset, this value should establish an upper bound on the quality of models computed with smaller training sets. This is indeed the case for small training sets (rightmost values in the figure), for which the predictivity estimated by LOO is above the majority of the predictivity indexes calculated. But at the same time, we know that LOO cross-validation tends to overestimate the errors, which explains the higher predictivity estimated by some other methods when 𝑚 = 𝑁 -𝑛 is large enough.

Compare now the behavior of the two MMD-based algorithms of Section 3, 𝑄 2 𝑛 (un-weighted) and 𝑄 2 𝑛 * (weighted) are plotted using solid and dashed lines, respectively, for both kernel herding (in blue) and support points (in orange). FSSF test-sets are not considered, as the application of an iso-probabilistic transformation imposes knowledge of the input distribution, which is not known for this example. Compare first the unweighted versions of the two MMD-based estimators. For small values of the ratio 𝑟 𝑛 , 0.1 𝑟 𝑛 0.45, the relative behavior of support points and kernel herding coincides with what we observed in the previous section, support points (solid orange line) estimating a better performance than kernel herding (solid blue line), which, moreover, is close to the median of the empirical distribution of 𝑄 2 𝑅𝐶𝑉 . However, for 𝑟 𝑛 ≥ 0.5, the dominance is reversed, support points estimating a worse performance than kernel herding.

As 𝑟 𝑛 increases up to 𝑟 𝑛 0.7 the solid orange and blue curves crossover, and it is now 𝑄 2 𝑛 for kernel herding that approximates the RCV empirical median, while the value obtained with support points underestimates the predictivity index. Also, note that for (irrealistic) very large values of 𝑟 𝑛 both support points and kernel herding estimate lower 𝑄 2 values, which are smaller than the median of the RCV estimates.

Let us now focus on the effect of residual weighting, i.e., in estimators 𝑄 2 𝑛 * which use the weights computed by the method of Section 2.2, shown in dashed lines in Figure 8. First, note that while for kernel herding weighting leads, as in the previous section, to higher estimates of the predictivity (compare solid and dashed blue lines), this is not the case for support points (solid and dashed orange curves), which, for small split ratios, produces smaller estimates when weighting is introduced. In the large 𝑟 𝑛 region, the behavior is consistent with what we saw previously, weighting inducing an increase of the estimated predictivity. It is remarkable -and rather surprising -that 𝑄 2 𝑛 * for support points (the dashed orange line) does not present the discontinuity of the uncorrected curve.

The sum 𝑛 𝑖=1 𝑤 * 𝑖 of the optimal weights of support points and kernel herding (6) is shown in Figure 9 (orange and blue curves, respectively). The slow increase with 𝑛/𝑁 of the sum of kernel-herding weights (blue line) is consistent with the increase of the volume of the input region around each validation point when the size of the training set decreases. The behavior of the sum of weights is more difficult to interpret for support points (orange line) but is consistent with the behavior of 𝑄 2 𝑛 * on Figure 8. Note that the energy-distance kernel (10) used for support points cannot be used for the weighting method of Section 2.2 as 𝐾 𝐸 is not positive definite but only conditionally positive definite. A full understanding of the observed curves would require a deeper analysis of the geometric characteristics of the designs generated by the two MMD methods, in particular of their interleaving with the training designs, which is not compatible with the space constraints of this manuscript.

While a number of unanswered points remain, in particular how deeply the behaviours observed may be affected by the poor predictivity resulting from the chosen PLS modeling methodology, the example presented in this section shows that the construction of test sets via MMD minimization and estimation of the predictivity index using the weighted estimator 𝑄 2 𝑛 * is promising as an efficient alternative to RCV: at a much lower computational cost, it builds performance estimates based on independent data the model developers may not have access to. Moreover, kernel herding proved, in the examples studied in this manuscript, to be a more reliable option for designing the test set, exhibiting a behavior that is consistent with what is expected, and very good estimation quality when the residuals over the design points are appropriately weighted.

Conclusion

Our study shows that ideas and tools from the design of experiments framework can be transposed to the problem of test-set selection. This paper explored approaches based on support points, kernel herding and FSSF, considering the incremental construction of a test set (i) either as a particular space-filling design problem, where design points should populate the holes left in the design space by the training set, or (i) from the point of view of partitioning a given dataset into a training set and a test set. A numerical benchmark has been performed for a panel of test-cases of different dimensions and complexity. Additionally to the usual predictivity coefficient, a new weighted metric (see [START_REF] Pronzato | Validation design I: construction of validation designs via kernel herding[END_REF]) has been proposed and shown to improve assessment of the predictivity of a given model for a given test set.

This weighting procedure appears very efficient for interpolators, like Gaussian process regression models, as it corrects the bias when the points in the test set used to predict the errors are far from the training points. For the first three testcases (Section 4), pairing one iterative design method with the weight-corrected estimator of the predictivity coefficient 𝑄 2 shows promising results as the estimated 𝑄 2 characteristic is close to the true one even for test-sets of moderate size.

Weighting can also be applied to models that do not interpolate the training data. For the industrial test-case of Section 5, the true 𝑄 2 value is unknown, but the weightcorrected estimation 𝑄 2 𝑛 * of 𝑄 2 is close to the value estimated by Leave-One-Out cross validation and to the median of the empirical distribution of 𝑄 2 values obtained by random 𝑘-fold cross-validation. At the same time, estimation by 𝑄 2 𝑛 * involves a much smaller computational cost than cross-validation methods, and uses a dataset fully independent from the one used to construct the model.

To each of the design methods considered to select a test set a downside can be attached. FSSF requires knowledge of the input distribution to be able to apply an iso-probabilistic transformation if necessary; it tends to select many points along the boundary of the candidate set considered. Support points require the computation of the 𝑁 (𝑁 -1)/2 distances between all pairs of candidate points, which implies important memory requirements for large 𝑁; the energy-distance kernel on which the method relies cannot be used for the weighting procedure. Finally, the efficient implementation of kernel herding relies on analytical expressions for the potentials 𝑃 𝐾 , 𝜇 , see Appendices A and B, which are available for particular distributions (like the uniform and the normal) and kernels (like Matérn) only. The great freedom in the choice of the kernel 𝐾 gives a lot of flexibility, but at the same time implies that some non-trivial decisions have to be made; also, the internal parameters of 𝐾, such as its correlation lengths, must to be specified. Future work should go beyond empirical rules of thumb and study the influence of these choices.

We have only computed numerical tests with independent inputs. Kernel herding and support points are both well suited for probability measures not being equal to the product of their marginals, which is a frequent case in real datasets. We have also only considered incremental constructions, as they allow to stop the validation procedure as soon as the estimation of the model predictivity is deemed sufficiently accurate, but it is also possible to select several points at once, using support points [START_REF] Mak | Support points[END_REF], or MMD minimization in general [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF].

Further developments around this work could be as follows. Firstly, the incremental construction of a test set could be coupled with the definition of an appropriate stopping rule, in order to decide when it is necessary to continue improving the model (possibly by supplementing the initial design with the test set, which seems well suited to this). The MMD 𝑑 𝐾 |𝑚 (𝜁 * 𝑛 , 𝜇) of Section 2.2 could play an important role in the derivation of such a rule. Secondly, the approach presented gives equal importance to all the 𝑑 inputs. However, it seems that inputs with a negligible influence on the output should receive less attention when selecting a test set. A preliminary screening step that identifies the important inputs would allow the test-set selection algorithm to be applied on these variables only. For example, when a X 𝑁 ⊂ R 𝑑 dataset is to be partitioned into X 𝑚 ∪ X 𝑛 , one could use only 𝑑 < 𝑑 components to define the partition, but still use all 𝑑 components to build the model and estimate its (weighted) 𝑄 2 . Note, however, that this would imply a slight violation of the conditions mentioned in introduction, as it renders the test set dependent on the function observations.

Finally, in some cases the probability measure 𝜇 is known up to a normalizing constant. The use of a Stein kernel then makes the potential 𝑃 𝐾 ,𝜇 identically zero [START_REF] Chen | Stein points[END_REF][START_REF] Chen | Stein point Markov Chain Monte Carlo[END_REF], which would facilitate the application of kernel herding. Also, more complex problems involve functional inputs, like temporal signals or images, or categorical variables; the application of the methods presented to kernels specifically designed for such situations raises challenging issues. 𝜃 𝑥 5

√ 5𝑥 2 + 25𝜃𝑥 + 8 √ 5𝜃 2 .
The expressions 𝑃 𝐾 𝜈, 𝜃 , 𝜇 (𝑥) for 𝜈 = 1/2 and 𝜈 = 3/2 can be found in [START_REF] Pronzato | Bayesian quadrature and energy minimization for spacefilling design[END_REF].

When 𝜇 is the standard normal distribution N (0, 1), the potential 𝑃 𝐾 5/2, 𝜃 , N (0,1) is 𝑃 𝐾 5/2, 𝜃 ,N (0,1) (𝑥) = 𝑇 𝜃 (𝑥) + 𝑇 𝜃 (-𝑥), where 𝑇 𝜃 (𝑥) = .

Fig. 1 Fig. 2

 12 Fig.[START_REF] Baudin | Open TURNS: An industrial software for uncertainty quantification in simulation[END_REF] Additional points (ordered, green) complementing an initial design (red crosses), 𝜇 is uniform on [0, 1], the candidate points are in gray.

2 Fig. 3

 23 Fig.3Left: 𝑓 1 (x) (test-case 1); right:𝑓 2 (x) (test-case 2); x ∈ X = [0, 1] 2 .

Fig. 4 Fig. 5

 45 Fig.[START_REF] Chen | Stein point Markov Chain Monte Carlo[END_REF] Test-case 1: predictivity assessment of a poor (left), good (center) and very good (right) model with kernel herding, support points and FSSF test sets.

Fig. 6

 6 Fig. 6Test-case 3: predictivity assessment of a poor (left), good (center) and very good (right) model with kernel herding, support points and FSSF test sets.

YFig. 7

 7 Fig. 7 Test-case CATHARE: inputs output scatter plots (𝑁 = 10 3 )

Figure 8

 8 Figure 8 compares various ways of extracting an 𝑛-point test set from an 𝑁point dataset to estimate model predictivity, for different splitting ratios 𝑛/𝑁 ∈ {0.1, 0.15, 0.2, . . . , 0.9}.Consider RCV first. For each value of 𝑟 𝑛 = 𝑛/𝑁, the empirical distribution of 𝑄2 𝑅𝐶𝑉 obtained from 𝑅 = 10 3 random splittings of X 𝑁 into X 𝑚 ∪ X 𝑛 is summarized by a boxplot. Depending on 𝑟 𝑛 , we can roughly distinguish three behaviors. For 0.1 ≤ 𝑟 𝑛 0.3 the distribution is bi-modal, with the lower mode corresponding to unlucky test-set selections leading to poor performance evaluations. When 0.3 𝑛/𝑁 0.7, the distribution looks uni-modal, revealing a more stable performance evaluation. Note that this is (partly) in line with the recommendations discussed in section 1. For 𝑟 𝑛 0.7, the variance of the distribution increases with 𝑟 𝑛 : many unlucky training sets lead to poor models. Note that the median of the empirical distribution slowly decreases as 𝑟 𝑛 increases, which is consistent with the intuition that the model predictivity should decrease when the size of the training set decreases.For completeness, we also show by a red diamond on the left of Figure8the value of 𝑄 2 𝐿𝑂𝑂 computed by LOO cross-validation. In principle, being computed using the entire dataset, this value should establish an upper bound on the quality

Fig. 8

 8 Fig.8Test-case CATHARE: estimated 𝑄 2 . The box-plots are for random cross-validation, the red diamond (left) is for 𝑄 2 𝐿𝑂𝑂 .

Fig. 9

 9 Fig.9Test-case CATHARE: sum of the weights[START_REF] Chen | Super-samples from kernel herding[END_REF].

  (𝑚) )] . The model above allows us to study how well ISE 𝜇 (X 𝑚 , y 𝑚 ) is estimated using a given test set. Denote by Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚 ) the expected squared error when estimating ISE 𝜇 (X 𝑚 , y 𝑚 ) by ISE 𝜉 𝑛 (X 𝑚 , y 𝑚 ), Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚

  𝜉 𝑛 , uniform on X 𝑛 , by a nonuniform measure 𝜁 𝑛 supported on X 𝑛 , 𝜁 𝑛 = 𝑛 𝑖=1 𝑤 𝑖 𝛿 x (𝑚+𝑖) with weights w 𝑛 = (𝑤 1 , . . . , 𝑤 𝑛 ) chosen such that the estimation error Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚 ), and thus 𝑑 2 𝐾 |𝑚 , 𝜇 (X 𝑛 ) = 𝑃 𝐾 |𝑚 ,𝜇 (x (𝑚+1) ), . . . , 𝑃 𝐾 |𝑚 ,𝜇 (x (𝑚+𝑛) ) , with 𝑃 𝐾 |𝑚 ,𝜇 (x) defined by[START_REF] Goodfellow | Deep learning[END_REF] in Appendix A, {K |𝑚 (X 𝑛 )} 𝑖, 𝑗 = 𝐾 |𝑚 (x (𝑚+𝑖) , x (𝑚+ 𝑗) ) for any 𝑚 ∩ X 𝑛 = ∅, the 𝑛 × 𝑛 matrix K |𝑚 (X 𝑛 ), whose element 𝑖, 𝑗 equals 𝐾 |𝑚 (x (𝑚+𝑖) , x (𝑚+ 𝑗) ), is positive definite. The elementwise (Hadamard) product K |𝑚 (X 𝑛 ) • K |𝑚 (X 𝑛 ) is thus positive definite too, implying that K |𝑚 (X 𝑛 ) is positive definite. The optimal weights w * 𝑛 minimizing 𝑑 2

	where p
	𝑖, 𝑗 = 1, . . . , 𝑛, and ℰ 𝐾 |𝑚 (𝜇) =	∫	X 2 𝐾 |𝑚 (x, x ) d𝜇(x)d𝜇(x ).
	When X 𝐾 |𝑚	(𝜁 𝑛 , 𝜇) are thus
			2	(𝜉 𝑛 , 𝜇) as the squared Maximum-Mean-Discrepancy (MMD)
			𝐾 |𝑚
	between 𝜉 𝑛 and 𝜇 for the kernel 𝐾 |𝑚 ; see (18) in Appendix A. Note that 𝜎 2 only
	appears as a multiplying factor in (4), with the consequence that 𝜎 2 does not impact
	the choice of a suitable 𝜉 𝑛 .
	When 𝜂 𝑚 does not interpolate y 𝑚 , similar developments still give Δ 2 (𝜉 𝑛 , 𝜇; X 𝑚 , y 𝑚 ) =
	𝜎 2 𝑑 2	(𝜉 𝑛 , 𝜇), with now
	𝐾 |𝑚
		𝐾 |𝑚 (x, x ) = 2 𝐾 |𝑚 (x, x ) + 2 𝛿 𝑚 (x) 𝛿 𝑚 (x ) 𝐾 |𝑚 (x, x )
				+ 𝛿 2 𝑚 (x) + 𝐾 |𝑚 (x, x) 𝛿 2 𝑚 (x ) + 𝐾 |𝑚 (x , x ) ,
	where 𝛿 𝑚 (x) is given by (3).
	The idea is to replace 𝐾 |𝑚	(𝜁 𝑛 , 𝜇), is minimized. Direct
	calculation gives
		𝑑 2 𝐾 |𝑚	(𝜁 𝑛 , 𝜇) = ℰ 𝐾 |𝑚 (𝜇) -2 w 𝑛 p	𝐾 |𝑚 ,𝜇 (X 𝑛 ) + w 𝑛 K |𝑚 (X 𝑛 )w 𝑛 ,

  is the product of its marginals 𝜇 [𝑖 ] on the X [𝑖 ] , 𝐾 is the product of one-dimensional kernels 𝐾 [𝑖 ] , and the one-dimensional integral in 𝑃 𝐾 [𝑖 ] , 𝜇 [𝑖 ] (𝑥) is known explicitly for each 𝑖 ∈ {1, . . . , 𝑑}. Indeed, for x = (𝑥 1 , . . . , 𝑥 𝑑 ) ∈ X, we then have 𝑃 𝐾 ,𝜇 (x) = 𝑑 𝑖=1 𝑃 𝐾 [𝑖 ] ,𝜇 [𝑖 ] (𝑥 𝑖 ); see [40]. When 𝐾 is the product of Matérn kernels with regularity parameter 5/2 and correlation lengths 𝜃 𝑖 , 𝐾 (x, x ) = 𝑑 𝑖=1 𝐾 5/2, 𝜃 𝑖 (𝑥 𝑖 -𝑥

𝑖 ), with

Acknowledgements

This work was supported by project INDEX (INcremental Design of EXperiments) ANR-18-CE91-0007 of the French National Research Agency (ANR). The authors are grateful to Guillaume Levillain and Thomas Bittar for their code development during their work at EDF. Thanks also to Sébastien Da Veiga for fruitful discussions.

Appendix Appendix A: Maximum Mean Discrepancy

Let 𝐾 be a positive definite kernel on X × X, defining a reproducing kernel Hilbert space (RKHS) H 𝐾 of functions on X, with scalar product 𝑓 , 𝑔 H 𝐾 and norm 𝑓 H 𝐾 ; see, e.g., [START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF]. For any 𝑓 ∈ H 𝐾 and any probability measures 𝜇 and 𝜉 on X, we have

where we have denoted 𝐾 x (•) = 𝐾 (x, •) and used the reproducing property 𝑓 (x) = 𝑓 , 𝐾 x H 𝐾 for all x ∈ X, and where, for any probability measure 𝜈 on X and x ∈ X,

is the potential of 𝜈 at x. 𝑃 𝐾 ,𝜈 ∈ H 𝐾 and is called kernel embedding of 𝜈 in ML. In some cases, the potential can be expressed analytically (see. Appendix 6), otherwise it can be estimated by numerical quadrature (Quasi Monte Carlo). Cauchy-Schwartz inequality applied to [START_REF] Geffraye | CATHARE2 V2.5_2: A single version for various applications[END_REF] gives

The Maximum Mean Discrepancy (MMD) between 𝜉 and 𝜇 (for the kernel 𝐾 and set X) is 𝑑 𝐾 (𝜉, 𝜇) = 𝑃 𝐾 , 𝜉 -𝑃 𝐾 ,𝜇 H 𝐾 . Direct calculation gives

where the random variables 𝜁 and 𝜁 in ( 19) are independent, see [START_REF] Smola | A Hilbert space embedding for distributions[END_REF]. When 𝐾 is the energy distance kernel [START_REF] Veiga | Basics and Trends in Sensitivity Analysis[END_REF], one recovers the expression [START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF] for the corresponding MMD. One may refer to [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] for an illuminating exposition on MMD, kernel embedding, and conditions on 𝐾 (the notion of characteristic kernel) that make 𝑑 𝐾 a metric on the space of probability measures on X. The distance and Matérn kernels considered in this paper are characteristic.

Appendix B: Analytical computation of potentials for Matérn kernels

As for tensor-product kernels, the potential is the product of the one-dimensional potentials, we only consider one-dimensional input spaces.

For 𝜇 the uniform distribution on [0, 1] and 𝐾 the Matérn kernel 𝐾 5/2, 𝜃 with smoothness 𝜈 = 5/2 and correlation length 𝜃, see [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF], we get 𝑃 𝐾 5/2, 𝜃 , 𝜇 (𝑥) =