Rachele Allena
email: rachele.allena@ensam.eu

Christophe Cluzel

Identication of anisotropic tensile strength of cortical bone using Brazilian test

Keywords: Cortical bone, Anisotropy, Brazilian test, Brittle strength

For a proper analysis of cortical bone behaviour, it is essential to take into account both the elastic stiness and the failure criteria. While ultrasound methods allow complete identication of the elastic orthotropic coecients, tests used to characterise the various failure mechanisms and to identify the brittle tensile strength in all directions are currently inadequate. In the present work we propose the Brazilian test as a complement to conventional tensile tests. In fact, this experimental technique, rarely employed in the biomechanics eld, has the potential to provide an accurate description of the anisotropic strength of cortical bone. Additionally, it allows to assess the scale inuence on failure behaviour which may be attributed to an intrinsic length in correlation with the cortical bone microstructure. In order

to correctly set up the Brazilian test, several aspects such as the machining, the geometrical parameters of the specimen and the loading conditions were determined. The nite element method was used to evaluate the maximal 1. Introduction 1.1. Bone's structure and behaviour Bone presents a hierarchical structure [START_REF] Currey | Bones: structure and mechanics[END_REF] [START_REF] Rho | Mechanical properties and the hierarchical structure of bone[END_REF]) [START_REF] Vashishth | Hierarchy of bone microdamage at multiple length scales[END_REF] which is organised in dierent levels as follows: i) the macrostructure: cancellous and cortical bone, ii) the mesostructure (from 10 to 500 µm): haversian system, osteons, trabeculae and iii) the microstructure (1-10 µm): the lamellae and the osteocytes iv) the nanostructure (from a few hundred nanometers to 1 µm): brillar collagen and embedded mineral v) the sub-nanostructure (below a few hundred nanometers): collagen, molecules and proteins.

The complex structure of the bone has been the object of many studies during the last decades in order to decipher the inuence of each level on both the mechanical and failure behaviour [START_REF] Currey | Bones: structure and mechanics[END_REF] . At the nanoscale, the orientation of the collagen brils and their degree of mineralisation (Turner-Walker and Parry, 1995) may aect the Young's modulus leading to a failure stress in the bres direction. At the microscale, the stacking of successive lamellae, each composed by collagen bres oriented in a single direction, provides an isotropic mechanical behaviour in the lamellae plan, while weak properties are observed along the perpendicular direction. At the mesoscale, the osteons structure supplies a transverse isotropy for both the stiness and the failure stresses [START_REF] Rho | Mechanical properties and the hierarchical structure of bone[END_REF] [START_REF] Ascenzi | Individual-specic multi-scale nite element simulation of cortical bone of human proximal femur[END_REF]. Such a behaviour is maintained at the macroscale due to the main orientation of the osteons along the longitudinal axis of the bone. Finally, at this level the interface between the osteons and the interstitial lamellae (the cement line) brings a further weakness to the failure behaviour.

Conventional mechanical tests in traction, compression and torsion on specimens obtained from cortical bone of the femur diaphysis were carried out by [START_REF] Reilly | The elastic and ultimate properties of compact bone tissue[END_REF]. They actually observed that the Young's modulus along the longitudinal direction is double that measured along the circumferential or radial directions. Therefore, the anisotropy of the elastic behaviour is clearly marked and complies with the geometrical organisation of the bone at the mesoscopic scale. Several experimental tests may be used to evaluate the strength for a brittle and anisotropic material like bone. Tensile testing is one of the classical methods to measure bone's mechanical properties. Nevertheless, specimens must have relatively large dimensions (15-20 mm in length, 4-8 mm in width) and they must be specically designed to obtain the majority of the strain in the central region [START_REF] Reilly | The elastic modulus for bone[END_REF] [START_REF] Ashman | Elastic properties of cancellous bone: Measurement by an ultrasonic technique[END_REF].

If one assumes that the external force is applied without inducing a bending moment, the tensile test provides a good assessment of bone's strength, but is limited in its ability to evaluate the eects of anisotropy due to the constraints on the dimensions of specimens.

Bending tests are usually employed for testing the bones of small animals, for which a tensile test is dicult to set up. In such a test, the entire bone is loaded until failure leading to tensile stresses on one side of the bone and compressive stresses on the other side. Additionally, tensile or compressive stresses increase from the neutral axis to the external boundaries of the specimen. Thus, failure commonly occurs on the tensile side since bone is weaker in tension than in compression [START_REF] Reilly | The elastic and ultimate properties of compact bone tissue[END_REF]) and may also be highly sensitive to surface defects due to the machining of the specimens for instance. Bending may be applied to the bone through either a three-point or a four-point loading. The former is very simple to set up, but it may cause high shear stress around the middle section of the bone. The latter induces pure bending and ensures zero transverse shear stress between the two upper loading points. Nevertheless, if the specimen is rather small in length and the bending moment is maximum under the loading point, the stress state is not easy to determine. Furthermore, in both three-point and four-point bending tests the total length of the specimen should be about sixteen times the thickness of the specimen to guarantee that 85-90 % of the bone exion is actually due to bending. Unfortunately, this length-width ratio cannot be acquired in whole bones such as femora or tibias.

For compression testing, relatively small specimens (7-10 mm long) can be used and therefore machined along the three directions, but the measurement tends to be less accurate than those for tensile tests because of edge eects. In those regions in fact, the strain is likely to be higher than in the central region, possibly due to the misalignment of the specimen faces or other problems associated with specimen machining. Then, because of friction between the contact surfaces of the bone specimen and the plates of the testing machine, one may have a unidirectional strain at the boundaries and a stress static state in the central region, such that the specimen acquires a barrel-like shape. Although an extensometer is usually employed during tensile tests to determine the axial strain in the specimen, this is not possible in compression due to the small dimensions of the specimen. In this case image correlation represents an alternative method to evaluate the stress-strain relationship.

Despite a lower accuracy of the results compared to tensile tests, compressive testing presents some major advantages. First, specimens do not have to be as large as tensile specimens. Second, machining of compressive specimens is easier than for tensile specimens and may be done in dierent directions to investigate the anisotropic behaviour of the bone. Nevertheless, compression tests do not initiate the same failure modes as tensile tests (for which failure mode and crack shape show a specic brittle mechanism).

In recent years, shear tests have been developed to determine the shear modulus of elasticity of the bone. Among them we mention the rail shear test, the torsion tube, cross-beam specimens and tension-compression of notched specimens, including the Iosipescu (ASTM D5379) [START_REF] Iosipescu | New accurate procedure for single shear testing of metals[END_REF] [START_REF] Funk | Eect of cement modulus on the shear properties of the bone-cement interface[END_REF] [START_REF] Sharma | Studies on locational variation of shear properties in cortical bone with iosipescu shear test[END_REF] and the Arcan tests [START_REF] Arcan | A method to produce uniform plane-stress states with applications to ber-reinforced materials[END_REF].

Although the previous resistance tests allow partial assessment of the anisotropic characteristics of cortical bone's behaviour and identication of some fracture modes, they fail in evaluating the anisotropy in traction. For this reason, here we propose the Brazilian test as an alternative experimental approach to characterise the bone failure responses along the longitudinal, circumferential and radial axes. Such a test presents interesting features, which appear to be decidedly appropriate to study bone's mechanical be-haviour and to obtain a complete predictive model. In the Brazilian test, a cylindrical specimen is loaded in compression until failure over a short strip along the specimen length at each end of the vertical diameter. Compression induces tensile stresses normal to the loading direction, which are approximately constant within a region around the centre. Therefore, for a brittle material, a crack appears perpendicular to the maximum traction stress direction, leading to the splitting of the cylinder into two halves.

The Brazilian test has some interesting characteristics. Firstly, it greatly simplies the traction loading of a brittle material. Secondly, it permits reduction of the size of the specimen down to that limited by testing a representative volume of the material. For the specic case of cortical bone, such a reduction in dimensions (e.g. some millimetres in diameter) leads to three further benets: i) it decreases the probability of nding very large defects that may induce macroscopic rupture, ii) it provides information on the correlation between specimen size and defect distribution and iii) it enables the analysis of the traction fracture along the three main axes of the bone. Therefore, the Brazilian test may be employed to provide an accurate identication of the anisotropic maximal traction stresses in cortical bone.

Materials and methods

Sample preparation

Specimens were obtained from a bovine tibia sourced from a local butcher and conserved at -18 • . Once the tibia was defrosted, the internal marrow and spongy bone were removed and the bone was cleaned with water. The three main local axes of the bone were chosen as follows (Fig. 2):

• the longitudinal axis x 1 corresponds to the main direction of the tibia;

• the circumferential axis x 2 coincides with the azimuthal direction;

• the radial axis x 3 is aligned with the outward radius of the bone's section.

First, 25 bone cylinders were machined using diamond-tipped tubular drills of internal diameters φ 10, 8, 6 and 4 mm. For the sake of convenience, the machining was performed along the x 1 and x 3 directions, which maintains the ability to obtain the three fracture stresses σ f 11 , σ f 22 and σ f 33 ((Fig. 2), the superscript f indicates failure). Second, the cylinders were sectioned perpendicular to the cylinder axis using a diamond disc saw. Furthermore, for those machined along the x 1 direction, more than one specimen was obtained.

Finally, 29 specimens were acquired. The length L of the samples was set to 6.5 mm, 5.2 mm, 3.9 mm and 2.6 mm respectively for φ = 10 mm, φ = 8 mm, φ = 6 mm and φ = 4 mm (Fig. 3). Such values provide a minimal average ratio φ/L equal to 1.54. Before sectioning, the three main axes x 1 , x 2 and

x 3 were identied on each specimen which allows to classify the specimens as follows: x i _F j , with x i and F j indicating the cutting axis and the loading direction, respectively (Fig. 2). During cutting, water was used in order to reduce both friction and temperature rise.

Brazilian test for cortical bone

The Brazilian tests were performed at room temperature right after the cutting, using a universal traction-compression machine INSTRON 5500-R equipped with a 5 KN sensor. We have assumed that the room humidity does not inuence the specimens behaviour. The machine was controlled by xing the displacement rate of the upper plate at 0.2 mm/min. The positioning of the specimen between the two plates of the machine as shown in Fig. 4 must During a regular test, the crack was generally initiated at the centre of the cylinder along the vertical axis (Fig. 5). Nevertheless, abnormal splitting might be observed due to i) shear stress (Fig. 6a), ii) crushing issues (Fig. 6b) or iii) a non centred crack. Defects such as those presented in Fig. 6a were mainly found during a preliminary series of tests with specimens having a ratio φ/L < 1.54. To limit the crushing of the contact surface (Fig. 6b), a cushion can be inserted between the specimen and each load plate as described in the standard for Brazilian tests applied to rocks (ISRM 1978, ASTM 2008). In our case, a 0.52 mm thick square of cardboard was used (Fig. 4). The imprint was measured after each test to estimate the contact area and we found that it can actually be dened independently of the specimen diameter φ through the angle α (Fig. 4a) as described in Wang et al. (2004).

Structural analysis of Brazilian test for anisotropic elastic behaviour

Through a structural analysis, we may be able to evaluate the maximal tensile stress σ xx,max along the x direction at the centre of each specimen.

For isotropic materials, an analytical solution was proposed by [START_REF] Peltier | Theoretical investigation of the brazilian test[END_REF] giving the tensile stress in the centre of the disc as follows

σ xx,max = 2F φLπ ( 1
)
where F is the applied load.

To account for the eect of a soft cushion between the specimen and the loading plates, a factor of correction k was introduced by Hondros (1959)

and [START_REF] Wang | The attened brazilian disc specimen used for testing elastic modulus, tensile strength and fracture toughness of brittle rocks: analytical and numerical results[END_REF] as a function of the angle α (Fig. 7). Thus, Eq. [1] becomes : overview of anisotropy eects, the denition of the maximum tensile stresses (Fig. 2) is similar to that proposed in Eq. [2], but the coecient of correction is now expressed as a function of both the direction of the failure stress and of the loading. Thus, we have:

σ ii,max = β ii-j 2F j φLπ (3)
where β ii-j is the correction factor and F j is the applied vertical load.

The subscripts ii and j indicate the principal stresses and loading direction, respectively. The main objective of the structural analysis is to nd the coecient β ii-j for the dierent directions independently of the specimen diameter φ.

In the present study, the analysis was performed using the nite elements The problem was solved using the plane stress hypothesis. The displacement of the point A was constrained along the x direction while the point B was totally constrained to prevent rigid body motion (Fig. 7). Finally, the vertical load was applied along the upper and lower boundaries l (blue lines in Fig. 7), which, as previously explained (Sec. 2.2), have been calculated

3. Results

Stress state in the loaded specimen

In this section we present the numerical results and in particular we discuss the stress eld inside the specimen. As it is possible to observe in Fig. 8a and 8b, for a load per length unit F = 1400 N/mm (which is the same for each tested diameter), the compressive (σ yy ) and tensile (σ xx ) stresses are heterogeneously distributed. Their pattern is very similar to that of the isotropic case as reported in Wang et al. ( 2004) and specically σ yy and σ xx are maximal along the loading surfaces and at the centre, respectively.

Actually, there exists a relationship between such stresses and the failure mechanism. In fact, as shown in Fig. 5, the crack is distinctly open at the centre of the disc ((x,y) = (0,0)) where the stress state is plane and given by

σ =      σ xx 0 0 0 σ yy 0 0 0 0      (4)
with σ xx = 55 MPa and σ yy = -147 MPa (blue line in Fig. 8d and8b, respectively).

Let n and t be respectively the normal and the tangent vectors to the failure plane dened as where θ is the angle between n and the x axis. Then, the normal (σ n ) and the shear (τ ) stresses read

n =      cosθ sinθ 0      and t =      -sinθ cosθ 0      (
σ n = n tr σ n (6) τ = t tr σ n (7)
with n tr the transposition of n.

It is interesting to evaluate the evolution of σ n and τ for i) θ varying between 0 • and 180 • and ii) the axial coordinate x of the point of interest (x,y) varying between ±0.73 mm from the centre of the disc (Fig. 9). For θ = 0 • , we nd σ n = 55 MPa and no shear stress, while τ is maximal (± 100 MPa) for θ = 45 • and 135 • . Finally, for θ = 90 • , σ n is equal to -147 MPa showing a compressive stress state (Fig. 9). It can be noticed that for all these stresses, the maximal values are found at the centre of the disc (blue line in Fig. 8b,8d,9). For a brittle material, the failure plane is a useful parameter to evaluate the cracking mechanism and the corresponding stress.

0°45°90°135°180°-150

Here, failure is not activated at θ = 90 • nor at θ = 45 • . On the contrary, the traction stress σ xx is assumed to be responsible for the failure each time the crack occurs parallel to the loading axis.

The main objective of the numerical simulations was to evaluate the correction factor β ii-j dened in Eq. [START_REF] Ashman | Elastic properties of cancellous bone: Measurement by an ultrasonic technique[END_REF], which is independent of the diameter φ of the disc. For an isotropic material, we found that such a coecient is equal to 1 in the case of a concentrated load F j and to 0.912 in the case of a distributed load as described in Sec. 2.3, which is very close to 0.92, the coecient analytically calculated from [START_REF] Wang | The attened brazilian disc specimen used for testing elastic modulus, tensile strength and fracture toughness of brittle rocks: analytical and numerical results[END_REF].

In order to use β ii-j as a consistent indicator, the variation of the stress state must be low with respect to the cracking position. In Fig. 8b and8d σ yy and σ xx respectively are plotted for a plane placed at x = 0, 0.24, 0.48 and 0.73 mm for a disc with a diameter of 6mm. We notice that if the crack occurs between ±0.4 mm from the vertical axis of the disc, the maximum stress only varies by about ±5.5 %. To keep such a low variability, the corresponding spatial tolerances for φ = 4, 8 and 10 mm are ±0.27, ±0.4 and ±0.67 mm, respectively. As an example, in Fig. 5, the diameter φ of the specimen is equal to 6 mm and the position of the crack is at 0.16 mm from the centre with an error of -1 % for the coecient β ii-j .

Finally, as mentioned in Sec. 2.4, a source of uncertainty for the correction factor β ii-j is related to the variations of the elastic coecients. According to the sensitivity analysis that has been carried out, the results for the four loading cases are reported in Table 2.

Specimen σ f 11 σ f 22 σ f 33 x 1 _F 2 - - β 33-2 = 1.007± 4% x 1 _F 3 - β 22-3 = 1.007 ± 4% - x 2 _F 1 - - β 33-1 = 0.802 ± 5%
x 2 _F 3 β 11-3 = 1.044 ± 3.5% --

x 3 _F 1 - β 22-1 = 0.802 ± 5% - x 3 _F 2 β 11-2 = 1.045 ±3.5% - -
Table 2: Results of the sensitivity analysis and values of correction factors β ii-j for the four loading cases.

Experimental data

The experimental tests were exploited to assess the failure force as well as the crack direction, which must be vertical, and shape, which must be sharp-cut. Furthermore, by using the correction coecients β ii-j derived from the numerical analysis (Table 2), the values of the tensile failure stress were determined depending on the specimen diameter φ for each direction of failure stress tested σ f 11 , σ f 22 and σ f 33 (Fig. 10). Among the 29 tests carried out, 4 were stopped due to a crushing problem on the loading area (Sec. 2.2) and 4 presented a cracking mechanism outside the admissible region (Sec. 3.1) (hollow arrows in Fig. 10). For these specic cases, stress leading to failure was not usable as a value to rupture, but as an underestimation of the failure stress. As the traction stress σ xx is not homogeneous within the sample (Fig. 8d), it may be of interest to identify a failure region for each specimen diameter rather than simply determining the relationship between the failure stress and the sample dimensions. Thus, a rectangular area S f ailure of height h f ailure

and width e f ailure can be dened for each diameter φ such that 0.9 σ xx,max < σ xx < σ xx,max . We notice that the dimensions and consequently the area of the failure region decrease with the specimen diameter (Table 3).

Discussion

The Brazilian test is suitable for brittle materials only, but the experimental validation of the failure mechanism is very easy to achieve because the crack must be unique and in a vertical plane as described in Tavallali and Vervoort (2010). Additionally, if the φ/L ratio is controlled and optimised, the rare faulty tests may be attributed to machining or positioning defects. In the present work, although the bovine cortical bone we tested seemed rather young with a marked microstructure, the experimental dispersion was quite reasonable and the anisotropy of brittle fracture clearly appeared leading to a ratio σ ii,max /σ ii,min of the order of 2.

For elastic isotropic materials, the fairly simple geometry of the specimen used for the Brazilian test allows the existence of analytical descriptions of the stress eld either for a concentrated or a distributed load. In this case, the analytical solution and our numerical simulation were in very good agreement. Specically, for a concentrated load, the correction coecient β ii-j dened in Eq. [START_REF] Ashman | Elastic properties of cancellous bone: Measurement by an ultrasonic technique[END_REF] is exactly equal to 1, while for a distributed load as described in Sec. 2.3, β ii-j is equal to 0.912.

For an anisotropic material such as cortical bone, the elastic coecients deduced from [START_REF] Bernard | Accurate measurement of cortical bone elasticity tensor with resonant ultrasound spectroscopy[END_REF] were used to run the numerical simulations for specimens of dierent diameters. We were able to determine the correction factors β ii-j associated to each failure stress and we found that all the coecients are between 0.802 and 1.05 or in a range of 0.92 ±14 %. This results in a variation of the maximum stress of the order of 14 %. Furthermore, according to the sensitivity analysis we performed, the uncertainties on β ii-j due to the variation of the elastic parameters are not higher than 5 %, which is quite low. Therefore, the coecients can be directly used or, for better accuracy, recalculated after verication of the rigidity by, for example, an ultrasonic method.

The Brazilian test also allowed us to assess the scale inuence on failure mechanism. The areas of the failure regions for the dierent specimens reported in Table 3 are very small for a tensile test on a brittle material, which results in failure stresses for specimens with a diameter of 4 mm higher than those for larger samples (Fig. 10). Previous works have focused on this specic aspect and have used either a Weibull distribution of the defect size [START_REF] Fok | A numerical study on the application of the weibull theory to brittle materials[END_REF] or a cohesive crack model [START_REF] Guinea | Assessment of the tensile strength through size eect curves[END_REF] to describe such a behaviour. In both cases, the size eect is attributed to an intrinsic length in correlation with the microstructure of the material, below which the failure stress increases. This might also be the case for cortical bone. In fact, we can see that as the specimen diameter φ decreases, the dimensions h f ailure and e f ailure of the failure region decrease too (Table 3) and approach the dimensions of a portion of the cement line (Sec. According to the previous remarks, it would be interesting to perform the Brazilian test on a large number of specimens within a range of small dimensions. In fact, this would allow to consistently investigate the scale inuence and the statistical dispersion and to characterise a suitable nonlocal model to be adopted for numerical simulations.

Figure 1 :

 1 Figure 1: Bovine bone microstructure: sections perpendicular to the longitudinal axis

1. 3 .

 3 Brazilian test for brittle materials The Brazilian test was rst introduced by Carneiro (1943) and Akazawa (1943) to determine the tensile strength of brittle materials such as rock, concrete or ceramic, which is dicult to evaluate by performing a direct uniaxial tensile test. It is widely used in the eld of civil engineering and has been the object of numerous works for both the calculation of stresses and the identication of material properties (Li et al., 2013). In the biomechanics eld, it has been employed to determine the tensile strength of archeological cortical bone (Turner-Walker and Parry, 1995) and articially aged bone (Turner-Walker, 2011). Additionally, (Huang et al., 2012) proposed a numerical analysis of the Brazilian test of heterogeneous specimens in order to analyse the tensile strength of dental amalgams.

Figure 2 :

 2 Figure 2: Coordinates system, traction stresses (σ f ii ) and loading directions (F j ).

Figure 3 :Figure 4 :

 34 Figure 3: Specimens diameters: 4, 6, 8 and 10 mm.

Figure 5 : 2 )Figure 6 :

 526 Figure 5: Appearance of a vertical crack at the centre of a specimen x 1 _F 3 with a 8.12 mm diameter.

(

  FE) method, which provides a better validation and simplify the management of various input and output data. The FE software COMSOL 3.5a was used to run two dimensional (2D) simulations and to evaluate the linear elastic stress eld within the samples along x 1 , x 2 and x 3 with two loading directions each. The cylindrical specimens were represented as circles with an anisotropic elastic behaviour. The elastic material parameters were deduced from Bernard et al. (2013) (

Figure 8 :

 8 Figure 8: (a, c) Plot of σ yy and σ xx , respectively for a x 1 _F 2 specimen of diameter 6mm. (b, d) Outline of σ yy and σ xx respectively along the vertical diameter (blue line, x = 0) and along vertical lines placed at x = 0.24 mm (red line), x = 0.48 mm (purple line) and x = 0.73 mm (green line).

Figure 9 :

 9 Figure 9: Normal and shear stress distribution at the centre and at 0.24 mm (red lines), 0.48 mm (purple lines) and 0.73 mm (green lines) from the centre along the x axis.

Figure 10 :

 10 Figure 10: Maximum tensile stress versus diameter for three directions of loading.

  1.1), which may constitute a weakness for failure behaviour as mentioned in (Norman and Wang, 1997) (O'Brien et al., 2007) (MFeerick et al., 2013).

  we have proposed the Brazilian test as an alternative technique to investigate both the anisotropic strength and failure mechanism of cortical bone. In fact, although this test has rarely been employed in the eld of bone biomechanics (Turner-Walker and Parry, 1995) (Turner-Walker, 2011) (Huang et al., 2012), it presents some interesting features. Firstly, it allows testing of brittle materials in traction through the use of a compressive load. Secondly, it allows to reduce the specimen dimensions down to those of the representative volume of the material. Then, for specic case of cortical bone it has been possible to assess the tensile failure along its three main axes and its anisotropy.

6. Bibliography

Table 1

 1

		).
	Modulus	Value	Reference or formula
	E 1	20.3 GPa (Bernard et al., 2013)
	E 2	12.8 GPa (Bernard et al., 2013)
	E 3	12.8 GPa (Bernard et al., 2013)
	G 12	6.38 GPa (Bernard et al., 2013)
	G 13	6.32 GPa (Bernard et al., 2013)
	G 23	6.38 GPa (Bernard et al., 2013)
	ν 12	0.421	(Bernard et al., 2013)
	ν 13	0.434	(Bernard et al., 2013)
	ν 23	0.348	(Bernard et al., 2013)
	ν 21	0.265	ν 12 E 2 E 1
	ν 31	0.273	ν 13 E 3 E 1
	ν 32	0.348	ν 23 E 3 E 2

Table 1 :

 1 Elastic parameters sourced and deduced from Bernard et al. (2013).

  Failure region area S f ailure = h • e (mm 2 ) 1.2 2.64 4.64 7.4

	Specimen diameter φ (mm)	4	6	8	10
	Failure region height h f ailure (mm)	1.5 2.2	2.9 3.7
	Failure region width e f ailure (mm)	0.8 1.2	1.6	2

Table 3 :

 3 Values of the failure region area according to the specimen diameter φ.

using the angle α = 14 • . Thus, l is equal to 2.4 mm, 1.9 mm, 1.4 mm and 1 mm for φ = 10 mm, φ = 8 mm φ = 6 mm and φ = 4 mm, respectively.

Sensitivity analysis

The correction factor β ii-j may change with respect to the elasticity coecients. Therefore, a sensitivity analysis was performed for each specimen by varying the Young's moduli and Poisson's ratios by ±10 % relative to the 'benchmark' values.


Magali Champion

Gaelle Chastaing

Sébastien Gadat

Clémentine Prieur

L 2 -Boosting for sensitivity analysis with dependent inputs

Keywords: L2-boosting, convergence, dependent variables, generalized ANOVA decomposition, sensitivity analysis

This paper is dedicated to the study of an estimator of the generalized Hoeffding decomposition. We build such an estimator using an empirical Gram-Schmidt approach and derive a consistency rate in a large dimensional setting. We then apply a greedy algorithm with these previous estimators to a sensitivity analysis.

We also establish the consistency of this L2-boosting under sparsity assumptions of the signal to be analyzed. The paper concludes with numerical experiments, that demonstrate the low computational cost of our method, as well as its efficiency on the standard benchmark of sensitivity analysis.

Introduction

In many scientific fields, it is desirable to extend a multivariate regression model as a specific sum of increasing dimension functions. Functional ANOVA decompositions and High Dimensional Representation Models (HDMR) [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF]Li et al., 2010) are well known expansions that make it possible to understand model behavior and to detect how inputs interact with each other.

When input variables are independent, Hoeffding establishes the uniqueness of the decomposition provided that the summands are mutually orthogonal [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF]. However, in practice, this assumption is sometimes difficult to justify, or may even be wrong (see Li and Rabitz (2010) for an application to correlated ionosonde data, or [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF], who studied an adjusted neutron spectrum inferred from a correlated dependent nuclear dataset).

When inputs are correlated, the orthogonality properties of the classical Sobol decomposition [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] are no longer satisfied. As pointed out by several authors [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF][START_REF] Da Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF], a global sensitivity analysis based on this decomposition may lead to erroneous conclusions. Following the work of [START_REF] Stone | The use of polynomial splines and their tensor products in multivariate function estimation[END_REF], later applied to Machine Learning by [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF], and to sensitivity analysis by [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -Application to sensitivity analysis[END_REF], we consider a hierarchically orthogonal decomposition in this paper whose uniqueness has been proven under mild conditions on the dependence structure of the inputs [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -Application to sensitivity analysis[END_REF]. In other words, any model function can be uniquely decomposed as a sum of hierarchically orthogonal component functions. Two summands are considered to be hierarchically orthogonal whenever all of the variables included in one of them are also involved in the other. For a better understanding of the paper, this generalized ANOVA expansion will be referred to as a Hierarchically Orthogonal Functional Decomposition (HOFD).

It is of great importance to develop estimation procedures since the analytical formulation for HOFD is rarely available. In this paper, we focus on an alternative method proposed in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF] to estimate the HOFD components. It consists in constructing a hierarchically orthogonal basis from a suitable Hilbert orthonormal basis. Inspired by the usual Gram-Schmidt algorithm, the procedure recursively builds a multidimensional basis for each component that satisfies the identifiability constraints imposed on this summand. Each component is then well approximated on a truncated basis where the unknown coefficients are deduced by solving an ordinary least-squares regression. Nevertheless, in a high-dimensional paradigm, this procedure suffers from a curse of dimensionality. Moreover, it is numerically observed that only a few of the coefficients are not close to zero, meaning that only a small number of predictors restore the major part of the information contained in the components. Thus, it is important to be able to select the most relevant representative functions and to then identify the HOFD with a limited computational budget.

With this in mind, we suggest in this article to transform the ordinary leastsquares regression into a penalized regression, as has been proposed in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF]. In the present paper, we focus on the L 2 -boosting developed by [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF] to deal with the ℓ 0 penalization. The L 2 -boosting is a greedy strategy that performs variable selection and shrinkage. The choice of such an algorithm is motivated by the fact that the L 2 -boosting is very intuitive and easy to implement. It is also closely related (from the practical point of view)

to the LARS algorithm proposed by [START_REF] Efron | Least angle regression[END_REF], which solves the Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Bühlmann | Statistics for high-dimensional data[END_REF]. The L 2 -boosting and the LARS both select predictors using the maximal correlation with the current residuals.

The question that naturally arises now is the following: provided that the theoretical procedure of component reconstruction is well tailored, do the estimators obtained by the L 2 -boosting converge to the theoretical true sparse parameters when the number of observations tends to infinity? The goal of this paper is to provide an overall consistent estimation of a signal spanned into a large dimensional dictionary derived from a HOFD. Hence, our work significantly improves the results of [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF]: we first address the convergence rate of the empirical HOFD and then use this result to obtain a sparse estimator of the unknown signal. We will need to manage sufficient theoretical conditions to ensure the consistency of our estimator. In addition, we discuss these conditions and provide some numerical examples in which such conditions are fulfilled.

The article is organized as follows. The notation used in the paper is presented in Section 2.1. Section 2.2 provides the HOFD representation of the model function. In Section 2.3, we review the procedure detailed in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF] that consists in constructing well-tailored hierarchically orthogonal bases to represent the components of the HOFD. Finally, we highlight the curse of dimensionality that we are exposed to, and present the L 2 -boosting. Section 3 describes our main theoretical results on the proposed algorithms. One interesting application of the general theory is the global sensitivity analysis (SA).

In Section 4, we apply the L 2 -boosting to estimate the generalized sensitivity indices defined in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -Application to sensitivity analysis[END_REF][START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF]. After recalling the form of these indices, we numerically compare the L 2 -boosting performance with a Lasso strategy and the Forward-Backward algorithm proposed by [START_REF] Zhang | Adaptive forward-backward algorithm for learning sparse representations[END_REF]. Section 5 contains the conclusion, and the proofs of the two main theorems are given in the Supplementary Material.

2. Estimation of the generalized Hoeffding decomposition components

Notation

We consider a measurable function f of a random real vector

X = (X 1 , • • • , X p ) of R p , p ≥ 1.
The response variable Y is a real-valued random variable defined as:

Y = f (X) + ε, (2.1)
where ε stands for a centered random variable independent of X and models the variability of the response around its theoretical unknown value f . We denote the distribution law of X by P X , which is unknown in our setting, and we assume that P X admits a density function p X with respect to the Lebesgue measure on R p . Note that P X is not necessarily a tensor product of univariate distributions since the components of X may be correlated.

Furthermore, we suppose that

f ∈ L 2 R (R p , B(R p ), P X ), where B(R p ) denotes the Borel set of R p . The Hilbert space L 2 R (R p , B(R p ), P X ) is denoted by L 2 R
, for which we use the inner product •, • , and the norm • as follows:

h, g = h(x)g(x)p X dx = E(h(X)g(X)) h 2 = h, h = E h(X) 2 , ∀h, g ∈ L 2 R ,
where E(•) stands for the expected value,

V (•) = E[(• -E(•)) 2 ] denotes the vari- ance, and Cov(•, * ) = E[(• -E(•))( * -E( * ))] the covariance.
For any 1 ≤ i ≤ p, we denote the marginal distribution of X i by P X i and naturally extend the former notation to

L 2 R (R, B(R), P X i ) := L 2 R,i .

The generalized Hoeffding decomposition

Let us denote [1 : k] := {1, 2, • • • , k}, with k ∈ N * , and let S be the collection of all subsets of [1 : p]. We also define S * := S \ {∅}. For u ∈ S, the subvector

X u of X is defined as X u := (X i ) i∈u . Conventionally, for u = ∅, X u = 1.
The marginal distribution (resp. density) of X u is denoted by P Xu (resp. p Xu ).

A functional ANOVA decomposition consists in expanding f as a sum of increasing dimension functions:

f (X) = f ∅ + p i=1 f i (X i ) + 1≤i<j≤p f ij (X i , X j ) + • • • + f 1,••• ,p (X) = u∈S f u (X u ), (2.2) where f ∅ is a constant term, f i , i ∈ [1 : p] are the main effects, f ij , f ijk , • • • , i, j, k ∈ [1 : p]
:= L 2 R (R u , B(R u ), P Xu ).
We then define H 0 u , u ∈ S \ ∅ as follows:

H 0 u = h u ∈ H u , h u , h v = 0, ∀ v ⊂ u, ∀ h v ∈ H 0 v ,
where ⊂ denotes the strict inclusion.

Definition 1 (Hierarchical Orthogonal Functional Decomposition -HOFD). Under Assumptions (C.1) and (C.2) in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -Application to sensitivity analysis[END_REF], the decomposition (2.2) is unique as soon as we assume f u ∈ H 0 u for all u ∈ S.

Remark 1. The components of the HOFD (2.2) are assumed to be hierarchically orthogonal, that is, f u , f v = 0 ∀v ⊂ u.

To obtain more information about the HOFD, the reader is referred to [START_REF] Hooker | Generalized functional anova diagnostics for high-dimensional functions of dependent variables[END_REF] and [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -Application to sensitivity analysis[END_REF]. In this paper, we are interested in estimating the summands in (2.2). As underlined in [START_REF] Huang | Projection estimation in multiple regression with application to functional anova models[END_REF], estimating all components of (2.2) suffers from a curse of dimensionality, leading to an intractable problem in practice. To bypass this issue, we will assume (without loss of generality) that f is centered, so that f ∅ = 0. Furthermore, most of the models are only governed by low-order interaction effects, as pointed out in [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF], [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF] and Li et al. (2010). We thus suppose that f is well approximated by:

f (X) ≃ u∈S * |u|≤d f u (X u ), d ≪ p, (2.3)
so that interactions of order ≥ d + 1 can be neglected. The choice of d, which is directly related to the notion of effective dimension in the superposition sense (see Definition 1 in [START_REF] Wang | The effective dimension and quasi-Monte Carlo integration[END_REF]), is addressed in [START_REF] Zuniga | Metamodelling with independent and dependent inputs[END_REF], but is not of great interest in the present article since it is assumed to be fixed by the user. Even by choosing a small d, the number of components in (2.3) can become prohibitive if the number of inputs p is high. We are therefore interested in estimation procedures under sparse assumptions when the number of variables p is large.

In the next section, we recall the procedure to identify components of (2.3).

As a result of this strategy, we highlight the curse of dimensionality when p becomes large, and we propose to use a greedy L 2 -boosting to tackle this issue.

Practical determination of the Sparse HOFD

General description of the procedure

We propose a two-step estimation procedure in this section to identify the components in (2.3): the first one is a simplified version of the Hierarchical Orthogonal Gram-Schmidt (HOGS) procedure developed in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF], and the second consists of a sparse estimation in the dictionary learned by the empirical HOGS.

In the following, we have chosen to use the so-called L 2 -boosting procedure instead of the widely used Lasso estimator. This choice is motivated by two reasons.

• First, from a technical point of view, the empirical HOGS will produce a noisy estimation of the theoretical dictionary in which the true signal f is expanded. Hence, the arguments produced for the Lasso estimation would have to be completely adjusted to this situation with errors in the variables.

Moreover, as an M-estimator, such a modification is far from being trivial (see [START_REF] Cavalier | Adaptive estimation for inverse problems with noisy operators[END_REF] for an example of oracle inequalities derived from M estimators with noise in the variables). In contrast, the approximation obtained in the empirical HOGS can be easily handled with the boosting algorithm since we just have to quantify how the empirical inner products built with noisy variables are close to theoretical ones. Our proofs rely precisely on this strategy: we obtain a uniform bound on our statistical estimation of the HOGS dictionary, and then take advantage of the sequential description of the boosting with empirical inner products.

• Second, in order to obtain consistent estimation with the boosting procedure, we do not need to make any coherence assumption on the dictionary (such as the RIP assumption of [START_REF] Candes | The Dantzig selector: statistical estimation when p is much larger than n[END_REF] or the weakest RE(s, c 0 ) assumption of [START_REF] Bickel | Simultaneous analysis of lasso and Dantzig selector[END_REF]). Such assumptions are generally necessary to assert some consistency results for the Dantzig and Lasso procedures, such as Sparse Oracle Inequalities (SOI). Nevertheless, it would be only reasonable here to impose these latter assumptions on the theoretical version of the HOGS although it seems difficult to deduce coherence results on the empirical HOGS from coherence results on the theoretical version of the HOGS. Our Theorem 2 below will not produce a SOI in expectation and our results will instead be expressed in probability. We will discuss the asymptotics involved in our Theorem 2 after its statement, and underline the differences with the state of the art results on the Lasso estimator.

To carry out this two-step procedure, we assume that we observe two independent and identically distributed samples, (y r ,

x r ) r=1,••• ,n 1 and (y s , x s ) s=1,••• ,n 2 ,
from the distribution of (Y, X) (the initial sample can be split into such two samples). We define the empirical inner product •, • n and the empirical norm • n associated with an n-sample as:

h, g n = 1 n n s=1 h(x s )g(x s ), h n = h, h n . Also, for u = (u 1 , • • • , u t ) ∈ S, we define the multi-index l u = (l u 1 , • • • , l ut ) ∈ N t .
We use the notation Span {B} to define the set of all finite linear combinations of elements of B, also referred to as the linear span of B.

Step 1 and Step 2 of our sparse HOFD procedure will be described in detail below.

Remark 2. The procedure could be extended to any higher order approximation, but we think that the description of the methodology for d = 2 provides a better understanding. We have thus chosen to only describe this situation for the sake of clarity.

Step 1: Hierarchically Orthogonal Gram-Schmidt Procedure

For each i ∈ [1 : p], let {1 , ψ i l i , l i ∈ N * } denote an orthonormal basis of H i := L 2 (R, B(R), P X i ). For L ∈ N * , for i = j ∈ [1 : p],
we set:

H L ∅ = Span {1} and H L i = Span 1, ψ i 1 , • • • , ψ i L ,
as well as:

H L ij = Span 1, ψ i 1 , • • • , ψ i L , ψ j 1 , • • • , ψ j L , ψ i 1 ⊗ ψ j 1 , • • • , ψ i L ⊗ ψ j L ,
where ⊗ denotes the tensor product between two elements of the basis. We define H L,0 u , the approximation of H 0 u , as:

H L,0 u = h u ∈ H L u , h u , h v = 0, ∀ v ⊂ u, ∀ h v ∈ H L,0 v ,
The recursive procedure below aims at constructing a basis for H L,0 i and a basis for

H L,0 ij for any i = j ∈ [1 : p]. Initialization For any 1 ≤ i ≤ p, define φ i l i := ψ i l i , l i ∈ [1 : L]. Then, as a result of the orthogonality of {ψ i l i , l i ∈ N}, we obtain H L,0 i := Span φ i 1 , • • • φ i L .
For this step, we just need the orthogonality of the constant function equal to 1 with each of the ψ i l i , l i ∈ N * . However, orthogonality is needed for the proof of the consistency of the L 2 -boosting procedure (see the proof of Lemma 4 in the Supplementary Material).

Second order interactions

Let u = {i, j}, with i = j ∈ [1 : p]. Since the dimension of H L ij is equal to L 2 + 2L + 1
, and the approximation space H L,0 ij is subject to 2L + 1 constraints, its dimension is then equal to L 2 . We want to construct a basis for H L,0 ij , which satisfies the hierarchical orthogonal constraints. We are looking for such a basis in the form:

φ ij l ij (X i , X j ) = φ i l i (X i ) × φ j l j (X j ) + L k=1 λ i k,l ij φ i k (X i ) + L k=1 λ j k,l ij φ j k (X j ) + C l ij , (2.4) with l ij = (l i , l j ) ∈ [1 : L] 2 . The constants (C l ij , (λ i k,l ij ) L k=1 , (λ j k,l ij ) L k=1
) are determined by resolving the following constraints:

φ ij l ij , φ i k = 0, ∀ k ∈ [1 : L] φ ij l ij , φ j k = 0, ∀ k ∈ [1 : L] φ ij l ij , 1 = 0.
(2.5)

We first solve the linear system:

A ij λ l ij = D l ij , (2.6) with λ l ij = t λ i 1,l ij • • • λ i L,l ij λ j 1,l ij • • • λ j L,l ij , and
D l ij = -              φ i l i × φ j l j , φ i 1 . . . φ i l i × φ j l j , φ i L φ i l i × φ j l j , φ j 1 . . . φ i l i × φ j l j , φ j L              , A ij = B ii B ij t B ij B jj , with B ij =     φ i 1 , φ j 1 • • • φ i 1 , φ j L . . . . . . φ i L , φ j 1 • • • φ i L , φ j L     .
As shown in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF], A l ij is a definite positive Gramian matrix and (2.6) provides a unique solution in λ l ij . C l ij is then deduced with:

C l ij = -E φ i l i ⊗ φ j l j (X i , X j ) + L k=1 λ i k,l ij φ i k (X i ) + L k=1 λ j k,l ij φ j k (X j ) .
(2.7)

Higher interactions This construction can be extended to any |u| ≥ 3. We refer the interested reader to [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF]. Just note that the dimension of the approximation space H L,0 u is given by L u = L |u| , where |u| denotes the cardinality of u.

Empirical procedure Algorithm 1 below proposes an empirical version of the HOGS procedure. It consists in substituting the inner product •, • with its

empirical version •, • n 1 obtained with the first dataset (y r , x r ) r=1,••• ,n 1 .
Algorithm 1: Empirical HOFD (EHOFD)

Input: Orthonormal system (ψ i l i ) L l i =0 of H i , i ∈ [1 : p], i.i.d. observations O 1 := (y r , x r ) r=1,••• ,n 1 of (2.1), threshold |u max | Initialization: for any i ∈ [1 : p] and l i ∈ [1 : L], first define φi l i ,n 1 = ψ i l i .
• For any u such that 2 ≤ |u| ≤ |u max |, write the matrix Âu n 1 as well as Dlu n 1 obtained using the former expressions with •, • n 1 .

• Solve (2.6) with the empirical inner product

•, • n 1 . Compute λl ij n 1
and Ĉn 1 l ij with (2.7).

• The empirical version of the basis given by (2.4) is then:

∀u ∈ [2 : |u max |] ĤL,0,n 1 u = Span φu 1,n 1 , • • • , φu L |u| ,n 1 .
Step 2: Greedy selection of Sparse HOFD Each component f u of the HOFD defined in Definition 1 is a projection onto H 0 u . Since for u ∈ S * , the space ĤL,0,n 1 u well approximates H 0 u , it is then natural to approximate f by:

f (x) ≃ f (x) = u∈S * |u|≤d fu (x u ), with fu (x u ) = lu β u lu φu lu,n 1 (x u ),
where l u is the multi-index l u = (l i ) i∈u ∈ [1 : L] |u| . For the sake of clarity (since there is no ambiguity), we will omit the summation support of l u in the sequel.

We now consider the second sample (y s , x s ) s=1,••• ,n 2 and we attempt to recover the unknown coefficients β u lu lu,|u|≤d on the regression problem:

y s = f (x s ) + ε s , s = 1, • • • , n 2 .
However, the number of coefficients is equal to d k=1 p k L k . When p becomes large, the usual least-squares estimator is not adapted to estimate the coefficients (β u lu ) lu,u . We then use the penalized regression:

( βu lu ) ∈ Argmin β u lu ∈R 1 n 2 n 2 s=1 y s - u∈S * |u|≤d lu β u lu φu lu,n 1 (x s u ) 2 + λJ(β 1 1 , • • • , β u lu , • • • ), (2.8)
where J(•) is the ℓ 0 -penalty, i.e.,

J(β 1 1 , • • • , β u lu , • • • ) = u∈S * |u|≤d lu 1(β u lu = 0).
Of course, such an optimization procedure is not tractable and we instead consider the relaxed L 2 -boosting [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF] to solve this penalized problem. Mimicking the notation of [START_REF] Temlyakov | Weak Greedy Algorithms[END_REF] and [START_REF] Champion | Sparse regression and support recovery with L 2 -boosting algorithm[END_REF],

we define the dictionary D of functions as:

D = { φ1 1,n 1 , • • • φ1 L,n 1 , • • • , φu 1,n 1 , • • • , φu Lu,n 1 , • • • }.
The quantity G k ( f ) denotes the approximation of f at step k as a linear combination of elements of D. At the end of the algorithm, the estimation of f is denoted by f . The L 2 -boosting is described in Algorithm 2. Initialization: G 0 ( f ) = 0.

for k = 1 to k up do 1. Select φu k lu k ,n 1 ∈ D such that: Y -G k-1 ( f ), φu k lu k ,n 1 n 2 = max φu lu,n 1 ∈D Y -G k-1 ( f ), φu lu,n 1 n 2 . (2.9) 2.
Compute the new approximation of f as:

G k ( f ) = G k-1 ( f ) + γ Y -G k-1 ( f ), φu k lu k ,n 1 n 2 • φu k lu k ,n 1 . (2.10) end Output: f = G kup ( f ).
For any step k, Algorithm 2 selects a function from D that provides sufficient information about the residual Y -G k-1 ( f ). The shrinkage parameter γ is the standard step-length parameter of the boosting algorithm. It actually smoothly inserts the next predictor into the model, making a refinement of the greedy algorithm possible and statistically guaranteeing its convergence rate.

Remark 3. In a deterministic setting, the shrinkage parameter is not really useful and may be set to 1 (see [START_REF] Temlyakov | Weak Greedy Algorithms[END_REF] for further details). It is particularly useful from a practical point of view to smooth the boosting iterations.

An algorithm for our new sparse HOFD procedure Algorithm 3 below now provides a simplified description of our sparse HOFD procedure, whose steps have been described above.

Algorithm 3: Greedy Hierarchically Orthogonal Functional Decomposition

Input: Orthonormal system (ψ i l i ) L l i =0 of L 2 (R, B(R), P X i ), i ∈ [1 : p], i.i.d. observations O := (y j , x j ) j=1...n of (2.1) Initialization: Split O in a partition O 1 ∪ O 2 of size (n 1 , n 2 ). • For any u ∈ S, use Step 1 with observations O 1 to construct the approximation ĤL,0,n 1 u := Span φu 1,n 1 , • • • , φu Lu,n 1 of H L,0 u (see Algorithm 1).
• Use an L 2 -boosting algorithm on O 2 with the random dictionary

D = { φ1 1,n 1 , • • • φ1 L,n 1 , • • • , φu 1,n 1 , • • • , φu Lu,n 1 , • • • } to obtain the Sparse Hierarchically Orthogonal Decomposition (see Algorithm 2).
We now obtain a strategy to estimate the components of the decomposition (2.3) in a high-dimensional paradigm. We aim to show that the obtained estimators are consistent, and that the two-step procedure (summarized in Algorithm 3) is numerically convincing. The next section is devoted to the asymptotic properties of the estimators.

Consistency of the estimator

In this section, we study the asymptotic properties of the estimator f obtained from Algorithm 3 described in Section 2. To do this, we restrict our study to the case of d = 2 and assume that f is well approximated by first and second order interaction components (see Remark 4 below). Hence, the observed signal Y may be represented as:

Y = u∈S * |u|≤2 lu β u,0 lu φ u lu (X u ) + ε, E(ε) = 0, E(ε 2 ) = σ 2 ,
where

β 0 = β u,0 lu lu,u
is the true parameter, and the functions φ u lu lu , |u| ≤ 2 are constructed according to the HOFD described in Section 2.3. We assume that we have an n-sample of observations divided into two samples, O 1 , and

O 2 . Samples in O 1 (resp. in O 2 ) of size n 1 = n/2 (resp. of size n 2 = n/2) are
used for the construction of φu lu,n 1 lu,u described in Algorithm 1 (resp. for the L 2 -boosting Algorithm 2 to estimate β u lu lu,u . The goal of this section is to study the consistency of f = G kn ( f ) when the sample size n tends to infinity. Its objective is also to determine an optimal number of steps k n necessary to obtain a consistent estimator from Algorithm 2.

Remark 4. We choose the truncature order d = 2 in order to simplify the presentation, but it may be extended to arbitrary larger thresholds independent of the sample size n. This choice is legitimate as soon as the function f is well approximated by low interaction components and this assumption is well suited for many practical situations [START_REF] Rabitz | Efficient input-output model representations[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF]. Indeed, a datadependent choice of d n (with d n -→ +∞ as n → +∞) would rely on a smoothness assumption on the signal f with respect to the order of the considered interactions by using the size of the bias term induced by the truncature given in Theorem 5 of [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF]. However, this challenging task is far beyond the scope of this paper and we have chosen to leave this problem open.

Assumptions

We first briefly recall some notation: for all sequences (a n ) n≥0 , (b n ) n≥0 , we

write a n = O n→+∞ (b n ) when a n /b n is a bounded sequence for large enough n. Now, for any random sequence (X n ) n≥0 , X n = O P (a n ) means that |X n /a n | is bounded in probability.
We have chosen to present our assumptions in three parts to deal with the dimension, the noise and the sparseness of the entries.

Bounded Assumptions (H b ) The first set of hypotheses matches the bounded case and is adapted to the special situation of bounded support for the random variable X, for example, when each X j follows a uniform law on a compact set

K j ⊂ K where K is a compact set of R independent of j ∈ [1 : p].
It is referred to as (H b ) in the sequel and corresponds to the following three conditions:

(H 1 b ) M := sup i∈[1:p] l i ∈[1:L] φ i l i (X i ) ∞ < +∞, (H 2 b )
The number of variables p n satisfies:

p n = O n→+∞ (exp(Cn 1-ξ ))
, where 0 < ξ ≤ 1 and C > 0.

(H 3,ϑ b ) The Gram matrices A ij introduced in (2.6) satisfies:

∃C > 0 ∀(i, j) ∈ [1 : p n ] 2 det(A ij ) ≥ Cn -ϑ ,
where det referes to the determinant of a matrix.

Roughly speaking, this will be the favorable situation from a technical point of view since it will be possible to apply a matricial Hoeffding type inequality. It may be possible to slightly relax such a hypothesis using a sub-exponential tail argument. For the sake of simplicity, we have chosen to only restrict our work to the settings of (H b ).

Regardless of the joint law of the random variables (X 1 , . . . , X p ), it is always possible to build an orthonormal basis (φ i l i ) 1≤l i ≤L from a bounded (frequency truncated) Fourier basis and, therefore, (H 1 b ) is not as restrictive in practice. Assumption (H 2 b ) deals with the high dimensional situation. We are in fact interested in practical situations where the number of variables can be much larger than the number of observations n. Hence, in our mathematical study, the number of variables p n can grow exponentially fast with the number of observations n. This obviously implies that the collection of subsets u also depends on n and will now be designated as S * n . As a consequence, S * n also increases rapidly and is much larger than n.

Note that Hypothesis (H 3,ϑ b ) stands for a lower bound of the determinant of the Gram matrices involved in the HOFD. It is shown in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF] that each of these Gram matrices is invertible and, as a result, each det(A ij ) is positive. Nevertheless, if ϑ = 0, this hypothesis assumes that such an invertibility is uniform over all choices of tensor (i, j). This hypothesis may be too strong for a large number of variables p n → +∞ when ϑ = 0. However, when ϑ > 0, Hypothesis (H 3,ϑ b ) drastically relaxes the case ϑ = 0 and becomes very weak. The verification of (H 3,ϑ b ) requires the computation of an order of p 2 n determinants of size L 2 × L 2 . We have checked this assumption in our experiments. However, for very large values of n, this may become impossible from a numerical point of view.

In the following, the parameters ϑ and ξ will be related to each other and we will obtain a consistency result of the sparse HOFD up to the condition ϑ < ξ/2. This constraint implicitly limits the size of p n since log

p n = O n→+∞ (n 1-ξ ).
Noise Assumption (H ε,q ) We will assume the noise measurement ε to obtain some bounded moments of sufficiently high order, which is true for Gaussian or bounded noise. This assumption is given by: (H ε,q ) E(|ε| q ) < ∞, for one q ∈ R + .

Sparsity Assumption (H s,α ) The last assumption concerns the sparse representation of the unknown signal described by Y in the basis (φ u lu (X u )) u . Such a hypothesis will be useful to assess the statistical performance of the L 2 -boosting and will be referred to as (H s,α ) below. It is legitimate due to our high dimension setting and our motivation to identify the main interactions X u .

(H s,α ) There exists α > 0 such that the parameter β 0 satisfies:

β 0 ℓ 1 := u∈S * n |u|≤d lu β u,0 lu = O n→+∞ (n α ).

Main results

We recall below that . is the L 2 norm on functions decomposed in the orthonormal basis (φ u lu ) u . We first provide our main result on the efficiency of the EHOFD (Algorithm 1).

Theorem 1. Assume that (H b ) holds with ξ (resp. ϑ) given by (H 2 b ) (resp. (H 3,ϑ b )) and that a constant Λ exists such that λ l ij 2 ≤ Λ for any couple (i, j).

Then, if ϑ < ξ/2, the sequence of estimators φu lu,n 1 u satisfies:

sup u∈S * n ,|u|≤d lu φu lu,n 1 -φ u lu = ζ n,0 = O P (n ϑ-ξ/2 ).
The proof of this theorem can be found in the Supplementary Material.

Let us mention the contribution of Theorem 1 compared to the results obtained in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF]. Proposition 5.1 of [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF] leads to an almost sure convergence of their estimator without any quantitative rate when the number of functions in the HOFD is kept fixed and does not grow with the number of observations n. In contrast, in our high dimensional paradigm, we allow S * n to grow with n and also obtain an almost sure result associated with a convergence rate. This will be essential for the derivation of our next result.

Our second main result concerns the L 2 -boosting that recovers the unknown f up to a preprocessing estimation of ( φu lu,n 1 ) lu,u on a first sample O 1 . Such a result is satisfied provided the sparsity assumption (H s,α ) holds. We assume that

Y = f (X) + ε, f (X) = u∈S * n |u|≤d lu β u,0 lu φ u lu (X u ) ∈ H L u ,
where β 0 = (β u,0 lu ) lu,u is the true parameter that expands f . To the best of our knowledge, such a high dimensional inference with noise in the variables appears to be novel. As already pointed out above, the greedy boosting seems to be a well-tailored approach to handle noisy dictionaries in comparison to a penalized regression strategy, which relies on a somewhat unverifiable "RIP-type" hypothesis on the learned dictionary.

Theorem 2 (Consistency of the L 2 -boosting). Consider an estimation f of f from an i.i.d. n-sample broken down into O 1 ∪ O 2 . Assume that functions φu lu,n 1 lu,u are estimated from the first sample O 1 under (H b ) with ϑ < ξ/2, and that a constant Λ exists such that λ l ij 2 ≤ Λ for any couple (i, j). f is then defined by (2.10) of Algorithm 2 on O 2 as:

f (X) = G kn ( f ), with f = u∈S * n |u|≤d lu β u,0 lu φu lu,n 1 (X u ).
If we assume that (H s,α ) and (H ε,q ) are satisfied with q > 4/ξ and α < ξ/4-ϑ/2, then a sequence k n := C log n exists, where C < ξ/2-ϑ-2α 2•log 3 such that:

f -f P -→ 0, when n → +∞.
In particular, for Gaussian noises that possess arbitrary large moments, the constraint on q disappears and Theorem 2 can be applied as soon as ξ < 1.

Let us discuss the asymptotic setting involved in our theorem. First, our result is a result in probability rather than in expectation. It is a frequently encountered fact that SOI in expectation are derived with additional assumptions on the coherence of the dictionary; some detailed discussions can be found in [START_REF] Bickel | Simultaneous analysis of lasso and Dantzig selector[END_REF] and [START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF]. With some coherence and boundedness assumptions, [START_REF] Bickel | Simultaneous analysis of lasso and Dantzig selector[END_REF] deduced convergence rates of the Lasso estimator in expectation as soon as:

β 0 ℓ 0 log(p) n -→ 0. (3.1)
Later, [START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF] extended the study of the Lasso behavior with a result on the Lasso estimator on bounded variables without any coherence assumption and showed a consistency result in probability when:

β 0 ℓ 1 log(p) n -→ 0. (3.2)
Hence, the rate is damaged by the appearance of the √ . in (3.2) in comparison with (3.1). Concerning the boosting algorithm, [START_REF] Champion | Sparse regression and support recovery with L 2 -boosting algorithm[END_REF] also obtained consistency results in probability under the asymptotic setting given in (3.2) without a coherence assumption. It should be observed that our results with a noisy dictionary requires that:

inf i,j det(A ij ) -1 β 0 2 ℓ 1 log p n -→ 0 as n -→ +∞, (3.3)
which is a stronger assumption in comparison with (3.2). From a technical point of view, the asymptotic setting is due to inequality (S4.10) where β 0 2 ℓ 1 ζ n appears instead of β 0 ℓ 1 ζ n for boosting algorithms without noise on the variables (see the proof of Theorem 2 in the Supplementary Material).

In favorable cases where all linear systems defined through the Gram matrices A ij are well conditioned, ϑ = 0 and the condition becomes β 0 2 ℓ 1 log p n -→ 0, and there is still a price to pay for the preliminary estimation of the elements of the HOGS. Theorem 2 can be applied only for sequences of coefficients such that β u,0 lu L 1 n 1/4 . Note also that the degeneracy of the Gram determinants must be strictly larger than n -1/2 . For example, when ϑ = 1/4, the norm β u,0 lu L 1 cannot be larger than n 1/8 . We briefly describe the proof below and provide the technical details in the Supplementary Material.

Sketch of Proof of Theorem 2. Mimicking the scheme of [START_REF] Bühlmann | Boosting for high-dimensional linear models[END_REF] and [START_REF] Champion | Sparse regression and support recovery with L 2 -boosting algorithm[END_REF], the proof first consists in defining the theoretical residual of Algorithm 2 at step k as:

R k ( f ) = f -G k ( f ) = f -G k-1 ( f ) -γ Y -G k-1 ( f ), φu k lu k ,n 1 n 2 • φu k lu k ,n 1 (3.4)
Furthermore, following the work of [START_REF] Champion | Sparse regression and support recovery with L 2 -boosting algorithm[END_REF], we introduce a phantom residual in order to reproduce the behavior of a deterministic boosting, studied in [START_REF] Temlyakov | Weak Greedy Algorithms[END_REF]. This phantom algorithm is the theoretical L 2boosting, performed using the randomly chosen elements of the dictionary by Equations (2.9) and (2.10), but updated using the deterministic inner product.

The phantom residuals Rk ( f ), k ≥ 0, are defined as follows:

R0 ( f ) = f Rk ( f ) = Rk-1 ( f ) -γ Rk-1 ( f ), φu k lu k ,n 1 φu k lu k ,n 1 , (3.5)
where φu k lu k ,n 1 has been selected with Equation (2.9) of Algorithm 2. The aim is to decompose the quantity ff to introduce the theoretical residuals and the phantom ones:

f -f = G kn ( f ) -f ≤ f -f + R kn ( f ) -Rkn ( f ) + Rkn ( f ) . (3.6)
We then have to show that each term on the right-hand side of (3.6) converges towards zero in probability.

Numerical Applications

This section is devoted to the numerical efficiency of the two-step procedure given in Section 2, and primarily focuses on the practical use of the HOFD through sensitivity analysis (SA). SA aims to identify the most contributive variables to the variability of a regression model [START_REF] Saltelli | Sensitivity Analysis[END_REF][START_REF] Cacuci | Sensitivity and Uncertainty Analysis, Volume II: Applications to Large-Scale Systems[END_REF]. The most common quantification is a variance-based index, known as the Sobol index [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. This measure relies on the Hoeffding decomposition that provides an elegant and meaningful theoretical framework when inputs are known to be independent. However, as mentioned in the introduction, the interpretation of such indices may be irrelevant when strong dependencies arise.

The HOFD presented in Section 2.2 is of great interest in this situation because it provides a general and rigorous multivariate regression extension that can be used to define sensitivity indices well-tailored to dependent inputs. As detailed in [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables -Application to sensitivity analysis[END_REF], the model variance can be expanded as follows:

V (Y ) = u∈S * n   V (f u (X u )) + u∩v =u,v Cov(f u (X u ), f v (X v ))  
Therefore, to measure the contribution of X u , for |u| ≥ 1, in terms of model variability, it is then quite natural to define a sensitivity index S u as follows:

S u = V (f u (X u )) + u∩v =u,v Cov(f u (X u ), f v (X v )) V (Y ) . (4.1)
Furthermore, we deduce the empirical estimation of (4.1) once we have applied the procedure described in Algorithm 3 to obtain ( fu , fv , u ∩ v = u, v).

Description

We end this paper with a short simulation study, focused primarily on the performance of the greedy selection algorithm for the prediction of generalized sensitivity indices. Since the estimation of these indices consists in estimating the summands of the generalized functional ANOVA decomposition (referred to as HOFD), we begin by constructing a hierarchically orthogonal system of functions to approximate the components. As pointed out above (see Assumption (H 3,ϑ b ) in Theorem 1 and 2), the invertibility of each linear system plays an important role in our theoretical study. For each situation, we therefore measured the degeneracy of the matrices involved, given by:

d(A) = inf i,j∈[1:p] det(A ij ).
We then use a variable selection method to select a sparse number of predictors. The goal is to numerically compare three variable selection methods: the L 2 -boosting, the Forward-Backward greedy algorithm (referred to as FoBa below), and the Lasso estimator. As pointed out above, we have an n-sample of i.i.d.

observations (y s , x s ) s=1,••• ,n broken down into two samples of size n 1 = n 2 = n/2.
The first sample is used to construct the system of functions according to Algorithm 1. The second sample is used to solve the penalized regression problem given by (2.8) and illustrated here:

( βu lu ) lu,u ∈ Argmin β u lu ∈R 1 n 2 n 2 s=1 y s - u∈S |u|≤d lu β u lu φu lu,n 1 (x s u ) 2 + λJ(β 1 1 , • • • , β u lu , • • • ).
We will now briefly describe how we use the Lasso, the FoBa and the Boosting.

Feature selection Algorithms

FoBa procedure The FoBa algorithm, as well as the L 2 -boosting, use a greedy exploration to minimize the previous criterion when J(•) is a ℓ 0 penalty, i.e.,

J(β 1 1 , • • • , β u lu , • • • ) = u∈S * n |u|≤d lu 1(β u lu = 0).
This algorithm is an iterative scheme that sequentially selects or deletes an element of D that has the least impact on the fit, i.e., that significantly reduces the model residual. This algorithm is described in [START_REF] Zhang | Adaptive forward-backward algorithm for learning sparse representations[END_REF] and used for HOFD in [START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: Numerical methods[END_REF]. We refer to these references for a more indepth description of this algorithm. This procedure depends on two shrinkage parameters, ǫ and δ. The parameter ǫ is the stopping criterion that predetermines if a large number of predictors is going to be introduced into the model.

The second parameter, δ ∈]0, 1], offers a flexibility in the backward step since it allows the algorithm to smoothly eliminate a predictor at each step.

In our numerical experiments, we found a well-suited behavior of the FoBa procedure with ǫ = 10 -2 and δ = 1/2.

Calibration of the Boosting As previously reported by [START_REF] Champion | Sparse regression and support recovery with L 2 -boosting algorithm[END_REF], we fixed the shrinkage parameter to γ = 0.7 since it provides a suitable value for high dimensional regression, even though we did not find any extreme differences when γ varies in [0.5; 1[. Since the optimal value for k up is unknown in practice, we use a C p Mallows-type criterion to fix the optimal number of iterations. This stopping criterion is much more important than the choice of the shrinkage parameter. It is, of course, induced by γ since it depends on the sequence of the boosting iterations.

Like in the LARS algorithm, we follow the recommendations of [START_REF] Efron | Least angle regression[END_REF] to select the best solution. First, we define a large number of iterations, say K. For each step k ∈ {1, • • • , K}, the boosting algorithm computes an estimation of the solution β(k). On the basis of this, we compute the following quantity:

E Boost k = 1 n n 2 s=1 y s - φu lu,n 1 ∈D βu lu (k) φu lu,n 1 (x s u ) 2 -n 2 + 2k,
where the implied set of functions φu lu,n 1 has been selected through the first k steps of the algorithm. Finally, we choose the optimal number of selected functions kup such that: kup = Argmin

k=1,••• ,K E Boost k .
Lasso algorithm Since the ℓ 0 strategy is very difficult to handle and may suffer from a lack of robustness, the ℓ 0 penalty is often replaced by the λ × ℓ 1 strategy that yields the Lasso estimator for a given penalization parameter λ > 0, i.e.,

J(β 1 1 , • • • , β u lu , • • • ) = u∈S * n |u|≤d lu β u lu = 0 .
Several algorithms have been proposed in the literature to solve the Lasso regression. One of the most popular is the LARS method, described in [START_REF] Efron | Least angle regression[END_REF], because it performs a solution that coincides with the theoretical regularization path { β(λ), λ ∈ R + }. However, the LARS strategy is very expensive in large Lasso problems. To make a good numerical comparison with the greedy algorithms, we choose to perform a coordinate descent algorithm proposed by [START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF][START_REF] Friedman | Pathwise coordinate optimization[END_REF] because of its low computational cost compared to the LARS implementation. The tuning parameter λ is first selected by generalized cross-validation, and the Lasso Coordinate Descent (LCD) algorithm is performed with the R lassoshooting package.

Datasets

Each experiment on each dataset was randomly reproduced 50 times to compute the Monte-Carlo errors. Since each dataset has very few instances, the size L of the initial orthonormal systems has to be small. Here, we arbitrarily choose 5 ≤ L ≤ 8 and the approximation performance does not suffer from the sensitivity of L in these models.

First Dataset: the Ishigami function Well known in sensitivity analysis, the analytical form of the Ishigami model is given by: Y = sin(X 1 ) + a sin 2 (X 2 ) + bX 4 3 sin(X 1 ),

where we set a = 7 and b = 0.1, and where it is assumed that the inputs are independent. In the numerical experiment, we consider the following cases:

1. For all i = 1, 2, 3, the inputs are uniformly distributed on [-π, π]. We choose n = 300 observations, with the first eight Legendre basis functions (L = 8).

2. For all i = 1, 2, 3, the inputs are uniformly distributed on [-π, π]. We choose n = 300 observations, with the first eight Fourier basis functions.

Each time, the number of predictors is m n = pL + p 2 L 2 = 408 ≥ n.

Second Dataset: the g-Sobol function This function is referred to in [START_REF] Saltelli | Sensitivity Analysis[END_REF], and 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 ⊆ [1 : p].
Here, we take p = 25 and a = (0, 0, 0, 1, 1, 2, 3, 4.5, 4.5, 4.5, 9, 9, 9, 9, 9, 99, • • • , 99).

For the construction of the hierarchical basis functions, we choose the first five Legendre polynomials (L = 5). We use n = 2000 evaluations of the model and the number of predictors m n = pL + p 2 L 2 = 7625, which clearly exceeds the sample size n.

Third dataset: dependent inputs The third data set stands for a rarely investigated situation where the inputs are correlated. As proposed by [START_REF] Mara | Variance-based sensitivity analysis of computer models with dependent inputs[END_REF], we generate a sample set according to the following distribution: X 1 and X 2 are uniformly sampled in the set S:

S := (x 1 , x 2 ) ∈ [-1, 1] 2 | 2x 2 1 -1 ≤ x 2 ≤ 2x 2 1 .
Furthermore, X 3 is also sampled uniformly in [-1; 1]. Then, Y is built following

Y = X 1 + X 2 + X 3 .
The inputs X 1 and X 2 are clearly not independent and we do not exactly know the analytical Sobol indices. We choose n = 100 observations, with the first six Legendre basis functions (L = 6).

The tank pressure model

This real case study concerns a shell closed by a cap and subject to an internal pressure. E cap ) and the yield strength (σ y,shell and σ y,cap ) of the shell and the cap. The last parameter is the internal pressure (P int ) applied to the shell. Some strong correlations exist between some of the inputs in the system as a result of the constraints of manufacturing processes, for example, between the shell radius and its thickness. The system is modeled by a 2D finite element ASTER code.

Input distributions are provided in Table 4.1. The large variability of P int in the model corresponds to the different internal pressure values that could be applied to the shell by the user.

To measure the contribution of the correlated inputs to the output variability, we estimate the generalized sensitivity indices. We do n = 1000 simulations. We use the first Hermite basis functions, whose maximum degree is 5 for every parameter.

Inputs Distribution

R int U([1800; 2200]), γ(R int , T shell ) = 0.85 T shell U([360; 440]), γ(T shell , T cap ) = 0.3 T cap U([180; 220]), γ(T cap , R int ) = 0.3 E cap αN (µ, Σ) + (1 -α)N (µ, Ω) σ y,cap α = 0.02, µ =    210 500   , Σ =    350 0 0 29   , Ω =    175

Results

We consider both the estimation of the sensitivity indices, the ability to select the good representation of the different signals, and the computation time needed to obtain the sparse representation. "Greedy" refers to the Foba procedure and "LCD" refers to the Lasso coordinate descent method. Our method is, of course, referred to as "Boosting". This may be due to the influence of the dependency among the inputs X 1 and X 2 in this dataset on the Lasso estimator.

Furthermore, it clearly appears in Table 4.2 that our proposed L 2 -boosting is the fastest method. This is particularly true on the 25-dimension g-Sobol function where the fraction of additional time required by the LCD algorithm in comparison to the L 2 -boosting is about 100. Although we do not have access to the theoretical support recovery β 0 , we can observe that the results of the L 2 -boosting are equivalent to those of other algorithms in terms of its feature selection ability. Hence, for the same degree of accuracy, our method seems to be much faster. Note that we have computed the maximal "degeneracy" that is involved in the resolution of the linear systems and quantified by Assumption (H 3,ϑ b ) in column 2 of Table 4.3. In many cases, we obtain a significantly larger value than 0. The third column of Table 4.3 shows the admissible size of the parameter ϑ, and we can check that the number of variables p n allowed by (H 2 b ) and the balance between ξ and ϑ (ξ should be greater than 2ϑ in our theoretical results) is not restrictive since n 1-2ϑ is always significantly greater than log(m n ) in Table 4

Conclusions and Perspectives

This paper provides a rigorous framework for the hierarchically orthogonal Gram-Schmidt procedure in a high-dimensional paradigm, with the use of the greedy L 2 -boosting. Overall, the procedure falls into the category of sparse estimation with a noisy dictionary, and we demonstrate its consistency up to some mild assumptions about the structure of the real underlying basis. From a mathematical point of view, assumption (H 1 b ) presents a restrictive condition, and to relax it would open a wider class of basis functions for applications. We leave this development open for a future study, which could be based either on the development of a concentration inequality for unbounded random matrices or on a truncating argument. It also appears that our algorithm produces very satisfactory numerical results through our three datasets as a result of its very low computational cost. It can also be extended with some further numerical work to a larger truncation order of d ≥ 3. Such an improvement may also be of interest from a theoretical point of view when dealing with a function that smoothly depends on the interaction order. In particular, a data-driven adaptive choice of d may be of practical interest in the future.

Algorithm 2 :

 2 The L 2 -boostingInput: Observations O 2 := (y s , x s ) s=1,••• ,n 2 , shrinkage parameter γ ∈]0, 1]and number of iterations k up ∈ N * .

  Figure 4.1 illustrates a simulation of tank distortion. We are interested in the von Mises stress, detailed in von Mises (1913) on the point y indicated in Figure 4.1. The von Mises stress makes it possible to predict materialyielding that occurs when the material yield strength is reached. The selected point y corresponds to the point for which the von Mises stress is maximal in the tank. We therefore want to prevent the tank from material damage induced by plastic deformations. In order to provide a large panel of tanks able to resist the internal pressure, a manufacturer wants to know the parameters that contribute the most to the von Mises criterion variability. In the model that we propose, the von Mises criterion depends on three geometrical parameters: the shell internal radius (R int ), the shell thickness (T shell ), and the cap thickness (T cap ). It also depends on five physical parameters concerning Young's modulus (E shell and

Figure 4 . 1 :

 41 Figure 4.1: Tank distortion at point y

Figure 4 . 3 :

 43 Figure 4.3: Representation of the first-order components on the first dataset (Ishigami function) described through the Fourier basis.

Figure 4 . 4 :

 44 Figure 4.4: Representation of the first-order components on the second dataset (g-Sobol function).

Figure 4 . 5 :

 45 Figure 4.5: Representation of the first-order components on the third dataset (dependent inputs).

  Figure 4.6: Dispersion of the first order sensitivity indices of the tank model parameters for the three methods.

  are the interaction effects, and the last component f 1,••• ,p is the

	residual.
	Decomposition (2.2) is generally not unique. However, under mild assump-
	tions on the joint density p X (see Assumptions (C.1) and (C.2) in Chastaing
	et al. (2012)), the decomposition is unique under some additional orthogonality
	assumptions.

More precisely, let us introduce H ∅ = H 0 ∅ as the set of constant functions, and for all u ∈ S * , H u

Table 4

 4

.1: Description of inputs of the shell model

Table 4 .

 4 .3. 3: Degeneracy of the linear systems and admissible size of m n (n 1-2ϑ should be greater than log(m n ).

	Dataset Ishigami function Case1	log(n) Degeneracy d(A) ϑ ≥ log(1/d(A)) 0.6388 [0.0786, +∞[	n 1-2ϑ 122.3821	log(m n ) 6.0113
	Ishigami function Case1	0.76	[0.0481, +∞[	173.3094	6.0113
	g-Sobol function	0.9745	[0.0034, +∞[	1899	8.9392
	Dependent inputs	0.628	[0.101, +∞[	39.4457	4.8363

Acknowledgment The authors are indebted to Fabrice Gamboa for his stimulating discussions and his numerous suggestions on the subject. The authors also thank the two anonymous referees and the associate editor for their helpful comments and suggestions on an earlier version of the manuscript.

are obviously unknown, but we obtain similar results for the three methods.

Finally, as illustrated in Figure 4.6, the most contributive parameter to the von Mises criterion variability is the internal pressure P int , which is not surprising. Concerning the geometric characteristics, the main parameters of the three methods are cap thickness, T cap , and shell thickness, T shell , using their expensive code, although the shell internal radius does not seem to be that important.

Computation time and accuracy The performances of the three methods are illustrated in Table 4.2 on the basis of their computational cost and the accuracy of the feature selection.

Regarding the statistical accuracy, it should be noted that each estimator of high dimensional regression possesses a comparable dispersion on all the datasets and performs quite similarly on the first dataset. The Lasso estimator seems a little bit imprecise in the third dataset in comparison with the FoBa and Boosting methods. Finally, the LCD method is also outperformed on the third dataset (with dependent inputs): it selects a significantly larger number of sensitivity