Sufficient dimension reduction (SDR) methods aim to identify the central subspace that preserves all the information in the predictors. In this paper, we consider the SDR problem of a multivariate response $\mathbf{Y}\in\mathbb{R}^{q}$ on a multivariate predictor $\mathbf{X}\in\mathbb{R}^{p}$. The central subspace $\cal S_{\mathbf{Y}\mid\mathbf{X}}$ is defined as the intersection of all the subspaces $\cal S$ such that $\mathbf{Y}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\mathbf{X}\mid\mathbf{P}_{\cal S}\mathbf{X}$, where $\mathbf{P}_{\cal S}$ is the projection matrix onto $\cal S$. By construction, the central subspace $\cal S_{\mathbf{Y}\mid\mathbf{X}}$ is the smallest dimension reduction subspace that contains all the information in the conditional distribution of $\mathbf{Y}$ given $\mathbf{X}$. Many methods have been proposed for the recovery of the central subspace or a portion of the central subspace (Li, 1991; Cook and Weisberg, 1991; Bura and Cook, 2001; Chiaromonte et al., 2002; Yin and Cook, 2003; Cook and Ni, 2005; Li and Wang, 2007; Zhou and He, 2008). See Li (2018) for a comprehensive review. Although the central subspace is well defined for both univariate and multivariate response, most existing SDR methods consider the case with univariate response, while extension to multivariate response is non-trivial.

The definition of central subspace is not very constructive, as it requires taking intersection of all subspace $\mathcal{S}\subseteq\mathbb{R}^{p}$ such that $\mathbf{Y}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\mathbf{X}\mid\mathbf{P}_{\cal S}\mathbf{X}$. It is indeed a very ambitious goal to estimate the central subspace without specifying a model between $\mathbf{Y}$ and $\mathbf{X}$. To achieve this, we often need additional assumptions such as the linearity and the coverage conditions. The linearity condition requires that, for any basis of the central subspace $\boldsymbol{\beta}$, we must have that $\mathrm{E}(\mathbf{X}\mid\boldsymbol{\beta}^{\mathrm{\tiny T}}\mathbf{X})$ is linear in $\boldsymbol{\beta}^{\mathrm{\tiny T}}\mathbf{X}$. The linearity condition is guaranteed if $\mathbf{X}$ is elliptically contoured, and allows us to connect the central subspace to the conditional expectation $\mathrm{E}(\mathbf{X}\mid\mathbf{Y})$. Define $\boldsymbol{\Sigma}_{\mathbf{X}}$ as the covariance of $\mathbf{X}$ and the inverse regression subspace

Then following property is well-known and often adopted in developing SDR methods.

The coverage condition further assumes that $\mathcal{S}_{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})}=\boldsymbol{\Sigma}_{\mathbf{X}}\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}$. It follows that we can estimate the central subspace by modeling the conditional expectation of $\mathbf{X}$. Indeed, many SDR methods approximate $\mathrm{E}(\mathbf{X}\mid\mathbf{Y})$. For example, the most classical SDR method, sliced inverse regression (SIR), slices univariate $Y$ into several categories and estimate the mean of $\mathbf{X}$ within each slice. Most later methods also follow this slice-and-estimate procedure. Apparently, the slicing scheme is important to the estimation. If there are too few slices, we may not be able to fully capture the dependence of $\mathbf{X}$ on $Y$; however, if there are too many slices, there is a lack of enough samples within each slice to allow accurate estimation.

In this section, we lay the foundation for the application of MDDM in SDR. We show that the subspace spanned by MDDM coincides with the inverse regression subspace in (3.3). In particular, we have the following Proposition 2, which was used in Zhang et al. (2020), without a proof, in the context of multivariate linear regression.

Therefore, the rank of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$ is the dimensionality of the inverse regression subspace; and the non-trivial eigenvectors of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$ contain all the information for $\mathcal{S}_{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})}$. Combining Proposition1& 2, we immediately have that (i) under the linearity condition, $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}\subseteq\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}$; and (ii) under the linearity and coverage conditions, $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}=\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}$.

Henceforth, we assume both the linearity and coverage conditions, which are assumed either explicitly or implicitly in inverse regression type dimension reduction methods (e.g., Li, 1991; Cook and Ni, 2005; Zhu et al., 2010; Cook and Zhang, 2014). Then the central subspace is related to the eigen-decomposition of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$. Specifically, we have the following scenarios.

If $\mathrm{cov}(\mathbf{X})=\sigma^{2}\mathbf{I}_{p}$ for some $\sigma^{2}>0$, then obviously $\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}=\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}$. This includes single index and multiple index models with uncorrelated predictors. Let $K$ be the rank of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$, then the dimension of the central subspace is $K$; and the first $K$ eigenvectors of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$ span the central subspace.

If $\mathrm{cov}(\mathbf{X}\mid\mathbf{Y})=\sigma^{2}\mathbf{I}_{p}$ for some $\sigma^{2}>0$, then we have $\boldsymbol{\Sigma}_{\mathbf{X}}=\sigma^{2}\mathbf{I}_{p}+\mathrm{cov}\{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})\}$. Because $\mathrm{span}[\mathrm{cov}\{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})\}]=\mathcal{S}_{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})}$, we can show that $\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}=\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}=\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}$. To see this, let $\mathrm{cov}\{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})\}=\mathbf{U}\mathbf{U}^{\mathrm{\tiny T}}$ for some $\mathbf{U}\in\mathbb{R}^{p\times K}$, then $\mathrm{span}(\mathbf{U})=\mathrm{span}[\mathrm{cov}\{\mathrm{E}(\mathbf{X}\mid\mathbf{Y})\}]=\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}$ and we may also write $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})=\mathbf{U}\boldsymbol{\Psi}\mathbf{U}^{{\mathrm{\tiny T}}}$ for some symmetric positive definite matrix $\boldsymbol{\Psi}\in\mathbb{R}^{K\times K}$. Then the result follows by applying the Woodbury matrix identity to $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}=(\sigma^{2}\mathbf{I}_{p}+\mathbf{U}\mathbf{U}^{\mathrm{\tiny T}})^{-1}=\sigma^{-2}\mathbf{I}_{p}-\sigma^{-2}\mathbf{U}(\sigma^{2}\mathbf{I}_{K}+\mathbf{U}^{{\mathrm{\tiny T}}}\mathbf{U})^{-1}\mathbf{U}^{{\mathrm{\tiny T}}}$. The non-trivial eigenvectors of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$ again span the central subspace.

For general covariance structure, the $d$-dimensional central subspace $\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}=\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}\mathrm{span}\{\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})\}$ can be obtained via generalized eigen-decomposition. Specifically, consider the generalized eigenvalue problem

where $\mathbf{v}_{i}^{\mathrm{\tiny T}}\boldsymbol{\Sigma}_{\mathbf{X}}\mathbf{v}_{j}=0$ for $i\neq j$. Then, similar to (Li, 2007; Chen et al., 2010), it is straightforward to show that the generalized eigenvector spans the central subspace, $\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}=\mathrm{span}(\mathbf{v}_{1},\dots,\mathbf{v}_{K})$.

Existing works in SDR often focus on the eigen-decomposition or the generalized eigen-decomposition of $\mathrm{cov}\{\mathrm{E}(\mathbf{X}\mid Y=y)\}$, where non-parametric estimates of $\mathrm{E}(\mathbf{X}\mid Y=y)$ are obtained from slicing the support of the univariate response $Y$. Comparing to these approaches, the MDDM approach requires no tuning parameter selection (i.e. specifying slicing schemes). Moreover, high-dimensional theoretical study of MDDM is easier and does not require additional assumptions on the conditional mean function $\mathrm{E}(\mathbf{X}\mid\mathbf{Y})$ such as smoothness in the empirical mean function of $\mathbf{X}$ given $Y$ (e.g. sliced stable condition in Lin et al. (2018)).

So far, we have discussed model-free SDR. Another important research area in SDR is model-based methods, which provide invaluable intuition for the use of inverse regression estimation under the assumption that the conditional distribution of $\mathbf{X}\mid\mathbf{Y}$ is normal. In this section, we consider the principal fitted component (PFC) model, which was discussed in details in Cook and Forzani (2009) and Cook (2007), and generalize it from univariate response to multivariate response. We argue that (generalized) eigen-decomposition of MDDM is potentially advantageous to the likelihood-based approaches under PFC model. This is somewhat surprising but reasonable, considering that the advantages of MDDM over least squares and likelihood-based estimation was demonstrated in Zhang et al. (2020) under the multivariate linear model.

Let $\mathbf{X}_{\mathbf{y}}\sim\mathbf{X}\mid(\mathbf{Y}=\mathbf{y})$ denote the conditional variable, then the PFC model is

where $\boldsymbol{\Gamma}\in\mathbb{R}^{p\times K}$, $K<p$, is a non-stochastic orthogonal matrix, $\boldsymbol{\nu}_{\mathbf{y}}\in\mathbb{R}^{K}$ is the latent variable that depends on $\mathbf{y}$. Then the latent variable $\boldsymbol{\nu}_{\mathbf{y}}$ is fitted as $\boldsymbol{\nu}_{\mathbf{y}}=\boldsymbol{\alpha}\mathbf{f}_{\mathbf{y}}$ with some user-specified functions $\mathbf{f}_{\mathbf{y}}=(f_{1}(\mathbf{y}),\ldots,f_{m}(\mathbf{y}))^{{\mathrm{\tiny T}}}\in\mathbb{R}^{m}$, $m\geq K$, that maps $q$-dimensional response to $m$-dimensional. In the univariate PFC model, $q=1$, so the $m$ functions can be viewed as an expansion of the response (similar to slicing). For our multivariate extensions of the PFC model, there is no requirement of $m\geq q$. The PFC model can be written as

where $\boldsymbol{\Gamma}$ and $\boldsymbol{\alpha}$ are estimated similarly to the multivariate reduced-rank regression with $\mathbf{X}\in\mathbb{R}^{p}$ being the response and $\mathbf{f}_{\mathbf{y}}\in\mathbb{R}^{m}$ being the predictor. Finally, the central subspace under this PFC model is $\boldsymbol{\Delta}^{-1}\mathrm{span}(\boldsymbol{\Gamma})$, which simplifies to $\mathrm{span}(\boldsymbol{\Gamma})$ if we further assume isotropic error (i.e. isotropic PFC model) $\boldsymbol{\Delta}=\mathrm{cov}(\mathbf{X}\mid\mathbf{Y})=\sigma^{2}\mathbf{I}_{p}$.

For the PFC model, our MDDM approach is the same as the model-free MDDM counterpart, and has two main advantages over the likelihood-based PFC estimation: (i) there is no need to specify the functions $\mathbf{f}_{\mathbf{y}}$, and thus no risk of mis-specification; (ii) extensions to high-dimensional setting is much more straightforward. Moreover, under the isotropic PFC model, the central subspace $\mathcal{S}_{\mathbf{Y}\mid\mathbf{X}}=\mathrm{span}(\boldsymbol{\Gamma})$ is exactly the first $K$ eigenvectors of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$.

Based on the results in the last section, penalized eigen-decomposition of MDDM can be used for estimating the central subspace in high dimension when the covariance $\boldsymbol{\Sigma}_{\mathbf{X}}$ or the conditional covariance $\mathrm{cov}(\mathbf{X}\mid\mathbf{Y})$ is proportional to the identity matrix $\mathbf{I}_{p}$. We consider the construction of such an estimate. It is worth mentioning that the penalized decomposition of MDDM we develop here is immediately applicable to the dimension reduction of multivariate stationary time series in (Lee and Shao, 2018), which is beyond the scope of this article. Moreover, it is well-known that $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}$ is not easy to estimate in high dimensions. Then, even when for general covariance structure, the eigen-decomposition of MDDM provides estimate of the inverse regression subspace (though may differ from the central subspace) that is useful for exploratory data analysis (e.g. detecting and visualizing non-linear mean function).

As such, we consider the estimation of the eigenvectors of $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$. We assume that $\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$ has $K$ nontrivial eigenvectors, denoted by $\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{K}$. We use the shorthand notation $\mathbf{M}=\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y})$. Also, we note that, given the first $k-1$ eigenvectors, $\boldsymbol{\beta}_{k}$ is the top eigenvector of $\mathbf{M}_{k}$, where $\mathbf{M}_{k}=\mathbf{M}-\sum_{l<k}(\boldsymbol{\beta}_{l}^{\mathrm{\tiny T}}\mathbf{M}\boldsymbol{\beta}_{l})\boldsymbol{\beta}_{l}\boldsymbol{\beta}_{l}^{\mathrm{\tiny T}}$.

It is well-known that the eigenvectors cannot be accurately estimated in high dimensions without additional assumptions. We adopt the popular sparsity assumption that many entries in $\boldsymbol{\beta}_{k}$ are zero. To estimate these sparse eigenvectors, denote $\widehat{\mathbf{M}}_{1}=\mathrm{MDDM}_{n}(\mathbf{X}\mid\mathbf{Y})$, where the sample $\mathrm{MDDM}_{n}$ is defined in (2.1). We find $\widehat{\boldsymbol{\beta}}_{k},k=1,\ldots,K$ as follows:

where $\widehat{\mathbf{M}}_{1}=\mathrm{MDDM}_{n}(\mathbf{X}\mid\mathbf{Y})$, $\widehat{\mathbf{M}}_{k}=\widehat{\mathbf{M}}_{1}-\sum_{l<k}\delta_{l}\widehat{\boldsymbol{\beta}}_{l}\widehat{\boldsymbol{\beta}}_{l}^{\mathrm{\tiny T}}$ for $k>1$ with $\delta_{l}=\widehat{\boldsymbol{\beta}}_{l}^{\mathrm{\tiny T}}\widehat{\mathbf{M}}_{1}\widehat{\boldsymbol{\beta}}_{l}$, and $s$ is a tuning parameter.

We solve the above problem by combining the truncated power method with hard thresholding. For a vector $\mathbf{v}\in\mathbb{R}^{p}$ and a positive integer $s$, denote $v_{s}^{*}$ as the $s$-th largest value of $|v_{j}|,j=1,\ldots,p$. The hard-thresholding operator is $\mathrm{HT}(\mathbf{v},s)=(v_{1}I(|v_{1}|\geq v_{s}^{*}),\ldots,v_{p}I(|v_{p}|\geq v_{s}^{*}))^{\mathrm{\tiny T}}$, which sets the $p-s$ elements in $\mathbf{v}$ to zero. We solve (4.7) by Algorithm 1, where the initialization $\widehat{\boldsymbol{\beta}}_{1}^{(0)}$ may be randomly generated. Note that Yuan and Zhang (2013) proposed Algorithm 1 to perform principal component analysis through penalized eigen-decomposition on the sample covariance.

In our algorithm, we require a pre-specified sparsity level $s$ and subspace dimension $K$. In theory, we show that our estimators for $\boldsymbol{\beta}_{k}$, $k=1,\dots,K$, are all consistent for their population counterparts when the sparsity $s$ is sufficiently large (i.e. larger than the population sparsity level) and the number of directions $K$ is no bigger than the true dimension of the central subspace. Therefore, our method is flexible in the sense that the pre-specified $s$ and $K$ do not have to be exactly correct. In practice, especially in exploratory data analysis, the number of sequentially extracted directions are often set to be small (i.e. $K=1,2\leavevmode\nobreak\ \mbox{or}\leavevmode\nobreak\ 3$), while the determination of true central subspace dimension is a separate and important research topic in SDR (e.g. Bura and Yang, 2011; Luo and Li, 2016) and is beyond the scope of this paper. Moreover, the pre-specified sparsity level $s$ combined with $\ell_{0}$ regularization is potentially convenient for post-dimension reduction inference (Kim et al., 2020), as seen in the post-selection inference of canonical correlation analysis that is done over subsets of variables of pre-specified cardinality (McKeague and Zhang, 2020).

As pointed out by a referee, other sparse principal component analysis (PCA) methods can potentially be applied to decompose MDDM. We choose to extend the algorithm in Yuan and Zhang (2013) to facilitate computation and theoretical development. For computationally efficient sparse PCA methods such as Zou et al. (2006); Witten et al. (2009), their theoretical properties are unfortunately unknown. Hence, we expect the theoretical study of their MDDM-variants to be very challenging. On the other hand, for the theoretically justified sparse PCA methods such as Vu and Lei (2013); Cai et al. (2013), the computation is less efficient.

Now we consider the general (arbitrary) covariance structure $\boldsymbol{\Sigma}_{\mathbf{X}}$. We continue to use $\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{K}$ to denote the nontrivial eigenvectors of $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}\mathrm{span}(\mathrm{MDDM}(\mathbf{X}\mid\mathbf{Y}))$ so that the central subspace is spanned by the $\boldsymbol{\beta}$’s. Again, we assume that these eigenvectors are sparse. In principle, we could assume that $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}$ is also sparse and construct its estimate accordingly. However, $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}$ is a nuisance parameter for our ultimate goal, and additional assumptions on it may unnecessarily limit the applicability of our method. Hence, we take a different approach as follows.

To avoid estimating $\boldsymbol{\Sigma}_{\mathbf{X}}^{-1}$, we note that $\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{K}$ can also be viewed as the generalized eigenvectors defined as follows, which is equivalent to (3.4),

Directly solving the generalized eigen-decomposition problem in (4.8) is not easy if we want to further impose penalties, because it is difficult to satisfy the orthogonality constraints. Therefore, we further consider another form for (4.8) that does not involve the orthogonal constraints. This alternative form is based on the following lemma.

Motivated by Lemma 1, we consider the penalized problem that $\boldsymbol{\beta}_{k}=\arg\max_{\boldsymbol{\beta}}\boldsymbol{\beta}^{\mathrm{\tiny T}}\widehat{\mathbf{M}}_{k}\boldsymbol{\beta}$ such that $\boldsymbol{\beta}^{\mathrm{\tiny T}}\boldsymbol{\Sigma}_{\mathbf{X}}\boldsymbol{\beta}=1,\|\boldsymbol{\beta}\|_{0}\leq s$, where $\widehat{\mathbf{M}}_{1}=\mathrm{MDDM}_{n}(\mathbf{X}\mid\mathbf{Y})$ and $\widehat{\mathbf{M}}_{k}=\widehat{\mathbf{M}}_{1}-\widehat{\boldsymbol{\Sigma}}_{\mathbf{X}}\left(\sum_{l<k}\delta_{l}\widehat{\boldsymbol{\beta}}_{l}\widehat{\boldsymbol{\beta}}_{l}^{\mathrm{\tiny T}}\right)\widehat{\boldsymbol{\Sigma}}_{\mathbf{X}}$ for $k>1$ with $\delta_{l}=\widehat{\boldsymbol{\beta}}_{l}^{\mathrm{\tiny T}}\widehat{\mathbf{M}}\widehat{\boldsymbol{\beta}}_{l}$, and $s$ is a tuning parameter. We adopt the RIFLE algorithm in Tan, Wang, Liu and Zhang (2018) to solve this problem. See the details in Algorithm 2. In our simulation studies, we considered randomly generated initial value $\widehat{\boldsymbol{\beta}}_{1}^{(0)}$ and fixed step size $\eta=1$, and observed reasonably good performance.

Although Algorithm 2 is a generalization of the RIFLE Algorithm in Tan, Wang, Liu and Zhang (2018), there are important differences between these two. On one hand, the RIFLE Algorithm only extracts the first generalized eigenvector, whereas Algorithm 2 is capable of estimating multiple generalized eigenvectors by properly deflating the MDDM. In sufficient dimension reduction problems, the central subspace often has a structural dimension greater than 1, and it is necessary to find more than one generalized eigenvector. Hence, Algorithm 2 is potentially more useful than the RIFLE algorithm in practice. On the other hand, the usefulness of RIFLE Algorithm has been demonstrated in several statistical applications, including sparse sliced inverse regression. Here Algorithm 2 decomposes the MDDM, which is the first time the penalized generalized eigenvector problem is used to perform sufficient dimension reduction in a slicing-free manner in high dimensions. A brief analysis of the computation complexity is included in Section S3 in the Supplementary Materials.

We compare our slicing-free approaches to the state-of-the-art high-dimensional extensions of sliced inverse regression estimators. Both univariate response and multivariate response settings are considered. Specifically, for univariate response simulations, we include Rifle-SIR (Tan, Wang, Liu and Zhang, 2018) and Lasso-SIR (Lin et al., 2019) as two main competitors; for multivariate response simulations, we mainly compare our method with the projective resampling approach to SIR (PR-SIR, Li et al., 2008), which is a computationally expensive method that repeatedly projects the multivariate response to one-dimensional subspaces. For Rifle-SIR, we adopt the Rifle algorithm to estimate the leading eigenvector of the sample matrix $\mathrm{cov}\{\mathrm{E}(\mathbf{X}\mid Y)\}$ based on slicing. In addition, we include the oracle-SIR as a benchmark method, where we perform SIR on the subset of truly relevant variables (hence a low-dimensional estimation problem). For all these SIR-based methods, we include two different slicing schemes by setting the number of slices to be $3$ and $10$, where $3$ is the minimal number of slices required to obtain our two-dimensional central subspace and $10$ is a typical choice used in the literature. To evaluate the performances of these SDR methods, we use the subspace estimation error defined as $\mathcal{D}(\widehat{\boldsymbol{\beta}},\boldsymbol{\beta})=\|\mathbf{P}_{\widehat{\boldsymbol{\beta}}}-\mathbf{P}_{\boldsymbol{\beta}}\|_{F}/\sqrt{2K}$, where $\widehat{\boldsymbol{\beta}},\boldsymbol{\beta}\in\mathbb{R}^{p\times K}$ are the estimated and the true basis matrices of the central subspace and $\mathbf{P}_{\widehat{\boldsymbol{\beta}}},\mathbf{P}_{\boldsymbol{\beta}}\in\mathbb{R}^{p\times p}$ are the corresponding projection matrices. This subspace estimation error is always between $0$ and $1$, and a small value indicates a good estimation.

First, we consider the following six models for univariate response regression: $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are single index models (i.e. $K=1$), $\mathcal{M}_{3}$–$\mathcal{M}_{5}$ are multiple index models (i.e. $K=2$), that are widely used in SDR literature; $\mathcal{M}_{6}$ is isotropic PFC model with $K=1$. Specifically,

where $\mathbf{X}\sim N_{p}(0,\boldsymbol{\Sigma}_{\mathbf{X}})$ and $\epsilon\sim N(0,1)$ for $\mathcal{M}_{1}$–$\mathcal{M}_{5}$, and $Y\sim N(0,1),\boldsymbol{\epsilon}\sim N_{p}(0,\mathbf{I}_{p})$ for the isotropic PFC model ($\mathcal{M}_{6}$). The sparse directions in the central subspace $\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2}\in\mathbb{R}^{p}$ are orthogonal as we let the first $s=6$ elements in $\boldsymbol{\beta}_{1}$ and the $6$-th to $12$-th elements in $\boldsymbol{\beta}_{2}$ to be $1/\sqrt{6}$ (while all other elements are zero). For $\mathcal{M}_{1}$–$\mathcal{M}_{5}$, we consider both the independent predictor setting with $\boldsymbol{\Sigma}_{\mathbf{X}}=\mathbf{I}_{p}$ and the correlated predictor setting with auto-regressive correlation that $\Sigma_{X}(i,j)=0.5^{|i-j|}$ for $i,j=1,2,...,p$. For each of model settings, we vary the sample size $n\in\{200,500,800\}$ and predictor dimension $p\in\{200,500,800,1200,2000\}$ and simulate 1000 independent data sets.

For our method, we applied the generalized eigen-decomposition algorithm (Algorithm 2) in all these six models (even when the covariance $\mathbf{X}$ is identity matrix). In the single index models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, we use the random initialization ($\widehat{\boldsymbol{\beta}}^{(0)}$ is randomly generated from $p$-dimensional standard normal) for our algorithm and Rifle-SIR to demonstrate the robustness to initialization. The step size in the algorithm is simply fixed as $\eta=1$. For more challenging multiple index models, $\mathcal{M}_{3}-\mathcal{M}_{6}$, we consider the best case scenarios for each method, therefore true parameter $\boldsymbol{\beta}$ is used as the initial value and an optimal $\eta\in\{0.1,0.2,\dots,1.0\}$ is selected from a separate training sample with 400 observations. The results based on 1000 replications for $n=200$ and $p=800$ are summarized in Table 1, while the rest of the results can be found in the Supplemental Materials. Overall, the slicing-free MDDM approach is much more accurate than existing SIR-based methods. It is almost as accurate as the oracle-SIR. Moreover, it is clear that SIR-type methods are rather sensitive to the choice of the number of slices.