Coordinate-independent sparse sufficient dimension reduction and variable selection

Introduction

where \mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}} stands for independence and $P_{(\cdot)}$ represents the projection matrix with respect to the standard inner product, then $\mathcal{S}$ is called a dimension reduction space. The central subspace $\mathcal{S}_{y|{\mathbf{x}}}$, which is the intersection of all dimension reduction spaces, is an essential concept of SDR. Under mild conditions, it can be shown that $\mathcal{S}_{y|{\mathbf{x}}}$ is itself a dimension reduction subspace [Cook (1994, 1998a)], which we assume throughout this article, and then it is taken as the parameter of interest. The dimension $d$ of $\mathcal{S}_{y|{\mathbf{x}}}$, usually far less than $p$, is assumed to be known in this article. We also assume throughout that $n>p$.

There has been considerable interest in dimension reduction methods since the introduction of sliced inverse regression [SIR; Li (1991)] and sliced average variance estimation [SAVE; Cook and Weisberg (1991)]. Li (1992) and Cook (1998b) proposed and studied the method of principal Hessian directions (PHD), and the related method of iterative Hessian transformations was proposed by Cook and Li (2002). Chiaromonte, Cook and Li (2002) proposed partial sliced inverse regression for estimating a partial central subspace. Yin and Cook (2002) introduced a covariance method for estimating the central $k$th moment subspace. Most of these and many other dimension reduction methods are based on the first two conditional moments and as a class are called F2M methods [Cook and Forzani (2009)]. They provide exhaustive estimation of $\mathcal{S}_{y|{\mathbf{x}}}$ under mild conditions. Recently, Li and Wang (2007) proposed another F2M method called directional regression (DR). They argued that DR is more accurate than or competitive with all of the previous F2M dimension reduction proposals. In contrast to these and other moment-based SDR approaches, Cook (2007) introduced a likelihood-based paradigm for SDR that requires a model for the inverse regression of ${\mathbf{x}}$ on $y$. This paradigm, which is broadly referred to as principal fitted components (PFC), was developed further by Cook and Forzani (2009). Likelihood-based SDR inherits properties and methods from general likelihood theory and can be very efficient in estimating the central subspace.

All of the aforementioned dimension reduction methods suffer because the estimated linear reductions usually involve all of the original predictors ${\mathbf{x}}$. As a consequence, the results can be hard to interpret, the important variables may be difficult to identify and the efficiency gain may be less than that possible with variable selection. These limitations can be overcome by screening irrelevant and redundant predictors while still estimating a few linear combinations of the active predictors. Some attempts have been made to address this problem in dimension reduction generally and SDR in particular. For example, Li, Cook and Nachtsheim (2005) proposed a model-free variable selection method based on SDR. Zou, Hastie and Tibshirani (2006) proposed a sparse principal component analysis. Ni, Cook and Tsai (2005) introduced a shrinkage version of SIR, while Li and Nachtsheim (2006) suggested a sparse version of SIR. Li (2007) studied sparse SDR by adapting the approach of Zou, Hastie and Tibshirani (2006). Zhou and He (2008) proposed a constrained canonical correlation procedure ($C^{3}$) based on imposing the $L_{1}$-norm constraint on the effective dimension reduction estimates in CANCOR [Fung et al. (2002)], followed by a simple variable filtering method. Their procedure is attractive because they showed that it has the oracle property [Donoho and Johnstone (1994), Fan and Li (2001)]. More recently, Leng and Wang (2009) proposed a general adaptive sparse principal component analysis and Johnstone and Lu (2009) studied the large $p$ theory in sparse principal components analysis.

However, most existing sparse dimension reduction methods are conducted stepwise, estimating a sparse solution for a basis matrix of the central subspace column by column. Instead, in this article, we propose a unified one-step approach to reduce the number of variables appearing in the estimate of $\mathcal{S}_{y|{\mathbf{x}}}$. Our approach, which hinges operationally on Grassmann manifold optimization, is able to achieve dimension reduction and variable selection simultaneously. Additionally, our proposed method has the oracle property: under mild conditions the proposed estimator would perform asymptotically as well as if the true irrelevant predictors were known.

We start in Section 2.1 by reviewing the link between many SDR methods and a generalized eigenvalue problem disclosed by Li (2007). In Section 2.2, we describe a new SDR penalty function that is invariant under orthogonal transformations and targets the removal of row vectors from the basis matrix. Based on this penalty function, in Section 2.3, a coordinate-independent penalized procedure is proposed which enables us to incorporate many model-free and model-based SDR approaches into a simple and unified framework to implement variable selection within SDR. A fast algorithm, which combines a local quadratic approximation [Fan and Li (2001)] and an eigensystem analysis in each iteration step, is suggested in Section 2.4 to handle our Grassmann manifold optimization problem with its nondifferentiable penalty function. In Section 2.5, we describe the oracle property of our estimator. Its proof differs significantly from those in the context of variable selection in single-index models [e.g., Fan and Li (2001), Zou (2006)] because the focus here is on subspaces rather than on coordinates. Results of simulation studies are reported in Section 3, and the Boston housing data, is analyzed in Section 4. Concluding remarks about the proposed method can be found in Section 5. Technical details are given in the Appendix.

Theory and methodology

Li (2007)showed that many moment based sufficient dimension reduction methods can be formulated as a generalized eigenvalue problem in the following form:

Following Cook (2004), Li (2007) showed that the eigenvectors $\{\bolds\delta_{n1},\ldots,\break\bolds\delta_{nd}\}$ from (1) can be obtained by minimizing a least square objective function. Let

subject to $\bolds\alpha^{T}\mathbf{N}_{n}\bolds\alpha=\mathbf{I}_{d}$, where $\operatorname{tr}(\cdot)$ stands for the trace operator, $\|\cdot\|_{r}$ denotes the $L_{r}$ norm, $\tau_{2}$ is some positive constant and $\tau_{1,j}\geq 0$ for $j=1,\ldots,d$ are the lasso shrinkage parameters that need to be determined by some method like cross validation (CV). The solution $\widehat{\mathbf{V}}_{s}$ is called the sparse sufficient dimension reduction estimator. As a result of the lasso constraint, $\widehat{\mathbf{V}}_{s}$ is expected to have some elements shrunk to zero.

We can see that Li’s sparsity method is coordinate dependent because the $L_{1}$ penalty term is not invariant under the orthogonal transformation of the basis and it forces individual elements of the basis matrix $\widehat{\mathbf{V}}_{s}$ to zero. However, variable screening requires that entire rows of $\widehat{\mathbf{V}}_{s}$ be zero, which is not the explicit goal of Li’s method. To see this more clearly, partition ${\mathbf{x}}$ as $({\mathbf{x}}_{1}^{T},{\mathbf{x}}_{2}^{T})^{T}$, where ${\mathbf{x}}_{1}$ corresponds to $q$ elements of ${\mathbf{x}}$ and ${\mathbf{x}}_{2}$ to the remaining elements. If

then ${\mathbf{x}}_{2}$ can be removed, as given ${\mathbf{x}}_{1}$, ${\mathbf{x}}_{2}$ contains no further information about $y$. Let the $p\times d$ matrix $\bolds\eta$ be a basis for $\mathcal{S}_{y|{\mathbf{x}}}$ and partition $\bolds\eta=(\bolds\eta_{1}^{T},\bolds\eta_{2}^{T})^{T}$ in accordance with the partition of ${\mathbf{x}}$. Then the condition (3) is equivalent to $\bolds\eta_{2}=0$ [Cook (2004)], so the corresponding rows of the basis are zero vector.

In effect, Li’s method is designed for element screening, not variable screening. Our experience reflects this limitation and reinforces the notion that $\widehat{\mathbf{V}}_{s}$ may not be sufficiently effective at variable screening. Inspired by Li’s method, we propose a new variable screening method—called coordinate-independent sparse estimation (CISE)—in the next subsection. We will show that CISE is simpler and more effective than Li’s method at variable screening.

CISE can be applied not only to moment based SDR approaches but also model based approaches. Cook (2007) and Cook and Forzani (2008) developed several powerful model-based dimension reduction approaches, collectively referred to as principal fitted components (PFC). PFC-based SDR methods can also be formulated in the same way as (1), as summarized in the next proposition. In preparation, consider the following model for the conditional distribution of ${\mathbf{x}}$ given $y$,

2 A coordinate-independent penalty function

We define a general coordinate-independent penalty function as

Given $h_{1}=\cdots=h_{p}=\sqrt{(\cdot)}$, we have a special coordinate-independent penalty function:

A method for selecting the tuning parameters will be discussed in Section 2.6. We can see that the penalty function $\rho$ has the same form as the group lasso proposed by Yuan and Lin (2006) but their concepts and usages are essentially different. Through this article, we shall use only $\rho$ in application and theory to demonstrate our ideas.

Penalty (5) is appealing for variable selection because it is independent of the basis used to represent the span of $\mathbf{V}$, $\rho(\mathbf{V})=\rho(\mathbf{VO})$ for any orthogonal matrix $\mathbf{O}$, and because it groups the row vector coefficients of $\mathbf{V}$. This motivated us to consider the regularized function (5) that can shrink the corresponding row vectors of irrelevant variables to zero. Another appealing feature of using this penalty is its oracle property, which is discussed in Section 2.5.

3 Coordinate-independent sparse estimation

Recall the generalized eigenvalue problem (1) and the associated notation. Formally,

where $\mathbf{G}_{n}=\mathbf{N}_{n}^{-1/2}\mathbf{M}_{n}\mathbf{N}_{n}^{-1/2}$ and we use $\mathbf{G}$ to denote its population analog in what follows. Hence, the ordinary sufficient dimension reduction estimation (OSDRE) given in (2) is

By using the coordinate independent penalty function given in last subsection, we propose the following coordinate-independent sparse sufficient dimension reduction estimator (CISE):

where $\rho(\mathbf{V})$ is defined in (5).

The minimizer (6) is equivalent to $\widehat{\mathbf{V}}=\mathbf{N}_{n}^{-1/2}\widehat{\bolds\Gamma}$ where

4 Algorithm

To overcome the nondifferentiability of $\rho(\cdot)$, we adopt the local quadratic approximation of Fan and Li (2001); that is, we approximate the penalty function locally with a quadratic function at every step of the iteration as follows.

By using the second-order Taylor expansion and some algebraic manipulation, we have

where $C_{0}$ stands for a constant with respect to $\mathbf{V}$.

5 Oracle property

In what follows, without loss of generality, we assume that only the first $q$ predictors are relevant to the regression, where $d\leq q<p$. Given a $p\times d$ matrix $\mathbf{K}$, $\mathbf{K}_{(q)}$ and $\mathbf{K}_{(p-q)}$ indicate the sub-matrices consisting of its first $q$ and remaining $p-q$ rows. If $\mathbf{K}$ is $p\times p$, then the notation indicates its first $q$ and the last $p-q$ block sub-matrices. In the context of the single-index model, Fan and Li (2001) and Zou (2006) have shown that, with the proper choice of the penalty functions and regularization parameters, the penalized likelihood estimators have the oracle property. With continuous penalty functions, the coefficient estimates that correspond to insignificant predictors must shrink toward 0 as the penalty parameter increases, and these estimates will be exactly 0 if that parameter is sufficiently large. In this section, we present theorems which establish the oracle property of CISE.

The distance between the subspaces spanned by the columns of $\mathbf{V}_{n}$ and $\mathbf{V}$, denoted as $D(\mathbf{V}_{n},\mathbf{V})$, is defined as the square root of the largest eigenvalue of

This distance criterion was first used by Li, Zha and Chiaromonte (2005) in the sufficient dimension reduction setting. See Gohberg, Lancaster and Rodman (2006) for more details. We use the following assumptions to establish the oracle property.

Let $\mathbf{V}_{0}$ denote the minimizer of (6) when the population matrices $\mathbf{M}$ and $\mathbf{N}$ are used in place of $\mathbf{M}_{n}$ and $\mathbf{N}_{n}$. Then ${\mathbf{V}_{0}}_{(p-q)}=0$.

$\mathbf{M}_{n}=\mathbf{M}+O_{p}(n^{-1/2})$ and $\mathbf{N}_{n}=\mathbf{N}+O_{p}(n^{-1/2})$.

Given some mild method-specified conditions, the minimizer of (6) $\widehat{\mathbf{V}}$ is a consistent estimator of a basis of the central subspace. For example, SIR provides the consistent estimate of the central subspace given that the linearity and coverage conditions hold [Cook (1998a), Chiaromonte, Cook and Li (2002)]. Consequently, the population version $\mathbf{V}_{0}$ will be a basis of the central subspace. Therefore, Assumption 1 is a reasonable one which facilitates our following presentations. Assumption 2 is mild and typically holds. These two assumptions suffice for our main results.

We state our theorems here, but their proofs are relegated to the Appendix. The constrained objective function in the minimization problem (7) is denoted as $Q(\mathbf{V};\mathbf{M}_{n}):=f(\mathbf{V};\mathbf{M}_{n})+\rho(\mathbf{V})$ where $f(\mathbf{V};\mathbf{M}_{n})=-\operatorname{tr}(\mathbf{V}^{T}\mathbf{M}_{n}\mathbf{V})$. The first theorem establishes existence of CISE.

It is clear from Theorem 1 that by choosing the $\theta_{i}$’s properly, there exists a root-$n$ consistent CISE. The next transition theorem states an oracle-like property of CISE.

The second part of Theorem 2 is actually valid in a generalized sense. The OSDRE in the exact oracle property, denoted as $\dot{\mathbf{V}}_{n(O)}$, is obtained by using the $q\times q$ $\mathbf{M}_{n}$ and $\mathbf{N}_{n}$ formed with the first $q$ variables (denoted as $\mathbf{M}_{n(O)}$ and $\mathbf{N}_{n(O)}$). Usually, $\mathbf{N}_{n(O)}=\mathbf{N}_{n(q)}$. From the definition, it is straightforward to see that $\mathbf{M}_{n(O)}=\mathbf{M}_{n(q)}$ for the PCA, SIR and PFC methods. Thus, in these cases, Theorem 2 establishes the exact oracle property. We conjecture that $\mathbf{M}_{n(O)}$ should be very close to $\mathbf{M}_{n(q)}$ for any SDR method that satisfies Assumptions 1 and 2. From the proof of Theorem 2(ii), we can conclude that if

the exact oracle property still holds. The next result establishes that the condition above holds for DR and SAVE under certain conditions.

Suppose the linearity and constant variance conditions [Li and Wang (2007)] hold and $(na_{n})^{-1}=O_{p}(1)$. Then condition (10) is satisfied for the DR and SAVE methods.

By this proposition, Theorem 2 and the discussion above, we know that from asymptotic viewpoints CISE is effective for all of the commonly used SDR methods. We summarize this major result in the following theorem.

In this paper, we make no attempt to further analysis general conditions for the validity of (10), but we think that such studies certainly warrant future research.

6 Choice of tuning parameters

where $\widehat{\mathbf{v}}_{i}$ is the $i$th row vector of the OSDRE $\widehat{\mathbf{V}}$ defined in (6), and $r>0$ is some pre-specified parameter. Following the suggestions of Zou (2006), $r=0.5$ is used in both the simulation study and the illustration in Section 4. Such a strategy effectively transforms the original $p$-dimensional tuning parameter selection problem into a univariate one. By Lemma 2 in the Appendix, $\widehat{\mathbf{v}}_{i}$ is root-$n$ consistent. Thus, it is easily to verify that the tuning parameter defined in (11) satisfies the conditions on $a_{n}$ and $b_{n}$ needed by Theorem 2 as long as $\sqrt{n}\theta\rightarrow 0$ and ${n}^{(1+r)/2}\theta\rightarrow\infty$. Hence, it suffices to select $\theta\in[0,+\infty)$ only.

To choose the tuning parameter $\theta$, we use the following criterion which has a form similar to ones used by Li (2007) and Leng and Wang (2009):

Simulation studies

We report the results of four simulation studies in this section, three of which were conducted using forward regression models and one was conducted using an inverse regression model. We compared our method with the $C^{3}$ method [Zhou and He (2008)] and the SSIR method [Ni, Cook and Tsai (2005)]. BIC and RIC [Shi and Tsai (2002)] were used in SSIR to select the tuning parameters, and two $\alpha$ levels (0.01 and 0.005) were used in the $C^{3}$ method. We used SIR and PFC to generate $\mathbf{M}_{n}$ and $\mathbf{N}_{n}$ for CISE selection. For these methods, denoted CIS-SIR and CIS-PFC, we report only the results using the BIC criterion to select tuning parameters as we tend to believe that BIC has consistency property. Unreported simulations using the RIC criterion show slightly better performance in some cases though.

In each study, we generated 2500 datasets with the sample size $n=60$ and $n=120$. For the $C^{3}$ method, the quadratic spline with four internal knots was used, as suggested by Zhou and He (2008). Six slices were used for the SSIR method. We calculated $\mathbf{M}_{n}$ in the PFC model setting using $f(y)=(|y|,y,y^{2})^{T}$ for all simulation studies.

where $\epsilon\sim N(0,1)$, ${\mathbf{x}}=(x_{1},\ldots,x_{24})^{T}\sim N(0,\bolds\Sigma)$ with $\Sigma_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq 24$, and ${\mathbf{x}}$ and $\epsilon$ are independent. In this study, the central subspace is spanned by the direction ${\bolds\beta}_{1}=(1,1,1,0,\ldots,0)^{T}$ with twenty-one zero coefficients.

where $\epsilon\sim N(0,1)$, ${\mathbf{x}}=(x_{1},\ldots,x_{24})^{T}\sim N(0,\bolds\Sigma)$ with $\Sigma_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq 24$, and $x$ and $\epsilon$ are independent. In this study, the central subspace is spanned by the direction ${\bolds\beta}_{1}=(1,1,1,0,\ldots,0)^{T}$ with twenty-one zero coefficients. In short, this study was identical to the first, except the error was increased by a factor of $4$.

where $\epsilon\sim N(0,1)$, ${\mathbf{x}}=(x_{1},\ldots,x_{24})^{T}\sim N(0,\bolds\Sigma)$ with $\Sigma_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq 24$, and $x$ and $\epsilon$ are independent. In this study, the central subspace is spanned by the directions ${\bolds\beta}_{1}=(1,0,\ldots,0)^{T}$ and $\bolds\beta_{2}=(0,1,\ldots,0)^{T}$.

where $\bolds\epsilon\sim N(0,\mathbf{I}_{24})$, $y\sim N(0,1)$, $\Delta_{ij}=0.5^{|i-j|}$ for $1\leq i,j\leq 24$, and $y$ and $\bolds\epsilon$ are independent. The first column of $\bolds\Gamma$ is $(0.5,0.5,0.5,0.5,0,\ldots,0)^{T}$ and the second column of $\bolds\Gamma$ is $(0.5,-0.5,0.5,-0.5,0,\ldots,0)^{T}$. In this study, the central subspace is the column space of $\bolds\Delta^{-1}\bolds\Gamma$.

The simulation results from these four studies are summarized in Tables 2–5, respectively. The standard errors of the $r_{k}$’s, $\sqrt{r_{k}(1-r_{k})}/50$, are typically less than 0.01 throughout this section. In Study 1, the signal-to-noise ratio is close to 5 (the ratio of the stand deviation of $x_{1}+x_{2}+x_{3}$ to 0.5). Because of the large signal-to-noise ratio, all the considered methods show very good performance, but CIS-SIR, CIS-PFC and $C^{3}$ perform slightly better than SSIR. In Study 2, we decreased the signal-to-noise ratio to about 1.2 and now CIS-SIR and CIS-PFC perform much better than $C^{3}$ and SSIR. In both Studies 3 and 4, CISE is generally superior to the other two methods, especially for CIS-PFC and the rate $r_{3}$. It should be pointed out that the superiority of CISE becomes more significant when $n$ gets larger. When $n=120$, $C^{3}$ still cannot perform exact identifications well, while SSIR rarely identifies all relevant and irrelevant variables correctly.

While both CISE and $C^{3}$ have the oracle property, they differ in many aspects. CISE is a unified method that can be applied to many popular sufficient dimension reduction methods, including PCA, PFC, SIR, SAVE and DR. On the other hand, $C^{3}$ is based on one specified sufficient dimension reduction method, canonical correlation [Fung et al. (2002)]. We regard $r_{3}$, the estimated probability all relevant and irrelevant variables are identified correctly, as the most important aspect of a method. On that measure CISE typically dominates $C^{3}$. There was only one case (Table 1, $n=60$) in which $C^{3}$ did slightly better than CISE. Additionally, CISE seems conceptually simpler and is easily implemented.

Boston housing data

We applied our method to the Boston housing data, which has been widely studied in the literature. The Boston housing data contains 506 observations, and can be downloaded from the web site http://lib.stat.cmu.edu/datasets/boston_corrected.txt. The re-sponse variable $y$ is the median value of owner-occupied homes in each of the 506 census tracts in the Boston Standard Metropolitan Statistical Areas. The 13 predictor variables are per capita crime rate by town ($x_{1}$); proportion of residential land zoned for lots over 25,000 sq.ft ($x_{2}$); proportion of nonretail business acres per town ($x_{3}$); Charles River dummy variable ($x_{4}$); nitric oxides concentration ($x_{5}$); average number of rooms per dwelling ($x_{6}$); proportion of owner-occupied units built prior to 1940 ($x_{7}$); weighted distances to five Boston employment centers ($x_{8}$); index of accessibility to radial highways ($x_{9}$); full-value property-tax rate ($x_{10}$); pupil–teacher ratio by town ($x_{11}$); proportion of blacks by town ($x_{12}$); percentage of lower status of the population ($x_{13}$).

Previous studies suggested that we remove those observation with crime rate greater than 3.2, as a few predictors remain constant except for 3 observations in this case [Li (1991)]. So we used the 374 observations with crime rate smaller than 3.2 in this analysis. All the methods considered in Section 3 were applied to this dataset. Scatter-plotting of each predictor against $y$, we concluded that it would be sufficient to use $\mathbf{f}=(\sqrt{y},y,y^{2})^{T}$ in the PFC model. Since PFC is a scale-invariant method, we did not standardize the data as many other methods do. Similar to the previous studies in the literature, we pick up two directions to estimate the central subspace. The estimated bases of the central subspace for all the considered methods are summarized in Table 4.1.

The coefficients in Table 4.1 from CIS-SIR, CIS-PFC and SSIR are based on the original dataset, while the coefficients of $C^{3}$ is based on a data-specific weighted version [Zhou and He (2008)]. As suggested by CIS-PFC, explanatory variables $x_{6}$, $x_{7}$, $x_{10}$, $x_{11}$, $x_{12}$ and $x_{13}$ would be important in explaining $y$.

In Table 7, we used the bootstrap to assess the accuracy of variable selection for all methods except $C^{3}$, as it is not clear how the weighting procedure used by Zhou and He should be automated. Without weighting we encountered serious convergence problems in the $C^{3}$ algorithm. This bootstrap study can be considered as another simulation study.

The bootstrap procedure was conducted as follows. First, we randomly chose with replacement 374 observations for $y$ jointly with $x_{6}$, $x_{7}$, $x_{10}$, $x_{11}$, $x_{12}$ and $x_{13}$. Secondly, we separately randomly selected 374 observations for $x_{1}$, $x_{2}$, $x_{3}$, $x_{4}$, $x_{5}$, $x_{8}$ and $x_{9}$. Then we combine them to make one complete bootstrap dataset. In this way, we mimic the results of the analysis of original data, forcing $x_{1}$, $x_{2}$, $x_{3}$, $x_{4}$, $x_{5}$, $x_{8}$ and $x_{9}$ to be irrelevant. This procedure was repeated 2500 times. The resulting rates $r_{1}$, $r_{2}$ and $r_{3}$ are shown in Table 7. The results show a pattern similar to those in simulation studies and again CISE performed quite well.

The establishment of the oracle property in this paper takes advantage of the simple trace form of the objective function: $-\operatorname{tr}(\mathbf{V}^{T}\mathbf{M}_{n}\mathbf{V})$. However we believe that the proof in the Appendix can be extended to more general objective functions. Moreover, it is also of great interests to see whether CISE and its oracle property are still valid in high-dimensional settings in which $p>n$.

We have seen that $\mathbf{N}_{n}$ usually takes the form of the marginal sample covariance matrix of ${\mathbf{x}}$, while $\mathbf{M}_{n}$ depends on the specific method. In practice, how to choose $\mathbf{M}_{n}$ for variable selection is an important issue and merits thorough investigation. In addition, it is well demonstrated that for the multiple regression model, the BIC criterion tends to identify the true sparse model well if the true model is included in the candidate set [Wang, Li and Tsai (2007)]. The consistency of the BIC criterion proposed in Section 2.6 deserves further study as well.

In order to prove the theorems, we first state a few necessary lemmas. For notation convenience, we need the following additional definitions. Define the Stiefel manifold $\operatorname{St}(p,d)$ as

The tangent space $T_{\bolds\Gamma}(p,d)$ of $\bolds\Gamma\in\operatorname{St}(p,d)$ is defined by

If $\mathbf{Z}\in T_{\bolds\Gamma}(p,d),{\bolds\Gamma}\in\operatorname{St}(p,d)$, we have: {longlist}

$R(\bolds\Gamma+t\mathbf{Z})=\bolds\Gamma+t\mathbf{Z}-(1/2)t^{2}\bolds\Gamma\mathbf{Z}^{T}\mathbf{Z}+O(t^{3})$.

This lemma comes from Lemma 10 and Proposition 12 of Manton (2002).

where ${\bolds\Gamma_{0}}$ denotes any minimizer of (8) when $\mathbf{G}_{n}$ is taken as the population matrix $\mathbf{G}$.

Note that $\bolds\Gamma_{0}=\mathbf{N}^{1/2}\mathbf{V}_{0}$, and thus it is equivalent to show that

since $D(\bolds{\Gamma}_{*},{\bolds\Gamma}_{0})=O_{p}(n^{-1/2})$ and $D(\cdot,\cdot)$ satisfies the triangle inequality.

From Proposition 20 of Manton (2002), it is straightforward to see

For any small $\epsilon$, if we can show that there exits a sufficiently large constant $C$, such that

By using Lemma 1, for $\mathbf{Z}\in T_{\bolds\Gamma_{*}}(p,d)$ we have

where the second inequality holds because ${\mathbf{1}}_{j}\mathbf{N}_{n}^{-1/2}{\bolds\Gamma}_{*}=0$ for any $j>q$ by Assumption 1, and the last inequality comes from first-order Taylor expansion and the definition of $a_{n}$. In addition, according to the theorem’s condition $\sqrt{n}a_{n}\stackrel{{\scriptstyle p}}{{\rightarrow}}0$, we known that $\Delta_{2}$ is $o_{p}(1)$. Furthermore, based on Lemma 1 and Assumption 2, we have

where $\bolds\Lambda=\operatorname{diag}\{\bolds\Lambda_{1},\bolds\Lambda_{2}\}$ is the diagonal eigenvalue matrix of $\mathbf{G}$ with the first $d\times d$ sub-matrix $\bolds\Lambda_{1}$. By using the definition of $\mathbf{Z}$ in (5), we get

where we use the fact $\operatorname{tr}(\mathbf{A}^{T}\mathbf{A}\bolds\Lambda_{1})-\operatorname{tr}(\mathbf{A}\mathbf{A}^{T}\bolds\Lambda_{1})=0$ because $\mathbf{A}$ is skew-symmetric. Here the last inequality follows from basic properties of trace operator for semi-positive definite matrix. As a consequence, by the Cauchy–Schwarz inequality for trace operator, the third term in $\Delta_{1}$ is uniformly bounded by $\|\mathbf{B}\|_{s}\times\|\sqrt{n}({\mathbf{G}}_{n}-\mathbf{G}){\bolds{\Gamma}}_{0}\|_{s}$. Therefore, as long as the constant $C$ is sufficiently large, the first two terms in $\Delta_{1}$ will always dominate the third term and $\Delta_{2}$ with arbitrarily large probabilities. This implies inequality (13), and the proof is complete. {pf*}Proof of Theorem 2 (i) To prove this part, we need represent (7) as vector forms. Define

where $\mathbf{t}_{i}$ denotes the $i$th column vector of $\mathbf{V}$, $\mathbf{C}_{l}$’s are $pd\times pd$ block-diagonal matrices, $\mathbf{C}_{kl}$’s $pd\times pd$ block matrices, $\mathbf{C}_{l}$ and $\mathbf{C}_{kl}$ contain $\mathbf{N}_{n}$ in the $l$th diagonal block and in the $(k,l)$ as well as $(l,k)$ blocks, respectively. The $pd\times pd$ symmetric matrices $\mathbf{C}_{kl}$ are defined for all the pairs of different indices belonging to $\mathcal{J}$, given by the $d(d-1)/2$ combinations of the indices $1,\ldots,d$.

where $\mathbf{A}$ is a $pd\times pd$ block-diagonal matrix with all diagonal blocks $\mathbf{M}_{n}$. Of course, in the above equation each $\mathbf{v}_{i}$ is regarded as a function of $\mathbf{t}$.

By using the equality representation of the compact Stiefel manifolds $\operatorname{St}(p,d)$, (7) is equivalent to

As a consequence, this enables us to apply an improved global lagrange multiplier rule proposed by Rapcsák (1997).

is a $(pd\times[d(d+1)/2])$-dimensional matrix, and $\partial^{g}f(\mathbf{V}_{n})/\partial\mathbf{t}$ and $\partial^{g}\rho(\mathbf{V}_{n})/\partial\mathbf{t}$ are defined in a similar form of $\partial^{g}Q^{*}(\mathbf{t})/\partial\mathbf{t}$ by replacing $Q^{*}$ with $f$ and $\rho$, respectively. By Theorem 1 and noting that $\partial f(\mathbf{V}_{n})/\partial\mathbf{t}$ is linear in $\mathbf{t}$,

where $\kappa_{1},\ldots,\kappa_{d-1d}$ are a sequence of constants satisfy they are not all the zeros. Define a sequence of $pd$-dimensional vectors $\mathbf{z}_{ij}$’s,

Similarly, $\kappa_{i}=o_{p}(1)$. Consequently, we can conclude all the $\kappa_{i}$ and $\kappa_{ij}$ equal to zero in probability which yields contradiction. As a result, with probability tending to 1 (w.p.1), (5) cannot hold, which implies there exists $j>q$ so that

(ii) For convenience purposes, first decompose the matrix $\mathbf{M}_{n}$ and $\mathbf{N}_{n}$ into the following block form:

where $\mathbf{M}_{n(q)}$ and $\mathbf{N}_{n(q)}$ are the first $q\times q$ sub-matrices. It then follows that

and $\mathbf{G}_{n(q)}=\mathbf{N}_{n(q)}^{-1/2}\mathbf{M}_{n(q)}\mathbf{N}_{n(q)}^{-1/2}$. Note that

where $2a_{n}^{-1}\operatorname{tr}(\mathbf{Z}^{T}\mathbf{G}_{n(q)}{\widehat{\bolds\Gamma}}_{n(O)})=0$ by using Lemma 2, and

Now we first deal with $\mathbf{M}_{n1}$. Rewrite it as

Let ${\mathbf{P}}_{\bolds\Omega}(\bolds\Sigma_{n})=\bolds\Omega(\bolds\Omega^{T}\bolds\Sigma_{n}\bolds\Omega)^{-1}\bolds\Omega^{T}\bolds\Sigma_{n}$ and let ${\mathbf{Q}}_{\bolds\Omega}(\bolds\Sigma_{n})=\mathbf{I}_{p}-{\mathbf{P}}_{\bolds\Omega}(\bolds\Sigma_{n})$. Then

By construction, $\mathcal{S}_{y|{\mathbf{x}}}\subseteq\operatorname{span}(\bolds\Omega)$. Under certain conditions [Cook (1998a)], we know $\operatorname{span}\{\bolds\Sigma^{-1}[\bolds\Sigma-\operatorname{Var}({\mathbf{x}}|y)]\}\subseteq\mathcal{S}_{y|{\mathbf{x}}}$. Hence,

By (21) and (22), $\mathbf{M}_{n(O)i}=\mathbf{M}_{n(q)i}+O_{p}(n^{-1})$ for $i=3,\ldots,6$, can be proved in a similar fashion to the foregoing proofs. We omit the details here for saving some space. It follows that for the DR method,

Thus, condition (10) is satisfied as long as $(na_{n})^{-1}=O_{p}(1)$.

Note that for SAVE, $\mathbf{M}_{n}$ takes the form of $\mathbf{M}_{n1}$ for DR. Thus, condition (10) is also satisfied for SAVE.

Acknowledgments

The authors thank the Associate Editor and two anonymous referees for their many helpful comments that have resulted in significant improvements in the article. Specially, we are grateful to the Associate Editor for helping us complete the proof of Proposition 3. The authors would also like to thank Dr. Jianhui Zhou and Dr. Liqiang Ni for providing us the codes for computing the $C^{3}$ and SSIR estimators.

References