Sparse kernel sufficient dimension reduction

Abstract

The sufficient dimension reduction (SDR) with sparsity has received much attention for analysing high-dimensional data. We study a nonparametric sparse kernel sufficient dimension reduction (KSDR) based on the reproducing kernel Hilbert space, which extends the methodology of the sparse SDR based on inverse moment methods. We establish the statistical consistency and efficient estimation of the sparse KSDR under the high-dimensional setting where the dimension diverges as the sample size increases. Computationally, we introduce a new nonconvex alternating directional method of multipliers (ADMM) to solve the challenging sparse SDR and propose the nonconvex linearised ADMM to solve the more challenging sparse KSDR. We study the computational guarantees of the proposed ADMMs and show an explicit iteration complexity bound to reach the stationary solution. We demonstrate the finite-sample properties in simulation studies and a real application.

KEYWORDS:

• Sufficient dimension reduction
• kernel dimension reduction
• nonconvex optimisation
• alternating direction method of multipliers

1. Introduction

Consider the sufficient dimension reduction (SDR) of high-dimensional random vector $\mathit{x}\in {\mathbb{R}}^p$ in the presence of response $y\in \mathbb{R}$. We aim to find the low-dimensional central subspace ${\mathcal{S}}_{y|\mathit{x}}\subset {\mathbb{R}}^p$ such that $\dim \left({\mathcal{S}}_{y|\mathit{x}}\right)=d\ll p$ and $y\perp \perp \mathit{x}|P_{{\mathcal{S}}_{y|\mathit{x}}}\left(\mathit{x}\right)$, where the dimension p may diverge as the sample size increases, and $P_{{\mathcal{S}}_{y|\mathit{x}}}(\cdot )$ is the projection of a vector onto the subspace ${\mathcal{S}}_{y|\mathit{x}}$. The inverse moment methods are extensively studied under the fixed-dimension setting, including the sliced inverse regression (SIR) (Li 1991), sliced average variance estimation (SAVE) (Cook and Weisberg 1991), and directional regression (DR) (Li and Wang 2007). Fan, Xue, and Yao (2017), Luo, Xue, Yao, and Yu (2022), and Yu, Yao, and Xue (2022) extended these inverse moment methods and introduced sufficient forecasting using factor models for high-dimensional predictors when both $n,p\to \mathrm{\infty }$. Ying and Yu (2022), Zhang, Xue, and Li (2024), and Zhang, Li, and Xue (2024) studied the sufficient dimension reduction for Fréchet regression and distribution-on-distribution regression. (Lin, Zhao, and Liu 2018; Lin, Zhao, and Liu 2019) introduced the Lasso-SIR under the high-dimensional setting when $p/n\to 0$ and the diagonal thresholding SIR (DT-SIR) introduced a new diagonal thresholding screening under the ultra high-dimensional setting when $p/n\to \mathrm{\infty }$. However, the restricted linearity condition on the conditional expectation of $\mathit{x}$ given $P_{{\mathcal{S}}_{y|\mathit{x}}}\left(\mathit{x}\right)$ is required for inverse moment methods. To remove the linearity condition, the kernel sufficient dimension reduction (KSDR) was proposed by Fukumizu, Bach, and Jordan (2009). The nonconvex KSDR does not require the linearity condition nor the elliptical distribution of $\mathit{x}$, but it poses a significant computational challenge. Wu, Miller, Chang, Sznaier, and Dy (2019) solved this problem using an iterative spectral method, but such a method is hard to extend when a sparsity assumption is assumed in high-dimensional settings. Shi, Huang, Feng, and Suykens (2019) studies nonconvex penalty in a kernel regression problem. Although both the theoretical property and optimisation algorithm cannot be applied in SDR settings, it shows the promises of nonconvex penalised KSDR in high-dimensional settings.

The sparse estimation of central subspace has received considerable attention in the past decade. The main focus of the literature is on sparse SDR based on inverse moment methods, which is equivalent to solving the generalised eigenvectors from the estimated kernel matrix (Li 2007). For example, the kernel matrix for SIR is the conditional covariance matrix of $\mathit{x}$ given y. In view of such equivalence, Li (2007) followed the similar spirit of the sparse principal component analysis (SPCA) (Zou, Hastie, and Tibshirani 2006; Zou and Xue 2018) to propose the sparse SDR using ${\ell }_1$ penalisation. Bondell and Li (2009) proposed a shrinkage inverse regression estimation and studied the asymptotic consistency and asymptotic normality. Chen, Zou, and Cook (2010) proposed a coordinate-independent sparse estimation method and established the oracle property when the dimension p is fixed. Lin et al. (2019) established the consistency of the Lasso-SIR by assuming the number of slices also diverges. Neykov, Lin, and Liu (2016) followed the convex semidefinite programming (d'Aspremont, Ghaoui, Jordan, and Lanckriet 2005; Amini and Wainwright 2008) to give a high-dimensional analysis for the sign consistency of the DT-SIR in a single-index model. Tan, Wang, Zhang, Liu, and Cook (2018) proposed a sparse SIR estimator based on convex formulation and established the Frobenius norm convergence of projection subspace. Lin et al. (2018) established the consistency of DT-SIR by assuming a stability condition for $E\left(\mathit{x}|y\right)$. Tan, Shi, and Yu (2020) proved the optimality of ${\ell }_0$-constrained sparse solution and provided a three-steps adaptive estimation to approximate $\ell -0$ solution with near-optimal property. The sparse estimation in the KSDR area has not been well studied.

In this paper, we propose the nonconvex estimation scheme for sparse KSDR, motivated by the nonconvex estimation of sparse SDR (Li 2007). The proposed methods target the ideal ${\ell }_0$-constrained estimator and can be efficiently solved with statistical and computational guarantees after approximation. We follow d'Aspremont et al. (2005) to reformulate sparse SDR and sparse KSDR as the ${\ell }_0$-constrained trace maximisation problem. It is worth pointing out that KSDR is also equivalent to solving the eigenvectors from the nonconvex residual covariance estimator in RKHS (Fukumizu et al. 2009). The residual covariance estimator in RKHS plays a similar role as the kernel matrix in inverse moment methods. Thus, we propose that sparse SDR solves the ${\ell }_0$-constrained convex M-estimation problem, whereas sparse KSDR solves the ${\ell }_0$-constrained nonconvex M-estimation problem. We provide a high-dimensional analysis for the global solution of the ${\ell }_0$-constrained SDR and prove its asymptotic properties, including both variable selection consistency and convergence rate. We demonstrate the power of this general theory by obtaining the explicit convergence rates for the ${\ell }_0$-constrained SDR under the high-dimensional setting when the dimension p diverges. Furthermore, we also provide a high-dimensional analysis of the global solution of the ${\ell }_0$-constrained KSDR and prove its asymptotic properties, such as the consistency in variable selection and estimation, without assuming the linearity condition.

Computationally, the ${\ell }_0$-constrained trace maximisation is NP-hard. There has been some work showing success in adjusting the weight of coefficients adaptively, using convex penalties such as adaptive-lasso (Zou 2006) and adaptive-ridge regression (Frommlet and Nuel 2016). The performance of adaptive-based approaches heavily relies on the accuracy of initial estimation. Seeking for a more robust ${\ell }_0$ penalisation approximation without adaptive procedure, we use the folded concave penalisation (Fan and Li 2001; Fan, Xue, and Zou 2014) as an alternative to the ${\ell }_1$ penalisation in Li (2007), Lin et al. (2019), and Lin et al. (2018). We propose the novel nonconvex alternating direction method of multipliers (ADMM) to solve the folded concave penalised SDR and the nonconvex linearised ADMM to solve the folded concave penalised KSDR. In particular, sparse KSDR has not only a nonconvex loss function but also a folded concave penalty and a positive definite constraint. We provide the important convergence guarantees of the proposed nonconvex ADMMs for solving two challenging nonconvex and nonsmooth problems. We also establish the explicit iteration complexity bound for the proposed nonconvex ADMM and its linearised variant to reach the stationary solution of sparse SDR and sparse KSDR, respectively. To the best of our knowledge, our work is the first in the literature to study the challenging nonconvex and nonsmooth optimisation problem for sparse KSDR.

The rest of this paper is organised as follows. We will revisit the methods for sparse SDR in Section 2 and then present the proposed methods for sparse KSDR in Section 3. The proposed nonconvex ADMM and its linearised variant are presented in Sectio 4, and statistical and computational guarantees are presented in Section 5. Section 6 evaluates the numerical performance in several different simulation settings. Section 7 applies sparse KDR in a real data analysis. Section 7 includes a few concluding remarks.

Before proceeding, we define some useful notations. Let $\Vert \mathit{\theta }{\Vert }_0=\mathrm{card}\{{\theta }_j\ne 0,j=1,\dots d\}$, and $|\mathcal{S}|$ is the cardinality of a set $\mathcal{S}$. For a matrix $\mathit{A}$, let $\Vert \mathit{A}{\Vert }_{\mathrm{\infty }}$ be the maximum norm, $\Vert \mathit{A}{\Vert }_1$ the matrix ${\ell }_1$ norm, $\Vert \mathit{A}{\Vert }_2$ the spectral norm, $\Vert \mathit{A}{\Vert }_F$ the Frobenius norm, $\Vert \mathit{A}{\Vert }_{\mathrm{tr}}$ the trace norm, defined as $\Vert \mathit{A}{\Vert }_{\mathrm{tr}}=\mathrm{Tr}\left[\right({\mathit{A}}^T\mathit{A}{)}^{1/2}]$. Let $\Vert \cdot {\Vert }_{{\mathcal{H}}_d}$ denote the norm in d-dimensional Hilbert space. To estimate the central subspace, it is equivalent to finding a low-rank matrix $\mathbf{\Theta }\in {\mathbb{R}}^{p\times d}$ of full column rank such that $y\perp \perp \mathit{x}|{\mathbf{\Theta }}^T\mathit{x}$. We may also denote the central subspace by ${\mathit{S}}_{y|\mathit{x}}=\mathrm{span}\left(\mathbf{\Theta }\right)$, as the column space of matrix $\mathbf{\Theta }$.

2. Preliminaries

This section first revisits the sparse SDR methods in Li (2007) and then follows d'Aspremont et al. (2005) to reformulate sparse SDR via the trace maximisation problem. Note that SDR aims to estimate the linear combination $\mathbf{\Theta }$ such that $y\perp \perp \mathit{x}|{\mathbf{\Theta }}^T\mathit{x}$ given independently and identically distributed data $({\mathit{x}}_1,y_1),\dots ,({\mathit{x}}_n,y_n)$. Li (2007) showed that SDR methods based on inverse moments, such as SIR, DR, SAVE, and their variants, essentially solve the generalised eigenvector of a certain kernel matrix. When we assume the product of the inverse covariance matrix and kernel matrix is a normal matrix, the problem is also equivalent to solving an eigenvector of a certain kernel matrix. More specifically, at population level, the kernel matrix of SIR is ${\mathit{M}}_{\mathrm{sir}}^{\star }={\mathbf{\Sigma }}_x^{-1}\mathrm{cov}\left[E\right\{\mathit{z}|y\left\}\right]$; the kernel matrix of DR is: $\begin{array}{rl}{\mathit{M}}_{\mathrm{dr}}^{\star }& =2{\mathbf{\Sigma }}_x^{-1}\left(E\right[E^2\left(\mathit{z}{\mathit{z}}^T|y\right)]+E^2[E\left(\mathit{z}|y\right)E\left({\mathit{z}}^T|y\right)]\\ & +E\left[E\right({\mathit{z}}^T|y\left)E\right(\mathit{z}|y\left)\right]E\left[E\right(\mathit{z}|y\left)E\right({\mathit{z}}^T|y\left)\right]-{\mathit{I}}_p),\end{array}$where $\mathit{z}=\mathit{x}-E\left(\mathit{x}\right)$ and ${\mathbf{\Sigma }}_x$ the covariance matrix of all predictors $\mathit{x}$; the kernel matrix of DT-SIR (Lin et al. 2018) is ${\mathit{M}}_{\mathrm{DT}-\mathrm{sir}}^{\star }={\mathbf{\Sigma }}_{x_{\mathrm{DT}}}^{-1}\mathrm{cov}\left[E\right\{({\mathit{x}}_{\mathrm{DT}}-E{\mathit{x}}_{\mathrm{DT}})|y\left\}\right],$ where ${\mathit{x}}_{\mathrm{DT}}$ is the subset of random variables in $\mathit{x}$ after the diagonal thresholding.

Given the finite sample size n, we have an empirical version of such kernel matrices $\hat{\mathit{M}}$, by replacing population covariance matrix ${\mathbf{\Sigma }}_x$ with sample covariance matrix ${\hat{\mathbf{\Sigma }}}_x$, and all means $E(\cdot )$ with empirical mean $E_n(\cdot )$. For example, ${\hat{\mathit{M}}}_{\mathrm{sir}}={\hat{\mathbf{\Sigma }}}_x^{-1}\mathrm{cov}\left[E_n\right\{\mathit{z}|y\left\}\right]$. Li (2007) shows these SDR methods are equivalent to solving the leading eigenvector of the empirical version of such kernel matrices.

After estimating the kernel matrix $\hat{\mathit{M}}$, we study the sparse estimation of the leading SDR direction, i.e. the largest eigenvector of $\hat{\mathit{M}}$. We use the well-known fact: solving the largest eigenvector ${\mathit{\theta }}_1$ is equivalent to trace optimisation (Wold, Esbensen, and Geladi 1987) as follows: $\underset{\mathit{\theta }}{max}\mathrm{tr}\left({\mathit{\theta }}^T\hat{\mathit{M}}\mathit{\theta }\right)\mathrm{s}.\mathrm{t}.\Vert \mathit{\theta }{\Vert }_2^2\le 1.$To impose the sparsity, we introduce ${\ell }_0$ constraint in the following trace maximisation to penalise the cardinality of non-zero elements in vector $\mathit{\theta }$: ${\hat{\mathit{\theta }}}_1=\underset{\mathit{\theta }}{\arg max}\mathrm{tr}\left({\mathit{\theta }}^T\hat{\mathit{M}}\mathit{\theta }\right)\mathrm{s}.\mathrm{t}.\Vert \mathit{\theta }{\Vert }_2^2\le 1,\Vert \mathit{\theta }{\Vert }_0\le t$Note that $\mathrm{span}\left({\mathit{\theta }}^{T}x\right)=\mathrm{span}\left(\mathit{\theta }{\mathit{\theta }}^{T}x\right)$. With $\mathbf{\Theta }=\mathit{\theta }{\mathit{\theta }}^T$, then it is equivalent to solving the positive semidefinite estimate ${\hat{\mathbf{\Theta }}}_1$ from $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;{\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\le t^2;\mathbf{\Theta }⪰0;\mathrm{rank}\left(\mathbf{\Theta }\right)=1$where $\Vert \mathbf{\Theta }{\Vert }_0$ is defined as the element-wise ${\ell }_0$ norm. When $\mathbf{\Theta }⪰0$ and $\mathrm{rank}\left(\mathbf{\Theta }\right)=1$, $\mathrm{tr}\left(\mathbf{\Theta }\right)=1$ implies $\Vert \mathit{\theta }{\Vert }_2^2=1$, and ${\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\le t^2$ is equivalent to $\Vert \mathit{\theta }{\Vert }_0\le t$. However, the rank-one constraint that $\mathrm{rank}\left(\mathbf{\Theta }\right)=1$ is hard to solve efficiently in practice. Therefore, we consider a relaxation of removing this constraint. Instead, we solve the dominant eigenvector of ${\hat{\mathbf{\Theta }}}_1$ as ${\hat{\mathit{\theta }}}_1$, which is the SDR direction of our interest. Generally, this relaxation tends to have a solution very close to rank one matrix (d'Aspremont et al. 2005).

After removing the rank-one constraint, it is equivalent to solve ${\hat{\mathbf{\Theta }}}_1$ from the following ${\ell }_0$ penalised estimation problem for some sufficiently large γ: (1) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\gamma {\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(1) where given ${\hat{\mathbf{\Theta }}}_1$, we follow d'Aspremont et al. (2005) to truncate it to obtain the dominant (sparse) eigenvector of ${\hat{\mathit{\theta }}}_1$, and it will be the leading sparse SDR direction of our interest.

In the sequel, given ${\hat{\mathit{\theta }}}_1,\dots ,{\hat{\mathit{\theta }}}_k$, we solve the following sparse SDR direction ${\hat{\mathit{\theta }}}_{k+1}$ in a sequential fashion. We denote ${\hat{\mathit{M}}}_1=\hat{\mathit{M}}$, and for $k\ge 1$, let ${\hat{\mathit{M}}}_{k+1}$ be the projection of ${\hat{\mathit{M}}}_k$ on the space orthogonal to previous SDR directions. To remove the effect of previous eigenvectors, ${\hat{\mathit{\theta }}}_{k+1}$ is estimated as the largest (sparse) eigenvector of matrix ${\hat{\mathit{M}}}_{k+1}$, which is a Hotelling's deflation from previous vectors: ${\hat{\mathit{M}}}_{k+1}={\hat{\mathit{M}}}_k-\left({\hat{\mathit{\theta }}}_k^T{\hat{\mathit{M}}}_k{\hat{\mathit{\theta }}}_k\right){\hat{\mathit{\theta }}}_k{\hat{\mathit{\theta }}}_k^T.$Mackey (2009) studied the property of Hotelling's deflation in this type of sparse eigenvector estimation problem. Especially when we estimate the sparse eigenvector on the true support, this deflation helps eliminate the influence of a given eigenvector by setting the associated eigenvalue to zero and achieve ${\hat{\mathit{\theta }}}_1\perp \perp {\hat{\mathit{\theta }}}_2\perp \perp \dots \perp \perp {\hat{\mathit{\theta }}}_k$. For the sparse estimation of ${\hat{\mathit{\theta }}}_{k+1}$, we plug in ${\hat{\mathit{M}}}_{k+1}$ to (1(1) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\gamma {\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(1) ), then solve ${\mathit{\theta }}_{k+1}$ from: $\underset{\mathit{\theta }\in {\mathbb{R}}^p}{max}\mathrm{tr}\left({\mathit{\theta }}^T{\hat{\mathit{M}}}_{k+1}\mathit{\theta }\right)\mathrm{s}.\mathrm{t}.\Vert \mathit{\theta }{\Vert }_2^2\le 1,\Vert \mathit{\theta }{\Vert }_0\le t,$which is further equivalent to the following constrained problem following the same rank relaxation: (2) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left({\hat{\mathit{M}}}_{k+1}\mathbf{\Theta }\right)-\gamma {\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0.$(2) (2(2) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left({\hat{\mathit{M}}}_{k+1}\mathbf{\Theta }\right)-\gamma {\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0.$(2) ) is still a nonconvex problem since ${\ell }_0$ constraint is included. To solve the computational challenge, one existing relaxation is using folded concave penalty to approximate ${\ell }_0$ constraint. Let $P_{\lambda }\left(|\mathbf{\Theta }|\right)$ be a folded concave penalty function such as the SCAD (Fan and Li 2001) and MCP (Zhang 2010), where $|\mathbf{\Theta }|$ is the matrix such that taking the absolute value of matrix $\mathbf{\Theta }$ element-wisely.

Specifically, in this paper, we consider the element-wise SCAD penalty $P_{\delta }\left(t\right)$ as: $P_{\delta }\left(t\right)=\{\begin{array}{ll}{\displaystyle \frac{2\mathrm{at}}{(a+1)\delta }}& {\displaystyle |t|\le \frac{\delta }{a},}\\ {\displaystyle 1-\frac{(\delta -|t|{)}^2}{(1-\frac{1}{a^2}){\delta }^2}}& {\displaystyle \frac{\delta }{a}\le |t|\le \delta ,}\\ 1& |t|\ge \delta .\end{array}$where a is a constant, usually chosen as 3.7 in the SCAD penalty. In Section 5, we discuss how such penalisation can well approximate ${\ell }_0$ penalisation. In this paper, We use SCAD penalty in this paper for simplicity. It is worth mentioning that we can also apply other nonconvex penalties and approximations here, such as MCP (Zhang 2010) and coefficient thresholding (Liu, Zhang, Xue, Song, and Kang 2024) to achieve the same theoretical property and similar numerical performance.

With this approximation, we solve ${\hat{\mathbf{\Theta }}}_1$ from the following folded concave penalised estimation of sparse SDR: (3) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(3) and then for $k=1,2,\dots ,d-1$, sequentially solve ${\hat{\mathbf{\Theta }}}_{k+1}$, given ${\hat{\mathit{M}}}_{k+1}$, from $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left({\hat{\mathit{M}}}_{k+1}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0.$Without loss of generality, this idea of estimating leading sparse eigenvalues can also be applied to the sparse Principal Component Analysis (PCA) problem.

Theorem 4 will provide the computational guarantee for the above approximation to ${\ell }_0$ penalisation. Section 6 will also present the numerical performance.

3. Methodology

The methodology and theory of the sparse SDR depend on the linearity condition and constant variance condition, which usually hold when predictors follow elliptical distributions, but elliptical conditions can be violated in many cases. In this section, following the framework built by Fukumizu et al. (2009) and Li, Chun, and Zhao (2012), we extend the framework to achieve sparse KSDR, which does not require linearity and constant variance conditions. The key concept is that we construct a conditional covariance operator in reproducing kernel Hilbert space (RKHS), which evaluates the amount of information y that cannot be explained by the function of linear combinations of $\mathit{x}$. Then, we derive the empirical estimator of such conditional covariance operator with finite samples and finally minimise the trace of such estimator with sparsity constraints.

First, we provide the necessary definitions in RKHS to construct the conditional covariance operator. Let $k_{\mathit{x}}:{\mathrm{\Omega }}_{\mathit{x}}\times {\mathrm{\Omega }}_{\mathit{x}}\to \mathbb{R}$ and $k_y:{\mathrm{\Omega }}_y\times {\mathrm{\Omega }}_y\to \mathbb{R}$ be positive definite kernel functions and ${\mathcal{H}}_{\mathit{x}}$ and ${\mathcal{H}}_y$ be the RKHS generated by $k_{\mathit{x}}$ and $k_y$, such that ${\mathcal{H}}_{\mathit{x}}$ satisfies reproducing property: $\begin{array}{rl}& \mathrm{\forall }\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}},k(\cdot ,\mathit{x})\in {\mathcal{H}}_{\mathit{x}};\\ & \mathrm{\forall }\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}},\mathrm{\forall f}\in {\mathcal{H}}_{\mathit{x}},〈f,k(\cdot ,\mathit{x}){〉}_{{\mathcal{H}}_{\mathit{x}}}=f\left(\mathit{x}\right).\end{array}$And ${\mathcal{H}}_y$ satisfies such reproducing property accordingly in space ${\mathcal{H}}_y$. Then define ${\mathbf{\Sigma }}_{\mathit{xx}}:{\mathcal{H}}_{\mathit{x}}\to {\mathcal{H}}_{\mathit{x}}$ be the variance operator of $\mathit{x}$, which satisfies the relation $〈f,{\mathbf{\Sigma }}_{\mathit{xx}}g{〉}_{{\mathcal{H}}_{\mathit{x}}}=\mathrm{cov}\left(f\right(\mathit{x}),g(\mathit{x}\left)\right)$ for all $f,g\in {\mathcal{H}}_{\mathit{x}}$, and so as ${\mathbf{\Sigma }}_{\mathrm{yy}}$. Define ${\mathbf{\Sigma }}_{y\mathit{x}}:{\mathcal{H}}_{\mathit{x}}\to {\mathcal{H}}_y$ as the covariance operator satisfying $〈g,{\mathbf{\Sigma }}_{y\mathit{x}}f{〉}_{{\mathcal{H}}_y}=\mathrm{cov}\left(g\right(y),f(x\left)\right)$ for all $g\in {\mathcal{H}}_y$ and $f\in {\mathcal{H}}_{\mathit{x}}$, and ${\mathbf{\Sigma }}_{\mathit{x}y}={\mathbf{\Sigma }}_{y\mathit{x}}^{\ast }$ be the adjoint operator of ${\mathbf{\Sigma }}_{\mathrm{yx}}$. Based on Baker (1973), we know there exists unique bounded operators ${\mathit{V}}_{y\mathit{x}}$ and ${\mathit{V}}_{\mathit{x}y}$ satisfy: ${\mathbf{\Sigma }}_{y\mathit{x}}={\mathbf{\Sigma }}_{\mathrm{yy}}^{1/2}{\mathit{V}}_{y\mathit{x}}{\mathbf{\Sigma }}_{\mathit{xx}}^{1/2};{\mathbf{\Sigma }}_{\mathit{x}y}={\mathbf{\Sigma }}_{\mathit{xx}}^{1/2}{\mathit{V}}_{\mathit{x}y}{\mathbf{\Sigma }}_{\mathrm{yy}}^{1/2}.$Then, the conditional covariance operator can be defined as: ${\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{x}}={\mathbf{\Sigma }}_{\mathrm{yy}}-{\mathbf{\Sigma }}_{\mathrm{yy}}^{1/2}{\mathit{V}}_{y\mathit{x}}{\mathit{V}}_{\mathit{x}y}{\mathbf{\Sigma }}_{\mathrm{yy}}^{1/2}.$For convenience, we sometime write it as: ${\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{x}}={\mathbf{\Sigma }}_{\mathrm{yy}}-{\mathbf{\Sigma }}_{y\mathit{x}}{\mathbf{\Sigma }}_{\mathit{xx}}^{†}{\mathbf{\Sigma }}_{\mathit{x}y}.$ where ${\mathbf{\Sigma }}_{\mathit{xx}}^{†}$ is the pseudo inverse of ${\mathbf{\Sigma }}_{\mathit{xx}}$. Specifically, ${\mathbf{\Sigma }}_{\mathit{xx}}^{†}=({\mathbf{\Sigma }}_{\mathit{xx}}+{ϵ}_{\mathit{s}}\mathit{I}{)}^{-1}$, as the Tychonoff regularised inverse matrix of ${\mathbf{\Sigma }}_{\mathit{xx}}^{†}$, given some sufficient small constant ${ϵ}_{\mathit{s}}$.

Also, define the sparse orthogonal space ${\mathbb{S}}_d^p\left(\mathbb{R}\right)=\{\mathbf{\Theta }\in M(p\times d;\mathbb{R}\left)\right)|{\mathbf{\Theta }}^T\mathbf{\Theta }={\mathbf{I}}_d;\Vert {\mathbf{\Theta }}_{\cdot i}{\Vert }_0\le k$ for i=1,…,d$\}$, representing the leading d sparse SDR directions. When d = 1, it represents the leading sparse SDR vector.

In what follows, we study the sparse estimation of $\mathbf{\Theta }\in {\mathbb{S}}_d^p\left(\mathbb{R}\right)$ such that $y\perp \perp \mathit{x}|{\mathbf{\Theta }}^T\mathit{x}$. To simplify notation, we introduce $\mathit{S}={\mathbf{\Theta }}^T\mathit{x}$, and $k_{\mathit{s}}$, ${\mathbf{\Sigma }}_{\mathit{ss}}$, ${\mathbf{\Sigma }}_{y\mathit{s}}$ and ${\mathbf{\Sigma }}_{\mathit{s}y}$ are defined in the same manner as $k_{\mathit{x}}$, ${\mathbf{\Sigma }}_{\mathit{xx}}$, ${\mathbf{\Sigma }}_{y\mathit{x}}$ and ${\mathbf{\Sigma }}_{\mathit{x}y}$. Then we have: ${\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{s}}={\mathbf{\Sigma }}_{\mathrm{yy}}-{\mathbf{\Sigma }}_{y\mathit{s}}{\mathbf{\Sigma }}_{\mathit{ss}}^{†}{\mathbf{\Sigma }}_{\mathit{s}y}$As shown in Theorem 4 of Fukumizu et al. (2009), if ${\mathcal{H}}_{\mathit{x}}$ and ${\mathcal{H}}_y$ are rich enough RKHS, $y\perp \perp \mathit{x}|{\mathbf{\Theta }}^T\mathit{x}⟺{\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{s}}={\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{x}},$ and also for any $\mathbf{\Theta }\in {\mathbb{S}}_d^p\left(\mathbb{R}\right)$, we have ${\mathbf{\Sigma }}_{\mathrm{yy}|{\mathbf{\Theta }}^T\mathit{x}}\ge {\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{x}}.$ This conclusion holds for our sparse orthogonal space as well, under the sparsity assumption of SDR vectors.

This result indicates that such a conditional covariance operator is a good estimator for sufficient dimension reduction at the population level. However, the population-level operator is not available in practice. At the sample level, we need to derive the empirical version of such an operator in coordinate representation and minimise it. To start with, we have the coordinate representation of the kernel matrix as: ${\mathit{K}}_{\mathit{x}}=\left(\begin{array}{ccc}{k}_{\mathit{x}}({\mathit{x}}_1,{\mathit{x}}_1)& \cdots & k_{\mathit{x}}({\mathit{x}}_1,{\mathit{x}}_n)\\ ⋮& \ddots & ⋮\\ k_{\mathit{x}}({\mathit{x}}_n,{\mathit{x}}_1)& \cdots & k_{\mathit{x}}({\mathit{x}}_n,{\mathit{x}}_n)\end{array}\right).$Let ${\mathit{G}}_{\mathit{x}}=\mathit{Q}{\mathit{K}}_{\mathit{x}}\mathit{Q}$ be the centralised gram matrix, where $\mathit{Q}={\mathit{I}}_n-{\mathbf{1}}_n{\mathbf{1}}_n^T/n$. Similarly, we define ${\mathit{K}}_y,{\mathit{K}}_{\mathit{s}},{\mathit{G}}_y,{\mathit{G}}_{\mathit{s}}$. Following Lemma 1 of Li et al. (2012), under the spanning systems: $\begin{array}{rl}{\hat{\mathcal{H}}}_X& =\mathrm{span}\left\{k_x\right(\cdot ,{\mathit{x}}_i)-{\mathbb{E}}_n{k}_x(\cdot ,\mathit{x}):i=1,\dots ,n\},\\ {\hat{\mathcal{H}}}_Y& =\mathrm{span}\left\{k_y\right(\cdot ,y_i)-{\mathbb{E}}_n{k}_y(\cdot ,y):i=1,\dots ,n\}.\end{array}$The variance and covariance operators can be written as: $\begin{array}{rl}& {\hat{\mathbf{\Sigma }}}_{\mathit{xx}}={\mathit{G}}_{\mathit{x}},{\hat{\mathbf{\Sigma }}}_{y\mathit{x}}={\mathit{G}}_{\mathit{x}},{\hat{\mathbf{\Sigma }}}_{\mathit{x}y}={\mathit{G}}_y,{\hat{\mathbf{\Sigma }}}_{\mathrm{yy}}={\mathit{G}}_y\\ & {\hat{\mathbf{\Sigma }}}_{\mathit{ss}}={\mathit{G}}_{\mathit{s}},{\hat{\mathbf{\Sigma }}}_{\mathit{s}y}={\mathit{G}}_y,{\hat{\mathbf{\Sigma }}}_{y\mathit{s}}={\mathit{G}}_{\mathit{s}}.\end{array}$Placing these coordinate representations into ${\mathbf{\Sigma }}_{\mathrm{yy}|\mathit{s}}$, we derive the coordinate representation of ${\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathit{s}}$ as: ${\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathit{s}}={\mathit{G}}_y-{\mathit{G}}_{\mathit{s}}{\mathit{G}}_{\mathit{s}}^{\text{‡}}{\mathit{G}}_y,$where ${\mathit{G}}_{\mathit{s}}^{\text{‡}}=({\mathit{G}}_{\mathit{s}}+{ϵ}_{\mathit{s}}\mathit{I}{)}^{-1}$ is the Tychonoff regularised inverse matrix of ${\mathit{G}}_{\mathit{s}}$ given some sufficient small constant ${ϵ}_{\mathit{s}}$. Then we minimise the trace of $\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathbf{\Theta }}^T\mathit{x}}\right)$ to estimate sparse SDR directions from $\mathbf{\Theta }\in {\mathbb{S}}_d^p\left(\mathbb{R}\right)$.

To deal with the orthogonality constraint in ${\mathbb{S}}_d^p\left(\mathbb{R}\right)$, we use a similar sequential procedure as in Section 2. To this end, we solve ${\hat{\mathit{\theta }}}_1$ from the following trace minimisation to minimise the sum of conditional variance: $\underset{\mathit{\theta }\in {\mathbb{SP}}_1^p\left(\mathbb{R}\right)}{min}\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathit{\theta }}^T\mathit{x}}\right)$where ${\mathbb{SP}}_1^p\left(\mathbb{R}\right)=\{\mathit{\theta }|\mathit{\theta }\in {\mathbb{R}}^p,\Vert \mathit{\theta }{\Vert }_2=1;\Vert \mathit{\theta }{\Vert }_0\le t\}.$ It is worth pointing out the equivalence between the above trace minimisation problem and the following regularised nonparametric regression problem: $\begin{array}{rl}& \underset{\mathbf{\Theta }\in {\mathbb{S}}_1^p\left(\mathbb{R}\right)}{min}\underset{f\in {\mathcal{H}}_{\mathcal{X}}^{\mathcal{S}}}{min}\frac{1}{n}\sum _{i=1}^n\sum _{a=1}^{\mathrm{\infty }}E|{\xi }_a\left(Y_i\right)-\frac{1}{n}\sum _{j=1}^n\left({\xi }_a\right(Y_j\left)\right)-\left(f\right(X_i)\\ & -\frac{1}{n}\sum _{j=1}^n\left[f\right(X_j\left)\right]\left){|}^2\right\}+\frac{{ϵ}_n}{n}\Vert f{\Vert }_{{\mathcal{H}}_{\mathcal{X}}^{\mathcal{S}}}^2\end{array}$This equivalence provides insight into sparse KSDR, and it is also useful to establish the theoretical properties in Section 5.

In sparse KSDR, we assume y is the function of ${\mathit{\theta }}_1^{\star T}\mathit{x},\dots ,{\mathit{\theta }}_d^{\star T}\mathit{x}$. In this case, ${\mathit{\theta }}_1^{\star },\dots ,{\mathit{\theta }}_d^{\star }$ are the true sparse SDR vectors of our interest. Previously, we have introduced how to estimate ${\mathit{\theta }}_1^{\star }$. Here, we introduce how to solve these non-leading SDR directions sequentially.

Given that we have already estimated orthogonal sufficient directions ${\hat{\mathit{\theta }}}_1,\dots ,{\hat{\mathit{\theta }}}_k$ for $k\ge 1$, we sequentially solve the following KSDR direction ${\mathit{\theta }}^{k+1}$ from the minimisation of $\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\hat{\mathit{\theta }}}_1^T\mathit{x},\dots ,{\hat{\mathit{\theta }}}_k^T\mathit{x},{\mathit{\theta }}^T\mathit{x}}\right)$. By orthogonality of all sufficient directions, $\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\hat{\mathit{\theta }}}_1^T\mathit{x},\dots ,{\hat{\mathit{\theta }}}_k^T\mathit{x},{\mathit{\theta }}^T\mathit{x}}\right)=\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|({\hat{\mathit{\theta }}}_1{\hat{\mathit{\theta }}}_1^T+\cdots +{\hat{\mathit{\theta }}}_k{\hat{\mathit{\theta }}}_k^T+\mathit{\theta }{\mathit{\theta }}^T)\mathit{x}}\right)$. Hence, we only need to solve ${\hat{\mathit{\theta }}}_{k+1}^{\mathrm{init}}=\underset{\mathit{\theta }\in {\mathbb{SP}}_1^p\left(\mathbb{R}\right)}{min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|({\hat{\mathit{\theta }}}_1{\hat{\mathit{\theta }}}_1^T+\cdots +{\hat{\mathit{\theta }}}_k{\hat{\mathit{\theta }}}_k^T+\mathit{\theta }{\mathit{\theta }}^T)\mathit{x}}\right),$Then we project ${\hat{\mathit{\theta }}}_{k+1}^{\mathrm{init}}$ to the space orthogonal to all previous KSDR vectors, denoted as the final estimator ${\hat{\mathit{\theta }}}_{k+1}$.

Here, we justify why this projection procedure is correct. Define ${\mathit{\theta }}_p={\mathcal{P}}_{{\hat{\mathit{\theta }}}_1,\dots .{\hat{\mathit{\theta }}}_k}^{\mathrm{\perp }}\mathit{\theta }$ as the projection of $\mathit{\theta }$ on the direction orthogonal to all previous k KSDR vectors, such that $\mathrm{span}({\hat{\mathit{\theta }}}_1,\dots ,{\hat{\mathit{\theta }}}_k,\mathit{\theta })=\mathrm{span}({\hat{\mathit{\theta }}}_1,\dots ,{\hat{\mathit{\theta }}}_k,{\mathit{\theta }}_p)$. By definition of the conditional covariance operator, we have: $\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\hat{\mathit{\theta }}}_1^T\mathit{x},\dots ,{\hat{\mathit{\theta }}}_k^T\mathit{x},{\mathit{\theta }}^T\mathit{x}}\right)=\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\hat{\mathit{\theta }}}_1^T\mathit{x},\dots ,{\hat{\mathit{\theta }}}_k^T\mathit{x},{\mathit{\theta }}_p^T\mathit{x}}\right).$Since the non-orthogonal direction to previous KSDR directions plays no effect in reducing the trace, we only need to consider the projection ${\mathit{\theta }}_p={\mathcal{P}}_{{\hat{\mathit{\theta }}}_1,\dots .{\hat{\mathit{\theta }}}_k}^{\mathrm{\perp }}\mathit{\theta }$ for the sparse estimation of the following KSDR direction ${\mathit{\theta }}^{k+1}$. By the projection procedure, we guarantee that ${\hat{\mathit{\theta }}}_1\perp \perp {\hat{\mathit{\theta }}}_2\perp \perp \dots \perp \perp {\hat{\mathit{\theta }}}_k$.

Finally, we discuss solving ${\ell }_0$ constraint and rank one constraint. Recall that $P_{\lambda }\left(|\mathbf{\Theta }|\right)$ is a folded concave penalty. In practice, we solve the following folded concave relaxation of sparse KSDR to approximate the ${\ell }_0$ constraint inside ${\mathbb{SP}}_1^p\left(\mathbb{R}\right)$, and solve ${\hat{\mathbf{\Theta }}}_1$ from $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\mathbf{\Theta }}^{T}x}\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;P_{\delta }\left(|\mathbf{\Theta }|\right)\le t^{\prime };\mathbf{\Theta }⪰0,$which can be rewritten in terms of a folded concave penalised estimation problem: $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\mathbf{\Theta }}^{T}x}\right)+\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0.$Given ${\hat{\mathbf{\Theta }}}_1$, we follow Section 2 to truncate it to obtain the first KSDR direction ${\hat{\mathit{\theta }}}_1$. Next, in the sequential fashion, given ${\hat{\mathit{\theta }}}_1,\dots ,{\hat{\mathit{\theta }}}_k$ for $k\ge 1$, we solve ${\hat{\mathbf{\Theta }}}_{k+1}$ from $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|({\hat{\mathit{\theta }}}_1{\hat{\mathit{\theta }}}_1^T+\cdots +{\hat{\mathit{\theta }}}_k{\hat{\mathit{\theta }}}_k^T+{\mathbf{\Theta }}^T)\mathit{x}}\right)+P_{\lambda }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$where we used the fact that $\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\hat{\mathit{\theta }}}_1^T\mathit{x},\dots ,{\hat{\mathit{\theta }}}_k^T\mathit{x},{\mathit{\theta }}^T\mathit{x}}\right)=\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\hat{\mathit{\theta }}}_1^T\mathit{x},\dots ,{\hat{\mathit{\theta }}}_k^T\mathit{x},{\mathit{\theta }}_p^T\mathit{x}}\right)$. Now, we compute the projection ${\mathbf{\Theta }}_{k+1}^{\mathrm{proj}}={\mathbb{P}}_{{\mathit{\theta }}_1\oplus \cdots \oplus {\mathit{\theta }}_k}^{\mathrm{\perp }}{\hat{\mathbf{\Theta }}}_{k+1}$ and solve the dominant (sparse) eigenvector as the desired sparse KSDR direction ${\hat{\mathit{\theta }}}_{k+1}$.

4. Nonconvex optimisation

Section 4 first proposes an efficient nonconvex ADMM to solve the nonconvex optimisation of sparse SDR and then presents an efficient nonconvex linearised ADMM to solve the even more challenging nonconvex optimisation of sparse KSDR. Both nonconvex ADMMs enjoy the important computational guarantees, which will be presented in Section 5.

4.1. Nonconvex ADMM for sparse SDR

In this section, we propose a novel nonconvex ADMM to solve the challenging nonconvex and nonsmooth optimisation in the sparse SDR. The proposed nonconvex ADMM enjoys the convergence guarantees to the ‘ϵ-stationary solution’ (see the definition in Theorem 5) with an explicit iteration complexity bound. The convergence guarantee of nonconvex algorithms can be defined as varying from algorithms. The ϵ-stationary solution introduced in Jiang, Lin, Ma, and Zhang (2019) is an applicable one to depict the convergence in the aspect of the Lagrangian function. More details about the ϵ-stationary solution and computational guarantees can be found in Theorem 5 of Section 5.1.

To begin with, we introduce a new variable $\mathbf{\Omega }$ and also a slacking variable $\mathbf{\Phi }$. By adding the equality constraint $\mathbf{\Theta }-\mathbf{\Omega }-\mathbf{\Phi }=0$, (3(3) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(3) ) is equivalent to: $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{min}-\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)+P_{\lambda }\left(|\mathbf{\Theta }|\right)+\frac{\mu }{2}\Vert \mathbf{\Phi }{\Vert }_F^2\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Omega }\right)=1;\mathbf{\Omega }⪰0,\mathbf{\Theta }-\mathbf{\Omega }-\mathbf{\Phi }=0.$The slacking variable $\mathbf{\Phi }$ plays an essential role in the convergence guarantee of our proposed nonconvex ADMM. Accordingly, the augmented Lagrange function is $\begin{array}{rl}{\mathcal{L}}_p(\mathbf{\Theta },\mathbf{\Omega },\mathbf{\Phi },\mathbf{\Lambda })& =-\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)+P_{\lambda }\left(|\mathbf{\Theta }|\right)+\frac{\mu }{2}\Vert \mathbf{\Phi }{\Vert }_F^2+{\mathcal{I}}_{\left\{\mathrm{tr}\right(\mathbf{\Omega })=1;\mathbf{\Omega }⪰0\}}\\ & +〈\mathbf{\Lambda },\mathbf{\Theta }-\mathbf{\Omega }-\mathbf{\Phi }〉+\frac{\rho }{2}\Vert \mathbf{\Theta }-\mathbf{\Omega }-\mathbf{\Phi }{\Vert }_F^2.\end{array}$where $\mathbf{\Lambda }$ is the Lagrange variable and $\rho >0$. Given the initial value $({\mathbf{\Omega }}^0,{\mathbf{\Phi }}^0,{\mathbf{\Lambda }}^0)$, we solve $({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}^{q+1},{\mathbf{\Phi }}^{q+1},{\mathbf{\Lambda }}^{q+1})$ iteratively for $q=0,1,2,\dots$ as follows: $\begin{array}{rl}{\mathbf{\Theta }}^{q+1}& =\underset{\mathbf{\Theta }}{\arg min}{\mathcal{L}}_p(\mathbf{\Theta },{\mathbf{\Omega }}^q,{\mathbf{\Phi }}^q,{\mathbf{\Lambda }}^q)+\frac{h}{2}\Vert {\mathbf{\Theta }}^q-\mathbf{\Theta }{\Vert }_F^2\\ {\mathbf{\Omega }}^{q+1}& =\underset{\mathbf{\Omega }}{\arg min}{\mathcal{L}}_p({\mathbf{\Theta }}^{q+1},\mathbf{\Omega },{\mathbf{\Phi }}^q,{\mathbf{\Lambda }}^q)\\ {\mathbf{\Phi }}^{q+1}& ={\mathbf{\Phi }}^q-\gamma \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathbf{\Phi }}^q}{\mathcal{L}}_p({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}^{q+1},\mathbf{\Phi },{\mathbf{\Lambda }}^q)\\ {\mathbf{\Lambda }}^{q+1}& ={\mathbf{\Lambda }}^q-\rho ({\mathbf{\Theta }}^{q+1}-{\mathbf{\Omega }}^{q+1}),\end{array}$where h>0 is chosen to guarantee the convexity of $\mathbf{\Theta }$ step and γ & ρ are chosen to guarantee the convergence. Their choices will be discussed in Theorem 5. Note that the $\mathbf{\Phi }$ step and the $\mathbf{\Lambda }$ step are straightforward. Fortunately, we also have the closed-form solutions in both $\mathbf{\Theta }$ step and $\mathbf{\Omega }$ step. In the $\mathbf{\Theta }$ step, we know that $\begin{array}{rl}& \underset{\mathbf{\Theta }}{min}{\mathcal{L}}_p(\mathbf{\Theta },{\mathbf{\Omega }}^q,{\mathbf{\Phi }}^q,{\mathbf{\Lambda }}^q)+\frac{h}{2}\Vert {\mathbf{\Theta }}^q-\mathbf{\Theta }{\Vert }_F^2\\ & =\underset{\mathbf{\Theta }}{min}\frac{\rho +h}{2}{\Vert \mathbf{\Theta }-\frac{1}{\rho +h}(\rho {\mathbf{\Omega }}^q+\rho {\mathbf{\Phi }}^q+h{\mathbf{\Theta }}^q-{\mathbf{\Lambda }}^q-\hat{\mathit{M}})\Vert }_F^2+P_{\lambda }\left(|\mathbf{\Theta }|\right),\end{array}$which can be analytically solved by using the univariate folded concave thresholding rule $S_{\lambda }(b,c)=\arg \underset{u}{min}\frac{1}{2}\Vert u-b{\Vert }^2+c{P}_{\lambda }\left(|u|\right).$ More specifically, we have ${\mathbf{\Theta }}^{q+1}=S_{\lambda }\left(\frac{1}{\rho +h}\right(\rho {\mathbf{\Omega }}^q+\rho {\mathbf{\Phi }}^q+h{\mathbf{\Theta }}^q-{\mathbf{\Lambda }}^q-\hat{\mathit{M}}),\frac{\rho +h}{2}).$In the $\mathbf{\Omega }$ step, it is easy to see that $\underset{\mathbf{\Omega }}{min}{\mathcal{L}}_p({\mathbf{\Theta }}^{q+1},\mathbf{\Omega },{\mathbf{\Phi }}^q,{\mathbf{\Lambda }}^q)=\underset{\mathbf{\Omega }⪰0:\mathrm{tr}\left(\mathbf{\Omega }\right)=1}{min}\Vert \mathbf{\Omega }+{\mathbf{\Phi }}^q-{\mathbf{\Theta }}^{q+1}-\frac{1}{\rho }{\mathbf{\Lambda }}^q{\Vert }_F^2,$which can be solved by using the projection onto the simplex in the Euclidean space (See Algorithm 1 of Ma 2013). To solve ${\mathbf{\Omega }}^{q+1}$, we take the singular value decomposition of ${\mathit{W}}^q={\mathbf{\Theta }}^{q+1}-{\mathbf{\Phi }}^q+\frac{1}{\rho }{\mathbf{\Lambda }}^q$ as ${\mathit{W}}^q=\mathit{U}\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_p){\mathit{U}}^T$. After analytically solving the projection of $\mathit{\sigma }$ onto the simplex in the Euclidean space. Denoting $\hat{\mathit{\xi }}=({\hat{\xi }}_1,\dots ,{\hat{\xi }}_p)$, it has closed-form solution: $\hat{\mathit{\xi }}=\arg \underset{\mathit{\xi }}{min}\frac{1}{2}\Vert \mathit{\xi }-\mathit{\sigma }{\Vert }_2^2\mathrm{s}.\mathrm{t}.\sum _{j=1}^p{\xi }_j=t,{\xi }_j\ge 0.$Then we obtain the closed-form solution ${\mathbf{\Omega }}^{q+1}=\mathit{U}\mathrm{diag}({\hat{\xi }}_1,\dots ,{\hat{\xi }}_p){\mathit{U}}^T$.

In the $\mathbf{\Phi }$ step, we solve the derivative as: $\frac{\mathrm{\partial }}{\mathrm{\partial }{\mathbf{\Phi }}^q}{\mathcal{L}}_p({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}^{q+1},\mathbf{\Phi },{\mathbf{\Lambda }}^q)=\mu \mathbf{\Phi }-\mathbf{\Lambda }+\rho (\mathbf{\Phi }+\mathbf{\Omega }-\mathbf{\Theta }).$Now, we summarise the proposed nonconvex ADMM in Algorithm 1.

4.2. Nonconvex linearised ADMM for sparse KSDR

In this section, we provide a detailed derivation of the nonconvex linearised ADMM algorithm for sparse KSDR. We present the algorithm for solving the leading KSDR vector, and the algorithm for the following vectors follows similarly. To begin with, we introduce new variables ${\mathbf{\Omega }}_1$ and ${\mathbf{\Omega }}_2$ and slacking variables ${\mathbf{\Phi }}_1$ and ${\mathbf{\Phi }}_2$. By adding two equality constraints that $\mathbf{\Theta }-{\mathbf{\Omega }}_1-{\mathbf{\Phi }}_1=0,\mathbf{\Theta }-{\mathbf{\Omega }}_2-{\mathbf{\Phi }}_2=0$, (2.7) is equivalent to: $\begin{array}{rl}& \underset{(\mathbf{\Theta },{\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2,{\mathbf{\Phi }}_1,{\mathbf{\Phi }}_2)}{min}\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathbf{\Theta }}^{T}x}\right)+P_{\lambda }\left(|{\mathbf{\Omega }}_1|\right)+\frac{\mu }{2}(\Vert {\mathbf{\Phi }}_1{\Vert }_F^2+\Vert {\mathbf{\Phi }}_2{\Vert }_F^2)\\ & s.t\mathbf{\Theta }-{\mathbf{\Omega }}_1-{\mathbf{\Phi }}_1=0,\mathbf{\Theta }-{\mathbf{\Omega }}_2-{\mathbf{\Phi }}_2=0,\mathrm{tr}\left({\mathbf{\Omega }}_2\right)=1,{\mathbf{\Omega }}_2⪰0.\end{array}$The slacking variables ${\mathbf{\Phi }}_1$ and ${\mathbf{\Phi }}_2$ are essential for the convergence guarantee of our proposed algorithm. Accordingly, the augmented Lagrange function is $\begin{array}{rl}{\mathcal{L}}_{\rho }(\mathbf{\Theta },{\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2,{\mathbf{\Lambda }}_1,{\mathbf{\Lambda }}_2)& =\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathbf{\Theta x}}\right)+{\mathcal{I}}_{\left\{\mathrm{tr}\right({\mathbf{\Omega }}_2)=1;{\mathbf{\Omega }}_2⪰0\}}+\frac{\mu }{2}(\Vert {\mathbf{\Phi }}_1{\Vert }_F^2+\Vert {\mathbf{\Phi }}_2{\Vert }_F^2)\\ & +P_{\lambda }\left(|{\mathbf{\Omega }}_1|\right)+〈{\mathbf{\Lambda }}_1,\mathbf{\Theta }-{\mathbf{\Omega }}_1-{\mathbf{\Phi }}_1〉+〈{\mathbf{\Lambda }}_2,\mathbf{\Theta }-{\mathbf{\Omega }}_2-{\mathbf{\Phi }}_2〉\\ & +\frac{\rho }{2}\Vert \mathbf{\Theta }-{\mathbf{\Omega }}_1-{\mathbf{\Phi }}_1{\Vert }_F^2+\frac{\rho }{2}\Vert \mathbf{\Theta }-{\mathbf{\Omega }}_2-{\mathbf{\Phi }}_2{\Vert }_F^2,\end{array}$where ${\mathbf{\Lambda }}_1$ and ${\mathbf{\Lambda }}_2$ are the Lagrange variable and $\rho >0$.

It is important to point out that $\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathbf{\Theta }}^{T}x}\right)$ introduces the additional computational challenges since the $\mathbf{\Theta }$ step does not lead to a closed-form solution. To this end, we introduce a linearised counterpart of this trace function as $\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathbf{\Theta }}^{\left(q\right)T}x}\right)+(\mathbf{\Theta }-{\mathbf{\Theta }}^q{)}^T\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{\Theta }}\mathrm{tr}({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathbf{\Theta }}^{\left(q\right)T}x}).$After plugging this linearised term, we denote the new Lagrange function as ${\mathcal{L}}_{\rho }^{\ast }(\mathbf{\Theta },{\mathbf{\Theta }}^q,{\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2,{\mathbf{\Phi }}_1,{\mathbf{\Phi }}_2,{\mathbf{\Lambda }}_1,{\mathbf{\Lambda }}_2).$ Given the initial solution $({\mathbf{\Omega }}_1^0,{\mathbf{\Omega }}_2^0,{\mathbf{\Phi }}_1^0,{\mathbf{\Phi }}_2^0,{\mathbf{\Lambda }}_1^0,{\mathbf{\Lambda }}_2^0)$, we solve $({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}_1^{q+1},{\mathbf{\Omega }}_2^{q+1},{\mathbf{\Phi }}_1^{q+1},{\mathbf{\Phi }}_2^{q+1},{\mathbf{\Lambda }}_1^{q+1},{\mathbf{\Lambda }}_2^{q+1})$ iteratively for $q=0,1,2,\dots$ as follows: $\begin{array}{rl}{\mathbf{\Theta }}^{q+1}& =\underset{\mathbf{\Theta }}{\arg min}{\mathcal{L}}_{\rho }^{\ast }(\mathbf{\Theta },{\mathbf{\Omega }}_1^q,{\mathbf{\Omega }}_2^q,{\mathbf{\Phi }}_1^q,{\mathbf{\Phi }}_2^q,{\mathbf{\Lambda }}_1^q,{\mathbf{\Lambda }}_2^q)+\frac{h}{2}\Vert \mathbf{\Theta }-{\mathbf{\Theta }}^q{\Vert }_F^2\\ {\mathbf{\Omega }}_1^{q+1}& =\underset{{\mathbf{\Omega }}_1}{\arg min}{\mathcal{L}}_{\rho }({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2^q,{\mathbf{\Phi }}_1^q,{\mathbf{\Phi }}_2^q,{\mathbf{\Lambda }}_1^q,{\mathbf{\Lambda }}_2^q)\\ {\mathbf{\Omega }}_2^{q+1}& =\underset{{\mathbf{\Omega }}_2}{\arg min}{\mathcal{L}}_{\rho }({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}_1^{q+1},{\mathbf{\Omega }}_2,{\mathbf{\Phi }}_1^q,{\mathbf{\Phi }}_2^q,{\mathbf{\Lambda }}_1^q,{\mathbf{\Lambda }}_2^q)\\ {\mathbf{\Phi }}_1^{q+1}& ={\mathbf{\Phi }}_1^q-\gamma \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathbf{\Phi }}_1^q}{\mathcal{L}}_{\rho }({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}_1^{q+1},{\mathbf{\Omega }}_2^{q+1},{\mathbf{\Phi }}_1^q,{\mathbf{\Phi }}_2^q,{\mathbf{\Lambda }}_1^q,{\mathbf{\Lambda }}_2^q)\\ {\mathbf{\Phi }}_2^{q+1}& ={\mathbf{\Phi }}_2^q-\gamma \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathbf{\Phi }}_2^q}{\mathcal{L}}_{\rho }({\mathbf{\Theta }}^{q+1},{\mathbf{\Omega }}_1^{q+1},{\mathbf{\Omega }}_2^{q+1},{\mathbf{\Phi }}_1^q,{\mathbf{\Phi }}_2^q,{\mathbf{\Lambda }}_1^q,{\mathbf{\Lambda }}_2^q)\\ {\mathbf{\Lambda }}_1^{q+1}& ={\mathbf{\Lambda }}_1^q-\rho ({\mathbf{\Theta }}^{q+1}-{\mathbf{\Omega }}_1^{q+1})\\ {\mathbf{\Lambda }}_2^{q+1}& ={\mathbf{\Lambda }}_2^q-\rho ({\mathbf{\Theta }}^{q+1}-{\mathbf{\Omega }}_2^{q+1}),\end{array}$where h>0 is chosen to guarantee the convexity of $\mathbf{\Theta }$ step and $\gamma ,\rho$ are chosen to guarantee the algorithmic convergence (see Theorem 10). In what follows, we present the closed-form solution for each subproblem.

For $\mathbf{\Theta }$ step, with the linearisation, ${\mathbf{\Theta }}^{q+1}$ can be solved by taking derivatives as: $\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{\Theta }}{\mathcal{L}}_{\rho }^{\ast }(\mathbf{\Theta },{\mathbf{\Theta }}^q,{\mathbf{\Omega }}_1^q,{\mathbf{\Omega }}_2^q,{\mathbf{\Phi }}_1^q,{\mathbf{\Phi }}_2^q,{\mathbf{\Lambda }}_1^q,{\mathbf{\Lambda }}_2^q)+h(\mathbf{\Theta }-{\mathbf{\Theta }}^q)=0.$Defining $g\left(\mathbf{\Theta }\right)={\mathit{G}}_s({\mathit{G}}_s+{ϵ}_{n}I{)}^{-1}$, then $\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathbf{\Theta x}}\right)=\mathrm{tr}\left(g\right(\mathbf{\Theta }\left){\mathit{G}}_y\right)$, then the derivatives can be solved by using the chain rule of matrix derivatives: $\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{\Theta }}\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathbf{\Theta x}}\right)=\mathrm{vec}{\left(\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{\Theta }}g\right(\mathbf{\Theta }\left)\right)}^T\mathrm{vec}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }g\left(\mathbf{\Theta }\right)}\mathrm{tr}\right({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathbf{\Theta x}}\left)\right),$which has a complex closed form. For space consideration, the detailed calculation of closed-form derivatives and its numerical approximation is deferred to the supplement. Based on $\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{\Theta }}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|\mathbf{\Theta x}}\right)$, we can obtain the close-form derivative. By solving the linearisation, we obtain the minimiser: $\begin{array}{rl}& q_1(\mathbf{\Theta },{\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2,{\mathbf{\Phi }}_1,{\mathbf{\Phi }}_2,{\mathbf{\Lambda }}_1,{\mathbf{\Lambda }}_2)\\ & =\frac{\rho ({\mathbf{\Omega }}_1+{\mathbf{\Phi }}_1+{\mathbf{\Omega }}_2+{\mathbf{\Phi }}_2)-{\mathbf{\Lambda }}_1-{\mathbf{\Lambda }}_2-\frac{\mathrm{\partial }}{\mathrm{\partial }g\left(\mathbf{\Theta }\right)}\mathrm{tr}\left({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|\mathbf{\Theta x}}\right)+h{\mathbf{\Theta }}^q}{2\rho +h}.\end{array}$For ${\mathbf{\Omega }}_1$ step, it is equivalent to solve the subproblem $\begin{array}{rl}& \underset{{\mathbf{\Omega }}_1}{min}{P}_{\lambda }\left(\mathbf{\Theta }\right)+〈{\mathbf{\Lambda }}_1,\mathbf{\Theta }-{\mathbf{\Omega }}_1-{\mathbf{\Phi }}_1〉+\frac{\rho }{2}\Vert \mathbf{\Theta }-{\mathbf{\Omega }}_1-{\mathbf{\Phi }}_1{\Vert }_F^2\\ & =\underset{{\mathbf{\Omega }}_1}{min}\Vert {\mathbf{\Omega }}_1-(\mathbf{\Theta }-{\mathbf{\Phi }}_1+\frac{{\mathbf{\Lambda }}_1}{\rho }){\Vert }_F^2+P_{\lambda }\left(|{\mathbf{\Omega }}_1|\right).\end{array}$This can be efficiently solved by using the univariate penalised least squares solution: ${\mathbf{\Omega }}_1^{q+1}=S_{\lambda }(\mathbf{\Theta }-{\mathbf{\Phi }}_1+\frac{{\mathbf{\Lambda }}_1}{\rho },2),$where $S_{\lambda }$ has been defined in Section 4.1. For ${\mathbf{\Omega }}_2$ step, we are solving the optimisation problem: ${\mathbf{\Omega }}_2=\underset{\mathrm{tr}\left({\mathbf{\Omega }}_2\right)=1}{\arg min}-{\mathbf{\Lambda }}_2{\mathbf{\Omega }}_2+\frac{\rho }{2}\Vert \mathbf{\Theta }-{\mathbf{\Omega }}_2-{\mathbf{\Phi }}_2{\Vert }_F^2=\underset{\mathrm{tr}\left({\mathbf{\Omega }}_2\right)=1}{\arg min}\frac{\rho }{2}\Vert {\mathbf{\Omega }}_2+{\mathbf{\Phi }}_2-\mathbf{\Theta }-{\mathbf{\Lambda }}_2/\rho {\Vert }_F^2.$It can be solved by projection onto the simplex in the Euclidean space. We have presented exactly the same algorithm in Section 4.1. Thus, we omit it here. For ${\mathbf{\Phi }}_1$ and ${\mathbf{\Phi }}_2$ steps, we can solve them straightforwardly: $\begin{array}{rl}{\mathbf{\Phi }}_1^{q+1}& ={\mathbf{\Phi }}_1^q+\gamma (\mu {\mathbf{\Phi }}_1^q+{\mathbf{\Lambda }}_1+\rho ({\mathbf{\Phi }}_1^q+{\mathbf{\Omega }}_1^{q+1}-{\mathbf{\Theta }}^{q+1}\left)\right)\\ {\mathbf{\Phi }}_2^{q+1}& ={\mathbf{\Phi }}_2^q+\gamma (\mu {\mathbf{\Phi }}_2^q+{\mathbf{\Lambda }}_2+\rho ({\mathbf{\Phi }}_2^q+{\mathbf{\Omega }}_2^{q+1}-{\mathbf{\Theta }}^{q+1}\left)\right).\end{array}$Now we can summarise the proposed nonconvex linearised ADMM in Algorithm 2:

5. Theoretical properties

In this section, we establish a theoretical guarantee for both sparse SDR and sparse KSDR estimators. For both SDR and KSDR, first we study the statistical guarantee of ${\ell }_0$-constrained estimator, showing their consistency, allowing p diverging as $n\to \mathrm{\infty }$. Then, we study how our folded concave penalised estimator can approximate ${\ell }_0$-constrained estimator and prove the computation guarantee of our proposed nonconvex ADMM and its linearised version.

5.1. Asymptotic properties of sparse SDR

We study the statistical and computational properties of ${\ell }_0$-constrained SDR, such as the consistency and convergence rate in this subsection and the computational guarantees of folded concave penalised SDR in this subsection. The ${\ell }_0$-constrained SDR ${\hat{\mathit{\theta }}}_1=\underset{\mathit{\theta }}{\arg max}\mathrm{tr}\left({\mathit{\theta }}^T\hat{\mathit{M}}\mathit{\theta }\right)\mathrm{s}.\mathrm{t}.\Vert \mathit{\theta }{\Vert }_2^2\le 1,\Vert \mathit{\theta }{\Vert }_0\le t$achieves the desired theoretical properties under the high-dimensional setting, while folded concave penalised SDR is computationally feasible and it asymptotically converges to ${\ell }_0$-constrained SDR. First, we define the true sufficient dimension vectors. For the selected method, d is the smallest integer such that, the top d eigenvectors of ${\mathit{M}}^{\star }$ are ${\mathit{\theta }}_i^{\star }$ for $i=1,\dots ,d$, such that $y\text{╨}\mathit{x}|{\mathit{\theta }}_i^{\star }$ for $i=1,\dots ,d$. Then by definition, we know ${\mathit{\theta }}_1^{\star },\dots ,{\mathit{\theta }}_d^{\star }$ are true sufficient vectors. Now we introduce assumptions thm(A1)–thm(A3).

(A1)

$E({\xi }^T\mathit{x}|{\mathit{\theta }}_1^T\mathit{x},\dots ,{\mathit{\theta }}_d^T\mathit{x})$ is a linear function of ${\mathit{\theta }}_1^T\mathit{x}$, $\dots ,{\mathit{\theta }}_d^T\mathit{x}$ for any $\xi \in {\mathbb{R}}^p$.

(A2)

$\mathrm{cov}(\mathit{x}|{\mathit{\theta }}_1^T\mathit{x},\dots ,{\mathit{\theta }}_d^T\mathit{x})$ is degenerate.

(A3)

True sufficient dimension vectors ${\mathit{\theta }}_i^{\star }$ satisfy $\Vert {\mathit{\theta }}_i^{\star }{\Vert }_0\le t$ for $i=1,\dots ,d$.

Note that thm(A1), thm(A2), and thm(A3) are also called the linearity condition, the constant variance condition, and the sparsity condition, respectively. thm(A1) and thm(A2) hold when $\mathit{x}$ follows elliptical distributions. thm(A3) is a standard assumption in sparse SDR literature. See Li (1991), Li and Wang (2007), Li (2007) and Luo et al. (2022) for justifications of these conditions.

Let $\sin \mathrm{\angle }({\mathit{\theta }}_a,{\mathit{\theta }}_b)=\sqrt{1-({\mathit{\theta }}_a^T{\mathit{\theta }}_b{)}^2}$. We first study the convergence rate of the first estimated sparse SDR direction ${\hat{\mathit{\theta }}}_1$. Theorem 1 derives this estimation bound based on the general estimation bound that $\Vert \hat{\mathit{M}}-{\mathit{M}}^{\star }{\Vert }_2=O_p\left(\tau \right(n,p\left)\right)$.

Theorem 5.1

Suppose (A1)–(A3), $\Vert \hat{\mathit{M}}-{\mathit{M}}^{\star }{\Vert }_2=O_p\left(\tau \right(n,p\left)\right)$, ${\lambda }_1>{\lambda }_2$, $\Vert {\mathit{\theta }}_1^{\star }{\Vert }_0\le t$, and the spectral decomposition for ${\mathit{M}}^{\star }$ exists. Then we have: $\sin \mathrm{\angle }({\hat{\mathit{\theta }}}_1,{\mathit{\theta }}_1^{\star })=O_p\left(\tau \right(n,p\left)\right).$

In practice, $O_p\left(\tau \right(n,p\left)\right)$ has been derived for each specific SDR method. Corollary 1 discusses the explicit convergence rate, given the explicit convergence rate $O_p\left(\tau \right(n,p\left)\right)$ for different methods.

Corollary 5.2

1. Suppose (A1)–(A3), ${\lambda }_1>{\lambda }_2$, $\Vert {\mathit{\theta }}_1^{\star }{\Vert }_0\le t$, covariance matrix ${\mathbf{\Sigma }}_x={\mathit{I}}_{p\ast p}$, and $p=o\left(n^{1/2}\right)$ as $n\to \mathrm{\infty }$, then for sparse SIR and sparse DR, we have: $\sin \mathrm{\angle }({\hat{\mathit{\theta }}}_1,{\mathit{\theta }}_1^{\star })=O_p\left(p{n}^{-1/2}\right)$

2. Let $v\in [0,1]$ corresponds the stability of the central curve $E\left[\mathit{x}|y\right]$. When $p=o\left(n\right)$, with proper choice of number of slices $H=(\frac{n}{p}{)}^{\frac{1}{2v+2}}$ and (A3)–(A4) in Lin et al. (2018) are satisfied, for sparse SIR, we have: $\sin \mathrm{\angle }({\hat{\mathit{\theta }}}_1,{\mathit{\theta }}_1^{\star })\to 0\mathrm{as}n\to \mathrm{\infty }.$Further, if we have ${\mathbf{\Sigma }}_x={\mathit{I}}_{p\ast p}$, we have convergence rate: $\sin \mathrm{\angle }({\hat{\mathit{\theta }}}_1,{\mathit{\theta }}_1^{\star })=O_p\left({\left(\frac{p}{n}\right)}^{\frac{v}{2v+2}}\right)$

3. Let $s=|\mathcal{S}|\ll p$ where $\mathcal{S}=\left\{i|{\mathit{\theta }}_j\right(i)\ne 0\mathrm{for}\mathrm{some}j,1\le j\le d\}$. Under the high-dimensional setting that $p\gg n$, by assuming (A3)–(A6) and (S1) in Lin et al. (2018), and ${\mathbf{\Sigma }}_x={\mathit{I}}_{p\ast p}$, then for the sparse DT-SIR, with proper choice of $H=(\frac{n}{p}{)}^{\frac{1}{2v+2}}$, we have convergence rate: $\sin \mathrm{\angle }({\hat{\mathit{\theta }}}_1,{\mathit{\theta }}_1^{\star })=O_p\left({\left(\frac{s}{n}\right)}^{\frac{v}{2v+2}}\right).$

where the stability v of the central curve $E\left[\mathit{x}|y\right]$ is defined in definition 1 of Lin et al. (2019). We also provide the complete definition of it in the supplemental material.

It is worth mentioning that our framework is very flexible. By applying different SDR methods to estimate the the $\mathit{M}$ matrix, we can derive the convergence rate under different settings. Let $J^{\star }=\{j:{\theta }_{1j}^{\star }\ne 0\}$ is the support set of ${\mathit{\theta }}_1^{\star }$. We further prove the variable selection consistency in Theorem 2.

Theorem 5.3

Under the assumption of Theorem 1, further if ${min}_{j\in J^{\star }}|{\theta }_{1j}^{\star }|\ge C_0\tau (n,p)$ for some constant $C_0>0$, and $\Vert {\mathit{\theta }}_1^{\star }{\Vert }_0=t$, then $\mathrm{supp}\left({\hat{\mathit{\theta }}}_1\right)=\mathrm{supp}\left({\mathit{\theta }}_1^{\star }\right)$ with probability tending to 1 as $n\to \mathrm{\infty }$.

Next, recall that in Section 2, for $k\ge 2$, by the same optimisation procedure, we can further estimate ${\hat{\mathit{\theta }}}_k$ from ${\hat{\mathit{M}}}_k$. Next, we study the estimation bound and consistency for ${\hat{\mathit{\theta }}}_k$ ($k\ge 2$) in Theorem 3.

Theorem 5.4

Suppose that the conditions in Theorems 1 and 2 are satisfied, $\Vert {\mathit{\theta }}_k^{\star }{\Vert }_0\le t$, ${\lambda }_k>{\lambda }_{k+1}$ and ${\lambda }_k=\mathcal{O}\left(1\right)$ for $1\le k\le d-1$. Then for $1\le k\le d-1$, we have: $\begin{array}{rl}\Vert {\hat{\mathit{M}}}_{k+1}-{\mathit{M}}_{k+1}^{\star }{\Vert }_2& =O_p\left(\tau \right(n,p\left)\right)\\ \sin \mathrm{\angle }({\hat{\mathit{\theta }}}_{k+1},{\mathit{\theta }}_{k+1}^{\star })& =O_p\left(\tau \right(n,p\left)\right),\end{array}$where ${\hat{\mathit{M}}}_{k+1}$ and ${\mathit{M}}_{k+1}^{\star }$ are defined in Section 2.

Now, we provide the computational guarantees of our proposed nonconvex ADMM. First, we show that the folded concave penalised estimation well approximates the ${\ell }_0$ penalisation. To this end, we use the rescaled folded concave penalty function, as we defined in Section 2. It is easy to see ${lim}_{\delta \to 0}{P}_{\delta }\left(t\right)=\Vert t{\Vert }_0$ for any $t\in \mathbb{R}$. Thus, $P_{\delta }\left(t\right)$ converges to $|t{|}_0$ in the ${\ell }_0$ penalisation. For matrix $\mathbf{\Theta }$, we define the penalisation $P_{\delta }\left(|\mathbf{\Theta }|\right)=\sum _{i=1}^p\sum _{j=1}^{d}P\left({\mathrm{\Theta }}_{\mathrm{ij}}\right)$, which is the summation of univariate rescaled penalisation for all entries. Here, we reveal the connection between ${\ell }_0$ penalised programming (1(1) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\gamma {\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(1) ) and scaled folded-concave penalised programming (3(3) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(3) ) in the following theorem.

Theorem 5.5

Suppose we choose $\left\{{\delta }_u\right\}$ to be sequence such that ${lim}_{u\to \mathrm{\infty }}{\delta }_u=0^{+}$. For the programming (3(3) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(3) ), we choose penalisation function $P_{\delta }\left(t\right)$ as we defined. Let ${\mathbf{\Theta }}^u$ be the global minimiser of (3(3) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(3) ) when we choose $\delta ={\delta }_u$, and assume ${\mathbf{\Theta }}^{\mathrm{\infty }}$ is a limit point of $\left\{{\mathbf{\Theta }}^u\right\}$. Then ${\mathbf{\Theta }}^{\mathrm{\infty }}$ is also a global minimiser to (1(1) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\gamma {\mathbf{1}}^T\Vert \mathbf{\Theta }{\Vert }_0\mathbf{1}\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(1) ).

Since the nonconvex problem (3(3) $\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}}{max}\mathrm{tr}\left(\hat{\mathit{M}}\mathbf{\Theta }\right)-\lambda P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{s}.\mathrm{t}.\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,$(3) ) has widely varying curvature, it is nontrivial to prove the algorithmic convergence. By taking advantage of the structured nonconvexity, the next theorem shows that our proposed nonconvex ADMM converges to an ϵ-stationary solution with an explicit iteration complexity bound. In the following theorem, we use $\mathrm{vec}\left(\mathbf{\Theta }\right)$ to represent the vectorisation of a matrix from ${\mathbb{R}}^{p\times d}$ to a vector ${\mathbb{R}}^{\mathrm{pd}}$.

Theorem 5.6

Let L be the Lipschitz bound of loss function, $\rho >max(\frac{18\sqrt{3}+6}{13}L,6\mathrm{hL})$ and $\gamma \in \{\begin{array}{ll}{\displaystyle (\frac{\sqrt{13{\rho }^2-12\mathrm{\rho L}-72L^2}-2\rho }{9{\rho }^2-12\mathrm{\rho L}-72L^2},+\mathrm{\infty })}& {\displaystyle \rho \in (\frac{2+2\sqrt{19}}{3}L,+\mathrm{\infty })}\\ {\displaystyle (\frac{2\rho -\sqrt{13{\rho }^2-12\mathrm{\rho L}-72L^2}}{72L^2+12\mathrm{\rho L}-9{\rho }^2},}\\ \frac{2\rho +\sqrt{13{\rho }^2-12\mathrm{\rho L}-72L^2}}{72L^2+12\mathrm{\rho L}-9{\rho }^2})& {\displaystyle \rho \in (\frac{18\sqrt{3}+6}{13}L,\frac{2+2\sqrt{19}}{3}L]}\end{array}$For any $0<ϵ<min\{\frac{1}{L},\frac{1}{3h}\}$, with no more than $\mathcal{O}\left(\frac{1}{{ϵ}^4}\right)$ number of iterations, our proposed algorithm obtains an ϵ-stationary solution $(\hat{\mathbf{\Theta }},\hat{\mathbf{\Omega }},\hat{\mathbf{\Lambda }})$, as defined in Jiang et al. (2019), for any matrix $\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}$ with finite Frobenius norm satisfying that: $\begin{array}{rl}& \left(\mathrm{vec}\right(\mathbf{\Theta }-\hat{\mathbf{\Theta }}\left){)}^T\right[\mathrm{vec}(\hat{\mathit{M}}+\mathrm{\partial }{P}_{\lambda }(|\mathbf{\Theta }|)-\hat{\mathbf{\Lambda }}\mathbf{I})]\ge -ϵ\\ & \Vert \hat{\mathbf{\Theta }}-\hat{\mathbf{\Omega }}\Vert \le ϵ\\ & \mathrm{tr}\left(\hat{\mathbf{\Omega }}\right)=1,\hat{\mathbf{\Omega }}⪰0.\end{array}$

Remark: When $ϵ=0$, the constraint exactly implicates the KKT-condition. Thus, the solution will be a stationary solution. We will show numerical performance of such a near-stationary solution in the next section.

5.2. Asymptotic properties of sparse KSDR

In this section, we establish the consistency of kernel sparse dimension reduction under milder assumptions. Based on Fukumizu et al. (2009), we assume following assumptions thm(A4)–thm(A7).

(A4)

For any bounded continuous function g on $\mathcal{Y}$, the mapping $\mathbf{\Theta }\to E_x\left[E_{y|{\mathbf{\Theta }}^T\mathit{X}}\right[g\left(y\right)|{\mathbf{\Theta }}^T\mathit{X}{]}^2]$ is continuous on $\mathbf{\Theta }\in {\mathbb{S}}_k^p\left(\mathbb{R}\right)$ for k=1,2,…,d.

(A5)

For any $k=1,2,\dots d$, let $\mathbf{\Theta }\in {\mathbb{S}}_k^p\left(\mathbb{R}\right)$, and $P_{\mathbf{\Theta }}$ be the probability distribution of the random variable ${\mathbf{\Theta }}^T\mathit{x}$ on $\mathcal{X}$. The Hilbert space ${\mathcal{H}}_{\mathcal{X}}^S+\mathbb{R}$ is dense in $L^2\left(P_{\mathbf{\Theta }}\right)$ for any $\mathbf{\Theta }\in {\mathbb{S}}_k^p\left(\mathbb{R}\right)$.

(A6)

There is a measurable function $\varphi :\mathcal{X}\to \mathbb{R}$ such that $E|\varphi \left(\mathit{x}\right){|}^2<\mathrm{\infty }$ and the Lipschitz condition: $\Vert k_s({\mathbf{\Theta }}^T\mathit{x},\cdot )-k_s({\tilde{\mathbf{\Theta }}}^T\mathit{x},\cdot ){\Vert }_{{\mathcal{H}}_d}\le \varphi \left(\mathit{x}\right)D(\mathbf{\Theta },\tilde{\mathbf{\Theta }})$ holds for all $\mathbf{\Theta },\tilde{\mathbf{\Theta }}\in {\mathbb{SP}}_k^p\left(\mathbb{R}\right)$ and $x\in \mathcal{X}$ for any k=1,2,…,d, and the kernel $k_s(x,\cdot )$ is defined in Section 3. D is a distance which is compatible with the topology of ${\mathbb{SP}}_k^p\left(\mathbb{R}\right)$. And we assume $\Vert {\mathbf{\Sigma }}_{\mathrm{yy}}{\Vert }_{\mathrm{tr}},\Vert {\mathbf{\Sigma }}_{\mathrm{yx}}{\Vert }_{\mathrm{HS}}\le \mathrm{\infty }$.

(A7)

The true sparse SDR directions ${\mathbf{\Theta }}^{\star }=({\mathit{\theta }}_1^{\star },\dots ,{\mathit{\theta }}_d^{\star })$ satisfies ${\mathbf{\Theta }}^{\star }\in {\mathbb{SP}}_d^p\left(\mathbb{R}\right)$.

Note that thm(A1) and thm(A2) are not required for the sparse KSDR, and we rewrite the sparsity assumption thm(A3) in our defined sparse orthogonal space as thm(A7). In thm(A4) to thm(A6), we introduce some mild regularity conditions in RKHS, which are well justified in Fukumizu et al. (2009).

Theorem 5.7

Under assumption (A4) to (A7), and we also suppose true variable set $\mathrm{supp}\left({\mathit{\theta }}_i^{\star }\right)$ are correlated with a finite number of variables; Suppose our choice of kernel function is continuous and bounded, and suppose the regularisation parameter ${ϵ}_n$ satisfies: ${ϵ}_n\to 0;{\left(\frac{n}{p^t}\right)}^{\frac{1}{t+2}}{ϵ}_n\to \mathrm{\infty };p\lesssim n^{\alpha }\mathrm{such}\mathrm{that}\alpha <1/t\mathrm{as}n\to \mathrm{\infty }.$Then the functions $\mathrm{tr}\left[{\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathit{\theta }}^T\mathit{x}}\right]$ and $\mathrm{tr}\left[{\mathbf{\Sigma }}_{\mathrm{yy}|{\mathit{\theta }}^T\mathit{x}}\right]$ are continuous on ${\mathbb{SP}}_1^p\left(\mathbb{R}\right)$, and it converges in probability that: $\underset{\mathit{\theta }\in {\mathbb{SP}}_1^p\left(\mathbb{R}\right)}{sup}|\mathrm{tr}\left[{\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\mathit{\theta }}^T\mathit{x}}^{\left(n\right)}\right]-\mathrm{tr}\left[{\mathbf{\Sigma }}_{\mathrm{yy}|{\mathit{\theta }}^T\mathit{x}}\right]|\to 0.$

Theorem 6 establishes the consistency of the first KDR vector. The growth rate assumption in the theorem says that the growth rate of p is negatively affected by the sparsity level t. When t is smaller, p can diverge faster as $n\to \mathrm{\infty }$. Because of this uniform convergence, we establish the variable selection consistency by following Theorem 7, allowing a diverging dimension of p.

Theorem 5.8

Under all conditions in Theorem 6, for $p\lesssim n^{\alpha }$, such that $\alpha <1/t$ as $n\to \mathrm{\infty }$, and $k=1,\dots ,d$, we have: $\mathrm{Pr}\left(\mathrm{supp}\right({\hat{\mathit{\theta }}}_k)=\mathrm{supp}({\mathit{\theta }}_k^{\star }\left)\right)\to 1.$

When variable selection consistency is achieved, we can further construct estimation consistency of $\hat{\mathit{\theta }}$:

Theorem 5.9

If all conditions in Theorem 7 hold, and we assume $p\lesssim n^{\alpha }$, such that $\alpha <1/t$ as $n\to \mathrm{\infty }$. Then we denote: $\hat{\mathit{\theta }}=\underset{\mathit{\theta }}{\arg min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\mathit{\theta }}^T\mathit{x}}\right)\mathrm{s}.\mathrm{t}.\Vert \mathit{\theta }{\Vert }_2=1,\Vert \mathit{\theta }{\Vert }_0\le t$And define the corresponding matrix ${\tilde{\mathbf{\Theta }}}_1={\hat{\mathit{\theta }}}_1{\hat{\mathit{\theta }}}_1^T$. When ${\tilde{\mathbf{\Theta }}}_1$ is nonempty, for any arbitrary open set U in ${\mathbb{S}}_1^p(\mathbb{R}$) with true direction ${\mathit{\theta }}^{\star }{\mathit{\theta }}^{\star T}\in U$, we have (4) $\underset{n\to \mathrm{\infty }}{lim}Pr({\tilde{\mathbf{\Theta }}}_1\in U)=1.$(4) For following directions, let ${\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_k$ correspond to the true first k KDR vectors, and ${\hat{\mathit{\theta }}}_{k+1}=\underset{\mathit{\theta }}{\arg min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|(\sum _{i=1}^k{\hat{\mathit{\theta }}}_i{\hat{\mathit{\theta }}}_i^T+\mathit{\theta }{\mathit{\theta }}^T{)}^T\mathit{x}}\right)\mathrm{s}.\mathrm{t}.\Vert \mathit{\theta }{\Vert }_2=1,\Vert \mathit{\theta }{\Vert }_0\le t.$Denote ${\tilde{\mathbf{\Theta }}}_{k+1}=\sum _{i=1}^k{\hat{\mathit{\theta }}}_i{\hat{\mathit{\theta }}}_i^T+{\hat{\mathit{\theta }}}_{k+1}{\hat{\mathit{\theta }}}_{k+1}^T$, and the true space ${\mathbf{\Theta }}_{k+1}=\sum _{i=1}^{k+1}({\mathit{\theta }}_i^{\star }{\mathit{\theta }}_i^{\star }{)}^T$. When ${\hat{\mathit{\theta }}}_i$ are all not nonempty for $i=1,\dots ,k$, for any arbitrary open set U in ${\mathbb{S}}_{k+1}^p\left(\mathbb{R}\right)$ with ${\mathbf{\Theta }}_{k+1}\in U$, we have (5) $\underset{n\to \mathrm{\infty }}{lim}Pr({\tilde{\mathbf{\Theta }}}_{k+1}\in U)=1.$(5)

Theorem 8 extends the KSDR consistency results from Fukumizu et al. (2009) to a diverging p case. Making use of the ${\ell }_0$-constrained parameter space, we have both variable selection and estimation consistency under certain conditions.

Following the statistical property, Theorem 9 shows the asymptotic equivalency of folded concave constraint and ${\ell }_0$ constraint. Theorem 10 shows that our nonconvex linearised ADMM guarantees the convergence to an ϵ-stationary solution.

Similar with Theorem 4, for Theorem 9, we consider $\begin{array}{r}\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times d}}{min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\mathbf{\Theta }}^{T}x}\right)+\gamma {\mathbf{1}}^T|\mathbf{\Theta }{|}_0\mathbf{1}\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0\\ \underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times d}}{min}\mathrm{tr}\left({\mathbf{\Sigma }}_{\mathrm{yy}|{\mathbf{\Theta }}^{T}x}\right)+\gamma P_{\delta }\left(|\mathbf{\Theta }|\right)\mathrm{tr}\left(\mathbf{\Theta }\right)=1;\mathbf{\Theta }⪰0,\end{array}$where $P_{\delta }\left(|\mathbf{\Theta }|\right)$ is defined in Section 2. Considering d = 1, we can also solve sufficient directions iteratively, as introduced in Section 2.

Theorem 5.10

Suppose we choose $\left\{{\delta }_u\right\}$ to be sequence such that ${lim}_{u\to \mathrm{\infty }}{\delta }_u=0^{+}$. For the programming (4.12), we choose penalisation function $P_{\delta }\left(t\right)$ as we defined. Let ${\mathbf{\Theta }}^u$ be the global minimiser of (4.12) when we choose $\delta ={\delta }^u$, and ${\mathbf{\Theta }}^{\mathrm{\infty }}$ is a limit point of $\left\{{\mathbf{\Theta }}^u\right\}$. Then ${\mathbf{\Theta }}^{\mathrm{\infty }}$ is also a global minimiser to (4.11).

Parallel to Theorem 5, we show the property for the ϵ-stationary solution of our proposed nonconvex ADMM in Theorem 10.

Theorem 5.11

Let L be the Lipschitz constant of the loss function. Let $\rho >max(\frac{18\sqrt{3}+6}{13}L,6\mathrm{hL})$and $\gamma \in \{\begin{array}{ll}{\displaystyle (\frac{\sqrt{13{\rho }^2-12\mathrm{\rho L}-72L^2}-2\rho }{9{\rho }^2-12\mathrm{\rho L}-72L^2},+\mathrm{\infty })}& {\displaystyle \rho \in (\frac{2+2\sqrt{19}}{3}L,+\mathrm{\infty })}\\ {\displaystyle (\frac{2\rho -\sqrt{13{\rho }^2-12\mathrm{\rho L}-72L^2}}{72L^2+12\mathrm{\rho L}-9{\rho }^2},}\\ \frac{2\rho +\sqrt{13{\rho }^2-12\mathrm{\rho L}-72L^2}}{72L{}^2+12\mathrm{\rho L}-9{\rho }^2})& {\displaystyle \rho \in (\frac{18\sqrt{3}+6}{13}L,\frac{2+2\sqrt{19}}{3}L]}\end{array}.$For any $0<ϵ<min\{\frac{1}{L},\frac{1}{3h}\}$, with no more than $\mathcal{O}\left(\frac{1}{{ϵ}^4}\right)$ number of iterations, our proposed algorithm obtains an ϵ-stationary solution $(\hat{\mathbf{\Theta }},{\hat{\mathbf{\Omega }}}_1,{\hat{\mathbf{\Omega }}}_2,{\hat{\mathbf{\Lambda }}}_1,{\hat{\mathbf{\Lambda }}}_2)$, as defined in Jiang et al. (2019), for any $\mathbf{\Theta }\in {\mathbb{R}}^{p\times p}$ with finite Frobenius norm,satisfying that $\begin{array}{rl}& \left(\mathrm{vec}\right(\mathbf{\Theta }-\hat{\mathbf{\Theta }}\left){)}^T\right[\mathrm{vec}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{\Theta }}\mathrm{tr}\right({\hat{\mathbf{\Sigma }}}_{\mathrm{yy}|{\hat{\mathbf{\Theta }}}^{T}x})-({\hat{\mathbf{\Lambda }}}_1+{\hat{\mathbf{\Lambda }}}_2\left)\mathbf{I}\right)]\ge -ϵ\\ & ({\mathbf{\Omega }}_1-{\hat{\mathbf{\Omega }}}_1{)}^T[P_{\lambda }^{\mathrm{\prime }}\left({\mathbf{\Omega }}_1\right)-({\hat{\mathbf{\Lambda }}}_1+{\hat{\mathbf{\Lambda }}}_2)\mathbf{I}]\ge -ϵ\\ & \Vert \hat{\mathbf{\Theta }}-{\hat{\mathbf{\Omega }}}_1\Vert \le ϵ\\ & \Vert \hat{\mathbf{\Theta }}-{\hat{\mathbf{\Omega }}}_2\Vert \le ϵ\\ & \mathrm{tr}\left({\hat{\mathbf{\Omega }}}_2\right)=1,{\hat{\mathbf{\Omega }}}_2⪰0.\end{array}$

6. Numerical properties

In this section, we evaluate the performance of the sparse SIR/DR/KSDR estimator in our semidefinite relaxation framework, the DT-SIR estimator proposed by Lin et al. (2019), the variable selection of penalised linear regression with SCAD penalisation and marginal correlation screening.

Here we compare the numerical performance of these methods in four data generation models as follows with sample size n = 100.

• Model 1: $y=\frac{0.5(x_1+x_2)}{0.5+\left(0.5\right(x_1-x_2)+1.5{)}^2}+(1+0.5(x_1-x_2){)}^2+{\sigma }^{2}e$, where $e\sim N(0,1)$, $\mathit{x}\sim N(0,{\mathit{I}}_{10})$, and $\sigma =0.2,0.4,0.8$ for Model 1a–1c.

• Model 2: $y=\frac{1}{2}(x_1^2+x_2^2)+{\sigma }^{2}e$, where $e\sim N(0,1)$, $\mathit{x}\sim N(0,{\mathit{I}}_{20})$, and $\sigma =0.1,0.2,0.3$ for Model 2a–2c.

• Model 3: $y=\frac{1}{2}(x_1+1{)}^2+{\sigma }^{2}e$, where $e\sim N(0,1)$, $\mathit{x}\sim N(0,{\mathit{I}}_{40})$, and $\sigma =0.2,0.4,0.8$ for Model 3a–3c.

• Model 4: $y={\sin }^2(\pi x_1+1)+{\sigma }^{2}e$, where $\sigma =0.2,0.4,0.8$ for Model 4a-4c, $e\sim N(0,1)$, and $\mathit{x}$ is uniformly distributed on $\{\mathit{x}\in {\mathbb{R}}^{10}|\mathrm{\exists }{x}_i>0.7,\mathrm{for}i=1,2,\dots 10\}$

To compare the numerical performance of these methods, we replicate our data generation procedure 100 times and apply all the methods to these 100 replications. We compute the estimated sparse SDR vector ${\hat{\mathit{\theta }}}_k$ for $k=1,2,\dots ,d$ for each method each time. Then, we evaluate the performance of different methods using multiple $R^2$, True Positive Rate (TPR) and True Negative Rate (TNR). The multiple $R^2$ is defined as $R^2=\frac{1}{d}\sum _{k=1}^d\Vert P_{\mathrm{span}\left({\mathbf{\Theta }}^{\star }\right)}{\hat{\mathit{\theta }}}_k{\Vert }_2^2.$ True Positive Rate(TPR) is defined as the fraction of variables correctly specified as non-zero in at least one out of SDR vectors, and True Negative Rate(TNR) is the fraction of variables correctly specified as zero in all SDR vectors.

The simulation results are summarised in the following Table 1–3. In the following Tables 1 and 2, we report the mean and standard deviation for these three indicators, with sparse SIR/DR/KSDR, DT-SIR, and SCAD penalised regression. The marginal correlation method doesn't have $R^2$, and we report its true positive/negative rate in Table 3. The tuning procedure is as follows. For the proposed methods, we choose the penalisation parameter λ via the oracle tuning procedure. Specifically, we construct two independent datasets with the same sample size n = 100 for training and tuning respectively. For each λ, we estimate the parameter $\mathbf{\Theta }$ on the training dataset and evaluate the corresponding loss function on the tuning dataset. We choose the best penalisation parameter to minimise the loss function on the tuning dataset. In practice, we may use the cross-validation to find the proper penalisation parameter.

Table 1. Summary of $R^2$, TP and TN for model 1 and 2.

	SpKSDR	DT-SIR	SpSIR	SpDR	SCAD	SpKSDR	DT-SIR	SpSIR	SpDR	SCAD
	Model 1a					Model 2a
$R^2$	0.833	0.843	0.980	0.995	0.656	0.975	0.224	0.215	0.965	0.172
	(0.303)	(0.153)	(0.098)	(0.050)	(0.063)	(0.112)	(0.250)	(0.278)	(0.146)	(0.355)
TPR	0.917	0.834	0.965	0.916	1	0.975	0.570	0.220	0.965	0.060
	(0.096)	(0.223)	(0.067)	(0.013)	(0)	(0.224)	(0.341)	(0.278)	(0)	(0.192)
TNR	0.975	0.423	0.976	0.969	0.660	0.988	0.429	0.805	0.990	0.905
	(0.069)	(0.255)	(0.066)	(0.091)	(0.191)	(0.038)	(0.270)	(0.070)	(0.042)	(0.103)
	Model 1b					Model 2b
$R^2$	0.917	0.834	0.965	0.916	0.654	0.950	0.235	0.180	0.929	0.183
	(0.203)	(0.158)	(0.151)	(0.278)	(0.064)	(0.154)	(0.274)	(0.261)	(0.199)	(0.363)
TPR	0.917	1	0.970	0.916	1	0.950	0.490	0.180	0.940	0.060
	(0.231)	(0)	(0.119)	(0.187)	(0)	(0.230)	(0.362)	(0.261)	(0.191)	(0.192)
TNR	0.975	0.423	0.976	0.969	0.626	0.981	0.464	0.795	0.949	0.895
	(0.069)	(0.255)	(0.066)	(0.091)	(0.229)	(0.205)	(0.270)	(0.065)	(0.181)	(0.136)
	Model 1c					Model 2c
$R^2$	0.900	0.757	0.950	0.720	0.649	0.900	0.233	0.185	0.920	0.196
	(0.203)	(0.158)	(0.151)	(0.278)	(0.073)	(0.230)	(0.266)	(0.263)	(0.197)	(0.371)
TPR	0.917	1	0.955	0.720	1	0.950	0.525	0.195	0.930	0.060
	(0.190)	(0)	(0.144)	(0.278)	(0)	(0.230)	(0.392)	(0.274)	(0.173)	(0.192)
TNR	0.967	0.375	0.978	0.930	0.610	0.975	0.454	0.799	0.981	0.893
	(0.080)	(0.228)	(0.060)	(0.070)	(0.236)	(0.052)	(0.281)	(0.069)	(0.048)	(0.138)

Table 2. Summary of $R^2$, TP and TN for model 3 and 4.

	SpKSDR	DT-SIR	SpSIR	SpDR	SCAD	SpKSDR	DT-SIR	SpSIR	SpDR	SCAD
	Model 3a					Model 4a
$R^2$	0.889	0.953	1	0.990	0.995	0.667	0.449	0.600	0.128	0.608
	(0.333)	(0.060)	(0)	(0.099)	(0.017)	(0.479)	(0.337)	(0.492)	(0.333)	(0.445)
TPR	0.889	1	1	1	1	0.667	0.900	0.600	0.130	0.700
	(0.333)	(0)	(0)	(0)	(0)	(0.479)	(0.302)	(0.462)	(0.338)	(0.461)
TNR	0.997	0.850	1	0.990	0.983	0.963	0.502	0.950	0.873	0.881
	(0.009)	(0.138)	(0)	(0.239)	(0.039)	(0.052)	(0.297)	(0.062)	(0.132)	(0.172)
	Model 3b					Model 4b
$R^2$	0.850	0.913	1	0.940	0.996	0.567	0.381	0.360	0.100	0.438
	(0.362)	(0.099)	(0)	(0.238)	(0.012)	(0.504)	(0.332)	(0.482)	(0.293)	(0.451)
TPR	0.850	1	1	1	1	0.567	0.840	0.360	0.110	0.550
	(0.365)	(0)	(0)	(0)	(0)	(0.504)	(0.368)	(0.482)	(0.314)	(0.500)
TNR	0.996	0.840	1	0.940	0.985	0.952	0.542	0.920	0.889	0.873
	(0.009)	(0.138)	(0)	(0.239)	(0.030)	(0.055)	(0.268)	(0.060)	(0.052)	(0.127)
	Model 3c					Model 4c
$R^2$	0.778	0.573	1	0.500	0.992	0.533	0.257	0.320	0.164	0.335
	(0.440)	(0.348)	(0)	(0.503)	(0.021)	(0.507)	(0.285)	(0.469)	(0.367)	(0.429)
TPR	0.778	0.900	1	0.500	1	0.533	0.730	0.320	0.210	0.450
	(0.440)	(0.302)	(0)	(0.503)	0	(0.507)	(0.447)	(0.469)	(0.409)	(0.500)
TNR	0.994	0.781	1	0.987	0.984	0.948	0.520	0.915	0.862	0.863
	(0.011)	(0.178)	(0)	(0.013)	(0.035)	(0.055)	(0.287)	(0.059)	(0.156)	(0.136)

Table 3. Summary of TP and TN for marginal correlation method.

	Model 1c	Model 2c	Model 3c	Model 4c
TPR	0.980	0.425	1	0.550
	(0.098)	(0.329)	(0)	(0.500)
TNR	0.873	0.913	1	0.950
	(0.012)	(0.018)	(0)	(0.056)

For implementation, Sparse SIR/DR/KSDR are implemented using our proposed nonconvex ADMM, DT-SIR is implemented using R package ‘LassoSIR’, and SCAD penalised regression is implemented using R package ‘ncvreg’. For the marginal correlation method, we compute the correlation between y and each variable $x_i$, then rank the absolute value of the correlation and select the first d variables, where d is the true dimension.

We also compare the performance of $R^2$ for kernel dimension reduction in Fukumizu et al. (2009) in Model 4(a, b, c) as a benchmark. The average $R^2$ are 0.642, 0.546 and 0.513 for Model 4(a), 4(b) and 4(c) separately. The results are comparable with SKDR results. SKDR performs slightly better given it eliminates some noise from irrelevant variables, thanks to the penalisation. Given kernel dimension reduction doesn't provide a sparse solution, and doesn't provide a meaningful TPR and TNR, we didn't include them in the main table.

We draw the following conclusions from the results: (1) sparse SIR and DT-SIR are compatible in most cases, while sparse SIR outperforms DT-SIR in $R^2$ in general. In model 1 and model 3, when the assumption for SIR holds, our sparse SIR is almost perfect. DT-SIR also performs well in these two settings but sometimes includes many irrelevant variables and leads to an unsatisfying TNR and $R^2$. (2) The sparse DR is also compatible with SIR in general. Although it may not outperform SIR in model 1 and model 3, when the function is symmetric over $\mathit{x}$ in model 2, SIR fails to recover the true space, while sparse DR is very stable in this case. (3) sparse KSDR performs best in model 4 when linearity conditions thm(A1) and thm(A2) don't hold, and both SIR and DR lose consistency. While sparse KSDR also performs fairly well compared to SIR and DR in models 1–3, when the linearity condition holds and no matter whether the function is symmetric over $\mathit{x}$. (4) The penalised regression fails in model 2 and model 4 when the relationship between $\mathit{x}$ and y cannot be approximated by linear regression. Also, regression cannot identify the linear combination $\mathit{\theta }$, which is also of interest. (5) Similarly, the marginal correlation method fails in model 2 when the correlation between $\mathit{x}$ and y is completely non-linear, and they have population correlation zero. What's more, using correlation-based methods also cannot identify the linear combination $\mathit{\theta }$ of our interest.

7. Application to microbiome data analysis

In this section, we apply our method to analyse the host transcriptomic data and 16S rRNA gene sequencing data from paired biopsies from IPAA patients with UC and familial adenomatous polyposis (Morgan et al. 2015). Morgan et al. (2015) extracted 5 principal components of microbiome 16s rRNA gene sequencing data, and analysed the relationship between these PCs and the PCs from host transcriptomic data. However, the result is hard to interpret since PCs are linear combinations of all host transcriptomics. Thus we apply our proposed sparse KSDR to explore the more interpretable sparse linear combinations which are relevant to these microbiome PCs.

First we select 20 representative genes of host transcriptomic data from two groups. The first group of genes are those most different from two locations (PPI and Porch) of a patient by Kolmogorov–Smirnov test. The second group of genes are important genes pointed out in Morgan et al. (2015), which are important to the pathogenesis of inflammatory bowel disease (IBD). And we treat 9 principal components of microbiome 16s rRNA gene sequencing as y in our model separately, and apply our kernel sparse diemsnion reduction estimation on the host transcriptomic data of PPI for each patient, with sample size n = 196 and p = 20. The selected genes and their coefficients are summarised in the following Table 4.

Table 4. Selected genes and corresponding coefficients using KSDR.

PCs	Gene 1	Gene 2	Gene 3	Gene 4
PC1	MMEL1(0.52)	IL1RN(−0.50)	TLR1(−0.30)	IL10(−0.27)
PC2	TNF (−0.42)	MUC6(0.37)	TLR1(−0.36)	NFKB1(−0.29)
PC3	MUC6(0.58)	IL1RN(0.50)	SOX14(0.45)	FUT2(0.24)
PC4	MUC6(0.74)	IL1RN7096 (0.36)	TLR1(0.26)	NFKB1(0.23)
PC5	TLR1(−0.49)	MUC6(0.41)	IFNG(−0.40)	NFKB1(0.39)

In the result, leading PCs (PC1-PC5) are mainly correlated with protein-related genes such as TLR1 and MMEL1 It reflects that these PCs may be mainly composed of antibiotic-signature microbiome. Several genes that have been widely studied that are important to the pathogenesis of IBD is also selected. For example, IL10 is selected in PC1, MUC6 is selected in PC2, 3, 4, and 5, and IL1RN is selected in PC1, 3, 4. Among them, IL10 is well studied that its effectiveness in signalling a subgroup of patients with IBD (Begue et al. 2011). MUC6 may have a role in epithelial wound healing after mucosal injury in IBD in addition to mucosal protection (Buisine et al. 2001). IL1RN has a well-established pathological role in IBD (Stokkers et al. 1998). Our result emphasises the effects of these genes may be actively related to microbiome RNA environment of IBD patients.

8. Conclusion

Motivated by the desired property of sparse SDR, we propose the nonconvex estimation scheme for sparse KSDR. In aspect of statistical property, we prove the asymptotic consistency for both estimators. In the aspect of optimisation, we show that they can both be solved by nonconvex ADMM algorithm efficiently with a provable convergence guarantee. We also show the practical usefulness of our proposed methods in various simulation settings and real data analysis.

Acknowledgments

The authors were grateful for the insightful comments and suggestions from the editor and reviewers.

Disclosure statement

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

Additional information

Funding

The authors were partially supported by the National Institutes of Health (NIH) grant 1R01GM152812 and National Science Foundation (NSF) grants DMS-1811552, DMS-1953189, CCF-2007823, and DMS-2210775.