Ikram Jomaa

Noureddine Issaoui

Thierry Roisnel

Houda Marouani
email: houdamarouani2015@gmail.com

Insight into non-covalent interactions in a tetrachlorocadmate salt with promising NLO properties: Experimental and computational analysis

Keywords: Cd(II) complex: X-ray diffraction, DFT calculations, Hirshfeld surface, AIM, NLO

This paper deals with the crystal structure of the new non-centrosymmetric organic-inorganic hybrid material, (C 6 H 14 N) 2 [CdCl 4 ]. Single crystal X-ray diffraction analysis shows that this compound crystallizes in the orthorhombic system, with the space group Cmc2 1 and the following parameters a = 27.257 (2) Å, b = 8.3560 (6) Å, c = 7.8872 (5) Å, V = 1796.4 (2) Å 3 and Z = 4. The structure provides a new interesting example of infinite inorganic chains of [CdCl 5 ] n following the c crystallographic direction. The monoprotonated cyclohexylammonium cations are linked to the anions via electrostatic, multiple bifurcated N-H…Cl hydrogen bonds and van der Waals interactions. To support experimental results, DFT calculations have been accomplished via the B3LYP method with 6-311++G(d,p) and LANL2DZ basis set on molecular geometry, vibrational and electronic properties. The noncovalent interactions were studied through AIM and RDG analysis and quantitatively using the Hirshfeld surfaces (HS) associated with 2D fingerprint plots. The NLO properties have been also investigated by DFT and compared to the urea reference

ACCEPTED MANUSCRIPT

Introduction

Organo-metal halides assembly investigation represents for a quite while one of the most active areas of chemical research and material science due to their designable and tuneable characteristics as well as their potential applications in numerous fields including In particular, Halogenocadmates(II) crystal chemistry is exceptionally diverse due to the fact that Cd 2+ is a d 10 metal ion. It demonstrated a propensity for forming a variety of crystalline structures by self-assembling from suitable solution mixtures, their structural flexibility enables a great variety of structural archetypes to be constructed in the form of onedimensional linear chains [START_REF] Chu | An ABX3 organic-inorganic perovskitetype material with the formula (C 5 N 2 H 9 )CdCl 3 : Application for detection of volatile organic solvent molecules[END_REF], two-dimensional layers [START_REF] Jin | Structural diversity of a series of chlorocadmate(II) and chlorocuprate(II) complexes based on benzylamine and its N-methylatedderivatives[END_REF], and ribbon-like structures [START_REF] Hajji | Structural and computational analyses ACCEPTED MANUSCRIPT of a 2-propanolammonium-chlorocadmate(II) assembly: Pivotal role of hydrogen bonding and H-H interactions[END_REF].

Else, cyclohexylammonium organic entities have been found to be a useful reagent for metal ions, due to their reactivity and adaptability in stabilizing various structures [START_REF] Yun | Synthesis and crystal structure of bis(dicyclohexylammonium) Tetracyanonickel(II)[END_REF][START_REF] Sarr | Bis(cyclohexylammonium) tetra-chlorido-diphenyl-stannate(IV)[END_REF].

Cyclohexylammonium organic molecules can be associated with chlorocadmates to produce systems with a potentially forceful network of hydrogen bond interactions.

From a structural point of view, two categories of organo-metal halides hybrids are distinguished [10]: Class I where both organic and inorganic constituents are linked by weak non covalent interactions while the components of class II hybrids are combined through stronger covalent bonds. The weakness of the bonds in Class I is responsible for the flexibility and dynamics of structures. Moreover, the assembly of organo-metal halides can create noncentrosymmetric materials with important nonlinear optical coefficients and specific physical and chemical properties. Developments in the field of non-linear optics play a major role in the field of photonics including optical information processing, sensor protector applications, data storage, etc [START_REF] Suresh | Review on theoretical aspect of nonlinear optics[END_REF].

In the present contribution to the above materials, generally expressed as (R-NH 3 ) 2 MX 4 (where R = is an organic cation, M = divalent metal, and X = halogen), and that consists of an extended network of metal halide tetrahedra or octahedra (typically being observed), sharing faces, edges, or vertex and alternating with a monolayer or multilayer organic group, we report the crystal structure study in detail of (C 6 H 14 N) 2 [CdCl 4 ] compound, Hirshfeld surface analysis, and fingerprints plots calculations. Furthermore, even though computational analysis has become an impressive characterization approach only a modest number of papers report some theoretical ACCEPTED MANUSCRIPT calculations of organic-chlorocadmates compounds [START_REF] Jomaa | C 8 H 14 N 2 )[END_REF][START_REF] Yang | Crystal structure,vibrational spectra, optical properties and density functional theoretical approach of[END_REF]. So here a state-of-the-art computational study was undertaken on the occurrence and nature of weak non covalent interactions in the assembly of the title structure, and have been as well combined with experimental results throughout this work for structural, vibratory, and electronic characterization. Topological properties and RDG analysis were also computed.

Experimental details

Synthesis and crystallization

Metal/amine molar ratio 1:2 solution of: CdCl 2 .H 2 O (0.201 g, 1 mmol, Aldrich, purity 98%) mixed in dilute HCl (5 mL, 3 M) was added drop wise to a solution of cyclohexylamine (0.198 g, 2 mmol, Aldrich, purity 99%) dissolved in ethanol. The mixture is left under stirring at room temperature for 30 minutes. One week later well-defined single crystals of (C 6 H 14 N) 2 [CdCl 4 ] were formed. Schematically the reaction can be written as follow:

2 C 6 H 13 N + CdCl 2 . H 2 O + 2 HCl EtOH / H 2 O C 6 H 14 N 2 CdCl 4 + H 2 O
Elemental analysis, calc. (found): C, 31.68% (31.52); H, 6.16% (6.12); N, 6.16% (6.10); Cd, 24.73% (24.79).

X-ray data collection and physical measurements

Data were collected at 150 K using an APEXII Bruker-AXS diffractometer equipped with Graphite monochromator and MoKα radiation (λ = 0.71073 Å). Absorption corrections were performed using the multi-scan technique using the SADABS program [13]. The total number of measured reflections was 7344 among which 2013, were independent and 1985 had the intensity I >2σ(I). The structure was solved with SHELXT [14], which revealed the position of all non-hydrogen atoms, and then refined with full-matrix least-square methods based on F 2 (SHELXL) [START_REF] Sheldrick | Crystal structure refinement with SHELXL[END_REF] included in the WINGX program [START_REF] Farrugia | WinGX and ORTEP for windows: an update[END_REF]. A final refinement on F 2 converged at R(F 2 ) = 0.044 and wR(F 2 ) = 0.114. A summary of the crystal data and the structure refinement parameters are included in Table 1. An ORTEP drawing of the asymmetric unit is exposed in Fig. 1a.

IR spectrum was registered at room temperature in the frequency range 4000-400 cm -1 with a Perkin-Elmer FT-IR 1000 spectrometerequipped with a diamond micro-ATR access.

Computational details

ACCEPTED MANUSCRIPT

Initially, thex-ray structure of (C 6 H 14 N) 2 [CdCl 4 ] compound imported on a CIF format was optimized in the gas phase with the Gaussian 09 software package [17].

All the quantum-chemical calculations have been performed via the hybrid B3LYP (Beckethree-parameter hybrid exchange functional combined with the Lee-Yang-Parr correlation functional) [18, 19] levels, the 6-311++G(d,p) basis set was used for all atoms except for the cadmium atom, the LANL2DZ basis set and effective core potentials (ECPs) (LANL2DZ basis and ECP built-in) have been used to represent the metal. Vibrational frequencies were calculated at the same level to ensure that each stationary point was a true minimum. We choose this functional because experimental and theoretical values are in good agreement. We have used the same method in similar compounds and have found good results [START_REF] Jomaa | C 8 H 14 N 2 )[END_REF][START_REF] Tahenti | A tetrachlorocobaltate(II) salt with 2-amino-5-picolinium: synthesis, theoretical and experimental characterization[END_REF].

On the other hand, the hybrid exchange-correlation functional named CAM-B3LYP combines the hybrid qualities of B3LYP and the long-range correction. Many theoretical studies show that CAM-B3LYP yields atomization energies of similar quality to those from B3LYP, while also performing well for charge transfer excitations in some models like in the case of organic compounds, which B3LYP underestimates enormously. In our case, in this paper, we study a hybrid compound which is formed by an organic and inorganic entity. For this type of these compounds, in our previous works [21,22], we have tried many methods, studied the dispersion effect and we have found that this functional haven't any effect in the energies or have very negligible effect.

Results and discussion

X-ray diffraction and geometry optimization

Single crystal X-ray diffraction analysis reveals that the asymmetric unit of the title material is formed by a monoprotonated cationic entity which is the cyclohexylammonium (C 6 H 11 NH 3 + ) and an anionic entity which is [CdCl 3 ] -, where Cd, Cl2 and Cl3 are in special positions in the mirror x = 0 and x = 1/2. By symmetry, the pyramid with a square base is

ACCEPTED MANUSCRIPT

formed. Otherwise, the organic entity exhibits a disorder manifested by the thermal agitation of all hydrogen, nitrogen, and carbon atoms.

The optimized molecular structure of (C 6 H 14 N) 2 CdCl 4 compound is given in Fig. 1b.

A comparison between selected optimized geometrical parameters, with experimental X-ray crystallographic parameters, is given in Table 2. The slight discrepancies observed between the computed and the experimental bond lengths and angles can be explained by the fact that the calculations relate to an isolated molecule, whereas the experimental results correspond to interacting molecules in the crystal lattice.

Interatomic
Further, we can propose a simple angular parameter used for penta-coordinate materials and defined by Alvarez

ACCEPTED MANUSCRIPT

So we can admit that the square pyramid CdCl 5 is slightly distorted.

The inorganic network consists of 1-D polymeric chains of square-pyramidal CdCl 5 vertex-sharing units developing parallel to the (b , c ) plane. The cationic layers are regularly intercalated between these anionic chains, the cohesion of the cations with the [CdCl 5 ] n chains is ensured through N-H…Cl hydrogen bonds resulting in the stabilization of the hybrid structure (Fig. 2a). A perspective view of the inorganic part (Fig. 2b), shows that the onedimensional polymeric chains of square-pyramidal CdCl 5 units are located around x = 0 and 1/2. These entities are linked to each other via Cl2 and Cl2 ii chlorine atoms (ii:

x, -y+1, z-1/2)
by vertex-sharing to form zigzag chains (Fig. 2c). Similar vertex-sharing example was also observed in (C 5 N 2 H 9 )CdCl 3 [START_REF] Chu | An ABX3 organic-inorganic perovskitetype material with the formula (C 5 N 2 H 9 )CdCl 3 : Application for detection of volatile organic solvent molecules[END_REF].

Two neighboring chains are further connected into a 2D supramolecular network by means of N-H… Cl type bonds with donor-acceptor distances ranging from 3.199 (5) to 3.372 ( 6) Å (Details in Table 3), three of which are bifurcated: N1A-H1A1…(Cl2 ii , Cl3 i ), N1B-

H1B2…(Cl1 iii , Cl3 iii ) (iii: x, -y+1, z+1/2) and N1B-H1B3…(Cl1, Cl2 ii ) (ii: x, -y+1, z-1/2), forming 𝑅 6 2 (8)

Hirshfeld surface analysis

Here, the Hirshfeld surface (HS) and their associated 2D fingerprint plots analysis of the asymmetric unit was undertaken as a systematic computational investigation method, that provides insight, quantification and exploration of intermolecular close contacts present within the crystal packing.The HS mapped over d norm of (C 6 H 14 N) 2 CdCl 4 illustrated in Fig. 3, allows a surface with a color scheme to be displayed, where dark red spots highlight the shortest intermolecular contacts that are attributed to N-H…Cl interactions, the white areas mark the contacts around the van der Waals separation and the blue areas illustrate the longest intermolecular contacts in the structure.

Further, the intermolecular interactions involved within the structures are also visible in the full 2D fingerprint (Fig. 4a) that can be decomposed (Fig. 4b,c,d a low contribution on the HS (0.4% and 1.1%). Despite the lowness of these percentages, we can suggest animprovement of the crystalline stability owing to their cumulative effect.

Such visual analysis of intermolecular interactions is consistent with those found by X-ray diffraction analysis.

Vibrational analysis

The vibrational properties of the (C 6 H 14 N) 2 [CdCl 4 ]hybrid material were discussed based on a vibrational study using experimental and calculated infrared spectroscopy. The experimental IR spectrum measured between 400 and 4000 cm -1 at room temperature and simulated using B3LYP/LanL2DZ basis set are shown in Fig. 5. For vibrational frequency, uniform scaling factors were adopted: frequencies below 1700 are multiplied by 0.983, so that greater than 1700 are multiplied by 0.958 [START_REF] Jomaa | C 8 H 14 N 2 )[END_REF][28][29][30][31]. This correction is made to find a good

ACCEPTED MANUSCRIPT

agreement with the experimental frequencies. The correlation between observed and theoretical spectra is reasonable; the differences observed especially in the region of higher wavenumbers can be easily attributed to the calculations because these were performed in the gas phase where the crystalline packing forces were not considered. The detailed assignment of the characteristic bands is essentially based on comparisons with data previously reported for similar materials. Brief discussions of the assignments for the most important functional groups are presented below.

[CdCl 4 2-] vibrations

According to the theoretical calculations and by reference to other compounds including the chlorocadmate anion [11,12], it is found that the inorganic entity leads to vibrations below 500 cm -1 . The asymmetric and symmetric stretching of Cd-Cl appeared at 189 and 174 cm -1 , respectively. The band corresponding to the bending mode of Cl-Cd-Cl appeared in the 74 cm -1 .

C 6 H 11 NH 3 + vibrations

The N-H asymmetric and symmetric stretching vibrations are observed in the highfrequency region between 3444 and 3119 cm -1 . The DFT computations give the frequency of these bands between 3495 and 3104 cm -1 .The bands observed at 1598 and 1565 cm -1 are respectively assigned to the asymmetric and symmetric bending modeof N-H.The corresponding bands are calculated to be found at 1605 and 1563 cm -1 . The bandsof medium intensities between 1486-1454 cm -1 are attributed to the characteristic NH 3 + asymmetric and symmetric stretching vibrations and predicted at 1503 and 1422 cm -1 . In cyclohexylammonium perchlorate [21] these modes are observed nearly in the same region.

The asymmetric and symmetric stretching vibrations of (-CH 2 ) are observed at 2936 and 2858 cm -1 respectively. Whereas, these vibrations are calculated at ν as = 2921 cm -1 and ν s = 2880 cm -1 . The bands that lie between1121 and 1007 cm -1 in the experimental spectrum correspond to the C-N stretching. The same vibration is calculated between 1145and 1003 cm - 1 by DFT level. Withal, the bands that appear around 547 and 455 cm -1 are attributed to the out of plane vibration modes of the δ(C-C) et δ(C-N) groups and predicted to be found in region 565-476 cm -1 .

Quantum theory of AIMs and RDG analysis

ACCEPTED MANUSCRIPT

R. Bader [START_REF] Bader | Atoms in Molecules, A Quantum Theory[END_REF] suggests that the presence of BCPs and the bonding paths between two atoms confirm the interaction between them. The graphical representation of the AIM analysis of the (C 6 H 14 N) 2 CdCl 4 is illustrated in Fig. 6, here the BCPs are located along the Cl… H bonds.

AIM analysis reveals that the title material is stabilized by eight hydrogen bonds that meet Koch and Popelier criteria [32,33] and that the interactions Cl47… H9, Cl47… H2, Cl44… H21, Cl44… H4 exhibit topological properties similar to the interactions Cl47… H30, Cl47… H23, Cl46… H42, Cl46… H25. The corresponding geometric, energetic, and topological parameters of the BCPs are collated in Table S2. With reference to Rozas criteria indicate as well the depletion of electronic charge along the binding path. The observed density for the Cl44… H4 or Cl46… H25 interactions (0.019515-0.019514 a.u.) is the highest probably because the distance between these two atoms is the smallest. Thus these last two interactions are considered to be the strongest according to the value of the calculated hydrogen bond energy which is of the order of -13.91126 kJ.mol -1 .

Further, the non-covalent interactions (NCIs) of a stabilizing or destabilizing nature existing in the (C 6 H 14 N) 2 CdCl 4 compound were also characterized by the NCI plot computational tool (RDG) that permits to assess the self-assembled structure. The noncovalent interactions are represented by isosurfaces instead of critical points.The illustration of various interactions types in the (C 6 H 14 N) 2 CdCl 4 structure is presented in Figs. 7a andb.

According to Fig. 7a, the attractive, van der Waals and repulsive interactions appear as blue, green, and red scattered spikes, respectively. Each individual peak corresponds to a single non-covalent interaction. The RDG peaks towards -0.02, -0.025 and -0.045 u.a with sign (λ2)ρ correspond to a strong attractive interaction confirmed also by the low-density values.

Regarding the isosurface color scheme in Fig. 7b, light blue spots are observed between the hydrogen and chlorine atoms, which justify the strong attractive interaction NH ... Cl. The differences in hydrogen bond strength can be seen from the depth of the blue, the absence of the strong blue colored isosurface clearly indicates that the structure of (C 6 H 14 N) 2 CdCl 4 does not exhibit strong hydrogen bonds.Interactions that appear as green patches around the anion are attributed to van der Waals interactions. The elliptical red plaque

ACCEPTED MANUSCRIPT

located at the center of the cycle is linked to repulsive interactions that show a strong steric effect. These results are consistent with those of the AIM analysis.

Molecular electrostatic potential (MEP)

The electrostaticpotential surface is relevant to electrondensitytherebyitisusedto predictthe reaction sites or attackswithpotentialelectrophiles or nucleophiles as well as hydrogen bonding interactions. The electrostatic potential at any point 𝑟 in the space surrounding a molecule can be expressed as [START_REF] Politzer | The fundamental nature and role of the electrostatic potential in atoms and molecules[END_REF]:

V(r)= 𝑍 𝐴 𝑅 𝐴 -𝑟 𝐴 - 𝝆 (𝒓 ′ )d𝑟 ′ 𝑟 ′ -𝑟
Where Z A is the charge of the nucleus A located at 𝑅 𝐴 and 𝜌(𝑟′) is the molecular electron density function. The sign of V(𝑟 ) at a particular region depends upon whether the effect of the nucleus or the electronsis dominant there :

-If V(𝑟 ) > 0 : region of strong attraction, electrophilic sites (blue color).

-If V(𝑟 ) < 0 : region of strong repulsion, nucleophilic sites (red color).

-Neutral and slightlyneutralpotential (green and yellowcolorrespectively).

The electronic distribution can be easily seen from the different colorations of their mapped surfaces, Fig. 8 illustrates the mapped MEP surfaces of (C 6 H 14 N) 2 [CdCl 4 ]compound.

The area of negative potentials (excess in electrons) defined by a red color is associated with [CdCl 4 ] 2-group as a nucleophilic site while the most positive potential (deficient in electrons) defined by a blue color is observed at the level of the hydrogen atoms, as expected due to its higher donor capacity, especially near the -NH 3 + group of the cyclohexylammonium cation so candidate for nucleophilic attacks.

On the cations side chain a notable yellow color typical of slightly inert region is observed. This confirms the existence of intermolecular N-H•••Cl interactions between the two organic cations and the inorganic anion in the compound understudy. The electrophilic and nucleophilic sites explain the formation of Hydrogen bonds, which is defined as being a donor-acceptor interaction. It demonstrates that the electrostatic potential of the donor becomes less negative as the acceptor becomes more negative [36].

Frontier orbitals and chemical quantum descriptors

The 3D plots of HOMO-1, HOMO, LUMO, and LUMO+1 of (C 6 H 14 N) 2 [CdCl 4 ] are provided in Fig. 9, the red and the green colors of the plots correspondto the positive and the negative phase, respectively. Orbital analysis revealed that the HOMO components are ACCEPTED MANUSCRIPT located on the anionic part, while the LUMO components are mainly located on the NH 3 + group of the organic cation and on the [CdCl 4 ] group. The energy level of the HOMO orbital is equal to-5.14eV and the energy level of the LUMO orbital is equal to -1.51eV, which gives an energy gap of about 3.63 eV. This value agrees very well with the corresponding energy gap reported earlier of other organic chlorocadmate(II) compounds [START_REF] Jomaa | C 8 H 14 N 2 )[END_REF].

Through the results of the calculated HOMO and LUMO energies and by means of equations reported in the literature [6,8], the global chemical reactivity descriptors (The electronegativity ( χ), hardness ( η), softness (S), electrophilicity index (ψ) and the chemical potential (μ)) can be defined. The results of this surveyusing the TD-DFT method are gathered in Table 4.

The energies of HOMO and LUMO orbitals of the studied compound are related to the ionization potential I and to the electronic affinity A, respectivelyby the following relationships: I = -E HOMO and A = -E LUMO . All the determined energies in the gas phase of (C 6 H 14 N) 2 [CdCl 4 ] material were found to be negative as well as the chemical potential.

Hardness and chemical softness are also good indicators of chemical stability; here the material has hardness equal to 1.81eV and chemical softness equal to 0.27 eV. The electrophile index ψ = 3.05 eV is relatively high, so we can conclude that these findings are an index of high kinetic stability and low chemical reactivity of the crystal structure.

Non-linear optical properties

The need for the investigation of hybrids was to promote new sophisticated materials with targeted physical and chemical features [START_REF] Jomaa | C 8 H 14 N 2 )[END_REF], in the NLO application field, organicinorganic complex crystals have been developed to combine the advantages of inorganic crystals, such as good stability and resistance to physical andchemical attack, with the advantages of organic crystals, in order to retain the high NLO efficiency of organic materials and overcome its shortcoming.

The electric dipole moment, polarizability, and first-order hyperpolarizability values of the title materialusing DFT/B3LYP/LANL2DZ method are grouped in Table 5. (α) and (β) values of the Gaussian output are in atomic units (a.u.), so they have been converted into electronic units (esu) (α; 1 a.u. = 0.1482×10 -24 esu, β; 1 a.u. = 8.6393×10 -33 esu).

The quantities cited (Table 5) are obtained from the following equations:

μ= (μ x 2 + μ y 2 + μ z 2 ) 1/2 α tot = 1 3 (α xx + α yy + α zz ) ACCEPTED MANUSCRIPT β= [(β xxx + β xyy + β xzz ) 2 + (β yyy + β yzz + β yxx ) 2 + (β zzz + β zxx + β zyy ) 2 ] 1/2
The calculated values of the dipole moment, the polarizability and the first-order hyperpolarizability of (C 6 H 14 N) 2 CdCl 4 compound are respectively 10.15 Debye, 35.75×10 -24 esu, and 2776.41×10 -33 esu. The comparison of these values with the urea which is a standard NLO reference material (μ= 1.3732 Debye, α= 5.0477×10 -24 esu, β= 372.8×10 -33 esu) shows a large value of β tot which is a bit higher than7x β urea indicating important second-order optical effects. According to the literature [12], hydrogen bonds play an important role not only in the generation of the crystal structure and its stability, but also they contribute to the improvement of the hyperpolarizability of molecular systems and therefore the enhancement of the optical properties of the material.

The dipole moment is 7 times higher than μ urea which can be explained by the noncentrosymmetric nature of this molecule. The highest absolute value of the dipole moment is observed for the component μ z in this direction the value is equal to 8.61 D, while the polarizability is mainly along the x-axis (α xx = 286.790 a.u.). The calculated polarizability α tot is found to be higher than 7α urea . All these results nominate the (C 6 H 14 N) 2 CdCl 4 as a good material for NLO application.

Conclusion

This contribution is intended to cast light on the non-covalent weak interactions and on the ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

D-H•••A D-H (Å) H•••A (Å) D•••A (Å) D-H•••A (°) N1A-H1A1

catalysis [ 1 ]

 1 , biochemistry [2], magnetism, electrical conductivity [3], and nonlinear optical (NLO) devices [4].

  Hirshfeld surfaces mapped with d norm and their associated 2D fingerprint were plotted using Crystal Explorer version 3.1 software [23]. The topological analysis was carried out with the Multiwfn program [24] while the reduced density gradient of the title compound is graphed by Multiwfn and plotted by the VMD program [25]. The optimized molecular structure is used to perform HOMO-LUMO distributions and MEP analysis.

  β are the two largest Cl-Cd-Cl bond angles, with [β<α].Thus, any structure intermediate between trigonal bipyramid and square pyramid must have an angular parameter in the range0 ≤τ≤ 1, a perfect trigonal bipyramid presents (α = 180 and β = 120) resulting in τ= 1 while an ideal square pyramid has α = β and its angular parameter is τ = 0.

Fig. 5 f

 5 Fig.5 f, shows the pairs (d i , d e ) of chlorine atoms. It can be seen that (d e + d i ~ 3.6 Å) > 3.5 Å, the sum of the VdW radii of the chlorine atoms. This reveals the absence of close Cl ... Cl contacts.

[ 34 ]

 34 , all of the eight hydrogen bonds are considered weak since the Laplacian and the energy density values are positive. Positive Laplacian values that vary from 0.05617 to 0.00874 a.u.

  NLO properties of the (C 6 H 14 N) 2 CdCl 4 material. Structural analysis has shown that noncovalent assembly is mainly controlled by N-H…Cl H-bonds. The entire supramolecular self-assembly and consequent cooperatively of weak noncovalent interactions have been analyzed in detail (AIM, RDG, and MEP analysis) and further quantified by theoretical Hirshfeld surface analysis. The theoretical investigations support the experimental findings.The HOMO and LUMO energies and other calculated quantum parameters reveal the hardness and the great stability of the material. The assignment of the vibrational bands wasprecisely performed by comparison with the vibration modes frequencies of homologous compounds and with the theoretical results. The compound is a good nonlinear optical material candidate.

Fig. 1 .

 1 Fig. 1. ORTEP representation of the asymmetric unit of (C 6 H 14 N) 2 [CdCl 4 ] with atom-labeling scheme. Displacement ellipsoids are drawn at the 50% probability level. H atoms are represented as small spheres of arbitrary radii. Hydrogen bonds are denoted as dashed lines (a) and the optimized molecular structure (b).

Fig. 2 .

 2 Fig. 2. The packing of the title compound viewed along the c-axis (a), while (b,c,d,e) represents a perspective view showing the disposition of the 1-D polymeric chains of the inorganic part and the patterns of hydrogen bonding in (C 6 H 14 N) 2 [CdCl 4 ], CdCl 5 is given in polyhedral representation. The dotted lines indicate the hydrogen bonds.

Fig. 3 .

 3 Fig. 3. 3D Hirshfeld surfaces mapped with d norm in the range [-0.806 to 1.173] of (C 6 H 14 N) 2 [CdCl 4 ], red dotted lines represent hydrogen bonds. (iii) x, -y, z-1/2.

Fig. 4 .

 4 Fig. 4. 2D fingerprint plots of the full (a) and the main intercontacts showing percentage of various intermolecular contacts contributed to the Hirshfeld surface (C 6 H 14 N) 2 [CdCl 4 ](b,c,d,e,f). Surface to the side (g) highlight the relevant surface patch associated with H…Cl/Cl…H contacts.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Theoretical and experimental IR spectrum of (C 6 H 14 N) 2 CdCl 4 .

Fig. 8 .

 8 Fig. 8. Molecular electrostatic potential surface (MEPs) of (C 6 H 14 N) 2 [CdCl 4 ].

Fig. 9 .

 9 Fig. 9. The frontier molecular orbitals of (C 6 H 14 N) 2 CdCl 4 molecule at DFT/LANL2DZ level.

Fig. 1 .

 1 Fig. 1. ORTEP representation of the asymmetric unit of (C 6 H 14 N) 2 [CdCl 4 ] with atom-labeling scheme. Displacement ellipsoids are drawn at the 50% probability level. H atoms are represented as small spheres of arbitrary radii. Hydrogen bonds are denoted as dashed lines (a) and the optimized molecular structure (b).

Fig. 2 .

 2 Fig. 2. The packing of the title compound viewed along the c-axis (a), while (b,c,d,e) represents a perspective view showing the disposition of the 1-D polymeric chains of the inorganic part and the patterns of hydrogen bonding in (C 6 H 14 N) 2 [CdCl 4 ], CdCl 5 is given in polyhedral representation. The dotted lines indicate the hydrogen bonds.

Fig. 3 .

 3 Fig. 3. 3D Hirshfeld surfaces mapped with d norm in the range [-0.806 to 1.173] of (C 6 H 14 N) 2 [CdCl 4 ], red dotted lines represent hydrogen bonds. (iii) x, -y, z-1/2.

Fig. 4 .

 4 Fig. 4. 2D fingerprint plots of the full (a) and the main intercontacts showing percentage of various intermolecular contacts contributed to the Hirshfeld surface (C 6 H 14 N) 2 [CdCl 4 ](b,c,d,e,f). Surface to the side (g) highlight the relevant surface patch associated with H…Cl/Cl…H contacts.

Fig. 5 .

 5 Fig. 5. Theoretical and experimental IR spectrum of (C 6 H 14 N) 2 CdCl 4 .

Fig. 6 .

 6 Fig. 6. AIM molecular graph screening the different bond critical points (BCPs) of (C 6 H 14 N) 2 CdCl 4 . The BCPs are represented by small orange balls and the connecting paths by orange lines.

Fig. 7 .

 7 Fig. 7. Plots of the reduced density gradient scatter (a), and color-filled isosurface of the (C 6 H 14 N) 2 [CdCl 4 ] complex illustrating the non-bonded interactions (b).

Fig. 8 .

 8 Fig. 8. Molecular electrostatic potential surface (MEPs) of (C 6 H 14 N) 2 [CdCl 4 ].

Fig. 9 .

 9 Fig. 9. The frontier molecular orbitals of (C 6 H 14 N) 2 CdCl 4 molecule at DFT/LANL2DZ level.

  single broad peak in the middle of the scattered points, which point towards the origin and correspond to (d e + d i ~ 1.3 Å), a value less than the sum of the VdW radii of the atoms of hydrogen (1.09 + 1.09 = 2.18 Å). The Cd… Cl / Cl… Cd contacts which occupy 3.7% of the entire HS show on its 2D graph the presence of a pair of symmetrical wings (d e = d i ~ 1.4 Å) value less than the sum of the VdW radii of the atoms involved (Cd: 1.58 Å Cl: 1.75 Å).

	These findings affirm that the last three intercontact are considered to be close
	contacts.
	The H… Cd / Cd… H contacts represent 1.1% of the total HS and which correspond
	to d e + d i ~ 3.6 Å> (2.67 Å = 1.09 + 1.58). These contacts are considered non-close contact.
	,e,f) to quantify
	individual contributions and highlight pairs of atoms in close contact.
	In fact, As expected, the major contribution (55.1%) comes from H…Cl/Cl…H

intercontacts that appear as two pointed symmetrical spikes with a maximum sum of d e + d i ~ 2.2 Å less than the sum of the VdW radii of the chlorine atom and the hydrogen atom (1.75 + 1.09 = 2.84 Å). H… H contacts comprise 39.7% of the entire HS. They are characterized by ACCEPTED MANUSCRIPT a

Table S1 .

 S1 The value (ER

HCl = 1.36) greater than unity, of H… Cl / Cl… H contacts show the propensity to form N-H ... Cl hydrogen bonds. H ... H contacts appear with an enrichment ratio close to unity (0.86), this value agrees well with Jelch's expectation [27]. These contacts occupy more than half of HS with S H = 67.8%. On the other hand, the Cd ... Cl / Cl ... Cd contacts are extremely enriched (ER = 2.59), they represent a major attraction in the crystal despite their small surface area. The Cl… Cl, and H… Cd / Cd… H contacts represent as well

Table captions Table 1 .

 captions1 Crystal data and structure refinement parameters for (C 6 H 14 N) 2 CdCl 4 .

Table 2 .

 2 Principal intermolecular distances (Å) and bond angles (°) in (C 6 H 14 N) 2 CdCl 4 by Xray data (with estimated standard deviation in parentheses) and by theoretical calculations.

Table 3 .

 3 Hydrogen-bonds geometry (Å, °) of (C 6 H 14 N) 2 CdCl 4 .

Table 4 .

 4 Global reactivity descriptors and calculated energy values of (C 6 H 14 N) 2 CdCl 4 at DFT/LANL2DZ level.

Table 5 .

 5 The electric dipole moment, polarizability, and first order hyperpolarizability values of (C 6 H 14 N) 2 CdCl 4 .

Table 1 .

 1 Crystal data and structure refinement parameters for (C 6 H 14 N) 2 CdCl 4 .

	CCDC

Table 3 .

 3 Hydrogen-bonds geometry (Å, °) of (C 6 H 14 N) 2 CdCl 4 .

Table 4 .

 4 Global reactivity descriptors and calculated frontier molecular orbital parameters of (C 6 H 14 N) 2 CdCl 4 using DFT/LANL2DZ level.

	•••Cl3 i	0.91	2.51	3.199 (5)	133
	N1A-H1A1•••Cl2 ii	0.91	2.82	3.372 (6)	121
	N1A-H1A2•••Cl1 iii	0.91	2.34	3.248 (6)	176
	N1A-H1A3•••Cl1	0.91	2.34	3.223 (5)	165
	N1B-H1B1•••Cl3 i	0.91	2.31	3.199 (5)	164
	N1B-H1B2•••Cl1 iii	0.91	2.75	3.248 (6)	116
	N1B-H1B2•••Cl3 iii	0.91	2.79	3.289 (5)	116
	N1B-H1B3•••Cl1	0.91	2.58	3.223 (5)	128
	N1B-H1B3•••Cl2 ii	0.91	2.65	3.372 (6)	137
	Symmetry codes: (i) x, y, z-1; (ii) x, -y+1, z-1/2; (iii) x, -y, z-1/2.
	Energies (eV)		Ionization potential A	1.51594704
	Frontier orbitals		Electron affinity I		5.14975726
	E HOMO	-5.14975726	Descriptors
	E LUMO	-1.51594704	Global electrophilicity (ω)	3.05681992
	|E HOMO -E LUMO |	3.63381022	Electronegativity (χ)	3.33285215
	E HOMO-1	-5.32663135	Global chemical potential (μ) -3.33285215
	E LUMO+1	-1.16764113	Global hardness (η)	1.81690511
	|E HOMO-1 -E LUMO+1 |	4.15899022	Global softness (S)		0.27519324

Acknowledgements

We are grateful to the Tunisian Ministry of Higher Education Scientific Research for the provided support. Supplementary Information: Tables S1-S2. Crystallographic data for the structural analysis have been deposited at the Cambridge Crystallographic Data Centre, CCDC No 2068977. These data can be obtained free of charge via https: //www.ccdc.cam.ac.uk/conts /retri eving.html, or from the CCDC, 12 Union Road, Cambridge, CB2 1EZ, UK: fax: (+44) 01223-336-033; e-mail: deposit@ccdc.cam.ac.

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT


Thierry Gonon

Céline Helbert

Christophette Blanchet-Scalliet

Bruno Demory

Gaussian process regression on nested spaces

Keywords: Metamodel, Kriging, Gaussian process regression, High dimension, Infinite conditioning, Multifidelity, Nested spaces, Variable-size design space problems

Metamodels are widely used in industry to predict the output of an expensive computer code. As industrial computer codes involve a large amount of input variables, creating directly one big metamodel depending on the whole set of inputs may be a very challenging problem. Industrialists choose instead to proceed sequentially. They build metamodels depending on nested sets of variables (the variables that are set aside are fixed to nominal values), i.e. the dimension of the input space is progressively increased. However, at each step, the previous piece of information is lost as a new Design of Experiment (DoE) is generated to learn the new metamodel. In this paper, an alternative approach will be introduced, based on all the DoEs rather than just the last one. This metamodel uses Gaussian process regression and is called seqGPR (sequential Gaussian process regression). At each step n, the output is supposed to be the realization of the sum of two independent Gaussian processes Y n-1 + Z n . The first one Y n-1 models the output at step n -1. It is defined on the input space of step n -1 which is a subspace of the one of step n. The second Gaussian process Z n is a correction term defined on the input space of step n. It represents the additional information provided by the newly released variables. Z n has the particularity of being null on the subspace where Y n-1 is defined so that there is a coherence between the steps. Firstly, some candidate Gaussian processes for (Z n ) n 2 are suggested, which have the property of being null on an infinite continuous set of points. Then, an EM (Expectation-Maximization) algorithm is implemented to estimate the parameters of the processes. Finally, the metamodel seqGPR is compared to a classic kriging metamodel where the output is assumed to be the realization of one second order stationary Gaussian process. The comparison is made on two analytic examples, a first one with two steps, up to dimension 4, and a second one with three steps, up to dimension 15. The introduced methodology is also tested on an industrial example which goes from dimension 11 to dimension 15. In all these test cases, seqGPR performs better than, or at least as well as kriging.

Introduction

In industry, numerical codes are widely used to model physical phenomena involved in manufacturing (see section 5.3 for an example of industrial product). Computer experiments make the search for performance easier and cheaper than physical experiments. However, computer codes are confronted to some limitations. As they model complex physical phenomena, they are often computationally expensive. This problem is solved by the use of metamodels [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF], [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. A metamodel is a simpler statistical model, like radial basis neural network [START_REF] Cheng | Neural networks: A review from a statistical perspective[END_REF], kriging model [START_REF] Santner | The design and analysis of computer experiments[END_REF], [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF], support vector regression [START_REF] Clarke | Analysis of support vector regression for approximation of complex engineering analyses[END_REF]. This simpler model is fit on a sample of well-chosen runs, also called design of experiments (DoE) [START_REF] Pronzato | Design of computer experiments: space filling and beyond[END_REF], [START_REF] Dupuy | Dicedesign and diceeval: Two r packages for design and analysis of computer experiments[END_REF]. Metamodels provide fast approximations of the code outputs. Besides, they have a proactive role as they help to select relevant simulations to be run [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]]. The computer code is then used more efficiently.

Building surrogate models can be difficult when numerical codes involve lots of variables: geometrical parameters, environmental ones etc. Complex engineering systems can be found in [START_REF] Makowski | Global sensitivity analysis for calculating the contribution of genetic parameters to the variance of crop model prediction[END_REF], [START_REF] Lefebvre | A methodological approach for statistical evaluation of aircraft infrared signature[END_REF], [START_REF] Auder | Screening and metamodeling of computer experiments with functional outputs. application to thermal-hydraulic computations[END_REF]. Dimension reduction is a common approach to deal with it. It consists in considering only the influential inputs. The influential inputs are often detected by sensitivity analysis methods ( [Saltelli, 2000], [Sobol, 1993], [START_REF] Iooss | Shapley effects for sensitivity analysis with correlated inputs: comparisons with sobol'indices, numerical estimation and applications[END_REF]). However, this approach leads to a vicious circle. Indeed, on one side, the variables selected by sensitivity analysis are relevant only if the metamodel is sufficiently accurate, and on the other side, the dimension reduction is used to increase the metamodel accuracy. An alternative method that is usually chosen in industry is to act step by step. The first studies carried out on the code have focused on a small amount of variables deemed as important based on expert knowledge. The other variables are set to nominal values. A metamodel is created and used only for these variables until the understanding of their influence is sufficient. Then, to obtain better performances, finer studies are proceeded with, in which new variables are added progressively. The classical approach is to build a metamodel on a new independent space-filling design in the wider input space. Yet, this approach may be criticized as the information provided by the previous DoE's and metamodels is disregarded.

In the present work, an original metamodel is introduced, which takes into account the information from all the steps. The chosen approach consists in creating dependent metamodels from one step to another. The metamodels are based on Gaussian process regression [START_REF] Santner | The design and analysis of computer experiments[END_REF] [ [START_REF] Williams | Gaussian processes for machine learning[END_REF]. One way of taking into account the different designs could be to use a multifidelity model for which the levels are not defined on the same input space. That is what is done in [START_REF] Hebbal | Multifidelity modeling with different input domain definitions using deep gaussian processes[END_REF] with deep Gaussian processes. Nevertheless, this model is really expensive to learn. Furthermore, the multifidelity context is not suitable since the runs obtained at the previous steps have the same accuracy level. Another way could be to define a virtual categorical input, equal to the number of the step, that influences the choice of the input variables to consider. [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF] defines a Gaussian process whose input space varies dynamically with this virtual categorical input. This modelization does not take into account the fact that the input sets are entwined from one step to another. This particularity can lead to a simpler modelization. The approach in this work, called seqGPR (sequential Gaussian process regression), is based on a recursive statistical model inspired by the autoregressive multifidelity model [Kennedy and O'Hagan, 2000] but with the same accuracy for all the runs and with an input space of increasing dimension. The output of the code at a given step is modeled as the realization of a Gaussian process being the sum of the Gaussian process that models the previous step and a correction term. The correction term stands for the information provided by the new variables.

As regards this probabilistic model, one first difficulty is to build an appropriate correction term to represent the gap between the same numerical model observed on two nested subspaces. Thus, the correction process must be null on the subspace corresponding to the previous step. As this subspace is composed of an infinite continuous set of points, the question is then how to build a Gaussian process null on an infinite continuous set of points with a numerically computable covariance kernel. One idea could be to enjoin the prediction to verify the nullity property. For example, relevant points from the concerned subspace can be added to the training sample. This technique is used in [Da Veiga and Marrel, 2012] in the context of monotonicity, boundedness, convexity constraints. One drawback is the greediness of that method. Besides, the nullity cannot be verified in the whole subspace. Instead of establishing the nullity in retrospect, it can be verified at first glance, as an intrinsic property of the correction term. The correction process could be defined in a finite dimensional way, with adequate basis functions, as it is done in [START_REF] Maatouk | Gaussian process emulators for computer experiments with inequality constraints[END_REF], [Lopera, 2019], [START_REF] Bachoc | Sequential construction and dimension reduction of gaussian processes under inequality constraints[END_REF] to ensure some properties of monotonicity, boundedness, convexity. Yet, this approach is greedy as the number of basis functions increases with the dimension. In this paper, the Gaussian process introduced in [START_REF] Gauthier | Spectral approach for kernel-based interpolation[END_REF] which is an extension of the conditional expectation to an infinite continuous set of points is retained. A tractable kernel is sought for this process which is then compared to another candidate verifying the nullity condition.

Once a tractable candidate for the correction term is proposed, another difficulty lies in the estimation of the hyperparameters. One way is to optimize the likelihood. This task can be numerically difficult, since the likelihood involves several sets of hyperparameters for the different processes. One way to reduce the dimension of the optimization space is to propose nested designs as it is done in [Le Gratiet and Garnier, 2014]. In this paper an alternative based on the expectation maximization algorithm [START_REF] Friedman | The elements of statistical learning[END_REF] is introduced. It allows the reduction of the dimension with limited constraints on the designs.

The formulation of the metamodel seqGPR (sequential Gaussian process regression) is detailed in section 2. Two candidates for the correction term are proposed in section 3. An EM (expectation-maximization) algorithm is suggested in section 4 to estimate the model parameters. Finally, the method seqGPR is tested on two analytic examples and one industrial case in section 5.

Metamodel

Underlying processes Let f : [0, 1] d → R be an output of an expensive simulation code depending on d variables. A response surface of f has to be created taking into account that N -1 previous studies of increasing dimension have already been completed. Each of these steps is a focus on the relation between the output and a subset of free variables. The other variables are temporarily fixed but then progressively released in the further steps. This paper deals with the handling of the last step, denoted by N , where all inputs are free. Let I n denote the index set of the variables that are released at step n ∈ 1, N and d n = card(I n ). For every index set I, let x I denote the associated subvector of x. The following framework, based on Gaussian process regression, is introduced to model the successive studies:

= m + Z 1 of mean m ∈ R and with Z 1 : Ω × [0, 1] d1 → R, a centered Gaussian Process of covariance kernel σ 2 1 ρ 1 . Let y 1 = f 1 (X 1
) represent the vector of observations of the output on DoE X 1 .

• At step 2, a second range of simulations is then launched on a subspace of a higher dimension [0, 1] d1+d2 .

The

: Ω × [0, 1] d1+d2 → R,
which is the sum of the process at step 1 Y 1 , and a correction term Z 2 , which represents the additional information provided by the released variables. As f 1 and f 2 are restrictions of the same function f , Y 1 and Y 2 have to coincide on the subspace defined at step 1. This results in the following definition of

Y 2 : Y 2 (x I1 , x I2 ) = Y 1 (x I1 ) + Z 2 (x I1 , x I2 ), ∀(x I1 , x I2 ) ∈ [0, 1] d1+d2 ,
with Z 2 a centered Gaussian process independent from Z 1 , of covariance kernel σ 2 2 ρ 2 such that for all x I1 in [0, 1] d1 , Z 2 (x I1 , xI2 ) = 0. Let y 2 = f 2 (X 2 ) denote the vector of observations for this step.

• In the same way as previous steps, at step n ∈ 3, N , a DoE X n is created in the higher space [0, 1]

Y n (x I1∪•••∪In-1 , x In ) = Y n-1 (x I1∪•••∪In-1 ) + Z n (x I1∪•••∪In-1 , x In ), ∀(x I1∪•••In-1 , x In ) ∈ [0, 1] d1+•••+dn , with Z n a centered Gaussian process independent from (Z 1 , • • • , Z n-1 ) of covariance kernel σ 2 n ρ n and such that Z n (x I1∪•••∪In-1 , xIn ) = 0, ∀x I1∪•••∪In-1 ∈ [0, 1] d1+•••+dn-1 .
The coincidence of the (Y n ) 1 n N on the subspaces is illustrated on figure 1 for N = 2.

Metamodel seqGPR Finally, the problem is modeled by the following statistical model:

       Y 1 (x I1 ) = m + Z 1 (x I1 ), ∀x I1 ∈ [0, 1] d1 , Y n (x I1∪•••∪In-1 , x In ) = Y n-1 (x I1∪•••∪In-1 ) + Z n (x I1∪•••∪In-1 , x In ), ∀n ∈ 2, N , ∀x I1∪•••∪In ∈ [0, 1] d1+•••+dn , (1) where:
• The processes (Z n ) 1 n N are independent Gaussian processes of law:

Z n ∼ GP(0, σ 2 n ρ n (x I1∪•••∪In , t I1∪•••∪In )) ∀n ∈ 1, N .
(2) • The processes (Z n ) 2 n N verify the following property:

Z n (x I1∪•••In-1 , x In ) = 0, ∀n ∈ 2, N , ∀x I1∪•••In-1 ∈ [0, 1] d1+•••+dn-1 . (3)
The same formulas of prediction as for a classic kriging metamodel apply. For every x ∈ [0, 1] d , the prediction mean ŷ(x) and variance v(x) are defined by

ŷ(x) = E [Y N (x) | Y 1 (X 1 ) = y 1 , • • • , Y N (X N ) = y N ] , v(x) = V ar (Y N (x) | Y 1 (X 1 ) = y 1 , • • • , Y N (X N ) = y N ) .
Issues This model raises two issues. The first difficulty is to build (Z n ) 2 n N such that it verifies the nullity property (cf (3)). This issue is tackled in section 3. Secondly, the estimation of the parameters of Y 1 and (Z n ) 2 n N at step N , in the presence of training samples from the different steps, stands ahead as a thorny hurdle to overcome. It is detailed in section 4.

Candidates for the correction process

This section answers the first issue, which is to build a process that is null on an infinite continuous set of points. More precisely, the goal is to build a process Z :

[0, 1] d J +d I × Ω → R such that: Z(x J , g(x J )) = 0 ∀x J ∈ [0, 1] d J , (4) with g : [0, 1] d J → [0, 1] d I a deterministic function.
In this section, two candidates are suggested for the process Z, both based on a latent Gaussian process

∼ Z ∼ GP 0, σ 2 r((x J , x I ), (t J , t I )) of covariance kernel σ 2 r.

Red (reduced) process

A natural idea to obtain the nullity property (4) is to subtract from

∼ Z its value on the subspace of [0, 1] d J +d I equal to {(x J , g(x J )) | x J ∈ [0, 1] d J }. The first candidate is: Z Red (x J , x I ) = ∼ Z(x J , x I ) - ∼ Z(x J , g(x J )).
It is called the Red process. It is a centered Gaussian process of covariance kernel σ 2 ρ Red with:

ρ Red ((x J , x I ), (t J , t I )) = r ((x J , x I ), (t J , t I )) -r ((x J , x I ), (t J , g(t J )))
-r ((x J , g(x J )), (t J , t I )) + r ((x J , g(x J )), (t J , g(t J )) .

P (preconditioned) process

Another way to obtain the nullity property ( 4) is to use the "conditional expectation" [START_REF] Gauthier | Spectral approach for kernel-based interpolation[END_REF].

The second candidate is defined by

Z P (x J , x I ) = ∼ Z(x J , x I ) -E ∼ Z(x J , x I ) | ∼ Z(s J , g(s J )), ∀s J ∈ [0, 1] d J , ∀(x J , x I ) ∈ [0, 1] d J +d I .
(5)

It is called the P process. The term

E ∼ Z(x J , x I ) | ∼ Z(s J , g(s J )), ∀s J ∈ [0, 1] d J is the orthogonal projection of ∼ Z(x J , x I ) in the sub Gaussian space engendered by the ∼ Z(s J , g(s J )) s J ∈[0,1] d J
. Its expression is given by [START_REF] Gauthier | Spectral approach for kernel-based interpolation[END_REF]:

E ∼ Z(x J , x I ) | ∼ Z(s J , g(s J )), ∀s J ∈ [0, 1] d J = +∞ n=1 φ n (x J , x I ) [0,1] d J ∼ φ n (s J ) ∼ Z(s J , g(s J ))ds J , (6)
with:

φ n (x J , x I ) = 1 λn [0,1] d J σ 2 r((x J , x I ), (s J , g(s J ))) ∼ φ n (s J )ds J ,
and (λ n , ∼ φ n ) n 1 are solutions of the eigen problem:

[0,1] d J σ 2 r ((x J , g(x J )), (s J , g(s J ))) ∼ φ n (s J )ds J = λ n ∼ φ n (x J ), ∀x J ∈ [0, 1] d J ,
such that:

[0,1] d J ∼ φ n (s J ) ∼ φ m (s J )ds J = δ nm , ∀n, m 1,
where δ nm is the Kronecker symbol. Z P is a centered Gaussian process of covariance kernel:

σ 2 ρ P ((x J , x I ), (t J , t I )) = σ 2 r((x J , x I ), (t J , t I )) - +∞ n=1 λ n φ n (x J , x I )φ n (t J , t I ), ∀(x J , x I ) ∈ [0, 1] d J × [0, 1] d I , ∀(t J , t I ) ∈ [0, 1] d J × [0, 1] d I .
This kernel cannot be used just as it is in practice, because the solutions of the eigen problem are not explicit in general and the sums are infinite. In this paper, two ways of computing this kernel are proposed. The first is based on the discretization of the spectral decomposition. The second is based on a wise choice of r (the correlation kernel of the latent process ∼ Z) for which an explicit formula is known.

Numerical approximation

In this paragraph, a first way of computing the kernel is given by an approximation based on the discretization of the spectral decomposition that reduces the functional eigen problem to a finite dimensional one. In that case, the terms from the spectral decomposition vanish and the formula of the resulting approximate kernel only depends on the kernel of the latent process.

Proposition 1 Let D = t (i) J , g(t (i) J ) 1 i L be a uniform sample in the subspace {(t J , g(t J )), t J ∈ [0, 1] d J }.
Then, discretizing the spectral decomposition of the P process Z by a Monte-Carlo method using D is equivalent to approximating the P process by the process ∼ Z conditioned to be null on the points of D:

Z D = ∼ Z | ∼ Z(D) = 0 . Z D is a centered Gaussian process of covariance kernel σ 2 ρ D : ρ D (x, t) = r(x, t) -r(x, D)r(D, D) -1 r(D, t) ∀x, t ∈ [0, 1] d J +d I .
See appendix A for the proof of this proposition.

This process does not match the original expectations as the nullity is not fully respected on the continous set of points. The size of D must be high to correctly approach the full nullity. It is more and more difficult as the dimension of the input space increases, and it leads to very expensive computations (with inversion of huge matrices). Therefore, a second way of computing this kernel is proposed. It consists in finding forms of the kernel of ∼ Z that enable the kernel of Z P to be tractable. Analytical expression in the case of a product kernel Proposition 2 Let r J and r I be two stationary correlation kernels. If r is of the form:

r((x J , x I ), (t J , t I )) = r J (x J , t J )r I (x I -g(x J ), t I -g(t J )),
then the P process is a centered Gaussian process equal to:

Z P (x J , x I ) = ∼ Z(x J , x I ) -r I (x I -g(x J ), 0) ∼ Z(x J , g(x J )),
whose covariance kernel is σ 2 ρ P with:

ρ P ((x J , x I ), (t J , t I )) = r J (x J , t J ) [r I (x I -g(x J ), t I -g(t J )) -r I (x I -g(x J ), 0)r I (0, t I -g(t J ))] .
See appendix B for the proof of this proposition.

Corollary 1 If g is constant equal to c ∈ [0, 1] d I and if r is a stationary kernel of the form r((x J , x I ), (t J , t I )) = r J (x J , t J )r I (x I , t I ),

then the P process is a centered Gaussian process equal to:

Z P (x J , x I ) = ∼ Z(x J , x I ) -r I (x I , c) ∼ Z(x J , c),
whose covariance kernel is σ 2 ρ P with:

ρ P ((x J , x I ), (t J , t I )) = r J (x J , t J ) [r I (x I , t I ) -r I (x I , c)r I (c, t I )] .
The model presented in the corollary (where the released variables were previously fixed at a constant value) is illustrated in figure 2.

In the following, the P processes are assumed to be built as in proposition 2 or corollary 1.

Interpretation of the candidates

An illustration of the processes at stake in the construction of P and Red processes is shown in figure 3. They are used to build Z(x 1 , x 2 ) such that Z(x 1 , 0) = 0. One main difference between Red and P process is that the initial latent process is disturbed locally in the case of the P process, whereas it is disturbed globally in the case of the Red process. Indeed, ∼ Z(x 1 , 0) (see panel 3a), which is used to build the Red process, is constant in x 2 and consequently generates a disruption of ∼ Z(x 1 , x 2 ) (see panel 3b) in the whole input space [0, 1] 2 . On the contrary, the conditional expectation

E ∼ Z(x 1 , x 2 ) | ∼ Z(t 1 , 0)∀t 1 ∈ [0, 1] (see panel 3c), used
to build the P process, depends on x 2 , especially near the line x 2 = 0, and tends to 0 far from the line. The modification it implies on the latent process is therefore located almost only along the line. that should be null on the line x 2 = 0. On left panel, the process which is removed from the latent process to create the Red process. On right panel, the process (denoted by E for expectation) which is removed from the latent process to create the P process.

Estimation of the parameters

This section answers the second issue which is to estimate the parameters of the metamodel seqGPR defined in section 2. An EM algorithm is presented to ease the maximization of the likelihood.

Notations At step N , the following training samples (also called observed data) are available: (X 1 , y 1 ), ..., (X N , y N ). The training data can be written as:

         Y 1 (X 1 ) = y 1 , . . . Y N -1 (X N -1 ) = y N-1 , Y N (X N ) = y N ,
or more simply: Y = y, with:

Y =    Y 1 (X 1 ) . . . Y N (X N )    and y =    y 1 . . . y N    .
The law of Y 1 (see equations ( 1) and ( 2)) depends on the parameters: m (scalar mean parameter), σ 2 1 (scalar variance parameter), and θ 1 (vector of covariance parameters). For all n ∈ 2, N , Z n (see equation ( 2)) is either a Red or P process built on the latent process

∼ Z n ∼ GP(0, σ 2 n r n ) (see section 3) of parameters: σ 2
n (scalar variance parameter), θ n (vector of covariance parameters). In order to emphasize the dependence of ρ n (the covariance kernel of Z n , n ∈ 1, N ) on the parameter θ n , this kernel is denoted by ρ θn . The purpose of this section is to estimate the parameters η = (m,

σ 2 1 , θ 1 η1 , σ 2 2 , θ 2 η2 , • • • , σ 2 n , θ n ηn , • • • , σ 2 N , θ N η N ) based
on the observed data. The maximum likelihood estimator is used. The following loss function1 has to be minimized:

l(η; y) = log |Cov η (Y, Y)| + (y -E η [Y]) Cov η (Y, Y) -1 (y -E η [Y]),
with for any matrix M , M denoting the transpose of M . This optimization problem is complex as η can reach big dimensions.

In what follows, the notation X n ⊂ X n-1 means that the submatrix of X n composed of the columns corresponding to x I1∪•••∪In-1 is included in X n-1 . In this case, the designs X n and X n-1 are said nested.

l(η; y) = l 1 (η 1 ; y 1 ) + N n=2 l n (η n ; z n ),
with l 1 (η 1 ; y 1 ) (respectively l n (η n ; z n ), n ∈ 1, N ) the loss function associated with the training data Y 1 (X 1 ) = y 1 (respectively

Z n (X n ) = z n ).
Proof (z n ) 1 n N is observed because the designs are nested. The loss function can be decoupled because

(Y 1 , Z 2 , • • • , Z N ) are independent.
If all the designs are nested, the parameters η n can be estimated separately by optimizing the loss functions l n .

Non nested designs

If the designs are not nested, the loss function is not decoupled. An EM (Expectation-Maximization, see [START_REF] Friedman | The elements of statistical learning[END_REF]) algorithm is then used to reduce the dimension of the optimization problem. To define this algorithm, the notion of complete data is introduced. In the case of not nested designs

Definition 1 (Complete data) Let ∼ X n (∀1 n N ) denote the union of all training samples that concern Z n . ∼ X n is a matrix of d 1 + • • • + d n columns, composed of the concatenation of the parts of X n , • • • , X N corre- sponding to the input variables x I1∪•••∪In . The complete data is the random vector: (Y 1 ( ∼ X 1 ), Z 2 ( ∼ X 2 ), • • • , Z N ( ∼ X N )).
∼ X 1 , • • • , ∼ X N play the role of X 1 , • • • , X N in
l c (η; z 1 , • • • , z n ) = N n=1 l n c (η n ; z n ),
where l 1 c (η 1 ; z 1 ) is the loss function associated to the training data Y 1 ( ∼ X 1 ) = z 1 , and

l n c (η n ; z n ) (n ∈ 2, N ) is the loss function associated to the training data Z n ( ∼ X n ) = z n .
As some of the complete data is not observed (missing data), an EM algorithm seems adequate to optimize the loss function.

Definition 2 (EM algorithm) The EM algorithm is defined by • Expectation Instead of the observed data loss function, the expectation of the complete data loss function conditioned by the observed data is considered. This quantity can be decomposed in N terms:

Q(η, η * ) = E η * [l c (η; T ) | Y = y] , = N n=1 Q n (η n ; η * ), with ∀n ∈ 1, N Q n (η n , η * ) = E η * l n c (η n ; Z n ( ∼ X n )) | Y = y .
• Maximization The EM algorithm consists in building the sequence (η (i) ) i 0 = (η

(i) 1 , • • • , η (i)
N ) i 0 such that η (0) is user defined and ∀i 0, η (i+1) is solution of the optimization problem

min η Q(η, η (i) ). For all n ∈ 1, N , η (i+1) n
is solution of the following optimization problem:

min ηn Q n (η n , η (i) ). (7)
See appendix C for the explicit formulas of the Q n . It is known that the sequence (η (i) ) i 0 verifies that (l(η (i) ; y)) i 0 is a decreasing sequence (see appendix D). Finally, η is estimated by a term of the sequence of sufficiently high rank, the choice of the rank being done by imposing a maximum number of iterations and a threshold on the variation of likelihood between iterations.

Remark To optimize

Q n , the matrices cov Z n ( ∼ X n ), Z n ( ∼ X n ) = σ 2 n ρ θn ( ∼ X n , ∼ X n ) (n ∈ 1, N ) need to be inverted. First, for all n ∈ 1, N , ∼
X n is assumed not having any redundant points such that the matrices have no identical rows. This condition is sufficient for the first matrix to be invertible as long as ρ θ1 is a positive-definite kernel. As it is not the case of the ρ θn (n 2), to ensure that the matrices are not singular, each ∼ X n (n ∈ 2, N ) is also supposed to not contain any point on which Z n is null so there are no rows full of zeros in the matrices.

Examples of application

In what follows, three metamodels are compared on three examples, two analytic functions and one output from an industrial code:

• seqGPR: the metamodel introduced in this paper. The kernel of Y 1 and the ones on which are based (Z n ) N n=2 are stationary tensor product matern 5 2 covariance kernels. • K_N: a classic kriging metamodel with a stationary tensor product matern 5 2 covariance kernel and a constant mean which is built on the last training sample (X N , y N ).

• K_tot: a classic kriging metamodel similar to K_N, but built on all training samples

(X 1 , y 1 ), • • • , (X N , y N ).
The philosophy of the method seqGPR is to explain a maximum of the data by Y 1 , and to correct it thanks to the (Z n ) N n=2 , so all the parameters for Y 1 are estimated individually and on the contrary some parameters of the Z n are grouped. Different parameter sparsities for the Z n are tried, to seek balance between fitting and robustness:

• Full: θ n = (θ 1 n , • • • , θ d1+•••+dn n ), one component θ i n for each input variable x i . • Robust: θ n = (θ I1∪•••∪In-1 n , • • • , θ I1∪•••∪In-1 n I1∪•••∪In-1 , θ 1 n , • • • , θ dn n In ), with θ I1∪•••∪In-1 n
common to all components of x I1∪•••∪In-1 , and one parameter θ i n for each component of x In .

Four versions of seqGPR are compared. P_full and P_rob (resp. Red_full and Red_rob) use P (resp. Red) processes defined in subsection 3.2 (resp. 3.1), respectively with parameter sparsity of type Full and Robust for the (Z n ) N n=2 . The metamodels are compared in terms of RMSE on a test sample. The training samples are not nested so the EM algorithm is used to estimate the parameters of the 4 versions of seqGPR. A simple likelihood estimation is used for the two kriging models.

Analytic example in dimension 4

The output f considered in this example is a function of 4 inputs

x = (x 1 , x 2 , x 3 , x 4 ): f (x 1 , x 2 , x 3 , x 4 ) = f 1 (x 1 , x 2 ) + f 2 (x 1 , x 2 , x 3 , x 4 ), with:        f 1 (x 1 , x 2 ) = 4 -2.1(4x 1 -2) 2 + (4x1-2) 4 3 (4x 1 -2) 2 + (4x 1 -2)(2x 2 -1) + -4 + 4(2x 2 -1) 2 (2x 2 -1) 2 , f 2 (x 1 , x 2 , x 3 , x 4 ) = 4 exp -x -0.3 2 .
The study is composed of two steps:

• At step 1, computer code evalutations are run at DoE X 1 in dimension 2, such that f 1 (X 1 ) = y 1 , with

f 1 defined by f 1 (x 1 , x 2 ) = f (x 1 , x 2 , x3 , x4 ).
Only the first two variables x 1 and x 2 are free and the other two variables are fixed:

x3 = x1+x2 2 , x4 = 0.2x 1 + 0.7.
• At step 2, new simulations are launched at points X 2 , a design in dimension 4. The last two variables (x 3 , x 4 ) are now released.

The total number of variables is d = 4, the number of steps is N = 2, the index set of variables released at step 1 is I 1 = {1, 2} and the index set of variables released at step 2 is I 2 = {3, 4}.

Figure 5 shows RMSE of the different methods computed on a Sobol sequence of size 10000 used as a test set. The RMSE is computed on 100 different training samples (X 1 , y 1 ) and (X 2 , y 2 ). The 100 RMSE's are represented by a boxplot for each metamodel. Different sizes of training sample are tested. The standard deviation of the output on the test set is represented by the black horizontal line.

All metamodel performances improve with the size of the training sample. Results show that including the previous information of (X 1 , y 1 ) is useful, as it greatly improves the performance of the kriging metamodel (K_tot is better than K_2). The robust seqGPR metamodels are better than the full, they are equivalent to or better than K_tot. It seems that the cases with a small or high number of training points are more discriminating than in the middle cases. With a small number of points (10pts-5pts), K_tot seems more destabilized than seqGPR. With a high number of points (20pts-40pts), seqGPR is more accurate than K_tot. Finally, seqGPR using Red are better than seqGPR using P. In the following examples, only the robust versions of seqGPR, which are more performing and whose parameter estimation is less complex, are compared to K_tot.

Analytic example in dimension 15

The following example should be less favorable to seqGPR. The objective function considered is in higher dimension and not decomposed as a sum of functions that respects the order in which the variables are released

f : [-3, 3] 15 → R x → a 1 x + a 2 sin x + a 3 cos x + x M x. (8)
The function f was first used in [Oakley and O'Hagan, 2004]. The values of its coefficients a 1 , a 2 , a 3 and M can be found in www.sheffield.ac.uk/st1jeo. The inputs are rescaled in [0, 1] and are rearranged by decreasing order of Sobol index.

Figure 5: Boxplots of the metamodel RMSE's on 100 studies, for the analytic example in 4D. In abscissa are the sizes of the training samples.The following sizes are tested from left to right: X 1 of size 10 and X 2 of size 5 (10pts-5pts), X 1 of size 10 and X 2 of size 10 (10pts-10pts), X 1 of size 10 and X 2 of size 20 (10pts-20pts), X 1 of size 20 and X 2 of size 10 (20pts-10pts), X 1 of size 20 and X 2 of size 20 (20pts-20pts), X 1 of size 20 and X 2 of size 40 (20pts-40pts) K_tot P_rob Red_rob Median (•10 -4 ) 46.039 37.135 36.903 Interquartile range (•10 -4 ) 7.258 6.660 6.570

Table 1: Performance of the metamodels for the 100 studies made on the analytic 15D example. The performances are shown in terms of median of RMSE's and interquartile range (q 75% -q 25% ).

(a) Car engine cooling system (b) Fan system In practice, the first studies are done on the most influential inputs, chosen according to experts' physical knowledge. Here, as the function is analytical, the choice of the releasing order for the inputs is made by sensitivity analysis, which is directly performed on the function. Following the distribution of the Sobol indices (see figure 6), a study of N = 3 steps is made, adding respectively the groups of variables I 1 = {1, 2, 3}, I 2 = {4, 5, 6, 7, 8, 9}, and I 3 = {10, 11, 12, 13, 14, 15}. The variables that are fixed are set to 0.5. Three training samples, one for each step, are generated with 5 points per dimension considered at the step: The performances of the metamodels are shown in table 1. All versions of seqGPR are better than the classic kriging, with a slight improvement for the Red version in comparison with the P version. Red and P are almost equivalent as opposed to the previous example. This can be due to the fact that the P version disturbs less the initial stationary kernel because the variables which are fixed at a given step are set to a constant instead of varying in function of the free variables, as it is the case in the previous example.

X 1 ⊂ [0, 1] 3 of size 15, X 2 ⊂ [0, 1]

Industrial example: fan system in dimension 15

The industrial product which is studied in this section is the fan system which is part of a car engine cooling system from Valeo company (see figure 7). This cooling system is composed of a cold fluid circulating in part in the car engine to regulate its temperature and in part in a radiator where it is itself cooled down. When the car moves, the wind generated by the car speed and reaching the radiator is sufficient to evacuate the heat from the fluid. When the car is motionless, the fan is activated to replace the wind.

The output considered here is the Pressure difference (∆P ) between the upstream and the downstream of the air flow crossing the fan. It is function of 15 input variables. One input is the flowrate (Q). The 14 others are geometric. The fan blade geometries are supposed identical with each other. The geometry is monitored at 5 different sections (see figure 8a). At each section, the stagger angle and the chord length are controlled (see figure 8b). The chord is the line formed by the two borders of the blade at the considered section. The stagger angle is the angle between the chord and the rotation axis. The sweep is manipulated at each section. It is defined as the distance between the line formed by the rotation axis and the right border of the first section and the right border of the considered section (see figure 8c). At section 1, it is always equal to 0 and therefore not kept as an input.

To summarize, the variables are Table 2: Performance of the metamodels for the 100 studies made on the industrial 15D example. The performances are shown in terms of median of RMSE's and interquartile range (q 75% -q 25% )

• The chord lengths at the five sections: Chord1, Chord2, Chord3, Chord4, Chord5

• The sweeps at sections 2 to 5: Swe2, Swe3, Swe4, Swe5

In the rest of this example, all variables are adimensionalized in [0, 1] and the output ∆P is adimensionalized in [-1, 1].

Industrial experts at Valeo have carried out a two-step (N = 2) study in the context of this work.

• At step 1, all the sweeps are fixed to constants: Swe2 = 0.517645, Swe3 = 0.82, Swe4 = 1, Swe5 = 1.

An OLH of size 126 has been created on the 11 free inputs (I 1 = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}).

• At step 2, all the inputs are free (I 2 = {12, 13, 14, 15}). An OLH of size 299 has been created for the 15 inputs.

In order to obtain robust results, a cross-validation procedure is implemented. Each OLH is decomposed in 30 subsamples of size 50 optimized for the maximin criterion (they can share common points). The 30 complementary subsamples of size 249 from the OLH of step 2 are used as test samples.

Table 2 shows performances of the different metamodels on the test samples. The metamodel seqGPR is better than K_tot. Its version with P is this time slightly better than with Red, again probably because the fixed inputs are constant as in the previous example.

Conclusion

In the framework of Gaussian process regression, the problem of building metamodels in nested spaces of increasing dimension is studied. The variables are first fixed and then released progressively. The suggested method, called seqGPR, consists in creating a kriging metamodel with a non-stationary kernel on the reunion of all training samples. This kernel is the covariance kernel of a process defined recursively. At a given step, the process modeling the output is equal to the sum of the process modeling the output at the previous step and an independent correction term.

The correction term must be a Gaussian process which is null on a part of the input space. Two candidates are proposed: the Red (Reduced) process, and the P (Preconditioned) process. The P process, found in the literature, is intractable. Two ways are tried to make it tractable. The first way is to discretize its spectral decomposition. But this technique leads to a Gaussian process conditioned to be null on some points of the space. This process is never exactly null everywhere in the part we want. To approach this nullity, a lot of points are needed in the conditioning, and that is very computationally expensive. In a second approach, a tractable exact expression is directly sought, in exchange for a modification of the latent process it is built upon.

Then, the issue of the parameter estimation is addressed. Instead of directly optimizing the likelihood, an EM algorithm is implemented, which is adapted to the additive structure composed of independent processes. Finally, the metamodel seqGPR is compared to a classic kriging metamodel, based on a stationary kernel, on two analytic and one industrial examples. Different levels of parameter sparsity for the correction processes are tried. The level of sparsity that seems to always be better than classic kriging is the Robust one with either Red or P processes. The Red version seems to achieve higher performances than the P version most of the time, so the Red_rob version is recommended.

Software and Acknowledgments

Implementation have been done in R, where all the kernels and methods of estimation and prediction have been implemented using the package RcppArmadillo.

This work, carried out at Ecole Centrale de Lyon, was financially subsidized by Valeo. The authors wish to show their particular gratitude to Mr. Manuel Henner and Mr. Nicolas Yoan François for their support and fruitful discussion on the industrial point of view. The authors are grateful to the Members of GdRMascotnum for their useful advice, and in particular to Mr. Olivier Roustant and Mr. Xavier Bay for constructive discussions on the subject on the theoretical point of view. The authors would also like to acknowledge Mr. Philippe du Bouays and Ms. Laura Gay for their careful proofreading.

Appendix A Proof of proposition 1

This section deals with the proof of proposition 1, which gives an approximation of the process Z P (defined in equation ( 5)) by discretizing its spectral decomposition.

Proof

• Approximation of (λ n , ∼ φ n ):

The following eigenvalue problem is considered:

[0,1] d J σ 2 r((x J , g(x J )), (s J , g(s J ))) ∼ φ(s J )ds J = λ ∼ φ(x J ) ∀x J ∈ [0, 1] d J .
It is discretized using a Monte-Carlo method. A sample (s

(i) J ) 1 i L is generated uniformly in [0, 1] d J to build D = s (i) J , g(s (i) J ) 1 i L which is a discretization of the subspace {(s J , g(s J )), s J ∈ [0, 1] d J }. The Monte-Carlo approximation of the integral is: 1 L L i=1 σ 2 r((x J , g(x J )), (s (i) J , g(s (i) J ))) ∼ φ(s (i) J ).
Discretizing in x J with D, the eigen problem becomes a finite dimensional one:

1 L σ 2 r(D, D) ∼ Φ = γ ∼ Φ. (9)
The solutions of (9) are noted (γ n , V n ) 1 n L . The V n are taken such that they verify:

V n V m = δ nm . The discretization of ∼ φ n (noted ∼ Φ n ) must verify the discrete equivalent of [0,1] d J ∼ φ n (t J ) ∼ φ m (t J )dt J = δ nm , which is 1 L ∼ Φ n ∼ Φ m = δ nm .
So the following relation between V n and

∼ Φ n is verified: V n = 1 √ L ∼ Φ n ⇔ ∼ Φ n = √ LV n .
• Approximation of φ n :

As

φ n (x J , x I ) = 1 λn [0,1] d J σ 2 r ((x J , x I ), (s J , g(s J ))) ∼ φ n (s J )ds J ∀(x J , x I ) ∈ [0, 1] d J × [0, 1] d I ,
using the same Monte-Carlo approximation and replacing (λ n , ∼ φ n ) by (γ n , ∼ Φ n ), the following approximation of φ n is obtained:

φ D n (x) = 1 γn 1 L σ 2 r(x, D) ∼ Φ n ∀x ∈ [0, 1] d J +d I .
• Decomposition of the matrix σ 2 r(D, D)

-1

(γ n , V n ) n 1 are the eigenvalues and (orthonormal) eigenvectors of 1 L σ 2 r(D, D). So ( 1 γn , V n ) n 1 are the eigenvalues and (orthonormal) eigenvectors of L σ 2 r(D, D)

-1 , and:

L σ 2 r(D, D) -1 = L n=1 1 γ n V n V n , = L n=1 1 γ n 1 √ L ∼ Φ n 1 √ L ∼ Φ n , = 1 L L n=1 γ n ∼ Φ n ∼ Φ n . So σ 2 r(D, D) -1 = 1 L 2 L n=1 γ n ∼ Φ n ∼ Φ n .
• Approximation of the process Z P :

The integral [0,1] d J ∼ φ n (t J ) ∼ Z(t J , g(t J ))dt J , involved in the formula of Z P (see (5)), is discretized the same way as before. It becomes:

1

L ∼ Φ n ∼ Z(D).
The approximation of the process is:

Z D (x) = ∼ Z(x) - L n=1 φ D n (x) 1 L ∼ Φ n ∼ Z(D) , = ∼ Z(x) - L n=1 1 L 1 γ n σ 2 r(x, D) ∼ Φ n 1 L ∼ Φ n ∼ Z (

B Proof of proposition 2

This section deals with the proof of proposition 2 which gives a closed form formula of the process Z P (see equation ( 5)) for a particular choice of the kernel of ∼ Z.

Proof

• Eigenvalue problem:

The eigenvalue problem can be rewritten as:

[0,1] d J σ 2 r ((x J , g(x J )), (s J , g(s J )))

∼ φ n (s J )ds J = λ n ∼ φ n (x J ), ∀x J ∈ [0, 1] d J , ⇔ [0,1] d J σ 2 r J (x J , s J ) r I (g(x J ) -g(x J ), g(s J ) -g(s J ))

∼ φ n (s J )ds J = λ n ∼ φ n (x J ), ⇔ [0,1] d J σ 2 r J (x J , s J ) r I (0, 0)

=1 ∼ φ n (s J )ds J = λ n ∼ φ n (x J ), ⇔ [0,1] d J σ 2 r J (x J , s J ) ∼ φ n (s J )ds J = λ n ∼ φ n (x J ). (10)
• Expression of φ n :

φ n can be rewritten as:

φ n (x J , x I ) = 1 λn [0,1] d J σ 2 r ((x J , x I ), (s J , g(s J )))

∼ φ n (s J )ds J , ∀(x J , x I ) ∈ [0, 1] d J × [0, 1] d I , = 1 λn [0,1] d J σ 2 r J (x J , s J ) r I (x I -g(x J ), g(s J ) -g(s J ))

∼ φ n (s J )ds J , = 1 λn [0,1] d J σ 2 r J (x J , s J ) ∼ φ n (s J )ds J r I (x I -g(x J ), 0)) ,

= 1 λn λ n ∼ φ n (x J ) r I (x I -g(x J ), 0)) , = ∼ φ n (x J )r I (x I -g(x J ), 0)) .
The second last equality is due to the fact that ∼ φ n is solution of the eigenvalue problem (10).

• The P process can be rewritten as: The second equality is the formula of the conditional expectation (see (6)).

So

Z P (x) = ∼ Z(x) -r I (x I -g(x J ), 0) ∼ Z(x J , g(x J )).

Figure 1 :

 1 Figure 1: Illustration of the model for N = 2. f is the realization of Y 2 (see panels on the left and center). On the cross section of equation x I2 = xI2 , Y 2 = Y 1 (see panels on the left and right)

Figure 2 :

 2 Figure 2: Illustration of the model for N = 2 (see figure 1), when xI2 is constant.

Figure 3 :

 3 Figure3: On middle panel, a realization of the latent process from which are built the two candidates for Z 1 that should be null on the line x 2 = 0. On left panel, the process which is removed from the latent process to create the Red process. On right panel, the process (denoted by E for expectation) which is removed from the latent process to create the P process.

Figure 4 :

 4 Figure 4: Definition of the

Figure 4

 4 Figure 4 gives an illustration of the different

  the case of nested designs. According to proposition 3, if the values of Y 1 on ∼ X 1 , denoted by z 1 , and the values of the Z n on the ∼ X n (n ∈ 2, N ), denoted by z n , were observed, the loss function associated to all those training samples, denoted by l c , could be decomposed as

Figure 6 :

 6 Figure 6: Sobol indices of the objective function (8).

Figure 7 :

 7 Figure 7: On left panel, diagram of the car engine cooling system. On right panel, diagram of fan system. [Valeo, ]

  9 of size 45, and X 3 ⊂ [0, 1] 15 of size 75. A Sobol sequence of size 10000 in dimension 15 is used as a test sample. The study is done 100 times, each time with new training samples.

Figure 8 :

 8 Figure 8: On the left panel, visualization of the five sections of the blade. On the middle panel, definition of the stagger angle and the chord length at a given section. On the right panel, definition of the sweep at section 5. K_tot P_rob Red_rob Median (•10 -3 ) 44.585 41.207 42.088 Interquartile range (•10 -3 ) 5.892 4.657 4.609

  ) = 0 .The process approximating Z P is the conditioned (on a finite set of points) Gaussian process Z D . It is a centered Gaussian process of covariance kernel σ 2 ρ D with:ρ D (x, x ) = r(x, x ) -r(x, D)r(D, D) -1 r(D, x ) ∀x, x ∈ [0, 1] d J +d I .

Z

  P (x J , x I ) = ∼ Z(x J , x I ) -+∞ n=1 φ n (x J , x I ) ∼ φ n (s J ) ∼ Z(s J , g(s J ))ds J , ∀(x J , x I ) ∈ [0, 1] d J +d I , = ∼ Z(x J , x I ) -+∞ n=1 ∼ φ n (x J )r I (x I -g(x J ), 0)) ∼ φ n (s J ) ∼ Z(s J , g(s J ))ds J , = ∼ Z(x J , x I ) -+∞ n=1 ∼ φ n (x J ) ∼ φ n (s J ) ∼ Z(s J , g(s J ))ds J r I (x I -g(x J ), 0)) ,and, because ∼ Z(x J , g(x J )) belongs to the sub Gaussian space engendered by the{ ∼ Z(s J , g(s J )) ∀s J ∈ [0, 1] d J }, its projection in this subspace is equal to itself: ∼ Z(x J , g(x J )) = E[ ∼ Z(x J , g(x J )) | ∼ Z(s J , g(s J ))∀s J ], ∀x J ∈ [0, 1] d J , x J ) r I (g(x J ) -g(x J ), 0)) J , g(t J ))ds J , = +∞ n=1 ∼ φ n (x J ) ∼ φ n (s J ) ∼Z(s J , g(s J ))ds J .

•

  At step 1, a first series of computer code evaluations is run. The set x I1 of the first d 1 components of x, is free in [0, 1] d1 . Different values are explored. They are stored in DoE X 1 ⊂ [0, 1] d1 . The other components are fixed to preset values xI2∪•••∪I N (entirely determined by x I1 ). Let f 1 denote the corresponding restriction of f on the subspace [0, 1] d1 . The function f 1 is assumed to be the realization of a stationary Gaussian process Y 1

  variables x I2 , fixed at step 1, are released. Different values of x I1∪I2 , stored in DoE X 2 , are explored. The other variables x I3∪•••∪I N are fixed to the new set of values xI3∪•••∪I N ∈ [0, 1]

  d3+

  ••+d N . Let f 2 be the corresponding restriction of f on the subspace [0, 1] d1+d2 . The function f 2 is assumed to be the realization of a Gaussian process Y 2

  d1+

  ••+dn . The variables x I1∪•••∪In-1 and x In are free and the others x In+1∪•••∪I N are set to xIn+1∪•••∪I N ∈ [0, 1]

  dn+1+

  ••+d N . The corresponding restriction of f on this subspace, f n , is evaluated on X n . The values are stored in y n = f n (X n ). The function f n is modeled as the realization of a Gaussian process Y

n : Ω × [0, 1]

d1+

  ••+dn → R. Following the same arguments as in step 2, Y n is defined as:

equal to twice the negative log-likelihood up to a constant

C Formulas of the EM procedure

This section gives the formulas of the Q n involved in the optimization problems (7) used to estimate the parameters η n in the EM algorithm. 1∼ X1 is the unit column vector of size nrows ∼ X 1 . Denoting by X 1 and X 2 two DoE's, 0 X 1 ,X 2 is the null matrix of size nrows(X 1 ) × nrows(X 2 ).

Using the expectation of a quadratic form formula to Gaussian vectors:

Partial analytical solution The optima m (i+1) and σ (i+1) n

(n ∈ 1, N ) have analytical forms obtained by solving the system formed when the corresponding partial derivatives of the Q n vanish. Finally, at each new iteration i + 1, the goal is to find θ

, where E n (η (i) ) and C n (η (i) ) (n ∈ 1, N ) are given in equation ( 11).

D Proof of the monotonicity of the sequence from the EM algorithm

. Let y, u, and t be the realizations of respectively Y, U, and T. By definition of the conditional density:

So (up to a constant) l(η; y) c = l c (η; t) + 2 log(h∼

Replacing the observations by the corresponding random variables

Taking the expectation assuming η and conditioning by Y = y l(η; y)

The difference between the loss function taken at the updated and previous terms of the EM sequence is the sum of two quantities.

The first quantity is negative by definition of the EM algorithm. The second quantity can be rewritten as: Indeed the inequality l(η (i+1) )-l(η (i) ) 0 is verified by the sequence η (i) i built following the EM algorithm.