Discussion on ‘Review of sparse sufficient dimension reduction’

I congratulate the authors on an excellent overview of an important research area. Sufficient dimension reduction methods are based on the model-free driving condition that $Y\text{⫫}\mathbf{\text{X}}\mid {\mathbf{\text{P}}}_{\mathcal{S}}\mathbf{\text{X}}$, where $\mathbf{\text{X}}\in {\mathbb{R}}^p$ is multivariate and potentially high-dimensional, ${\mathbf{\text{P}}}_{\mathcal{S}}$ is the projection onto the dimension reduction subspace $\mathcal{S}\subseteq {\mathbb{R}}^p$. Equipped with variable selection and variable screening techniques, many modern sparse sufficient dimension reduction methods have been developed in the past few years, and they can work really well in the model development stage of high-dimensional data analysis. This review paper is very timely and provides a thorough overview of sparse sufficient dimension reduction methods and sheds lights on future research directions.

Specific contributions of this paper include the following. First of all, it recasts many moment-based sufficient dimension reduction methods as a generalised eigen-decomposition problem and also as a constrained trace optimisation (Equation (2.1) in the paper). These two formulations are crucial for studying sparse sufficient dimension reduction in high-dimensional settings where p>n or $p\gg n$. This paper also provides a comprehensive overview of the key developments in sparse sufficient dimension reduction literature, from the first paper by Ni et al. (2005), to the newest theoretical and computational breakthrough in Tan et al. (2020). Various techniques for sparse sufficient dimension reduction are discussed, including the shrinkage and regularised estimation, variable screening and the trace pursuit, and the sequential procedures. Advantages and possible pitfalls of different approaches are well explained by the authors. Finally, the paper provides the theoretical foundations for establishing the minimax rates of convergence for estimating a sparse dimension reduction subspace. Further calculations are also included to provide insights about second-order dimension reduction methods such as SAVE and DR.

This paper is inspiring for substantive research in high-dimensional multivariate statistics, and encouraging for current and future studies in dimension reduction. In what follows, I draw two connections which may be suggestive of future directions.

1. Least squares formulations

In the seminal work of Li and Duan (1989), a well-known connection between the ordinary least square estimator in the regression of Y on $\mathbf{\text{X}}\in {\mathbb{R}}^p$ and the one-dimensional dimension reduction subspace can be derived. Specifically, consider the following bivariate loss function, (1) $\begin{array}{rl}& L(a+{\mathbf{b}}^T\mathbf{X},Y),L(u,v)\mathrm{i}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{x}\mathrm{i}\mathrm{n}u,\\ & \mathrm{e}.\mathrm{g}.L(u,v)=(u-v{)}^2;\end{array}$(1) and define the unique minimiser as follows, (2) $(\alpha ,\text{β})\equiv \arg \underset{a,\mathbf{b}\in {\mathbb{R}}^p}{min}\mathrm{E}\left\{L\right(a+{\mathbf{b}}^T\mathbf{X},Y\left)\right\}.$(2) Then under the linearity condition of dimension reduction subspace, it was shown that $\text{β}$ is contained in the central subspace (or any dimension reduction subspace).

A direct consequence of this somewhat surprising result is that one can simply use OLS to extract the direction in nonlinear models of the following form, (3) $Y=g({\mathit{\gamma }}^T\mathbf{X},\epsilon ),$(3) where g is some unknown function, ϵ is the error term, and $\mathit{\gamma }\in {\mathbb{R}}^{p\times d}$ spans the central subspace under this model. The model (3(3) $Y=g({\mathit{\gamma }}^T\mathbf{X},\epsilon ),$(3) ) is known as the single-index model when d = 1, and the multiple-index model when d>1. Then Li–Duan Theorem (Li & Duan, 1989, Theorem 2.1) implies that the solution $\text{β}$ from ordinary least squares estimation, or more generally from (2(2) $(\alpha ,\text{β})\equiv \arg \underset{a,\mathbf{b}\in {\mathbb{R}}^p}{min}\mathrm{E}\left\{L\right(a+{\mathbf{b}}^T\mathbf{X},Y\left)\right\}.$(2) ), lies within the central subspace: $\text{β}\in \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\left(\mathit{\gamma }\right)={\mathcal{S}}_{Y\mid \mathbf{\text{X}}}$.

For single-index models, this means that $\text{β}$ from OLS fitting, $(\alpha ,\text{β})=\arg \underset{a,\mathbf{\text{b}}}{min}\mathrm{E}\left\{\right(Y-a-{\mathbf{\text{b}}}^T\mathbf{\text{X}}{)}^2\}$, is the same as $\mathit{\gamma }$ in population up to a scalar multiplication (i.e. $\text{β}=c\cdot \mathit{\gamma }$ for some constant $c\ne 0$). As such, we might still use OLS to study regression graphics even when there is a nonlinear relationship between Y and $\mathbf{\text{X}}$ (Cook, 1998). In high-dimensional setting, one may replace OLS with its penalised version such as LASSO regression and achieve consistent variable selection and directional estimation (Neykov, Lin, et al., 2016; Neykov, Liu, et al., 2016).

For multiple-index models and model-free sufficient dimension reduction, the (sparse) estimation of the central subspace is much more challenging. In high-dimensional sparse sufficient dimension reduction, it is thus desirable to have a penalised least squares formulation. Indeed, as this paper points out, the computationally tractable and rate optimal sparse SIR is eventually obtained by the adaptive estimation scheme based on a least squares formulation (Tan et al., 2020). Finally, the paper concludes that ‘Since SAVE and DR can not be rewritten as a least-square formulation, we do not define refined sparse SAVE and DR estimator’. Not surprisingly, developing theoretical solid and computationally feasible methods for high-dimensional SAVE and DR is very challenging and requires substantial efforts.

2. Sparse/constrained canonical correlation

As discussed in Section 3.1 of this paper, sparse sufficient dimension reduction subspace can be obtained by the C${}^3$ method: constrained canonical correlation (Zhou and He, 2008). The idea is to estimate the constrained canonical variates between the B-spline basis functions of response transformation, $\pi \left(\mathbf{\text{Y}}\right)\in {\mathbb{R}}^{m+k_n}$, where m is the spline order and $k_n$ is the number of interval knots in B-spline transformation, and the predictor $\mathbf{\text{X}}\in {\mathbb{R}}^p$. Then this procedure can be viewed as an estimation for sparse SIR directions. However, as the authors noted, this method may not be directly applicable to very high-dimensional settings, at least in theory. In the past few years, there are some advances in both the theoretical and computational aspects of high-dimensional canonical correlation analysis. In Mai and Zhang (2019), they solved the sparse canonical correlation analysis (CCA) problem using an iterative penalised least squares approach. This new iterative penalised least squares approach can be directly applied to estimate the sparse sufficient dimension reduction directions when combined with B-spline transformations of the response.

To illustrate the idea, we consider the estimation of the leading CCA directions. For a multivariate $\mathbf{\text{X}}\in {\mathbb{R}}^p$ and a multivariate $\mathbf{\text{Y}}\in {\mathbb{R}}^q$, the leading CCA directions are defined through a pair of linear combinations ${\mathit{\alpha }}_1^T\mathbf{\text{Y}}$ and ${\text{β}}_1^T\mathbf{\text{X}}$ such that the correlation between them are maximised. When $max(p,q)\gg n$, the sparse CCA problem assumes that the population solution ${\mathit{\alpha }}_1^{\star }\in {\mathbb{R}}^q$ and ${\text{β}}_1^{\star }\in {\mathbb{R}}^p$ are both sparse so that we can estimate them with a limited sample size. Then the leading sparse CCA directions can be obtained by solving the following constrained optimisation problem, (4) $\begin{array}{rl}(\hat{\mathit{\alpha }},\widehat{\text{β}})& =\arg \underset{\mathit{\alpha }\in {\mathbb{R}}^q,\text{β}\in {\mathbb{R}}^p}{min}\sum _{i=1}^n({\mathit{\alpha }}^T{\mathbf{Y}}_i-{\text{β}}^T{\mathbf{X}}_i{)}^2\\ & +{\lambda }_{\alpha }\parallel \mathit{\alpha }{\parallel }_1+{\lambda }_{\alpha }\parallel \text{β}{\parallel }_1,\end{array}$(4) subject to constraints that ${\mathit{\alpha }}^T{\hat{\mathbf{\Sigma }}}_{\mathbf{\text{Y}}}\mathit{\alpha }=1$ and ${\text{β}}^T{\hat{\mathbf{\Sigma }}}_{\mathbf{\text{X}}}\text{β}=1$. Note that the data are centred so that $\sum _{i=1}^n{\mathbf{\text{X}}}_i=0$ and $\sum _{i=1}^n{\mathbf{\text{Y}}}_i=0$, also that the sample covariance ${\hat{\mathbf{\Sigma }}}_{\mathbf{\text{Y}}}$ and ${\hat{\mathbf{\Sigma }}}_{\mathbf{\text{X}}}$ do not need to be positive definite. Then it can be shown that, the above optimisation can be solved by iteratively solving the following two LASSO regression problems. Specifically, the sparse CCA solutions can be obtained as follows, (5) $\begin{array}{rl}\tilde{\mathit{\alpha }}& =\arg \underset{\mathit{\alpha }\in {\mathbb{R}}^q}{min}\sum _{i=1}^n({\widehat{\text{β}}}_1^T{\mathbf{X}}_i-{\mathit{\alpha }}^T{\mathbf{Y}}_i{)}^2+{\lambda }_{\alpha }\parallel \mathit{\alpha }{\parallel }_1,\\ {\hat{\mathit{\alpha }}}_1& =\frac{\tilde{\mathit{\alpha }}}{\sqrt{{\tilde{\mathit{\alpha }}}^T{\hat{\mathbf{\Sigma }}}_{\mathbf{Y}}\tilde{\mathit{\alpha }}}};\end{array}$(5) (6) $\begin{array}{rl}\widetilde{\text{β}}& =\arg \underset{\text{β}\in {\mathbb{R}}^p}{min}\sum _{i=1}^n({\hat{\mathit{\alpha }}}_1^T{\mathbf{Y}}_i-{\text{β}}^T{\mathbf{X}}_i{)}^2+{\lambda }_{\text{β}}\parallel \text{β}{\parallel }_1,\\ {\widehat{\text{β}}}_1& =\frac{\widetilde{\text{β}}}{\sqrt{{\widetilde{\text{β}}}^T{\hat{\mathbf{\Sigma }}}_{\mathbf{X}}\widetilde{\text{β}}}}.\end{array}$(6) For sequential sparse CCA directions, ${\hat{\mathit{\alpha }}}_k\in {\mathbb{R}}^q$ and ${\widehat{\text{β}}}_k\in {\mathbb{R}}^p$, $k=1,2,\dots ,$ we can use the similar iterative penalised least squares formulation after deflation of the data from the previous $(k-1)$ estimated directions. Theoretical results show that this approach can consistently estimate the population directions (any fix number of pairs $k=1,\dots ,K$) with an overwhelming probability in ultra-high dimensions $\log (p+q)=o\left(n\right)$. More importantly, due to its simplicity, the iterative penalised least squares approach can be extremely fast and scalable – even much faster than some existing convex formulations. This approach might be useful for sparse sufficient dimension reduction when the response is multivariate and even high-dimensional.

Disclosure statement

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

Additional information

Notes on contributors

Xin Zhang

Xin Zhang is an Associate Professor in Statistics at the Florida State University.