Loïc Le Gratiet

Stefano Marelli

Bruno Sudret

Metamodel-based sensitivity analysis: Polynomial chaos expansions and Gaussian processes

Keywords: Polynomial Chaos Expansions, Gaussian Processes, Kriging, Error estimation, Sobol' indices

Global sensitivity analysis is now established as a powerful approach for determining the key random input parameters that drive the uncertainty of model output predictions. Yet the classical computation of the so-called Sobol' indices is based on Monte Carlo simulation, which is not affordable when computationally expensive models are used, as it is the case in most applications in engineering and applied sciences. In this respect metamodels such as polynomial chaos expansions (PCE) and Gaussian processes (GP) have received tremendous attention in the last few years, as they allow one to replace the original, taxing model by a surrogate which is built from an experimental design of limited size. Then the surrogate can be used to compute the sensitivity indices in negligible time. In this chapter an introduction to each technique is given, with an emphasis on their strengths and limitations in the context of global sensitivity analysis. In particular, Sobol' (resp. total Sobol') indices can be computed analytically from the PCE coefficients. In contrast, confidence intervals on sensitivity indices can be derived straightforwardly from the properties of GPs. The performance of the two techniques is finally compared on three well-known analytical benchmarks (Ishigami, G-Sobol and Morris functions) as well as on a realistic engineering application (deflection of a truss structure).

Introduction

In modern engineering sciences computational models are used to simulate and predict the behavior of complex systems. The governing equations of the system are usually discretized so as to be solved by dedicated algorithms. In the end a computational model (a.k.a. simulator) is built up, which can be considered as a mapping from the space of input parameters to the space of quantities of interest that are computed by the model. However, in many situations the values of the parameters describing the properties of the system, its environment and the various initial and boundary conditions are not perfectly well-known. To account for such uncertainty, they are typically described by possible variation ranges or probability distribution functions.

In this context global sensitivity analysis aims at determining which input parameters of the model influence the most the predictions, i.e. how the variability of the model response is affected by the uncertainty of the various input parameters. A popular technique is based on the decomposition of the response variance as a sum of contributions that can be associated to each single input parameter or to combinations thereof, leading to the computation of the so-called Sobol' indices.

As presented earlier in this book (see Variance-based sensitivity analysis: Theory and estimation algorithms), the use of Monte Carlo simulation to compute Sobol' indices requires a large number of samples (typically, thousands to hundreds of thousands), which may be an impossible requirement when the underlying computational model is expensive-to-evaluate. To bypass this difficulty, surrogate models may be built. Generally speaking, a surrogate model (a.k.a. metamodel or emulator) is an approximation of the original computational model:

x ∈ D X ⊂ R d → y = G(x) (1)
which is constructed based on a limited number of runs of the true model, the so-called experimental design:

X = x (1) , . . . , x (n) . (2)
Once a type of surrogate model is selected, the parameters have to be fitted based on the information contained in the experimental design X and associated runs of the original computational model Y = y i = G(x (i) ), i = 1, . . . , n . Then the accuracy of the surrogate shall be estimated by some kind of validation technique. For a general introduction to surrogate modelling the reader is referred to [START_REF] Storlie | Implementation and evaluation of nonparametric regression procedures for sensitivity analysis of computationally demanding models[END_REF] and to the recent review by Iooss and Lemaître [START_REF] Iooss | Uncertainty management in simulation-optimization of complex systems: algorithms and applications[END_REF].

In this chapter we discuss two classes of surrogate models that are commonly used for sensitivity analysis, namely polynomial chaos expansions (PCE) and Gaussian processes (GP). The use of polynomial chaos expansions in the context of sensitivity analysis has been originally presented in Sudret [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] using a non intrusive least-square method. Other non-intrusive strategies for the calculation of PCE coefficients include spectral projection through sparse grids (e.g. Crestaux et al. [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]; Buzzard and Xiu [START_REF] Buzzard | Variance-based global sensitivity analysis via sparse-grid interpolation and cubature[END_REF]; Buzzard [START_REF] Buzzard | Global sensitivity analysis using sparse grid interpolation and polynomial chaos[END_REF]) and sparse polynomial expansions (e.g. Blatman and Sudret [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF]). In the last five years numerous application examples have been developed using PCE for sensitivity analysis, e.g. Fajraoui et al. [START_REF] Fajraoui | Use of global sensitivity analysis and polynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media[END_REF];

Younes et al. [START_REF] Younes | Use of global sensitivity analysis to help assess unsaturated soil hydraulic parameters[END_REF]; Brown et al. [START_REF] Brown | Global sensitivity analysis of the impact of impurities on CO2 pipeline failure[END_REF]; Sandoval et al. [START_REF] Sandoval | Sensitivity study of dynamic systems using polynomial chaos[END_REF]. Recent extensions to problems with dependent input parameters can be found in Sudret and Caniou [START_REF] Sudret | Analysis of covariance (ancova) using polynomial chaos expansions[END_REF]; Munoz Zuniga et al. [START_REF] Munoz Zuniga | Metamodelling with independent and dependent inputs[END_REF].

In parallel, Gaussian process modeling has been introduced in the context of sensitivity analysis by Welch et al. [START_REF] Welch | Screening, predicting, and computer experiments[END_REF]; Oakley and O'Hagan [START_REF] Oakley | Probabilistic sensitivity analysis of complex models a Bayesian approach[END_REF]; Marrel et al. [START_REF] Marrel | An efficient methodology for modeling complex computer codes with gaussian processes[END_REF][START_REF] Marrel | Calculations of Sobol indices for the Gaussian process metamodel[END_REF]. Recent developments in which metamodeling errors are taken into account in the analysis have been proposed by Le Gratiet et al. [START_REF] Gratiet | A bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF]; Chastaing and Le Gratiet [START_REF] Chastaing | Anova decomposition of conditional gaussian processes for sensitivity analysis with dependent inputs[END_REF].

The chapter first recalls the basics of the two approaches and details how they can be used to compute sensitivity indices. The two approaches are then compared on different benchmark examples as well as on an application in structural mechanics.

2 Polynomial chaos expansions

Mathematical setup

Let us consider a computational model G : x ∈ D X ⊂ R d → y = G(x) ∈ R. Suppose that the uncertainty in the input parameters is modeled by a random vector X with prescribed joint probability density function (PDF) f X (x). The resulting (random) quantity of interest Y = G(X) is obtained by propagating the uncertainty in X through G. Assuming that Y has a finite variance (which is a physically meaningful assumption when dealing with physical systems), it belongs to the so-called Hilbert space of second order random variables, which allows for the following spectral representation to hold [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]:

Y = ∞ j=0 y j Z j . (3)
The random variable Y is therefore cast as an infinite series, in which {Z j } ∞ j=0 is a numerable set of random variables (which form a basis of the Hilbert space), and {y j } ∞ j=0 are coefficients. The latter may be interpreted as the coordinates of Y in this basis. In the sequel we focus on polynomial chaos expansions, in which the basis terms {Z j } ∞ j=0 are multivariate orthonormal polynomials in the input vector X, i.e. Z j = Ψ j (X).

Polynomial chaos basis

In the sequel we assume that the input variables are statistically independent, so that the joint PDF is the product of the d marginal distributions: f X (x) = d i=1 f Xi (x i ), where the f Xi (x i ) are the marginal distributions of each variable {X i , i = 1, . . . , d} defined on D Xi . For each single variable X i and any two functions φ 1 , φ 2 : x ∈ D Xi → R, we define the functional inner product by the following integral (provided it exists):

φ 1 , φ 2 i def = E [φ 1 (X i ) φ 2 (X i )] = D X i φ 1 (x) φ 2 (x) f Xi (x) dx. (4)
Using the above notation, classical algebra allows one to build a family of orthogonal polynomials {P

(i) k , k ∈ N} satisfying P (i) j , P (i) k i def = E P (i) j (X i ) P (i) k (X i ) = D X i P (i) j (x) P (i) k (x) f Xi (x) dx = a (i) j δ jk , (5)
see e.g. Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]. In the above equation subscript k denotes the degree of the polynomial

P (i)
k , δ jk is the Kronecker symbol equal to 1 when j = k and 0 otherwise and a (i) j corresponds to the squared norm of

P (i) j : a (i) j def = P (i) j 2 i = P (i) j , P (i) j i . (6)
In general orthogonal bases may be obtained by applying the Gram-Schmidt orthogonalization procedure, e.g. to the canonical family of monomials 1, x, x 2 , . . . . For standard distributions, the associated families of orthogonal polynomials are well-known [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. For instance, if X i ∼ U(-1, 1) has a uniform distribution over [-1, 1], the resulting family is that of the so-called Legendre polynomials. When X i ∼ N (0, 1) has a standard normal distribution with zero mean value and unit standard deviation, the resulting family is that of the Hermite polynomials. The families associated to standard distributions are summarized in Table 1 (taken from Sudret [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -contributions to structural reliability and stochastic spectral methods[END_REF]).

Table 1: Classical families of orthogonal polynomials (taken from Sudret [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -contributions to structural reliability and stochastic spectral methods[END_REF])

Type of variable Distribution Orthogonal polynomials Hilbertian basis

ψ k (x) Uniform U(-1, 1) 1 [-1,1] (x)/2 Legendre P k (x) P k (x)/ 1 2k+1
Gaussian N (0, 1)

1 √ 2π e -x 2 /2
Hermite H e k (x)

H e k (x)/ √ k! Gamma Γ(a, λ = 1) x a e -x 1 R + (x) Laguerre L a k (x) L a k (x)/ Γ(k+a+1) k! Beta B(a, b) 1 [-1,1] (x) (1-x) a (1+x) b B(a) B(b) Jacobi J a,b k (x) J a,b k (x)/J a,b,k J 2 a,b,k = 2 a+b+1 2k+a+b+1 Γ(k+a+1)Γ(k+b+1) Γ(k+a+b+1)Γ(k+1)
Note that the obtained family is usually not orthonormal. By enforcing normalization, an orthonormal family

ψ (i) j ∞ j=0
is obtained from Eqs.( 5), [START_REF] Berveiller | Presentation of two methods for computing the response coefficients in stochastic finite element analysis[END_REF] as follows (see Table 1):

ψ (i) j = P (i) j / a (i) j i = 1, . . . , d, j ∈ N. (7)
From the sets of univariate orthonormal polynomials one can now build multivariate orthonormal polynomials with a tensor product construction. For this purpose let us define the multi-indices α ∈ N d , which are ordered lists of integers:

α = (α 1 , . . . , α d ) , α i ∈ N. (8)
One can associate a multivariate polynomial Ψ α to any multi-index α by

Ψ α (x) def = d i=1 ψ (i) αi (x i ), (9)
where the univariate polynomials ψ

(i)
k , k ∈ N are defined above, see Eqs.( 5), [START_REF] Berveiller | Stochastic finite elements: a non intrusive approach by regression[END_REF]. By virtue of Eq.( 5) and the above tensor product construction, the multivariate polynomials in the input vector X are also orthonormal, i.e.

E [Ψ α (X) Ψ β (X)] def = D X Ψ α (x)Ψ β (x) f X (x) dx = δ αβ ∀ α, β ∈ N d , (10)
where δ αβ is the Kronecker symbol which is equal to 1 if α = β and zero otherwise. With this notation, it can be proven that the set of all multivariate polynomials in the input random vector X forms a basis of the Hilbert space in which Y = G(X) is to be represented [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]:

Y = α∈N d y α Ψ α (X). ( 11
)

Non standard variables and truncation scheme

In practical sensitivity analysis problems the input variables may not necessarily have standardized distributions as the ones described in Table 1. Thus reduced variables U with standardized distributions are introduced first through an isoprobabilistic transform:

X = T (U ). ( 12
)
For instance, when dealing with independent uniform distributions with support D Xi = [a i , b i ], i = 1, . . . , d, the isoprobabilistic transform reads:

X i = a i + b i 2 + b i -a i 2 U i U i ∼ U([-1, 1]). ( 13
)
In the case of Gaussian independent variables {X i ∼ N (µ i , σ i ) , i = 1, . . . , d}, the one-to-one mapping reads:

X i = µ i + σ i U i , U i ∼ N (0, 1) (14)
In the general case when the input variables are non Gaussian (e.g. Gumbel distributions, see application in Section 4.4), the one-to-one mapping may be obtained as follows:

X i = F -1 Xi (Φ(U i )) U i ∼ N (0, 1) (15)
where F Xi (resp. Φ) is the cumulative distribution function (CDF) of variable X i (resp. the standard normal CDF).

This isoprobabilistic transform approach also allows one to address problems with dependent variables.

For instance, if the input vector X is defined by a set of marginal distributions and a Gaussian copula, it can be transformed into a set of independent standard normal variables using the Nataf transform [20; 35].

The representation of the random response in Eq.( 11) is exact when the infinite series is considered.

However, in practice, only a finite number of terms may be computed. For this purpose a truncation scheme has to be selected. Since the polynomial chaos basis consists of multivariate polynomials, it is natural to consider as a truncated series all the polynomials up to a given maximum degree. Let us define the total degree of a multivariate polynomial Ψ α by:

|α| def = d i=1 α i . (16)
The standard truncation scheme consists in selecting all polynomials such that |α| is smaller than or equal to a given p. This leads to a set of polynomials denoted by

A d,p = α ∈ N d : |α| ≤ p of cardinality: card A d,p = d + p p = (d + p)! d! p! . ( 17
)
The maximal polynomial degree p may typically be equal to 3 -5 in practical applications. Note that the cardinality of A d,p increases exponentially with d and p. Thus the number of terms in the series, i.e. the number of coefficients to be computed, increases dramatically when d is large, say d > 10. This complexity is referred to as the curse of dimensionality. This issue may be solved using specific algorithms to compute sparse PCE, see e.g. Blatman and Sudret [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]; Doostan and Owhadi [START_REF] Doostan | A non-adapted sparse approximation of pdes with stochastic inputs[END_REF].

Computation of the coefficients and error estimation

The use of polynomial chaos expansions has emerged in the late eighties in uncertainty quantification problems under the form of stochastic finite element methods [START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF]. In this setup the constitutive equations of the physical problem are discretized both in the physical space (using standard finite element techniques) and in the random space using polynomial chaos expansion. This results in coupled systems of equations which require ad-hoc solvers, thus the term "intrusive approach".

Non intrusive techniques such as projection or stochastic collocation have emerged in the last decade as a means to compute the coefficients of PC expansions from repeated evaluations of the existing model G considered as a black-box function. In this section we focus on a particular non intrusive approach based on least-square analysis.

Following Berveiller et al. [START_REF] Berveiller | Presentation of two methods for computing the response coefficients in stochastic finite element analysis[END_REF][START_REF] Berveiller | Stochastic finite elements: a non intrusive approach by regression[END_REF], the computation of the PCE coefficients may be cast as a least-square minimization problem (originally termed "regression" problem) as follows: once a truncation scheme

A ⊂ N d is chosen (for instance, A = A d,p
), the infinite series is recast as the sum of the truncated series and a residual:

Y = G(X) = α∈A y α Ψ α (X) + ε, (18)
in which ε corresponds to all those PC polynomials whose index α is not in the truncation set A. The least-square minimization approach consists in finding the set of coefficients y = {y α , α ∈ A} which minimizes the mean square error

E ε 2 def = E   G(X) - α∈A y α Ψ α (X) 2   . ( 19
)
The set of coefficients y is computed at once by solving:

y = arg min y∈R cardA E   G(X) - α∈A y α Ψ α (X) 2   . (20)
In practice the discretized version of the problem is obtained by replacing the expectation operator in Eq.( 20) by the empirical mean over a sample set: ŷ = arg min

y∈R cardA 1 N N i=1 G(x (i) ) - α∈A y α Ψ α (x (i) ) 2 . (21)
In this expression, X = x (i) , i = 1, . . . , n is a sample set of points called experimental design (ED) that is typically obtained by Monte Carlo simulation of the input random vector X. To solve the leastsquare minimization problem in Eq.( 21) the computational model G is first run for each point in the ED, and the results are stored in a vector Y = y (1) = G(x (1) ), . . . , y

(n) = G(x (n) ) T .
Then the so-called information matrix is calculated from the evaluation of the basis polynomials onto each point in the ED:

A = A ij def = Ψ j (x (i) ) , i = 1, . . . , n, j = 1, . . . , card A . ( 22
)
The solution of the least-square minimization problem finally reads:

ŷ = A T A -1 A T Y. (23)
The points used in the experimental design may be obtained from crude Monte Carlo simulation. However other types of designs are of common use, especially Latin Hypercube sampling (LHS), see McKay et al.

[39], or quasi-random sequences such as the Sobol' or Halton sequence [START_REF] Niederreiter | Random number generation and quasi-Monte Carlo methods[END_REF]. The size of the experimental design is of crucial importance: it must be larger than the number of unknowns cardA for the problem to be well-posed. In practice we use the thumb rule n ≈ 2 -3 card A [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF].

The simple least-square approach summarized above does not allow one to cope with the curse of dimensionality. Indeed the standard truncation scheme requires approximately 3 • d+p p runs of the original model G(x), which is in the order of 10 4 when e.g. d ≥ 15, p ≥ 5. However, in practice most of the problems lead eventually to sparse expansions, i.e. PCE in which most of the coefficients are zero or negligible. In order to find directly the significant polynomials and associated coefficients, sparse PCE have been introduced recently by Blatman and Sudret [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF][START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]; Bieri and Schwab [START_REF] Bieri | Sparse high order FEM for elliptic sPDEs[END_REF]. The recent developments make use of specific selection algorithms which, by solving a penalized least-square problem, lead by construction to sparse expansions. Of interest in this chapter is the use of the least-angle regression algorithm (LAR, Efron et al. [START_REF] Efron | Least angle regression[END_REF]), which was introduced in the field of uncertainty quantification by Blatman and Sudret [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]. Details can be found in Sudret [START_REF] Sudret | Polynomial chaos expansions and stochastic finite element methods[END_REF]. Note that other techniques based on compressive sensing have also been developed recently, see e.g. Doostan and Owhadi [START_REF] Doostan | A non-adapted sparse approximation of pdes with stochastic inputs[END_REF]; Sargsyan et al.

[47]; Jakeman et al. [START_REF] Jakeman | Enhancing 1 -minimization estimates of polynomial chaos expansions using basis selection[END_REF].

Error estimation

The truncation of the polynomial chaos expansion introduces an approximation error which may be computed a posteriori. Based on the data contained in the experimental design, the empirical error may be computed from Eq.( 21) once least-square minimization problem has been solved:

ε emp = 1 N N i=1 G(x (i) ) - α∈A ŷα Ψ α (x (i) ) 2 . (24)
However, this estimator usually underestimates severely the mean square error in Eq. [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]. In particular, if the size N of the experimental design is close to the number of unknown coefficients card A, the empirical error tends to zero whereas the true mean square error does not.

A more robust error estimator can be derived based on the cross-validation technique. The experimental design is split into a training set and a validation set: the coefficients of the expansion are computed using the training set (Eq.( 21)) whereas the error is estimated using the validation set. The leave-one-out cross-validation is a particular case in which all points but one are used to compute the coefficients.

Setting aside x (i) ∈ X , a PCE denoted by G PC\i (X) is built up using the experimental design X \x (i) def =

x (1) , . . . , x (i-1) , x (i+1) , . . . , x (n) . Then the error is computed at point x (i) :

∆ i def = G(x (i) ) -G PC\i (x (i) ). ( 25
)
The LOO error is defined by:

ε LOO = 1 n n i=1 ∆ 2 i = 1 n n i=1 G(x (i) ) -G PC\i (x (i) ) 2 . ( 26
)
After some algebra this reduces to:

ε LOO = 1 n n i=1 G(x (i) ) -G PC (x (i) ) 1 -h i 2 , (27)
where h i is the i-th diagonal term of matrix A(A T A) -1 A T (matrix A is defined in Eq.( 22)) and G PC (•) is now the PC expansion built up from the full experimental design X .

As a conclusion, when using a least-square minimization technique to compute the coefficients of a PC expansion, an a posteriori estimator of the mean-square error is readily available. This allows one to compare PCEs obtained from different truncation schemes and select the best one according to the leave-one-out error estimate.

2.6 Post-processing for sensitivity analysis

Statistical moments

The truncated PC expansion Ŷ = G P C (X) = α∈A ŷα Ψ α (X) contains all the information about the statistical properties of the random output Y = G(X). Due to the orthogonality of the PC basis, mean and standard deviation of Ŷ may be computed directly from the coefficients ŷ. Indeed, since

Ψ 0 ≡ 1, we get E [Ψ α (X)] = 0 ∀ α = 0.
Thus the mean value of Ŷ is the first term of the series:

E Ŷ = E α∈A ŷα Ψ α (X) = ŷ0 . (28)
Similarly, due to Eq.( 10) the variance of Ŷ may be cast as:

σ 2 Ŷ def = Var Ŷ = E Ŷ -ŷ0 2 = α∈A α =0 ŷ2 α . (29)
In other words the mean and variance of the random response may be obtained by a mere combination of the PCE coefficients once the latter have been computed.

Sobol' decomposition and indices

As already discussed in Chapter 4, global sensitivity analysis is based on Sobol' decomposition of the computational model G (a.k.a generalized ANOVA decomposition), which reads [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF]:

G(x) = G 0 + d i=1 G i (x i ) + 1≤i<j≤d G ij (x i , x j ) + • • • + G 12...d (x), (30)
that is, as a sum of a constant G 0 , univariate functions

{G i (x i ) , 1 ≤ i ≤ d}, bivariate functions {G ij (x i , x j ) , 1 ≤ i < j ≤ d}, etc.
A recursive construction is obtained by the following recurrence relationship:

G 0 = E [G(X)] G i (x i ) = E [G(X)|X i = x i ] -G 0 G ij (x i , x j ) = E [G(X)|X i , X j = x i , x j ] -G i (x i ) -G j (x j ) -G 0 . (31)
Using the set notation for indices

A def = {i 1 , . . . , i s } ⊂ {1, . . . , d} , (32)
the Sobol' decomposition in Eq.( 30) reads:

G(x) = G 0 + A⊂{1, ... ,d} A =∅ G A (x A ), (33)
where x A is a subvector of x which only contains the components that belong to the index set A. It can be proven that the summands are orthogonal with each other:

E [G A (x A ) G B (x B )] = 0 ∀ A, B ⊂ {1, . . . , d} , A = B. (34)
Using this orthogonality property, one can decompose the variance of the model output

V def = Var [Y ] = Var     A⊂{1, ... ,d} A =∅ G A (x A )     = A⊂{1, ... ,d} A =∅ Var [G A (X A )] (35)
as the sum of so-called partial variances defined by:

V A def = Var [G A (X A )] = E G 2 A (X A ) . (36)
The Sobol' index attached to each subset of variables

A def = {i 1 , . . . , i s } ⊂ {1, .
. . , d} is finally defined by:

S A = V A V = Var [G A (X A )] Var [Y ] . (37)
First-order Sobol' indices quantify the portion of the total variance V that can be apportioned to the sole input variable X i :

S i = V i V = Var [G i (X i )] Var [Y ] . (38)
Second-order indices quantify the joint effect of variables (X i , X j ) that cannot be explained by each single variable separately:

S ij = V ij V = Var [G ij (X i , X j )] Var [Y ] . (39)
Finally, total Sobol' indices S tot i quantify the total impact of a given parameter X i including all of its interactions with other variables. They may be computed by the sum of the Sobol' indices of any order that contain X i :

S tot i = A i S A . (40)
Amongst other methods, Monte Carlo estimators of the various indices are available in the literature and thoroughly discussed in (see Variance-based sensitivity analysis: Theory and estimation algorithms).

Their computation usually requires 10 3-4 runs of the model G for each index, which leads to a global computational cost that is not affordable when G is expensive-to-evaluate.

Sobol' indices and PC expansions

As can be seen by comparing Eqs.( 18) and ( 33), both polynomial chaos expansions and Sobol' decomposition are sums of orthogonal functions. Taking advantage of this property, it is possible to derive analytic expressions for Sobol' indices based on a PC expansion, as originally shown in Sudret [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF].

For this purpose let us consider the set of multivariate polynomials Ψ α which depend only on a subset of variables A = {i 1 , . . . , i s } ⊂ {1, . . . , d}:

A A = {α ∈ A : α k = 0 if and only if k ∈ A} . (41)
The union of all these sets is by construction equal to A. Thus we can reorder the terms of the truncated PC expansion so as to exhibit the Sobol' decomposition:

G PC (x) = y 0 + A⊂{1, ... ,d} A =∅ G PC A (x A ) where G PC A (x A ) def = α∈A A y α Ψ α (x). ( 42
)
Consequently, due to the orthogonality of the PC basis, the partial variance V A reduces to:

V A = Var G PC A (X A ) = α∈A A y 2 α . (43)
In other words, from a given PC expansion, the Sobol' indices at any order may be obtained by a mere combination of the squares of the coefficients. More specifically, the PC-based estimator of the first-order Sobol' indices read:

Ŝi = α∈Ai ŷ2 α α∈A , α =0 ŷ2 α where A i = {α ∈ A : α i > 0 , α j =i = 0} . ( 44
)
and the total PC-based Sobol' indices read:

Ŝtot i = α∈A tot i ŷ2 α α∈A , α =0 ŷ2 α A tot i = {α ∈ A : α i > 0} . (45)

Summary

Polynomial chaos expansions allow one to cast the random response G(X) as a truncated series expansion.

By selecting an orthonormal basis w.r.t. the input parameter distributions, the corresponding coefficients can be given a straightforward interpretation: the first coefficient y 0 is the mean value of the model output whereas the variance is the sum of the squares of the remaining coefficients. Similarly, the Sobol' indices are obtained by summing up the squares of suitable coefficients. Note that in low dimension (d < 10)

the coefficients can be computed by solving a mere ordinary least-square problem. In higher dimensions advanced techniques leading to sparse expansions must be used to keep the total computational cost (measured in terms of the size N of the experimental design) affordable. Yet the post-processing to get the Sobol' indices from the PCE coefficients is independent of the technique used.

3 Gaussian process-based sensitivity analysis

A short introduction to Gaussian processes

Let us consider a probability space (Ω Z , F Z , P Z ), a measurable space (S, B((S)) and an arbitrary set T .

A stochastic process Z(x), x ∈ T , is Gaussian if and only if for any finite subset C ⊂ T , the collection of random variables Z(C) has a Gaussian joint distribution. In our framework, T and S represent the input and the output spaces. Therefore, we have T = R d and S = R.

A Gaussian process is entirely specified by its mean m(x

) = E Z [Z(x)] and covariance function k(x, x ) = cov Z (Z(x), Z(x )
) where E Z and cov Z denote the expectation and the covariance with respect to (Ω Z , F Z , P Z ).

The covariance function k(x, x ) is a positive definite kernel. It is often considered stationary i.e. k(x, x ) is a function of x-x . The covariance kernel is the most important term of a Gaussian process regression.

Indeed, it controls the smoothness and the scale of the approximation. A popular choice for k(x, x ) is the stationary isotropic squared exponential kernel defined as :

k(x, x ) = σ 2 exp - 1 2θ 2 ||x -x || 2 .
It is parametrized by the parameter θ -also called characteristic length scale or correlation lengthand the variance parameter σ 2 . We give in Figure 1 examples of realizations of Gaussian processes with stationary isotropic squared exponential kernels.

We observe that m(x) is the trend around which the realizations vary, σ 2 controls the range of their variation and θ controls their oscillation frequencies. We highligh that Gaussian processes with squared exponential covariance kernels are infinitely differentiable almost surely. As mentioned in [START_REF] Stein | Interpolation of Spatial Data[END_REF], this choice of kernel can be unrealistic due to its strong regularity.

Gaussian process regression models

The principle of Gaussian process regression is to consider that the prior knowledge about the computational model G(x), x ∈ R d , can be modeled by a Gaussian process Z(x) with a mean denoted by m(x)

and a covariance kernel denoted by k(x, x ). Roughly speaking, we consider that the true response is a realization of Z(x). Usually, the mean and the covariance are parametrized as follows:

m(x) = f T (x)β, (46)
and

k(x, x ) = σ 2 r(x, x ; θ), (47)
where f T (x) is a vector of p prescribed functions and β, σ The shaded areas represent the point-wise 95% confidence intervals.

Predictive distribution

Consider an experimental design X = x (1) , . . . , x (n) , x (i) ∈ R d , and the corresponding model responses

Y = G(X ).
The predictive distribution of G(x) is given by:

[Z(x)|Z(X ) = Y, σ 2 , θ] ∼ GP (m n (x), k n (x, x )) , (48)
where In these expressions R = [r(x i , x j ; θ)] i,j=1,...,n , r(x) = [r(x,

m n (x) = f T (x) β + r T (x)R -1 Y -F β , (49)
k n (x, x ) = σ 2   1 -f T (x) r T (x)   0 F T F R     f (x ) r(x )     , (50
x (i) ; θ)] i=1,...,n , F = [f T (x (i) )] i=1,...,n and β = F T R -1 F -1 F T R -1 Y. ( 51
)
The term β denotes the posterior distribution mode of β obtained from the improper non-informative prior distribution π(β) ∝ 1 [START_REF] Robert | The Bayesian choice: from decision-theoretic foundations to computational implementation[END_REF].

Remark. The predictive distribution is given by the Gaussian process Z(x) conditioned by the known observations Y. The Gaussian process regression metamodel is given by the conditional expectation m n (x) and its mean squared error is given by the conditional variance k n (x, x). An illustration of m n (x) and

k n (x, x) is given in Figure 2.
The reader can note that the predictive distribution [START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF] integrates the posterior distribution of β.

However, the hyper-parameters σ 2 and θ are not known in practice and shall be estimated with the maximum likelihood method [28; 46] or a cross-validation strategy [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of gaussian processes with model misspecification[END_REF]. Then, their estimates are plugged in the predictive distribution. The restricted maximum likelihood estimate of σ 2 is given by:

σ2 = (Y -F β) T R -1 (Y -F β) n -p . (52)
Unfortunately, such a closed form expression does not exist for θ and it has to be numerically estimated.

Remark. Gaussian process regression can easily be extended to the case of noisy observations. Let us suppose that Y is tainted by a white Gaussian noise ε :

Y obs = Y + σ ε (X )ε.
The term σ ε (X ) represents the standard deviation of the observation noise. The mean and the covariance of the predictive distribution [Z(x) obs |Z(X ) = Y obs , σ 2 , θ] is then obtained by replacing in Equations ( 49), ( 50) and ( 51) the correlation matrix R by σ 2 R + ∆ ε where ∆ ε is a diagonal matrix given by : 2) ) . . .

∆ ε =         σ ε (x (1) ) σ ε (x (
σ ε (x (n) )        
.

We emphasize that the closed form expression for the restricted maximum likelihood estimate of σ 2 does not exist anymore. Therefore, this parameter has to be numerically estimated.

Sequential design

To improve the global accuracy of the GP model, it is usual to augment the initial design set X with new points. An important feature of Gaussian process regression is that it provides an estimate of the model mean-square error through the term k n (x, x ) (50) which can be used to select these new points.

The most common but not efficient sequential criterion consists in adding the point x (n+1) where the mean-square error is the largest:

x (n+1) = arg max x k n (x, x). (53)
More efficient criteria can be found in Bates et al. [START_REF] Bates | Experimental design and observation for large systems[END_REF]; van Beers and Kleijnen [START_REF] Van Beers | Customized sequential designs for random simulation experiments: Kriging metamodelling and bootstrapping[END_REF]; Le Gratiet and Cannamela [START_REF] Le Gratiet | Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes[END_REF].

Model selection

To build up a GP model, the user has to make several choices. Indeed, the vector of functions f (x) and the class of the correlation kernel r(x, x ; θ) need to be set (see Rasmussen and Williams [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for different examples of correlation kernels). These choices and the relevance of the model are tested a posteriori with a validation procedure. If the number n of observations is large, an external validation may be performed on a test set. Otherwise, a cross-validation procedure may be used. An interesting property of GP models is that a closed form expression exists for the cross-validation predictive distribution, see for instance Dubrule [START_REF] Dubrule | Cross validation of kriging in a unique neighborhood[END_REF]. It allows for deriving efficient methods of parameter estimation [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of gaussian processes with model misspecification[END_REF] or sequential design [START_REF] Le Gratiet | Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes[END_REF].

Some usual stationary covariance kernel are listed below.

The squared exponential covariance function. The form of this kernel is given by:

k(x, x ) = σ 2 exp - 1 2θ 2 ||x -x || 2 .
This covariance function corresponds to Gaussian processes which are infinitely differentiable in mean square and almost surely. We illustrate in Figure 3 The ν-Matérn covariance function. This covariance kernel is defined as follow (see [START_REF] Stein | Interpolation of Spatial Data[END_REF]):

k ν (h) = 2 1-ν Γ(ν) √ 2||h|| θ ν K ν √ 2ν||h|| θ ,
where ν is the regularity parameter, K ν is a modified Bessel Function and Γ is the Euler Gamma function. A Gaussian process with a ν-Matérn covariance kernel is ν-Hölder continuous in mean square and ν -Hölder continuous almost surely with ν < ν. Three popular choice of ν-Matérn covariance kernels are the ones for ν = 1/2, ν = 3/2 and ν = 5/2 :

k ν=1/2 (h) = exp - ||h|| θ , k ν=3/2 (h) = 1 + √ 3||h|| θ exp - √ 3||h|| θ , and
k ν=5/2 (h) = 1 + √ 5||h|| θ + 5||h|| 2 3θ 2 exp - √ 5||h|| θ .
We illustrate in Figure 4 the 1-dimensional ν-Matérn kernel for different values of ν.

The γ-exponential covariance function. This kernel is defined as follow:

k γ (h) = exp - ||h|| θ γ .
For γ < 2 the corresponding Gaussian process are not differentiable in mean square whereas for γ = 2 is is infinitely differentiable (it corresponds to the squared exponential kernel). We illustrate in Figure 5

Sensitivity analysis

To perform a sensitivity analysis from a GP model, two approaches are possible. The first one consists in substituting the true model G(x) with the mean of the conditional Gaussian process m n (x) in [START_REF] Schoebi | Polynomial-Chaos-based Kriging[END_REF].

However, it may provide biased sensitivity index estimates. Furthermore it does not allow one to quantify the error on the sensitivity indices due to the metamodel approximation. The second one consists in substituting G(x) by a Gaussian process Z n (x) having the predictive distribution

[Z(x)|Z(X ) = Y, σ 2 , θ]
shown in [START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF]. This approach makes it possible to quantify the uncertainty due to the metamodel approximation and allows for building unbiased index estimates.

Main effects visualization

From now on, the input parameter x ∈ R d is considered as a random input vector X = (X 1 , . . . , X d ) with independent components. Before focusing on variance-based sensitivity indices, the inference about the main effects is studied in this section. Main effects are a powerful tool to visualize the impact of each input variable on the model output (see e.g. Oakley and O'Hagan [START_REF] Oakley | Probabilistic sensitivity analysis of complex models a Bayesian approach[END_REF]). The main effect of the group of input variables X A , A ⊂ {1, . . . , d} is defined by

E [G(X)|X A ].
Since the original model G may be time-consuming to evaluate, it is substituted for by its approximation, i.e. G(X)

≈ E [Z n (X)|X A ],
where

Z n (x) ∼ [Z(x)|Z(X ) = Y, σ 2 , θ]. Since E [Z n (X)|X A ]
E [Z n (X)|X A ] ∼ GP (E [m n (X)|X A ] , E [E [k n (X, X )|X A ] |X A ]) . (54)
The

term E [m n (X)] represents the approximation of E [G(X)|X A ] and E [E [k n (X, X )|X A ] |X A ]
is the mean-square error due to the metamodel approximation. Therefore, with this method one can quantify the error on the main effects due to the metamodel approximation. For more detail about this approach, the reader is referred to Oakley and O'Hagan [START_REF] Oakley | Probabilistic sensitivity analysis of complex models a Bayesian approach[END_REF]; Marrel et al. [START_REF] Marrel | Calculations of Sobol indices for the Gaussian process metamodel[END_REF].

Variance of the main effects

Although the main effect enables one to visualize the impact of a group of variables on the model output, it does not quantify it. To perform such an analysis, consider the variance of the main effect:

V A = Var (E [Z n (X)|X A ]) , (55)
or its normalized version which corresponds to the Sobol' index:

S A = V A V = Var (E [Z n (X)|X A ]) Var (Z n (X)) . ( 56
)
Sobol' indices are the most popular measures to carry out a sensitivity analysis since their value can easily be interpreted as the part of the total variance due to a group of variables. However, in contrary to the partial variance V A , it does not provide information about the order of magnitude of the contribution to the model output variance of variable group X A .

Analytic formulae

The above indices are studied in Oakley and O'Hagan [START_REF] Oakley | Probabilistic sensitivity analysis of complex models a Bayesian approach[END_REF] where the estimation of V A and V is performed separately. Indeed, computing the Sobol' index S A requires considering the joint distribution of V A and V , which makes it impossible to derive analytic formulae. According to Oakley and O'Hagan [START_REF] Oakley | Probabilistic sensitivity analysis of complex models a Bayesian approach[END_REF], closed form expressions in terms of integrals can be obtained for the two quantities E Z [V A ] and Var Z (V A ). The quantity E Z [V A ] is the sensitivity measure and Var Z (V A ) represents the error due to the metamodel approximation. Nevertheless, V A is not a linear transform of Z n (X) and its full distribution cannot be established.

Variance estimates with Monte-Carlo integration

To evaluate the Sobol' index S A , it is possible to use the pick-freeze approaches presented in Chapter 4 and in Sobol [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF]; Sobol et al. [START_REF] Sobol | Estimating the approximation error when fixing unessential factors in global sensitivity analysis[END_REF]; Janon et al. [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF]. By considering the formula given in Sobol [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF], S A can be approximated by:

S A,N = 1 N N i=1 Z n (X (i) )Z n (X (i) A ) -1 2N N i=1 Z n (X (i) ) + Z n (X (i) A ) 2 1 N N i=1 Z n (X (i) ) 2 -1 2N N i=1 Z n (X (i) ) + Z n (X (i) A ) 2 , (57)
where (X (i) , X

A ) i=1,...,N is a N -sample from the random variable (X, X ∼A ).

In particular, this approach avoids to compute the integrals presented in Oakley and O'Hagan [START_REF] Oakley | Probabilistic sensitivity analysis of complex models a Bayesian approach[END_REF] and thus simplify the estimation of V A and V . Furthermore, it takes into account their joint distribution.

Remark. This result can easily be extended to the total Sobol' index

S tot i = A⊃ i S A .
The reader is referred to Sobol et al. [START_REF] Sobol | Estimating the approximation error when fixing unessential factors in global sensitivity analysis[END_REF] and Variance-based sensitivity analysis: Theory and estimation algorithms in this handbook for examples of pick-freeze estimates of S tot A .

Numerical estimates of Sobol' indices by Gaussian process sampling

The sensitivity index S A,N (57) is obtained after substituting the Gaussian process Z n (x) for the original computational model G(x). Therefore, it is a random variable defined on the same probability space as Z n (x). The aim of this section is to present a simple methodology to get a sample S A,N of S A . From this sample, an estimate of S A (56) and a quantification of its uncertainty can be deduced.

Sampling from the Gaussian predictive distribution

To obtain a realization of S A,N , one has to obtain a sample of Z n (x) on (X (i) , X

A ) i=1,...,N and then use Eq. ( 57). To deal with large N , an efficient strategy is to sample Z n (x) using the Kriging conditioning method, see for example Chilès and Delfiner [START_REF] Chilès | Geostatistics: modeling spatial uncertainty[END_REF]. Consider first the unconditioned, zero-mean Gaussian process:

Z(x) = GP (0, k(x, x )) . ( 58
)
Then, the Gaussian process:

Zn (x) = m n (x) -mn (x) + Z(x), ( 59
)
where mn (x

) = f T (x) β + r T (x)R -1 Z(X ) -F β and β = F T R -1 F -1 F T R -1 Z(X )
has the same distribution as Z n (x). Therefore, one can compute realizations of Z n (x) from realizations of Z(x). Since Z(x) is not conditioned, the problem is numerically easier. Among the available Gaussian process sampling methods, several can be mentioned: Cholesky decomposition [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], Fourier spectral decomposition [START_REF] Stein | Interpolation of Spatial Data[END_REF], Karhunen-Loeve spectral decomposition [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] and the propagative version of the Gibbs sampler [START_REF] Lantuéjoul | Simulation of a Gaussian random vector: A propagative version of the Gibbs sampler[END_REF].

Remark. Let suppose that a new point x (n+1) is added to the experimental design X . A classical result of conditional probability implies that the new predictive distribution

[Z(x)|Z(X ) = Y, Z(x (n+1) ) = G(x (n+1) ), σ 2 , θ] is identical to [Z n (x)|Z n (x (n+1) ) = G(x (n+1) ), σ 2 , θ].
Therefore, Z n (x) can be viewed as an unconditioned Gaussian process and, using the Kriging conditioning method, realizations of [Z(x)|Z(X ) = Y, Z(x (n+1) ) = G(x (n+1) ), σ 2 , θ] can be derived from realizations of Z n (x) using the following equation:

Z n+1 (x) = k n (x (n+1) , X ) k n (x (n+1) , x (n+1) ) G(x (n+1) ) -Z n (x (n+1) ) + Z n (x). ( 60
)
Therefore, it is easy to calculate a new sample of S A,N after adding a new point x (n+1) to the experimental design set X . This result is used in the R CRAN package "sensitivity" to perform sequential design for sensitivity analysis using a Stepwise Uncertainty Reduction (SUR) strategy [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF].

Meta-model and Monte-Carlo sampling errors

Let us denote by S N A,i , i = 1, . . . , m a sample set of S A,N [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -contributions to structural reliability and stochastic spectral methods[END_REF] where of size m > 0. From this sample set, the following unbiased estimate of S A can be deduced:

ŜA = 1 m m i=1 S N A,i . (61)
with variance:

σ2 ŜA = 1 m -1 m i=1 S N A,i -ŜA 2 . (62)
The term σ2

ŜA represents the uncertainty on the estimate of S A (56) due to the metamodel approximation.

Therefore, with the presented strategy, one can both obtain an unbiased estimate of the sensitivity index S A and a quantification of its uncertainty.

Finally it may be of interest to evaluate the error due to the pick-freeze approximation and to compare it to the error due to the metamodel. To do so, one can use the central limit theorem [31; 17] or a bootstrap procedure [START_REF] Gratiet | A bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF]. In particular, a methodology to evaluate the uncertainty on the sensitivity index due to both the Gaussian process and to the pick-freeze approximations is presented in Le Gratiet et al. [START_REF] Gratiet | A bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF]. It makes it possible to determine the value of N such that the pick-freeze approximation error is negligible compared to that of the metamodel.

Summary

Gaussian Process regression makes it possible to perform sensitivity analysis on complex computational models using a limited number of model evaluations. An important feature of this method is that one can propagate the Gaussian process approximation error to the sensitivity index estimates. This allows the construction of sequential design strategies optimized for sensitivity analysis. It also provides a powerful tool to visualize the main effect of a group of variables and the uncertainty of its estimate.

Another advantage of this approach is that Gaussian process regression has been thoroughly investigated in the literature and can be used in various problems. For example, the method can be adapted for non-stationary numerical models by using a treed Gaussian process as in Gramacy and Taddy [START_REF] Gramacy | Categorical inputs, sensitivity analysis, optimization and importance tempering with tgp version 2, an r package for treed gaussian process models[END_REF].

Furthermore, it can also be used for multifidelity computer codes, i.e. codes which can be run at multiple level of accuracy (see Le Gratiet et al. [START_REF] Gratiet | A bayesian approach for global sensitivity analysis of (multifidelity) computer codes[END_REF]).

Applications

In this section, metamodel-based sensitivity analysis is illustrated on several academic and engineering examples.

Ishigami function

The Ishigami function is given by:

G(x 1 , x 2 , x 3 ) = sin(x 1 ) + 7 sin(x 2 ) 2 + 0.1x 4 3 sin(x 1 ). (63)
The input distributions of X 1 , X 2 and X 3 are uniform over the interval [-π, π] 3 . This is a classical academic benchmark for sensitivity analysis, with first-order Sobol' indices:

S 1 = 0.3138 S 2 = 0.4424 S 3 = 0. ( 64
)
To compare polynomial chaos expansions and Gaussian process modeling on this example, experimental designs of different sizes n are considered. For each size n, 100 Latin Hypercube Sampling sets (LHS) are computed so as to replicate the procedure and assess statistical uncertainty.

For the polynomial chaos approach, the coefficients are calculated based on a degree-adaptive LARS strategy (for details, see Blatman and Sudret [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]), resulting in a sparse basis set. The maximum polynomial degree is adaptively selected in the interval 3 ≤ p ≤ 15 based on LOO cross-validation error estimates (see Eq. ( 27)).

For the Gaussian process approach, a tensorized Matérn-5/2 covariance kernel is chosen (see Rasmussen

and Williams [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]) with trend functions given by:

f T (x) = 1 x 2 x 2 2 x 3 1 x 3 2 x 4 1 x 4 2 . (65)
The hyper-parameters θ are estimated with a Leave-One-Out cross validation procedure while the parameters β and σ 2 are estimated with a restricted maximum likelihood method.

First we illustrate in Figure 6 the accuracy of the models with respect to the sample size n. The Nash-Sutcliffe model efficiency coefficient (also called predictivity coefficient) is defined as follows:

Q 2 = 1 - ntest i=1 (G(x (i) ) -Ĝ(x (i) )) 2 ntest i=1 (G(x (i) ) -Ḡ) 2 , Ḡ = 1 n test ntest i=1 G(x (i) ), (66) 20
where Ĝ(x (i) ) is the prediction given by the polynomial chaos or the Gaussian process regression model on the i th point of a test sample of size n test = 10, 000. This test sample set is randomly generated from a uniform distribution. The closer Q 2 is to 1, the more accurate the metamodel is. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 40 60 80 100 120 0.2 0.4 0.6 0.8 1.0

Polynomial chaos n Q2 q q q q q q q q q q 40 60 80 100 120 We emphasize that checking the metamodel accuracy (see Figure 6) is very important since a metamodelbased sensitivity analysis provides sensitivity indices for the metamodel and not for the true model G(x).

Therefore, the estimated indices are relevant only if the considered surrogate model is accurate.

Figure 7 shows the Sobol' index estimates with respect to the sample size n. For the Gaussian process regression approach, the convergence for is reached for n = 100. It corresponds to a Q 2 coefficient greater than 90%. Convergence of the PCE approach is somewhat faster, with comparable accuracy achieved 21 with n = 60 and almost perfect accuracy for n = 100. Therefore, the convergence of the estimators of the Sobol' indices in Eqs. [START_REF] Marelli | UQLab: A framework for uncertainty quantification in Matlab[END_REF] to ( 39) is expected to be comparable to that of Q 2 . Note that the PCE approach also provides second order-and total Sobol' indices for free, as shown in Sudret [START_REF] Sudret | Polynomial chaos expansions and stochastic finite element methods[END_REF].

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Polynomial chaos n Sobol q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 40 60 80 100 120 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 40 60 80 100 120 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q 40 60 80 Gaussian process regression n Sobol q q q q q q q 40 60 80 100 120 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q 40 60 80

G-Sobol function

The G-Sobol function is given by :

G(x) = d i=1 |4x i -2| + a i 1 + a i , a i ≥ 0. ( 67
)
To benchmark the described metamodel-based sensitivity analysis methods in higher dimension, we select d = 15. The exact first-order Sobol' indices S i are given by the following equations: with n = 100 samples in the experimental design. This is expected due to the good accuracy of both metamodels for all the n considered (see Figure 8).

V i = 1 3(1 + a i ) 2 , i = 1, . . . , d, V = d i=1 (1 + V i ) -1, S i = V i /V.
Finally, Table 2 provides the Sobol' index estimates median and root mean square error for n = 100 and n = 500. As presented in Figure 9, the estimates of the largest Sobol' indices are very accurate. Note that the remaining first order indices are insignificant. One can observe that the RMS error over the 100 LHS replications is slightly smaller when using PCE for both n = 100 and n = 500 ED points. Note that the second order-and total Sobol' indices are also available for free when using PCE.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 100 200 300 400 500 0.80 0.85 0.90 0.95 1.00

Polynomial chaos

n Q2 q q q q q q q q q q q q 100 200 300 400 500 q q q q q q q q 100 200 300 400 500 Polynomial chaos n Sobol q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q 100 200 300 400 500 Gaussian process regression n Sobol q q q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q 100 200 300 400 500

Morris function

The Morris function is given by [START_REF] Schoebi | Polynomial-Chaos-based Kriging[END_REF]:

G(x) = 20 i=1 β i w i + 20 i<j β ij w i w j + 20 i<j<l β ijl w i w j w l + 5w 1 w 2 w 3 w 4 (70)
where X i ∼ U(0, 1), i = 1, . . . , 20 and w i = 2(x i -1/2) for all i except for i = 3, 5, 7 where w i = 2(1.1x i /(x i + 0.1) -1/2). The coefficients are defined as β i = 20, i = 1, . . . , 10; β ij = -15, i, j = 1, . . . , 6; Finally, Table 3 provides the a detailed breakdown of the Sobol' index estimates, including median and root mean square error (RMSE), for n = 100 and n = 500. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 100 200 300 400 500 q q q q q q q q q q q q q q q q q q q q q q q q q 100 200 300 400 500 Polynomial chaos n Sobol q q q q q q q q q q q q q q q q q q q q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q 100 200 300 400 500 Gaussian process regression n Sobol q q q q q q q q q q q q 100 200 300 400 500 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q 100 200 300 400 500

β ijl = -

Maximum deflection of a truss structure

Sensitivity analysis is also of great interest for engineering models whose input parameters may have different distributions. As an example consider the elastic truss structure depicted in Figure 12 (see e.g. Blatman and Sudret [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF]). This truss is made of two types of bars, namely horizontal bars with cross-section A 1 and Young's modulus (stiffness) E 1 on the one hand oblique bars with cross-section A 2 and Young's modulus (stiffness) E 2 on the other hand. The truss is loaded with six vertical loads applied on the top chord. Of interest is the maximum vertical displacement (called deflection) at mid-span. This quantity is computed using a finite element model comprising elastic bar elements. The various parameters describing the behavior of this truss structure are modeled by independent random variables that account for the uncertainty in both the physical properties of the structure and the applied loads. Their distributions are gathered in Table 4. These input variables are collected in the random vector

X = {E 1 , E 2 , A 1 , A 2 , P 1 , . . . , P 6 } . ( 71
)
Using this notation, the maximal deflection of interest is cast as:

u 4 = G FE (X). ( 72
)
Different sparse polynomial chaos expansions are calculated assuming a maximal degree 3 < p < 10 using LARS and the best expansion (in terms of smallest LOO error) is retained. For the Gaussian process regression model, a tensorized Matérn-5/2 covariance kernel is considered with a constant trend function f (x) = 1. The hyper-parameter θ is estimated with a Leave-One-Out cross validation procedure and the parameters β and σ 2 are estimated with the maximum likelihood method.

The first order Sobol' indices obtained from PCE and GP metamodels are reported in Table 5 in the case when the experimental design is of size 100. In decreasing importance order, the important variables are the properties of the chords (horizontal bars), then the loads close to mid-span, namely P 3 and P 4 . The Sobol' indices of the latter are identical due to the symmetry of the model. Then come the loads P 2 and P 5 . The other variables (the loads P 1 and P 6 and the properties of the oblique bars) appear unimportant. The estimates of the three largest first-order Sobol' indices which correspond to variables E 1 , P 3 and P 5 obtained for various sizes n of the LHS experimental design are plotted in Figure 13 as a function of n.

The reference solution is obtained by Monte-Carlo sampling with a sample set of size 6,000,000. Both PCE-and GP-based Sobol' indices converge to stable estimates as soon as n ≥ 60.

Conclusions

Sobol' indices are recognized as good descriptors of the sensitivity of the output of a computational model to its various input parameters. Classical estimation methods based on Monte Carlo simulation are computationally expensive though. The required costs, in the order of 10 3 -10 4 model evaluations, are often not compatible with the advanced simulation models encountered in engineering applications.

For this reason, surrogate models may be first built up from a limited number of runs of the computational model (the so-called experimental design), and the sensitivity analysis is then carried out by substituting the surrogate model for the original one.

Polynomial chaos expansions and Gaussian processes are two popular methods that can be used for this purpose. The advantage of the PCE approach is that the Sobol' indices at any order may be computed analytically once the expansion is available. In this contribution, least-square minimization techniques are presented to compute the PCE coefficients, yet any intrusive or non intrusive method could be used as an alternative. q q q q q q q q q q q q q q q Polynomial chaos expansion q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Gaussian process regression q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 13: Truss structure -First-order Sobol' index estimates as a function of the sample size n for the truss model. The horizontal solid lines represent reference values of input variables E 1 , P 3 and P 5 from a Monte Carlo estimate on 6,000,000 samples.

In contrast Gaussian process surrogate models are used together with Monte Carlo simulation for estimating the Sobol' indices. The advantage of this approach is that the metamodel error can be included in the estimators. Note that bootstrap techniques can be used similarly to calculate and include metamodeling error also for PCE-based sensitivity analysis, as demonstrated by Dubreuil et al. [START_REF] Dubreuil | Construction of bootstrap confidence intervals on sensitivity indices computed by polynomial chaos expansion[END_REF].

As shown in the various comparisons, PCE and GP give similar accuracy (measured in terms of the Q 2 validation coefficient) for a given size of the experimental design in a broad range of applications. The replication of the analyses with different random designs of the same size show a smaller scatter using GP for extremely small designs, whereas PCE becomes more stable for medium-size designs. Selecting the

2 σ = 1 ;Figure 1 :

 211 Figure 1: Examples of Gaussian process realizations with squared exponential kernels and different means.

Figure 2 :

 2 Figure 2: Examples of predictive distribution. The solid line represents the mean of the predictive distribution, the non-solid lines represent some of its realizations and the shaded area represents the 95% confidence intervals based on the variance of the predictive distribution.

3 Figure 3 :

 33 Figure 3: The squared exponential kernel in function of h = xx with different correlation lengths θ and examples of resulting Gaussian process realizations.

Figure 4 : 2 Figure 5 :

 425 Figure 4: The ν-Matérn kernel in function of h = xx with different regularity parameters ν and examples of resulting Gaussian process realizations.

Figure 6 :

 6 Figure 6: Q 2 coefficient as a function of the sample size n for the Ishigami function. For each n, the box-plots represent the variations of Q 2 obtained over 100 LHS replications.

Figure 7 :

 7 Figure 7: First-order Sobol' index estimates as a function of the sample size n for the Ishigami function. The horizontal solid lines represent the exact values of S 1 , S 2 and S 3 . For each n, the box-plot represents the variations obtained from 100 LHS replications. The validation set comprises n test = 10, 000 samples.

  this example, vector a = {a 1 , a 2 , . . . , a d } is equal to: a = {1, 2, 5, 10, 20, 50, 100, 500, 1000, 1000, 1000, 1000, 1000, 1000, 1000} . (69) As in the previous section, different sample sizes n are considered and 100 LHS replications are computed for each n. Sparse polynomial chaos expansions are obtained with the same strategy as for the Ishigami function: adaptive polynomial degree selection with 3 < p < 15 and LARS-based calculation of the coefficients. For the Gaussian process regression model, a tensorized Matérn-5/2 covariance kernel is considered with a constant trend function f (x) = 1. The hyper-parameter θ is estimated with a Leave-One-Out cross validation procedure and the parameters β and σ 2 are estimated with the maximum likelihood method. The accuracy of the metamodels with respect to n is presented in Figure 8. It is computed from a test sample set of size n test = 10, 000. The convergence of the estimates of the first four first-order Sobol' indices is represented in Figure 9. Both metamodel-based estimations yield excellent results already

Figure 8 :

 8 Figure 8: Q 2 coefficient as a function of the sample size n for G-Sobol academic example. For each n, the box-plot represents the variations of Q 2 obtained from 100 LHS.

Figure 9 :

 9 Figure 9: Sobol' index estimates with respect to the sample size n for G-Sobol function. The horizontal solid lines represent the true values of S 1 , S 2 , S 3 and S 4 . For each n, the box-plot represents the variations obtained from 100 LHS.

  10, i, j, l = 1, . . . , 5. The remaining coefficients are set equal to β i = (-1) i and β ij = (-1) i+j and all the rest are zero. The reference values of the first-order Sobol' indices of the Morris function are calculated by a large Monte Carlo-based sensitivity analysis (n = 10 6 ). As in the previous section different sample sizes n are considered and 100 LHS replications are computed for each n. Sparse polynomial chaos expansions are obtained by adaptive polynomial degree selection 5 < p < 13 and LARS-based calculation of the coefficients.The accuracy of the metamodels with respect to n is presented in Figure10. It is computed from a test sample of size n test = 10, 000. As expected due to the complexity and dimensionality of the Morris function, both metamodels show a slower overall convergence rate with the number of samples with respect to the previous examples. Polynomial chaos expansions show in this case remarkably more scattering in their performance for smaller experimental designs with respect to Gaussian process regression. This is likely due to the comparatively large amount of prior information in the form of trend functions provided to the Gaussian process, not used for PCE.The convergence of the estimates of three selected first-order Sobol' indices (the largest S 9 and two intermediate ones S 3 and S 8 ) is represented in Figure11. Both methods perform very well with as few as 250 samples. PCE, however, shows a more standard convergence behavior both in mean value in dispersion. Gaussian process regression retrieves the Sobol' estimates very accurately even with extremely small experimental designs, but no clear convergence pattern can be seen for larger datasets.

Figure 10 :

 10 Figure 10: Q 2 coefficient as a function of the sample size n for the Morris function example. For each n, the box-plot represents the variations of Q 2 obtained from 100 LHS.

Figure 11 :

 11 Figure 11: First-order Sobol' index estimates as a function of the sample size n for the Morris function. The horizontal solid lines represent the exact values of S 3 , S 8 and S 9 . For each n, the box-plot represents the variations obtained from 100 LHS replications.

Figure 12 :

 12 Figure 12: Model of a truss structure with 23 members. The quantity of interest is the maximum displacement at mid-span u 4 .

A 1 ( 4 A 2 4 P 1 -

 14241 m 2 ) Lognormal 2.0 × 10 -3 2.0 × 10 -(m 2 ) Lognormal 1.0 × 10 -3 1.0 × 10 -P 6 (N) Gumbel 5.0 × 10 4 7.5 × 10 3

Table 4 :

 4 Probabilistic input model of the truss structure

	Variable	Distribution	Mean	Standard Deviation
	E 1 , E 2 (Pa) Lognormal	2.1 × 10 11	2.1 × 10 10

Table 5 :

 5 Truss structure -First order Sobol' indices Variable Reference PCE Gaussian Process

	A 1	0.365	0.366	0.384
	E 1	0.365	0.369	0.362
	P 3	0.075	0.078	0.075
	P 4	0.074	0.076	0.069
	P 5	0.035	0.036	0.029
	P 2	0.035	0.036	0.028
	A 2	0.011	0.012	0.015
	E 2	0.011	0.012	0.008
	P 6	0.003	0.005	0.002
	P 1	0.002	0.005	0.000

best technique is in the end problem-dependent, and it is worth comparing the two approaches using the same experimental design, as it can be done in recent sensitivity analysis toolboxes such as OpenTURNS [2] and UQLab [36].

Finally it is worth mentioning that the so-called derivative-based global sensitivity measures (DGSM) originally introduced by Sobol' and Kucherenko [51] can also be computed using surrogate models. In particular, polynomial chaos expansions may be used to compute the DGSM analytically, as shown in Sudret and Mai [61]. The recent combination of polynomial chaos expansions and Gaussian processes into PC-Kriging [48] also appears promising for estimating sensitivity indices from extremely small experimental designs.

0.000 0.000 1.1×10 -4 1.1×10 -3 1.4×10 -4 5.3×10 -4 3.0×10 -4 0.013 0.013 S 8 0.000 0.000 0.000 3.3×10 -4 1.7×10 -5 6.5×10 -4 7.1×10 -4 0.013 0.013 S 9 0.000 0.000 0.000 4.1×10 -4 1.1×10 -5 8.5×10 -4 4.4×10 -4 0.14 0.013 S 10 0.000 0.000 0.000 2.4×10 -4 1.1×10 -5 2.2×10 -4 1.7×10 -4 0.013 0.013 S 11 0.000 0.000 0.000 9.5 ×10 -4 1.2×10 -5 5.5×10 -4 -9.9×10 -5 0.013 0.013 S 12 0.000 0.000 0.000 5.2×10 -4 2.1 ×10 -5 2.6×10 -4 4.1×10 -4 0.013 0.013 S 13 0.000 0.000 0.000 5.1×10 -4 5.9 ×10 -6 9.8×10 -4 4.7×10 -4 0.013 0.013 S 14 0.000 0.000 0.000 8.8×10 -4 1.9 ×10 -5 1.8×10 -4 6.9×10 -4 0.013 0.013 S 15 0.000 0.000 0.000 8.6×10 -4 9.7×10 -6 7.2×10 -4 3.1×10 -4 0.013 0.013 0.000 0.000 0.000 0.083 0.004 0.000 0.000 0.064 0.077 S 13 0.000 0.000 0.000 0.081 0.003 0.000 0.000 0.064 0.074 S 14 0.000 0.000 0.000 0.020 0.003 0.000 0.000 0.070 0.078 S 15 0.000 0.000 0.000 0.140 0.003 0.000 0.000 0.065 0.075 S 16 0.000 0.000 0.000 0.040 0.005 0.001 0.000 0.077 0.074 S 17 0.000 0.000 0.000 0.264 0.004 0.000 0.000 0.065 0.075 S 18 0.000 0.000 0.000 0.084 0.004 0.000 0.000 0.064 0.075 S 19 0.000 0.000 0.000 0.083 0.004 0.000 0.000 0.064 0.076 S 20 0.000 0.000 0.000 0.049 0.004 0.000 0.000 0.064 0.075