Variable screening in multivariate linear regression with high-dimensional covariates

Abstract

We propose two variable selection methods in multivariate linear regression with high-dimensional covariates. The first method uses a multiple correlation coefficient to fast reduce the dimension of the relevant predictors to a moderate or low level. The second method extends the univariate forward regression of Wang [(2009). Forward regression for ultra-high dimensional variable screening. Journal of the American Statistical Association, 104(488), 1512–1524. https://doi.org/10.1198/jasa.2008.tm08516] in a unified way such that the variable selection and model estimation can be obtained simultaneously. We establish the sure screening property for both methods. Simulation and real data applications are presented to show the finite sample performance of the proposed methods in comparison with some naive method.

Keywords:

• Dimension reduction
• forward regression
• multiple correlation coefficient
• multivariate regression
• variable selection

1. Introduction

High-dimensional multivariate regression has been widely applied in bioinformatics, chemometrics, and medical image analysis where many of the response variables are highly correlated (Cai et al., 2013; Ferte et al., 2013; Jia et al., 2017; Peng et al., 2010; Smith & Fahrmeir, 2007). For instance, in genetics study, we are interested in the association between correlated phenotypes (involved in biological pathways) and genotypes, as genetic effects and their possible interaction have been recognized as an important component for the genetic architecture of each complex phenotype (Yi, 2010). For this kind of problem, the number of covariates or explanatory variables is much larger than the number of observations or samples. Traditional methods of subset selection and stepwise procedure become infeasible when confronted with high dimensionality (Breiman, 1995).

Statistical methods and theories have been developed to solve this problem through various approaches such as network-based regularization method (C. Li & Li, 2008; Ren et al., 2019, 2017), graphical model (B. Li, Chuns et al., 2012; Yin & Li, 2011), correlation-based screening (B. Li, Chuns et al., 2012; Song et al., 2016) and group lasso (Y. Li et al., 2015; J. Wang et al., 2019; Yang & Zou, 2015).

Variable selection methods for regression models with a univariate response have been proposed in the past. Some popular methods include the bridge regression (Frank & Friedman, 1993; Fu, 1998), LASSO (Tibshirani, 1996), SCAD (Fan & Li, 2001), LARS (Efron et al., 2004), elastic net (Zou & Hastie, 2005), adaptive LASSO (H. H. Zhang & Lu, 2007; Zou, 2006), and Dantzig selector (Candes & Tao, 2007; Y. Kong et al., 2016), among others. On the other hand, variable screening procedures have been developed to reduce the dimensionality from an ultrahigh dimension to a lower dimension which is smaller than the sample size (Fan & Lv, 2008; X. Kong et al., 2017; G. Li, Peng et al., 2012; H. Wang, 2009; Zhu et al., 2011).

For variable selection under multivariate regression models, one simple approach is to apply some variable selection method to univariate regression of each response separately. Such an approach may produce sub-optimal results since it does not utilize the joint information among the responses (Breiman & Friedman, 1997; Kim et al., 2009). To improve the estimation, various attempts have been made. One approach is to use dimension reduction techniques such as the reduced rank regression (Chen & Huang, 2012; He et al., 2018; Zhao et al., 2017) and the sliced inverse regression (Setdji & Cook, 2004; N. Zhang et al., 2019). Another approach is to use a block-structured regularization method to select a subset which can be used as predictors for all outcome variables (Obozinski et al., 2011; Peng et al., 2010; Turlach et al., 2005). The latter approach assumes that a covariate affects either all or none of the responses. However, this assumption may be too strong when each response variable is affected by different sets of predictors. Rothman et al. (2010) proposed a penalized framework to estimate multivariate regression coefficient and covariance matrix simultaneously under ${\ell }_1$ penalty. Lee and Liu (2012) further improved Rothman et al.'s (2010) work by using a weighted ${\ell }_1$ regularization. Cai et al. (2013) proposed a method to first estimate the regression coefficients in a column-wise fashion with Dantzig selector and then to estimate the precision matrix by solving a constrained ${\ell }_1$ minimization problem.

In high-dimensional setting, most of the aforementioned multivariate regression methods use the technique of regularization to estimate the regression coefficient matrix (Obozinski et al., 2011; Peng et al., 2010; Turlach et al., 2005). However, a well-chosen penalty requires an efficient exploration of the correlation structure of the responses. It is reported that simultaneously estimating covariance and selecting variables via joint optimization can be numerically unstable in high-dimensional cases (Deshpande et al., 2019; Pecanka et al., 2019; Ren et al., 2019).

In this study, we propose two methods in parallel for variable screening and variable selection, namely the multiple correlation coefficient (MCC) screening (Section 3) and the unified forward regression (UFR) (Section 4). The first method is for dimension reduction which filters out covariates that have weak correlation with the response variables. It significantly reduces the feature space to a moderate or low dimension that covers the set of relevant predictors almost certainly. The second method is for variable selection which uses an extended forward regression (FR) (H. Wang, 2009) to identify all relevant predictors consistently under mild conditions. By MCC all relevant predictors are identified or screened, whereas by UFR both variable selection and model estimation are obtained. We illustrate the finite sample performance of the proposed methods in comparison with a naive method by simulation (Section 5) and a real data application (Section 6). We conclude the paper in Section 7 and defer the technical proofs in Appendix.

2. Notation and assumptions

Let $\mathbf{y}=(y_1,y_2,\dots ,y_q{)}^{\mathrm{\top }}$ denote the q-dimensional response vector of interest. Let $\mathbf{x}=(x_1,x_2,\dots ,x_p{)}^{\mathrm{\top }}$ denote the p-dimensional covariates or predictors. Denote the covariance matrices of $\mathbf{y}$ and $\mathbf{x}$ by ${\mathrm{\Sigma }}_y$ and ${\mathrm{\Sigma }}_x=\left({\sigma }_{ij}\right)$, respectively. Without loss of generality, assume that $\mathrm{E}\left(x_k\right)=0$ and $\text{var}\left(x_k\right)=1$ for $k=1,\dots ,p$ and that $\mathrm{E}\left(y_j\right)=0$ for $j=1,\dots ,q$. In practice, these can be achieved by standardization and centralization.

Consider the multivariate linear regression model (1) $\mathbf{y}=B^{\mathrm{\top }}\mathbf{x}+\overrightarrow{ϵ},$(1) where B is a $p\times q$ matrix of coefficients and $\overrightarrow{ϵ}$ is the random error vector which is independent with $\mathbf{x}$. For $j=1,\dots ,q$ and $k=1,\dots ,p$, denote ${\overrightarrow{\text{β}}}_j$ as the jth column vector of B and ${\overrightarrow{\text{β}}}_{\left(k\right)}$ as the kth row vector of B. If ${\overrightarrow{\text{β}}}_{\left(k\right)}\ne \mathbf{0}$, $x_k$ is referred to as a relevant predictor.

Let $F=\{1,\dots ,p\}$ denote the full model of predictors. Let $S=\{k:{\overrightarrow{\text{β}}}_{\left(k\right)}\ne \mathbf{0}\}$ denote the true model. Denote the compliment of S by $S^c$. Denote the cardinalities of F and S as $\left|F\right|=p$ and $\left|S\right|=p_0$, respectively. Throughout, let $\Vert \cdot \Vert$ denote the Euclidean norm of a vector.

Let $\left\{\right({\mathbf{y}}_i,{\mathbf{x}}_i):i=1,\dots ,n\}$ denote independent and identically distributed samples of $(\mathbf{y},\mathbf{x})$. Denote $X_{n\times p}=({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n{)}^{\mathrm{\top }}$ and $Y_{n\times q}=({\mathbf{y}}_1,\dots ,{\mathbf{y}}_n{)}^{\mathrm{\top }}$. For $j=1,\dots ,q$, let ${\mathbf{y}}_{\left(j\right)}$ denote the jth column of Y.

Assume that $\mathbf{x}$ is high dimensional with p being much larger than the sample size n (in the sense of Cai & Lv, 2007). Assume that the response vector is associated with only a small portion of predictors, i.e., $p_0/p$ is small and $p_0$ is $\mathrm{O}\left(n\right)$ (Fan & Lv, 2008). This sparsity principle is frequently adopted and deemed useful in analysis.

3. Multiple correlation coefficient

We first propose to use a multiple correlation coefficient (MCC) to identify S. It is known that the multiple correlation coefficient between $\mathbf{y}$ and $x_k$ is defined as ${\rho }_k={max}_{\overrightarrow{\alpha }\in {\mathbf{R}}^q}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}({\overrightarrow{\alpha }}^{\mathrm{\top }}\mathbf{y},x_k)$ and its square can be further expressed as (2) ${\rho }_k^2=\mathrm{E}\left({\overrightarrow{\gamma }}_k^{\mathrm{\top }}\mathbf{y}{x}_k\right),$(2) where ${\overrightarrow{\gamma }}_k={\mathrm{\Sigma }}_y^{-1}\mathrm{E}\left(\mathbf{y}{x}_k\right)$ (Anderson, 2003, Section 12.2).

Given the standardized samples, we estimate ${\rho }_k^2$ by (3) ${\hat{\rho }}_k^2=\frac{1}{n}\sum _{i=1}^n{{\hat{\overrightarrow{\gamma }}}_k}^{\mathrm{\top }}{\mathbf{y}}_i{x}_{ik},$(3) where ${\hat{\overrightarrow{\gamma }}}_k=(\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{\mathbf{y}}_i^{\mathrm{\top }}{)}^{-1}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}$. Note that the computation of ${\hat{\rho }}_k^2$ is simple and fast through matrix algebra and does not involve any iteration. Then, we estimate S by ${\hat{S}}_{\mathrm{M}\mathrm{C}\mathrm{C}}=\{k:{\hat{\rho }}_k^2\ge \tau \}$, where τ is the threshold which determines the size of the estimated predictors. Here we adopt the threshold of Fan and Lv (2008) by choosing $\tau ={\hat{\rho }}_{(p-d_n+1)}^2$, where ${\hat{\rho }}_{\left(1\right)}^2\le \cdots \le {\hat{\rho }}_{\left(p\right)}^2$ are the order statistics and $d_n=\lceil n/\log \left(n\right)\rceil$ ($\lceil \cdot \rceil$ is the ceiling function), so that $d_n$ predictors with the largest values of ${\hat{\rho }}_k^2$ are retained. The naive correlation coefficient (NCC) method of Fan and Lv (2008) estimates S by ${\hat{S}}_{\mathrm{N}\mathrm{C}\mathrm{C}}={\cup }_{j=1}^p\{k:{\hat{\rho }}_{k,j}^2\ge {\tau }_j\},$ where ${\hat{\rho }}_{k,j}$ is the sample correlation coefficient between $y_j$ and $x_k$ and ${\tau }_j$ is determined in the same way as in MCC with respect to the jth response.

We now show that the MCC-based screening procedure has the sure screening property (i.e., the probability of selecting all true relevant predictors tends to one) and reduces the dimensionality of predictors below the sample size.

We state some assumptions first.

Assumption 3.1

Let ${\lambda }_{min}\left(A\right)$ and ${\lambda }_{max}\left(A\right)$ denote the smallest and largest eigenvalue of a positive definite matrix A, respectively. Assume that there exist two positive constants ${\tau }_{min}<{\tau }_{max}$ such that $2{\tau }_{min}<{\lambda }_{min}\left({\mathrm{\Sigma }}_y^{-1}\right)\le {\lambda }_{max}\left({\mathrm{\Sigma }}_y^{-1}\right)<2^{-1}{\tau }_{max}$and $2{\tau }_{min}<{\lambda }_{min}\left({\mathrm{\Sigma }}_x\right)\le {\lambda }_{max}\left({\mathrm{\Sigma }}_x\right)<2^{-1}{\tau }_{max}.$

Assumption 3.2

Assume that (i) for $j=1,\dots ,q$, $\Vert {\overrightarrow{\text{β}}}_j\Vert \le C_B$ for some positive constant $C_B$ and that (ii) for $k=1,\dots ,p$, ${\text{β}}_{min}=\underset{k\in S}{min}\underset{j}{min}\left|{\text{β}}_{kj}\right|\ge {\nu }_B{n}^{-{\xi }_{min}}$ for some positive constants ${\xi }_{min}$ and ${\nu }_B$.

Assumption 3.3

Assume that there exist positive constants $0<\eta <4^{-1}$, K such that (i) $n^{-1}\log \left(pq\right)\le \eta$, and (ii) $\mathrm{E}\left(e^{t{x}_k^2}\right)\le K$ for $\left|t\right|\le \eta$ and all $k=1,\dots ,p$.

Assumption 3.1 requires the matrix X to be well behaved. Assumption 3.2 requires the smallest nonzero regression coefficient does not converge too fast. Otherwise, it cannot be consistently identified. (See Fan & Peng, 2004 for more discussions.) Assumption 3.3 ensures the exponential convergence rate of arbitrary order moments of $\mathbf{x}$ and $\overrightarrow{ϵ}$ (Cai et al., 2011) which is superior to the polynomial type counterpart (Ravikumar et al., 2010).

Theorem 3.1

Under Assumptions 3.1–3.3, if ${\rho }_k^2\ge \tau$ for all $k\in S$, then $\mathrm{P}(S\subset {\hat{S}}_{\mathrm{M}\mathrm{C}\mathrm{C}})\to 1$ as $n\to \mathrm{\infty }$.

Theorem 3.1 reveals that for a properly chosen threshold τ, the probability that MCC detects all relevant predictors tends to one.

4. Unified forward regression

In this section, we propose a unified forward regression (UFR) for variable selection. It extends Wang's (2009) forward regression method for the multivariate response case.

Let $M=\{k_1,\dots ,k_t\}$ denote a generic subset of F with $\left|M\right|=t$. Denote ${\mathbf{x}}_{\left(M\right)}=(x_{k_1},\dots ,x_{k_t}{)}^{\mathrm{\top }}$ and denote $X_{\left(M\right)}=({\mathbf{x}}_{1\left(M\right)},\dots ,{\mathbf{x}}_{n\left(M\right)}{)}^{\mathrm{\top }}$ as the subset of X corresponding to M. We first describe a naive forward regression (NFR) method that combines the selected variables obtained by repeatedly applying Wang's (2009) forward regression method to univariate regressions with respect to every response. The procedure is summarized as follows. Initially, set $S_{\left(j\right)}^{\left(0\right)}=\mathrm{\varnothing }$ for $j=1,\dots ,q$. Perform forward regression with respect to the jth response by iterating the following two steps for $\ell =1,\dots ,n$.

1. For every $k\in F\mathrm{\setminus }{S}_{\left(j\right)}^{(\ell -1)}$, let $M_{kj}^{(\ell -1)}=S_{\left(j\right)}^{(\ell -1)}\bigcup \left\{k\right\}$. Compute the sum square of residuals ${\mathrm{R}\mathrm{S}\mathrm{S}}_{kj}^{(\ell -1)}={\mathbf{y}}_{\left(j\right)}^{\mathrm{\top }}(I_n-{\tilde{H}}_k^{(\ell -1)}){\mathbf{y}}_{\left(j\right)}$, where ${\tilde{H}}_{kj}^{(\ell -1)}=X_{\left(M_{kj}^{(\ell -1)}\right)}(X_{\left(M_{kj}^{(\ell -1)}\right)}^{\mathrm{\top }}{X}_{\left(M_{kj}^{(\ell -1)}\right)}{)}^{-1}{X}_{\left(M_{kj}^{(\ell -1)}\right)}^{\mathrm{\top }}$. Let $a_{\ell j}={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{k\in F\mathrm{\setminus }{S}_{\left(j\right)}^{(\ell -1)}}{\mathrm{R}\mathrm{S}\mathrm{S}}_{kj}^{(\ell -1)}$.

2. Update $S_{\left(j\right)}^{\left(\ell \right)}=S_{\left(j\right)}^{(\ell -1)}\bigcup \left\{a_{\ell j}\right\}$.

The solution path of NFR is obtained by $\{S_{\mathrm{N}\mathrm{F}\mathrm{R}}^{\left(\ell \right)}={\cup }_{j=1}^q{S}_{\left(j\right)}^{\left(\ell \right)}:\ell =1,\dots ,n\}$.

Next, we propose the unified forward regression to select predictors by applying a modified forward regression algorithm that makes use of all response variables simultaneously. The procedure is modified from the previous one as follows. Initially, set $S^{\left(0\right)}=\mathrm{\varnothing }$. Perform a modified forward regression by iterating the following two steps for $\ell =1,\dots ,n$.

1. For every $k\in F\mathrm{\setminus }{S}^{(\ell -1)}$, let $M_k^{(\ell -1)}=S^{(\ell -1)}\bigcup \left\{k\right\}$. Compute the sum square of residuals ${\mathrm{R}\mathrm{S}\mathrm{S}}_k^{(\ell -1)}=\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-{\tilde{H}}_k^{(\ell -1)}\left)Y\right\}$, where ${\tilde{H}}_k^{(\ell -1)}=X_{M_k^{(\ell -1)}}(X_{M_k^{(\ell -1)}}^{\mathrm{\top }}{X}_{M_k^{(\ell -1)}}{)}^{-1}{X}_{M_k^{(\ell -1)}}^{\mathrm{\top }}$. Let $a_{\ell }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{k\in F\mathrm{\setminus }{S}^{(\ell -1)}}{\mathrm{R}\mathrm{S}\mathrm{S}}_k^{(\ell -1)}$.

2. Update $S^{\left(\ell \right)}=S^{(\ell -1)}\bigcup \left\{a_{\ell }\right\}$.

The solution path of UFR is obtained by $\{S_{\mathrm{U}\mathrm{F}\mathrm{R}}^{\left(\ell \right)}=S^{\left(\ell \right)}:\ell =1,\dots ,n\}$. Notice that both NFR and UFR terminate automatically after n iterations. It is seen that the UFR algorithm makes use of all response variables simultaneously by the trace operator. It has nearly one qth computation cost of NFR.

We show that the proposed UFR method also possesses the sure screening property. Also, we add a few more assumptions to facilitate the development of the theory.

Assumption 4.1

Assume that (i) $\mathbf{x}$ follows elliptically contoured distribution, whose density admits the form $\left|{\mathrm{\Sigma }}_x{|}^{-1/2}g\right\{(\mathbf{x}-\overrightarrow{\mu }{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_x^{-1}(\mathbf{x}-\overrightarrow{\mu }\left)\right\}$ with $\overrightarrow{\mu }=\mathrm{E}\mathbf{x}$ and $g(\cdot )>0$, denoted by $\mathrm{E}\mathrm{C}(\overrightarrow{\mu },{\mathrm{\Sigma }}_x,g)$, and that (ii) the distribution of $\overrightarrow{ϵ}$ is normal.

Assumption 4.2

There exist positive constants ξ, ${\xi }_0$ and ν such that (i) $\log \left(p\right)\le \nu n^{\xi }$, (ii) $p_0\le \nu n^{{\xi }_0}$, and (iii) $\xi +6{\xi }_0+12{\xi }_{min}<1$.

Assumption 4.3

The row vectors of B, i.e, ${\overrightarrow{\text{β}}}_{\left(k\right)}$, $k=1,\dots ,p$, have the same ‘all-or-nothing’ structure, i.e., the entries of ${\overrightarrow{\text{β}}}_{\left(k\right)}$ are either all zero or none-zero.

Usually, the normality assumption of $\mathbf{x}$ is imposed to facilitate theory development (Fan & Lv, 2008; H. Wang, 2009). Here in Assumption 4.1, we relax it to elliptically contoured distribution and show its sufficiency to obtain Lemma 1 of H. Wang (2009) in Appendix. Assumption 4.1, together with Assumption 3.1, ensures the sparse Riesz assumption (C. Zhang & Huang, 2008) to derive some key inequalities in proving Theorem 4.1. Assumption 4.2 has been popularly assumed in the literature of ultra-high dimensional inference (Fan & Lv, 2008; H. Wang, 2009). It implies that the dimension of the covariates diverges to infinity at an exponential rate (Fan & Lv, 2008). Assumption 4.3 implies that all responses are associated with the same covariates (Turlach et al., 2005). It warrants the row-wise selection of B by UFR in contrast to the element-wise selection by NFR, which enables UFR to reach the sure screening property in fewer steps than NFR.

Define $\mathcal{K}=2{\tau }_{max}\nu C_B^2{\tau }_{min}^{-2}{\nu }_B^{-4}$, where the factors are defined in Assumptions 3.1, 3.2 and 4.2. Applying Theorem 1 of H. Wang (2009), we can readily get $\mathrm{P}(S\subset S_{\mathrm{N}\mathrm{F}\mathrm{R}}^{\left(q\mathcal{K}\nu n^{2{\xi }_0+4{\xi }_{min}}\right)})\to 1$, i.e., the NFR selects all relevant predictors with high probability after $q\mathcal{K}\nu n^{2{\xi }_0+4{\xi }_{min}}$ steps for the multivariate regression setting. While the following theorem shows that the UFR can do the job in much fewer (one qth) steps.

Theorem 4.1

Under Assumptions 3.1–4.3, $\mathrm{P}(S\subset S_{\mathrm{U}\mathrm{F}\mathrm{R}}^{\left(K\nu n^{2{\xi }_0+4{\xi }_{min}}\right)})\to 1$ as $n\to \mathrm{\infty }$.

We adopt the BIC criteria to select the best subset of variables from a solution path (Liang et al., 2012; H. Wang, 2009). Let (4) $\begin{array}{rl}\mathrm{B}\mathrm{I}\mathrm{C}\left(M\right)& =\log \left[n^{-1}\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(M\right)}\left)Y\right\}\right]\\ & +n^{-1}\left|M\right|(\log n+2\log p),\end{array}$(4) where $H_{\left(M\right)}=X_{\left(M\right)}\{X_{\left(M\right)}^{\mathrm{\top }}{X}_{\left(M\right)}{\}}^{-1}{X}_{\left(M\right)}^{\mathrm{\top }}$. We then choose the subset of the variables from the solution path which minimizes $\mathrm{B}\mathrm{I}\mathrm{C}\left(M\right)$. The selection consistency was showed by H. Wang (2009) and Sofer et al.(2014).

Note that the UFR method is consistent if $p=O\left(n^{\alpha }\right)$ for some $\alpha >0$, while the MCC method works with $\log \left(p\right)=O\left(n^{\alpha }\right)$. In this sense, the MCC method can handle higher dimensional variable screening than the UFR method. Secondly, the UFR method is computationally more expensive than the MCC method as the former involves n−1 more times of matrix inversion operation than the latter. On the other hand, when p and n are of the same order, both MCC and UFR perform well in terms of coverage probability and UFR performs better in yielding a parsimonious model with high specificity in terms of model size and correct fit (defined later) as seen in simulation.

5. Simulation

We conduct numerical studies to investigate the finite sample performance of the proposed methods, i.e., MCC and UFR, in comparison with the naive correlation coefficient (NCC) method and the naive forward regression (NFR).

5.1. Models

Consider five models for generating the p-dimensional covariates $\mathbf{x}$ in Table 1, which are adopted from Examples 1 and 2 of Fan and Lv (2008), Example 1 of Tibshirani (1996), and Examples 4 and 5 of H. Wang (2009), respectively.

Table 1. Five models.

Model	n	q	$p_0$	${\mathrm{\Sigma }}_x$	Entries of $B_0$
1	200	4	8	$I_p$	$N(4\log (n)/\sqrt{n},1)$
2	75	5	3	$0.5I_p+0.5{\mathbf{1}}_p{\mathbf{1}}_p^{\mathrm{\top }}$	$N(4\log (n)/\sqrt{n},1)$
3	200	3	3	${\sigma }_{kr}={0.5}^{|k-r|},1\le k,r\le p$	$\mathrm{\Gamma }(2,1)$
4	300	2	5	$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{I_{p_0},4^{-1}(I_{p-p_0}+p_0{\mathbf{1}}_{(p-p_0)}{\mathbf{1}}_{(p-p_0)}^{\mathrm{\top }}\left)\right\}$	$\exp \left(9\right)$
5	200	6	10	$I_p$	$\exp \left(9\right)$

For models 1 to 3, $\mathbf{x}$ follows a multivariate normal distribution with zero mean vector and covariance matrix ${\mathrm{\Sigma }}_x$ of the structure of identity, autoregressive and compound symmetry, respectively. In model 4, $\mathbf{x}$ is generated by $x_r=(z_r+w_r)/\sqrt{2}$ for $r=1,\dots ,p_0$ and $x_r=(z_r+\sum _{r^{\prime }=1}^{p_0}{z}_{r^{\prime }})/2$ for $r=p_0+1,\dots ,p$, where $z_r$ and $w_r$ are independent standard normal variables (H. Wang, 2009). Note that model 4 is a challenging case as the correlation coefficient of the relevant predictors and the response variables are much smaller than the correlation coefficient of irrelevant predictors and the response variables. The details of ${\mathrm{\Sigma }}_x$ are provided in Table 1. In model 5, we generate independent components of both $\mathbf{x}$ and $\overrightarrow{ϵ}$ to be e−1 where e is an exponential random variable with parameter one. This model is used to examine the robustness of the proposed methods against the departure from the normality assumption. Consider the number of predictors p to be 1000, 5000, and 10,000, respectively, which are all much larger than the sample sizes considered in the five models. Recall $p_0$ is the number of relevant predictors. Denote the first $p_0$ rows of B by $B_0$. We generate independent entries of $B_0$ from distributions given in the last column of Table 1, where $N(4\log (n)/\sqrt{n},1)$ is a normal random variable with mean $4\log \left(n\right)/\sqrt{n}$ and variance 1, $\mathrm{\Gamma }(2,1)$ denotes a random variable of gamma distribution with shape parameter 2 and scale parameter 1, and $\exp \left(9\right)$ is an exponential random variable with parameter 9. They are all independent with $\mathbf{x}$. Set the remaining entries (the last $p-p_0$ rows) of B to be zero.

For the multivariate response case, the signal-to-noise ratio is given by $R^2=\frac{\mathrm{t}\mathrm{r}\left\{\mathrm{v}\mathrm{a}\mathrm{r}\right(B^{\mathrm{\top }}\mathbf{x}\left)\right\}}{\mathrm{t}\mathrm{r}\left\{\mathrm{v}\mathrm{a}\mathrm{r}\right(\mathbf{y}\left)\right\}}=\frac{\mathrm{t}\mathrm{r}\left(B^{\mathrm{\top }}{\mathrm{\Sigma }}_{x}B\right)}{p{\sigma }^2+\mathrm{t}\mathrm{r}\left(B^{\mathrm{\top }}{\mathrm{\Sigma }}_{x}B\right)}.$We chose the values of ${\sigma }^2$ such that the signal-to-noise ratios are $30\mathrm{\%}$, $60\mathrm{\%}$, and $90\mathrm{\%}$, respectively.

Throughout, set the number of replications N to be 1000.

5.2. Evaluation criteria

For MCC screening, we use Fan and Lv's (2008) hard threshold method to retain the relevant predictors. For both NFR and UFR, we use the BIC criterion (4(4) $\begin{array}{rl}\mathrm{B}\mathrm{I}\mathrm{C}\left(M\right)& =\log \left[n^{-1}\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(M\right)}\left)Y\right\}\right]\\ & +n^{-1}\left|M\right|(\log n+2\log p),\end{array}$(4) ) to determine the relevant predictors.

We adopt five measures as described in Table 2 to evaluate the finite sample performance of the proposed methods, where the model size (MS) is the number of the selected relevant predictors, the coverage probability (CP) measures how likely all the relevant predictors are identified, the percentage of correctly fitted (CF) measures the capability in identifying the true model correctly, the correct zero (CZ) characterizes the capability in producing sparse solution, and the incorrect zero (IZ) characterizes the method's under-fitting effects. Ideally, we wish a method to have MS close to $p_0$, CP, CF, CZ all close to 100% and IZ close to zero.

Table 2. Measures for the finite sample performance of variable selection.

Model size	$\mathrm{M}\mathrm{S}=\left|\hat{S}\right|$
Coverage probability	$\mathrm{C}\mathrm{P}=\mathrm{P}(S\subset \hat{S})$
% of correctly fitted model	$\mathrm{C}\mathrm{F}=\mathrm{P}(S=\hat{S})$
% of correct zero	$\mathrm{C}\mathrm{Z}=(p-p_0{)}^{-1}|S^{\mathrm{c}}\cap {\hat{S}}^{\mathrm{c}}|$
% of incorrect zero	$\mathrm{I}\mathrm{Z}=p_0^{-1}|S\cap {\hat{S}}^{\mathrm{c}}|$

For $b=1,\dots ,N$, let ${\hat{B}}^{\left(b\right)}$ denote the estimate of B under the bth replication. The corresponding selected model is denoted by ${\hat{S}}^{\left(b\right)}=\{k:{\widehat{\overrightarrow{\text{β}}}}_{\left(k\right)}^{\left(b\right)}\ne \mathbf{0},k=1,\dots p\}$. The empirical MS is computed as $\mathrm{M}\mathrm{S}=N^{-1}\sum _{b=1}^N\left|{\hat{S}}^{\left(b\right)}\right|$ and the empirical values of the other measures are similarly computed.

5.3. Results

Tables 3–7 report the finite sample performance of the four competing methods in terms of the measures given in Table 2 under various numbers of covariates p and signal strength $R^2$. We summarize the findings as follows. (i) The MCC method is uniformly superior to the NCC method with larger coverage probability (CP), better estimation of sparsity (with larger CZ and smaller IZ), as expected. (ii) As we adopted the fixed threshold procedure for MCC and NCC, these two methods produce conservatively large coverage of predictors at the cost of large model size. For the same reason, the percentage of incorrect zero is larger than the other two regression-based methods (UFR and NFR). So the resulting percentages of correctly fitted models for MCC and NCC are zero. (iii) When comparing UFR with NFR, the UFR demonstrates its superiority over NFR uniformly in all five measures across all five models (including Model 5 with the non-normal distribution). This corroborates the advantage of UFR in utilizing the correlation within responses over NFR. When comparing UFR with PWL, both methods perform comparably when the signal strength is as small as $30\mathrm{\%}$. When the signal strength is as large as $60\mathrm{\%}$ or $90\mathrm{\%}$, UFR outperforms PWL in all five measures in general. (iv) The UFR method performs inferior to the MCC method in cases of ultra-high dimensional covariates especially under lower signal strength, as pointed out earlier. For instance, in Model 1 of Table 3, the coverage probability of UFR reduces by $83\mathrm{\%},$ while the counterpart of MCC reduces by $33\mathrm{\%}$ when the dimension of predictors p increases from 1000 to 10,000 at the signal strength of $30\mathrm{\%}$. (v) As for the impact of the signal strength, the percentage of incorrect zeros rises under the weak signal strength cases from those under the strong signal strength cases. It is consistent with the findings for the univariate case in H. Wang (2009) and Y. Li et al. (2017). However, as the signal strength increases (e.g., from $30\mathrm{\%}$ to $90\mathrm{\%}$), the percentages of coverage probability (CP) and probability of correct fit (CF) increase significantly (e.g., $61.9\mathrm{\%}$ to $98.3\mathrm{\%}$ and $28.8\mathrm{\%}$ to $58.8\mathrm{\%}$, respectively, with p = 5000) and the percentage of incorrect zeros (IZ) drops quickly (e.g., from $53.7\mathrm{\%}$ to $2.35\mathrm{\%}$ with p = 5000) by both NFR and UFR as seen in Table 3. (vi) To examine the impact of the sample size, Table 8 reports the performance of the proposed methods under Model 1 with a number of covariates p fixed at 5000 and varying sample size n to be 100, 200, and 400, respectively. It is seen that the measures of model size (MS), coverage probability (CP), probability of correct fit (CF) and probability of incorrect zero (IZ) are sensitive to sample size. The improvement of performance is significant. For instance, when the sample size increases from n = 100 to n = 200 with signal strength ${\mathrm{R}}^2=60\mathrm{\%}$, the CP increases from $52.2\mathrm{\%}$ to $80.4\mathrm{\%}$ on average and the percentage of incorrect zero drops from $47.8\mathrm{\%}$ to $19.7\mathrm{\%}$ on average.

Table 3. Five measures of the performance of variable selection defined in Table 2 obtained by the four competing methods under various numbers of covariates (p) and signal-to-noise ratio ($R^2$) for Model 1 in Table 1 with $(n,q,p_0)=(200,4,8)$.

Method	p	$R^2$ (%)	MS	CP (%)	CF (%)	CZ (%)	IZ (%)
MCC	1000	30	38	75.9	0	98.4	24.8
		60	38	96.7	0	98.6	3.3
		90	38	99.5	0	98.6	0.5
	5000	30	38	58.2	0	99.6	41.8
		60	38	93.2	0	99.7	6.8
		90	38	98.5	0	99.7	1.5
	10,000	30	38	50.8	0	99.8	49.2
		60	38	90.6	0	99.9	9.4
		90	38	98.1	0	99.9	1.9
NCC	1000	30	38	26.2	0	97.9	73.8
		60	38	40.6	0	98.9	59.2
		90	38	50.5	0	98.2	49.5
	5000	30	38	13.9	0	99.6	86.1
		60	38	30.7	0	99.6	69.3
		90	38	41.8	0	99.6	58.2
	10,000	30	38	10.9	0	99.8	89.1
		60	38	25.6	0	99.8	74.4
		90	38	37.9	0	99.8	62.0
NFR	1000	30	1.1	1.5	0	99.9	98.5
		60	8.1	12.9	0	99.9	87.1
		90	8.3	99.7	90.9	99.9	0.3
	5000	30	2.1	0.5	0	99.9	99.4
		60	6.5	6.2	0	99.9	93.8
		90	7.8	95.3	76.9	100	4.7
	10,000	30	3.1	0.4	0	100	99.6
		60	7.4	4.6	0	100	95.4
		90	10.3	87.9	65.8	100	12.1
UFR	1000	30	0.3	4.1	0	100	95.9
		60	7.6	94.9	71.9	100	5.1
		90	8.0	100	100	100	0
	5000	30	0.1	1.1	0	100	98.9
		60	6.9	86.4	49.9	100	13.6
		90	8.0	100	99.9	100	0
	10,000	30	0.1	0.7	0	100	99.3
		60	5.9	73.8	34.2	100	26.1
		90	8.0	100	99.9	100	0

Table 4. Five measures of the performance of variable selection defined in Table 2 obtained by the four competing methods under various numbers of covariates (p) and signal to noise ratio ($R^2$) for Model 2 in Table 1 with $(n,q,p_0)=(75,5,3)$.

Method	p	$R^2\left(\mathrm{\%}\right)$	MS	CP (%)	CF (%)	CZ (%)	IZ (%)
MCC	1000	30	18	80.7	0	98.4	19.3
		60	18	98.1	0	98.5	1.9
		90	18	99.7	0	98.5	0.3
	5000	30	18	68.1	0	99.7	31.9
		60	18	95.5	0	99.7	4.5
		90	18	99.6	0	99.7	0.4
	10,000	30	18	61.3	0	99.8	38.7
		60	18	95.1	0	99.8	4.9
		90	18	99.7	0	99.9	0.3
NCC	1000	30	18	79.7	0	96.5	23.3
		60	18	95.8	0	96.5	4.7
		90	18	98.9	0	96.5	1.4
	5000	30	18	68.9	0	99.3	39.1
		60	18	94.3	0	99.3	7.3
		90	18	98.7	0	99.3	3.4
	10,000	30	18	68.1	0	99.6	51.9
		60	18	89.9	0	99.6	10.1
		90	18	95.3	0	99.7	2.1
NFR	1000	30	1.6	7.9	0.1	99.9	92.1
		60	3.2	51.6	23.7	99.9	48.4
		90	3.4	99.2	93.2	99.9	0.8
	5000	30	1.2	3.5	0	99.9	96.5
		60	3.9	35.6	11.3	99.9	64.4
		90	4.6	97.8	88.1	99.9	2.2
	10,000	30	1.6	2.3	0	99.9	97.4
		60	3.8	28.6	6.7	99.9	71.3
		90	4.7	96.4	84.7	99.9	3.6
UFR	1000	30	1.5	45.6	11.7	99.9	54.4
		60	3.0	99.8	99.1	100	0.2
		90	3.0	100	100	100	0
	5000	30	1.0	28.8	2.4	99.9	71.2
		60	2.9	98.9	97.5	100	1.0
		90	3.0	100	100	100	0
	10,000	30	0.9	24.4	2.0	99.9	75.6
		60	3.0	98.6	95.6	100	1.4
		90	3.0	100	100	100	0

Table 5. Five measures of the performance of variable selection defined in Table 2 obtained by the four competing methods under various numbers of covariates (p) and signal-to-noise ratio ($R^2$) for Model 3 in Table 1 with $(n,q,p_0)=(200,3,3)$.

Method	p	$R^2\left(\mathrm{\%}\right)$	MS	CP (%)	CF (%)	CZ (%)	IZ (%)
MCC	1000	30	38	100	0	96.5	0
		60	38	100	0	96.5	0
		90	38	100	0	96.5	0
	5000	30	38	100	0	99.3	0
		60	38	100	0	99.3	0
		90	38	100	0	96.3	0
	10,000	30	38	99.9	0	99.6	0.1
		60	38	100	0	99.6	0
		90	38	100	0	99.6	0
NCC	1000	30	38	29.5	0	98.3	70.5
		60	38	47.8	0	98.3	52.2
		90	38	54.8	0	98.4	45.2
	5000	30	38	18.2	0	99.7	81.8
		60	38	38.1	0	99.6	61.9
		90	38	52.0	0	99.7	48.0
	10,000	30	38	13.3	0	99.8	86.7
		60	38	36.0	0	99.8	63.9
		90	38	50.6	0	99.8	49.4
NFR	1000	30	2.2	99.9	0.2	100	60.1
		60	3.0	67.8	17.8	99.9	32.2
		90	4.7	90.6	71.6	99.9	9.4
	5000	30	2.2	37.0	0	100	63.0
		60	3.0	65.0	13.7	100	35.0
		90	3.7	89.8	69.2	100	10.2
	10,000	30	2.1	36.2	0.1	100	63.8
		60	3.9	63.1	11.7	100	36.9
		90	5.7	89.2	67.4	100	10.8
UFR	1000	30	2.1	71.5	22.6	100	28.5
		60	2.9	96.4	89.1	100	3.6
		90	3.0	99.9	99.6	100	0
	5000	30	2.0	67.2	14.4	100	32.8
		60	2.8	94.1	82.3	100	5.9
		90	3.0	99.9	99.7	100	0.1
	10,000	30	1.9	63.7	10.9	100	36.3
		60	2.8	93.6	80.6	100	6.4
		90	3.0	99.9	99.5	100	0

Table 6. Five measures of the performance of variable selection defined in Table 2 obtained by the four competing methods under various numbers of covariates (p) and signal to noise ratio ($R^2$) for Model 4 in Table 1 with $(n,q,p_0)=(300,2,5)$.

Method	p	$R^2\left(\mathrm{\%}\right)$	MS	CP (%)	CF (%)	CZ (%)	IZ (%)
MCC	1000	30	53	21.7	0	97.1	78.3
		60	53	27.8	0	97.1	72.2
		90	53	37.7	0	97.2	79.7
	5000	30	53	19.8	0	99.4	80.2
		60	53	26.7	0	99.4	73.1
		90	53	36.1	0	99.4	63.9
	10,000	30	53	17.7	0	99.7	82.3
		60	53	25.2	0	99.7	74.8
		90	53	35.1	0	99.7	64.9
NCC	1000	30	53	19.8	0	97.1	80.2
		60	53	20.7	0	97.1	79.2
		90	53	21.8	0	97.1	78.1
	5000	30	53	16.8	0	99.4	83.1
		60	53	18.6	0	99.4	81.3
		90	53	19.7	0	99.4	80.3
	10,000	30	53	16.2	0	99.7	83.8
		60	53	17.3	0	99.7	82.7
		90	53	20.0	0	99.7	80.0
NFR	1000	30	2.7	15.6	0	99.9	84.4
		60	4.9	34.3	0	99.9	65.7
		90	5.0	65.1	5.0	99.9	34.9
	5000	30	2.6	13.2	0	99.9	86.8
		60	4.8	27.5	0	99.9	72.5
		90	5.5	57.8	3.1	99.9	42.3
	10,000	30	2.0	12.1	0	99.9	87.9
		60	4.8	26.6	0	99.9	73.4
		90	6.5	54.7	1.9	99.9	45.3
UFR	1000	30	2.4	35.9	0	99.9	64.1
		60	4.4	75.8	10.3	99.9	24.2
		90	5.0	96.8	70.8	99.9	3.2
	5000	30	2.1	28.9	0	99.9	71.1
		60	4.2	70.9	7.2	99.9	29.1
		90	5.0	95.5	63.2	99.9	4.5
	10,000	30	2.0	24.7	0	99.9	75.3
		60	4.1	69.3	4.9	99.9	30.7
		90	5.0	95.6	61.9	99.9	4.4

Table 7. Five measures of the performance of variable selection defined in Table 2 obtained by the four competing methods under various numbers of covariates (p), and signal-to-noise ratio ($R^2$) for Model 5 in Table 1 with $(n,q,p_0)=(200,6,10)$.

Method	p	$R^2\left(\mathrm{\%}\right)$	MS	CP (%)	CF (%)	CZ (%)	IZ (%)
MCC	1000	30	38	85.2	0	97.0	14.8
		60	38	97.2	0	97.1	2.8
		90	38	99.8	0	97.2	0.2
	5000	30	38	70.1	0	99.4	29.9
		60	38	93.3	0	99.4	6.7
		90	38	99.5	0	99.4	0.6
	10,000	30	38	63.0	0	99.7	37.0
		60	38	90.6	0	99.7	9.4
		90	38	99.3	0	99.7	0.8
NCC	1000	30	38	61.2	0	96.8	38.8
		60	38	78.8	0	96.9	21.2
		90	38	85.8	0	97	14.2
	5000	30	38	40.8	0	99.3	59.2
		60	38	64.9	0	99.4	24.3
		90	38	75.7	0	99.4	24.3
	10,000	30	38	33.0	0	99.6	66.9
		60	38	60.0	0	99.7	39.9
		90	38	72.5	0	99.7	27.5
NFR	1000	30	3.8	7.6	0	99.9	92.4
		60	10.2	50.9	0.2	99.9	49.0
		90	10.7	96.5	63.9	99.9	3.5
	5000	30	2.5	4.3	0	99.9	95.7
		60	6.4	40.2	0.1	99.9	59.8
		90	10.6	94.6	49.9	99.9	5.4
	10,000	30	2.4	3.3	0	99.9	96.7
		60	10.6	35.3	0.1	99.9	64.7
		90	11.6	93.9	44.3	99.9	6.1
UFR	1000	30	2.7	27.3	0	100	72.7
		60	9.9	99.3	92.8	100	0.7
		90	10.0	100	99.9	100	0
	5000	30	1.3	13.3	0	100	86.7
		60	4.1	98.2	84.5	100	1.8
		90	10.0	100	99.9	100	0
	10,000	30	1.0	10.1	0	100	89.9
		60	9.8	97.8	81.3	100	2.2
		90	10.0	100	99.9	100	0

Table 8. Five measures of the performance of variable selection defined in Table 2 obtained by the four competing methods under various sample sizes (n) and signal-to-noise ratio ($R^2$) for Model 1 in Table 1 with $(p,q,p_0)=(5000,4,8)$.

Method	n	$R^2\left(\mathrm{\%}\right)$	MS	CP (%)	CF (%)	CZ (%)	IZ (%)
MCC	100	30	22	39.3	0	99.6	60.3
		60	22	72.6	0	99.7	27.4
		90	22	92.3	0	99.7	7.7
	200	30	38	81.2	0	99.4	18.8
		60	38	95.5	0	99.4	4.5
		90	38	98.9	0	99.4	1.1
	400	30	67	96.7	0	98.8	3.3
		60	67	99.3	0	98.8	0.7
		90	67	99.8	0	98.8	0.1
NCC	100	30	22	24.7	0	99.6	75.3
		60	22	45.9	0	99.6	54.0
		90	22	57.4	0	99.6	42.6
	200	30	38	48.3	0	99.3	51.8
		60	38	65.2	0	99.3	34.8
		90	38	70.9	0	99.4	29.1
	400	30	67	62.1	0	98.8	37.9
		60	67	72.2	0	98.8	27.8
		90	67	76.7	0	98.8	23.3
NFR	100	30	2.1	1.4	0	99.9	98.6
		60	6.2	14.8	0	99.9	85.2
		90	8.3	78.8	14.2	99.9	21.2
	200	30	3.0	10.6	0	99.9	89.4
		60	8.0	50.0	0.2	99.9	49.9
		90	8.9	84.8	24.0	100	15.2
	400	30	7.7	33.3	30.0	100	66.6
		60	8.0	65.9	2.2	100	34.1
		90	9.0	86.8	30.7	100	13.1
UFR	100	30	1.0	3.2	0	100	96.8
		60	5.6	69.2	12.0	100	30.8
		90	8.0	99.9	99.1	100	0.1
	200	30	2.8	35.5	0.1	100	64.5
		60	7.6	95.6	69.8	100	4.4
		90	8.0	99.9	99.7	100	0
	400	30	6.4	79.5	12.3	100	20.5
		60	7.9	98.4	88.1	100	1.5
		90	8.0	100	100	100	0

In conclusion, the MCC method performs better when the dimension of covariates is ultra-high ($\log \left(p\right)=O\left(n^{\alpha }\right)$) with respect to the sample size and the UFR method outperforms the MCC method when the dimension of covariates is of polynomial order ($p=O\left(n^{\alpha }\right)$).

6. Real data application

We apply the proposed methods to a real data set regarding bone mineral density (BMD) (Reppe et al., 2010). The data were collected from 84 post-menopausal Caucasian women aged from 50 to 86. For each subject, there are two responses, namely the body mass index and total hip z-score (a measure of how strong the bone in the hip), and 8649 gene expression levels in trans-iliacal bone biopsies served as covariates. It is known that low bone mineral density is usually related to fragile bone and osteoporosis and progressive reduction of bone strength which leads to increasing susceptibility of bone fractures (Cooper, 1997; Reppe et al., 2010). The goal of the study is to identify the genes that are related to BMD.

Table 9 reports the genes identified by the five competing methods. The MCC method identified 19 genes which include all 13 genes identified by NFR except gene TNK2. The PWL method identified 12 genes which all identified by NFR except PAIP1. And the UFR found 10 significant genes which are all contained in the set identified by NFR.

Table 9. Selected genes for the BMD data.

method	genes
MCC	ACSL3, NIPSNAP3B, DLEU2, C1ORF61, DKK1, SOST, ABCA8,
	AFFX-M27830-M-at, RRNF216, PLIN5, PACS2, MUM1, CRYGS,
	RABEP2, PEX14, USP2, FKBP14, FAM55C, MAPK8
NCC	ACSL3, NIPSNAP3B, DLEU2, C1ORF61, DKK1, SOST, ABCA8,
	AFFX-M27830-M-at, RRNF216, PLIN5, MEPE, COPS4, CRYGS
	RABEP2, CTSE, PEX14, FKBP14, NARG1, PAIP1
NFR	ACSL3, NIPSNAP3B, DLEU2, C1ORF61, DKK1, SOST, ABCA8,
	AFFX-M27830-M-at, RNF216, PLIN5, FKBP14, FAM55C, TNK2
UFR	ACSL3, NIPSNAP3B, DLEU2, C1ORF61, DKK1, SOST, ABCA8,
	AFFX-M27830-M-at, RNF216, PLIN5

To examine the quality of variable selection of these methods, we compare the prediction mean square error ($\mathrm{E}\Vert \mathbf{y}-B_S^{\mathrm{\top }}{\mathbf{x}}_S{\Vert }^2$) obtained by the three methods. To this end, we randomly split the data into a training set of 60 samples and a testing set of the remaining 24 samples. The average prediction mean square errors over the 100 replications for MCC, NCC, NFR and UFR are 273.7, 293.5, 271.0 and 241.6, respectively. Clearly, the UFR method is the winner. All the eight genes (ACSL3, NIPSNAP3B, DLEU2, C1ORF61, DKK1, SOST, ABCA8, and AFFX-M27830-M-at) identified by Reppe et al. (2010) were selected by the four competing methods. The UFR method discovered two more genes, RNF216 and PLIN5, with the smallest prediction mean square errors. Similar to the findings in simulation, both MCC and NCC selected more genes than NFR and UFR with larger prediction error.

7. Conclusion

We propose two methods for variable screening in high-dimensional multivariate linear regression. The MCC method has the advantage of computational ease and can provide fast variable screening to obtain an accurate subset with a dimension below the ample size. The proposed UFR method has the feature of discovering all relevant predictors consistently at nearly the same computational cost as the univariate forward regression. The performance of UFR is sensitive to the dimensionality and signal strength. Our theory assumes Gaussian distribution for the response variables. The numerical study also shows the robustness of the proposed methods against non-normality. It is of interest to investigate the problem under more general non-homogeneously sparse assumption and nonlinear models.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendix

A.1. Proof of Theorem 3.1

Let $C_K=\tau +{\eta }^{-1}K$ and $\vartheta =3{\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)}$. When ${\rho }_k^2\ge \tau$ for all $k\in S$, it suffices to show that under Assumptions 3.1–3.3 (A1) $\begin{array}{rl}& \mathrm{P}(\underset{1\le k\le p}{max}|{\hat{\rho }}_k^2-{\rho }_k^2|\ge \vartheta )\\ & =\mathrm{P}(\underset{1\le k\le p}{max}\left|\frac{1}{n}\sum _{i=1}^n{\hat{\overrightarrow{\gamma }}}_k^{\mathrm{\top }}{\mathbf{y}}_i{x}_{ik}-\mathrm{E}{{\overrightarrow{\gamma }}_k}^{\mathrm{\top }}\mathbf{y}{x}_k\right|\ge \vartheta )\\ & <3p^{-\tau }{q}^{-\tau }.\end{array}$(A1) Observe that (A2) $\begin{array}{rl}& |{\hat{\rho }}_k^2-{\rho }_k^2|\\ & =\left|({\hat{\overrightarrow{\gamma }}}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}+{\overrightarrow{\gamma }}_k^{\mathrm{\top }}(\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}-\mathrm{E}\mathbf{y}{x}_k)\right|\\ & \le \left|({\hat{\overrightarrow{\gamma }}}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|+\left|{\overrightarrow{\gamma }}_k^{\mathrm{\top }}(\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}-\mathrm{E}\mathbf{y}{x}_k)\right|\\ \\ & \triangleq A_1+A_2.\end{array}$(A2) First, we compute the rate of convergence for $A_2$. For $k=1,\dots ,p$, let ${\omega }_k=n^{-1}\sum _{i=1}^n{\overrightarrow{\gamma }}_k^{\mathrm{\top }}({\mathbf{y}}_i{x}_{ik}-\mathrm{E}\mathbf{y}{x}_k)$. (So $A_2=\left|{\omega }_k\right|$.) Let $t_1=\eta \sqrt{n^{-1}\log \left(pq\right)}$. Applying the inequalities $\mathrm{P}\left(\right|U|\ge V)\le e^{-t_{1}V}\mathrm{E}{e}^{t_1\left|U\right|}$ for any V >0 and $|e^s-1-s|\le s^2{e}^{max(s,0)}$ for any $s\in \mathbb{R}$, we obtain (A3) $\begin{array}{rl}& \mathrm{P}(A_2\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le e^{-C_K\log \left(pq\right)}\mathrm{E}\exp \left(t_1\right|{\omega }_k\left|\right)\\ & \le \exp \{-C_K\log (pq)+t_1^2\mathrm{E}{\omega }_k^2\exp (t_1\left|{\omega }_k\right|\left)\right\}\\ & \le \exp \{-C_K\log (pq)+{\eta }^{-1}K\log (pq\left)\right\}\\ & =p^{-\tau }{q}^{-\tau }.\end{array}$(A3) Second, we compute the rate of convergence for $A_1$. Notice that (A4) $\begin{array}{rl}\left|\right({\hat{\gamma }}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}|& \le \left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\\ & +\left|{\gamma }_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|.\end{array}$(A4) Parallel to (A3(A3) $\begin{array}{rl}& \mathrm{P}(A_2\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le e^{-C_K\log \left(pq\right)}\mathrm{E}\exp \left(t_1\right|{\omega }_k\left|\right)\\ & \le \exp \{-C_K\log (pq)+t_1^2\mathrm{E}{\omega }_k^2\exp (t_1\left|{\omega }_k\right|\left)\right\}\\ & \le \exp \{-C_K\log (pq)+{\eta }^{-1}K\log (pq\left)\right\}\\ & =p^{-\tau }{q}^{-\tau }.\end{array}$(A3) ), we have (A5) $\begin{array}{rl}& \mathrm{P}(\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le e^{-C_K\log \left(pq\right)}\mathrm{E}\exp \left\{t_1\left(\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\right)\right\}\\ & \le \exp \{-C_K\log (pq)+t_1^2\mathrm{E}{\left({\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right)}^2\\ & \times \exp \left(t_1\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\right)\}\\ & \le \exp \{-C_K\log (pq)+{\eta }^{-1}K\log (pq\left)\right\}\\ & =p^{-\tau }{q}^{-\tau }.\end{array}$(A5) Similar to (A5(A5) $\begin{array}{rl}& \mathrm{P}(\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le e^{-C_K\log \left(pq\right)}\mathrm{E}\exp \left\{t_1\left(\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\right)\right\}\\ & \le \exp \{-C_K\log (pq)+t_1^2\mathrm{E}{\left({\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right)}^2\\ & \times \exp \left(t_1\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\right)\}\\ & \le \exp \{-C_K\log (pq)+{\eta }^{-1}K\log (pq\left)\right\}\\ & =p^{-\tau }{q}^{-\tau }.\end{array}$(A5) ), (A6) $\mathrm{P}(\left|{\gamma }_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\le p^{-\tau }{q}^{-\tau }.$(A6) Combining (A4(A4) $\begin{array}{rl}\left|\right({\hat{\gamma }}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}|& \le \left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\\ & +\left|{\gamma }_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|.\end{array}$(A4) ), (A5(A5) $\begin{array}{rl}& \mathrm{P}(\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le e^{-C_K\log \left(pq\right)}\mathrm{E}\exp \left\{t_1\left(\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\right)\right\}\\ & \le \exp \{-C_K\log (pq)+t_1^2\mathrm{E}{\left({\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right)}^2\\ & \times \exp \left(t_1\left|{\hat{\gamma }}_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\right)\}\\ & \le \exp \{-C_K\log (pq)+{\eta }^{-1}K\log (pq\left)\right\}\\ & =p^{-\tau }{q}^{-\tau }.\end{array}$(A5) ) and (A6(A6) $\mathrm{P}(\left|{\gamma }_k^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\le p^{-\tau }{q}^{-\tau }.$(A6) ), (A7) $\begin{array}{rl}& \mathrm{P}(A_1=\left|\right({\hat{\gamma }}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}|\ge 2{\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le 2p^{-\tau }{q}^{-\tau }.\end{array}$(A7) At last, combining (A3(A3) $\begin{array}{rl}& \mathrm{P}(A_2\ge {\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le e^{-C_K\log \left(pq\right)}\mathrm{E}\exp \left(t_1\right|{\omega }_k\left|\right)\\ & \le \exp \{-C_K\log (pq)+t_1^2\mathrm{E}{\omega }_k^2\exp (t_1\left|{\omega }_k\right|\left)\right\}\\ & \le \exp \{-C_K\log (pq)+{\eta }^{-1}K\log (pq\left)\right\}\\ & =p^{-\tau }{q}^{-\tau }.\end{array}$(A3) ) and (A7(A7) $\begin{array}{rl}& \mathrm{P}(A_1=\left|\right({\hat{\gamma }}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}|\ge 2{\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\\ & \le 2p^{-\tau }{q}^{-\tau }.\end{array}$(A7) ), we get $\mathrm{P}(A_1+A_2\ge 3{\eta }^{-1}{C}_K\sqrt{n\log \left(pq\right)})\le 3p^{-\tau }{q}^{-\tau },$which, coupled with (A2(A2) $\begin{array}{rl}& |{\hat{\rho }}_k^2-{\rho }_k^2|\\ & =\left|({\hat{\overrightarrow{\gamma }}}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}+{\overrightarrow{\gamma }}_k^{\mathrm{\top }}(\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}-\mathrm{E}\mathbf{y}{x}_k)\right|\\ & \le \left|({\hat{\overrightarrow{\gamma }}}_k-{\overrightarrow{\gamma }}_k{)}^{\mathrm{\top }}\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}\right|+\left|{\overrightarrow{\gamma }}_k^{\mathrm{\top }}(\frac{1}{n}\sum _{i=1}^n{\mathbf{y}}_i{x}_{ik}-\mathrm{E}\mathbf{y}{x}_k)\right|\\ \\ & \triangleq A_1+A_2.\end{array}$(A2) ), implies (A1(A1) $\begin{array}{rl}& \mathrm{P}(\underset{1\le k\le p}{max}|{\hat{\rho }}_k^2-{\rho }_k^2|\ge \vartheta )\\ & =\mathrm{P}(\underset{1\le k\le p}{max}\left|\frac{1}{n}\sum _{i=1}^n{\hat{\overrightarrow{\gamma }}}_k^{\mathrm{\top }}{\mathbf{y}}_i{x}_{ik}-\mathrm{E}{{\overrightarrow{\gamma }}_k}^{\mathrm{\top }}\mathbf{y}{x}_k\right|\ge \vartheta )\\ & <3p^{-\tau }{q}^{-\tau }.\end{array}$(A1) ). This completes the proof.

A.2. A lemma used to prove theorem 4.1

Lemma A.1

Suppose Assumptions 3.1, 4.1 and 4.2 hold. Let ${\mathrm{\Sigma }}_{x\left(M\right)}$ denote the submatrix of ${\mathrm{\Sigma }}_x$ after M. Let ${\hat{\mathrm{\Sigma }}}_{x\left(M\right)}$ and ${\hat{\mathrm{\Sigma }}}_x$ denote the corresponding estimators, respectively. Suppose $m=O(n^{2{\xi }_0}+4{\xi }_{min})$, where ${\xi }_0$ and ${\xi }_{min}$ are defined in Assumption 4.2. Then, with probability tending to one, we have (A8) ${\tau }_{min}\le \underset{\left|M\right|\le m}{min}{\lambda }_{min}\left({\hat{\mathrm{\Sigma }}}_{x\left(M\right)}\right)\le \underset{\left|M\right|\le m}{max}{\lambda }_{max}\left({\hat{\mathrm{\Sigma }}}_{x\left(M\right)}\right)\le {\tau }_{max}.$(A8)

Proof.

The proof is similar to that for Lemma 1 of H. Wang (2009). Here, we relax the normality assumption of $\mathbf{x}$ to the elliptically contoured distribution.

First, for $i=1,\dots ,n$, $j,k=1,\dots ,p$, let $U_i=(z_{ij}+z_{ik})/\sqrt{2(1+{\rho }_{jk})}$ and $V_i=(z_{ij}-z_{ik})/\sqrt{2(1-{\rho }_{jk})}$. By Assumption 4.1 and the additive property of elliptical contoured distribution (Fang et al., 2018), $(z_{ij},z_{ik})\sim {\mathrm{E}\mathrm{C}}_2(0,0,1,1,{\rho }_{jk},g)$, $U_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{E}\mathrm{C}(0,1,g)$ and $V_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{E}\mathrm{C}(0,1,g)$.

Second, observe that $\begin{array}{rl}\sum _{i=1}^n(z_{ij}{z}_{ik}-{\rho }_{jk})& =\frac{1}{4}[\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\\ & -\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}].\end{array}$Then, following Lemma A.3 of Bickel and Levina (2008), we have (A9) $\begin{array}{rl}& \mathrm{P}\left(\right|{\hat{\sigma }}_{ij}-{\sigma }_{ij}|\ge nv)\\ & =\mathrm{P}(\left|\sum _{i=1}^n\right(x_{ij}{x}_{ik}-{\sigma }_{jk}\left)\right|\ge nv)\\ & =\mathrm{P}(\left|\sum _{i=1}^n\right(z_{ij}{z}_{ik}-{\rho }_{jk}\left)\right|\ge \frac{nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}})\\ & =\mathrm{P}[|\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\\ & -\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & \le \mathrm{P}[\left|\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\right|\\ & +\left|\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}\right|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & \le \mathrm{P}[\left|\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\right|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & +\mathrm{P}[\left|\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}\right|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & =\mathrm{P}\{\left|\sum _{i=1}^n\right(U_i^2-1\left)\right|\ge \frac{2nv}{(1+{\rho }_{jk})({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\}\\ & +\mathrm{P}\{\left|\sum _{i=1}^n\right(V_i^2-1\left)\right|\ge \frac{2nv}{(1-{\rho }_{jk})({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\}.\end{array}$(A9) Further, for $i=1,\dots ,n$, let $W_i=U_i^2-1$ and $B_n^2=\sum _{i=1}^n\text{var}\left(W_i\right)$. Observe that by Jensen inequality there exist positive constants $c_1,\dots ,c_n$, $C_3$ and $C_4$ such that $\begin{array}{rl}& \left|\frac{\ln \{\mathrm{E}\exp (\zeta W_i\left)\right\}}{{\zeta }^2}\right|\le c_i^2\mathrm{f}\mathrm{o}\mathrm{r}\left|\zeta \right|<C_3\mathrm{a}\mathrm{n}\mathrm{d}\\ & {\overline{lim}}_{n\to \mathrm{\infty }}\frac{1}{B_n^2}\sum _{i=1}^n{c}_n^2\le C_4,\end{array}$satisfying condition (P) of Saulis and Statulevicius (1991). The same result holds when $W_i=V_i^2-1$. By Theorem 3.2 of Saulis and Statulevicius (1991), the first and second terms of (A9(A9) $\begin{array}{rl}& \mathrm{P}\left(\right|{\hat{\sigma }}_{ij}-{\sigma }_{ij}|\ge nv)\\ & =\mathrm{P}(\left|\sum _{i=1}^n\right(x_{ij}{x}_{ik}-{\sigma }_{jk}\left)\right|\ge nv)\\ & =\mathrm{P}(\left|\sum _{i=1}^n\right(z_{ij}{z}_{ik}-{\rho }_{jk}\left)\right|\ge \frac{nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}})\\ & =\mathrm{P}[|\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\\ & -\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & \le \mathrm{P}[\left|\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\right|\\ & +\left|\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}\right|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & \le \mathrm{P}[\left|\sum _{i=1}^n\left\{\right(z_{ij}+z_{ik}{)}^2-2(1+{\rho }_{jk})\}\right|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & +\mathrm{P}[\left|\sum _{i=1}^n\left\{\right(z_{ij}-z_{ik}{)}^2-2(1-{\rho }_{jk})\}\right|\ge \frac{4nv}{({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}]\\ & =\mathrm{P}\{\left|\sum _{i=1}^n\right(U_i^2-1\left)\right|\ge \frac{2nv}{(1+{\rho }_{jk})({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\}\\ & +\mathrm{P}\{\left|\sum _{i=1}^n\right(V_i^2-1\left)\right|\ge \frac{2nv}{(1-{\rho }_{jk})({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\}.\end{array}$(A9) ) are bounded, respectively, by $2\exp \{-\frac{2n{v}^2}{(1+{\rho }_{jk}{)}^2({\sigma }_{jj}{\sigma }_{kk})+2v(1+{\rho }_{jk}\left)\right({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\}$and $2\exp \{-\frac{2n{v}^2}{(1-{\rho }_{jk}{)}^2({\sigma }_{jj}{\sigma }_{kk})+2v(1-{\rho }_{jk}\left)\right({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\}.$Therefore, $\begin{array}{rl}& \mathrm{P}\left(\right|{\hat{\sigma }}_{ij}-{\sigma }_{ij}|\ge nv)\\ & \le 4max[\exp \left\{\frac{-2n{v}^2}{\begin{array}{c}(1+{\rho }_{jk}{)}^2({\sigma }_{jj}{\sigma }_{kk})\\ +2v(1+{\rho }_{jk})({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}\end{array}}\right\},\\ & \times \phantom{-2n{v}^2\begin{array}{c}(1+{\rho }_{jk}{)}^2({\sigma }_{jj}{\sigma }_{kk})\\ +2v(1+{\rho }_{jk})({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}\end{array}}\exp \left\{\frac{-2n{v}^2}{(1-{\rho }_{jk}{)}^2({\sigma }_{jj}{\sigma }_{kk})+2v(1-{\rho }_{jk}\left)\right({\sigma }_{jj}{\sigma }_{kk}{)}^{1/2}}\right\}].\end{array}$

This, together with ${\lambda }_{max}\left({\mathrm{\Sigma }}_x\right)<2^{-1}{\tau }_{max}$, implies that (A10) $\mathrm{P}\left(\right|{\hat{\sigma }}_{ij}-{\sigma }_{ij}|\ge nv)\le C_1\exp (-C_{2}n{v}^2)\mathrm{f}\mathrm{o}\mathrm{r}\left|v\right|\le \delta ,$(A10) where the positive constants $C_1$, $C_2$ and δ all depend on ${\tau }_{max}$ alone (Bickel & Levina, 2008, Lemma A.3).

The rest of the proof follows exactly the same as that for Lemma 1 of H. Wang (2009).

A.3. Proof of theorem 4.1

Our proof follows similar arguments as in the proof of Theorem 1 of H. Wang (2009).

Assume that no relevant predictor has been discovered in the first ℓ iterations, i.e., $S\not\subset S^{\left(\ell \right)}$. We evaluate the probability that at least one relevant will be identified in the $(\ell +1)$'s iteration or equivalently its complementary probability that the predictor selected by the $(\ell +1)$'s iteration is still an irrelevant one.

Let $X_{\left(S\right)}=({\mathbf{x}}_{1\left(S\right)},\dots ,{\mathbf{x}}_{n\left(S\right)})$ denote the subset of X corresponding to S. Let $B_{\left(S\right)}$ denote the coefficient matrix under the true model.

Denote $\begin{array}{rl}{H}_{\left(S^{\left(\ell \right)}\right)}& =X_{\left(S^{\left(\ell \right)}\right)}\{X_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{X}_{\left(S^{\left(\ell \right)}\right)}{\}}^{-1}{X}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }},\\ {\tilde{H}}_k^{\left(\ell \right)}& =X_{M_k^{\left(\ell \right)}}(X_{M_k^{\left(\ell \right)}}^{\mathrm{\top }}{X}_{M_k^{\left(\ell \right)}}{)}^{-1}{X}_{M_k^{\left(\ell \right)}}^{\mathrm{\top }},\\ {\mathbf{x}}_k^{\left(\ell \right)}& =(I_n-H_{S^{\left(\ell \right)}}){\mathbf{x}}_k,\\ H_k^{\left(\ell \right)}& ={\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^{-2}.\end{array}$Observe that (A11) $\begin{array}{rl}\mathrm{\Omega }\left(\ell \right)& =\mathrm{R}\mathrm{S}\mathrm{S}\left(S^{\left(\ell \right)}\right)-\mathrm{R}\mathrm{S}\mathrm{S}\left(S^{(\ell +1)}\right)\\ & =\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-{\tilde{H}}_k^{\left(\ell \right)}\left)Y\right\}-\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-{\tilde{H}}_k^{(\ell +1)}\left)Y\right\}\\ & =\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right({\tilde{H}}_k^{(\ell +1)}-{\tilde{H}}_k^{\left(\ell \right)}\left)Y\right\}\\ & =\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{a_{\ell +1}}^{(\ell {)}^{\mathrm{\top }}}{H}_{a_{\ell +1}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})Y\}.\end{array}$(A11) Assume $a_{\ell +1}\notin S$. We have (A12) $\begin{array}{rl}\mathrm{\Omega }\left(\ell \right)& \ge \underset{j\in S}{max}\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})Y\}\\ & \ge \mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})Y\},\end{array}$(A12) where $\begin{array}{rl}\hat{k}& ={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{k\in S}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}\\ & H_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}.\end{array}$Further, observe that the last inequality of (A12(A12) $\begin{array}{rl}\mathrm{\Omega }\left(\ell \right)& \ge \underset{j\in S}{max}\mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})Y\}\\ & \ge \mathrm{t}\mathrm{r}\left\{Y^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})Y\},\end{array}$(A12) ) is no less than (A13) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\}\\ & \ge \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}\\ & H_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\},\end{array}$(A13)

where $\mathit{ϵ}=({\overrightarrow{ϵ}}_1,\dots ,{\overrightarrow{ϵ}}_n{)}^{\mathrm{\top }}\in {\mathbb{R}}^{n\times q}$.

In what follows, we study the two terms in (A13(A13) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\}\\ & \ge \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}\\ & H_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\},\end{array}$(A13) ) separately.

Step 1: The first term of (A13(A13) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\}\\ & \ge \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}\\ & H_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\},\end{array}$(A13) ). Define $Q_{\left(S^{\left(k\right)}\right)}=I_n-H_{\left(S^{\left(k\right)}\right)}$. And denote ${\mathbf{x}}_k^{\left(\ell \right)}={\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}$. Then, the first term in (A13(A13) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\}\\ & \ge \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}\\ & H_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\},\end{array}$(A13) ) can be expressed as (A14) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\underset{j\in S}{max}\mathrm{t}\mathrm{r}\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}\frac{{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}}{\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^4}\\ & \times Q_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\phantom{\frac{4}{3}}\}\\ & =\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}\frac{{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}}{\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^2}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}.\end{array}$(A14) Define $\begin{array}{rl}{k}^{\ast }& ={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}_{k\in S}\mathrm{t}\mathrm{r}\left\{\right(X_{\left(S\right)}{B}_{\left(S\right)}{)}^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}\\ & Q_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}{B}_{\left(S\right)}\phantom{Q_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}}\}.\end{array}$Thus, the RHS of (A14(A14) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\underset{j\in S}{max}\mathrm{t}\mathrm{r}\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}\frac{{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}}{\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^4}\\ & \times Q_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\phantom{\frac{4}{3}}\}\\ & =\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}\frac{{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}}{\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^2}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}.\end{array}$(A14) ) is no less than (A15) $\begin{array}{rl}& \Vert {\mathbf{x}}_{k^{\ast }}^{\left(\ell \right)}{\Vert }^{-2}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_{k^{\ast }}^{\left(\ell \right)}{\mathbf{x}}_{k^{\ast }}^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge \underset{j\in S}{min}\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^{-2}\\ & \times \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_{k^{\ast }}^{\left(\ell \right)}{\mathbf{x}}_{k^{\ast }}^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =(\underset{k\in S}{max}\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^2{)}^{-1}\\ & \times \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge (\underset{k\in S}{max}\Vert {\mathbf{x}}_k{\Vert }^2{)}^{-1}\\ & \times \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\},\end{array}$(A15) where the last inequality follows after the fact $\Vert {\mathbf{x}}_k\Vert \ge \Vert {\mathbf{x}}_k^{\left(\ell \right)}\Vert$.

On the other hand, observe that (A16) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\sum _{j=1}^q\left(\sum _{k\in S}{\text{β}}_{kj}{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{\overrightarrow{\text{β}}}_j\right)\\ & \le \sum _{j=1}^q{\left(\sum _{k\in S}{\text{β}}_{kj}^2\right)}^{1/2}{\left\{\sum _{k\in S}{\left({\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{\overrightarrow{\text{β}}}_j\right)}^2\right\}}^{1/2}\\ & \le q\times \underset{j}{max}\Vert {\overrightarrow{\text{β}}}_j\Vert \mathrm{t}\mathrm{r}\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}\\ & {Q_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}}^{1/2}\times \sqrt{p_0}.\end{array}$(A16) Applying (A16(A16) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\sum _{j=1}^q\left(\sum _{k\in S}{\text{β}}_{kj}{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{\overrightarrow{\text{β}}}_j\right)\\ & \le \sum _{j=1}^q{\left(\sum _{k\in S}{\text{β}}_{kj}^2\right)}^{1/2}{\left\{\sum _{k\in S}{\left({\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{\overrightarrow{\text{β}}}_j\right)}^2\right\}}^{1/2}\\ & \le q\times \underset{j}{max}\Vert {\overrightarrow{\text{β}}}_j\Vert \mathrm{t}\mathrm{r}\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}\\ & {Q_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}}^{1/2}\times \sqrt{p_0}.\end{array}$(A16) ) to (A15(A15) $\begin{array}{rl}& \Vert {\mathbf{x}}_{k^{\ast }}^{\left(\ell \right)}{\Vert }^{-2}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_{k^{\ast }}^{\left(\ell \right)}{\mathbf{x}}_{k^{\ast }}^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge \underset{j\in S}{min}\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^{-2}\\ & \times \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_{k^{\ast }}^{\left(\ell \right)}{\mathbf{x}}_{k^{\ast }}^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =(\underset{k\in S}{max}\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^2{)}^{-1}\\ & \times \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge (\underset{k\in S}{max}\Vert {\mathbf{x}}_k{\Vert }^2{)}^{-1}\\ & \times \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\},\end{array}$(A15) ) and (A14(A14) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\underset{j\in S}{max}\mathrm{t}\mathrm{r}\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}\frac{{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}}{\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^4}\\ & \times Q_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\phantom{\frac{4}{3}}\}\\ & =\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}\frac{{\mathbf{x}}_k^{\left(\ell \right)}{\mathbf{x}}_k^{(\ell {)}^{\mathrm{\top }}}}{\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^2}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}.\end{array}$(A14) ), using the fact that $\underset{k\in S}{max}\Vert {\mathbf{x}}_k{\Vert }^2/n\le {\tau }_{max}$ with probability tending to one, and by Assumptions 3.1, 4.1 and Lemma 1 of H. Wang (2009), we obtain (A17) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge \frac{{\left[\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\right]}^2}{q^{2}n{\tau }_{max}{p}_0\times {max}_j\Vert {\overrightarrow{\text{β}}}_j{\Vert }^2}.\end{array}$(A17) Define ${\zeta }_{\left(S^{\left(\ell \right)}\right)}=(X_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{X}_{\left(S^{\left(\ell \right)}\right)}{)}^{-1}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}$. We have $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & =\mathrm{t}\mathrm{r}(B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}-{\zeta }_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{X}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{X}_{\left(S^{\left(\ell \right)}\right)}{\zeta }_{\left(S^{\left(\ell \right)}\right)}).\end{array}$Recall the assumption at the beginning of the proof that $S\not\subset S^{\left(k\right)}$. Then, by Assumptions 3.1, 3.2 and 4.1, and Lemma 1 of H. Wang (2009), we get (A18) $\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\ge qn{\tau }_{min}{\text{β}}_{min}^2$(A18) with probability tending to one. Applying (A18(A18) $\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\ge qn{\tau }_{min}{\text{β}}_{min}^2$(A18) ) to (A17(A17) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge \frac{{\left[\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\right]}^2}{q^{2}n{\tau }_{max}{p}_0\times {max}_j\Vert {\overrightarrow{\text{β}}}_j{\Vert }^2}.\end{array}$(A17) ) and using Assumptions 3.2, 4.1 and 4.2, we obtain (A19) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge n{\tau }_{max}^{-1}{p}_0^{-1}(\underset{j}{max}\Vert {\overrightarrow{\text{β}}}_j\Vert {)}^{-2}{\tau }_{min}^2{\text{β}}_{min}^4\\ & \ge {\tau }_{max}^{-1}{\nu }^{-1}{C}_B^{-2}{\tau }_{min}^2{\nu }_B^4{n}^{1-{\xi }_0-4{\xi }_{min}}.\end{array}$(A19) Step 2: The second term of (A13(A13) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\}\\ & \ge \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}\\ & H_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\},\end{array}$(A13) ). Recall ${\mathbf{x}}_k^{\left(\ell \right)}={\mathbf{x}}_k-H_{\left(S^{\left(\ell \right)}\right)}{\mathbf{x}}_k={\mathbf{x}}_k-X_{\left(S^{\left(\ell \right)}\right)}{\theta }_{k\left(S^{\left(\ell \right)}\right)}$where ${\theta }_{k\left(S^{\left(\ell \right)}\right)}={\left(X_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{X}_{\left(S^{\left(\ell \right)}\right)}\right)}^{-1}{X}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k.$Then, by Assumptions 3.1 and 3.2 and Lemma 1 of H. Wang (2009), we have $\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^2\ge n{\tau }_{min}$.

Moreover, since ${\mathbf{x}}_k^{\left(\ell \right)}=(I_n-H_{\left(S^{\left(\ell \right)}\right)}){\mathbf{x}}_k$, we have (A20) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{{\overrightarrow{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\overrightarrow{ϵ}\}\\ & =\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^{-2}\mathrm{t}\mathrm{r}({\overrightarrow{ϵ}}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}\overrightarrow{ϵ}-{\overrightarrow{ϵ}}^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}\overrightarrow{ϵ})\\ & \le {\tau }_{min}^{-1}{n}^{-1}\mathrm{t}\mathrm{r}({\overrightarrow{ϵ}}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}\overrightarrow{ϵ}-{\overrightarrow{ϵ}}^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}\overrightarrow{ϵ})\\ \\ & ={\tau }_{min}^{-1}{n}^{-1}\mathrm{t}\mathrm{r}\left({ϵ}^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}\overrightarrow{ϵ}\right)\\ & \le {\tau }_{min}^{-1}{n}^{-1}\underset{k\in S}{max}\underset{\left|M\right|\le m^{\ast }}{max}\mathrm{t}\mathrm{r}\left({\overrightarrow{ϵ}}^{\mathrm{\top }}{Q}_{\left(M\right)}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(M\right)}\overrightarrow{ϵ}\right)\\ & ={\tau }_{min}^{-1}{n}^{-1}\underset{k\in S}{max}\underset{\left|M\right|\le m^{\ast }}{max}\sum _{j=1}^q\left({\overrightarrow{ϵ}}_j^{\mathrm{\top }}{Q}_{\left(M\right)}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(M\right)}{\overrightarrow{ϵ}}_j\right),\end{array}$(A20) where $m^{\ast }=K\nu n^{2{\xi }_0+4{\xi }_{min}}$.

Notice that ${\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(M\right)}{\overrightarrow{ϵ}}_j$ is a normal random variable with mean zero and variance given by $\Vert Q_{\left(M\right)}{\mathbf{x}}_k{\Vert }^2\le \Vert {\mathbf{x}}_k{\Vert }^2$. Thus, the RHS of (A20(A20) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{{\overrightarrow{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\overrightarrow{ϵ}\}\\ & =\Vert {\mathbf{x}}_k^{\left(\ell \right)}{\Vert }^{-2}\mathrm{t}\mathrm{r}({\overrightarrow{ϵ}}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}\overrightarrow{ϵ}-{\overrightarrow{ϵ}}^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}\overrightarrow{ϵ})\\ & \le {\tau }_{min}^{-1}{n}^{-1}\mathrm{t}\mathrm{r}({\overrightarrow{ϵ}}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}\overrightarrow{ϵ}-{\overrightarrow{ϵ}}^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{H}_{\left(S^{\left(\ell \right)}\right)}\overrightarrow{ϵ})\\ \\ & ={\tau }_{min}^{-1}{n}^{-1}\mathrm{t}\mathrm{r}\left({ϵ}^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(S^{\left(\ell \right)}\right)}\overrightarrow{ϵ}\right)\\ & \le {\tau }_{min}^{-1}{n}^{-1}\underset{k\in S}{max}\underset{\left|M\right|\le m^{\ast }}{max}\mathrm{t}\mathrm{r}\left({\overrightarrow{ϵ}}^{\mathrm{\top }}{Q}_{\left(M\right)}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(M\right)}\overrightarrow{ϵ}\right)\\ & ={\tau }_{min}^{-1}{n}^{-1}\underset{k\in S}{max}\underset{\left|M\right|\le m^{\ast }}{max}\sum _{j=1}^q\left({\overrightarrow{ϵ}}_j^{\mathrm{\top }}{Q}_{\left(M\right)}{\mathbf{x}}_k{\mathbf{x}}_k^{\mathrm{\top }}{Q}_{\left(M\right)}{\overrightarrow{ϵ}}_j\right),\end{array}$(A20) ) is no greater than (A21) $q{\tau }_{min}^{-1}{n}^{-1}\underset{k\in S}{max}\Vert {\mathbf{x}}_k{\Vert }^2\underset{k\in S}{max}\underset{\left|M\right|\le m^{\ast }}{max}{\chi }_1^2,$(A21) where ${\chi }_1^2$ stands for a chi-square random variable with one degree of freedom. By Assumptions 3.1 and 4.1, and Lemma 1 of H. Wang (2009), we get that $n^{-1}\underset{k\in S}{max}\Vert {\mathbf{x}}_k{\Vert }^2\le {\tau }_{max}$ with probability tending to one.

On the other hand, the total number of combinations for $k\in S$ and $\left|M\right|\le m^{\ast }$ is no more than $d^{m^{\ast }+2}$. Then, by Assumption 4.2, we get $\begin{array}{rl}& \underset{k\in S}{max}\underset{\left|M\right|=\ell }{max}{\chi }_1^2\\ & \le 2(m^{\ast }+2)\log \left(p\right)\\ & \le 3K\nu n^{2{\xi }_0+4{\xi }_{min}}\times \nu n^{\xi }=3K{\nu }^2{n}^{\xi +2{\xi }_0+4{\xi }_{min}}\end{array}$with probability tending to one. Therefore, (A21(A21) $q{\tau }_{min}^{-1}{n}^{-1}\underset{k\in S}{max}\Vert {\mathbf{x}}_k{\Vert }^2\underset{k\in S}{max}\underset{\left|M\right|\le m^{\ast }}{max}{\chi }_1^2,$(A21) ) is bounded by (A22) $q{\tau }_{min}^{-1}{\tau }_{max}3K{\nu }^2{n}^{\xi +2{\xi }_0+4{\xi }_{min}-1}.$(A22) Combining (A13(A13) $\begin{array}{rl}& \mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_{\hat{k}}^{(\ell {)}^{\mathrm{\top }}}{H}_{\hat{k}}^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\}\\ & \ge \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}\\ & H_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})X_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\}\\ & -\underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{{\mathit{ϵ}}^{\mathrm{\top }}\right(I_n-H_{\left(S^{\left(\ell \right)}\right)}{)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}(I_n-H_{\left(S^{\left(\ell \right)}\right)})\mathit{ϵ}\},\end{array}$(A13) ), (A19(A19) $\begin{array}{rl}& \underset{k\in S}{max}\mathrm{t}\mathrm{r}\left\{B_{\left(S\right)}^{\mathrm{\top }}{X}_{\left(S\right)}{Q}_{\left(S^{\left(\ell \right)}\right)}^{\mathrm{\top }}{H}_k^{(\ell {)}^{\mathrm{\top }}}{H}_k^{\left(\ell \right)}{Q}_{\left(S^{\left(\ell \right)}\right)}{X}_{\left(S\right)}^{\mathrm{\top }}{B}_{\left(S\right)}\right\}\\ & \ge n{\tau }_{max}^{-1}{p}_0^{-1}(\underset{j}{max}\Vert {\overrightarrow{\text{β}}}_j\Vert {)}^{-2}{\tau }_{min}^2{\text{β}}_{min}^4\\ & \ge {\tau }_{max}^{-1}{\nu }^{-1}{C}_B^{-2}{\tau }_{min}^2{\nu }_B^4{n}^{1-{\xi }_0-4{\xi }_{min}}.\end{array}$(A19) ) and (A22(A22) $q{\tau }_{min}^{-1}{\tau }_{max}3K{\nu }^2{n}^{\xi +2{\xi }_0+4{\xi }_{min}-1}.$(A22) ), we have (A23) $\begin{array}{rl}{n}^{-1}\mathrm{\Omega }\left(\ell \right)& \ge {\tau }_{max}^{-1}{\nu }^{-1}{C}_B^{-2}{\tau }_{min}^2{\nu }_B^4{n}^{1-{\xi }_0-4{\xi }_{min}}\\ & -q{\tau }_{min}^{-1}{\tau }_{max}3K{\nu }^2{n}^{\xi +2{\xi }_0+4{\xi }_{min}-1}\\ & ={\tau }_{max}^{-1}{\nu }^{-1}{C}_B^{-2}{\tau }_{min}^2{\nu }_B^4{n}^{1-{\xi }_0-4{\xi }_{min}}\\ & \times (1-q{\tau }_{max}^2{\nu }^3{C}_B^2{\tau }_{max}^{-3}{\nu }_B^{-4}3K{n}^{\xi +3{\xi }_0+8{\xi }_{min}-1})\end{array}$(A23) uniformly for every $\ell \le K{n}^{{\xi }_0+4{\xi }_{min}}$. Recall $K=2{\tau }_{max}\nu C_B^2{\tau }_{min}^{-2}{\nu }_B^{-4}$ defined in Section 4. Then, by Assumption 4.2 and 4.3, we have (A24) $\begin{array}{rl}{n}^{-1}\Vert Y{\Vert }_F^2& \ge n^{-1}\sum _{\ell =1}^{K{n}^{{\xi }_0}+4{\xi }_{min}}\mathrm{\Omega }\left(\ell \right)\\ & \ge 2(1-q{\tau }_{max}^2{\nu }^3{C}_B^2{\tau }_{max}^{-3}{\nu }_B^{-4}3K{n}^{\xi +3{\xi }_0+8{\xi }_{min}-1})\\ & \stackrel{p}{\to }2,\end{array}$(A24) where $\Vert \cdot {\Vert }_F$ is the Frobenius norm.

Without loss of generality, we can assume $\text{var}\left(y_{i1}\right)+\cdots +\text{var}\left(y_{iq}\right)=1$, and we have $n^{-1}\Vert Y{\Vert }_F^2\stackrel{p}{\to }1$. (Otherwise, we can standardize it by letting $y_{ik}^{\ast }=y_{ik}/\sqrt{\text{var}\left(y_{i1}\right)+\cdots +\text{var}\left(y_{iq}\right)}$ for $i=1,\dots ,n$ and $k=1,\dots ,q$.) Thus, it is impossible to have $S^{\left(k\right)}\bigcup M_t=\mathrm{\varnothing }$ for every $1\le k\le K{n}^{{\xi }_0+4{\xi }_{min}}$, which implies that at least one relevant variable will be discovered within $K{n}^{{\xi }_0+4{\xi }_{min}}$ steps. This completes the proof.