G Chastaing
email: gal.chastaing@gmail.com

F Gamboa

C Prieur

Journal of Statistical Computation and Simulation Generalized Sobol sensitivity indices for dependent variables: numerical methods

Keywords: Sensitivity analysis, dependent variables, extended basis, functional decomposition, greedy algorithm, LARS AMS Subject Classification: 62G08, 62H99

The hierarchically orthogonal functional decomposition of any measurable function η of a random vector X = (X 1 , • • • , Xp) consists in decomposing η(X) into a sum of increasing dimension functions depending only on a subvector of X. Even when X 1 , • • • , Xp are assumed to be dependent, this decomposition is unique if the components are hierarchically orthogonal. That is, two of the components are orthogonal whenever all the variables involved in one of the summands are a subset of the variables involved in the other. Setting Y = η(X), this decomposition leads to the definition of generalized sensitivity indices able to quantify the uncertainty of Y due to each dependent input in X [1]. In this paper, a numerical method is developed to identify the component functions of the decomposition using the hierarchical orthogonality property. Furthermore, the asymptotic properties of the components estimation is studied, as well as the numerical estimation of the generalized sensitivity indices of a toy model. Lastly, the method is applied to a model arising from a real-world problem.

Introduction

In simulation models, input parameters can be affected by many sources of uncertainty. The objective of global sensitivity analysis is to identify and to rank the input variables that drive the uncertainty of the model output. The most popular methods are the variance-based ones [START_REF] Saltelli | Sensitivity analysis[END_REF]. Among them, the Sobol indices are widely used [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. This last method relies on the assumption that the inputs are independent. Under this assumption, Hoeffding [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF] shows that the model output can be uniquely cast as a sum of increasing dimension functions, where the integrals of every summand with respect to any of its own variables must be zero. A consequence of these conditions is that all summands of the decomposition are mutually orthogonal. Later, Sobol applies the latter decomposition to sensitivity analysis [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. It results that the global variance can be decomposed as a sum of partial variances that quantify the sensitivity of a set of inputs on the model response.

However, for models featuring dependent inputs, the use of Sobol indices may lead to a wrong interpretation because the sensitivity induced by the dependence between two factors is implicitly included in their Sobol indices. To handle this problem, a straightforward solution consists in computing Sobol sensitivity indices for independent groups of dependent variables. First introduced by Sobol [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], this idea is exploited in practice by Jacques et al. [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]. Nevertheless, this technique implies to work with models having several independent groups of inputs. Furthermore, it does not allow to quantify the individual contribution of each input. A different way to deal with this issue has been initiated by Borgonovo et al. [START_REF] Borgonovo | A new uncertainty importance measure[END_REF][START_REF] Borgonovo | Moment independent importance measures: New results and analytical test cases[END_REF]. These authors define a new measure based on the joint distribution of (Y, X). The new sensitivity indicator of an input X i measures the shift between the output distribution and the same distribution conditionally to X i . This moment free index has many properties and has been applied to some real applications [START_REF] Borgonovo | Moment-independent and variance-based sensitivity analysis with correlations: an application to the stability of a chemical reactor[END_REF][START_REF] Borgonovo | Model emulation and moment-independent sensitivity analysis: An application to environmental modelling[END_REF]. However, the dependence issue remains unsolved as we do not know how the conditional distribution is distorted by the dependence, and how it impacts the sensitivity index. Another idea is to use the Gram-Schmidt orthogonalization procedure. In an early work, Bedford [START_REF] Bedford | Sensitivity indices for (tree-)dependent variables[END_REF] suggests to orthogonalize the conditional expectations and then to use the usual variance decomposition on this new orthogonal collection. Further, the Monte Carlo simulation is used to compute the indices. Following this approach,the Gram-Schmidt process is exploited by Mara et al. [START_REF] Mara | Variance-based sensitivity analysis of computer models with dependent inputs[END_REF], before performing a polynomial regression to approximate the model. In both papers, the decorrelation method depends on the ordering of the variables, making the procedure computationally expensive and difficult to interpret. Following the construction of Sobol indices previously exposed, Xu et al. [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF] propose to decompose the partial variance of an input into a correlated and an uncorrelated contribution in the context of linear models. This last work has been later extended by Li et al. with the concept of HDMR [START_REF] Li | Global sensitivity analysis with independent and/or correlated inputs[END_REF][START_REF] Li | High dimensional model representations[END_REF]. In [START_REF] Li | Global sensitivity analysis with independent and/or correlated inputs[END_REF], the authors suggest to reconstruct the model function using classical bases (polynomials, splines,...), then to deduce the decomposition of the response variance as a sum of partial variances and covariances. Instead of a classical basis, Caniou et al. [START_REF] Caniou | Distribution-based global sensitivity analysis using polynomial chaos expansion[END_REF] use a polynomial chaos expansion to approximate the initial model, and the copula theory to model the dependence structure [START_REF]An introduction to copulas[END_REF]. Thus, in all these papers, the authors choose a type of model reconstruction before proceeding to the splitting of the response variance.

In a previous paper [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF], we revisited the Hoeffding decomposition in a different way, leading to a new definition of the functional decomposition in the case of dependent inputs. Inspired by the pioneering work of Stone [START_REF] Stone | The use of polynomial splines and their tensor products in multivariate function estimation[END_REF] and Hooker [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF], we showed, under a weak assumption on the inputs distribution, that any model function can be decomposed into a sum of hierarchically orthogonal component functions. Hierarchical orthogonality means that two of these summands are orthogonal whenever all variables included in one of the components are also involved in the other. The decomposition leads to generalized Sobol sensitivity indices able to quantify the uncertainty induced by the dependent model inputs.

The goal of this paper is to complement the work done in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF] by providing an efficient numerical method for the estimation of the generalized Sobol sensitivity indices. In our previous paper [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF], we have proposed a statistical procedure based on projection operators to identify the components of the hierarchically orthogonal functional decomposition (HOFD). The method consists in projecting the model output onto constrained spaces to obtain a functional linear system. The numerical resolution of these systems relies on an iterative scheme that requires to estimate conditional expectations at each step. This method is well tailored for independent pairs of dependent variables models. However, it is difficult to apply to more general models because of its computational cost. Hooker [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF] has also worked on the estimation of the HOFD components. This author studies the component estimation via a minimization problem under constraints using a sample grid. In general, this procedure is also quite computationally demanding. Moreover, it requires to get a prior on the inputs distribution at each evaluation point, or, at least, to be able to estimate them properly. In a recent article, Li et al. [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF] reconsider Hooker's work and also identify the HOFD components by a least-squares method. These last authors propose to approximate these components expanded on a suitable basis. They bypass some technical problem of degenerate design matrix by using a continuous descent technique [START_REF] Li | D-morph regression: application to modeling with unknown parameters more than observation data[END_REF].

In this paper, we propose an alternative to directly construct a hierarchical orthogonal basis. Inspired by the usual Gram-Schmidt algorithm, the procedure consists in recursively constructing for each component a multidimensional basis that satisfies the hierarchical orthogonal conditions. This procedure will be referred to as the Hierarchically Orthogonal Gram-Schmidt (HOGS) procedure. Then, each component of the decomposition can be properly estimated by a linear combination of this basis. The coefficients are then estimated by the usual least-squares method. Thanks to the HOGS Procedure, we show that the design matrix has full rank, so the minimization problem admits a unique and explicit solution. Furthermore, we study the asymptotic properties of the estimated components. Nevertheless, the practical estimation of the one-by-one component suffers from the curse of dimensionality when using the ordinary least-squares estimation. To handle this problem, we propose to estimate parameters of the model using variable selection methods. Two usual algorithms are briefly presented, and are adapted to our method. Moreover, the HOGS Procedure coupled with these algorithms is experimented on numerical examples.

The paper is organized as follows. In Section 2, we give and discuss the general results related to the HOFD. We remind Conditions (C.1) and (C.2) under which the HOFD is available. Further, a definition of generalized Sobol sensitivity indices is derived and discussed. Section 3 is devoted to the HOGS Procedure. We introduce the appropriate notation, and outline the procedure in detail. In Section 4, we adapt the least-squares estimation to our problem, and show the consistency of the HOGS Procedure in the least-squares estimation context. Further, we point out the curse of dimensionality, tackled by the use of a penalized minimization scheme. In Section 5, we present numerical applications. The first two examples are toy functions. The objective here is to show the efficiency of the HOGS Procedure that may be coupled with variable selection methods to estimate the sensitivity indices. The last example is an industrial case study. The objective is to detect the inputs that have the strongest impact on a tank distortion.

Generalized Sobol sensitivity indices

Functional ANOVA models are specified by a sum of functions depending on an increasing number of variables. A functional ANOVA model is said to be additive if only main effects are included in the model. It is said to be saturated if all interaction terms are included in the model. The existence and the uniqueness of such a decomposition is ensured by some identifiability constraints. When the inputs are independent, any squared-integrable model function can be exactly represented by a saturated ANOVA model with pairwise orthogonal components. As a result, the contribution of any group of variables to the model response is measured by the Sobol index, bounded between 0 and 1. Moreover, the sum of all the Sobol indices is equal to 1 [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The use of such an index is not excluded in the context of dependent inputs, but the information conveyed by the Sobol indices is redundant, and may lead to a wrong interpretation of the sensitivity in the model. In this section, we remind the main results established in Chastaing et al. [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF] for possibly non-independent inputs. In this case, the saturated ANOVA decomposition holds with weaker identifiability constraints than for the independent case. This leads to a generalization of the Sobol indices that is well suited to perform global sensitivity analysis when the inputs are not necessarily independent. First, we remind the general context and notation. The last part is dedicated to the generalization of the Hoeffding-Sobol decomposition when inputs are potentially dependent. The definition of the generalized sensitivity indices is introduced in the following.

First settings

Consider a measurable function η of a random vector X = (X 1 , • • • , X p ) ∈ R p , p ≥ 1, and let Y be the real-valued response variable defined as

Y : (R p , B(R p ), P X ) → (R, B(R)) X → η(X)
where the joint distribution of X is denoted by P X . For a σ-finite measure ν on (R p , B(R p )), we assume that P X << ν and that X admits a density p X with respect to ν, that is p X = dP X dν .

Also, we assume that η ∈ L 2 R (R p , B(R p ), P X ). As usual, we define the inner product •, • and the norm • of the Hilbert space L 2 R (R p , B(R p ), P X ) as

h 1 , h 2 = h 1 (x)h 2 (x)p X dν(x) = E(h 1 (X)h 2 (X)), h 1 , h 2 ∈ L 2 R (R p , B(R p ), P X ) h 2 = h, h = E(h(X) 2 ), h ∈ L 2 R (R p , B(R p ), P X )
Here E(•) denotes the expectation. Further,

V (•) = E[(• -E(•)) 2 ] denotes the variance, and Cov(•, * ) = E[(• -E(•))( * -E( * ))] the covariance. Let us denote [1 : k] := {1, 2, • • • , k}, ∀ k ∈ N * ,
u := (X u1 , • • • , X uk ).
Conventionally, if u = ∅, |u| = 0, and X ∅ = 1. Also, we denote by X -u the complementary vector of X u (that is, -u is the complementary set of u). The marginal density of X u (respectively X -u ) is denoted by p Xu (resp. p X-u ).

Further, the mathematical structure of the functional ANOVA models is defined through subspaces (H u ) u∈C and (H 0 u ) u∈C of L 2 R (R p , B(R p ), P X ). H ∅ ≡ H 0 ∅ denotes the space of constant functions. For u ∈ C \{∅}, H u is the space of square-integrable functions that depend only on X u . The space H 0 u is defined in a recursive way on u as:

H 0 u = h u ∈ H u , h u , h v = 0, ∀ v ⊂ u, ∀ h v ∈ H 0 v = H u ∩ v⊂u H 0 v ⊥ ,
where ⊂ denotes the strict inclusion, that is

A ⊂ B ⇒ A ∩ B = B.
Further, ⊆ will denote the inclusion when equality is possible.

Generalized Sobol sensitivity indices

Let us suppose that

P X << ν where ν(dx) = ν 1 (dx 1 ) ⊗ • • • ⊗ ν p (dx p ) (C.1)
Our main assumption is :

∃ 0 < M ≤ 1, ∀ u ∈ C \ {∅}, p X ≥ M • p Xu p X-u ν-a.e. (C.2)
Under these conditions, the following result states a general decomposition of η as a saturated functional ANOVA model, under the specific conditions of the spaces H 0 u (defined in Section 2.1), Theorem 2.1 Let η be any function in L 2 R (R p , B(R p ), P X ). Then, under (C.1) and (C.2), there exist unique functions

(η 0 , η 1 , • • • , η {1,••• ,p} ) ∈ H ∅ × H 0 1 × • • • H 0 {1,••
• ,p} such that the following equality holds :

η(X) = η 0 + p i=1 η i (X i ) + 1≤i<j≤p η ij (X i , X j ) + • • • + η {1,••• ,p} (X) = u∈C η u (X u ) a.e. (1)
To get a better understanding of Theorem 2.1, the reader could refer to its proof and further explanations in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF]. Notice that, unlike the Sobol decomposition with independent inputs, the component functions of (1) are hierarchically orthogonal, and no more mutually orthogonal. Thus, from now on, the obtained decomposition (1) will be abbreviated HOFD (for Hierarchically Orthogonal Functional Decomposition). Also, as mentioned in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF], the HOFD is said to be a generalized decomposition because it turns out to be the usual functional ANOVA decomposition when inputs are independent. The general decomposition of the output Y = η(X) given in Theorem 2.1 allows one to decompose the global variance as a simplified sum of covariance terms. Further below, we define the generalized sensitivity indices able to measure the contribution of any group of inputs in the model when inputs may be dependent : Definition 2.2 The sensitivity index S u of order |u| measuring the contribution of X u to the model response is given by :

S u = V (η u (X u )) + u,v∈C\{∅} u ⊆v,v ⊆u Cov(η u (X u ), η v (X v )) V (Y ) (2)
More specifically, the first order sensitivity index S i is given by :

S i = V (η i (X i )) + v∈C\{∅} i ∈v Cov(η i (X i ), η v (X v )) V (Y ) (3)
These indices are called generalized Sobol sensitivity indices because if all inputs are independent, it can be shown that Cov

(η u , η v ) = 0, ∀ u = v [1]. Thus, S u = V (η u (X u ))/V (Y ),

Interpretation of the sensitivity indices

The indices' interpretation is not obvious, as they are no more bounded and they can even be negative. We provide here an interpretation of the first order sensitivity index S i , split into two parts:

S i = V (η i (X i )) V (Y ) VSi + v =∅ i ∈v Cov(η i (X i ), η v (X v )) V (Y ) CoVSi .
The first part, VS i , could be identified as the full contribution of X i , whereas the second part, CoVS i , could be interpreted as the contribution induced by the dependence with the other terms of the decomposition. Thus, CoVS i would play the role of compensation. We detail here this interpretation, and we distinguish five cases, represented and explained further below.

VS i CoVS i (1)
S i > 0

(2)

S i > 0 (3)
S i > 0 (4)
S i < 0 S i = 0(5)
Case [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF]. The full contribution of X i is not important, but the uncertainty is induced by the dependence. Thus, X i has an influence through its dependence with the other variables.

Case [START_REF] Saltelli | Sensitivity analysis[END_REF]. The full contribution of X i is important, and the induced contribution is lower. In this case, X i has a strong influence.

Case [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The uncertainty of the sole X i is important, but weakened by the contribution induced by the dependence.

Case (4). The influence of X i is not direct, as it comes from its dependencies with other factors. This influence, obtained by negative covariances, is more important than the full contribution, so X i may not be so significant.

Case [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]. The case VS i = CoVS i = 0 means that X i is not contributive to the global variability. However, if VS i = -CoVS i = 0, the dependence makes a perfect compensation.

Thus, decisions on the influence of the model inputs can hardly be taken in Case (4) and the perfect compensation of Case [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF]. Also, the interpretation is subject to the index splitting, so it may have an impact on the conclusion drawn in sensitivity analysis. This impact has to be carefully considered and its study remains an open problem.

From now on, the next part is dedicated to the practical use of these indices. The analytical formulation of the HOFD components is rarely available in realistic applications. Indeed, their computation requires to know the mathematical form of η and the distribution of the input parameters. It also implies to search for components in a space constrained by very specific orthogonality conditions. Efficient numerical methods has then to be developed to estimate the generalized sensitivity indices. The following section is devoted to an original estimation scheme of the HOFD components based on two tools: extended bases [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF] and penalized regression [START_REF] Hastie | The elements of statistical learning[END_REF].

The hierarchically orthogonal Gram-Schmidt procedure

In Section 2, we show that each component η u belongs to a subspace H 0 u of L 2 (R p , B(R p ), P X ). Thus, to estimate η u , the most natural approach is to construct a good approximation space of H 0 u . In addition, we have seen that the generalized sensitivity indices are defined for any type of reference measures (ν i ) i∈ [1:p] . From now and until the end, we will assume that ν i , ∀ i ∈ [1 : p], are diffuse measures. Indeed, the non diffuse measures raise additional issues in the results developed further that we will not address in this paper.

In a Hilbert space, it is usual to call in an orthonormal basis to express any of the space element as a linear combination of the components of this basis. Further below, we will define the finite-dimensional spaces

H L u ⊂ H u and H 0,L u ⊂ H 0 u , ∀ u ∈ C
, as linear spans of some orthonormal systems that will be settled later. We denote by Span {B} the set of all finite linear combination of elements of B, also called the linear span of B.

Consider, for any

i ∈ [1 : p], a truncated orthonormal system (ψ i li ) Li li=0 of L 2 (R, B(R), P Xi ), with L i ≥ 1.
Without loss of generality, we simplify the notation, and we assume that

L i = L ≥ 1, for all i ∈ [1 : p]. Also, when there is no con- fusion, ψ i li is written ψ li . Moreover, we set ψ 0 = 1. For any u = {u 1 , • • • , u k } ∈ C, l u = (l u1 , • • • , l uk ) ∈ [1 : L] |u| := [1 : L] × • • • × [1 : L] is the multi-index associated with the tensor-product (⊗ k i=1 ψ lu i ).
To define properly the truncated spaces H L u ⊂ H u , we further assume that

∀u = {u i } k i=1 ∈ C, ∀ l u ∈ [1 : L] |u| , [ k i=1 ψ lu i (x ui )] 2 p X dν(x) < +∞ (C.3) Remark 3.1 A sufficient condition for (C.2) is to have 0 < M 1 ≤ p X ≤ M 2 (see Section 3 of [1]
). In this particular case, it is sufficient to assume that, ∀i ∈

[1 : p], ∀ l i ∈ [1 : L], ( i∈[1:p] ψ li (x i )) 2 dν(x) < +∞ to ensure (C.3). Under (C.3), we define, H L ∅ = Span {1}. Also, we set, ∀ i = j ∈ [1 : p], H L i = Span {1, ψ li , l i ∈ [1 : L]} .
Further, we write the multivariate spaces H L u , for u ∈ C, as

H L u = ⊗ i∈u H L i .
Then, the dimension of

H L u is dim(H L u ) = (L + 1) |u| .
Now, we focus on the construction of the theoretical finite-dimensional spaces (H 0,L u ) u∈C , that corresponds to the constrained subspaces of (H L u ) u∈C . Thus, for all u ∈ C, H 0,L u is defined as

H 0,L u = h u ∈ H L u , h u , h v = 0, ∀ v ⊂ u, ∀ h v ∈ H 0,L v Hence, dim(H 0,L u ) = dim(H L u ) -[ v⊂u v =∅ L |v| + 1] = L |u| .
Given an independent and identically distributed n-sample (y s , x s ) s=1,••• ,n from the distribution of (Y, X), the empirical version Ĥ0,L

u of H 0,L u is defined as Ĥ0,L ∅ = H L ∅ , and
Ĥ0,L u = g u ∈ H L u , g u , g v n = 0, ∀ v ⊂ u, ∀ g v ∈ Ĥ0,L v ,
where •, • n defines the empirical inner product associated with the n-sample. The space Ĥ0,L u varies with sample size n, but for notational convenience, we suppress the dependence on n.

The next procedure is an iterative scheme to construct (H 0,L u ) u∈C and ( Ĥ0,L u ) u∈C by taking into account their specific properties of orthogonality. This numerical method is referred to as the Hierarchically Orthogonal Gram-Schmidt (HOGS) procedure.

Hierarchically Orthogonal Gram-Schmidt Procedure

(1) Initialization. For any i ∈ [1 : p], we use the truncated orthonormal system

(ψ li ) L li=0 . Set φ li = ψ li , ∀ l i ≥ 1 and H 0,L i = Span {φ li , l i ∈ [1 : L]} .
As (φ li ) L li=1 is an orthonormal system, any function

h i ∈ H 0,L i satisfies E(h i (X i )) = 0.
(2) Second order interactions. Let u = {i, j} with i = j ∈ [1 : p]. The space H 0,L u is constructed in a recursive way. By Step (1), we have

H 0,L i = Span {φ li , l i ∈ [1 : L]} and H 0,L j = Span φ lj , l j ∈ [1 : L] . For all l ij = (l i , l j ) ∈ [1 : L] 2 , (a) set φ lij (X i , X j ) = φ li (X i ) × φ lj (X j ) + L k=1 λ i k φ i k (X i ) + L k=1 λ j k φ j k (X j ) + C lij (b) The constants (C lij , (λ i k ) L k=1 , (λ i k ) L k=1
) are determined by solving the hierarchical orthogonal constraints,

   φ lij , φ i k = 0, ∀ k ∈ [1 : L] φ lij , φ j k = 0, ∀ k ∈ [1 : L] φ lij , 1 = 0 Finally, H 0,L ij = Span φ lij , l ij ∈ [1 : L] 2 . Each function h ij ∈ H 0,L ij satis- fies the constraints imposed to H 0 ij .
(3) Higher order interactions. To build a basis

(φ lu ) lu∈[1:L] |u| of H 0,L u , with u = (u 1 , • • • , u k ),
we proceed recursively on k. By the previous steps of the Procedure, we have at hand, for any

v ∈ C such that 1 ≤ |v| ≤ k -1, H 0,L v = Span φ lv , l v ∈ [1 : L] |v| , dim(H 0,L v ) = L |v| .
Now, to obtain φ lu , for all l u = (l u1 , • • • , l uk ) ∈ [1 : L] |u| , we proceed as follows, (a) set

φ lu (X u ) = k i=1 φ lu i (X ui ) + v⊂u v =∅ lv∈[1:L] |v| λ lv,lu φ lv (X v ) + C lu (4) (b) compute the (1+ v⊂u v =∅ L |v| ) coefficients (C lu , (λ lv,lu ) lv∈[1:L] |v| ,v⊂u ) by solving φ lu , φ lv = 0, ∀v ⊂ u, ∀ l v ∈ [1 : L] |v| φ lu , 1 = 0. ( 5
)
The linear system (5), with (4), is equivalent to a sparse matrix system of the form A u φ Λ u = D lu , when C lu has been removed. The matrix

A u φ is a Gramian matrix involving terms E(Φ v1 (X v1 ) t Φ v2 (X v2 )) v1,v2⊂u , with (Φ vi (X vi )) lv i = φ lv i (X vi ), l vi ∈ [1 : L] |vi| , i = 1, 2. Λ u involves the coef- ficients (λ lv,lu ) lv∈[1:L] |v| ,v⊂u , and D lu involves -E( ui∈u φ • lu i Φ vi ) vi⊂u . In Lemma A.2, we show that A u
φ is a definite positive matrix, so the system (5) admits a unique solution. Finally, set

H 0,L u = Span φ lu , l u ∈ [1 : L] |u| .
The construction of ( Ĥ0,L u ) u∈C is very similar to the (H 0,L u ) u∈C one. However, as the spaces ( Ĥ0,L u ) u∈C depend on the observed n-sample, their construction requires to assume that the sample size n is larger than the size L. Thus, face to an expensive model with a limited budget, the number of observations n may be small. In this case, and in view of the HOGS procedure, L should be chosen accordingly. To build Ĥ0,L i , ∀ i ∈ [1 : p], we use the usual Gram-Schmidt procedure on (φ li ) L li=1 to get an orthonormal system ( φli ) L li=1 with respect to the empirical inner product •, • n . To build ( Ĥ0,L u ) u∈C,|u|≥2 , we can simply use the HOGS procedure while replacing

•, • with •, • n . Finally, we denote Ĥ0,L u = Span φlu , l u ∈ [1 : L] |u| ∀ u ∈ C \ {∅}.
In practice, polynomials or splines basis functions [START_REF] Droesbeke | Approches non paramétriques en régression[END_REF] will be considered. In the next section, we discuss the practical estimation of the generalized Sobol sensitivity indices using least-squares minimization, and its consistency. Further, we discuss the curse of dimensionality, and propose some variable selection methods to circumvent it.

Estimation of the generalized sensitivity indices 4.1 Least-Squares estimation

The effects (η u ) u∈C in the HOFD (1) satisfy

(η u ) u∈C = Argmin (hu)u∈C hu∈H 0 u E[(Y - u∈C h u (X u )) 2 ] (6)
Notice that η 0 , the expected value of Y , is not involved in the sensitivity indices estimation. Thus, Y is replaced with Ỹ := Y -E(Y ) in [START_REF] Borgonovo | A new uncertainty importance measure[END_REF]. Also, the residual term η {1,••• ,p} is removed from [START_REF] Borgonovo | A new uncertainty importance measure[END_REF] and is estimated afterwards. In Section 3, we defined the approximating spaces Ĥ0,L u of H 0 u , for u ∈ C \ {∅}. Thus, the minimization problem (6) may be replaced with its empirical version, min

β u lu ∈R, lu∈[1:L] |u| , u∈C 1 n n s=1     ỹs - u⊂[1:p] u =∅ lu∈[1:L] |u| β u lu φlu (x u s )     2 (7)
where ỹs := y sȳ, ȳ := 1 n n s=1 y s , and where every subspace Ĥ0,L u is spanned by the basis functions ( φlu ) lu∈[1:L] |u| constructed according to the HOGS Procedure of Section 3. The equivalent matrix form of ( 7) is

min β Y -X φβ 2 n ( 8
)
where

Y s = y s -ȳ, X φ = φ1 • • • φu • • • ∈ × u∈C M n,L |u| (R), where × u∈C M n,L |u| (R) denotes
the cartesian product of real entries matrices with n rows and L |u| columns.

For u ∈ C, ( φu ) s,lu = φlu (x u s ), and

β = (β u lu ) lu,u , ∀ s ∈ [1 : n], ∀ l u ∈ [1 : L] |u| .

Asymptotic results

The method exposed in Section 3 aims at estimating the ANOVA components, whose uniqueness is ensured by hierarchical orthogonality. However, we would like to make sure that our numerical procedure is robust, i.e. that the estimated summands converge to the theoretical ones. To do that, we suppose that the model function η is approximated by η L , where

η L (X) = u∈C η L u (X), with η L u = lu∈[1:L] |u| β u,0 lu φ lu (X u ),
where β u,0 lu stands for the true parameter. Further, we assume that the estimator ηL of

η L is given ηL (X) := u∈C ηL u (X u ), with ηL u (X u ) = lu∈[1:L] |u| βlu φlu (X u ),
where the coefficients ( βu lu ) lu∈[1:L] |u| ,u∈C are estimated by [START_REF] Borgonovo | Moment-independent and variance-based sensitivity analysis with correlations: an application to the stability of a chemical reactor[END_REF], and φlu ∈ Ĥ0,L u . Here, we are interested in the consistency of ηL when the dimension L |u| is fixed for all u ∈ C. This result is stated in Proposition 4.1.

Proposition 4.1 Assume that Y = η L (X) + ε, where η L (X) = u∈C lu∈[1:L] |u| β u,0 lu φ lu (X u ),
with E(ε) = 0, and

E(ε 2 ) = σ 2 * , E(ε • φ lu (X u )) = 0, ∀ l u ∈ [1 : L] |u| , ∀ u ∈ C. (β 0 = (β u,0 lu ) lu,u is the true parameter).
Further, let us consider the least squares estimation ηL of η L using the sample (y s , x s ) s∈[1:n] from the distribution of (Y, X), and the functions ( φlu ) lu , that is

ηL (X) = u∈C ηL u (X u ), where ηL u (X u ) = lu∈[1:L] |u| βu lu φlu (X u ),
where

β = Argmin β∈Θ , Y -X φβ 2 n
and Θ is a compact set of R m , m being the size of the vector β.

If we assume that (H.1) The distribution P X is equivalent to

⊗ p i=1 P Xi ; (H.2) For any u ∈ C, any l u ∈ [1 : L] |u| , φ lu = 1 and φlu n = 1 (H.3) For any i ∈ [1 : p], any l i ∈ [1 : L], φ2 li < +∞. Then, ηL -η L a.s. → 0 when n → +∞. ( 9
)
The proof of Proposition 4.1 is postponed to Appendix A.

Our aim here is to study how the approximating spaces Ĥ0,L u , constructed with the previous procedure, behave when n → +∞. However, we assume that the dimension of H 0,L u , L |u| , is fixed. By extending the work of Stone [START_REF] Stone | The use of polynomial splines and their tensor products in multivariate function estimation[END_REF], Huang [START_REF] Huang | Projection estimation in multiple regression with application to functional anova models[END_REF] explores the convergence properties of functional ANOVA models when L |u| is not fixed anymore. Nevertheless, the results are obtained for a general model space H u , and its approximating space Ĥu , ∀ u ∈ C. In [START_REF] Huang | Projection estimation in multiple regression with application to functional anova models[END_REF], the author states that if the basis functions are m-smooth and bounded, ηη converges in probability. For polynomials, Fourier transforms or splines, he specifically shows that ηη = O p (n -2m 2m+d ) (See [START_REF] Huang | Projection estimation in multiple regression with application to functional anova models[END_REF] p. 257), when d is the ANOVA order (i.e. η ≃ u∈C |u|≤d η u ), where d can be chosen by the notions of effective dimension [START_REF] Wang | The effective dimension and quasi-monte carlo integration[END_REF]. Thus, even if we show the convergence of ηL for d = p, where p is the model dimension, it is in our interest to have a small order d when p gets large to get a good rate of convergence. However, in practice, even for small d ≪ p, the number of terms blows up with the dimensionality of the problem, and so would the number of model evaluations when using an ordinary least-squares regression scheme. As an illustration, take d = 3, p = 8 and L = 5. In this case, m = 7740 parameters have to be estimated, which could be a difficult task in practice. To handle this kind of problem, many variable selection methods have been considered in the field of statistics. The next section aims at briefly exposing the variable selection methods based on a penalized regression. We particularly focus on the ℓ 0 penalty [START_REF] Temlyakov | Greedy approximation[END_REF] and on the Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

The variable selection methods

For simplicity, we denote by m the number of parameters in [START_REF] Borgonovo | Moment-independent and variance-based sensitivity analysis with correlations: an application to the stability of a chemical reactor[END_REF]. The variable selection methods usually deal with the penalized regression problem

min β Y -X φβ 2 n + λJ(β) (10)
where J(•) is positive valued for β = 0, and where λ ≥ 0 is a tuning parameter. The most intuitive approach is to consider the ℓ 0 -penalty J(β) = β 0 , where β 0 = m j=1 1(β j = 0). Indeed, the ℓ 0 regularization aims at selecting non-zero coefficients, thus at removing the useless parameters from the model. The greedy approximation [START_REF] Temlyakov | Greedy approximation[END_REF] offers a series of strategies to deal with the ℓ 0 -penalty. Nevertheless, the ℓ 0 regularization is a non convex function, and suffers from statistical instability, as mentioned in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Breiman | Better subset regression using the nonnegative garrote[END_REF]. The Lasso regression could be regarded as a convex relaxation of the optimization problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Indeed, the Lasso regression is based on ℓ 1 -penalty, i.e. [START_REF] Bedford | Sensitivity indices for (tree-)dependent variables[END_REF] with J(β) = β 1 , and β 1 = m j=1 |β j |. The Lasso offers a good compromise between a rough selection of non-zero elements, and a ridge regression (J(β) = m j=1 β 2 j ) that only shrinks coefficients, but is known to be stable [START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Buhlmann | Statistics for high-dimensional data[END_REF]. In the following, the proposed method will use either the ℓ 0 or the ℓ 1 regularization. The adaptive forward-backward greedy (FoBa) algorithm proposed in Zhang [START_REF] Zhang | Adaptive forward-backward algorithm for learning sparse representations[END_REF] is exploited here to deal with the ℓ 0 penalization. From a dictionary D that can be large and/or redundant, the FoBa algorithm is an iterative scheme that sequentially selects and discards the element of D that has the least impact on the fit. The aim of the algorithm is to efficiently select a limited number of predictors. The advantage of such an approach is that it is very intuitive, and easy to implement. Initiated by Efron et al. [START_REF] Efron | Least angle regression[END_REF], the modified LARS algorithm is adopted to deal with the Lasso regression. The LARS is a general iterative technique that builds up the regression function by successive steps. The adaptation of LARS to Lasso (the modified LARS) is inspired by the homotopy method proposed by Osborne et al. [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF]. The main advantage of the modified LARS algorithm is that it builds up the whole regularized solution path { β(λ), λ ∈ R}, exploiting the property of piecewise linearity of the solutions with respect to λ [START_REF] Buhlmann | Statistics for high-dimensional data[END_REF][START_REF] Osborne | On the lasso and its dual[END_REF]. Moreover, the modified LARS yields all lasso solutions, making this algorithm theoretically efficient.

In the next part, both the FoBa and the modified LARS algorithms are adapted to our problem and are then compared on numerical examples.

Summary of the estimation procedure

Provided an initial choice of orthonormal system (ψ li ) L li=0,i∈[1:p] , we first construct the approximation spaces Ĥ0,L u of H 0 u for |u| ≤ d, and d ≪ p, using the HOGS Procedure described in Section 3. A HOFD component η u is then a projection onto Ĥ0,L u , whose coefficients ( βu lu ) are defined by least-squares estimation. To bypass the curse of dimensionality, the FoBa algorithm or the modified LARS algorithm is used. Once the HOFD components are estimated, we derive the empirical estimation of the generalized Sobol sensitivity indices. Thus, Ŝi , ∀ i ∈ [1 : p], is given by Ŝi =

1 n n s=1 L li=1 βi li φli (x s i ) 2 V (Y ) + |v|≤d i ∈v 1 n n s=1 L li=1 βi li φli (x s i ) • lv βv lv φlv (x v s ) V (Y )
, where V (Y ) = 1/n n s=1 (y sȳ) 2 .

Application

In this section, we are interested by the numerical efficiency of the HOGS procedure introduced in Section 3, that may be coupled with a penalized regression, as done in Section 4.3. The goal of the following study is to show that our strategy gives a good estimation of the generalized sensitivity indices in an efficient way.

Description

In this study, we compare several numerical strategies summarized further below.

(1) The HOGS Procedure consists in constructing the basis functions that will be used to estimate the components of the functional ANOVA. Further, to estimate the m unknown coefficients, we may use (a) the usual least squares estimation when m < n, and n is the number of model evaluations. This technique is called LSEHOGS. (b) when m ≥ n, the HOGS is coupled with the adaptive greedy algorithm FoBa to solve the ℓ 0 -penalized problem. This is called FoBaHOGS. To relax the ℓ 0 penalization, the modified LARS algorithm may replace the greedy strategy, abbreviated LHOGS, where L stands for LARS. (2) The general method developed in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -application to sensitivity analysis[END_REF], based on projection operators, consists in solving a functional linear system, when the model depends on independent pairs of dependent inputs. This procedure is abbreviated POM for Projection Operators Method. (3) At last, we compare our strategy to a minimization under constraints detailed in [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF], and summarized as

min F Y -XF 2 n D n F = 0 ( 11
)
with

Y s = y s -ȳ, s = 1, • • • , n, X is a matrix composed of 1 and 0 elements, F s,u = η u (x u s )
, and D n is the matrix of hierarchical orthogonal constraints, where the inner product •, • has been replaced with its empirical version •, • n . However, the matrix X is not full rank, so the solution of ( 11) is not unique. This implies that the descent technique used to estimate F may give local solutions and lead to wrong results. Moreover, in a high-dimensional paradigm, the matrix of constraints D n becomes very large, leading to an intractable strategy in practice. To remedy to these issues, we consider that each component is parametrized, and constrained, i.e. we consider the following problem,

min Φ Y -XΦ 2 n D n Φ = 0 ( 12
)
where X s,lu = ψ lu (x u s ), s = 1, • • • , n, with (ψ lu ) lu,u the usual tensor basis of L 2 (R) (polynomial, splines, . . . ). The vector Φ is the set of unknown coefficients, and D n is the matrix of constraints given by [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF], where the inner product •, • has been replaced with its empirical version •, • n , on the parametrized functionals of the decomposition. The Lagrange function associated with ( 12) can be easily derived, and the linear system to be solved is the following

t XX -t D D 0 Φ λ = t XY 0 ,
where λ is the Lagrange multiplier. This procedure, substantially similar to the HOGS Procedure, is abbreviated MUC for Minimization Under Constraints.

In the following, the computational time cost of each model is negligible when compared with the procedures described above. In the next numerical examples, all these strategies are compared in terms of CPU time, as well as mean squared error, defined as

mse(η) = 1 |C| u∈C 1 n n s=1 [η u (x u s ) -η u (x u s )] 2 ,
where the functions ηu are estimated by one of the methodologies described above, and η u are the analytical functions.

Test cases and results

For every model, we consider that the functional ANOVA decomposition is truncated at order d = 2.

Test case 1: the g-Sobol function. Well known in the sensitivity literature [START_REF] Saltelli | Sensitivity analysis[END_REF], the g-Sobol function is given by

Y = p i=1 |4X i -2| + a i 1 + a i , a i ≥ 0,
where the inputs X i are independent and uniformly distributed over [0, 1]. The analytical Sobol indices are given by

S u = 1 D i∈u D i , D i = 1 3(1 + a i ) 2 , D = p i=1 (D i + 1) -1, ∀ u ∈ C.
We take p = 4, and a = (0, 1, 4.5, 99). We choose the Legendre polynomial basis of degree 5. Therefore, the number of parameters to be evaluated is m = 170. The purpose of this study is to show the efficiency of the HOGS strategy, and to compare it with the MUC one. Our aim is also to consolidate the asymptotic result given in Section 4 from a numerical viewpoint. To this end, we make two numerical tests, where The computational efficiency of both strategies are reported in Table 1. Test case 2: the Li function. This polynomial function has been introduced by Li et al. [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF], and is given by the following expression,

S1 S2 S3 S4 S12 S13 S14 S23 S24 S34
Y = g 1 (X 1 , X 2 ) + g 2 (X 2 ) + g 3 (X 3 ), X ∼ N (0, Σ), with g 1 (X 1 , X 2 ) = [a 1 X 1 + a 0 ][b 1 X 2 + b 0 ] g 2 (X 2 ) = c 2 X 2 2 + c 1 X 2 + c 0 g 3 (X 3 ) = d 3 X 3 3 + d 2 X 2 3 + d 1 X 3 + d 0 Σ =   σ 2 1 γσ 1 σ 2 0 γσ 1 σ 2 σ 2 2 0 0 0 σ 2 3   .
The normal distribution does not satisfy Condition (C.2). However, the Gaussian density makes it possible to compute a HOFD decomposition, as done in [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF]. Moreover, if the search of solutions is restricted to the polynomial spaces, the uniqueness of the HOFD components given in [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF] is ensured, whatever the type of distribution. Thus, the analytical form of the ANOVA components and the generalized Sobol indices can be derived in this case. To mimic the work done in [START_REF] Li | General formulation of hdmr component functions with independent and correlated variables[END_REF], we take

a 0 = c 1 = d 0 = 1, a 1 = b 0 = c 2 = d 1 = d 2 = 2 and b 1 = c 0 = d 3 = 3.
The variations are set equal to σ 1 = σ 2 = 0.2, σ 3 = 0.18 and γ = 0.6. For each component, we choose a Hermite basis of degree 10 to apply the HOGS Procedure and the MUC strategy. Thus, the number of parameters to be estimated is equal to m = 330. Further, we consider n = 300 model evaluations to estimate the parameters by the (L/FoBa)HOGS method. We repeat the test over 50 repetitions. We compare it to the MUC and the POM strategies in Table 2 on the estimated sensitivity indices. Table 3 shows the number of non-zero estimated coefficients for the FoBaHOGS and the LHOGS. The averaged elapsed time and the mse(η) computed for each method are also reported in Table 3. Test case 3: the tank pressure model. The case study is a shell closed by a cap and subjected to an internal pressure. Figure 3 illustrates a simulation of tank distortion. The calculation of interest is the von Mises stress [START_REF] Von Mises | Mechanik der festen krper im plastisch deformablen zustand[END_REF] at point y. The yielding of the material occurs as soon as the von Mises stress reaches the material yield strength. The selected point y corresponds to the point for which the von Mises stress is maximal in the tank. We want to prevent the tank from material damage induced by plastic strains. A 2D finite elements model the system using the code ASTER. In order to design reliable structures, a manufacturer wants to identify the most contributive parameters to the von Mises criterion variability. The von Mises criterion depends on three geometrical parameters: R int , T shell and T cap . It also depends on six material parameters, E shell , E cap , σ y,shell , σ y,cap , and P int . Table 4 gives the meaning and the distribution of the eight inputs.

The geometrical parameters are uniformly distributed because of the large choice left for the tank building. The correlation γ between the geometrical parameters is induced by the constraints of manufacturing processes. The physical inputs are normally distributed and their uncertainty are due to the manufacturing process and the properties of the elementary constituents variabilities. The large variability of P int in the model corresponds to the different internal pressure values which could be applied to the shell by the user. To measure the contribution of the correlated inputs to the output variability, we proceed to two estimation methods:

(1) As the FoBaHOGS technique gives very similar results with the LHOGS one, we perform the FoBaHOGS strategy with n = 1000 simulations for each of 50 replications. We use the 5-spline functions for the geometrical parameters and the Hermite basis functions of degree 7 for the physical parameters. Figure 5 displays the first order sensitivity indices S i , i = 1, • • • , 8, and their splits into VS i and CoVS i for an easier interpretation. (2) As we are faced to four independent groups of dependent variables, {R int , T shell , T cap }, {E cap , σ y,cap }, {E shell , σ y,shell } and {P int }, we estimate the usual Sobol indices of these four groups by the Monte Carlo procedure [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF], with two observations samples of respective size n = 10000, for each of 50 replications. The aim is to compare the results obtained by this approach to our procedure, and to illustrate the additional information one can get with the generalized sensitivity indices. Figure 4 plots the contribution of the independent groups into the model variability. Interpretation. Through the Case 1, we observe that MUC clearly underestimates the sensitivity indices even if n is large, whereas our strategy behaves well. Also, the MUC procedure gives a very large variability for S 1 and S 2 that does not appear in the LSE-HOGS one. Moreover, trough Figures 1(a)-2(a), one can observe that the variability of the estimation gets very small when n gets large, which numerically strengthens the convergence result of Section 4.2. In terms of numerical comparison, the CPU time is not significantly different from one technique to another, whereas the squared error is much smaller for the LSEHOGS. In Case 2, the POM shows its limitation, as only interaction terms involved in independent pairs of dependent inputs can be estimated. The MUC strategy behaves well in this situation, although the mean-squared error is the largest of the four methods. Nevertheless, in view of Case 1, MUC is not a robust technique, and the quality of the estimation might vary according to the model complexity. At last, (L/FoBa)HOGS give very similar results in terms of estimation. The LHOGS is much slower than FoBaHOGS because the LARS computes the whole regularization path { β(λ), λ ∈ R} before choos- ing one solution by selection criterion [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. Thus, although the LARS penalization offers a theoretical consistency, its time performance is low. The Case 3 compares the Jacques's method [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF] with our strategy. Through the usual approach, we notice that {R int , T shell , T cap } and {P int } are the most contributive groups in the model. However, we are unable to detect the influence of each input in this way.

Our strategy offers to do it. Moreover, we notice that the sensitivity index of any group in Figure 4 is equal to the sum of the individual contributions of Figure 5(a). This shows the relevance of our sensitivity indices. It also makes sense with the remark about the generalization made on page 6. Indeed, the model function can be decomposed as Y = η ∅ + η {Rint,Tshell,Tcap} + η {Ecap,σy,cap} + η {Eshell,σy,shell} + η {Pint} + interaction terms, where each summands are orthogonal. Thus, in this case, S Sobol {Rint,Tshell,Tcap} ≡ S {Rint,Tshell,Tcap} . Further, any summand can be decomposed as [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. For example, one can write, η {Rint,Tshell,Tcap} = η Rint + η Tshell + η Tcap + interaction terms. Thus, S Sobol {Rint,Tshell,Tcap} = S Rint + S Tshell + S Tcap . Moreover, the sensitivity can be explained thanks to the split of the indices, as showed in Figure 5(b)-5(c). The effects of the material parameters E shell , E cap , σ y,shell , σ y,cap are negligible, so we can conclude that they do not have any influence in the model.

The internal pressure P int has an influence on the model response, but the strongest contribution comes from the correlated set of geometrical inputs (R int , T shell , T cap ). More precisely, we deduce that T cap has an important full contribution, barely weakened by the contribution induced by the dependence. Thus, one can deduce that T cap is very significant in the model. In view of low full and induced contribution, we can reasonably deduce that R int has a very small influence in the model. The sensitivity index of T shell is quite small, but it should be noticed that the covariance part plays the role of compensation, and, if possible, one should work on the dependence with R int to increase its influence.

Conclusions and perspectives

This paper brings a new methodology to estimate the components of the generalized functional decomposition, when the latter satisfy hierarchical orthogonal constraints. Moreover, we show the consistency of the estimators when the usual least-squares estimation is used to estimate the unknown coefficients. From a practical point of view, it appears that the penalized regression should be often applied, and we observe that a selection variable strategy is numerically efficient. However, both FoBa and LARS suffer from limitations. On the one hand, FoBa is a very intuitive algorithm, but has no theoretical basis to guarantee the consistency of the estimators. On the other hand, LARS is computationally expensive when compared with a greedy strategy. It is also well-known that LARS suffers from numerical instability when predictors are strongly correlated. The future objective is to overcome these issues by exploring the numerical and theoretical properties of our methodology when the ℓ 0 penalization is relaxed by the L 2 -boosting [START_REF] Champion | Sparse regression and support recovery with L 2 -boosting algorithm[END_REF][START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF]. → 0, so we can conclude.

Let show (iv). We have,

k i=1 φlu i , φlv n - k i=1 φ lu i , φ lv ≤ k i=1 φlu i -k i=1 φ lu i , φlv n + k i=1 φ lu i , φlv n - k i=1 φ lu i , φ lv ≤ k i=1 φlu i -k i=1 φ lu i n + k i=1 φ lu i n φlv -φ lv n + k i=1 φ lu i , φ lv n - k i=1 φ lu i , φ lv .
By the strong law of large numbers,

k i=1 φ lu i , φ lv n - k i=1 φ lu i , φ lv a.s.
→ 0, and we can conclude with the previous arguments.

As (φ 1 ) L li=1 is an orthonormal system, we get

| φ 1 , 1 n | a.s. → E(φ 1 ) = 0 and T 1 n a.s. → φ 1 = 1. Therefore, φ1 -φ 1 a.s. → 0. Also, φ1 -φ 1 n ≤ 1-T 1 n T 1 n φ 1 n + | φ1,1 n| T 1 n .
Exactly with the same previous argument, we conclude that φ1φ 1 The resolution of (A4) leads to the resolution of a linear system, when removing C n lu , of the type

A u,n Λ u,n = D lu,n ,
where Λ u,n is the vector of unknown parameters (λ n lv,lu ) lv∈[1:L] |v| ,v⊂u , A u,n is the matrix whose block entries are ( φv1 , φv2 n ) v1,v2⊂u , and D lu,n involves block entries (-⊗ k i=1 φlu i , φvi n ) vi⊂u .

Also, the theoretical construction of the functions (φ lu ) lu consists in setting The resolution of (A5) leads to the resolution of a linear system, when removing C lu , of the type

φ lu = k i=1 φ lu i + v⊂u v =∅ lv∈[1:L] |v|
A u φ Λ u = D lu ,
where Λ u is the vector of unknown parameters (λ lv,lu ) lv∈[1:L] |v| ,v⊂u , A u φ is the matrix whose block entries are ( φ v1 , φ v2 ) v1,v2⊂u , and D lu involves block entries (-⊗ k i=1 φ lu i , φ vi ) vi⊂u .

We → 0, ∀ j, ∀ l uj .

Thus, we conclude that (A8) is satisfied.

Secondly, as k i=1 φlu i -k i=1 φ lu i n a.s.

→ 0, Assertion (iv) of Lemma A.3 can be applied. Assertion (iii) claims that A u,n tends to the theoretical matrix A u φ . Also, by (iv) of Lemma A.3, D lu,n a.s.

→ D lu . Hence, Λ u,n a.s. → Λ u . We also deduce that C n lu a.s.

→ C lu .

  which is the usual definition of Sobol indices for independent inputs. Proposition 2.3 Under (C.1) and (C.2), the sensitivity indices S u previously defined sum to 1, i.e. u∈C\{∅} S u = 1.

  (a) the model is evaluated n = 200 times over 50 runs. The estimated first and second order sensitivity indices are represented in Figure 1 when LSEHOGS and MUC methods are used. (b) the model is evaluated n = 1500 times over 50 runs. The estimated first and second order sensitivity indices are represented in Figure 2 in this condition.

Figure 1 .Figure 2 .

 12 Figure 1. Test case 1(a). Sensitivity indices estimation with n = 200 evaluations

Figure 3 .

 3 Figure 3. Tank distortion at point y Table 4. Description of inputs of the shell model Inputs Meaning Distribution R int shell internal radius U ([1800; 2200]), γ(R int , T shell ) = 0.85 T shell shell thickness U ([360; 440]), γ(T shell , Tcap) = 0.3 Tcap cap thickness U ([180; 220]), γ(Tcap, R int ) = 0.3 Ecap cap Young's modulus αN (µ, Σ) + (1α)N (µ, Ω), α = 0.02 σy,cap cap yield strength µ = 210 500 , Σ = 350 0 0 29 , Ω = 175 81 81 417 E shell shell Young's modulus αN (µ, Σ) + (1α)N (µ, Ω), α = 0.02 σ y,shell shell yield strength µ = 70 300 , Σ = 117 0 0 500 , Ω = 58 37 37 250 P int internal pressure N (80, 10)

Figure 4 .

 4 Figure 4. Sobol indices of the four groups of variables

Figure 5 .

 5 Figure 5. Test case 3. First order sensitivity indices with n = 1000

1 )Tφ. 2 n 1 / 2 (

 1212 is true.• Let l i ∈ [1 : L]. Suppose that (H ′ k ) is true for any k ≤ l i . Let us show (H ′ li+1 ) holds. By construction, we get, φli+1 = φ li+1 -li k=0 φ li+1 , φk n • φk φ li+1 -li k=0 φ li+1 , φk n • φk n li+1 , φk n • φk . For all k ≤ l i , φ li+1 , φk n = φ li+1 , φ k n + φ li+1 , φkφ k n ≤ | φ li+1 , φ k n | + φ li+1 , φkφ k n ≤ | φ li+1 , φ k n | + φ li+1 n • φkφ k n Byinduction, φkφ k a.s. → 0. By the usual law of large numbers, | φ li+1 , φ k n | a.s. → | φ li+1 , φ k | = 0 as the system (φ li ) L li=1 is orthonormal. As φ li+1 n a.s. → φ li+1 = 1, we deduce that φ li+1 , φk n a.s. → 0. And, φ li+1 , φk n φk ≤ φ li+1 , φk

λ

  lv,lu φ lv + C lu , where (C lu , (λ lv,lu ) lv,v⊂u ) are computed by the resolution of the following equationsφ lu , φ lv = 0, ∀v ⊂ u, ∀ l v ∈ [1 : L] |v| φ lu , 1 = 0. (A5)

a

  i1 • • • a is • b j1 • • • b jt where a i = φ lu i /T lu i n ,and b j = φ j lu j , φk n • φk /T lu j n . As already proved, for all i ∈ [1 : p], l ui ∈ [1 : φ lu j , φk n φk n a.s.

  and let C be the collection of all subsets of [1 : p]. As misuse of notation, we will denote the sets {i} by i, and {ij} by ij. For u ∈ C with u = {u 1 , • • • , u k }, we denote the cardinality of u by |u| = k and the corresponding random subvector by X

Table 1 .

 1 Test case 1. Numerical comparisons of the MUC and the LSEHOGS CPU time (in sec.) mse(η) × 10 -3

	Case 1(a)	MUC	0.1628	28.2
		LSEHOGS	0.35	6.3
	Case 1(b)	MUC	1.867	20.8
		LSEHOGS	0.658	2

Table 2 .

 2 Test case 2. Sensitivity indices estimation (with standard deviations) with n = 300

		S 1	S 2	S 3	S 12	S 13	S 23	S 123
	Analytical	0.4683	0.4652	0.0593	0.0072	0	0	0
	POM	0.4402 (0.021)	0.4718 (0.0401)	0.0810 (0.0012)	-0.0014 (0.001)	--	--	-
	FoBaHOGS	0.4488 (0.0216)	0.4699 (0.0190)	0.0714 (0.0233)	0.0041 (0.0028)	0 (0)	0 (0)	0.0058
	LHOGS	0.4536 (0.0216)	0.4733 (0.0193)	0.0745 (0.0227)	0.0065 (0.0017)	0.0013 (0.0017)	0.0006 (0.0009)	-0.0098
	MUC	0.4439 (0.0206)	0.4533 (0.0185)	0.0713 (0.0221)	0.0002 (0.0001)	0.0002 (0.0004)	0.0001 (0.0003)	0.0310

Table 3 .

 3 Test case 2. Numerical comparisons of the MUC, the POM and the FoBa/LHOGS with n = 300

		j 1( βj = 0) CPU time (in sec.) mse(η) × 10 -3
	POM	-	11.96	6.5
	MUC	-	11.97	7.9
	FoBaHOGS	5.9	3.38	5.1
	LHOGS	152.6	77.23	5.1

  Proof of Lemma A.[START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF].The first point (i) is trivial. Now, we have, by (H.2),φ lu , φlvφ lu , φ lv = φ lu , φlvφ lv lu , φlv nφ lu , φ lv ≤ φ lu , φlvφ lv n + | φ lu , φ lv nφ lu , φ lv | = φ lu n φlvφ lv n + | φ lu , φ lv nφ lu , φ lv | .By the usual strong law of large numbers, | φ lu , φ lv nφ lu , φ lv | φlvφ lu , φ lv = φluφ lu , φlv + φ lu , φlvφ lv ≤ φluφ lu φlv + φ lu φlvφ lv . φ lv = 1 (by (i)). Further, φlu , φlv nφ lu , φ lv ≤ φlu , φlvφ lv n + φlu , φ lv nφ lu , φ lv lv nφ lu , φ lv .

						≤ φ lu	φlv -φ lv	a.s. → 0.
						=1
	Further,
	φ a.s. → 0. By hypothesis, φlv -φ lv n a.s. → 0, and φ lu n a.s. → 1 by (i). Hence, (ii) holds.
	The point (iii) follows from
	φlu , By assumptions, φlu -φ lu satisfied, as φlv a.s.	a.s. → 0 and φlv -φ lv	a.s. → 0. Thus, the first point of (iii) is
						= φlu
	First, φlu , φ lv	φlu	n	=	1. By hypothesis,	φlv -φ lv	n	a.s. →	0. By (ii),

→ n φlvφ lv n + φlu , φ nφ lu , φ lv a.s.

  With the previous arguments, φ li+1 , φk n a.s. 1. Then, we conclude that (H ′ li+1 ) is true. Therefore, (H 1 ) is satisfied.• Let k ∈ [1 : p]. Suppose now that (H |ũ| ) is true for any 1 ≤ |ũ| ≤ k -1, and any l ũ ∈ [1 : L] |u| . Show that (H k ) is satisfied. Let u be such that u = {u 1 , • • • , u k }.First, as (H |ũ| ) is true for any 1 ≤ |ũ| ≤ k -1, results (i)-(ii)-(iii) of Lemma A.3 can be applied to any couple (φ lu , φlu ) such that |u| ≤ k -1.Further, we have seen that, for any l u ∈ [1 : L] |u| , (λ n lv,lu ) lv,v⊂u ) are computed by the resolution of the following equations φlu , φlv n = 0, ∀v ⊂ u, ∀ l v ∈ [1 : L] |v| φlu , 1 n = 0.

			n → 0, and T li+1	a.s.
	k
	φu lu =	φlu i +	λ n lv,lu φlv + C n lu ,
	i=1	v⊂u v =∅	lv∈[1:L] |v|
	where (C n lu , (A4)
			1 n T li+1	φ li+1 , φk n .

→

  have, φluφ lu = k i=1 φlu i -k i=1 φ lu i + v⊂u v =∅ lv∈[1:L] |v| (λ n lv,lu φlvλ lv,lu φ lv ) +C n lu -C lu ≤ k i=1 φlu i -k i=1 φ lu i + v⊂u v =∅ lv∈[1:L] |v| λ n lv,lu φlvλ lv,lu φ lv + C n lu -C lu , |v| λ n lv,lu φlvλ lv,lu φ lvEach univariate function φlu i is constructed from (φ k ) k≤lu i by the Gram-Schmidt procedure. Thus,

							(A6)
	and,
	φlu -φ lu	n	≤	k i=1 φlu i -k i=1 φ lu i n + C n lu -C lu .	+ v⊂u v =∅	lv∈[1:L] n (A7)
	First, we show that
				 	k i=1 φlu i -k i=1 φ lu i	a.s.	(A8)
				
	k i=1 φlu i -k i=1 φ lu i =	k i=1	φl u i -	lu i k=0 -1 T lu φl u i , φk n• φk i n	-k i=1 φ lu i
					≤	k i=1 φ lu i	k i=1	1 lu i n T	-1
					+	s+t=k (s,t)∈N×N *	lu j 1 -1 k=0
						1≤i1<•••<is≤k
						1≤j1<•••<jt≤k

→ 0, ∀ l ui ∈ [1 : L] k i=1 φlu i -k i=1 φ lu i n a.s. → 0, ∀ l ui ∈ [1 : L].

Acknowledgments

The authors would like to thank Loic Le Gratiet and Gilles Defaux for their help in providing and understanding the tank pressure model. This work has been supported by French National Research Agency (ANR) through the COSINUS program (project COSTA-BRAVA number ANR-09-COSI-015).

Appendix A. Convergence results

In this part, we restate and prove Proposition 4.1 of Section 4.2. For sake of clarity, we first define and recall some notation that will be used further.

Reminder

First, as mentioned in Section 4, we assume that Y is centered. Recall that, ∀ i ∈ [1 : p], L is the dimension of the spaces H 0,L i and Ĥ0,L i . Also, dim(H 0,L u ) = dim( Ĥ0,L

We refer (φ lu ) lu∈[1:L] |u| as the basis of H 0,L u and ( φlu ) lu∈[1:L] |u| as the basis of Ĥ0,L u constructed according to HOGS Procedure of Section 3. Thus, these functions all lie in L 2 R (R p , B(R p ), P X ).

•, • and • are used as the inner product and norm on L 2 R (R p , B(R p ), P X ),

• n and • n denote the empirical inner product and norm, that is

is the n-sample of observations from the distribution of (Y, X).

New settings

We set m := u∈C L |u| the number of parameters in the regression model. Denote, for all u ∈ C, Φ u (X u ) ∈ (L 2 (R, B(R), P X )) L |u| , with (Φ u (X u )) lu = φ lu (X u ), and by β any vector of Θ ⊂ R m , where (β

Recall that, for a, b ∈ N * , M a,b (R) denotes the set of all real matrices with a rows and b columns.

, where ( φu ) s,lu = φlu (x u s ), and we set

Denote by A φ be the m × m Gram matrix whose block entries are

The main convergence result is reminded further below.

Further, let us consider the least-squares estimation ηL of η L using the sample (y s , x s ) s∈[1:n] from the distribution of (Y, X), and the functions ( φlu ) lu , that is

, and Θ is a compact set of R m . If we assume that (H.1) The distribution P X is equivalent to

The proof of Proposition is broken up into Lemmas A.2-A.5. To prove (A1), we introduce ηL as the following approximation of η

and we write the triangular inequality,

Thus, it is enough to prove that ηL -ηL a.s. → 0, and that ηLη L a.s. → 0. Lemmas A.4 and A.5 deal with convergence results on ηLη L and on ηL -ηL , respectively. Lemmas A.2, A.3 are preliminary results to prove Lemmas A.4 and A.5.

A.1 Preliminary results

Lemma

First of all notice that when we consider a Gram matrix, by a classical argument on the associated quadratic form, the full rank of this matrix holds if and only if the associated functional vector has full rank in L 2 [START_REF] Holmes | On random correlation matrices[END_REF].

To begin with, set, for all i ∈ [1 : p],

As (φ li ) L li=1 is an orthonormal system, we obviously get G i = I (L+1)×(L+1) , where I denotes the identity matrix. Now we may rewrite the tensor product ⊗ p i=1 G i as

We obviously have ⊗ p i=1 G i = I. So that, using the remark of the beginning of the proof, the system

is linearly independent (⊗ p i=1 P Xi )-a.e.

As we assumed that ⊗ p i=1 P Xi and P X are equivalent by (H.1), we get that 1 (⊗ i∈u Φ i ) u⊆[1:p] u =∅ is linearly independent P Xa.e.. Now, we may conclude as in the classical Gram-Schmidt construction. Indeed, the construction of the system (Φ u ) u∈C involves an invertible triangular matrix.

A.2 Main convergence results

Lemma A.4 Remind that the true regression function is

Further, let ηL be the approximation of η L ,

Then, under (H.2)-(H.3), we have

→ 0 ∀ u ∈ C, and ηLη L a.s. → 0.

Proof of Lemma A.4.

For any u ∈ C,

Let us show that φluφ lu a.s.

→ 0. Actually, the proof of this convergence requires the use of Lemma A.3, so we also have to show that φluφ lu n a.s.

→ 0. These two results are going to be proved by a double induction on |u| ≥ 1 and on l u ∈ [1 : L] |u| . We set

Let us show that (H k ) is true for any k ≤ p :

• Let u = {i}, so k = 1. We used the Gram-Schimdt procedure on (φ li ) L li=1 to construct ( φli ) L li=1 . Let us show by induction on l i that φliφ li a.s.

Consequently, by induction, we deduce that every piece of the right-hand side of (A6) (respectively (A7)) tends to 0, so is φlu

Proof of Lemma A.5. First, we remind the true regression model,

with

, and β 0 = (β u,0 lu ) lu,u∈C the true parameter. Let

Due to Lemma A.2, ( t X φ X φ ) -1 is well defined. Thus,

By the law of large numbers, ( [START_REF] Golub | Matrix computations[END_REF]. We use this inequality to get

First, by (A11),

After computation, we get → A φ by (ii) of Lemma A.3.

Under (H.1), using the result of Lemma A.2, A φ is invertible. Then, n ( t X φ X φ ) -1t X φ -( t X φX φ) -1t X φ → 0.

Finally, as a consequence of Lemmas A.4-A.5, ηLη L ≤ ηL -ηL + ηLη L a.s. → 0.