Martina Flöser

Serge Florens

Thierry Champel

Diagrammatic Approach for the High-Temperature Regime of Quantum Hall Transitions

We use a general diagrammatic formalism based on a local conductivity approach to compute electronic transport in continuous media with long-range disorder, in the absence of quantum interference effects. The method allows us then to investigate the interplay of dissipative processes and random drifting of electronic trajectories in the high-temperature regime of quantum Hall transitions. We obtain that the longitudinal conductance σxx scales with an exponent κ = 0.767 ± 0.002 in agreement with the value κ = 10/13 conjectured from analogies to classical percolation. We also derive a microscopic expression for the temperature-dependent peak value of σxx, useful to extract κ from experiments.

Introduction.-The geometric concept of percolation is ubiquitous to electronic transport in strongly disordered media [1], in both the classical and quantum realm. Indeed, building on earlier studies in the context of metallic alloys and granular materials [2], recent advances have extended percolation ideas to the description of quantum phases in low-dimensional electron gases, ranging from metal/insulator transitions at low magnetic field to the high magnetic field regime associated to the quantum Hall effect [3][4][5]. Despite this very seductive geometrical analogy, difficulties arise for a microscopic description of transport because the electrical current does not just propagate on simple geometrical objects, such as the bulk or the boundaries of a percolation network. In fact, in a dissipative system the current density always spreads along extended structures, so that fractality of the transport network may be smeared in realistic situations [6]. While fully numerical simulations of transport models can account for such complexity [3,5], they bring finite size effects that give limitations for quantitative description of transport. For instance, an important question for metrological purposes [7] is the precise understanding of the accuracy of Hall conductance quantization, where percolation is known to play a role, both from theoretical grounds [6,8,9] and from local density of states [10] and transport measurements [11][12][13].

Our goal in this Letter is to show that percolation features of transport in continuous disordered media can be captured analytically by a diagrammatic approach, starting from local Ohm's law:

j(r) = σ(r)E(r), (1)
with j the local current density and E the local electrical field. This introduces σ(r) the local conductivity tensor, a spatially-dependent quantity due to inhomogeneities, that naturally encodes altogether dissipation, disorder and confinement [6,14,15]. The local conductivity model is expected to be accurate at high enough temperatures whenever phase-breaking processes, such as electron-phonon scattering, occur on length scales that are shorter than the typical variations of disorder. However, quantum mechanics may still be important to determine microscopically the quantitative behavior of the local conductivity tensor [16,17]. The main difficulty thus lies in solving the continuity equation ∇ • j = 0 in the presence of long-range random inhomogeneities in the sample, see Fig. 1. General formalism.-Our starting point follows early ideas proposed by several authors [18,19], where effective conductivity approaches were developed based on a local conductivity tensor σ(r). We consider the general situation of an arbitrary and continuous distribution of conductivity in a macroscopic d-dimensional sample of volume V , bounded by a surface S. The experimentally accessible quantity is the average current density j = σeff E 0 which is driven by applying a constant electric field E 0 at the boundary of the sample. This defines a position-independent effective conductivity tensor σeff , which is nothing but the macroscopic conductance tensor, up to a geometrical prefactor. Following Ref. 19, we decompose (arbitrarily at this stage) σ(r) = σ0 + δσ(r) into uniform and fluctuating parts, re-spectively. By expressing the electrical field by its scalar potential E(r) = -∇Φ(r), the continuity equation leads to the boundary value problem:

∇ • [σ 0 ∇Φ(r)] = -∇ • [δσ(r)∇Φ(r)] in V (2) Φ(r) ≡ Φ 0 (r) = -E 0 • r on S. (3)
By introducing the Green's function G(r, r ′ ) defined by

∇ • [σ 0 ∇G(r, r ′ )] = -δ(r -r ′ ) in V (4) G(r, r ′ ) = 0 for r on S, (5)
the scalar potential is formally given by

Φ(r) = Φ 0 (r) + V d d r ′ G(r, r ′ )∇ ′ • [δσ(r ′ )∇ ′ Φ(r ′ )] (6)
with the short-hand notation ∇ ′ = ∇ r ′ . Integrating by parts with ∇ ′ G(r, r ′ ) = -∇G(r, r ′ ) and taking the gradient on both sides of Eq. ( 6) leads to

E(r) = E 0 + V d d r ′ ∇ • [∇G(r, r ′ )δσ(r ′ )E(r ′ )] (7)
= E 0 + V d d r ′ Ĝ0 (r, r ′ )δσ(r ′ )E(r ′ ), (8)
where Ĝ0

ij = ∂ ∂ri ∂ ∂rj G(r, r ′ ).
Finally, multiplying Eq. ( 8) by δσ(r) and introducing a new local tensor χ such that δσ(r)E(r) = χ(r)E 0 , we obtain:

χ(r)E 0 = δσ(r)E 0 + δσ(r) V d d r ′ Ĝ0 (r, r ′ ) χ(r ′ )E 0 . ( 9
)
As Eq. ( 9) is valid for all possible choices of E 0 , the following tensorial equation also holds:

χ(r) = δσ(r) + δσ(r) V d d r ′ Ĝ0 (r, r ′ ) χ(r ′ ). ( 10
)
Spatial averaging of the current j(r) = [σ 0 + χ(r)]E 0 over conductivity fluctuations δσ(r) leads therefore to the effective conductivity σeff = σ0 + χ , where the spatial average on χ is performed while enforcing the integral equation (10). Although sample boundaries could be considered in principle, we now focus on an infinite sample, so that the Green's function [Eq. ( 4)] becomes translation-invariant

G(r, r ′ ) = d d p (2π) d e ip•(r-r ′ ) pσ 0 p + 0 + , (11)
where 0 + is a small positive quantity which ensures the correct boundary condition at infinity [Eq. (5)]. Systematic expansion at strong-dissipation.-Previous works either considered a mean-field solution of Eq. (10) in the peculiar case of binary randomness in the local conductivity tensor [19], or computed low order contributions for continuous disorder distribution [18,20]. Our aim is to present a systematic expansion controlled by weak fluctuations of the conductivity and to show that the nonperturbative regime of large conductivity fluctuations can be tackled by sufficient knowledge of the perturbative series. The spatial average on χ(r) can be obtained clearly after iterating Eq. ( 10) to all orders:

χ(r) = δσ(r) + d d r 1 δσ(r) Ĝ0 (r, r 1 )δσ(r 1 ) (12)
+ d d r 1 d d r 2 δσ(r) Ĝ0 (r, r 1 )δσ(r 1 ) Ĝ0 (r 1 , r 2 )δσ(r 2 ) + ... which can be expressed graphically as in Fig. 2. For in-

χ(r) = r + r r 1 + r r 1 r 2 + ... Ĝ0 (r -r 1 ) = r r 1 δσ(r) = r FIG.
2. Graphical representation of the strong-dissipation expansion (12) of the self-consistent transport equation (10).

coherent transport, self-averaging occurs and the spatial average over the local conductivity fluctuations may be replaced by an ensemble average. Let us first illustrate the method for a purely resistive and isotropic medium, so that σ0 = σ 0 1 and δσ(r) = δσ(r) 1, with δσ(r) = 0. In the limit of strongdissipation compared to the typical fluctuations of conductivity [σ 0 ≫ δσ 2 ], we get σeff = σ xx 1 with

σ xx = σ 0 - 1 σ 0 d d r d d p (2π) d p 2
x e ipr p 2 + 0 + δσ(r)δσ(0) (13)

= σ 0 - 1 σ 0 d d r δ(r) d δσ(r)δσ(0) = σ 0 - δσ 2 dσ 0 .
We thus recover previous results [20] obtained for weakly disordered media, which predict a reduction of the macroscopic conductance due to randomly distributed resistive barriers. Clearly, nontrivial geometrical aspects are absent at this order, because the dominant background of conductivity σ 0 prevents the percolating network to establish. This general formulation of transport [Eq. (12)] is immediately appealing because arbitrary orders of the strong-dissipation expansion can be generated in a compact fashion, fostering hope that the difficult limit of large conductivity fluctuations can be tackled by standard resummation methods.

Simplification for Gaussian randomness.-Under some microscopic assumptions, the conductivity tensor may follow a random Gaussian distribution, according to δσ(r) = 0 and δσ ij (r)δσ kl (r ′ ) = C ij;kl (rr ′ ), so that all moments of the local conductivity tensor are determined from Wick's theorem (in particular, all odd correlations vanish here). This hypothesis leads to a familiar-looking diagrammatic formulation for the strong-dissipation expansion, as shown in Fig. 3. An important technical point is that all particle reducible graphs (diagrams that can be split in two parts by cutting a single line of Ĝ0 ) are identically zero. This is because all such contributions contain the zero momentum limit of the Green's function [ Ĝ0 ] ij (p) = -p i p j G(p) which vanishes at zero momentum according to Eq. ( 11) (note the crucial role of the regularization parameter). Interestingly, the conductance correction χ now takes the precise form of a self-energy, in contrast to a fully quantum formulation of electronic transport [START_REF] Akkermans | Mesoscopic Physics of Electrons and Photons[END_REF] where vertex corrections associated to interference effects need to be accounted for. In what follows, we wish to use the method with the challenging regime of a strongly fluctuating local conductivity, that may lead to geometrical effects related to classical percolation. Clearly, the general perturbation series (12) in powers of δσ 2 /σ 2 0 then breaks down, so that high order terms will be needed.

χ(r) = r r 1 + + r r 1 r 2 r 3 r r 1 r 2 r 3 +...
Percolation regime of the semiclassical Hall effect.-We henceforth consider the semiclassical regime of the quantum Hall effect, which occurs in very high mobility twodimensional electron gases at large perpendicular magnetic field [12,13]. General physical arguments [6] as well as microscopic calculations [16,17] show that the electron dynamics can be described in this regime by a local Ohm's law with a randomly fluctuating Hall conductivity σ H (r) = σ H + δσ(r):

σ(r) = σ 0 -σ H -δσ(r) σ H + δσ(r) σ 0 with δσ(r) = 0. ( 14
)
According to the classical Hall's law, such purely offdiagonal fluctuations of the conductivity correspond to spatial modulations of the electron density brought by long-range random impurities [6,14]. The diagonal part in Eq. ( 14) accounts phenomenologically for dissipative processes, such as electron-phonon scattering, and is supposed for simplicity to be spatially uniform.

The explicit connection to geometrical percolation can now be made. At vanishing dissipation σ 0 → 0, drift currents follow from Hall's law and propagate along constant lines of Hall conductivity. Indeed, from Maxwell's equation ∇ × E = 0 and current conservation ∇ • j = 0, one gets the transport equation [∇δσ(r)] • j = 0. The lines of constant δσ(r) are typically closed, so that all electronic states are localized, except the ones living on the percolation cluster. However, the percolating state does not contribute to macroscopic transport either, as it must necessarily pass through saddle-points of the disordered landscape, where the transport equation becomes undetermined. Thus having finite σ 0 is required to establish a finite conductance in the sample, by connecting the different nearly localized states. This difficulty has led authors [6] to wonder whether purely geometric arguments are sufficient to understand the transport properties at small but finite dissipation, because the current carrying states become broad filaments that may smear the fractal structure of the percolation cluster. This question is now investigated in a controlled fashion.

At high temperature, the Hall conductivity fluctuations given by Eq. ( 14) follow the Gaussian distribution of disorder [START_REF][END_REF]. We also consider for simplicity Gaussian spatial correlations δσ(r)δσ(r ′ ) = δσ 2 e -|r-r ′ | 2 /ξ 2 , with correlation length ξ. Inspection of the diagrammatic series depicted in Fig. 3 shows that the effective conductivity obeys the following expansion:

σeff = 0 -σ H σ H 0 + σ 0 1 + ∞ n=1 a n δσ 2 n σ 2n 0 1 0 0 1 (15)
with dimensionless coefficients a n collecting all diagrams of order n in perturbation theory in δσ 2 /σ 2 0 . The Hall component is therefore not affected here, while the longitudinal conductance receives nontrivial corrections that encode the interplay of dissipation and percolation. The diagrammatic formulation of transport allowed us to compute this series up to sixth order [START_REF][END_REF].

As understood previously, the effective longitudinal conductivity σ xx must vanish when σ 0 → 0 for a continuous local conductivity model, and previous works [6,8,9] suggested a power-law dependence σ xx ∼ C δσ 2 κ/2 σ 1-κ 0 at small σ 0 , with nonuniversal dimensionless constant C and universal critical exponent κ characterizing the transport properties. While κ = 10/13 is often quoted as an exact value [1,6,8,9], Simon and Halperin [6] argued that one could not completely rule out the possibility that finite dissipation may spoil the connection to geometrical percolation and change the value of κ. In order to check that this is not the case, we performed careful Padé resummation [START_REF][END_REF] of the perturbative series (15) up to six loops, see Table I [START_REF][END_REF] built from the perturbative series (15). κ = 0.767 ± 0.002 seems to confirm the conjectured value κ = 10/13 ≃ 0.7692 based on the analogy to classical percolation [6,8,9]. We stress the good convergence of the Padé approximants for all values of the dissipation strength σ 0 , see Fig. 4. Note that partial resummation of perturbation theory in previous works [18] failed to recover the critical behavior associated to percolation in the strong coupling regime and this approximation led to an incorrect saturation of σ xx in the limit σ 0 → 0, which would apply only for transport model with discrete conductivity values [START_REF] Dykhne | [END_REF].

Microscopics of σ xx at plateau transitions.-We finally study the temperature behavior of transport in the percolation dominated regime. At high magnetic field, the local Hall conductivity is explicitly related to the Fermi distribution of Landau levels E m = ω c (m+1/2) with integer m, disorder landscape V (r) and chemical potential µ [16,17,[START_REF][END_REF]:

σ H (r) = e 2 h ∞ m=0 n F [E m + V (r) -µ] (16)
neglecting spin effects. We have introduced here the cyclotron energy ω c = |eB|/m * in terms of Planck's constant , electron charge e, applied perpendicular magnetic field B and effective mass m * . At temperatures T ≫ V 2 , the Fermi distribution n F (E) can be linearized, so that the random conductivity distribution (14) becomes Gaussian. Straightforward analysis [START_REF][END_REF] and our low-dissipation formula lead to a simple expression for the peak conductance measured at the transition region between two Landau levels (k B is Boltzmann's constant):

σ peak xx = σ bg (T, B) 1 + ∞ l=1 4π 2 lk B T ω c csch 2π 2 lk B T ω c κ . (17)
This expression plotted in Fig. 5 shows a sharp crossover at temperature k B T ⋆ = ω c /4 from a low-T power-law

σ bg (T, B) = C[σ 0 (T, B)] 1-κ e 2 h V 2 ω c κ . (18)
(17)-( 18), which combine microscopic parameters (such as the width of the disorder distribution) with geometrical effects through the exponent κ, should be useful for detailed analysis of transport measurement in quantum Hall samples.

Clearly, σ peak xx cannot diverge at T → 0 and is in fact expected to level off when reaching conductance values of the order of e 2 /2h [START_REF] Dykhne | [END_REF]24]. In this very-low-temperature regime, the linearization of the local Hall conductivity (16) breaks down, thereby putting a limit to the present diagrammatic calculation. Furthermore, quantum effects become important at low-T and lead [4,13] to a different exponent κ qu. ≃ 3/7 ≃ 0.42. The classical percolation exponent [6,8,9] κ = 10/13 ≃ 0.77 may be observin very high mobility samples dominated by smooth disorder [12,13]. Finally, at temperatures T > ω c /4, the leading magnetic field dependence of the longitudinal conductivity in Eq. ( 17) is provided by the ω -κ c ∝ B -κ term, as discussed previously [11,25].

Conclusion.-We have used a general diagrammatic method to compute fully microscopically the electronic transport in incoherent disordered conductors, leading to accurate determination of critical exponents for the conductivity in the classical percolation regime of the quantum Hall transition. This framework seems also well suited for efficient numerical implementations using the recently developed diagrammatic Monte Carlo methods [26], leading to envision progresses towards more realistic description of quantum Hall transport taking into account disorder effects.

We thank A. Freyn for precious help with symbolic

We consider here the problem of random Gaussian fluctuations of the local Hall conductivity in two dimensions (see Eq. ( 14) in the main text), split into an average Hall component σ H and a fluctuating term δσ(r), defined so that δσ(r) = 0. The dissipationless nature of the Hall component shows up by the fact that σ H exactly drops in the correlation function Ĝ0

ij = ∂ ∂ri ∂ ∂rj G(r, r ′ ): Ĝ0 ij (r) = - 1 σ 0 d 2 p (2π) 2 p i p j e ip•r p 2 + 0 + , (1)
with G(r, r ′ ) defined by Eq. ( 11) in the main text.

The first order diagram contributing to the conductivity is straightforwardly calculated in the case of Gaussian fluctuations of the Hall component in two dimensions (see Eq. ( 14) in the main text):

= d 2 p (2π) 2 K(p)ǫ Ĝ0 (p)ǫ = δσ 2 πσ 0 +∞ 0 dp pe -p 2 2π 0 dθ sin 2 (θ) -cos(θ) sin(θ) -cos(θ) sin(θ) cos 2 (θ) = δσ 2 2σ 0 1 0 0 1 (2)
with K(r) ≡ δσ(r)δσ(0) = δσ 2 e -|r| 2 /ξ 2 , and its Fourier transform K(p) = πξ 2 δσ 2 e -ξ 2 p 2 /4 . Here ǫ denotes the fully antisymmetric 2 × 2 matrix, ǫ = 0 1 -1 0 . Note that the conductivity correction [Eq. ( 2)] is positive and exactly opposite in sign to the one obtained in the case of pure longitudinal fluctuations of the conductivity in Eq. ( 13) of the main text.

All second and third order diagrams can be obtained analytically with the help of symbolic computation, see the results displayed in Table I. The method of computation for the second and third order contributions is to first express each of the several denominators appearing in a given graph using Feynman's identity:

1 x i = ∞ 0 dt i e -tixi . (3)
One can then perform the Gaussian integration over all momenta, and finally compute the remaining integrals over the auxiliary variables t i .

We have not managed to analytically obtain the diagrams of fourth order and beyond (except for the noncrossing ones, see below), and we had therefore recourse to a combination of analytical and numerical steps. First, Σ1 (p) = p p p + q q FIG. 1. Self-energy Σ1(p) entering the calculation of the noncrossing diagrams in Table II.

an automated script was used to generate all possible diagrams, discarding the particle reducible ones, which enables to output explicitely the corresponding functions that require full momentum integration. In order to avoid indefinite integrals, all two-dimensional momenta in an n th order diagram were combined into the hyperspherical coordinate K in dimension 2n, such that K 2 = n i=1 p i 2 . This allows analytical integration over |K|, leaving the bounded integration domain on the hypersphere in 2n dimensions. This numerical step was finally performed using the Vegas Monte Carlo integration routine from the GNU Scientific Library. Because only the complete sum of all diagrams at a given order matters, and since multidimensional integrals are time consuming, we have summed up all the contributions at a given order before performing the integration. The Monte Carlo evaluations were iterated until the relative error was below 0.1%, but we can also ascertain the good convergence of the numerics by benchmarking the routine on analytically tractable diagrams that have no crossings of the propagators, see Table II for comparison. The high (up to 6 th ) order non-crossing diagrams that we considered are obtained in the following way: we remark that these graphs are only composed of bare propagators and of the first order self-energy Σ1 appearing in Fig. 1. The momentum dependence of this self-energy is readily evalu-ated:

Σ1 (p) = δσ 2 σ 0 1 (p 2 x + p 2 y ) 2 a b b c , (4) a
= 1 2 (p 2 y -p 2 x ) e -p 2 x -p 2 y -1 + p 2 x p 2 y + p 4 y , b = -p x p y [e -p 2 x -p 2 y -1 + p 2 x + p 2 y ], c = 1 2 (p 2 x -p 2 y ) e -p 2 x -p 2 y -1 + p 2 x p 2 y + p 4 x .
We note that Σ1 (p = 0) recovers the first order contribution to the conductivity in Eq. ( 2). At finite momentum, the self-energy contains off-diagonal elements, although the final correction to the conductivity is purely diagonal. The analytical computation of the non-crossing diagrams then proceeds as previously described, using Feynman's trick and Gaussian integration. For instance the following fourth order contribution

= d 2 p (2π) 2 K(p)ǫ Ĝ0 (p)[ Σ1 (p) Ĝ0 (p)] 3 ǫ
only involves a single momentum integration, which can then be performed analytically. Its value is given in Table II.

Extrapolation to the weak dissipation regime

We present here the methodology to obtain the extrapolated behaviour of the effective diagonal conductivity in the limit σ 0 → 0, starting from the large-σ 0 expansion: with the first six coefficients a n given in Table III.

σ xx (σ 0 ) = σ 0 1 + ∞ n=1 a n δσ 2 n σ 2n 0 (5)
One standard method of extrapolation is the so-called DLog Padé approximants [1], which starts with the dimensionless logarithmic derivative of the function to extrapolate:

f (x) ≡ σ 0 σ xx (σ 0 ) dσ xx (σ 0 ) dσ 0 σ0/ √ δσ 2 →x . (6)
One then reexpands at small x the function f (x) to order N :

f N (x) = 1 + N n=1 b n x 2n (7)
with the coefficients b n given in Table IV. The DLog Order 1 2 3 4 5 6 Coefficient bn -1 log(4) -2.135 -6.919 13.823 TABLE IV. Coefficients bn used in the DLog Padé extrapolation, corresponding to the small-x series expansion (7) of the function f (x) defined in Eq. ( 6).

Padé method uses then an approximant for f (x) of the following form:

f N (x) = 1 + N n=1 c n x 2n 1 + N p=1 d n x 2n . (8)
The coefficients c n and d n are computed from the knowledge of the perturbative terms b n given in Table IV. ¿From the expected power-law behavior of the conductivity at small dissipation, σ xx ∝ δσ 2 κ/2 σ 1-κ 0 , one gets f (x) → (1 -κ) for x → ∞. The critical exponent κ is thus obtained by extrapolating the Padé approximant (8) to infinity, which simply reads 1-κ = c N /b N at the order N .

The corrections to the effective conductivity at second order require an order N = 2 DLog Padé approximant, which lead after integration of Eq. ( 6) to the formula:

σ xx ≃ σ 0 1 + 1 κ δσ 2 σ 2 0 κ/2 (9)
with κ = 0.72 ± 0.09. The error bar on κ is obtained here by expanding Eq. ( 9) to third order with κ arbitrary, and comparing the deviation from the resulting coefficient with the exact a 3 value. Eq. ( 9) captures the full crossover between the perturbative regime δσ 2 ≪ σ 2 0 (where strong dissipation controls transport) to the nonperturbative limit of vanishing dissipation σ 0 → 0 (where percolation effects dominate), see Fig. 4 in the main text.

In order to obtain a better estimate for the exponent, one must push the calculation of the effective conductivity to fourth order. Following the same strategy, the order N = 4 DLog Padé approximant provides the estimate κ = 0.779 ± 0.006, and the resulting formula for the effective conductivity takes the form:

σ xx (σ 0 ) = σ 0 1 + A δσ 2 σ 2 0 B 1 + C δσ 2 σ 2 0 D (10)
with dimensionless numbers A, B, C, D, leading to κ = 2B + 2D. Again, the error bar on κ is obtained from comparison to the next known coefficient, namely a 5 , expanding Eq. ( 10) to fifth order while keeping an arbitrary κ fixed (a small additional error due to the Monte Carlo evaluation of the coefficients was also taken into account). While our calculation of the sixth order corrections to the conductivity would allow us in principle to further refine the estimation of the exponent, we encounter in that case a spurious pole [2], that invalidates the method. One explanation why the Padé method becomes unstable at high orders can be understood already from the fourth order extrapolation (10), which leads to trivial sub-leading corrections to scaling at small dissipation:

σ xx (σ 0 ) ∝ δσ 2 κ/2 σ 1-κ 0 1 + E σ 2 0 δσ 2 + . . . . (11)
This shows that the DLog Padé method enforces a given value κ ′ ≃ 3 -κ for the sub-leading exponent κ ′ , which is unlikely to correspond with good precision to the right value. This lack of flexibility is the likely source of the instability of the Padé approximant, and authors [3] have used a generalized n-Fit method that circumvents this problem. For the case of the fourth order conductivity, the fitting formula has rather the following additive form:

σ xx (σ 0 ) = F σ 0 1 + G δσ 2 σ 2 0 H +(1-F )c 0 1 + I δσ 2 σ 2 0 J . (12)
The critical exponent is then given by κ = min[2H, 2J], while the independent subleading exponent reads κ ′ = max[2H, 2J]. All unknown numerical coefficients are obtained by expanding Eq. ( 12) at small x and fitting to the coefficients of Table IV. Estimating the error by comparison to the known a 5 coefficient, we find κ = 0.767 ± 0.002, in excellent agreement with the conjectured value κ = 10/13 ≃ 0.7692. Moreover, the Padé approximants show good convergence for all values of the dissipation strength σ 0 , see Fig. 4 in the main text.

High temperature microscopics of σxx at the plateau transition

The local Hall conductivity can be computed microscopically in the high magnetic field regime [4,5], and simply follows from Hall's law with Landau level quantization:

σ H (r) = e 2 h ∞ m=0 n F [E m + V (r) -µ] (13)
with standard Landau levels E m = ω c (m + 1/2), cyclotron frequency ω c = |eB|/m * , random disorder potential V (r), chemical potential µ and Fermi function n F (E) = 1/{exp[E/(k B T )] + 1}. In particular, microscopic calculations [5] show that deviations to the form (13) are small by the dimensionless parameter [l 2 B

V 2 ]/[ξ 2 ω c ] ≪ 1, with V 2 the width of the disorder distribution, l B = h/|eB| the magnetic length, and ξ the large correlation length of the disorder fluctuations. Note the smallness of l B ≃ 8nm at B = 10T, so that l B ≪ ξ for smooth disorder.

At temperatures such that T ≫ V 2 , the Fermi distribution in Eq. ( 13) can be linearized, so that Gaussian fluctuations of disorder provide Gaussian fluctuations for the Hall conductivity σ H (r) = σ H + δσ(r) with

σ H = e 2 h ∞ m=0 n F (E m -µ), (14)
δσ(r) = e 2 h ∞ m=0 n ′ F (E m -µ)V (r). (15)
The power-law behavior of the longitudinal conductivity at small dissipation, σ xx = Cσ 1-κ 0 δσ 2 κ/2 , leads to:

σ xx = Cσ 1-κ 0 e 2 h V 2 ∞ m=0 n ′ F (E m -µ) κ . ( 16
)
We re-express the sum over Landau levels in Eq. ( 16) by using Poisson summation formula:

+∞ m=0 f (m) = +∞ l=-∞ +∞ 0 dt e i2πlt f (t). ( 17
)
In the limit T < µ, one finds after standard manipulations [6]:

∞ m=0 n ′ F (E m -µ) = 1 ω c 1 + +∞ l=1 (-1) l cos 2πlµ ω c × 4π 2 lkB T ωc sinh 2π 2 lkB T ωc . (18)
Finally, by considering the plateau transition region between the filling factors ν and ν + 1, the chemical potential is pinned to µ = ω c (ν + 1/2), leading to expressions (17)- (18) in the main text.

FIG. 1 .

 1 FIG. 1. (color online) Two-dimensional sample with percolating random charge inhomogeneities: measurement of longitudinal Vxx and Hall Vxy voltages with applied current I.

FIG. 3 .

 3 FIG. 3. Diagrammatic expansion in the case of Gaussian fluctuations of the local conductivity. Wiggly lines denote the conductivity correlation functions.

σ 0 4 FIG. 4 .

 044 FIG. 4. (color online) Scaling function of the longitudinal conductance from the percolating (σ0 → 0) to the dissipative regime (σ0 → ∞). A comparison is made between first-order bare perturbation theory (controlled only at large dissipation) to the resummations of the n = 2 and n = 4 orders, showing good convergence for all values of σ0.

  FIG. 5. (color online) Temperature dependence of the peak longitudinal conductance from Eq. (17) in log-scale. A crossover occurs at kBT ⋆ = ωc/4 between a low-T powerlaw [σxx ∝ T -κ or T1-2κ ] and a high-T power-law [σxx ∝ cst or T 1-κ ], for elastic scattering σ0(T ) = σ0 or inelastic phonon scattering σ0(T ) = AT respectively.

  . Our most accurate result

	Order Method Exponent 1 -κ
	2	Padé	0.28 ± 0.09
	4	Padé	0.221 ± 0.006
	4	n-fit	0.233 ± 0.002
	∞ Conjecture 3/13 ≃ 0.2308
	TABLE I. Critical exponent 1 -κ obtained from Padé ap-
	proximants

TABLE I .

 I Diagrams to second and third order: multiplicity and analytical values. The resulting coefficients a4 and a6 are given in TableIII.

	Diagram	Multiplicity Analytical Value	Decimal Value
	second order
		1 -1 4 log(2)		-0.173287
		1 1 8 (1 -log(4))	-0.0482868
	third order
		1 1 96 3 -π 2 + 3 log[3](-3 + log[9]) + 12Polylog 2, 2 3	0.00504001
		2 1 32 log 27 16		0.0163515
		1 1 16 2 log[2] 2 -3 log[3] + log[8] + Polylog 2, 1 4	0.000760209
		2 1 384 (2 + 100 log[2] -63 log[3])	0.00547433
		1 1 8 log 32 27		0.0212374
		1 1 8 log 27 16		0.065406
		1 -1 48 -log[2] 6	+ 9 log[3] 64	0.0181345
		1 3 16 log 4 3		0.0539404

TABLE II .

 II Benchmarking the numerical Monte Carlo evaluation against analytically tractable non-crossing diagrams at fourth, fifth and sixth order respectively.

	Order	Diagram	Analytical value		Monte Carlo evaluation
	4		-44 log[2]+27 log[3] 32	≃ -0.02612	-0.02607
	5		162 log[3]+125 log[5]-544 log[2] 192	≃ 0.01084	0.01087
	6		-6496 log[2]-486 log[3]+3125 log[5] 1536	≃ -0.004632 -0.004630
	Order Method Coefficient an
	1 Analytical 1 2
	2 Analytical 1 8 -1 2 log(2)
	3 Analytical 0.2034560502
	4 Numerical -0.265 ± 0.001
	5 Numerical 0.405 ± 0.001
	6 Numerical -0.694 ± 0.001

TABLE III .

 III Coefficients an of the perturbative series (5) up to sixth loop order.

computation, and S. Bera, B. Piot, M. E. Raikh, V. Renard and F. Schoepfer for stimulating discussions.


Loïc Iapteff

Julien Jacques

Benoit Celse

Victor Lameiras

Franco Da Costa

Design of Experiment for Bayesian transferred model

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Design of Experiment for Bayesian transferred model

Loïc Iapteff 1 , Julien Jacques 2 , Benoit Celse 3 & Victor Costa 4 1 loic.iapteff @ifpen.fr, Rond-point de l'échangeur de Solaize, 69360 Solaize 2 julien.jacques @univ-lyon2.fr, 5 Avenue Pierre Mendès France, 69500 Bron 3 benoit.celse @ifpen.fr, Rond-point de l'échangeur de Solaize, 69360 Solaize 4 victor.costa @ifpen.fr, Rond-point de l'échangeur de Solaize, 69360 Solaize Résumé. Le groupe IFP commercialise des catalyseurs et doit s'engager sur leur performance. Il est donc nécessaire de disposer de modèles prédictifs fiables pour chaque nouvelle génération de catalyseurs. Ces modèles sont construits à partir de données expérimentales très couteuses. Afin d'optimiser les coûts, notre ambition est de réduire le nombre d'expérimentations nécessaires pour estimer un modèle associé à un nouveau type de catalyseur. De précédents travaux ont montré qu'une approche de transfert Bayésien permettait d'améliorer la qualité des modèles lorsque le nombre d'observations est réduit. Dans cet article, les plans d'expériences sont étudiés afin de déterminer comment sélectionner ce nombre réduit d'observations permettant d'obtenir les meilleurs modèles par transfert Bayésien. Cet article montre que l'algorithme Kennard and Stone peut, sous certaines conditions, offrir de meilleurs résultats que des plans optimaux.

Mots-clés. Plan d'expérience, Transfer learning, inférence bayésienne Abstract. IFP group develops catalysts and has to guarantee their performances. It is therefore crucial to have good predictive models for all new catalysts. These models are built upon very expensive experimental data. In order to minimize costs, we aim at reducing the number of new data points to measure to fit a model on the new catalyst. Previous work has shown that a Bayesian transfer approach can improve the quality of models when the number of observations is reduced. In this paper, experimental designs are studied in order to determine how to select this reduced number of observations to obtain the best models by Bayesian transfer. This article shows that the Kennard and Stone algorithm can, under certain conditions, offer better results than optimal designs.

Keywords. Design of Experiment, Transfer learning, Bayesian inference 1 The challenge IFP Group develops and sells catalysts to chemical and biochemical producers. Catalysts are solids that make the reaction feasible, faster, and/or at a lower temperature and pressure. Performance must be guaranteed and it is therefore crucial to have good predictive 1 models for all new catalysts. These models are built on very expensive experimental data, therefore the aim is to use a small number of data and select them carefully to build a good model. For the modeling of hydrocracking process, previous work has shown that the Bayesian transfer approach significantly reduces the number of points needed to obtain robust models [START_REF] Iapteff | Reducing the number of experiments required for modelling the hydrocracking process with kriging through Bayesian transfer learning[END_REF]. The method consists in estimating a model for a new catalyst using Bayesian models, whose priors depend on the old catalyst generation, for which important amounts of data are available. Celse, J.-J. D. Costa, and V. Costa 2016 shows the effectiveness of D-optimal designs for modeling the hydrocracking process using the kinetic model.

This paper focuses on hydrocracking modeling using Bayesian linear model (1):

y = βX + ϵ, (1) ϵ ∼ N (0, σ 2 I), β ∼ π(•),
where β = (β 0 , β 1 , . . . , β p ) T is the vector of model parameters, assumed to be a random variable of prior density π(•), X is the matrix of observations and y is the variable of interest. For the experiments, an industrial dataset from hydrocracking facilities is available. The prior π(•) is built with an old dataset that uses another catalyst.

The goal of this work is to show that when the sample size is not too small, basic geometric experimental design offer better results than optimal designs. This paper is composed of two parts: a first section presents the Design of Experiment (DoE) and a second section shows the results on our industrial dataset, using some of the most commonly employed DoE approaches. This section also explains why the results are not necessarily the best when using optimal designs.

Design of Experiment

The purpose of the Design of Experiments (DoE) is to organise the experiments carried out for scientific research or industrial studies in the best possible way. When experiments are used to study the relation between a quantity of interest and other features, DoE aims to maximise the information on this quantity from a reduced amount of data. The benefits of DoE are multiple, for instance it allows to control bias, to reduce random variation, to increase the precision of parameter estimates, to improve predictions of future observations, or to make a choice between competing models. In this section, we focus on most common DoE approaches for linear models.

Model independent approaches

The Kennard and Stone algorithm [START_REF] Kennard | Computer aided design of experiments[END_REF]) is an iterative algorithm for selecting the best observations regardless of the type of model. This algorithm is a sequential construction for a minimax space-filling design and is based on a set of possible observations in order to extract the optimal points. It is initialized with a training dataset composed of the two most distant points of the complete dataset and a candidate set composed of the remaining points. Then, at each iteration, one observation x i is moved from the candidate set to the training set such as x i = argmax

x i ∈Candidate ( min x i ′ ∈Training dist(x i , x i ′ ))
until the training dataset reaches a predefined size. Notice that in the original Kennard and Stone algorithm, the euclidean distance is considered. With such an approach, we try to obtain an uniform distribution over the features space for the training set.

Other DoE model-independent approaches are well studied, such as modified Kennard and Stone algorithm, Latin hypercube sampling or full factorial designs but are not developped in this paper. The model independent DoE approaches are also often used to initialise algorithms for obtaining optimal designs, allowing faster convergence.

Optimal DoE for Linear Model

In this section, DoE approaches for the linear model (1) is considered. For this model, the most commonly used criteria are the D/A/E/G-optimal designs [START_REF] Fedorov | Theory of Optimal Experiments Designs[END_REF][START_REF] De Aguiar | D-optimal designs[END_REF][START_REF] Wong | Comparing robust properties of A, D, E and G-optimal designs[END_REF]. When the objective is to improve the prediction's quality of the model, a popular criterion is the G-optimal criterion, which aims to minimize the maximum prediction uncertainty: the design must minimize max

x i (x T i (X T X) -1 )x i ).
The D, A and E optimal criteria use the dispersion matrix (X T X) -1 which is proportional to the covariance matrix of parameters. The idea is that minimizing the uncertainty in the parameters helps to minimize the uncertainty in prediction.

The D-optimal design is the design that minimizes det((X T X) -1 ). It is the most widely used criterion due to its efficiency and ease of implementation. The Fedorov's algorithm [START_REF] Fedorov | Theory of Optimal Experiments Designs[END_REF] gives an exact D-optimal design with an easy implementation. This algorithm ensures convergence to a local minimum, so different initializations can be considered. An usual initialization is to take a design obtained via the Kennard and Stone algorithm.

The A-optimal design and E-optimal design minimize the trace of (X T X) -1 and the largest eigenvalue of (X T X) -1 , respectively.

Bayesian Experimental Design

Since the objective of DoE is often to correctly estimate the model parameters, it can be interesting to provide a prior information on the value of these parameters. This approach is called Bayesian Experimental Design [START_REF] Chaloner | Bayesian experimental design: A review[END_REF] and also has the advantage that a single observation can be sufficient to build the model, the information matrix being always non-singular with a proper prior. Bayesian Experimental Design is commonly based on an utility function U (d, θ, y), defining the contribution of the choice of design d, yielding the y-values for the dependent variable and with the model's parameters θ. The optimal design d opt is thus obtained by maximizing the maximum expected utility:

d opt = argmax d Y Θ U (d, θ, y)p(y|θ, d)p(θ)dθdy, (2)
where p(y|θ, d) is the likelihood of the model and p(θ) the prior information, which is assumed not to be affected by the choice of the design.

In the Bayesian linear regression case, with a Gaussian prior for parameters β ∼ N (µ, Σ) and σ fixed, the posterior distribution of the parameters is known:

β|y, X ∼ N ((X T X + σ 2 Σ -1 ) -1 (X T y + σ 2 Σ -1 µ), σ 2 (X T X + σ 2 Σ -1 ) -1 ).
In our case of study, the parameters β are considered uncertain and σ is assumed to be known and identical to the parameter of the source model since we seek to model similar phenomena with identical tolerance.

Equivalently to the classical DoE framework, Bayesian alphabetical criteria are developed and use the matrix (X T X+σ 2 Σ -1 ) -1 , proportional to the covariance of the Bayesian linear model parameters, instead of the matrix (X T X) -1 . Then, Bayes D/A/E/G-optimal designs minimize respectively det((X T X + σ 2 Σ -1 ) -1 ), the trace of (X T X + σ 2 Σ -1 ) -1 , the largest eigenvalue of (X T X + σ 2 Σ -1 ) -1 and max

x i (x T i (X T X + σ 2 Σ -1 ) -1 x i ).
The Bayesian D-optimal design can be derived from the mutual information for the Gaussian linear model while the Bayesian A-optimal design can be obtained with the quadratic loss. Concerning the E-optimal design and G-optimal design, no utility function has yet been found whose maximisation leads to the construction of such a model.

Notice that Bayesian DoE is not only used if one wishes to build a Bayesian predictive model but can also be used in a non-Bayesian framework.

Use case: hydrocraking dataset

The dataset used for this experiment is obtained from two refineries using a new generation catalyst and is composed of 1004 observations described with 12 features. Let X be the matrix of observations and y the vector of the actual value of Diesel density to be predicted. In this experiment, different DoE are tested for a varying size of sample n. The aim is to build a design that allows to obtain the most efficient models with the fewest possible points for learning. The objective is to combined the DoE approach with the Bayesian transfer method [START_REF] Iapteff | Reducing the number of experiments required for modelling the hydrocracking process with kriging through Bayesian transfer learning[END_REF]), which has proven its efficiency in improving the quality of the model. The model fitted with the Bayesian transfer approach is compared to a classical "non-transfer" model, for which β is assumed to be no longer random and estimated by maximum likelihood. The "Bayesian transfer model" is the model learned with prior on β: β ∼ N (β prior , Σ prior ). As a reminder, for the Bayesian approach, the prior is built from a dataset from refineries using an old generation catalyst. To compare different methods, the target dataset is split: 600 observations are used as candidate for the design and the 404 remaining observations are used for testing the model. The RMSE score is thus always evaluated only on these 404 points whatever the size n of the sample.

For both models with and without transfer, three different designs are considered:

• Mean of random designs: samples are randomly selected from the candidates dataset, a thousand times, and scores are then averaged.

• Kennard and Stone,

• D-optimal design: the Fedorov's algorithm is used in order to build the D-optimal design for each sample size. For the model without transfer, classical D-optimal design is considered while Bayesian D-optimal design is considered for the model fit with Bayesian transfer.

Note that other alphabetical criteria have been tested and offer less good results than the D-optimal criterion for our prediction task and are therefore not presented here. The results are plotted in Figure 1. For all designs, the Bayesian transfer approach improves the results compared to a non-transfer approach. Nevertheless, the non-transfer approach with a D-optimal design provides very good results fairly quickly . The best results with a reduced number of points are obtained with a D-optimal design coupled with the Bayesian transfer approach. As the number of points increases, above 15, the Bayesian transfer approach coupled with the design obtained via Kennard and Stone is better.

In the present case study, the good performance of the Kennard and Stone method in comparison with the D-optimal approach can be explained by the fact that data came from two refineries. Indeed, the dataset is then not exactly homogeneous, and the discrete D-optimal design may be far from the theoretical optimum. Additional experiments on simulated dataset confirmed this behavior.

Another important point is that the design obtained with Kennard and Stone algorithm allows to reach the best possible score when number of observation is sufficient, around 25 observations for our application with 12 features.

Conclusion

This study shows that in some practical cases, building optimal designs is not the best solution. It is the case when the features space distribution is not homogeneous and the chosen model does not perfectly describes the phenomenon. Then, the Kennard and Stone algorithm allows to keep the important information when the number of observation is sufficient, and become more effective than optimal design. However, when the number of observations is low, D-optimality remains the most efficient.

Figure 1 :

 1 Figure 1: Comparison of the different designs for the Linear model. The red lines give the best possible design obtained by taking for each sample size n the best RMSE score among the thousand random samples.


Benoît Liquet
email: benoit.liquet@isped.u-bordeaux2.fr

Jérôme Saracco
email: jerome.saracco@u-bordeaux4.fr

Application of the bootstrap approach to the choice of dimension and the α parameter in the SIR α method

Keywords: Bootstrap, Dimension Reduction, Sliced Inverse Regression. 1

To reduce the dimensionality of regression problems, sliced inverse regression approaches make it possible to determine linear combinations of a set of explanatory variables X related to the response variable Y in general semiparametric regression context. From a practical point of view, the determination of a suitable dimension (number of the linear combination of X) is important. In the literature, statistical tests based on the nullity of some eigenvalues have been proposed. Another approach is to consider the quality of the estimation of the effective dimension reduction (EDR) space. The square trace correlation between the true EDR space and its estimate can be used as goodness of estimation. In this paper, we focus on the SIR α method and propose a naïve bootstrap estimation of the square trace correlation criterion. Moreover, this criterion could also select the α parameter in the SIR α method. We indicate how it can be used in practice. A simulation study is performed to illustrate the behaviour of this approach.

Introduction

Conventional parametric or nonparametric methods of curve fitting estimate the unknown relationship between an explanatory variable and a response variable. These methods work well when the dimensionality of the data is low. In high dimensions, all data are sparse. Thus the number of observations available to give information about the local behaviour of the regression function becomes very small with large dimensions. Moreover, the complexity of the possible underlying structure increases more than exponentially fast as the dimension increases, so enormous numbers of data are needed in order to decisively select one parametric model over another. This is the so-called curse of dimensionality which quickly defeats parametric or nonparametric regression analysis for fitting usefully predictive models. In the statistical literature, various strategies exist to challenge its effect. One of them is to consider dimension reduction models.

In the following, we consider a high-dimensional data set with one response variable Y and many predictor variables X 1 , . . . , X p . Let X = (X 1 , . . . , X p ) T denote the p-dimensional column vector of explanatory variables. Parametric approaches such as linear and polynomial regression, or more generally, transformation or generalized linear models, are examples of dimension-reduction methods. Other methods aim to relieve the problem of the curse of dimensionality without requiring assumptions as strong as those of the parametric ones. Such methods are nonparametric in nature with a parametric part and are called semiparametric. Many dimension-reduction tools are available. There are essentially based on two families of dimension-reduction models. In the first, the effects of the p predictor variables are assumed to be additive: the generalized additive models belong to this family, see for instance Hastie andTibshirani (1986, 1990), [START_REF] Stone | Additive regression and other nonparametric models[END_REF][START_REF] Stone | The dimensionality reduction principle for generalized additive models[END_REF] or [START_REF] Friedman | Flexible parsimonious smoothing and additive modeling (with discussion)[END_REF]. Among the well-known methods, there is the ACE (Alternating Conditional Expectation) introduced in the Brieman and [START_REF] Breiman | Estimating optimal transformations for multiple regression and correlation (with discussion)[END_REF] algorithm and the MARS (Multivariate Adaptative Regression Splines) method described by [START_REF] Friedman | Multivariate adaptative regression splines (with discussion)[END_REF].

The second family assumes that the effects of the explanatory variable X can be captured in a K (with K ≤ p) dimensional projection subspace, (X T β 1 , . . . , X T β K ) where β 1 , . . . , β K ∈ ℜ p . The popular PPR (Projection Pursuit Regression) approach belongs to those projection-based models.

It was introduced by [START_REF] Friedman | Projection pursuit regression[END_REF] and has been discussed by many authors, see for instance [START_REF] Chen | Estimation of a projection-pursuit type regression model[END_REF] or [START_REF] Hall | On projection pursuit regression[END_REF] among others. However, most of the dimension reduction methods mentioned so far involve a backfitting (iterative) algorithm which needs a smoothing curve-fitting at each step. These methods generally are computing-intensive and may not converge.

Several methods which make it possible to estimate the directions β k first, without a smoothing step of the link functions, have been proposed by [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] and [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF][START_REF] Li | On principal Hessian directions for data visualization and dimension reduction: another application of Stein's lemma[END_REF]. These are the SIR (Sliced Inverse Regression) and PHD (Principal Hessian Direction) methods. They are computationally simple, since they require only simple averaging and an eigenvalue decomposition.

In this paper, we will concentrate on the so-called SIR methods and more particulary on the SIR α method.

The dimension-reduction assumption is the following: Y ⊥ X|X T β 1 , . . . , X T β K where the notation U ⊥ V |W means that the random variable U and V are independent conditionally to the random variable W . The underlying semiparametric regression model can be written as :

Y = f (X T β 1 , . . . , X T β K , ǫ), (1)
where the response variable Y is associated with the p-dimensional regressor X only through the linear combinations X T β k and ǫ is a random error term independent of X. No assumption is made about the functional form of the unknown link function f or the distribution of ǫ. Since no structural conditions on f are imposed, the vectors β k are not identifiable, unlike the linear subspace spanned by the β k 's. [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] used this model to introduce the notion of EDR (Effective Dimension Reduction) space, namely the linear subspace E of R p spanned by the β k 's. Any vector in E is called an EDR direction. This semiparametric regression model appears to be a reasonable compromise between fully parametric and fully nonparametric modeling. Two estimation problems in (1) clearly arise. The first is to recover the EDR space. The second consists in estimating the link function f . Since {(X T β 1 , . . . , X T β K ) : X ∈ R p } theoretically captures all the information on the distribution of Y given X and can be estimated first, one can then easily apply standard nonparametric smoothing techniques (such as kernel, spline or wavelet regression) to estimate the link function f . The goal of dimension reduction is thus achieved.

In this paper our main objective concerns dimension K of the EDR space. Obviously, in most applications, dimension K is unknown and hence must be estimated from the data {(X i , Y i ), i = 1, . . . , n} with p < n. Several approaches have been proposed in the literature. For the original approach introduced by Li (1991), named SIR-I in the following, and for SIR-II (see Scott, 1994), the choice of dimension K is based on statistical nested test procedures. However, from a practical point of view, such test procedures present two main drawbacks. First, the validity of these test approaches rely on normality distribution or elliptically symmetric distribution assumptions for X. Secondly, the nested test procedures do not allow to control the overall level since knowledge of the level at each step does not imply knowledge of the overall level because of the problem of multiplicity of tests. Note that [START_REF] Bai | A chi-square test for dimensionality for non-Gaussian data[END_REF] studied the limiting distribution of the test statistic for dimensionality without the normality assumption, and obtained a necessary and sufficient condition for the chi-square limiting distribution to hold. However, this approach focused on a variant of SIR-I, called CANCOR and introduced by Fung et al. (2002).

There is an alternative to the test approach: this kind of approach is based on the trace criterion and was firstly developed by [START_REF] Ferré | Dimension choice for sliced inverse regression based on ranks[END_REF][START_REF] Ferré | Determining the dimension in Sliced Inverse Regression and related methods[END_REF] in the Sliced Inverse Regression and Principal Hessian direction context. In this paper, as in [START_REF] Ferré | Determining the dimension in Sliced Inverse Regression and related methods[END_REF], we prefer a model selection approach which does not require the estimation of the link function f (depending itself on the selected dimension). Thus, while a criterion based on the prediction of Y would be optimal, we only focus on the projections to choose the dimension on the EDR space, without including the estimation of f .

We consider a sample s = {(X i , Y i ), i = 1, . . . , n} from the model (1). Let Σ = V(X) be the covariance matrix of X, and let Σ be the empirical covariance matrix of the X i 's. Let B =

[β 1 , . . . , β K ] denote the p × K matrix of the true β k 's assumed to be linearly independent. Clearly, when the true dimension is K, we have Span(B) = E. Let B = [ b1 , . . . , bK ] denote the corresponding estimated EDR directions matrix. Let us define E = Span( B) the estimated EDR space. To study the closeness between two subspaces, the square trace correlation can be naturally used. The corresponding risk function is defined by considering the following expectation:

R k = E Trace(P k P k ) /k, (2)
where P k denotes the Σ-orthogonal projector onto the space spanned by the first k vector β l of B and P k is the Σ-orthogonal projector onto the space spanned the first k vector bl of B.

More precisely, let

B k = [β 1 , . . . , β k ] and B k = [ b1 , . . . , bk ], then P k = B k (B T k ΣB k ) -1 B T k Σ and P k = B k ( B T k Σ B k ) -1 B T k Σ.
Note that R k is defined for any dimension k lower than or equal to K. A value of R k close to one indicates that the set of the k estimated linear combinations of X is close to the ideal set. So in terms of dimensionality, k is a feasible solution. On the other hand, a value of R k perceptibly different from 1 means that this estimated set is slightly different from the ideal one, so the solution for the dimension is greater than k. Since R K will converge to one as n tends to infinity (for the true dimension K), then, for a fixed n, a reasonable way to assess whether an EDR direction is available is given by looking at how much R k departs from one. From a computational point of view, we require consistent estimates R k of R k , so the feasible solution for the dimension can be obtained by computing the values of R k for k = 1 to p and observing how much it departs from one. In [START_REF] Ferré | Dimension choice for sliced inverse regression based on ranks[END_REF][START_REF] Ferré | Determining the dimension in Sliced Inverse Regression and related methods[END_REF], consistent estimators of R k were proposed. They were obtained from asymptotic expansions of the loss function R k which require the normality or the elliptically symmetric distribution of X. Moreover, the estimate is not easily computable.To avoid complex computing, [START_REF] Ferré | Dimension choice for sliced inverse regression based on ranks[END_REF] suggested considering the SIR approach of [START_REF] Hsing | An asympotic theory for Sliced Inverse regression[END_REF] based on slices of cardinal 2 and to use a Jack-Knife (JK) method to estimate R k in this special case: the corresponding estimated loss function measuring how far the estimated subspace obtained by deleting two observations (that is a slice) differs from the sample estimated over the whole sample. However, Ferré mentioned that the Hsing and Carroll estimators, which use slices with two observations, are unable to capture some relevant information, see also [START_REF] Aragon | Sliced Inverse Regression (SIR): an appraisal of small sample alternatives to slicing[END_REF] for a simulation study exhibiting the undersmoothing effect of these estimators. Hence, the efficiency of the corresponding JK criterion to choose dimension K is not really clear from a practical point of view.

In this paper, we use this criterion to determine a suitable dimension K using the SIR α approach.

The SIR α method (described in the next section) depends on a parameter α which controls the combination of the SIR-I and the SIR-II methods. Therefore, the SIR α approach will be used to define P k , which depends on the choice of α and will be henceforth denoted P k,α . The practical choice of α can be based on the test approach proposed by [START_REF] Saracco | Pooled Slicing methods versus Slicing methods[END_REF] which does not require estimation of the link function. Two cross-validation criteria were also developed by [START_REF] Gannoun | Two Cross Validation Criteria for SIR α and PSIR α methods in view of prediction[END_REF] to select the parameter α, but these criteria require kernel smoothing estimation of the link function.

To determine the suitable α parameter and the suitable dimension K, we use the criterion defined in (2), which now becomes R k,α . Criterion R k,α is precisely defined in (5). To compute estimates for this risk function, computer intensive methods can be used. Note that R k,α is symmetrically defined since its value is invariant under any permutation of the observations (X i , Y i ) of s. We propose here to use a well-known resampling method, the bootstrap, which is a natural candidate for estimating R k,α . In Section 2, we give an overview of the SIR α approach. The bootstrap estimate of R k,α is described in Section 3. Section 4 is devoted to a simulation study which was conducted in order to show the efficiency of the proposed approach. A comparison with existing methods is also provided in Section 4.4. Finally, some concluding remarks are given in Section 5.

Overview of the SIR α approach

We give an overview of the univariate SIR approaches (that is when Y ∈ ℜ). While there are several possible variations, the basic principle of SIR methods (SIR-I, SIR-II or SIR α ) is to reverse the role of Y and X. Instead of regressing the univariate response variable Y on the multivariate covariable X, the explanatory variable X is regressed on the dependent variable Y . The SIR-I estimates, based on the first moment E(X|Y ), were introduced by [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] and [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF], and have been studied extensively by several authors: see for instance [START_REF] Carroll | Measurement error regression with unknown link: dimension reduction and data visualization[END_REF], [START_REF] Hsing | An asympotic theory for Sliced Inverse regression[END_REF], [START_REF] Zhu | Asymptotics of sliced inverse regression[END_REF], or [START_REF] Aragon | Sliced Inverse Regression (SIR): an appraisal of small sample alternatives to slicing[END_REF]. The estimation scheme and the application of SIR-I have been discussed in detail by [START_REF] Chen | Can SIR be as popular as multiple linear regression[END_REF]. Asymptotic properties of SIR have been investigated in several articles with emphasis on the convergence and the asymptotic distribution of the estimator of the EDR space, see for instance [START_REF] Hsing | Nearest neighbor inverse regression[END_REF]Carroll (1992), Koetter (1996), [START_REF] Zhu | Asymptotics for kernel estimate of Sliced Inverse Regression[END_REF], [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF], [START_REF] Hsing | Nearest neighbor inverse regression[END_REF], among others. However, this approach is "blind" for symmetric dependencies (see [START_REF] Cook | Discussion of "Sliced inverse regression[END_REF] or [START_REF] Kötter | Sliced Inverse Regression[END_REF]). Therefore, SIR-II estimates based on the inverse conditional second moment V(X|Y ) have been suggested, see for instance [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF], [START_REF] Cook | Discussion of "Sliced inverse regression[END_REF], [START_REF] Kötter | Sliced Inverse Regression[END_REF] or [START_REF] Yin | Asymptotic distributions for dimension reduction in the SIR-II method[END_REF]. Hence these two approaches concentrate on the use of the inverse conditional moments E(X|Y ) or V(X|Y ) to find the EDR space. To increase the chance of discovering all the EDR directions, the idea of the SIR α method is to conjugate the information provided by SIR-I and SIR-II methods: if an EDR direction can only be marginally detected by SIR-I or SIR-II, SIR α considers a mixture of these two methods (see the matrix M α defined below) which may provide all the EDR directions.

Let us now recall the geometric properties of model (1). In SIR α approach, [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] considered, for α ∈ [0, 1], the eigen-decomposition of the matrix

Σ -1 M α where M α = (1-α)M I Σ -1 M I +αM II .
The matrices M I and M II are respectively the matrices used in the usual SIR-I and SIR-II approaches. They are defined as follows:

M I = V(E(X|T (Y ))) and M II = E (V(X|T (Y )) -E(V(X|T (Y )))) Σ -1 (V(X|T (Y )) -E(V(X|T (Y ))))
T where T denotes a monotonic transformation of Y . Specific transformations T are ordinarily used in order to simplify the expression of the matrices M I and M II . The most usual one is the slicing function (defined below and used in the rest of the paper). Note that the monotonicity of T is necessary in order to avoid pathological case of SIR-I only due to a bad choice of T . It can be shown that, under the linearity condition (3) and the constant variance assumption defined in Remark 2, and for any monotonic transformation T , the eigenvectors associated with the largest K eigenvalues of Σ -1 M α are some EDR directions. Note also that, when α = 0 (resp. α = 1), SIR α is equivalent to SIR-I (resp. SIR-II). [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] proposed a transformation T , called a slicing, which categorizes the response Y into a new response with H > K levels. The support of Y is partitioned into H non-overlapping slices s 1 , . . . , s h , . . . , s H . With such a transformation T , the matrices of interest are now written as

M I = H h=1 p h (m h -µ)(m h -µ) T and M II = H h=1 p h V h -V Σ -1 V h -V , where µ = E(X), p h = P (Y ∈ s h ), m h = E(X|Y ∈ s h ), V h = V(X|Y ∈ s h ) and V = H h=1 p h V h .
Therefore, it is straightforward to estimate these matrices by substituting empirical versions of the moments for their theoretical counterparts, and therefore to obtain the estimation of the EDR directions. Each estimated EDR direction converges to an EDR direction at rate √ n, see for instance [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] or [START_REF] Saracco | Pooled Slicing methods versus Slicing methods[END_REF]. Asymptotic normality of the SIR α estimates has been studied

by Gannoun and Saracco (2003a).

Remark 1. The practical choice of the slicing function T is discussed in [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF], [START_REF] Kötter | Sliced Inverse Regression[END_REF] and [START_REF] Saracco | Pooled Slicing methods versus Slicing methods[END_REF]. Note that the user has to fix the slicing strategy and the number H of slices, then observations are assigned to slices by value. The SIR theory makes no assumption about the slicing strategy. In practice, there are naturally two possibilities: to fix the width of the slices or to fix the number of observations per slice. In their investigation of SIR-I, various researchers have preferred the second approach. From the sampling point of view, the slices are such that the number of observations in each slice is as close to each other as possible. To avoid artificial reduction of dimension, H must be greater than K. Also, in order to have at least two cases in each slice, H must be less than [n/2] where [a] denotes the integer part of a. [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] noticed that the choice of the slicing is less crucial than the choice of a bandwidth, as in kernel-based methods. Simulation studies (with p = 5 and 10) show that the influence of the "slicing parameter" is small when the sample size is greater than 100. Note that, in order to avoid the choice of a slicing (the transformation T is here the identity function and not the slicing function), the kernel-based estimate of SIR-I has been investigated, see [START_REF] Zhu | Asymptotics for kernel estimate of Sliced Inverse Regression[END_REF] or [START_REF] Aragon | Sliced Inverse Regression (SIR): an appraisal of small sample alternatives to slicing[END_REF]. However, these methods are hard to implement with regard to basic slicing one and are computationally slow. Moreover, [START_REF] Bura | Dimension reduction via parametric inverse regression. L 1 -statistical procedures and related topics (Neuchel[END_REF] and [START_REF] Bura | Estimating the structural dimension of regressions via parametric inverse regression[END_REF] proposed a parametric version of SIR-I.

Remark 2. Note that two crucial conditions for the theoretical succes of SIR α methods are the following: a linearity condition

E(b T X|β T 1 X, . . . , β T K X) is linear in β T 1 X, . . . , β T K X; ∀b ∈ ℜ p , (3)
and a constant variance condition

V(X|β T 1 X, . . . , β T K X) is non-random. ( 4
)
Condition ( 3) is satisfied when X has an elliptically symmetric distribution and ( 4) is satisfied when X follows a multivariate normal distribution (which is an elliptically one): in this case, the linearity condition is satisfied for any directions, even if they may not be in the EDR space. Note that it is not possible to verify (3) since this involves the unknown directions β k , contrary to ellipticity or normality condition which can be tested. Using the bayesian argument of [START_REF] Hall | On almost linearity of low dimensional projections from high dimensional data[END_REF], we can infer that (3) holds approximately for many high dimensional data sets. An interesting and detailed discussion on the linearity condition can be found in [START_REF] Chen | Can SIR be as popular as multiple linear regression[END_REF]. Note that when the distribution of X is far from being multinormal (or from conditions (3) and ( 4)), the dimension reduction may not be possible with SIR approaches: that is the choice of K (based on the eigenvalues or on R k ) will give K = p.

3 Bootstrap estimate of the risk function R k,α

The risk function is defined by:

R k,α = E Trace(P k P k,α ) /k, ∀k = 1 . . . , K (5)
where P k,α is the Σ-orthogonal projector onto the space spanned by the first k vector bl,α of

B α = [ b1,α , . . . , bK,α ],
with bl,α the eigenvector associated to the lth greater eigenvalues of Σ

-1 M α the empirical version of Σ -1 M α which is obtained from the sample s = {(X i , Y i ), i = 1, . . . , n}.
Let B be the number of bootstrap replications and let

s (b) = (X (b) i , Y (b)
i ), i = 1, . . . , n a nonparametric bootstrap sample replication. According to [START_REF] Efron | The jackknife, the bootstrap and other resampling plans[END_REF], a naïve bootstrap estimate of the mean square risk function is defined by:

R k,α = 1 B B b=1 R (b) k,α (6)
where

R (b) k,α = Trace P k,α P (b) k,α /k and P (b)
k,α is the projector onto the subspace spanned by the first k eigenvectors of the matrix of interest, which is obtained from the bootstrap replication sample s (b) . Note that the practical criterion R k,α could be computed for all k = 1 . . . , p whereas the R k,α is only defined for k = 1, . . . , K.

The aim is to find the dimension k of the model and to have a practical choice of α in the SIR α method thanks to the bootstrap estimate. The proposed method consists in evaluating the R k,α for all (k, α) ∈ {1, . . . , p} × [0, 1] and then in observing how much it departs from one. The best choice will be a couple ( K, α) which gives a value of R k,α close to one, such that K << p. In practice, there is no objective criteria in order to establish when a departure from one is close, but a visual expertise of the plot of the Rk,α versus k and α allows us to choose the best couple. This point will be illustrated with simulations in the next section. Note that we only use N α values in the interval [0, 1] for α: we choose α in the set S Nα = {α j = j/N α , j = 0, 1, . . . , N α -1}.

Simulation study

To evaluate the performance of the proposed method, we generate simulated data from the following regression model:

Y = (X T β 1 ) 2 exp(X T β 1 /θ) + γ(X T β 2 ) 2 exp(X T β 2 /θ) + ǫ, (7)
where X follows a p-dimensional standardized normal distribution and ǫ is standard normally distributed. We take β 1 = (1, 1, 1, 0, . . . , 0) T and β 2 = (0, . . . , 0, 1, 1, 1) T . To visualize the data, we simulate samples of n = 100 data points from this model for different levels of θ: 1, 5 and 100, and for γ = 0. Figure 1 shows the plots of the response variable Y versus the index X T β 1 for these different values of θ. Clearly, the parameter θ has an influence on the form of the dependence between the index X T β 1 and Y . For γ = 0, when the value of θ is large (resp. small or medium), model ( 8) is a symmetric (resp. non symmetric or partially symmetric) dependent model. Therefore, the parameter θ affects the choice of α in the SIR α method. The parameter γ clearly has an influence on the choice of the dimension K: when γ = 0, only one direction is used in model ( 8), whereas if γ is nonnull (for instance if γ is fixed at 1), we have two directions.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • -4 -3 -2 -1 0 1 2 3 0 50 100 150 200 θ θ = = 1 X t β β1 Y • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • -4 -2 0 2 4 0 10 20 30 40 θ θ = = 5 X t β β1 Y • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • -4 -2 0
The performance of the proposed method will be evaluated for different sample size (n = 300 or 500), various dimensions of the explanatory variable (p = 5 or 10), several choices of θ (=1, 5

or 100) and γ (=0 or 1). The number of slices used in the slicing step is given by the expression

H = max( √ n, p).
The number B of bootstrap replications is chosen to be equal to 50. In these simulations, the parameter α varies in S Nα for N α = 11.

In the next subsection, we detail our method applied on some simulated samples. Then, in subsection 4.2, we comment a complete simulation study. In subsection 4.3, we evaluate the robustness of the approach when X does not have a multivariate normal distribution. Simulations were performed with R. All the source codes are available from the authors by e-mail.

Simulated example

We consider here several samples corresponding to specific choices for the couple of parameters (k, α).

Symmetric dependent model (θ = 100) and one direction (γ = 0). We generate a sample data of size n = 300 from the model ( 8) with p = 5. For each value of α ∈ S Nα and k = 1, . . . , 5, we compute R k,α . For simplicity, we first investigate the effect of the parameter α on R k,α for fixed values of k, then we investigate the behavior of k on R k,α for fixed values of α. Finally, the influence of the couple (k, α) is simultaneously studied.

-Influence of α on SIR α estimates. k,α and the plots of the Rk,α versus the values of α for several values of k (1, 2 and 3). Note that for α = 0, the SIR α method fails ( R k,α ≃ 0.6 for every values of k). When α > 0, Rk,α provides greater values which are close to one for k = 1, whereas the Rk,α 's values are always lower than 0.8 when k = 2 or 3. 0 0.1 0.3 0.5 0.7 0.9 1 -Choice of dimension K in the SIR α method. Figure 3 gives the boxplots of the R (b) k,α 's and the plots of Rk,α versus the values of k for several values of α (0, 0.1 and 1). Note that for α = 0.1 and α = 1, the Rk,α 's are close to 1 when k = 1, and they decrease for k = 2 before slowly increasing again.

• • • • • • • • • • • • • • • • • • 0.
• • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=1 α α R ^kα α • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=2 α α R ^kα α • • • • • • • • • • • 0.
-Choice of dimension K and α in the SIR α method. The 3D-graphic in Figure 4 exhibits the Rk,α values versus α and k. Note that the plots at the bottom of Figure 2 (resp. 3) are contained in the 3D-plot of Figure 4 if we focus only on the α-axis for a fixed k (resp. on the k-axis for a fixed α).

In view of these figures, we clearly choose the dimension K = 1 and α = 0.1. Symmetric dependent model (θ = 100) and two directions (γ = 1). First, we generate a sample data of size n = 500 from the model ( 8) with p = 5. For each value of α and k = 1, . . . , 5, we compute R k,α .

-Influence of α on SIR α estimates. Figure 5 gives the boxplots of the R (b) k,α 's and the plots of Rk,α versus the values of α for several values of k (1, 2 and 3). Note that for α = 0, the SIR α method fails ( R k,α = 0.68). For α > 0, Rk,α gives greater values which are close to one for k = 1 and 2.

When k = 3, the values of R k,α are significantly lower than the previous ones. -Choice of dimension K in the SIR α method. Figure 6 gives the boxplots of the R (b) k,α 's and the plots of Rk,α versus the values of k for several values of α (0, 0.1 and 1).Note that for α = 0.1 and α = 1, the Rk,α 's are close to 1 when k = 1 and 2, then they decrease for k = 3 and increase after.

• • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=1 α α R ^kα α b 0 0.1 0.3 0.5 0.7 0.9 1 • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=2 α α R ^kα α b 0 0.1 0.3 0.5 0.7 0.9 1 • 0.0 0.2 0.4 0.6 0.8 1.0 k=3 α α R ^kα α b 0 0.1 0.3 0.5 0.7 0.9 1 • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=1 α α R ^kα α • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=2 α α R ^kα α • • • • • • • • • • • 0.0 0.2 0.
-Choice of dimension K and α in the SIR α method. The 3-D graphic in Figure 7 shows the Rk,α 's values versus α and k. In view of Figures 5, 6 and7, we can choose K = 2 and α = 0.1. In this case, we observe a plateau in Figure 7 for α = 0.1. The explanation of this phenomenon is that the estimated eigenvectors b1 (resp. b2 ) and b(b) 1 , corresponding to the s (b) samples, (resp. b(b) 2 ) are colinear (i.e. each direction appears to be individually identifiable); then the quality of the estimated EDR space for k = 1 and k = 2 is good (close to one). Nevertheless, if we only keep the first direction, we lose information provided by the second one (which is in the true two-dimensional

• • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 0 k R ^kα α b 1 2 3 4 5 • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 0.1 k R ^kα α b 1 2 3 4 5 • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 1 k R ^kα α b 1 2 3 4 5 • • • • • 1 2 3 4 5

Results of the simulation study

For each combination of θ (=1, 5 or 100), γ (=0 or 1), p (=5 or 10) and n(=300 or 500), we generate N = 500 samples {(X i , Y i ), i = 1, . . . , n} from model ( 8). The EDR space has been estimated by the SIR α method. For each sample, we evaluate the corresponding values of R α,k for k = 1, . . . , p and α varying in S Nα . For concision, we show here only few results for several combinations of θ, γ, p and n.

Single index (γ = 0) model. Table 1 shows the means of the R α,k 's over the N = 500 replications (denoted by R k,α ) for n = 300 and p = 5. When there is no symmetric dependence (θ = 1) in the model, the values 0 or 0.1 for α provide R α,k values close to 1 only for k = 1. When θ = 5 or 100 (symmetric dependent model), the SIR-I (α = 0) fails and we need to take α > 0 in order to get a R α,k value close to 1 for k = 1. When k > 1 and α > 0, the averages of the R α,k 's values are always smaller than those obtained with k = 1. values of θ when γ = 0, n = 300 and p = 5.

. R k,α θ = 1 θ = 5 θ = 100 α k = 1 k = 2 k = 3 k = 4 k = 1 k = 2 k = 3 k = 4 k = 1 k = 2 k = 3 k = 4 0 0.
Model with two directions (γ = 1). Let us consider the N = 500 samples of size n = 500 for p = 10 and θ = 100. Figure 8 shows the boxplots of the R k,α over the N = 500 replications and the plots of the mean of the Rk,α 's versus the values of α, for several values of k (1, 2 and 3). Note that for α = 0, the SIR α method fails (with ¯ R k,α ≤ 0.5). For α > 0, the Rk,α 's provide greater values that are close to one for k = 2. Figure 9 shows the boxplots of the R k,α 's values over the N = 500

• • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=1 R ^kα α 0 0.1 0.3 0.5 0.7 0.9 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=2 R ^kα α 0 0.1 0.3 0.5 0.7 0.9 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=3 R ^kα α 0 0.1 0.3 0.5 0.7 0.9 • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=1 α α R ^kα α • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=2 α α R ^kα α • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 k=3 α α R ^kα α Figure 8: Boxplots of the R k,α 's values (above) and plots of R k,α , mean of the R k,α 's, (below)
versus α, for k = 1, 2 and 3.

replications and the plots of the mean of the Rk,α 's versus the values of k, for several values of α (0, 0.1 and 1). Note that for α = 0.1 and 1, the Rk,α 's are close to 0.8 when k = 1, that they increase to 0.9 for k = 2, and then decrease for k = 3 and 4, before increasing slowly thereafter. We also show the results when p = 5 (with the same simulation parameters: γ = 1 , θ = 100, n = 500). Figure 10(b) shows the mean of the Rk,α 's values over the N = 500 replications versus α and k. Note that for α = 0, the SIR α method fails. For α > 0, Rk,α gives greater values that are close to one for k = 1 and 2. When k = 3, the values of R k,α are significantly lower than the previous ones. In this case, we observe the same plateau phenomenon as described in the second example of section 4.1. The two directions appear to be individually identifiable, and a straightforward consequence is that the quality of the estimated EDR space for k = 1 and k = 2 is good (close to one). Nevertheless, in order to keep all the EDR space information, we need to consider the two-dimensional subspace. Note that in the previous case (see Figure 10(a)), the first two directions were not individually identifiable, contrary to the two-dimensional EDR space; therefore the quality of the criterion for k = 1 is poor while it is close to one when k = 2. Hence, the best choice in mean is the dimension K = 2 and an α > 0.1.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 0 k R ^kα α 1 2 3 4 5 6 7 8 9 10 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 0.1 k R ^kα α 1 2 3 4 5 6 7 8 9 10 • • • • • • • • • • • • • • • • • • • • • • • • • 0.

Case of a non-multivariate X

We also consider the case when X does not come from a multivariate normal distribution in order to obtain some idea of the robustness of our approach: the distribution of X is

1 2 N p (µ 2 , Σ 2 ) + 1 2 N p (µ 1 , Σ 1 )
where µ 1 = (-2, . . . , -2) T , Σ 1 = 0.5I p , µ 2 = (2, . . . , 2) T , Σ 2 = I p . We present the results for the case where n = 500, p = 10, γ = 1 and θ = 100. Figure 11 shows the boxplots close to 1 when k = 2, then they decrease until k = 4, and that they slowly increase thereafter.

• • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=1 R ^kα α 0 0.1 0.3 0.5 0.7 0.9 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=2 R ^kα α 0 0.1 0.3 0.5 0.7 0.9 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 k=3 R ^kα α 0 0.1 0.3 0.5 0.7 0.9 • • • • • • • • • • • 0.
• • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 0 k R ^kα α 1 2 3 4 5 6 7 8 9 10 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 0.1 k R ^kα α 1 2 3 4 5 6 7 8 9 10 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.0 0.2 0.4 0.6 0.8 1.0 α α = = 1 k R ^kα
Figure 13 represents the mean of the Rk,α 's values over all the replications versus α and k. In view of figures 11, 12 and 13, the best choice in mean is the dimension K = 2 and an α = 0.1.

Comparison with other methods

In this subsection, we compare the proposed bootstrap method with the nested tests procedure introduced by [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] and with the method developed by [START_REF] Ferré | Determining the dimension in Sliced Inverse Regression and related methods[END_REF], which is based on an estimation of the risk function R k thanks to asymptotic expansions. These two methods only focus on the SIR-I approach (that is when α = 0 in SIR α ). Therefore, the choice of the parameter α is not necessary and is fixed at zero in our approach. This is a first limitation of these two existing methods for determining dimension K. In the following,we recall these two methods. Then we describe the simulated model. Finally we comment the results of this simulation study.

The nested tests procedure. [START_REF] Li | Sliced inverse regression for dimension reduction, with discussion[END_REF] follows a χ 2 distribution with (p -K)(H -K -1) degrees of freedom asymptotically. From this result, if the rescaled λ(p-k) is larger than the corresponding χ 2 value (say the 95th percentile), we may infer that there are at least k + 1 (significant) indices in the model.

Ferré's method. [START_REF] Ferré | Determining the dimension in Sliced Inverse Regression and related methods[END_REF] proposed to estimate the risk function R k defined in (2). Under technical assumptions and when X follows an elliptically distribution, Ferré showed that:

R k = 1 - 1 nk   k i=1 K j=k+1 (λ i + λ j + (-1 + κ)λ i λ j ) (λ i -λ j ) 2 + (p -K) q i=1 1 λ i   + O(n -3/2 )
for k = 1 to K where κ is the usual kurtosis parameter. From this asymptotic expansion, an estimator Rk is defined by substituting the eigenvalues and κ by their estimates. Nevertheless, this criterion depends on the true dimension K (unknown). Ferré proposed to compute Rk (in fact a best notation should be here Rk,K ) for K = 1 to p and k = 1 to p, and he mentioned that "the decision rule deduces from the observation of the p curves obtained by plotting Rk,K versus k, for k = 1 to p and for K = 1 to p". In practice, it is not really easy to determine the dimension from this graphical aspect (in particular when p is large). Final comments. Note that our bootstrap approach always provides a good decision concerning the choice of dimension, contrary to the other ones. Moreover, this approach also makes it possible the choice of α (which was fixed, α = 0, in the previous comparison) for all kinds of models (especially for symmetric dependent models).

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
The nested tests procedure and the Ferré approach rely on multinormal or elliptical assumption for the distribution of X. In the previous subsection, we show that our method still provides good performance in non elliptical case.

Concluding remarks

This article proposes a practical criterion based on bootstrap to choose the dimension and α in the SIR α method. The simulation study demonstrates good numerical behaviour of this approach, contrary to existing methods (which only focus on the choice of dimension K). We underline the pleasant graphical aspect of the criterion which allows the practitioner to choose both K and α.

YFigure 1 :

 1 Figure 1: Plots of Y versus the index X T β 1 for different values of θ with γ = 0.

  Figure 2 gives the boxplots of the R (b)

  Figure 2: Boxplots of the R (b) k,α 's values (above) and plots of Rk,α (below) versus α, for k = 1, 2 and 3.

Figure 3 :

 3 Figure 3: Boxplots of the R (b) k,α 's values (above) and plots of Rk,α (below) versus k, for α = 0, 0.1 and 1.

Figure 4 :

 4 Figure 4: Plot of Rk,α versus α and k.

Figure 5 :

 5 Figure 5: Boxplots of the R (b) k,α 's values (above) and plots of Rk,α (below) versus α, for k = 1, 2 and 3.

Figure 6 :

 6 Figure 6: Boxplots of the R (b) k,α 's values (above) and plots of Rk,α (below) versus k for α = 0, 0.1 and 1.

Figure 7 :

 7 Figure 7: Plot of Rk,α versus α and k.

Figure 10

 10 Figure 10(a) represents the mean of the Rk,α 's values over the N = 500 replications versus α and k. In view of Figures 8, 9 and 10(a), the best choice (in mean) is the dimension K = 2 and an α > 0.1.

Figure 9 :Figure 10 :

 910 Figure 9: Boxplots of the R k,α 's (above) and plots of R k,α (below) versus k for α = 0, 0.1 and 1.

Figure 11 :

 11 Figure 11: Boxplots of the R k,α 's values (above) and plots of R k,α (below) versus α, for k = 1, 2 and 3.

Figure 12 :

 12 Figure 12: Boxplots of the R k,α 's values (above) and plots of R k,α (below) versus k for α = 0, 0.1 and 1.

Figure 13 :

 13 Figure 13: Plot of R k,α versus α and k.

Figure 15 :

 15 Figure 15: Boxplots of the p-values of the nested tests for different values of (n, σ 2 )

Table 1 :

 1 Means of the R k,α 's values over the N = 500 replications, denoted by R k,α , for different

		989 0.706 0.743 0.843 0.922 0.692 0.738 0.843 0.393 0.554 0.694 0.831
	0.1	0.990 0.710 0.748 0.851 0.972 0.697 0.746 0.850 0.971 0.697 0.742 0.850
	0.2	0.989 0.713 0.750 0.851 0.965 0.703 0.749 0.849 0.971 0.697 0.743 0.850
	0.3	0.977 0.711 0.752 0.851 0.939 0.697 0.746 0.851 0.968 0.700 0.743 0.849
	0.4	0.890 0.706 0.749 0.852 0.898 0.691 0.744 0.850 0.967 0.700 0.742 0.849
	0.5	0.548 0.669 0.745 0.849 0.844 0.683 0.742 0.848 0.966 0.697 0.743 0.849
	0.6	0.250 0.537 0.717 0.847 0.798 0.674 0.740 0.848 0.964 0.698 0.744 0.849
	0.7	0.355 0.494 0.659 0.830 0.763 0.666 0.738 0.848 0.964 0.696 0.743 0.850
	0.8	0.409 0.580 0.703 0.825 0.731 0.656 0.736 0.847 0.963 0.697 0.744 0.849
	0.9	0.419 0.608 0.763 0.903 0.714 0.651 0.734 0.847 0.962 0.698 0.742 0.848
	1	0.420 0.613 0.780 0.957 0.698 0.642 0.730 0.846 0.962 0.696 0.742 0.848

Simulated model. We generate simulated data from the following regression model (in which the true dimension K is equal to 2):

where X follows a 10-dimensional standardized normal distribution and ǫ is normally distributed with σ 2 . We take β 1 = (1, 1, 0, . . . , 0) ′ and β 2 = (0, 0, 0, 1, 1, 0, . . . , 0) ′ . Note that we consider here the multinormal case corresponding with the theoretical assumption of the nested tests procedure and Ferré method. Moreover, this model is not a symmetric dependent model, then SIR-I should provide good performance.

The performances of the different methods are evaluated for different sample sizes (n = 100, 400

or 1000) and several choices of σ 2 (= 0.01 or 1). For each value of (n, σ 2 ), N = 500 samples have been generated. For each sample, we determined the dimension with the three methods (nested tests, Ferré's approach, bootstrap approach).

Results and comments. Using our bootstrap approach, Figure 14 shows the boxplots of the Rk,α 's values over the N = 500 replications and the plots of the mean of the Rk,α 's versus the values of k, in the case where σ 2 = 1, n = 100 and n = 400. Note that for n = 100, the Rk,α 's are close to 0.7 when k = 1, increase to 0.8 for k = 2, and then decrease for k = 3 and 4, before increasing slowly thereafter. In view of this figure, the best choice (in mean) is the dimension K = 2. Similar (and best) results have been observed for n = 400 and σ 2 = 1 (see Figure 14). When σ 2 = 0.01 and n = 100, 400 or 1000, we obtained similar graphics which are not shown in the paper.

Let us now comment the results from the nested tests procedure. Figure 15 shows the boxplots of the p-values of the nested tests for all values of (n, σ 2 ). On these graphics, if a p-value is lower than 5% for a value of k and is greater than 5% for k + 1, the dimension k + 1 is chosen. Note that the only case where this procedure had good performance (that is determine the true dimension) is for large sample size (n = 1000) and small variance (σ 2 = 0.01). The dimension K = 2 has been always chosen. In the other cases, the method failed. For instance, for medium sample size (n = 400) and small variance (σ 2 = 0.01), the procedure chose K = 2 for 23 simulated samples, and K = 1 for 477 samples. When the variance is σ 2 = 1, the procedure never find the true dimension ( K = 1 for n = 1000, no dimension reduction for n = 100).

Concerning Ferré approach, the estimate Rk,K of R K only gave reasonable values (in [0,1]) for K = 1 and for the true dimension K = 2. In the other cases, the values for Rk,K are always (in mean over the N replications) negative values. This phenomenon is particularly obvious for small sample size and large variance. This problem has been already mentioned by [START_REF] Ferré | Un critere de choix de la dimension dans la méthode SIR-II[END_REF] for a SIR-II version of this approach. An explanation of these "bad" values can be the proximity of the eigenvalues. From these remarks, it seems difficult to use this method to determine dimension, even if Rk,K appears to be well estimated for the true dimension (which is unknown in practice).

The method has been implemented in R, and the source codes are available from the author. They provide the three kinds of graphics in order to visualize the quality of the estimated EDR space when k (resp. α) varies for a fixed α (resp. k), and they give the 3D-plot of this quality when the couple (k, α) varies on a grid. We used a rough grid in our simulation study, but for a dataset, we recommend working with a thin grid.

Importantly the choice of parameter α seems to be less sensitive than the choice of dimension K. More precisely, for the choice of α, there are generally two scenarios: the case where α = 0 and the case where α > 0 (when a symmetric dependence occurs in the model). The 3D-graphics in Section 4 show that for each fixed k, the profiles of the criterion are nearly similar when α > 0.

Finally, this approach can be extended to the multivariate SIR α method named pooled marginal slicing PMS α , see [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF] or [START_REF] Barreda | Some extensions of multivariate sliced inverse regression[END_REF].