A selective overview of sparse sufficient dimension reduction

Abstract

High-dimensional data analysis has been a challenging issue in statistics. Sufficient dimension reduction aims to reduce the dimension of the predictors by replacing the original predictors with a minimal set of their linear combinations without loss of information. However, the estimated linear combinations generally consist of all of the variables, making it difficult to interpret. To circumvent this difficulty, sparse sufficient dimension reduction methods were proposed to conduct model-free variable selection or screening within the framework of sufficient dimension reduction. We review the current literature of sparse sufficient dimension reduction and do some further investigation in this paper.

Keywords:

• Minimax rate
• sparse sufficient dimension reduction
• variable selection
• variable screening

1. Introduction

The rapid development of data collection technology in areas, such as biology, financial econometrics and signal processing, has posed a great challenge for traditional multivariate analysis. High-dimensional data analysis becomes ubiquitous and increasingly important. Dimension reduction, and in particular sufficient dimension reduction for regression, offers an appealing avenue to tackle high-dimensional problems. It is often desirable to reduce the dimensionality of the problem by replacing the original high-dimensional data with a low-dimensional space composed of a few linear combinations of predictors, which are usually much smaller than the original dimension. Although sufficient dimension reduction is an effective way to extract relevant information from high-dimensional data sets, while grasping the important features or patterns in the data, the linear combinations usually consist of all original predictors which makes the interpretation difficult. This limitation can be overcome via variable selection, where a subset of relevant predictor variables is selected. The removal of the excess variables not only can reduce the noise to the precise estimation, alleviate the collinearity issue, but also help reduce the computational cost caused by high-dimensional data.

As one of the most important dimension reduction approaches, many variable selection methods have been developed. Some most popular variable selection approaches are developed under the linear model or the generalised linear model paradigm, such as nonnegative garrotte (Breiman, 1995), the least absolute shrinkage and selection operator (Lasso, hereafter) (Tibshirani, 1996), the smoothly clipped absolute deviation (SCAD, hereafter) (Fan & Li, 2001), adaptive Lasso (Zou, 2006), group Lasso (Yuan & Lin, 2006), Dantzig selector (Candes & Tao, 2007) and the minimax concave plus penalty (MCP, hereafter) (Zhang, 2010).

These model-based variable selection methods assume the underlying true model is known up to a finite dimensional parameter or the imposed working model is usefully similar to the true model. However, the true model might be in a complex form and it is usually unknown. If the underlying modelling assumption is violated, these variable selection methods might fail. Hence, model-free variable selection method, which does not require the full knowledge of the underlying true model, is called for. It has been shown that the general framework of sufficient dimension reduction is useful for variable selection (Bondell & Li, 2009) since no pre-specified underlying models between the response and the predictors are required. So model-free variable selection can be achieved through the framework of SDR (Cook, 1998; Li, 1991, 2000).

Let $\mathbf{X}=(X_1,\dots ,X_p{)}^{\mathrm{\top }}$ be the predictor and $Y$ be the scalar response. The goal of variable selection is to seek the smallest subset of the predictors ${\mathbf{X}}_{\mathcal{A}}$, with partition $\mathbf{X}=({\mathbf{X}}_{\mathcal{A}}^{\mathrm{\top }},{\mathbf{X}}_{{\mathcal{A}}^c}^{\mathrm{\top }}{)}^{\mathrm{\top }}$, such that (1) $Y\text{boxhU}{\mathbf{X}}_{{\mathcal{A}}^c}|{\mathbf{X}}_{\mathcal{A}}.$(1) Here $\mathcal{A}$ denotes a subset of indices of $\{1,\dots ,p\}$ corresponding to the relevant predictor set ${\mathbf{X}}_{\mathcal{A}}$, and ${\mathcal{A}}^c$ is the complement of $\mathcal{A}$, i.e., ${\mathbf{X}}_{\mathcal{A}}=\{x_i:i\in \mathcal{A}\}$ and ${\mathbf{X}}_{{\mathcal{A}}^c}=\{x_i:i\in {\mathcal{A}}^c\}$. Condition (1(1) $Y\text{boxhU}{\mathbf{X}}_{{\mathcal{A}}^c}|{\mathbf{X}}_{\mathcal{A}}.$(1) ) implies that ${\mathbf{X}}_{\mathcal{A}}$ contains all the active predictors in terms of predicting $Y$. The existence and uniqueness of $\mathcal{A}$ were discussed in details in Yin and Hilafu (2015). Ideally, we want to find the smallest index set $\mathcal{A}$ satisfying (1(1) $Y\text{boxhU}{\mathbf{X}}_{{\mathcal{A}}^c}|{\mathbf{X}}_{\mathcal{A}}.$(1) ), in which case no inactive predictors are included in ${\mathbf{X}}_{{\mathcal{A}}^c}$.

Model-free variable selection is closely related to sufficient dimension reduction, which aims to find $\text{β}\in {\mathbb{R}}^{p\times d}$ with $d\le p$, such that (2) $Y\text{boxhU}\mathbf{X}|{\text{β}}^{\mathrm{\top }}\mathbf{X},$(2) that is, $Y$ is independent of $\mathbf{X}$ conditioning on ${\text{β}}^{\mathrm{\top }}\mathbf{X}$. The column space of such $\text{β}$, $\mathbf{S}\mathbf{p}\mathbf{a}\mathbf{n}\left(\text{β}\right)$, is called a dimension reduction space. Under mild assumptions, such as given in Cook (1996) and Yin et al. (2008), the intersection of all such spaces is itself a dimension reduction space. In this case, we call the intersection the central subspace for the regression of $Y$ on $\mathbf{X}$, and denote it by ${\mathcal{S}}_{Y|\mathbf{X}}$. And its dimension, $d=\mathrm{d}\mathrm{i}\mathrm{m}\left({\mathcal{S}}_{Y|\mathbf{X}}\right)$, is usually much smaller than p, the dimension of the original predictor. Following the partition of $\mathbf{X}$, we can partition $\text{β}$ accordingly as $\text{β}=\left(\begin{array}{c}{\text{β}}_{\mathcal{A}}\\ {\text{β}}_{{\mathcal{A}}^c}\end{array}\right),{\text{β}}_{\mathcal{A}}\in {\mathbb{R}}^{|\mathcal{A}|\times d},{\text{β}}_{{\mathcal{A}}^c}\in {\mathbb{R}}^{(p-|\mathcal{A}|)\times d},$ where $|\mathcal{A}|$ is the cardinality of $\mathcal{A}$. Hence, (1(1) $Y\text{boxhU}{\mathbf{X}}_{{\mathcal{A}}^c}|{\mathbf{X}}_{\mathcal{A}}.$(1) ) is equivalent to ${\text{β}}_{{\mathcal{A}}^c}=0$.

Many methods have been proposed for estimating the basis of ${\mathcal{S}}_{Y|\mathbf{X}}$ in the literature, including sliced inverse regression (SIR, hereafter) (Li, 1991), sliced average variance estimation (SAVE, hereafter) (Cook & Weisberg, 1991), principal Hessian directions (PHD, hereafter) (Li, 1992), minimum average variance estimation (MAVE, hereafter) (Xia et al., 2002), directional regression (DR, hereafter) (Li & Wang, 2007), principal fitted component (PFC, hereafter) (Cook & Forzani, 2008), semiparametric approach (Ma & Zhu, 2012), etc. Several methods have been also suggested for simultaneously selecting the contributing predictors. These include shrinkage SIR (Ni et al., 2005), sparse SIR (Li, 2007; Li & Nachtsheim, 2006), sparse SAVE and sparse PHD (Li, 2007), constrained canonical correlation (Zhou & He, 2008), the general shrinkage strategy for inverse regression estimation (Bondell & Li, 2009), the regularised SIR estimator with SCAD penalty (Wu & Li, 2011) and coordinate independent sparse estimation (CISE, hereafter) (Chen et al., 2010), conditional covariance minimisation (Chen et al., 2017), etc.

Although these aforementioned methods can select the significant predictors without assuming an underlying parametric model, they are not designed for $p\gg n$ problems, in which the number of predictor variables is larger than the number of observations. The so-called large p small n problems are increasingly common with rapid technological advances in data collection and have attracted a lot of research interests. We hereby give a very brief review of model-free variable selections via sufficient dimension reduction approach under the $p\gg n$ setting. Li and Yin (2008) proposed sparse ridge SIR, which combined SIR with both ${\ell }_1$- and ${\ell }_2$-regularisation to achieve dimension reduction and variable selection simultaneously, even when p>n. Yu et al. (2013) suggested combining SIR with the Dantzig selector (Candes & Tao, 2007) to recover the central subspace in the general semiparametric models. A non-asymptotic error bound for the resulting estimator is derived and the error bound is of order $O_p\left(\right(\log p/n{)}^{1/2})$, which appears to be optimal. Moreover, they proposed another regularised version of SIR with the adaptive Dantzig selector. The resulting estimators defined from variable selection are asymptotically normal even when the predictor dimension diverges to infinity. It is worth mentioning that the $|\mathcal{A}|$ is fixed in Yu et al. (2013). Yu, Dong, Zhu (2016) proposed trace pursuit for model-free variable selection under the sufficient dimension reduction paradigm. Two distinct algorithms are proposed: stepwise trace pursuit (STP, hereafter) and forward trace pursuit (FTP, hereafter). Stepwise trace pursuit achieved selection consistency with fixed p and is applicable in the setting with p>n. Furthermore, forward trace pursuit can serve as an initial screening step to speed up the computation in the case of ultrahigh dimensionality. Li and Dong (2020) extended trace pursuit method to matrix-valued predictors based on Yu, Dong, Zhu (2016). To test the importance of rows, columns and submatrices of the predictor matrix in terms of predicting the response, three types of hypotheses are formulated under a unified framework. The asymptotic properties of the test statistics under the null hypothesis are established and a permutation testing algorithm is also introduced to approximate the distribution of the test statistics. Tan et al. (2018) developed a convex formulation for fitting sparse SIR in high dimensions. They solved the resulting convex optimisation problem via the linearised alternating direction methods of multiple algorithms and established an upper bound on the subspace distance between the estimated and the true subspaces. Unlike Yu et al. (2013), Lin et al. (2019) allowed $|\mathcal{A}|$ goes to infinity. By constructing artificial response variables made up from top eigenvectors of the estimated conditional covariance matrix, Lin et al. (2019) introduced a simple Lasso regression method to obtain an estimator of the sufficient dimension reduction space. The resulting algorithm, Lasso-SIR, is shown to be consistent and achieves the optimal convergence rate under certain sparsity conditions when p is of order $o\left(n^2{c}^2\right)$, where c is the generalised signal-to-noise ratio, which is only the first step of Tan et al. (2020). Moreover, Tan et al. (2020) discovered the possible trade-off between statistical guarantee and computational performance for sparse SIR and proposed an adaptive estimation scheme for sparse SIR which is computationally tractable and rate optimal under the condition that $\log p=o\left(n\right)$, which is weaker than Lin et al. (2019).

There is considerable literature on applying sufficient dimension reduction for model-free selection, but the study of developing screening consistency for the ultra-high dimensional setting is still lacking. To fulfil the aforementioned gaps, Zhu et al. (2011) proposed a variable screening procedure under a unified model framework, which contains a wide variety of commonly used parametric and semiparametric models. The new method does not require imposing a specific model structure on regression functions and thus is particularly appealing to ultrahigh-dimensional regressions. They also showed the proposed method achieves screening consistency even with the number of predictors growing at an exponential rate of the sample size. Yu, Dong, Shao (2016) proposed an approach called marginal SIR for model-free variable selection. Furthermore, marginal SIR with Dantzig selector exploits the sparsity structure in the marginal utility and achieves the desirable selection consistency property. Lin et al. (2017) first introduced a large class of models depending on the smallest non-zero eigenvalue of the kernel matrix of SIR, then the minimax rate for estimating the central space is derived, which is the first paper studied the minimax estimation of sparse SIR. However, they only considered the projection loss (Li & Wang, 2007). More importantly, their theoretical study is based on the assumption that the covariance matrix is diagonal. As far as we know, most of mentioned work mainly focus on SIR with consistency on variable selection. Qian et al. (2019) provided simultaneous analysis for PFC and SAVE. Furthermore, their approach allows many quantities such as the structural dimension, the number of important predictors and the number of slices to diverge with n. To deliver the most essential messages, in the following section, we focus our discussion on the papers mentioned above.

2. Review of sufficient dimension reduction

Sufficient dimension reduction aims to find the column space of $\text{β}$ with the smallest dimension d. In other words, sufficient dimension reduction is proposed as a problem of estimating a space, instead of the classic statistical problem of estimating parameters. As mentioned in the introduction, there are many approaches in the literature of sufficient dimension reduction for estimating the column space $\text{β}$: sliced inverse regression (SIR; Li, 1991), sliced average variance estimation (SAVE; Cook & Weisberg, 1991), minimum average variance estimation (MAVE; Xia et al., 2002), the kth moment estimation (Yin & Cook, 2002, 2003), inverse regression (Cook & Ni, 2005), directional regression (DR; Li & Wang, 2007), sliced regression (SR; Wang & Xia, 2008), likelihood acquired directions (LAD; Cook & Forzani, 2009), semiparametric approaches (Ma & Zhu, 2012, 2013a, 2013b, 2014), etc. We mainly review three inverse regression-based methods (SIR; SAVE and DR) for estimating ${\mathcal{S}}_{Y|\mathbf{X}}$ for our subsequent investigation.

Inverse regression methods constitute the oldest class of dimension reduction methods and are still under active development currently. The main idea of the inverse regression is to reverse the relation between the response and the predictors (Li, 1991). Instead of considering distributions or expectations of functions of Y conditional on $\mathbf{X}$, which suffers the curse of dimensionality when $\mathbf{X}$ is high dimensional, these inverse regression-based methods consider expectations of functions of $\mathbf{X}$ conditional on Y, which is suddenly a low dimensional problem because Y is univariate. The inverse regression-based methods often are based on some additional assumptions on the predictors to link the low dimensional problem and the original high dimensional problem. These additional assumptions are given as follows.

(W1)	linearity condition $E\left(\mathbf{X}|{\text{β}}^{\mathrm{\top }}\mathbf{X}\right)=\mathbf{\Sigma }\text{β}({\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }\text{β}{)}^{-1}{\text{β}}^{\mathrm{\top }}\mathbf{X}$.
(W2)	constant variance condition $\mathrm{c}\mathrm{o}\mathrm{v}\left(\mathbf{X}|{\text{β}}^{\mathrm{\top }}\mathbf{X}\right)=\mathbf{\Sigma }-\mathbf{\Sigma }\text{β}\left({\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }\text{β}{)}^{-1}{\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }\text{β}\right({\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }\text{β}{)}^{-1}{\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }$,

where $\mathbf{\Sigma }=\mathrm{c}\mathrm{o}\mathrm{v}\left(\mathbf{X}\right)$. As is known to all, SIR only requires the condition (W1) holds. However, SAVE and DR need both conditions.

When the linearity condition and the constant variance condition are satisfied, the inverse regression methods formulate the problem of estimating ${\mathcal{S}}_{Y|\mathbf{X}}$ into an eigen-decomposition problem. Let $\mathbf{M}$ be the kernel matrix of a specific inverse regression based dimension reduction method. For the sufficient dimension reduction methods that aim to estimate ${\mathcal{S}}_{Y|\mathbf{X}}$, the kernel matrices corresponding to the three most well-known inverse regression methods are summarised as below: $\begin{array}{rl}& \text{SIR: }{\mathbf{M}}_{\mathrm{S}\mathrm{I}\mathrm{R}}=\mathrm{v}\mathrm{a}\mathrm{r}\left\{E\right(\mathbf{X}|Y\left)\right\};\\ & \text{SAVE: }{\mathbf{M}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}=E\{\mathbf{\Sigma }-\mathrm{v}\mathrm{a}\mathrm{r}(\mathbf{X}|Y){\}}^2;\\ & \text{DR: }{\mathbf{M}}_{\mathrm{D}\mathrm{R}}\\ & =2E^2\left\{E\right(\mathbf{X}|Y\left)E\right({\mathbf{X}}^{\mathrm{\top }}|Y\left)\right\}\\ & +2E\left\{E\right({\mathbf{X}}^{\mathrm{\top }}|Y\left)E\right(\mathbf{X}|Y\left)\right\}E\left\{E\right(\mathbf{X}|Y\left)E\right({\mathbf{X}}^{\mathrm{\top }}|Y\left)\right\}\\ & +2E\left\{E^2\right(\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}\left)\right\}-2\mathbf{\Sigma }.\end{array}$ Assuming $d={\mathcal{S}}_{Y|\mathbf{X}}$ is known, the procedure for a generalised eigenvalue-decomposition of the kernel matrix $\mathbf{M}$, that is $\begin{array}{rl}\mathbf{M}{\text{β}}_i& ={\lambda }_i\mathbf{\Sigma }{\text{β}}_i,\text{with }{\text{β}}_i^{\mathrm{\top }}\mathbf{\Sigma }{\text{β}}_j=1\text{if }i=j,\\ & {\text{β}}_i^{\mathrm{\top }}\mathbf{\Sigma }{\text{β}}_j=0\text{else }i\ne j,\end{array}$ where $i=1,\dots ,p$, and ${\lambda }_1\ge \cdots \ge {\lambda }_d>0={\lambda }_{d+1}=\cdots ={\lambda }_p$ are the eigenvalues. Then the eigenvectors corresponding to the nonzero eigenvalues $\text{β}=({\text{β}}_1,\dots ,{\text{β}}_d)$ form a basis of ${\mathcal{S}}_{Y|\mathbf{X}}$. Thus the sufficient dimension reduction directions $\text{β}$ can also be identified through the following optimisation problem (Tan et al., 2020): (3) $\widehat{\text{β}}=\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}\mathbf{M}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\mathbf{\Sigma }\mathbf{B}={\mathbf{I}}_d.$(3)

3. The current literature of variable selection via sufficient dimension reduction

3.1. Oracle property under the setting p<n

In the general framework of condition (1(1) $Y\text{boxhU}{\mathbf{X}}_{{\mathcal{A}}^c}|{\mathbf{X}}_{\mathcal{A}}.$(1) ), the shrinkage SIR method is developed in Ni et al. (2005) by applying the Lasso approach to SIR. When a subset of predictors are irrelevant, then the corresponding row estimates of $\text{β}$ is equal to 0, and consequently to achieve variable selection. Let $\alpha =({\alpha }_1,\dots ,{\alpha }_p{)}^{\mathrm{\top }}$, with ${\alpha }_i\in \mathbb{R}$, $i=1,\dots ,p$, be the shrinkage vector. Then based on the expression (3(3) $\widehat{\text{β}}=\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}\mathbf{M}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\mathbf{\Sigma }\mathbf{B}={\mathbf{I}}_d.$(3) ), the estimation of the shrinkage vector can be rewritten to minimise the following function over α (Ni et al., 2005): (4) $\begin{array}{rl}\hat{\alpha }& =\underset{\alpha }{\arg min}n\left(\mathrm{V}\mathrm{e}\mathrm{c}\right(\widehat{\text{β}})\\ & -\mathrm{V}\mathrm{e}\mathrm{c}\left\{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right(\alpha \left)\widehat{\text{β}}\right\}{)}^{\mathrm{\top }}\hat{\mathbf{M}}\left(\mathrm{V}\mathrm{e}\mathrm{c}\right(\widehat{\text{β}})-\mathrm{V}\mathrm{e}\mathrm{c}\{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\alpha \right)\widehat{\text{β}}\left\}\right),\\ & \text{subject to }\sum _{i=1}^p|{\alpha }_i|\le t,t\ge 0,\end{array}$(4) where $\hat{\mathbf{M}}$ is the estimator of kernel matrix $\mathbf{M}$. To investigate the asymptotic behaviour, we consider the Lagrangian formulation of the constrained optimisation problem. Specially, the optimisation problem in expression (4(4) $\begin{array}{rl}\hat{\alpha }& =\underset{\alpha }{\arg min}n\left(\mathrm{V}\mathrm{e}\mathrm{c}\right(\widehat{\text{β}})\\ & -\mathrm{V}\mathrm{e}\mathrm{c}\left\{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right(\alpha \left)\widehat{\text{β}}\right\}{)}^{\mathrm{\top }}\hat{\mathbf{M}}\left(\mathrm{V}\mathrm{e}\mathrm{c}\right(\widehat{\text{β}})-\mathrm{V}\mathrm{e}\mathrm{c}\{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\alpha \right)\widehat{\text{β}}\left\}\right),\\ & \text{subject to }\sum _{i=1}^p|{\alpha }_i|\le t,t\ge 0,\end{array}$(4) ) can be reformulated as $\hat{\alpha }=\arg \underset{\alpha }{min}(||\mathbf{U}-\mathbf{W}\alpha |{|}^2+{\tau }_n\sum _{i=1}^p|{\alpha }_i|),$ for some non-negative penalty constant ${\tau }_n$. In which, $\mathbf{U}=n^{1/2}{\hat{\mathbf{M}}}^{1/2}\mathrm{V}\mathrm{e}\mathrm{c}\left(\widehat{\text{β}}\right),\mathbf{W}=n^{1/2}{\hat{\mathbf{M}}}^{1/2}\widehat{\text{β}}.$ Then the central dimension reduction subspace ${\mathcal{S}}_{Y|\mathbf{X}}$ is estimated by $\mathbf{S}\mathbf{p}\mathbf{a}\mathbf{n}\left\{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\right(\hat{\mathit{\alpha }}\left)\widehat{\text{β}}\right\}$. Li (2007) extended shrinkage SIR method to SAVE and PHD methods, where the central dimension subspace is estimated the same as Ni et al. (2005), and $\widehat{\text{β}}$ corresponds to the estimated central dimension reduction directions of SAVE and PHD methods, respectively. Bondell and Li (2009) proposed a general shrinkage estimation strategy for the entire inverse regression estimation family that is capable of simultaneous sufficient dimension reduction and variable selection. They considered the adaptive Lasso, $\hat{\alpha }=\arg \underset{\alpha }{min}\left(||\mathbf{U}-\mathbf{W}\alpha |{|}^2+{\tau }_n\sum _{i=1}^p{w}_i|{\alpha }_i|\right),$ where $\mathbf{w}=(w_1,\dots ,w_p{)}^{\mathrm{\top }}$ is a known weights vector. They also demonstrated that the proposed class of shrinkage estimators has the desirable oracle property of consistency in variable selection while retaining root n estimation consistency.

However, most existing sparse dimension reduction methods mentioned above are conducted stepwise, estimating a sparse solution for a basis matrix of the central subspace column by column. Instead, Chen et al. (2010) proposed a unified one-step approach to reduce the number of variables appearing in the estimate of ${\mathcal{S}}_{Y|\mathbf{X}}$. Their approach, which depends operationally on Grassmann manifold optimisation, can achieve dimension reduction and variable selection simultaneously. Additionally, their 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. More importantly, Chen et al. (2010) is an extension to Bondell and Li (2009), which combined SIR, SAVE, DR with adaptive Lasso to variable selection. 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, followed by a simple variable selection method. Using the B-spline basis functions generated for the response variable, the CANCOR method (Fung et al., 2002) is asymptotically equivalent to SIR. Suppose that the range of $Y$ is a bounded interval $[a,b]$, given $k_n$ interval knots in $[a,b]$ and the spline order m, we generate $m+k_n$ B-spline basis functions. Under the linearity condition, CANCOR estimates a set of effective dimension reduction directions by estimating the canonical variates between the B-spline basis functions and $\mathbf{X}$. Since the generated $m+k_n$ B-spline basis functions add to 1, we use in CANCOR the first $m+k_n-1$ basis function of $Y$, $\pi \left(Y\right)=\left({\pi }_1\right(Y),\dots ,{\pi }_{m+k_n-1}(Y){)}^{\mathrm{\top }}$. Let $\mathcal{X}=({\mathbf{X}}_1,\dots ,{\mathbf{X}}_n{)}^{\mathrm{\top }}$ and ${\mathrm{\Pi }}_{n\times (m+k_n-1)}=\left(\pi \right(Y_1),\dots ,\pi (Y_n){)}^{\mathrm{\top }}$ be the data matrices containing the predictor values and the B-spline basis function values. Then the CANCOR method is to estimate the canonical correlations between the columns of $\mathcal{X}$ and the columns of Π. The dimensionality of the central dimension reduction subspace is selected by performing the following sequential tests on the number of the non-zero canonical correlations, $H_0:{\gamma }_s>{\gamma }_{s+1}=0$ versus $H_1:{\gamma }_{s+1}>0$ for $s=0,1,\dots ,p-1$, where ${\gamma }_s$ are the asymptotic canonical correlations between $\pi \left(Y\right)$ and $\mathbf{X}$ in decreasing order. The dimensionality estimate for d is the smallest s such that $H_0$ is not rejected. The CANCOR method actually solves an optimisation problem that sequentially finds the directions $\text{β}$ with the maximum correlation between ${\text{β}}^{\mathrm{\top }}\mathbf{X}$ and some functions of $Y$. Their procedure is attractive because they demonstrated that it also has the oracle property.

Sparse sufficient dimension reduction methods mentioned above focus on the cases when p is fixed. For regressions with diverging p, estimation and variable selection methods are also developed in the framework of sufficient dimension reduction: Zhu et al. (2006) studied the asymptotic properties of SIR as p diverges, but their result is for SIR only, and variable selection is not studied at all. Zhu and Zhu (2009a) investigated weighted partial least squares with a diverging p, but again variable selection is not derived. Zhu and Zhu (2009b) investigated variable selection with a diverging number of predictors through inverse regression, but focused on single-index models only. By contrast, Wu and Li (2011) established asymptotic properties for a family of inverse regression estimators that includes SIR, studied simultaneous dimension reduction and variable selection with a particular emphasis on the latter and encompassed more general forms, while the number of predictors p is allowed to diverge as the sample size n approaches infinity. Wu and Li (2011) adopted the SCAD type penalty that was first introduced by Fan and Li (2001), and combined it with sufficient dimension reduction estimator, that is $\hat{\alpha }=\arg \underset{\alpha }{min}\left(||\mathbf{U}-\mathbf{W}\alpha |{|}^2+\sum _{i=1}^p{p}_{{\tau }_n}\left(|{\alpha }_i|\right)\right).$ The penalty $p_{{\tau }_n}(\cdot )$ are not necessarily the same for all i. Wu and Li (2011) also showed that the penalised estimator selects all truly contributing predictors and excludes all irrelevant ones with probability approaching one.

Based on the work in kernel dimension reduction, Chen et al. (2017) proposed a method to perform feature selection via a constrained optimisation problem. The corresponding SDR method can refer to Fukumizu et al. (2009); Fukumizu Leng (2014). Many previous kernel approaches are filter methods based on the Hilbert–Schmidt Independence Criterion (HSIC, Gretton et al., 2005). Chen et al. (2017) proposed to use the trace of the conditional covariance operator as a criterion for feature selection. Let $(H_1,k_1)$ denote an RKHS supported on $\mathbf{X}\subset {\mathbb{R}}^p$. Then the trace of the conditional covariance operator, $\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left({\mathbf{\Sigma }}_{YY|\mathbf{X}}\right)$ can be interpreted as a dependence measure, as long as the $H_1$ is large enough. Then the problem of supervised feature selection reduces to minimising the trace of the conditional covariance operator over subsets of features with controlled cardinality: $\underset{T:|T|=d}{min}Q\left(T\right):=\mathrm{T}\mathrm{r}\left({\mathbf{\Sigma }}_{YY|{\mathbf{X}}_T}\right).$ They also showed that empirical estimate of the criterion is consistent as the sample size increases. It is worth noting that kernel feature selection methods have the advantage of capturing nonlinear relationships between the features and the labels.

Theorem 3.1

Assume $n^{1/2}\left\{\mathrm{V}\mathrm{e}\mathrm{c}\right(\widehat{\text{β})}-\mathrm{V}\mathrm{e}\mathrm{c}\left(\text{β}\right)\}\to \mathcal{N}(0,\mathrm{\Gamma }),$ for some $\mathrm{\Gamma }>0,$ and that ${\mathbf{M}}_n^{1/2}={\mathbf{M}}^{1/2}+o\left(\frac{1}{\sqrt{n}}\right)$. Suppose that ${\tau }_n\to \mathrm{\infty }$ and $\frac{{\tau }_n}{\sqrt{n}}\to 0$ with $p<n,$ then the shrinkage estimator $\widehat{\text{β}}$ satisfies

1. consistency in variable selection, $Pr(\hat{\mathcal{A}}=\mathcal{A})\to 1,$ and

2. asymptotic normality, $n^{1/2}\left\{\mathrm{V}\mathrm{e}\mathrm{c}\right({\widehat{\text{β}}}_{\mathcal{A}})-\mathrm{V}\mathrm{e}\mathrm{c}({\text{β}}_{\mathcal{A}}\left)\right\}\to \mathcal{N}(0,\mathrm{\Lambda }),$ for some $\mathrm{\Lambda }>0$.

Remark 3.1

Theorem 3.1, part (a), indicates that the sparse sufficient dimension reduction estimator can select contributing predictors consistently, i.e., for all $i\notin \mathcal{A}$ we have $Pr({\hat{\alpha }}_i\ne 0)\to 0$, and for all $i\in \mathcal{A}$ we have $Pr({\hat{\alpha }}_i\ne 0)\to 1$. Theorem 3.1, part (b), further shows that the estimator for ${\text{β}}_{\mathcal{A}}$ that corresponds to the contributing predictors is root n consistent. The oracle property as shown in Theorem 3.1 is given in Bondell and Li (2009), Chen et al. (2010), Wu and Li (2011) and Zhou and He (2008). Most of the methods mentioned above cannot achieve the desired property with p>n, however, Wu and Li (2011) showed that their proposed method can obtain selection consistency when p diverge as the sample size n goes to infinity. Then we turn to investigate the oracle property with p>n.

3.2. Oracle property under the setting $p\gg n$

Large-p-small-n problems appear frequently in fields such as biology, economics and finance. While those variable selection methods have been successfully applied in many high-dimensional analyses, modern applications in areas such as genomics and high-frequency finance further push the dimensionality of data to an even larger scale, where p may grow exponentially with n. Such ultrahigh-dimensional data present simultaneous challenges of computational expediency, statistical accuracy and algorithm stability. It is difficult to directly apply the aforementioned variable selection methods to those ultrahigh-dimensional statistical learning problems due to the computational complexity inherent in those methods. To reduce the predictor dimension in semiparametric regressions, Yu et al. (2013) proposed a ${\rho }_1$-minimisation of SIR with the Dantzig selector (Candes & Tao, 2007), which is defined as (5) $\begin{array}{rl}& min||{\eta }_k|{|}_{{\ell }_1},(l=1,\dots ,k-1)\\ & \text{such that }||\hat{\mathbf{M}}{\eta }_k-{\nu }_k\hat{\mathbf{\Sigma }}{\eta }_k||\\ & \le {\zeta }_k,|{\eta }_k^{\mathrm{\top }}\hat{\mathbf{\Sigma }}{\eta }_k-1|\le {\zeta }_k,|{\eta }_k^{\mathrm{\top }}\hat{\mathbf{\Sigma }}{\eta }_l|\le {\zeta }_k,\end{array}$(5) where $k=1,\dots ,d$, ${\nu }_k={\eta }_k^{\mathrm{\top }}\hat{\mathbf{M}}{\eta }_k$, $|\eta {|}_{{\ell }_1}=\sum _{i=1}^p|{\eta }_i|$ and ${\eta }_0$ is a $p\times 1$ zero vector. Furthermore, they established a non-asymptotic error bound for the resulting estimator when $|\mathcal{A}|$ is fixed. Yu et al. (2013) also extended the regularisation concept to SIR with an adaptive Dantzig selector, which is defined by (6) $\begin{array}{rl}& min||W_k{\eta }_k|{|}_{{\ell }_1},\\ & \text{such that }||W_k^{-1}(\hat{\mathbf{M}}{\widehat{\text{β}}}_k^0-{\hat{\lambda }}_k^0\hat{\mathbf{\Sigma }}{\eta }_k)||\le {\zeta }_k,\end{array}$(6) where $W_k=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(w_{k1},\dots ,w_{kp})$ is the a known weight matrix and $w_{kj}$ is a specified positive value, which should vary inversely with the magnitude of ${\widehat{\text{β}}}_{kj}^0$. Yu et al. (2013) proposed a two-step estimation procedure to select the contributing predictors. In the first step, they screened out informative predictors based on (5(5) $\begin{array}{rl}& min||{\eta }_k|{|}_{{\ell }_1},(l=1,\dots ,k-1)\\ & \text{such that }||\hat{\mathbf{M}}{\eta }_k-{\nu }_k\hat{\mathbf{\Sigma }}{\eta }_k||\\ & \le {\zeta }_k,|{\eta }_k^{\mathrm{\top }}\hat{\mathbf{\Sigma }}{\eta }_k-1|\le {\zeta }_k,|{\eta }_k^{\mathrm{\top }}\hat{\mathbf{\Sigma }}{\eta }_l|\le {\zeta }_k,\end{array}$(5) ). This is called Dantzig selector based SIR. In the second step, they enhance the sparsity and the estimation efficiency with (6(6) $\begin{array}{rl}& min||W_k{\eta }_k|{|}_{{\ell }_1},\\ & \text{such that }||W_k^{-1}(\hat{\mathbf{M}}{\widehat{\text{β}}}_k^0-{\hat{\lambda }}_k^0\hat{\mathbf{\Sigma }}{\eta }_k)||\le {\zeta }_k,\end{array}$(6) ), based on the predictors selected in the first step, called iterative adaptive Dantzig selector based SIR. This ensures that all contributing predictors are selected with high probability and that the resulting estimator is asymptotically normal even when the predictor dimension diverges to infinity.

However, there is a gap between the optimisation problem and the theoretical results: there is no guarantee that the estimator obtained from solving the proposed biconvex optimisation problem is the global minimum. Most existing work in the high-dimensional sufficient dimension reduction literature involves nonconvex optimisation problems. Moreover, they seek to estimate a set of reduced predictors that are not identifiable by definition, rather than the central subspace. Yin and Hilafu (2015) proposed a sequential approach for estimating high-dimensional SIR. Both proposals are stepwise procedures that do not correspond to solving a convex optimisation problem. Moreover, as discussed in Yin and Hilafu (2015), theoretical properties for their proposed estimators are hard to establish due to the sequential procedure used to obtain the estimators. In the high-dimensional setting, Lin et al. (2018) proposed a screening approach to perform variable selection and established an error bound for the estimators, which allows $|\mathcal{A}|$ goes to infinity. The selected variables are then used to fit classic SIR. Furthermore, the resulting algorithm is shown to be consistent and achieved the optimal convergence rate under certain sparsity conditions when p is of order $o\left(n^2{c}^2\right)$, where c is the generalised signal-to-noise ratio. Tan et al. (2018) proposed a convex formulation for sparse SIR in the high-dimensional setting by adapting techniques from the sparse canonical correlation analysis. Their proposal estimates the central subspace directly and performs variable selection simultaneously. Moreover, the proposed method can be adapted for sufficient dimension reduction methods that can be formulated as generalised eigenvalue problems.

As mentioned in introduction, most literature mainly focus on SIR with consistency on variable selection. Qian et al. (2019) proposed methods under a unified minimum discrepancy framework with regularisation. Consistency results in both central subspace estimation and variable selection are established simultaneously for some famous SDR methods, including SIR, PFC and SAVE. More importantly, their approach allows many quantities such as the structural dimension, the number of important predictors and the number of slices to diverge with n. Unlike many high-dimensional SDR methods, their method did not necessarily require a sparsity condition on the predictor covariance matrix or the maximum eigenvalue of the predictor covariance matrix to be upper bounded. Furthermore, they developed a new algorithm that can efficiently solve a general class of high-dimensional sparse minimum discrepancy problems.

Many SDR methods can be rewritten as a minimisation problem using an objective function of the form (7) $L_{1n}(\mathrm{\Gamma },\mathbf{V})=\mathrm{t}\mathrm{r}\left(({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V}{)}^{\mathrm{\top }}{\mathrm{\Omega }}_n({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V})\right),$(7) where ${\gamma }_n$, ${\mathbf{\Sigma }}_n$ and ${\mathrm{\Omega }}_n$ are sample estimates for the population matrices $\tilde{\mathbf{M}}$, $\mathbf{W}$ and ${\mathbf{W}}^{-1}$. Here, $\tilde{\mathbf{M}}$ is a $p\times l$ kernel matrix associated with a particular SDR method, where $d\le l\le p$, $\mathbf{W}$ is some $p\times p$ positive definite matrix, $\mathrm{\Gamma }\in {\mathbb{R}}^{p\times d}$ and $\mathbf{V}\in {\mathbb{R}}^{d\times l}$ represent parameters to be estimated by minimisation of $L_{1n}$. The general form of (7(7) $L_{1n}(\mathrm{\Gamma },\mathbf{V})=\mathrm{t}\mathrm{r}\left(({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V}{)}^{\mathrm{\top }}{\mathrm{\Omega }}_n({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V})\right),$(7) ) is an adaptation of the minimum discrepancy approach proposed by Cook and Ni (2005). To identify the correct sparsity structure of ${\mathcal{S}}_{Y|\mathbf{X}}$ under $p\gg n$ scenarios, Qian et al. (2019) proposed to adopt coordinate-independent regularisation approach and imposed the penalty $P_{\mathbf{V}}\left(\mathrm{\Gamma }\right)$ with tuning parameter ${\lambda }_n$ on (7(7) $L_{1n}(\mathrm{\Gamma },\mathbf{V})=\mathrm{t}\mathrm{r}\left(({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V}{)}^{\mathrm{\top }}{\mathrm{\Omega }}_n({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V})\right),$(7) ) under the alternative constraint $\mathbf{V}{\mathbf{V}}^{\mathrm{\top }}={\mathbf{I}}_d$, given the objective function $\begin{array}{rl}{L}_{2n}(\mathrm{\Gamma },\mathbf{V})& =\frac{1}{2}\mathrm{t}\mathrm{r}\left(({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V}{)}^{\mathrm{\top }}{\mathrm{\Omega }}_n({\gamma }_n-{\mathbf{\Sigma }}_n\mathrm{\Gamma }\mathbf{V})\right)\\ & +{\lambda }_n{P}_{\mathbf{V}}\left(\mathrm{\Gamma }\right),\text{subject to }\mathbf{V}{\mathbf{V}}^{\mathrm{\top }}={\mathbf{I}}_d.\end{array}$ Given its minimiser $(\hat{\mathrm{\Gamma }},\hat{\mathbf{V}})=\arg \underset{\mathrm{\Gamma },\mathbf{V}}{min}{L}_{2n}(\mathrm{\Gamma },\mathbf{V})$, they simultaneously estimated ${\mathcal{S}}_{Y|\mathbf{X}}$ by $\mathbf{S}\mathbf{p}\mathbf{a}\mathbf{n}\left(\mathrm{\Gamma }\right)$ and estimated ${\mathcal{A}}_0$ by ${\hat{\mathcal{A}}}_0=\{1\le j\le p:e_j^{\mathrm{\top }}\hat{\mathrm{\Gamma }}{\hat{\mathrm{\Gamma }}}^{\mathrm{\top }}{e}_j>0\}$.

Tan et al. (2020) considered four loss functions

1. General loss. $L_G(\widehat{\text{β}},\text{β})=||\widehat{\text{β}}{\widehat{\text{β}}}^{\mathrm{\top }}-\text{β}{\text{β}}^{\mathrm{\top }}||$;

2. Projection loss. $L_P(\widehat{\text{β}},\text{β})=||\widehat{\text{β}}({\widehat{\text{β}}}^{\mathrm{\top }}\widehat{\text{β}}{)}^{-1}{\widehat{\text{β}}}^{\mathrm{\top }}-\text{β}({\text{β}}^{\mathrm{\top }}\text{β}{)}^{-1}{\text{β}}^{\mathrm{\top }}|{|}_F^2$;

3. Prediction loss. $L_X(\widehat{\text{β}},\text{β})=\underset{W\in {\mathbb{R}}^{p\times p}}{inf}||{\mathbf{\Sigma }}^{1/2}(\widehat{\text{β}}-\text{β}W)|{|}_F^2$;

4. Correlation loss. $L_C(\widehat{\text{β}},\text{β})=1-\frac{1}{d}\mathrm{T}\mathrm{r}\left[\right({\widehat{\text{β}}}^{\mathrm{\top }}\mathbf{\Sigma }\widehat{\text{β}}{)}^{-1}\left({\widehat{\text{β}}}^{\mathrm{\top }}\mathbf{\Sigma }\text{β}\right)\left({\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }\text{β}{)}^{-1}\right({\text{β}}^{\mathrm{\top }}\mathbf{\Sigma }\widehat{\text{β}}\left)\right]$,

where $||\cdot |{|}_F$ denotes the Frobenius norm of a matrix. Further, Tan et al. (2020) established the minimax lower bound for sparse SIR under general loss, projection loss and prediction loss. They proposed natural sparse SIR estimator and proved that the upper error bound associated with all four loss functions can match the minimax lower bound obtained, which implies that it is a rate-optimal estimator for sparse SIR. However, this optimal estimation is computational intractable. Then they developed the computational feasible counterpart for this natural sparse SIR estimator through convex relaxation. But their theoretical investigation suggested that such computational realisation for natural sparse SIR estimator cannot maintain the optimal estimation rate.

To further address this issue, they proposed a refined sparse SIR estimator. The refined sparse SIR estimator is also rate-optimal yet computational intractable. However, its computational feasible counterpart based on the adaptive estimation procedure is proven to be nearly rate-optimal. Compared to the Lasso-SIR (Lin et al., 2019), which was shown to be rate optimal only when $p=o\left(n^2\right)$, their sparse SIR approach is rate optimal even when $\log p=o\left(n\right)$. Therefore, their proposed sparse SIR estimator certainly enjoys a much wider range of applications. The reason why Lasso-SIR fails to work when $\log p=o\left(n\right)$ is that it requires the estimation of the eigenvalues and eigenvectors of the $p\times p$ non-sparse SIR kernel matrix. It is well known that the sample eigenvalues and eigenvectors are not even consistent when p/n has a nonzero limit as $n\to \mathrm{\infty }$. In summary, the minimax lower bound obtained, the two rate-optimal yet computational infeasible estimators, the two corresponding computational tractable counterparts, and the theoretical upper bound of the four estimators under four-loss functions together provide a thorough understanding of sparse SIR. It is also worth noting that Lin et al. (2019) is just the first step of Tan et al. (2020). Bondell and Li (2009) demonstrated that $\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\text{β}\right)=\mathcal{A}$, then the sparse representation of SIR relies on $|\mathcal{A}|$, the number of truly relevant predictors, where $\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\text{β}\right)$ denotes the support of $\text{β}$. Assuming $|\mathcal{A}|\le s$, sparse SIR is further defined through seeking $\text{β}$ such that (8) $\begin{array}{rl}\text{β}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}\mathbf{M}\mathbf{B}\right)\text{s.t. }{\mathbf{B}}^{\mathrm{\top }}\mathbf{\Sigma }\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s.\end{array}$(8) The above formulation of sparse SIR enjoys a similar fashion as that of sparse CCA (Gao et al., 2015). To get theoretical results, the following conditions are required.

(A1)	the conditional mean $E\left\{\right(\mathbf{X}-E\left(\mathbf{X}\right)|{\text{β}}^{\mathrm{\top }}\mathbf{X}\left)\right\}$ is linear in $\mathbf{X}$;
(A2)	$|{\mathcal{A}}_k|$ is bounded for $k=1,\dots ,d$;
(A3)	the nonzero eigenvalues ${\lambda }_1,\dots ,{\lambda }_d$ are distinct;
(A4)	there exists a positive constant $a_0$ such that $0<a_0<1/4$, $\log p/n\le a_0$ and $E(\exp [t\left\{{\mathbf{X}}_i-E\right({\mathbf{X}}_i\left){\}}^2\right])\le K<\mathrm{\infty }$, for $i=1,\dots ,p$, and all $|t|\le a_0$;
(A5)	$D_0=\underset{1\le i\le p}{max}\sum _{j=1}^p|{\sigma }_{ij}|$ is bounded constant as $p\to \mathrm{\infty }$;
(A6)	the restricted isometry and restricted orthogonality constants ${\delta }_{2S_k}^{A^k}$ and ${\theta }_{S_k,2S_k}^{A^K}$ satisfy ${\delta }_{2S_k}^{A^k}+{\theta }_{S_k,2S_k}^{A^K}<1$, where $A^k=\mathbf{M}-{\lambda }_k\mathbf{\Sigma }$.

See Yu et al. (2013) for more details.

Theorem 3.2

Suppose that Conditions A1–A6 are satisfied, and ${\zeta }_k=C_0(\log p/n{)}^{1/2}$. Then $||{\widehat{\text{β}}}_k-{\text{β}}_k|{|}^2\le \frac{4C^2{S}_K}{(1-{\delta }_{2S_k}^{A^k}-{\theta }_{S,2S}^{A^K}{)}^2}\frac{\log p}{n},$ with a probability greater than $1-58p^{-\tau }$ for some τ greater than $(\log p{)}^{-1}\log 58,$ where $A^k,$ ${\delta }_{2S_k}^{A^k}$ and ${\theta }_{S,2S}^{A^K}$ are defined in Condition A6.

Remark 3.2

Theorem 3.2 suggests that a small price can obtain a sparse solution, as the squared estimation error of the regularised estimation is optimal up to a factor of $\log p$. The consistency property of Lin et al. (2018), Tan et al. (2018) and Tan et al. (2020) are similar with Theorem 3.2, but the threshold for $||{\widehat{\text{β}}}_k-{\text{β}}_k|{|}^2$ can be different.

4. The current literature of variable screening

Although there is a vast literature of applying sufficient dimension reduction for model-free selection, the result of developing screening consistency for the ultra-high dimensional setting is scant. Therefore many scholars are concentrated on investigating methods to achieve screening consistency.

4.1. Marginal utility

Yu, Dong, Shao (2016) proposed an approach called marginal SIR for model-free variable selection. Since $\mathbf{M}$ contains all the regression information between $Y$ and $\mathbf{X}$, Yu, Dong, Shao (2016) considered the diagonal element of $\mathbf{M}$ as the marginal utility for the corresponding predictor. Specially, let $e_k$ be the standard unit vector in ${\mathbb{R}}^p$ with 1 being the kth element and 0 otherwise. They considered the following utility for ${\mathbf{X}}_k$: (9) $m_k=e_k^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathbf{M}{\mathbf{\Sigma }}^{-1}{e}_k.$(9) Yu, Dong, Shao (2016) refer to $m_k$ as the population level marginal SIR utility. To apply Dantzig selector for the estimation of the marginal SIR utility $m_k$, they defined $p_{\ell }=E\left\{\text{𝟙}\right(Y\in {\mathbf{J}}_{\ell }\left)\right\}$, $\ell =1,\dots ,H$. Let ${\mu }_{\ell }=E\left\{\mathbf{X}\text{𝟙}\right(Y\in {\mathbf{J}}_{\ell }\left)\right\}$. Then ${\mathbf{M}}_{\mathrm{S}\mathrm{I}\mathrm{R}}=\sum _{\ell =1}^H{p}_{\ell }E(\mathbf{X}|Y\in {\mathbf{J}}_{\ell })E({\mathbf{X}}^{\mathrm{\top }}|Y\in {\mathbf{J}}_{\ell })$ can be written as ${\mathbf{M}}_{\mathrm{S}\mathrm{I}\mathrm{R}}=\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }$. Therefore $m_k^{\mathrm{S}\mathrm{I}\mathrm{R}}=e_k^{\mathrm{\top }}\left(\sum _{\ell =1}^H{\nu }_{\ell }{\nu }_{\ell }^{\mathrm{\top }}/p_{\ell }\right)e_k,{\nu }_{\ell }={\mathbf{\Sigma }}^{-1}{\mu }_{\ell }.$ The marginal utility $m_k$ is estimated by ${\hat{m}}_k^{\mathrm{S}\mathrm{I}\mathrm{R}}=e_k^{\mathrm{\top }}\left(\sum _{\ell =1}^H{\hat{\nu }}_{\ell }{\hat{\nu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }\right)e_k,$ where $\hat{\mathbf{\Sigma }}=\sum _{i=1}^n{\mathbf{X}}_i{\mathbf{X}}_i^{\mathrm{\top }}/n$ and ${\hat{\mu }}_{\ell }=\sum _{i=1}^n{\mathbf{X}}_i\text{𝟙}(Y_i\in {\mathbf{J}}_{\ell })/n$. For a given threshold $b_n$, the active set $\mathcal{A}$ is estimated by including the predictors such that ${\hat{m}}_k^{\mathrm{S}\mathrm{I}\mathrm{R}}$ exceeds $b_n$ or $\hat{\mathcal{A}}=\{k:{\hat{m}}_k^{\mathrm{S}\mathrm{I}\mathrm{R}}\ge b_n\}$. Yu, Dong, Shao (2016) take an example that $\mathbf{X}=(X_1,\dots ,X_p{)}^{\mathrm{\top }}\sim \mathcal{N}(0,\mathbf{\Sigma })$. Let $\mathrm{v}\mathrm{a}\mathrm{r}\left(X_i\right)=1$, $\mathrm{c}\mathrm{o}\mathrm{v}(X_i,X_j)=0.6$ for $|i-j|=1$, and $\mathrm{c}\mathrm{o}\mathrm{v}(X_i,X_j)=0$ for $|i-j|>1$, $1\le i,j\le p$. Let $Y={\text{β}}^{\mathrm{\top }}\mathbf{X}+\epsilon$, where $\epsilon \sim \mathcal{N}(0,1)$ is independent of $\mathbf{X}$ and $\text{β}=(1.2,-2,0.\dots ,0{)}^{\mathrm{\top }}$. Then the active set for the linear regression models is $\mathcal{A}=\{1,2\}$. Consider five utilities for $X_1$: the marginal absolute Pearson correlation from Fan and Lv (2008), the marginal squared distance correlation utility from Li et al. (2012), the marginal fused Kolmogorov filter utility as defined in (5.3) of Yu, Dong, Shao (2016), the marginal independence SIR utility as defined in (5.1) of Yu, Dong, Shao (2016) and the marginal SIR utility as defined in (9(9) $m_k=e_k^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathbf{M}{\mathbf{\Sigma }}^{-1}{e}_k.$(9) ). Unfortunately, the first four independence screening methods will fail to recover the active predictor $X_1$, only marginal SIR achieves desired result.

4.2. Trace pursuit

Yu, Dong, Zhu (2016) proposed trace pursuit as a novel approach for model-free variable selection. They first extended the classical stepwise regression in linear models and proposed an STP algorithm for model-free variable selection. Furthermore, they proposed the FTP algorithm. After finding a solution path by adding one predictor into the model at a time, a modified Bayesian information criterion (BIC, hereafter) provides a chosen model that is guaranteed to include all important predictors. Finally, the two-stage trace pursuit algorithm uses FTP for initial variable screening.

For working index set $\mathcal{F}$ and index $j\in {\mathcal{F}}^c$, if we want to test (10) $\begin{array}{rl}& H_0:Y\text{boxhU}{X}_j|{\mathbf{X}}_{\mathcal{F}}\text{ vs. }{H}_a:\\ & Y\text{ is not independent of }{X}_j\text{ given }{\mathbf{X}}_{\mathcal{F}}.\end{array}$(10) For any index set $\mathcal{F}$, denote ${\mathbf{X}}_{\mathcal{F}}=\{X_i:i\in \mathcal{F}\}$, $\mathrm{v}\mathrm{a}\mathrm{r}\left({\mathbf{X}}_{\mathcal{F}}\right)={\mathbf{\Sigma }}_{\mathcal{F}}$. Taking SIR as an example, denote ${\mathbf{M}}^{\mathrm{S}\mathrm{I}\mathrm{R}}={\mathbf{\Sigma }}_{\mathcal{F}}^{-1/2}{\mathbf{M}}_{\mathrm{S}\mathrm{I}\mathrm{R}}{\mathbf{\Sigma }}_{\mathcal{F}}^{-1/2}$. Recall that $\mathcal{A}$ denotes the active index set satisfying $Y\text{boxhU}{\mathbf{X}}_{{\mathcal{A}}^c}|{\mathbf{X}}_{\mathcal{A}}$, and $\mathcal{I}=\{1,\dots ,p\}$ denotes the full index set. It is worth noting that, if the assumption (W1) holds true, then for any index set $\mathcal{F}$ such that $\mathcal{A}\subseteq \mathcal{F}\subseteq \mathcal{I}$, $\mathrm{t}\mathrm{r}\left({\mathbf{M}}_{\mathcal{A}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right)=\mathrm{t}\mathrm{r}\left({\mathbf{M}}_{\mathcal{F}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right)=\mathrm{t}\mathrm{r}\left({\mathbf{M}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right)$. It suggests that $\mathrm{t}\mathrm{r}\left({\mathbf{M}}_{\mathcal{F}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right)$ can be used to capture the strength of relationship between $Y$ and ${\mathbf{X}}_{\mathcal{F}}$. Denote $\mathcal{F}\cup j$ as the index set of j together with all the indices in $\mathcal{F}$. Given that ${\mathbf{X}}_{\mathcal{F}}$ is already in the model, then trace difference $\mathrm{t}\mathrm{r}\left({\mathbf{M}}_{\mathcal{F}\cup j}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right)-\mathrm{t}\mathrm{r}\left({\mathbf{M}}_{\mathcal{F}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right)$ can be used to test the contribution of the additional variable $X_j$ to $Y$. The idea of using trace difference is similar to the extra sums of squares test in the classical multiple linear regression setting. The following subset LCM assumption is required in Yu, Dong, Zhu (2016), (11) $\begin{array}{rl}& E\left(X_j|{\mathbf{X}}_{\mathcal{F}}\right)\text{ is a linear function of }{\mathbf{X}}_{\mathcal{F}}\text{ for any }\mathcal{F}\\ & \subset \mathcal{I}\text{ and }j\in {\mathcal{F}}^c.\end{array}$(11) Furthermore, they also provided the STP algorithm, that is

1. Initialisation. Set the initial working set to be $\mathcal{F}=\mathrm{\varnothing }$.

2. Forward addition. Find index $a_{\mathcal{F}}$ such that $a_{\mathcal{F}}=\arg \underset{j\in {\mathcal{F}}^c}{max}\mathrm{t}\mathrm{r}\left({\hat{\mathbf{M}}}_{\mathcal{F}\cup j}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right).$ If $T_{a_{\mathcal{F}}|\mathcal{F}}^{\mathrm{S}\mathrm{I}\mathrm{R}}=n\left\{\mathrm{t}\mathrm{r}\right({\hat{\mathbf{M}}}_{\mathcal{F}\cup a_{\mathcal{F}}}^{\mathrm{S}\mathrm{I}\mathrm{R}})-\mathrm{t}\mathrm{r}({\hat{\mathbf{M}}}_{\mathcal{F}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\left)\right\}>c^{\mathrm{S}\mathrm{I}\mathrm{R}}$, update $\mathcal{F}$ to be $\mathcal{F}\cup a_{\mathcal{F}}$.

3. Backward deletion. Find index $d_{\mathcal{F}}$ such that $d_{\mathcal{F}}=\arg \underset{j\in \mathcal{F}}{max}\mathrm{t}\mathrm{r}\left({\hat{\mathbf{M}}}_{\mathcal{F}\mathrm{\setminus }j}^{\mathrm{S}\mathrm{I}\mathrm{R}}\right).$ If $T_{d_{\mathcal{F}}|\mathcal{F}\mathrm{\setminus }{d}_{\mathcal{F}}}^{\mathrm{S}\mathrm{I}\mathrm{R}}=n\left\{\mathrm{t}\mathrm{r}\right({\hat{\mathbf{M}}}_{\mathcal{F}}^{\mathrm{S}\mathrm{I}\mathrm{R}})-\mathrm{t}\mathrm{r}({\hat{\mathbf{M}}}_{\mathcal{F}\mathrm{\setminus }{d}_{\mathcal{F}}}^{\mathrm{S}\mathrm{I}\mathrm{R}}\left)\right\}<c_{\mathrm{S}\mathrm{I}\mathrm{R}}$, update $\mathcal{F}$ to be $\mathcal{F}\mathrm{\setminus }{d}_{\mathcal{F}}$.

4. Repeat steps (b) and (c) until no predictors can be added or deleted.

The test for SAVE and DR can be defined in a parallel fashion if the following CCV assumption together with the subset LCM (11(11) $\begin{array}{rl}& E\left(X_j|{\mathbf{X}}_{\mathcal{F}}\right)\text{ is a linear function of }{\mathbf{X}}_{\mathcal{F}}\text{ for any }\mathcal{F}\\ & \subset \mathcal{I}\text{ and }j\in {\mathcal{F}}^c.\end{array}$(11) ) assumption holds true $\mathrm{v}\mathrm{a}\mathrm{r}\left(X_j|{\mathbf{X}}_{\mathcal{F}}\right)\text{ is nonrandom for any }\mathcal{F}\subseteq \mathcal{I}\text{ and }j\in {\mathcal{F}}^c.$ Li and Dong (2020) had a recent extension of trace pursuit to matrix-valued predictors. Suppose the response variable $Y\in \mathbb{R}$ and the predictor $\mathbf{X}\in {\mathbb{R}}^{p\times q}$ have the following general relationship: (12) $Y=g\left(\mathbf{X}\right)+\epsilon ,$(12) where g: ${\mathbb{R}}^{p\times q}\to \mathbb{R}$ is an unknown function, ϵ is independent of $\mathbf{X}$, and $E\left(\epsilon \right)=0$. Assume that $\mathbf{X}$ follows the matrix normal distribution, which is denoted as $\mathbf{X}\sim {\mathcal{N}}_{p,q}(\mu ,\mathbf{U},\mathbf{V})$ with $\mu \in {\mathbb{R}}^{p\times q}$, $\mathbf{U}\in {\mathbb{R}}^{p\times p}$ and $\mathbf{V}\in {\mathbb{R}}^{q\times q}$. Then, the row covariance matrix is $\mathbf{U}=E\left\{\right(\mathbf{X}-\mathit{\mu }\left)\right(\mathbf{X}-\mathit{\mu }{)}^{\mathrm{\top }}\left\}/\mathrm{t}\mathrm{r}\right(\mathbf{V})$, and the column covariance matrix is $\mathbf{V}=E\left\{\right(\mathbf{X}-\mathit{\mu }{)}^{\mathrm{\top }}(\mathbf{X}-\mathit{\mu })\left\}/\mathrm{t}\mathrm{r}\right(\mathbf{U})$.

Let ${\mathcal{I}}_{row}=\{1,\dots ,p\}$ be the full index set of rows and ${\mathbf{X}}_{j,\cdot }$ be the jth row of $\mathbf{X}$ for $j=1,\dots ,p$. Define the active row set $\mathcal{A}$ as $\mathcal{A}=\{j\in {\mathcal{I}}_{row}:Y\text{ depends on }{\mathbf{X}}_{j,\cdot }\text{ in model }(12\left)\right\}.$ Similarly, let ${\mathcal{I}}_{col}=\{1,\dots ,q\}$ be the full index set of columns and ${\mathbf{X}}_{\cdot ,k}$ be the kth column of $\mathbf{X}$ for $k=1,\dots ,q$. Define the active column set $\mathcal{B}$ as $\mathcal{B}=\{k\in {\mathcal{I}}_{col}:Y\text{ depends on }{\mathbf{X}}_{\cdot ,k}\text{ in model }(12\left)\right\}.$ Based on the active row and column predictors, model (12(12) $Y=g\left(\mathbf{X}\right)+\epsilon ,$(12) ) can be expressed as $Y=g^{\ast }\left({\mathbf{X}}_{\mathcal{A},\mathcal{B}}\right)+\epsilon ,$ where $g^{\ast }:{\mathbb{R}}^{|\mathcal{A}|\times |\mathcal{B}|}\to \mathbb{R}$ with $|\cdot |$ denoting the cardinality of a set, and ${\mathbf{X}}_{\mathcal{A},\mathcal{B}}$ denotes the submatrix of $\mathbf{X}$ that contains the active rows indexed by $\mathcal{A}$ and the active columns indexed by $\mathcal{B}$. Note that Y depends on $\mathbf{X}$ only through ${\mathbf{X}}_{\mathcal{A},\mathcal{B}}$. Li and Dong (2020) introduced procedures to recover the active row set $\mathcal{A}$ in detail. Let ${\mathbf{X}}_{j,\cdot }$, $j=1,\dots ,p$, be the jth row of $\mathbf{X}$ and ${\mathbf{X}}_{-j,\cdot }\in {\mathbb{R}}^{(p-1)\times q}$ be the matrix that includes all but the jth row of $\mathbf{X}$. To test the importance of ${\mathbf{X}}_{j,\cdot }$, they considered the following row hypotheses: (13) $\begin{array}{rl}& H_{0,\left\{j\right\}}^{row}:Y\text{boxhU}\mathbf{X}|{\mathbf{X}}_{-j,\cdot }\text{ vs. }{H}_{a,\left\{j\right\}}^{row}:\\ & Y\text{ is not independent of }\mathbf{X}\text{ given }{\mathbf{X}}_{-j,\cdot }.\end{array}$(13) Under the null hypothesis, $H_{0,\left\{j\right\}}^{row}$, the response Y depends on $\mathbf{X}$ only through ${\mathbf{X}}_{-j,\cdot }$. In the special case of q = 1, $\mathbf{X}$ becomes a p-dimensional vector, and (13(13) $\begin{array}{rl}& H_{0,\left\{j\right\}}^{row}:Y\text{boxhU}\mathbf{X}|{\mathbf{X}}_{-j,\cdot }\text{ vs. }{H}_{a,\left\{j\right\}}^{row}:\\ & Y\text{ is not independent of }\mathbf{X}\text{ given }{\mathbf{X}}_{-j,\cdot }.\end{array}$(13) ) is equivalent to testing the importance of one component of $\mathbf{X}$ given the other p−1 predictors. This special case is known as the marginal coordinate test (Cook, 2004). Let ${\mathbf{U}}_{-j,-j}\in {\mathbb{R}}^{(p-1)\times (p-1)}$ be the submatrix of $\mathbf{U}$ that excludes the jth row and the jth column of $\mathbf{U}$. Define the following quantity: ${\delta }_j^{row}=\mathrm{t}\mathrm{r}\left(\mathbf{M}\right)-\mathrm{t}\mathrm{r}\left({\mathbf{M}}_{-j,\cdot }\right),$ where $\mathbf{M}={\mathbf{U}}^{-1}E\left(\mathbf{X}Y\right){\mathbf{V}}^{-1}{E}^{\mathrm{\top }}\left(\mathbf{X}Y\right)$ and ${\mathbf{M}}_{-j,\cdot }={\mathbf{U}}_{-j,-j}^{-1}E\left({\mathbf{X}}_{-j,\cdot }Y\right){\mathbf{V}}^{-1}{E}^{\mathrm{\top }}\left({\mathbf{X}}_{-j,\cdot }Y\right)$. This trace difference ${\delta }_j^{row}$ is the key quantity to test the importance of the jth row of $\mathbf{X}$, which is same as Yu, Dong, Zhu (2016). Note that ${\delta }_j^{row}=0$ under $H_{0,\left\{j\right\}}^{row}$.

To develop the screening consistency for ultrahigh dimensional setting, Zhu et al. (2011) proposed a novel variable screening procedure under a unified model framework, which covers a wide variety of commonly used parametric and semiparametric models. They assumed that $E\left({\mathbf{X}}_i\right)=0$ and $\mathrm{v}\mathrm{a}\mathrm{r}\left(x_i\right)=1$ for $i=1,\dots ,p$ and $\mathrm{\Omega }\left(Y\right)=E\left\{\mathbf{X}F\right(Y|\mathbf{X}\left)\right\}$ for ease of explanation. It then follows by the law of iterated expectations that $\mathrm{\Omega }\left(Y\right)=\mathrm{c}\mathrm{o}\mathrm{v}\{\mathbf{X},\text{𝟙}(Y<y\left)\right\}$. Let ${\mathrm{\Omega }}_i\left(Y\right)$ be the ith element of $\mathrm{\Omega }\left(Y\right)$, and defined as ${\omega }_i=E\left\{{\mathrm{\Omega }}_i^2\right(Y\left)\right\},\dots ,i=1,\dots ,p.$ Then ${\omega }_i$ is to serve as the population quantity of our proposed marginal utility measure for predictor ranking. Intuitively, one can see that, if $x_i$ and $Y$ are independent, then $x_i$ and the indicator function $\text{𝟙}(Y\le y)$ change independently. Consequently, ${\omega }_i=0$. On the other hand, if $x_i$ and $Y$ are related, then ${\omega }_i$ must be positive. For ease of presentation, they assumed that the sample predictors are all standardised; that is, $n^{-1}\sum _{j=1}^n{\mathbf{X}}_{ji}=0$ and $n^{-1}\sum _{j=1}^n{\mathbf{X}}_{ji}^2=1$ for $i=1,\dots ,p$. A natural estimator of ${\omega }_i$ is ${\tilde{\omega }}_i=\frac{1}{n}\sum _{j=1}^n{\left\{\frac{1}{n}\sum _{k=1}^n{\mathbf{X}}_{ki}\text{𝟙}(Y_k<Y_j)\right\}}^2,i=1,\dots ,p,$ where ${\mathbf{X}}_{ki}$ denotes the kth element of $x_i$. The new method does not require imposing a specific model structure on regression functions, and thus is particularly appealing to ultrahigh-dimensional regressions. They showed that, with the number of predictors growing at an exponential rate of the sample size, the proposed procedure possesses consistency in ranking, which is both useful in its own right and can lead to consistency in selection. Lin et al. (2017) first introduced a large class of models depending on the smallest non-zero eigenvalue λ of the kernel matrix of SIR, then the determination of the minimax rate for estimating the central space over two classes is derived, which is the first paper that studied the minimax estimation of sparse SIR. Furthermore, they showed that the estimator based on the SIR procedure converges at rate $d\wedge \left(\right(sd+s\log \left(ep/s\right)\left)/\right(n\lambda \left)\right)$, which is the optimal rate for the single index models and multiple index models with fixed structural dimension d, fixed $s=|\mathcal{A}|$ and λ. However, Lin et al. (2017) only considered the projection loss (Li & Wang, 2007). More importantly, their theoretical study is actually based on the assumption that covariance matrix is diagonal.

Before discussing the consistency property, we need some conditions. Taking Yu, Dong, Shao (2016) as an example,

(C1)	The coverage condition: $\mathbf{S}\mathbf{p}\mathbf{a}\mathbf{n}\left\{\mathbf{\Sigma }E\right(\mathbf{X}|Y\in {\mathbf{J}}_{\ell })\ell =1,\dots ,H\}={\mathcal{S}}_{Y|\mathbf{X}}$.
(C2)	There exist $0<c<1/4$ and $0<q<\mathrm{\infty }$ such that $E\{\exp (t{\mathbf{X}}_k\left)\right\}\le q$ for all $|t|\le c$, $k=1,\dots ,p$. In addition, there exist positive constants ${\lambda }_{min}$ and ${\lambda }_{max}$ such that $0\le {\lambda }_{min}\le {\lambda }_{min}\left(\mathbf{\Sigma }\right)\le {\lambda }_{max}\left(\mathbf{\Sigma }\right)\le {\lambda }_{max}<\mathrm{\infty }$, where ${\lambda }_{min}\left(\mathbf{\Sigma }\right)$ and where ${\lambda }_{max}\left(\mathbf{\Sigma }\right)$ are the smallest and largest eigenvalue of $\mathbf{\Sigma }$, respectively.
(C3)	There exists $0<f<\mathrm{\infty }$ such that $||{\mathbf{\Sigma }}^{-1}|{|}_1\le f$.
(C4)	There exists $0<g<1-2\times r$ such that $f^2{s}^2\log p=O_p\left(n^g\right)$, where s is the cardinality of $\mathcal{A}$ and r is specified in condition (C5).
(C5)	There exists $0<a_2<\mathrm{\infty }$ and $r\le 1/2$ such that $\underset{k\in \mathcal{A}}{min}{m}_k>2a_2{n}^{-r}$.

More details please refer to Yu, Dong, Shao (2016).

Theorem 4.1

Assume above conditions hold, then the shrinkage estimator $\widehat{\text{β}}$ satisfies consistency in variable selection, $Pr(\hat{\mathcal{A}}\supseteq \mathcal{A})\to 1.$

Theorem 4.1 is given in Yu, Dong, Shao (2016), Yu, Dong, Zhu (2016), Lin et al. (2017), and Zhu et al. (2011).

5. Minimax rate

Recently, an impressive range of penalised SIR methods has been proposed to estimate the central subspace in a sparse fashion. However, few of them considered the sparse sufficient dimension reduction from a decision-theoretical point of view. To address this issue, Tan et al. (2020) established the minimax rates of convergence for estimating the sparse SIR directions under various commonly used loss functions in the literature of sufficient dimension reduction. Lin et al. (2019) introduced a simple Lasso regression method to obtain an estimator of the sufficient dimension reduction space, which is only the first step of Tan et al. (2020). Moreover, Tan et al. (2020) discovered the possible trade-off between statistical guarantee and computational performance for sparse SIR and proposed an adaptive estimation scheme for sparse SIR which is computationally tractable and rate optimal under the condition that $\log p=o\left(n\right)$, which is weaker than Lin et al. (2019).

As we can see that, the kernel matrix ${\mathbf{M}}_{\mathrm{S}\mathrm{I}\mathrm{R}}$ can be estimated as ${\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{I}\mathrm{R}}=\sum _{\ell =1}^H{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }.$ Then it is natural to estimate $\text{β}$ via replacing $\mathbf{M}$ and $\mathbf{\Sigma }$ in (8(8) $\begin{array}{rl}\text{β}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}\mathbf{M}\mathbf{B}\right)\text{s.t. }{\mathbf{B}}^{\mathrm{\top }}\mathbf{\Sigma }\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s.\end{array}$(8) ) by their sample estimators, which yields (14) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{I}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{I}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\end{array}$(14) The solution ${\widehat{\text{β}}}_{\mathrm{S}\mathrm{I}\mathrm{R}}$ in (14(14) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{I}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{I}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\end{array}$(14) ) is called the natural sparse SIR estimator. The following theorem establishes the lower bound and upper bound of the four loss functions for the natural sparse SIR estimator.

Theorem 5.1

Assume $n{\lambda }^2\ge C_0\log \frac{ep}{s}$ for some sufficiently large constant $C_0$. Then there exist positive constants C and $c_0$ such that $\begin{array}{rl}& \underset{\widehat{\text{β}}}{inf}\underset{\mathbb{P}\in \mathcal{P}}{sup}{E}_{\mathbb{P}}{L}_G(\widehat{\text{β}},\text{β})\ge C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\wedge c_0,\\ & \underset{\widehat{\text{β}}}{inf}\underset{\mathbb{P}\in \mathcal{P}}{sup}\mathbb{P}\left\{L_P(\widehat{\text{β}},\text{β})\ge C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\wedge c_0\right\}\ge 0.8,\\ & \underset{\widehat{\text{β}}}{inf}\underset{\mathbb{P}\in \mathcal{P}}{sup}\mathbb{P}\left\{L_X(\widehat{\text{β}},\text{β})\ge C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\wedge c_0\right\}\ge 0.8,\end{array}$ where $\mathcal{P}=\mathcal{P}(n,H,s,p,d,\lambda ;K,m)$.

Theorem 5.2

Assume that $\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\le c$ for some small constant $c\in (0,1)$. Then for any $C^{\mathrm{\prime }}>0,$ there exists a positive constant C such that $\begin{array}{rl}& L_G(\widehat{\text{β}},\text{β})\vee L_P(\widehat{\text{β}},\text{β})\vee L_X(\widehat{\text{β}},\text{β})\vee L_C(\widehat{\text{β}},\text{β})\\ & \le C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\end{array}$ with probability greater than $1-2\exp (-C^{\mathrm{\prime }}(s+\log \left(ep/s\right)\left)\right)$ uniformly over $\mathcal{P}\in P(n,H,s,p,d,\lambda ;K,m)$.

Since SIR can be rewritten as a least-square formulation, they finally proposed an adaptive estimation scheme for sparse SIR which is computationally tractable and rate optimal. More details about the adaptive sparse SIR estimator can refer to Tan et al. (2020).

6. Further investigation

6.1. Marginal utility

Motivated by Yu, Dong, Shao (2016), we can extend their method to SAVE and DR. Let ${\varphi }_{\ell }=E\left\{\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}\text{𝟙}\right(Y\in {\mathbf{J}}_{\ell }\left)\right\}$. Then ${\mathbf{M}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}=\sum _{\ell =1}^H{p}_{\ell }\{\mathbf{\Sigma }-\mathrm{v}\mathrm{a}\mathrm{r}(\mathbf{X}|Y\in {\mathbf{J}}_{\ell }){\}}^2$ can be written as ${\mathbf{M}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}=\sum _{\ell =1}^H{p}_{\ell }\{\mathbf{\Sigma }-{\varphi }_{\ell }/p_{\ell }+{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2{\}}^2$. Therefore $\begin{array}{rl}{m}_k^{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =e_k^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\left(\sum _{\ell =1}^H{p}_{\ell }\{\mathbf{\Sigma }-{\varphi }_{\ell }/p_{\ell }+{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2{\}}^2\right)\\ & \times {\mathbf{\Sigma }}^{-1}{e}_k.\end{array}$ The marginal utility $m_k$ is estimated by $\begin{array}{rl}{\hat{m}}_k^{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =e_k^{\mathrm{\top }}{\hat{\mathbf{\Sigma }}}^{-1}\left(\sum _{\ell =1}^H{\hat{p}}_{\ell }\{\hat{\mathbf{\Sigma }}-{\hat{\varphi }}_{\ell }/{\hat{p}}_{\ell }+{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }^2{\}}^2\right)\\ & \times {\hat{\mathbf{\Sigma }}}^{-1}{e}_k,\end{array}$ where ${\hat{\varphi }}_{\ell }=\sum _{i=1}^n{\mathbf{X}}_i{\mathbf{X}}_i^{\mathrm{\top }}\text{𝟙}(Y_i\in {\mathbf{J}}_{\ell })/n$. For a given threshold $b_n$, the active set $\mathcal{A}$ is estimated by including the predictors such that ${\hat{m}}_k^{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}$ exceeds $b_n$ or $\hat{\mathcal{A}}=\{k:{\hat{m}}_k^{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}\ge b_n\}$.

Next, we consider marginal DR with the Dantzig selector. Then $\begin{array}{rl}{\mathbf{M}}_{\mathrm{D}\mathrm{R}}& =\sum _{\ell =1}^{H}2p_{\ell }\left[E\right\{E^2\left(\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}\right)\}\\ & +E^2\left\{E\right(\mathbf{X}|Y\in {\mathbf{J}}_{\ell }\left)E\right({\mathbf{X}}^{\mathrm{\top }}|Y\in {\mathbf{J}}_{\ell }\left)\right\}\\ & +E\left\{E\right({\mathbf{X}}^{\mathrm{\top }}|Y\in {\mathbf{J}}_{\ell }\left)E\right(\mathbf{X}|Y\in {\mathbf{J}}_{\ell }\left)\right\}\\ & \times E\left\{E\right(\mathbf{X}|Y\in {\mathbf{J}}_{\ell }\left)E\right({\mathbf{X}}^{\mathrm{\top }}|Y\in {\mathbf{J}}_{\ell }\left)\right\}]-2\mathbf{\Sigma }\end{array}$ can be written as $\begin{array}{rl}{\mathbf{M}}_{\mathrm{D}\mathrm{R}}& =2\sum _{\ell =1}^H{p}_{\ell }{\left({\varphi }_{\ell }-\mathbf{\Sigma }\right)}^2+2{\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)}^2\\ & +2\left(\sum _{\ell =1}^H{\mu }_{\ell }^{\mathrm{\top }}{\mu }_{\ell }/p_{\ell }^2\right)\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right).\end{array}$ Therefore $\begin{array}{rl}{m}_k^{\mathrm{D}\mathrm{R}}& =2e_k^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\left\{\phantom{{\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)}^2}\sum _{\ell =1}^H{p}_{\ell }{\left({\varphi }_{\ell }-\mathbf{\Sigma }\right)}^2\right.\\ & +{\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)}^2\\ & +\left.\left(\sum _{\ell =1}^H{\mu }_{\ell }^{\mathrm{\top }}{\mu }_{\ell }/p_{\ell }^2\right)\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)\phantom{{\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)}^2}\right\}{\mathbf{\Sigma }}^{-1}{e}_k.\end{array}$ The marginal utility $m_k$ is estimated by $\begin{array}{rl}{\hat{m}}_k^{\mathrm{D}\mathrm{R}}& =2e_k^{\mathrm{\top }}{\hat{\mathbf{\Sigma }}}^{-1}\left\{\phantom{{\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)}^2}\sum _{\ell =1}^H{\hat{p}}_{\ell }{\left({\hat{\varphi }}_{\ell }-\hat{\mathbf{\Sigma }}\right)}^2\right.\\ & +{\left(\sum _{\ell =1}^H{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }^2\right)}^2\\ & +\left.\left(\sum _{\ell =1}^H{\hat{\mu }}_{\ell }^{\mathrm{\top }}{\hat{\mu }}_{\ell }/{\hat{p}}_{\ell }^2\right)\left(\sum _{\ell =1}^H{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }^2\right)\phantom{{\left(\sum _{\ell =1}^H{\mu }_{\ell }{\mu }_{\ell }^{\mathrm{\top }}/p_{\ell }^2\right)}^2}\right\}{\hat{\mathbf{\Sigma }}}^{-1}{e}_k.\end{array}$ For a given threshold $b_n$, the active set $\mathcal{A}$ is estimated by including the predictors such that ${\hat{m}}_k^{\mathrm{D}\mathrm{R}}$ exceeds $b_n$, or $\hat{\mathcal{A}}=\{k:{\hat{m}}_k^{\mathrm{D}\mathrm{R}}\ge b_n\}$. Following the proof of Yu, Dong, Shao (2016), we can expect the marginal SAVE and DR with the Dantzig selector to achieve selection consistency.

6.2. Minimax rate

Motivated by Tan et al. (2020), we can further investigate the natural sparse SAVE estimator and upper error bound. Let $E_n\mathbf{X}=\frac{1}{n}\sum _{i=1}^n{\mathbf{X}}_i$ and $\hat{\mathbf{\Sigma }}=\frac{1}{n-1}\sum _{i=1}^n({\mathbf{X}}_i-E_n\mathbf{X})({\mathbf{X}}_i-E_n\mathbf{X}{)}^{\mathrm{\top }}$ be the sample mean and sample covariance of $\mathbf{X}$, then the SAVE kernel matrix $\mathbf{M}$ is estimated as ${\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}=\sum _{\ell =1}^H{\hat{p}}_{\ell }\{\hat{\mathbf{\Sigma }}-{\hat{\varphi }}_{\ell }/{\hat{p}}_{\ell }+{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }^2{\}}^2.$ Similarly, the DR kernel matrix is estimated as $\begin{array}{rl}{\hat{\mathbf{M}}}_{\mathrm{D}\mathrm{R}}& =2\sum _{\ell =1}^H{\hat{p}}_{\ell }{\left({\hat{\varphi }}_{\ell }-\hat{\mathbf{\Sigma }}\right)}^2+2{\left(\sum _{\ell =1}^H{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }^2\right)}^2\\ & +2\left(\sum _{\ell =1}^H{\hat{\mu }}_{\ell }^{\mathrm{\top }}{\hat{\mu }}_{\ell }/{\hat{p}}_{\ell }^2\right)\left(\sum _{\ell =1}^H{\hat{\mu }}_{\ell }{\hat{\mu }}_{\ell }^{\mathrm{\top }}/{\hat{p}}_{\ell }^2\right).\end{array}$ Then it is natural to estimate $\text{β}$ via replacing $\mathbf{M}$ and $\mathbf{\Sigma }$ in (8(8) $\begin{array}{rl}\text{β}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}\mathbf{M}\mathbf{B}\right)\text{s.t. }{\mathbf{B}}^{\mathrm{\top }}\mathbf{\Sigma }\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s.\end{array}$(8) ) by their sample estimators, which yields (15) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\\ {\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{D}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s\end{array}$(15) The solution ${\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}$ and ${\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}$ in (15(15) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\\ {\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{D}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s\end{array}$(15) ) are called the natural sparse SAVE and DR estimator. The following theorem establishes the lower bound and upper bound of the four loss functions for the natural sparse SAVE and DR estimator.

Theorem 6.1

Assume $n{\lambda }^2\ge C_0\log \frac{ep}{s}$ for some sufficiently large constant $C_0$. Then there exist positive constants C and $c_0$ such that $\begin{array}{rl}& \underset{\widehat{\text{β}}}{inf}\underset{\mathbb{P}\in \mathcal{P}}{sup}{E}_{\mathbb{P}}{L}_G(\widehat{\text{β}},\text{β})\ge C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\wedge c_0,\\ & \underset{\widehat{\text{β}}}{inf}\underset{\mathbb{P}\in \mathcal{P}}{sup}\mathbb{P}\left\{L_P(\widehat{\text{β}},\text{β})\ge C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\wedge c_0\right\}\ge 0.8,\\ & \underset{\widehat{\text{β}}}{inf}\underset{\mathbb{P}\in \mathcal{P}}{sup}\mathbb{P}\left\{L_X(\widehat{\text{β}},\text{β})\ge C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\wedge c_0\right\}\ge 0.8,\end{array}$ where $\mathcal{P}=\mathcal{P}(n,H,s,p,d,\lambda ;K,m)$.

Theorem 6.2

Assume that $\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\le c$ for some small constant $c\in (0,1)$. Then for any $C^{\mathrm{\prime }}>0,$ there exists a positive constant C such that $\begin{array}{rl}& L_G(\widehat{\text{β}},\text{β})\vee L_P(\widehat{\text{β}},\text{β})\vee L_X(\widehat{\text{β}},\text{β})\vee L_C(\widehat{\text{β}},\text{β})\\ & \le C\frac{s\log \left(ep/s\right)}{n{\lambda }^2}\end{array}$ with probability greater than $1-2\exp (-C^{\mathrm{\prime }}(s+\log \left(ep/s\right)\left)\right)$ uniformly over $\mathcal{P}\in P(n,H,s,p,d,\lambda ;K,m)$. In which $\widehat{\text{β}}$ is constructed in (15(15) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\\ {\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{D}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s\end{array}$(15) ).

Following by Tan et al. (2020), $\widehat{\text{β}}$ in (15(15) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\\ {\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{D}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s\end{array}$(15) ) is rate optimal under general loss, projection loss and prediction loss. Moreover, the natural sparse SAVE estimator ${\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}$ and DR estimator ${\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}$ can also be regarded as one optimal estimator for the SAVE directions and DR directions. However, the estimation procedure (15(15) $\begin{array}{rl}{\widehat{\text{β}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{E}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s,\\ {\widehat{\text{β}}}_{\mathrm{D}\mathrm{R}}& =\underset{\mathbf{B}\in {\mathbb{R}}^{p\times d}}{\arg max}\mathrm{T}\mathrm{r}\left({\mathbf{B}}^{\mathrm{\top }}{\hat{\mathbf{M}}}_{\mathrm{D}\mathrm{R}}\mathbf{B}\right)\text{ s.t. }{\mathbf{B}}^{\mathrm{\top }}\hat{\mathbf{\Sigma }}\mathbf{B}\\ & ={\mathbf{I}}_d\text{ and }|\mathrm{S}\mathrm{u}\mathrm{p}\mathrm{p}\left(\mathbf{B}\right)|\le s\end{array}$(15) ) depends on the unknown sparsity parameter s and is computationally infeasible as it involves exhaustive search over all $\mathbf{B}\in {\mathbb{R}}^{p\times d}$ subject to the sparsity constraint. Tan et al. (2020) defined a refined sparse SIR estimator based on that SIR can be viewed as transformation-based projection pursuit. Since SAVE and DR cannot be rewritten as a least-square formulation, we do not define refined sparse SAVE and DR estimator.

Disclosure statement

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

Additional information

Funding

This research is supported by the National Natural Science Foundation of China Grant 11971170, the 111 project B14019 and the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning.

Notes on contributors

Lu Li

Lu Li is currently a Ph.D student at School of Statistics, East China Normal University.

Xuerong Meggie Wen

Dr Xuerong Meggie Wen is currently an associate professor of Statistics at Dept. of Mathematics and Statistics, Missouri University of Science and Technology.

Zhou Yu

Dr Zhou Yu is a Professor of Statistics at School of Statistics, East China Normal Univercity.