Alessandro Chiancone
email: al.chiancone@gmail.com

Florence Forbes

Stéphane Girard

Student Sliced Inverse Regression

Keywords: Dimension reduction, Inverse regression, Outliers, Robust estimation, Generalized Student distribution

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.

Introduction

Let us consider a regression setting where the goal is to estimate the relationship between a univariate response variable Y and a predictor X. When the dimension p of the predictor space is 1 or 2, a simple 2D or 3D plot can visually reveal the relationship and can be useful to determine the regression strategy to be used. If p becomes large such an approach is not feasible. A possibility to overcome problems arising in the context of regression is to make the assumption that the response variable does not depend on the whole predictor space but just on a projection of X onto a subspace of smaller dimension. Such a dimensionality reduction leads to the concept of sufficient dimension reduction and to that of central subspace [START_REF] Cook | Graphics for regressions with a binary response[END_REF]. The central subspace is the intersection of all dimension-reduction subspaces (d.r.s.). A subspace S is a d.r.s. if Y is independent of X given P S X, where P S is the orthogonal projection onto S. In other words, all the information carried by the predictors X on Y can be compressed in P S X. It has been shown under weak assumptions that the intersection of all d.r.s., and therefore the central subspace, is itself a d.r.s. [START_REF] Yin | Successive direction extraction for estimating the central subspace in a multiple-index regression[END_REF]. It is of particular interest to develop methods to estimate the central subspace as once it is identified, the regression problem can be solved equivalently using the lower-dimensional representation P S X of X in the subspace.

Among methods that lead to an estimation of the central subspace, Sliced Inverse Regression (SIR) [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] is one of the most popular. SIR is a semiparametric method assuming that the link function depends on d linear combinations of the predictors and a random error independent of X: Y = f (β T 1 X, . . . , β T d X, ). When this model holds, the projection of X onto the space spanned by the vectors {β i , i = 1, . . . , d} captures all the information about Y . In addition, [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] shows that a basis of this space can be recovered using an inverse regression strategy provided that the so called linearity condition holds. It has been shown that the linearity condition is satisfied as soon as X is elliptically distributed. Moreover, this condition approximately holds in high-dimensional datasets, see [START_REF] Hall | On almost linearity of low dimensional projections from high dimensional data[END_REF]. However, solutions have been proposed to deal with non elliptical distributed predictors and to overcome the linearity condition limitation [START_REF] Fukumizu | Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces[END_REF][START_REF] Li | Dimension reduction for nonelliptically distributed predictors[END_REF][START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF].

The inverse regression approach to dimensionality reduction gained then rapid attention [START_REF] Cook | Sliced inverse regression for dimension reduction: Comment[END_REF] and was generalized in [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] which shows the link between the axes spanning the central subspace and an inverse regression problem with Gaussian distributed errors. More specifically, in [START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF][START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF], it appears that, for a Gaussian error term and under appropriate conditions, the SIR estimator can be recovered as the maximum likelihood estimator of the parameters of an inverse regression model. In other words, although SIR is originally a model free method, the standard SIR estimates are shown to correspond to maximum likelihood estimators for a Gaussian inverse regression model. It is then not surprising that SIR has been observed, e.g. in [START_REF] Bura | Extending sliced inverse regression: The weighted chi-squared test[END_REF], to be at best under normality and that its performance may degrade otherwise. Indeed, the Gaussian distribution is known to have tails too light to properly accommodate extreme values. In particular, [START_REF] Sheather | A comparison of procedures based on inverse regression[END_REF] observes that SIR was highly sensitive to outliers, with additional studies, evidence and analysis given in [START_REF] Gather | A note on outlier sensitivity of sliced inverse regression[END_REF]. To downweight this sensitivity, robust versions of SIR have been proposed, mainly starting from the standard model free estimators and trying to make them more resistant to outliers. Typically, in [START_REF] Gather | A robustified version of sliced inverse regression[END_REF] classical estimators are replaced by high breakdown robust estimators and, recently in [START_REF] Dong | Robust inverse regression for dimension reduction[END_REF] two approaches are built: a weighted version of SIR and a solution based on the intra slice multivariate median estimator.

As an alternative, we propose to rather exploit the inverse regression formulation of SIR [START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF][START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF]. A new error term modeled by a multivariate Student distribution [START_REF] Kotz | Multivariate t Distributions and their Applications[END_REF] is introduced. Among the elliptically contoured distributions, the multivariate Student is a natural generalization of the multivariate Gaussian but its heavy tails can better accommodate outliers. The result in Proposition 6 of [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] is extended from Gaussian to Student errors showing that the inverse regression approach of SIR is still valid outside the Gaussian case, meaning that the central subspace can still be estimated by maximum likelihood estimation of the inverse regression parameters. It is then shown that the computation of the maximum likelihood estimators remains tractable in the Student case via an Expectation-Maximization (EM) algorithm which has a simple implementation and desirable properties.

The paper is organized as follows. In Section 2 general properties of the multivariate Student distribution and some of its variants are first recalled. The inverse regression model is introduced in Section 3 followed by the EM strategy to find the maximum likelihood estimator, the link with SIR and the resulting Student SIR algorithm. A simulation study is carried out in Section 4 and a real data application, showing the interest of this technique, is detailed in Section 5. The final section contains concluding remarks and perspectives. Proofs are postponed to the Appendix.

Multivariate generalized Student distributions

Multivariate Student, also called t-distributions, are useful when dealing with real-data because of their heavy tails. They are a robust alternative to the Gaussian distribution, which is known to be very sensitive to outliers. In contrast to the Gaussian case though, no closed-form solution exists for the maximum likelihood estimation of the parameters of the t-distribution. Tractability is, however, maintained both in the univariate and multivariate case, via the EM algorithm [START_REF] Mclachlan | Robust mixture modelling using the t distribution[END_REF] and thanks to a useful representation of the t-distribution as a so-called infinite mixture of scaled Gaussians or Gaussian scale mixture [START_REF] Andrews | Scale mixtures of normal distributions[END_REF]. A Gaussian scale mixture distribution has a probability density function of the form

P (x; µ, Σ, ψ) = ∞ 0 N p (x; µ, Σ/u) f U (u; ψ) du, (1)
where N p ( . ; µ, Σ/u) denotes the density function of the p-dimensional Gaussian distribution with mean µ and covariance Σ/u and f U is the probability distribution of a univariate positive variable U referred to hereafter as the weight variable. When f U is a Gamma distribution G(ν/2, ν/2) where ν denotes the degrees of freedom, expression (1) leads to the standard pdimensional t-distribution denoted by t p (x; µ, Σ, ν) with parameters µ (lo-cation vector), Σ (p×p positive definite scale matrix) and ν (positive degrees of freedom parameter). Its density is given by

t p (x; µ, Σ, ν) = ∞ 0 N p (x; µ, Σ/u) G(u; ν/2, ν/2) du = Γ((ν + p)/2) |Σ| 1/2 Γ(ν/2) (πν) p/2 [1 + δ(x, µ, Σ)/ν] -(ν+p)/2 , (2)
where δ(x, µ, Σ) = (x -µ) T Σ -1 (x -µ) is the Mahalanobis distance between x and µ. The Gamma distribution has probability density function G(u; α, γ) = u α-1 Γ(α) -1 exp(-γu)γ α , where Γ denotes the Gamma function.

If f U (u; ψ) is set equal to a Gamma distribution G(α, γ) without imposing α = γ, (1) results in a multivariate Pearson type VII distribution (see e.g. [START_REF] Johnson | Continuous Univariate Distributions[END_REF] vol.2 chap. 28) also referred to as the Arellano-Valle and Bolfarine's Generalized t distribution in [START_REF] Kotz | Multivariate t Distributions and their Applications[END_REF]. This generalized version is the multivariate version of the t-distribution considered in this work, its density is given by:

S p (x; µ, Σ, α, γ) = ∞ 0 N p (x; µ, Σ/u) G(u; α, γ) du (3) = Γ(α + p/2) |Σ| 1/2 Γ(α) (2πγ) p/2 [1 + δ(x, µ, Σ)/(2γ)] -(α+p/2) . (4)
For a random variable X following distribution (4), an equivalent representation useful for simulation is X = µ + U -1/2 X where U follows a G(α, γ) distribution and X follows a N (0, Σ) distribution.

Remark 1 (Identifiability). The expression (4) depends on γ and Σ only through the product γΣ which means that to make the parameterization unique, an additional constraint is required. One possibility is to impose that Σ is of determinant 1. It is easy to see that this is equivalent to have an unconstrained Σ with γ = 1.

Unconstrained parameters are easier to deal with in inference algorithms. Therefore, we will rather assume without loss of generality that γ = 1 with the notation S p (0, V, α, 1) ≡ S p (0, V, α) adopted in the next Section.

Student Sliced Inverse Regression

Let X ∈ R p be a random vector, Y ∈ R the real response variable and S Y |X the central subspace spanned by the columns of the matrix β ∈ R p×d . In the following, it is assumed that dim(S Y |X ) = d where d is known and d ≤ p. To address the estimation of the central subspace, we consider the inverse regression formulation of [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF], which models the link from Y to X. In addition to be a simpler regression problem, the inverse regression approach is of great interest because Proposition 6 in [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] states that in the Gaussian case, an estimation of the central subspace is provided by the estimation of the inverse regression parameters. In Subsection 3.1, the inverse regression model of [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] is extended by considering Student distributed errors. It is then shown in Subsection 3.2 that the estimation of the extended model is tractable via an Expectation-Maximization algorithm (EM). A link with SIR is presented in Subsection 3.3 and the resulting Student SIR algorithm is described in Subsection 3.4.

Student multi-index inverse regression model

In the spirit of [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF][START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF] the following regression model is considered

X = µ + VBc(Y ) + ε, (5)
where µ ∈ R p is a non random vector, B is a non random p × d matrix with B T B = I d , ε ∈ R p is a centered generalized Student random vector following the distribution given in (4), ε is assumed independent of Y , with scale matrix V, c : R → R d is a non random function. It directly follows from ( 5) that

E(X|Y = y) = µ + VBc(y), (6)
and thus, after translation by µ, the conditional expectation of X given Y is a random vector located in the space spanned by the columns of VB.

When ε is assumed to be Gaussian distributed, Proposition 6 in [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] states that B corresponds to the directions of the central subspace β. In [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF][START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF], it appears then that, under appropriate conditions, the maximum likelihood estimator of B is (up to a full rank linear transformation) the SIR estimator of β, i.e. Span{B} = Span{β}. Proposition 6 in [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] can be generalized to our Student setting, so that B still corresponds to the central subspace. The generalization of Proposition 6 of [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] is given below.

Proposition 1. Let X y be a random variable distributed as X|Y = y, let us assume that

X y = µ + VBc(y) + ε, (7)
with ε following a generalized Student distribution S p (0, V, α), c(y) ∈ R d is function of y and VB is a p × d matrix of rank d. Under model [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF], the distribution of Y |X = x is the same as the distribution of Y |B T X = B T x for all values x.

The proof is given in Appendix 7.1. According to this proposition, X can be replaced by B T X without loss of information on the regression of Y on X.

A procedure to estimate B is then proposed in the next Section

Maximum likelihood estimation via EM algorithm

Let (X i , Y i ), i = 1, . . . , n be a set of independent random variables distributed according to the distribution of (X, Y ) as defined in [START_REF] Fukumizu | Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces[END_REF]. The unknown quantities to be estimated in model [START_REF] Fukumizu | Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces[END_REF] are {µ, V, B, α} and the function c(.). Regarding c, we focus on projection estimators for each coordinate of c(.) = (c 1 (.), . . . , c d (.)). For k = 1, . . . , d, function c k (.) is expanded as a linear combination of h basis functions s j (.), j = 1, . . . , h as

c k (.) = h j=1 c jk s j (.), (8)
where the coefficients c jk , j = 1, . . . , h and k = 1, . . . , d are unknown and to be estimated while h is supposed to be known. Let C be a h × d matrix with the kth column given by (c 1k , . . . , c hk ) T and s(.) = (s 1 (.), . . . , s h (.)) T . Then, model ( 5) can be rewritten as

X = µ + VBC T s(Y ) + ε, with ε ∼ S p (0, V, α), (9)
where S p (0, V, α) is the multivariate centered generalized Student distribution with scale matrix V. For each i, it follows that conditionally to Y i , X i ∼ S p (µ + VBC T s i , V, α) where s i =s(Y i ). The density of the generalized Student distribution is available in closed form and given in (4). However to perform the estimation, a more useful representation of this distribution is given by its Gaussian scale mixture representation (3). Introducing an additional set of latent variables U = {U 1 , . . . , U n } with U i independent of Y i , one can equivalently write:

X i |U i = u i , Y i = y i ∼ N p (µ + VBC T s i , V/u i ), (10)
U i |Y i = y i ∼ G(α, 1). ( 11
)
Let us denote by θ = {µ, V, B, C, α} the parameters to estimate from realizations {x i , y i , i = 1, . . . , n}. In contrast to the Gaussian case, the maximum likelihood estimates are not available in closed-form for the tdistributions. However, they are reachable using an Expectation-Maximization (EM) algorithm. More specifically, at iteration (t) of the algorithm, θ is updated from a current value θ (t-1) to a new value θ (t) defined as θ (t) = arg max θ Q(θ, θ (t-1) ). Considering the scale mixture representation above, a natural choice for Q is the following expected value of the complete loglikelihood:

Q(θ, θ (t-1) ) = E U [ n i=1 log P (x i , U i |Y i = y i ; θ)|X i = x i , Y i = y i ; θ (t-1) ] (12) = n i=1 E Ui [log P (x i |U i , y i ; µ, V, B, C)|x i , y i ; θ (t-1) ] + E Ui [log P (U i ; α)|x i , y i ; θ (t-1) ] = - 1 2 n log det V + 1 2 p n i=1 E Ui [log(U i )|x i , y i ; θ (t-1) ] - 1 2 n i=1 E Ui [U i |x i , y i ; θ (t-1) ] (µ + VBC T s i -x i ) T V -1 (µ + VBC T s i -x i ) + n i=1 E Ui [log P (U i ; α)|x i , y i ; θ (t-1) ] .
Note that all computations are conditionally to the Y i 's and no assumption is made on the distribution of the Y i 's. The E-step therefore consists of computing the quantities

ūi (t) = E U i [U i |x i , y i ; θ (t-1) ] , (13)
ũi (t) = E U i [log U i |x i , y i ; θ (t-1) ] , (14)
while the M-step divides into two-independent M-steps involving separately parameters (µ, V, B, C) and α. The second quantity ( 14) is needed only in the estimation of α. The following notation is introduced for the next sections:

ū(t) = n i=1 ūi (t) n (15) ũ(t) = n i=1 ũi (t) n . (16)
E-step. The quantities ( 13) and ( 14) above require the posterior distribution of the U i 's. This distribution can be easily determined using the well known conjugacy of the Gamma and Gaussian distributions for the mean. It follows then from standard Bayesian computations that the posterior distribution is still a Gamma distribution with parameters specified below,

p(u i |X i = x i , Y i = y i ; θ (t-1) ) ∝ N p (x i ; µ (t-1) + V (t-1) B (t-1) C (t-1)T s i , V (t-1) /u i ) G(u i ; α (t-1) , 1) = G(u i ; α (t-1) + p 2 , 1 + 1 2 δ(x i , µ (t-1) + V (t-1) B (t-1) C (t-1)T s i , V (t-1) )), where δ(x i , µ+VBC T s i , V) = (µ+VBC T s i -x i ) T V -1 (µ+VBC T s i -x i ) is the Mahalanobis distance between x i and µ + VBC T s i when the covariance is V.
The required moments ( 13) and ( 14) are then well known for a Gamma distribution, so that it comes,

ū(t) i = α (t-1) + p 2 1 + 1 2 δ(x i , µ (t-1) + V (t-1) B (t-1) C (t-1)T s i , V (t-1) ) and ũ(t) i = Ψ(α (t-1) + p 2 ) -log(1 + 1 2 δ(x i , µ (t-1) + V (t-1) B (t-1) C (t-1)T s i , V (t-1) )) ,
where Ψ is the Digamma function. As it will become clear in the following M-step, ū(t) i acts as a weight for x i . Whenever the Mahalanobis distance of x i to µ (t-1) + V (t-1) B (t-1) C (t-1)T s i increases, the weight ū(t) i of x i decreases and the influence of x i in the estimation of the parameters will be downweighted in the next iteration. The idea of using weights to handle outliers is common in the literature, Weighted Inverse Regression (WIRE) [START_REF] Dong | Robust inverse regression for dimension reduction[END_REF] gives weights through a deterministic kernel function to ensure the existence of the first moment. Our approach does not require previous knowledge to select an appropriate kernel and refers to the wide range of t-distributions (the Cauchy distribution for which the first moment is not defined lies in this family). M-step. The M-step divides into the following two independent sub-steps. M-(µ, V, B, C) substep. Omitting terms that do not depend on the parameters in [START_REF] Sheather | A comparison of procedures based on inverse regression[END_REF], estimating (µ, V, B, C) by maximization of Q consists, at iteration (t), of minimizing with respect to (µ, V, B, C) the following G function,

G(µ, V, B, C) = log det V+ 1 n n i=1 ūi (t) (µ+VBC T s i -x i ) T V -1 (µ+VBC T s i -x i ) . ( 17
)
To this aim, let us introduce (omitting the index iteration (t) in the notation) the h × h weighted covariance matrix W of s(Y ) defined by:

W = 1 n n i=1 ūi (s i -s)(s i -s) T ,
the h × p weighted covariance matrix M of (s, X) defined by

M = 1 n n i=1 ūi (s i -s)(x i -x) T ,
and Σ the p × p weighted covariance matrix of

X Σ = 1 n n i=1 ūi (x i -x)(x i -x) T , (18)
where

x = 1 n i=1 ūi n i=1 ūi x i and (19)
s = 1 n i=1 ūi n i=1 ūi s i . ( 20
)
We derive then the following lemma.

Lemma 1. Using the above notations, G(µ, V, B, C) can be rewritten as

G(µ, V, B, C) = log det V + tr(ΣV -1 ) + tr(C T WCB T VB) -2tr(C T MB) + ū (µ -x + VBC T s) T V -1 (µ -x + VBC T s) .
The proof is given in Appendix 7.2. Thanks to this representation of G(.) it is possible to derive the following proposition which is a generalization to the multi-index case and Student setting of the result obtained in case of Gaussian error in [START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF].

Proposition 2. Under (9), if W and Σ are regular, then the M-step for (µ, V, B, C) leads to the updated estimations ( μ, V, B, Ĉ) given below

• B is made of the eigenvectors associated to the largest eigenvalues of

Σ -1 M T W -1 M, • V = Σ -(M T W -1 MB)(B T M T W -1 MB) -1 (M T W -1 MB) T , • Ĉ = W -1 M B( BT V B) -1 and • μ = x -V B ĈT s.
The proof is detailed in Appendix 7.3. Regarding parameter α it can be updated using an independent part of Q as detailed in the next M-step.

M-α substep.

Parameter α can be estimated by maximizing independently with regards to α,

n i=1 E U i [log P (U i ; α)|x i , y i ; θ (t-1) ] . (21)
Then, since

E U i [log p(U i ; α)|x i , y i ; θ (t-1) ] = -ūi (t) + (α -1) ũi (t) -log Γ(α) , (22)
setting the derivative with respect to α to zero, we obtain that α = Ψ -1 (ũ), where Ψ(.) is the Digamma function.

In practice, for the procedure to be complete, the choice of the h basis functions s j needs to be specified. Many possibilities for basis functions are available in the literature such as classical Fourier series, polynomials, etc. In the next section, we discuss a choice of basis functions which provides the connection with Sliced Inverse Regression (SIR) [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF].

Connection to Sliced Inverse Regression

As in the Gaussian case [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF][START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF], a clear connection with SIR can be established for a specific choice of the h basis functions. When Y is univariate a natural approach is to first partition the range of Y into h + 1 bins S j for j = 1, . . . , h + 1 also referred to as slices, and then defining h basis functions by considering the first h slices as follows,

s j (.) = 1I{. ∈ S j }, j = 1, . . . , h, (23)
where 1I is the indicator function. Note that it is important to remove one of the slices so that the basis functions remain independent. However, the following related quantities are defined for j = 1, . . . , h + 1:

n j = n i=1 ūi 1I{y i ∈ S j }, f j = n j n . ( 24
)
They represent respectively the number of y i in slice j weighted by the ūi and the weighted proportion in slice j. The following weighted mean of X given Y ∈ S j is then denoted by

xj = 1 n j n i=1 ūi 1I{y i ∈ S j }x i , (25)
and the p × p "between slices" covariance matrix by

Γ = h+1 j=1 f j (x j -x)(x j -x) T .
In this context, the following consequence of Proposition 2 can be established.

Corollary 1. Under ( 9) and ( 23), if Σ is regular, then the updated estimation B of B is given by the eigenvectors associated to the largest eigenvalues of Σ -1 Γ. In addition,

Γ = M T W -1 M.
The proof is given in Appendix 7.4. When all ūi = 1, the iterative EM algorithm reduces to one M-step and the quantities defined in this section correspond to the standard SIR estimators. The EM algorithm resulting from this choice of basis functions is referred to as the Student SIR algorithm. It is outlined in the next section.

Central subspace estimation via Student SIR algorithm

The EM algorithm can be outlined using Proposition 2 and Corollary 1. It relies on two additional features to be specified, initialization and stopping rule. As the algorithm alternates the E and M steps, it is equivalent to start with one of this step. It is convenient to start with the Maximization step since the initialization of quantities ūi , ũi can be better interpreted. If ūi is constant and ũi = 0, the first M-step of the algorithm results in performing standard SIR. Regarding an appropriate stopping rule of the algorithm, EM's fundamental property is to increase the log-likelihood at each iteration. A standard criteria is then the relative increase in log-likelihood, denoted by ∆(θ (t) , θ (t-1) ), between two iterations. At each iteration, for current parameters values, the log-likelihood is easy to compute using (4) and [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF]. Another natural criterion is to assess when parameter estimation stabilizes. Typically, focusing on the central subspace B, the following proximity measure [START_REF] Ferré | Determining the dimension in sliced inverse regression and related methods[END_REF][START_REF] Chavent | A sliced inverse regression approach for data stream[END_REF] can be considered:

r(B, B) = trace(BB T B BT ) d . (26)
The above quantity r ranges from 0 to 1 and evaluates the distance between the subspaces spanned by the columns of B and B. If d = 1, r is the squared cosine between the two spanning vectors. Although not directly related to the EM algorithm, in practice this criterion gave similar results in terms of parameter estimation. Experiments on simulated and real data are reported in the next two sections.

Determination of the central subspace dimension

Determining the dimension d of the central subspace is an important issue for which different solutions have been proposed in the literature. Most users rely on graphical considerations, e.g. [START_REF] Liquet | A graphical tool for selecting the number of slices and the dimension of the model in SIR and SAVE approaches[END_REF]. A more quantitative approach is to use cross validation after the link function is found. Although in that case, d may vary depending on the specific regression approach that the user selected. Other methods that can be easily used on real data, are mainly based on (sequential) tests [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF][START_REF] Schott | Determining the Dimensionality in Sliced Inverse Regression[END_REF][START_REF] Velilla | Assessing the number of linear components in a general regression problem[END_REF][START_REF] Ferré | Determining the dimension in sliced inverse regression and related methods[END_REF][START_REF] Barrios | A bootstrap method for assessing the dimension of a general regression problem[END_REF][START_REF] Bura | Extending sliced inverse regression: The weighted chi-squared test[END_REF]]. An alternative that uses a penalized likelihood criterion has been proposed in [START_REF] Zhu | On Sliced Inverse Regression With High-Dimensional Covariates[END_REF]. In our setting, formulated as a maximum likelihood problem, the penalized likelihood approach is the most natural. For a given value d of the central subspace dimension, we therefore propose to compute the Bayesian information criterion [START_REF] Kass | Bayes factors[END_REF] defined as

BIC(d) = -2L(d) + η log n , where η = p(p+3) 2 + 1 + d(2p-d-1+2h) 2
is the number of free parameters in the model and L(d) is the maximized loglikelihood computed at the parameters values obtained at convergence of the EM algorithm. Computing L(d) is a straightforward byproduct of the algorithm described above as this quantity is already used in our stopping

Algorithm 1 Student SIR algorithm

Set h and partition the Y range into h + 1 slices. Set the e.d.r. space dimension d and the desired tolerance value for convergence δ. Initialize the ū(0) i , ũ(0) i 's with ū(0) i = 1 and ũ(0) i = 0 for all i = 1, . . . , n (this first iteration of the algorithm gives the SIR estimation of Γ and B). while ∆(θ (t) , θ (t-1) ) < δ do M-step Compute:

• ū(t) and ũ(t) (eq. ( 15) and ( 16)),

f (t) = (f (t) 1 , . . . f (t) h ) T and f (t)
h+1 (eq. ( 24)),

• x(t) j and x(t) (eq. ( 25) and ( 19)),

• Σ (t) (eq. ( 18)),

• M (t) where each row is given by M

(t) j,. = f (t) j (x (t) j -x(t) ) T for j = 1, . . . , h, • W (t)-1 = diag 1 f (t) 1 , . . . , 1 f (t) h + 1 f (t) h+1 O, where O is the h × h matrix defined by O ij = 1, • Γ (t) = M (t)T W (t)-1 M (t) ,
• B (t) matrix of the d eigenvectors associated to the d largest eigenvalues of

Σ (t)-1 Γ (t) , • V (t) = Σ (t) -Γ (t) B (t) (B (t) T Γ (t) B (t) ) -1 (Γ (t) B (t) ) T , • C (t) = W (t)-1 M (t) B (t) (B (t)T V (t) B (t) ) -1 , • µ (t) = x(t) -V (t) B (t) C (t)T s(t) , • α (t) = Ψ -1 (ũ (t)
).

E-step

Update the ūi , ũi 's using the quantities estimated in the M-step:

ūi (t+1) = α (t) + p 2 1 + 1 2 δ(x i , µ (t) + V (t) B (t) C (t)T s i , V (t) ) , ũi (t+1) = Ψ(α (t) + p 2 ) -log(1 + 1 2 δ(x i , µ (t) + V (t) B (t) C (t)T s i , V (t) )) .
end while criterion. Following the BIC principle, an estimator of d can then be defined as the minimizer of BIC(d) over d ∈ {1, . . . , min(p, h)}. The performance of this criterion is investigated in the simulation study in Section 4 and used on the real data example of Section 5. The simulation study reveals that BIC can provide correct selections but requires large enough sample sizes. This limitation has been already pointed out in the literature (see e.g. [START_REF] Giraud | Introduction to high-dimensional statistics[END_REF]).

Simulation study

Student SIR is tested on simulated data under a variety of different models and distributions for the p-dimensional random variable X. The behavior of Student SIR is compared to SIR and four other techniques arising from the literature that claim some robustness. For comparison, the simulation setup described in [START_REF] Dong | Robust inverse regression for dimension reduction[END_REF][START_REF] Gather | A robustified version of sliced inverse regression[END_REF] is adopted.

Simulation setup

Three different regression models are considered:

I : Y = 1 + 0.6X 1 -0.4X 2 + 0.8X 3 + 0.2ε, II : Y = (1 + 0.1ε)X 1 , III : Y = X 1 /(0.5 + (X 2 + 1.5) 2 )) + 0.2ε,
where ε follows a standard normal distribution. The three models are combined with three possible distributions for the predictors X:

(i) X is multivariate normal distributed with mean vector 0 and covariance matrix defined by its entries as σ ij = 0.5 |i-j| ;

(ii) X is standard multivariate Cauchy distributed;

(iii) X = (X 1 , . . . , X p ) T , where each X i is generated independently from a mixture of normal and uniform distributions denoted by 0.8N (0, 1) + 0.2U(-ν, ν) where ν is a positive scalar value.

Models I, III are homoscedastic while model II is heteroscedastic. Case (ii) is built to test the sensitivity to outliers while the distribution of X is elliptical. In (iii) a non-elliptical distribution of X is considered. The dimension is set to p = 10, the dimension of the e.d.r. space is d = 1 for I, II and d = 2 for III. The nine different configurations of X and Y are simulated with a number of samples varying depending on the experiment. In all tables Student SIR is compared with standard SIR and four other approaches. Contour Projection (CP-SIR) [START_REF] Wang | Improving dimension reduction via contourprojection[END_REF][START_REF] Luo | Contour projected dimension reduction[END_REF] applies the SIR procedure on a rescaled version of the predictors. Weighted Canonical Correlation (WCAN) [START_REF] Zhou | Robust dimension reduction based on canonical correlation[END_REF] uses a basis of B-splines first estimating the dimension d of the central subspace and then the directions from the nonzero robustified version of the correlation matrices between the predictors and the B-splines basis functions. The idea of Weighted Inverse Regression (WIRE) [START_REF] Dong | Robust inverse regression for dimension reduction[END_REF] is to use a different weight function capable of dealing with both outliers and inliers. SIR is a particular case of WIRE with constant weighting function. Slice Inverse Median Estimation (SIME) replaces the intra slice mean estimator with the median which is well known to be more robust. All values referring to CP-SIR, WCAN, WIRE, SIME in the tables are directly extracted from [START_REF] Dong | Robust inverse regression for dimension reduction[END_REF]. Values relative to SIR have been recomputed using [START_REF] Cook | Ldr: a package for likelihoodbased sufficient dimension reduction[END_REF].

Results

To assess the sensitivity of the compared methods to different setting parameters, four sets of tests are carried out and reported respectively in Tables 1 and2. First, the 9 configurations of X and Y models are tested for fixed sample size n = 200, number of slices h = 5 and p = 10 (Table 1 (a)). Then, the effect of the sample size is illustrated for model I (Table 1

(b)).

The number of slices is varied to evaluate the sensitivity to the h value (Table reftb3 (a)) and at last, different values of ν are tested in the model (iii) case (Table 2 (b)). In all cases and tables, the different methods performance is assessed based on their ability to recover the central subspace which is measured via the value of the proximity measure r [START_REF] Zhu | On Sliced Inverse Regression With High-Dimensional Covariates[END_REF].

Student SIR shows its capability to deal with different configurations.

The proximity criterion [START_REF] Zhu | On Sliced Inverse Regression With High-Dimensional Covariates[END_REF] in Table 1 (a) is very close to one, for the first two regression models independently of the distributions of the predictors.

In the Gaussian case, Student SIR and SIR are performing equally well showing that our approach has no undesirable effects when dealing with simple cases. For configuration III -(iii), a slightly different value has been found for SIR compared to [START_REF] Dong | Robust inverse regression for dimension reduction[END_REF]. In this configuration however the trend is clear: standard SIR, Student SIR, WIRE and SIME show similar performance. In contrast, configurations I -(ii), II -(ii), III -(ii) illustrate that Student SIR can significantly outperform SIR. This is not surprising since the standard multivariate Cauchy has heavy tails and SIR is sensitive to outliers [START_REF] Gather | A robustified version of sliced inverse regression[END_REF]. Table 1 (b) illustrates on model I the effect of the sample size n: Student SIR exhibits the best performance among all methods. It is interesting to observe that, in case (ii), the smaller value of r for standard SIR does not depend on the sample size n. In contrast, adding observations results in a better estimation for Student SIR.

It is then known that SIR is not very sensitive to the number of slices h [START_REF] Liquet | A graphical tool for selecting the number of slices and the dimension of the model in SIR and SAVE approaches[END_REF]. In Table 2 (a), an analysis is performed with varying h. Student SIR appears to be as well not very sensitive to the number of slices.

Extra inliers as well as outliers can affect the estimation. In case (iii), parameter ν is controlling the extra observations magnitude. Under different values of ν = 0.5, 0.2, 0.1, 0.05 ,Table 2 (b) shows that both SIR and Student SIR are robust to inliers while CP-SIR and WCAN fail when ν is small and extra observations behave as inliers concentrated around the average.

In addition, a study on the behavior of SIR and Student SIR when X follows a standard multivariate Student distribution, with different degrees of freedom (df ), is shown in Table 3 (a). The multivariate Cauchy of model (ii) coincides with the multivariate Student with one degree of freedom. This setting is favorable to our model which is designed to handle heavy tails. Not surprisingly, Student SIR provides better results for small degrees of freedom but the difference with SIR is reduced as the degree of freedom increases and the multivariate Student gets closer to a Gaussian. The stan-dard deviation follows the same trend. In case III -(ii) the convergence of SIR becomes extremely slow. Regarding computational time, results are reported in Table 3 (b). Student SIR has multiple iterations, which increases computational time compared to SIR. It is interesting that, in the cases in which SIR fails (I -(ii), II -(ii), III -(ii) see Table 1), the convergence of Student SIR is fast, requiring less than a second on a standard laptop (Our Matlab code is available at https://hal.inria.fr/hal-01294982). All reported results have been obtained using a threshold of 0.01 for the relative increase of the Log-likelihood.

At last, the use of BIC as a selection criterion for the central subspace dimension d is investigated. As an illustration, last column of Table 3 (b) shows the number of times the criterion succeeded in selecting the correct dimension (i.e. d = 2 in this example) over 200 repetitions. BIC performs very well provided the sample size is large enough, this phenomenon being more critical as the number of outlying data increases. This is not surprising as this limitation of BIC has often been reported in the literature.

To summarize, through these simulations Student SIR shows good performance, outperforming SIR when the distribution of X is heavy-tailed (case (ii)) and preserving good properties such as insensitivity to the number of slices or robustness to inliers that are peculiar of SIR.

Real data application: The galaxy dataset

Data

The Galaxy dataset corresponds to n = 362, 887 different galaxies. This dataset has been already used in [START_REF] Chiancone | Collaborative sliced inverse regression[END_REF] with a preprocessing based on expert supervision to remove outliers. In this study all the original observations are considered, removing only points with missing values, which requires no expertise. The response variable Y is the stellar formation rate. The predictor X is made of spectral characteristics of the galaxies and is of dimension p = 46.

Evaluation setting

The number of samples n is very large and the proportion of outliers is very small compared to the whole dataset. The following strategy is adopted: 1000 random subsets of X of size n a = 3, 000, X a i , i = 1, . . . , 1000 and size n b = 30, 000, X b i , i = 1, . . . , 1000 are considered to compare the performance of SIR and Student SIR. First a reference result BSIR , Bst-SIR is obtained using the whole dataset X, using respectively SIR and Student SIR with the dimension of the e.d.r. space set to d = 3 and the number of slices to h = 1000. The value d = 3 was selected via BIC computed for d = 1, ..., 20, which is reliable for such a large sample size. The proximity measure r (26) between the two reference spaces is r( BSIR , Bst-SIR ) = 0.95. SIR and st-SIR are identifying approximately the same same e.d.r. space.

Results

Let BSIR

i , Bst-SIR i be the estimations of the basis of the e.d.r. space for the random subsets X a i , i = 1, . . . , 1000 using respectively SIR and Student SIR. The proximity measures r SIR i = r( BSIR , BSIR i ) and r st-SIR i = r( Bst-SIR , Bst-SIR i ) are considered. All results are obtained setting the number of slices to h = 10. The means (and standard deviations) of the resulting proximity measures r are respectively 0.86(0.09) for SIR and 0.87(0.09) for Student SIR. The experiment is better visualized in Figure 1 (a) where histograms show that Student SIR performs better than SIR most of the time. As expected SIR is less robust than Student SIR, obtaining with a higher frequency low values of r. The histograms show a difference between values around r = 0.96 (23.8% of random subsets for Student SIR, 17.2% for SIR).

In the second test, the sample size of the subsets is increased to n b = 30, 000. Accordingly, the number of slices is increased to h = 100. Not surprisingly, the means (and standard deviations) of r SIR i and r st-SIR i are increasing to 0.97(0.04) and 0.99(0.00). Student SIR however still performs better than SIR (Figure 1 (b)) with some low values of the proximity measure for SIR while Student SIR has almost all the values (93.4% of random subsets) concentrated around r = 0.98. The difference between the two approaches is then further emphasized in Figure 2 where the cloud of points in the upper left corner of the plot corresponds to datasets for which SIR was not able to estimate a correct basis of the e.d.r space while Student SIR shows good performance. Even if the true e.d.r space is unknown, this analysis suggests that Student SIR is robust to outliers and can be profitably used in real applications. Table 3: (a) Effect of the degree of freedom (df ) on the average of the proximity measure r (eq.( 26)) for sample size n = 200, the number of slices is h = 5 and the predictor dimension p = 10; and (b) Effect of the sample size on the computational time in seconds (standard deviations in brackets) and ratio of correct selections (d = 2) for BIC over 200 runs.

Model

Conclusion and future work

We proposed a new approach referred to as Student SIR to robustify SIR. In contrast to most existing approaches which aim at replacing the standard SIR estimators by robust versions, we considered the intrinsic characterization of SIR as a Gaussian inverse regression model [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] and modified it into a Student model with heavier tails. While SIR is not robust to outliers, Student SIR has shown to be able to deal with different kind of situations that depart from normality. As expected, when SIR provides good results, Student SIR is performing similarly but at a higher computational cost due to the need for an EM iterative algorithm for estimation.

Limitations of the approach include the difficulty in dealing with the case p > n or when there are strong correlations between variables. Student SIR as well as SIR still suffer from the need to inverse large covariance matrices. A regularization, to overcome this problem, has been proposed in [START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF] and could be extended to our Student setting. Another practical issue is how to set the dimension d of the central subspace. We have proposed the use of BIC as a natural tool in our maximum likelihood setting. It provided good results but may be not suited when the sample size is too small. A more complete study and comparison with other solutions would be interesting.

To conclude, Student SIR shows good performance in the presence of outliers and is performing equally well in case of Gaussian errors. In our experiments, the algorithm has shown fast convergence being a promising alternative to SIR since nowadays most datasets include outliers. Future work would be to extend this setting to a multivariate response following the lead of [START_REF] Barreda | Some extensions of multivariate SIR[END_REF][START_REF] Coudret | A new sliced inverse regression method for multivariate response[END_REF].

does not depend on y is independent on the distribution of ε. Indeed the proof uses only the properties of the conditional expectation seen as a projection operator. This means that in our case also, B corresponds to the mean central subspace as defined by E[X|B T X, Y = y] = E[X|B T X].

Proof of Lemma 1

The proof is adapted from the proof of lemma 1 in [START_REF] Bernard-Michel | Gaussian regularized sliced inverse regression[END_REF] taking into account the additional quantities ūi 's. Let us remark that

R def = G(µ, V, B, C) -log det V = 1 n n i=1 ūi Z T i V -1 Z i , (27)
where we have defined for i = 1, . . . , n,

Z i = µ + VBC T s i -x i (28)
= (µ -x + VBC T s) + VBC T (s i -s) -

(x i -x) ( 29
) def = Z 1 + Z 2,i -Z 3,i . (30)
Since Z 2,. and Z 3,. are centered, replacing the previous expansion in ( 27) yields

R = ū Z T 1 V -1 Z 1 + 1 n n i=1 ūi Z T 2,i V -1 Z 2,i + 1 n n i=1 ūi Z T 3,i V -1 Z 3,i -2 n n i=1 ūi Z T 2,i V -1 Z 3,i ,
where

Z T 1 V -1 Z 1 = (µ -x + VBC T s) T V -1 (µ -x + VBC T s), 1 n n i=1 ūi Z T 2,i V -1 Z 2,i = tr(C T WCB T VB), 1 n n i=1 ūi Z T 3,i V -1 Z 3,i = 1 n n i=1 ūi tr((x i -x) T V -1 (x i -x)) = 1 n n i=1
ūi tr(V -1 (x i -x)(x i -x) T ) = tr(V -1 Σ) and

1 n n i=1 ūi Z T 2,i V -1 Z 3,i = tr(C T MB),
and the conclusion follows.

G is then minimized when the λ k are the largest. As a consequence of (42), it also comes that V = Σ -V B( ĈT W Ĉ) BT V.

Replacing V B in (46) by the expression given in (37), it comes

V = Σ -M T Ĉ( ĈT W Ĉ) -1 ĈT M. ( 47
)
Grouping the results in (47), (39), (36) and the considerations after (45) gives the Proposition.

Proof of Corollary 1.

Let us remark that, under [START_REF] Schott | Determining the Dimensionality in Sliced Inverse Regression[END_REF], the coefficients W ij of W have an explicit form:

W = 1 n n i=1
ūi (s i -s)(s i -s) T

2f f t ū + f f t ū = 1 n n i=1 ūi s i s T i - f f t ū ,
where f = (f 1 , ..., f h ). Using [START_REF] Schott | Determining the Dimensionality in Sliced Inverse Regression[END_REF] the first sum corresponds to diag(f 1 , ..., f h ) leading to W = diag(f 1 , ..., f h ) -f f t ū . The inverse matrix of W can be calculated using Sherman-Morrison formula:

W -1 = diag 1 f 1 , . . . , 1 f h + 1 f h+1 O,
where O is the h × h matrix defined by O ij = 1 for all (i, j) ∈ {1, . . . , h} × {1, . . . , h}. Using (23) the jth row of M is given by: for all j = 1, . . . , h. Now taking into account that O 2 = hO, we have

M T W -1 M = h j=1 f j (x j -x)(x j -x) T + 1 f h+1 M T OM = h j=1 f j (x j -x)(x j -x) T + 1 hf h+1 (M T O)(M T O) T . (48)
Now, remarking that all the columns of M T O are equal to and thus replacing in (48) yields

M T W -1 M = h+1 j=1 f j (x j -x)(x j -x) T = Γ.
The result is then a consequence of Proposition 2.

Figure 1 :

 1 Figure 1: Histograms of the proximity measure (26) r SIR i

  i ∈ S j } -sj )(x i -x) xj T -f j xT -f j xT + f j xT = f j (x j -x) T ,

f

  j (x j -x) = h+1 j=1 f j (x j -x) -f h+1 (x h+1 -x) = -f h+1 (x h+1 -x),wheref h+1 = 1 n n i=1 ūi 1I{y i ∈ S h+1 } = ū -h j=1 f j it follows that (M T O)(M T O) T = hf 2 h+1 (x h+1 -x)(x h+1 -x) T

Table 1 :

 1 (a) Average of the proximity measure r (eq. (26)) for sample size n = 200; and (b) effect of sample size n on the average proximity measure r, both over 200 repetitions with standard deviation in brackets. Six methods are compared. SIR: sliced inverse regression; CP-SIR: contour projection for SIR; WCAN: weighted canonical correlation; WIRE: weighted sliced inverse regression estimation; SIME: sliced inverse multivariate median estimation and st-SIR: Student SIR. In all cases, the number of slices is h = 5 and the predictor dimension p = 10. Best r values are in bold.

		X	Method
				SIR	CP-SIR	WCAN	WIRE	SIME	st-SIR
		(i) .99(.01) .99(.01) .98(.01)	.98(.01) .99(.01) .99(.01)
	I	(ii)	.63(.18) .92(.04) .88(.06)	.87(.07)	.91(.04) .98(.01)
		(iii) .99(.01) .86(.12) .72(.27)	.98(.01)	.97(.01) .99(.01)
		(i) .99(.01) .98(.01) .98(.01)	.98(.01)	98(.01)	.99(.01)
	II	(ii)	.61(.18) .92(.04) .89(.06)	.87(.08)	.91(.05) .98(.01)
		(iii) .99(.01) .67(.25) .69(.28)	.98(.01)	.97(.02) .99(.01)
		(i)	.88(.06) .87(.06) .89(.05) .86(.06)	.87(.06)	.87(.06)
	III	(ii)	.40(.13) .78(.10) .78(.11)	.76(.11)	.78(.10) .85(.06)
		(iii) .84(.07) .63(.12) .67(.13) .85(.07) .85(.07) .84(.07)
							(a)
	Model	X	n	Method
				SIR		CP-SIR WCAN	WIRE	SIME	st-SIR
			50 .95(.03) .91(.09) .86(.11) .88(.11)	.90(.08) .95(.03)
		(i)	100 .98(.01) .96(.03) .96(.03) .95(.03) 200 .99(.01) .99(.01) .98(.01) .98 (.01) .99(.01) .99(.01) .96(.02) .98(.01)
			400	1(.00)		.99(.00) .99(.00) .99 (.01) .99(.00)	1(.00)
			50	.60(.22)	.66(.18) .57(.23) .49(.24)	.59(.21) .90(.07)
	I	(ii)	100 .62(.21) 200 .62(.20)	.85 (.08) .78(.11) .73(.15) .92(.04) .88(.06) .87(.07)	.81(.10) .96(.02) .91(.04) .98(.01)
			400 .62(.18)	.96(.02) .94(.03) .93(.03)	.96(.02) .99(.00)
			50 .95(.02) .45(.29) .18(.19) .73(.25)	.86(.09) .95(.02)
		(iii)	100 .98(.01) .66(.25) .35(.29) .94(.04) 200 .99(.01) .86(.12) .72(.27) .98(.01)	.94(.04) .98(.01) .97(.01) .99(.00)
			400 .99(.00) .96(.04) .96(.04) .93(.03) .99(.01) .99(.00)
							(b)

Table 2 :

 2 Effect of the number of slices (a) and of inlier magnitude ν (b) on the average proximity measure r (eq. (26)), over 200 repetitions with related standard deviation in brackets. Six methods are compared. SIR: sliced inverse regression; CP-SIR: contour projection for SIR; WCAN: weighted canonical correlation; WIRE: weighted sliced inverse regression estimation; SIME: sliced inverse multivariate median estimation and st-SIR: Student SIR. In all cases, the sample size is n = 200 and the predictor dimension p = 10. Best r values are in bold.

					Model -X	n	Method
	Model -X	df	Method				SIR	st-SIR	BIC
			SIR	st-SIR		200	.00(.00) .13(.05)	25/200
		3	.94(.05) .99(.00)		300	.01(.00) .09(.03) 109/200
	I -(ii)	5	.98(.02) .99(.00)	III-(i)	400	.04(.01) .33(.16) 156/200
		7	.98(.01) .99(.00)		500	.05(.01) .43(.17) 189/200
		10 .99(.01) .99(.00)		1000 .10(.02) .51(.17) 200/200
		3	.94(.05) .99(.00)		200	.00(.00) .13(.05)	21/200
	II -(ii)	5	.97(.02)	.99(.00)		300	.01(.00) .09(.05)	19/200
		7	.98(.01)	.99(.00)	III-(ii)	400	.04(.01) .33(.20)	39/200
		10 .99(.01) .99(.00)		500	.05(.01) .43(.18)	90/200
		3	.82(.08) .90(.04)		1000 .10(.02) .51(.20) 200/200
	III -(ii)	5	.88(.05) .92(.03)		200	.00(.00) .13(.05)	0/200
		7	.90(.04) .92(.03)		300	.01(.00) .13(.04)	12/200
		10 .90(.04) .92(.03)	III-(iii)	400	.04(.01) .38(.16)	22/200
		30 .91(.03) .92(.03)		500	.05(.01) .34(.13)	16/200
						1000 .10(.02) .51(.17

Acknowledgments

The authors would like to thank Didier Fraix-Burnet for his contribution to the data. This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).

, vertical axis r st-SIR i , i = 1, . . . , 1000, proximity measures computed using subsets of X of size n b = 30, 000. Almost all points are lying above the line y = x indicating that Student SIR improves SIR results and significantly so for the subsets in the upper left corner.

Appendix: Proofs

Proof of Proposition 1

The proof generalizes the proof of Proposition 6 in [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] to the generalized Student case. It comes from [START_REF] Fukumizu | Kernel dimension reduction in regression[END_REF] that X y follows a generalized Student distribution S p (µ y , V, α) where µ y = µ + VBc(y). Generalized Student distributions have similar properties to Gaussian distributions (see for instance section 5.5 in [START_REF] Kotz | Multivariate t Distributions and their Applications[END_REF]). In particular any affine transformation of a generalized Student distribution remains in this family. It follows that B T X|Y = y is distributed as S d (B T µ y , B T VB, α). Similarly, marginals and conditional distributions are retained in the family. It follows that X|B T X = B T x, Y = y is also a generalized Student distribution S p ( μ, Ṽ, α, γ) with

from which it is clear that Ṽ, α, γ and μ do not depend on y. It follows that

Note that for the proof of the proposition, it was necessary to show that the independence on y holds for each parameter of the distribution and not only for the mean. The independence on y of the mean is actually straightforward using [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] where it appears that the proof that E[X|B T X, Y = y]

Proof of Proposition 2

Cancelling the gradients of G(µ, V, B, C) yields the system of equations

From ( 31), we have μ

Replacing in ( 32) and ( 33) yields the simplified system of equations

It follows from the last equality that

Multiplying (37) by B T VB on the left, we get

and assuming W is regular, (38

Now, multiplying [START_REF] Barreda | Some extensions of multivariate SIR[END_REF] on the left and on the right by V and taking account of (36) entails

As a consequence of (42), it comes

and

Using this expression of V B above in (41), it comes

which means that the columns of B are stable by Σ -1 M T W -1 M and thus are eigenvectors of


Djeylan Aktas

Tao Wang

Olivier Alibart

Gian Luca Lippi
email: gian-luca.lippi@inln.cnrs.fr

Éric Picholle
email: eric.picholle@unice.fr

Sébastien Tanzilli

Quantum Coherence in Semiconductor Microlasers

Keywords: lasers, 270.5290 Photon statistics, 140.3430 Laser theory

The buildup of the coherence of a Class-B semiconductor microlaser in the subthreshold region is characterized by photon coincidence. Correlation values g (2) (0) > 2 confirm the key role of spatio-temporal dynamics in the photon statistics.

Introduction

Lasers are often perceived as threshold devices, with a rapid all-or-nothing transition around threshold. This approach usually holds for high quality factor ("Class A" [START_REF] Arecchi | Deterministic chaos in laser with injected signal[END_REF]) devices or macroscopic ones, but fails to correctly describe the transition in very small samples. The long carrier lifetime -compared to the cavity time -typical of semiconductor-based micro-and nanolasers ("Class B" devices) introduces an additional element of complexity to the problem due to the nontrivial dynamical evolution [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF][START_REF] Puccioni | Stochastic simulator for modeling the transition to lasing[END_REF] of the coherent fraction of electromagnetic field, and opens questions on coherence buildup.

Below threshold, the characterization of the latter is rendered difficult due to the small number of photons contained in the laser cavity at any given time, thus yielding very low optical intensity levels, close to or below the sensitivity limit of standard photodetectors. Fortunately, recent telecom-based quantum measurement techniques [START_REF] Ngah | Ultra-fast heralded single photon source based on telecom technology[END_REF] now allow the assessment of the laser's degree of coherence through its second order correlation function, g (2) (τ).

Experimental setup

We elected to characterize a commercial Class B laser, Thorlabs VCSEL-980 (λ = 980 nm). The attenuated beam reaches, through a standard 50/50 fiber coupler, two independent IdQuantique ID100 photon detectors. Detection events are correlated through a Ortec 567 time-to-amplitude converter (TAC). Coincidence events are accounted for and mapped. (Fig. 1) The jitter of the detectors is estimated at 40 ps, and the resolution of the TAC at 15 ps.

Results

This approach allows us to study the transition region at low pump values and to characterize the coherence buildup from a mostly incoherent regime well below threshold (Fig. 2a) towards the usual multimode regime above threshold, passing through different statistical mixtures of coherent and incoherent components. The growth of a signal at the center of the histograms, on top of the detector background, indicates the growth of short-term coherence. Its progressive enlargement with increasing pump currents concurs both with classical-based measurements above threshold [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF] and with numerical simulations of the transition region performed by taking into account intrinsic fluctuations at threshold [START_REF] Puccioni | Stochastic simulator for modeling the transition to lasing[END_REF].

The associated second order correlation function, g (2) (τ), also allows a characterization of the degree of coherence of the laser. Further increase in pump produces a flattening of the correlation (not shown), corresponding to the establishment of long-range temporal coherence shown in Fig. 3 (g (2) (0) = 1, right portion of the curve), in agreement with [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF][START_REF] Puccioni | Stochastic simulator for modeling the transition to lasing[END_REF]. Fig. 3. Experimental correlation function g (2) (0)

We also observe, for the first time in this context, correlation values g (2) (0) up to 2.14 ± 0.06 at I p = 1.30 mA (Fig. 3). Such values, significantly higher than the theoretical values for chaotic light (g (2) (0) = 2), can easily be interpreted by taking into account the temporal dynamics of the source in the extended threshold region, where bursts of shorts pulses can be classically observed [START_REF] Wang | Dynamical Buildup of Lasing in Mesoscale Devices[END_REF], in good agreement with our stochastic model [START_REF] Puccioni | Stochastic simulator for modeling the transition to lasing[END_REF], instead of the traditionally expected cw behaviour above threshold: the concentration of detection events during a pulse increases the probability of coincidence events, thus g (2) (0). The apparent decrease in correlation for very low values of the pump current (left side of Fig. 3) is an artefact stemming from the limited temporal resolution of the TAC.

Conclusion

Quantum measurement techniques based on telecom technology appear as a promising tool for the characterization of the coherence of micro-and nano-lasers at low pump levels, where classical measurements become difficult. Further developments, and the implementation of a better temporal resolution both for the photon detectors and the TDC, should also allow an indirect characterization of the temporal dynamics of these lasers. In addition, these experimental measurements validate the theoretical approach of the stochastic simulator [START_REF] Puccioni | Stochastic simulator for modeling the transition to lasing[END_REF].

Fig. 1 .

 1 Fig. 1. Experimental setup

Fig. 2 .

 2 Fig. 2. Progressive buildup of the laser coherence: Histograms of coincidence events (Δt = 0 : time-bin 78). From left to right, pumping intensities are (a) Ip = 1.15 mA; (b) Ip = 1. 20 mA & (c) Ip = 1.30 mA.