# Robust Sufficient Dimension Reduction via $\alpha$-Distance Covariance

Hsin-Hsiung Huanga and Feng Yub and Teng Zhangc CONTACT Hsin-Hsiung Huang.
Email<EMAIL_ADDRESS>Teng Zhang. Email<EMAIL_ADDRESS>a
Department of Statistics and Data Science, University of Central Florida, USA;
b Department of Mathematics, University of Minnesota Twin Cities, Minnesota,
USA; c Department of Mathematics, University of Central Florida, USA

###### Abstract

We introduce a novel sufficient dimension-reduction (SDR) method which is
robust against outliers using $\alpha$-distance covariance (dCov) in
dimension-reduction problems. Under very mild conditions on the predictors,
the central subspace is effectively estimated and model-free without
estimating link function based on the projection on the Stiefel manifold. We
establish the convergence property of the proposed estimation under some
regularity conditions. We compare the performance of our method with existing
SDR methods by simulation and real data analysis and show that our algorithm
improves the computational efficiency and effectiveness.

###### keywords:

$\alpha$-distance covariance; central subspace; sufficient dimension reduction
(SDR); manifold learning; robust statistics

††articletype: Journal of Nonparametric Statistics

##  1 Introduction

In regression analysis, sufficient dimension reduction (SDR) provides a useful
statistical framework to analyze a high-dimensional dataset without losing any
information. It finds the fewest linear combinations of predictors that
capture a full regression relationship. Let $Y$ be an univariate response and
$X=(x_{1},\ldots,x_{p})^{T}$ be a $p\times 1$ predictor vector, SDR aims to
find a $p\times d$ matrix $\bm{\beta}$ such that

$Y\perp\\!\\!\\!\perp X\mid\bm{\beta}^{T}X$

which denotes the statistical independence.

Sufficient dimension reduction (SDR) based on the conditional distribution of
the response (KC. Li, 19911; Cook Weisberg, 1991; Xia ., 2002; Yin Li, 2011)
provides the reduced predictors without loss of regression information.
Recently, SDR methods using distance covariance (dCov) have been developed
(Sheng Yin, 2013, 2016), and such methods do not need a constant covariance
condition, or distribution assumptions on $X$, $X\mid Y$ or $Y\mid X$.
Therefore, it has broad applications for continuous and discrete variables
from various distributions. Several robust sufficient dimension reduction
methods have proposed for coefficient estimation such as the robust sufficient
dimension reduction using the ball covariance (J. Zhang Chen, 2019) and the
expected likelihood based method that minimizes the Kullback-Leiblier distance
(Yin Cook, 2005; N. Zhang Yin, 2015). In this article, we propose a robust
estimation of sufficient dimension reduction the independence via the
$\alpha$-distance covariance ($\alpha$-dCov) between the response and the
predictors and develop a new algorithm for estimating directions in general
multiple-index models with a form

$Y=g(\bm{\beta}^{T}X,\epsilon),$

where $g$ is an unknown link function (Yin ., 2008; Xia, 2008; Sheng Yin,
2013). The rest of this article is organized as follows: Section 2 describes
our robust $\alpha$-dCov method and a corresponding outlier detection method,
including motivation, theoretical results, estimation algorithm, and testing
procedure. We introduce the consistency theorem in Section 3. Section 4
contains simulation and real data studies. We summarize our work in Section 5.

###  1.1 Generalized distance covariance

Distance covariance (Székely ., 2007) is a popular dependence measure for two
random vectors of possibly different dimensions and types. In recent years,
there have been concentrated efforts in the literature to understand the
distributional properties of the sample distance covariance in a high-
dimensional setting, with an exclusive emphasis on the null case that $X$ and
$Y$ are independent. Distance covariance can be generalized to include powers
of Euclidean distance. Define

$\displaystyle\nu^{2}(X,Y;\alpha)$
$\displaystyle:=\mbox{E}[\|X-X^{\prime}\|^{\alpha}\,\|Y-Y^{\prime}\|^{\alpha}]+\mbox{E}[\|X-X^{\prime}\|^{\alpha}]\,\mbox{E}[\|Y-Y^{\prime}\|^{\alpha}]$
$\displaystyle-2\mbox{E}[\|X-X^{\prime}\|^{\alpha}\,\|Y-Y^{\prime\prime}\|^{\alpha}],$
(1)

where $(X,Y),\;(X^{\prime},Y^{\prime}),\;(X^{\prime\prime},Y^{\prime\prime})$
are independent and identically distributed (i.i.d.) with respect to the joint
distribution of $(X,Y)$ (Székely Rizzo, 2014). As discussed in (Székely .,
2007), for every $0<\alpha<2$, $X$ and $Y$ are independent if and only if
$\nu^{2}(X,Y;\alpha)=0$. When $\alpha=1$, it reduces to the classical distance
covariance. When $0<\alpha<1$, it can be considered as a more robust version
of distance covariance as it reduces the influence of large values of
$\|X-X^{\prime}\|$, $\|Y-Y^{\prime}\|$, and $\|Y-Y^{\prime\prime}\|$ that
might be contributed to outliers.

###  1.2 Central Space Estimation via $\alpha$-dCov

Let $({\bf X},{\bf Y})=\\{(X_{i},Y_{i}):\;i=1,\ldots,n\\}$ be $n$ random
samples from random variables $(X,Y)$. In addition, ${\bf X}$ denotes a
$p\times n$ data matrix whose columns are $X_{1},\cdots,X_{n}$ and ${\bf
Y}=[Y_{1},\cdots,Y_{n}]$ denotes a $1\times n$ response data matrix. In this
article, we consider univariate responses. However, the method can naturally
be extended to multivariate responses without any issue due to the nature of
$\alpha$-dCov. The empirical solution of the SDR method based on $\alpha$-dCov
for these $n$ observations relies on solving the following objective function
(Székely ., 2007; Sheng Yin, 2016):

$\displaystyle\max_{\bm{\beta}\in\mathbb{R}^{p\times
d}}\nu^{2}_{n}(\bm{\beta}^{T}{\bf X},{\bf Y},\alpha).$ (2)

with constraint $\bm{\beta}^{T}\Sigma_{X}\bm{\beta}={\bf I}_{d}$ and $1\leq
d\leq p$, where $\nu_{n}$ is the empirical version of $\nu$ defined in
Equation (1.1).

The empirical distance dependence statistics $\nu_{n}$ is defined as follows.
For $k,l=1,\ldots,n$, we compute the Euclidean distance matrices
$(a_{kl})=(|X_{k}-X_{l}|^{\alpha}_{p})$ and
$(b_{kl})=(|Y_{k}-Y_{l}|^{\alpha})$ for $0<\alpha<2$ (Székely Rizzo, 2009).
Define

$A_{kl}=a_{kl}-\bar{a}_{k\cdot}-\bar{a}_{\cdot l}+\bar{a}_{\cdot\cdot},\quad
k,l=1,\ldots,n,$

where

$\bar{a}_{k\cdot}=\frac{1}{n}\sum^{n}_{l=1}a_{kl},\quad\bar{a}_{\cdot
l},=\frac{1}{n}\sum^{n}_{k=1}a_{kl},\quad\bar{a}_{\cdot\cdot}=\frac{1}{n^{2}}\sum^{n}_{k,l=1}a_{kl}.$

Similarly, define $B_{kl}=b_{kl}-\bar{b}_{k\cdot}-\bar{b}_{\cdot
l}+\bar{b}_{\cdot\cdot}$, for $k,l=1,\ldots,n$. The nonnegative sample
distance covariance $\nu_{n}({\bf X},{\bf Y})$ and sample distance correlation
$R_{n}({\bf X},{\bf Y})$ are defined by

$\nu^{2}_{n}(\bm{\beta}^{T}{\bf X},{\bf
Y},\alpha)=\frac{1}{n^{2}}\sum^{n}_{k,l=1}A_{kl}B_{kl}$ (3)

and

$R^{2}_{n}({\bf X},{\bf Y},\alpha)=\begin{cases}\frac{\nu^{2}_{n}({\bf X},{\bf
Y},\alpha)}{\nu^{2}_{n}({\bf X},\alpha)\nu^{2}_{n}({\bf Y},\alpha)},\;&\mbox{
if }\nu^{2}_{n}({\bf X},\alpha)\nu^{2}_{n}({\bf Y},\alpha)>0;\\\ 0,\;&\mbox{
if }\nu^{2}_{n}({\bf X},\alpha)\nu^{2}_{n}({\bf Y},\alpha)=0,\end{cases}$

respectively, where the sample distance variance is defined by

$\nu^{2}_{n}({\bf X},\alpha):=\nu^{2}_{n}({\bf X},{\bf
X},\alpha)=\frac{1}{n^{2}}\sum^{n}_{k,l=1}A^{2}_{kl}.$

Following Wu Chen (2021), we have the following equivalence. Let ${\bf
C}=\hat{\Sigma}^{\frac{1}{2}}_{X}\bm{\beta}$ and ${\bf
Z}=\hat{\Sigma}^{-\frac{1}{2}}_{X}{\bf X},$ the target function (2) can be
rewritten as

$\displaystyle\max_{{\bf C}}\nu^{2}_{n}({\bf C}^{T}{\bf Z},{\bf
Y},\alpha):=\frac{1}{n^{2}}\sum^{n}_{k,l=1}a_{kl}({\bf C})B_{kl},\mbox{ s.t.
}{\bf C}\in\mbox{St}(d,p),$ (4)

where $a_{kl}({\bf C})=\|{\bf C}^{T}Z_{k}-{\bf C}^{T}Z_{l}\|^{\alpha}$. We use
the same notation $\mbox{St}(d,p)=\\{{\bf C}\in\mathbb{R}^{p\times d}\mid{\bf
C}^{T}{\bf C}=I_{d}\\}$ with $d\leq p$ is referred to the Stiefel manifold and
$T_{{\bf C}}\mbox{St}(d,p)$ is the tangent space to $\mbox{St}(d,p)$ at a
point ${\bf C}\in\mbox{St}(d,p)$. We assume that $Y=g({\bf C}^{T}{\bf
Z},\epsilon)$, where ${\bf C}$ is a $p\times d$ matrix, $\epsilon$ is an
unknown random error independent of ${\bf Z}$, and $g$ is an unknown link
function. We propose a new method to estimate a basis of the central subspace
$S_{Y\mid{\bf Z}}=\mbox{Span}({\bf C})$ and denote $\nu^{2}_{n}({\bf
C}^{T}{\bf Z},{\bf Y},\alpha)$ as $F({\bf C})$.

##  2 Algorithm

We develop an iterative algorithm based on the gradient descent algorithm on
the Stiefel manifold. Here $P_{S}$ is a projection on the Stiefel manifold
(Dalmau-Cedeno Oviedo, 2017). By Proposition 3.4 (the projection onto Stiefel
manifolds) of PA. Absil Malick (2012), we let $\bar{{\bf C}}\in\mbox{St}(d,p)$
for any ${\bf C}$ such that $\|X-\bar{X}\|<\sigma_{d}(\bar{{\bf C}})$, where
$\sigma_{d}(\bar{{\bf C}})$ is the largest singular value, then the projection
of ${\bf C}$ onto $\mbox{St}(d,p)$ exists uniquely, and can be expressed as
$P_{S}({\bf C})=\sum^{d}_{i=1}u_{i}v_{i}^{T}$, given by a singular value
decomposition of ${\bf C}$. Alternatively, let the SVD of ${\bf
C}\in\mathbb{R}^{p\times d}$ be ${\bf C}={\bf U}\Sigma{\bf V}$, then
$P_{S}({\bf C})={\bf U}{\bf V}^{T}$.

Algorithm 1 rSDR: robust SDR

1: Input: The samples $\\{(y_{i},{\bf Z}_{i}),i=1,\cdots,n\\}$, initial ${\bf
C}^{(0)}$.

2: Initialization: ${\bf C}^{(0)}$.

3: for $\mathrm{iter}=0,1,\cdots$ do

4: Let ${\bf C}^{(\mathrm{iter}+1)}=P_{S}({\bf
C}^{(\mathrm{iter})}+\alpha^{(\mathrm{iter})}_{1}\partial_{\bf C}F({\bf
C}^{(\mathrm{iter})}))$ or ${\bf C}^{(\mathrm{iter}+1)}=P_{S}\Big{(}{\bf
C}^{(\mathrm{iter})}+\alpha^{(\mathrm{iter})}_{1}\partial_{\bf C}F({\bf
C}^{(\mathrm{iter})})({\bf I}-{\bf C}^{(\mathrm{iter})\,T}{\bf
C}^{(\mathrm{iter})})\Big{)}$, where $P_{S}(\cdot)$ is the projection on the
Stiefel manifold and $\frac{\partial}{\partial{\bf C}}F({\bf C})$, and
$\alpha^{(\mathrm{iter})}_{1}$ is chosen by a line search.

5: Repeat steps 4 until $\|F({\bf C}^{(\mathrm{iter})})-F({\bf
C}^{(\mathrm{iter}-1)})\|_{F}\leq{\epsilon}_{n}$ where ${\epsilon}_{n}$ is a
pre-specified threshold, or the number of iterations exceeds the upper limit:
$\mathrm{iter}>N^{(\max)}$.

6: end for

7: Output: Estimated coefficients $\hat{{\bf C}}$.

Now we derive the explicit formula for $\partial_{\bf C}F({\bf C})$, where
$F({\bf C})=\nu^{2}_{n}({\bf C}^{T}{\bf
Z},Y)=\frac{1}{n^{2}}\sum^{n}_{k,l=1}a_{kl}({\bf C})B_{kl}$. Recall that
$a_{kl}({\bf C})=\|{\bf C}^{T}{\bf Z}_{k}-{\bf C}^{T}{\bf Z}_{l}\|^{\alpha}$,
the gradient is

$\partial_{\bf C}F({\bf C})=\frac{\alpha}{n^{2}}\sum^{n}_{k,l=1}\frac{{\bf
C}^{T}({\bf Z}_{k}-{\bf Z}_{l})({\bf Z}_{k}-{\bf Z}_{l})^{T}}{\|{\bf
C}^{T}{\bf Z}_{k}-{\bf C}^{T}{\bf Z}_{l}\|^{2-\alpha}}B_{kl},$ (5)

and one may perform the manifold gradient descent algorithm as follows:

${\bf C}^{(\mathrm{iter}+1)}=P_{S}\Big{(}{\bf
C}^{(\mathrm{iter})}+\alpha^{(\mathrm{iter})}_{1}\partial_{\bf C}F({\bf
C}^{(\mathrm{iter})})({\bf I}-{\bf C}^{(\mathrm{iter})}{\bf
C}^{(\mathrm{iter})\,T})\Big{)}.$

We remark that while there are various advanced Stiefel manifold optimization
algorithms such as the ones based on the Cayley transform (Wen Yin, 2013; Zhu
., 2019) or geodesics (P. Absil ., 2009), we applied the standard projected
gradient descent algorithm as it is simpler to implementation and has the same
order of computational cost per iteration of $O(p^{2}d)$.

Implementation issues When implementing our approach, practical challenges may
arise due to the potential for an extremely small denominator in Equation (5),
disproportionately amplifying the influence of the $(k,l)$-th term. To
preemptively address this concern, we introduce a small positive
regularization parameter, denoted as $\eta$. Subsequently, we employ a
regularization technique on the objective function $F({\bf C})$ such that the
$(k,l)$-th term of the gradient in Equation (5) remains bounded. In
particular, we apply it to the regularized objective function denoted as
$F_{\eta}({\bf C})$:

$F_{\eta}({\bf C})=\frac{1}{n^{2}}(\|{\bf C}^{T}Z_{k}-{\bf
C}^{T}Z_{l}\|^{2}+\eta)^{\alpha/2}B_{kl},$

which leads to the regularized gradient formulation expressed as follows:

$\frac{\alpha}{n^{2}}\partial_{\bf C}F_{\eta}({\bf
C})=\sum^{n}_{k,l=1}\frac{{\bf C}^{T}({\bf Z}_{k}-{\bf Z}_{l})({\bf
Z}_{k}-{\bf Z}_{l})^{T}}{(\|{\bf C}^{T}{\bf Z}_{k}-{\bf C}^{T}{\bf
Z}_{l}\|^{2}+\eta)^{(2-\alpha)/2}}B_{kl}.$ (6)

##  3 Consistency Theory

We consider a model with a general noise term

${\bf Y}=g(\bm{\beta}_{0}^{T}{\bf X},\epsilon)=g({\bf C}_{0}^{T}{\bf
Z},\epsilon),$

where $\bm{\beta}_{0}$ is a $p\times d$ orthogonal matrix, $g(\cdot)$ is an
unknown link function, ${\bf
C}_{0}=\hat{\Sigma}^{\frac{1}{2}}_{X}\bm{\beta}_{0}$, and ${\bf
Z}=\hat{\Sigma}^{-\frac{1}{2}}_{X}{\bf X}$, and $\epsilon$ is independent of
${\bf Z}$. This model includes the model from Xia . (2002) that

${\bf Y}=g(\bm{\beta}_{0}^{T}{\bf X})+\epsilon$

as a special example.

Following Sheng Yin (2016), we have the asymptotic properties of the estimator
$\hat{{\bf C}}$ that is consistent. The statement and the proof is similar to
that of Sheng Yin (2013). It requires an additional assumption that depends on
the decomposition of $X$ into two independent components, and some discussions
on this condition are available in (Sheng Yin, 2013, Section 3.2). For
example, it is satisfied when $X$ is normal (N. Zhang Yin, 2015). In addition,
this assumption also holds asymptotically when $p$ is large (Hall Li, 1993).

The following proposition establishes the asymptotic properties of our
estimator ${\bf C}$ up to some rotation matrix ${\bf Q}$. This implies the
asymptotic property of the estimated central subspace as it is invariant to
the rotation matrix.

######  Proposition 3.1.

Let ${\bf C}\in\mathbb{R}^{d\times p}$ be a basis of the central subspace
$S_{Y\mid X}$ with ${\bf C}^{T}\Sigma_{X}{\bf C}=I_{d}$. Suppose $P^{T}_{{\bf
C}(\Sigma_{X})}X\perp\\!\\!\\!\perp Q^{T}_{{\bf C}(\Sigma_{X})}X$ and the
support of $X\in\mathbb{R}^{d\times p}$, say $S$, is a compact set. In
addition, assume that there exists ${\bf C}^{\prime}\in\mathbb{R}^{(p-d)\times
p}$ such that $[{\bf C},{\bf C}^{\prime}]^{T}\Sigma_{X}[{\bf C},{\bf
C}^{\prime}]=I_{p}$ and ${\bf C}^{T}X$ is independent of ${\bf C}^{\prime
T}X$. Let $\hat{{\bf C}}=\arg\min_{{\bf C}^{T}\Sigma_{X}{\bf
C}=I_{d}}\nu^{2}_{n}({\bf C}^{T}{\bf X},{\bf Y})$, then there exists a
rotation matrix ${\bf Q}$: ${\bf Q}^{T}{\bf Q}=I_{d}$ such that $\hat{{\bf
C}}\stackrel{{\scriptstyle P}}{{\to}}{\bf C}{\bf Q}$ (convergence in
probability) as $n\to\infty$.

###### Proof.

Following (Székely Rizzo, 2009, (4.1)), we have that for random variables $X$
and $Y$ from ${\mathbb{R}}^{p_{1}}$ and ${\mathbb{R}}^{p_{2}}$,

$\nu^{2}(X,Y,\alpha)=C\int_{t,s}\frac{|f_{X,Y}(t,s)-f_{X}(f)f_{Y}(t)|^{2}}{\|t\|^{p_{1}+\alpha}\|s\|^{p_{2}+\alpha}}{\mathrm{d}}t{\mathrm{d}}s,$

where $f_{X}$, $f_{Y}$, $f_{X,Y}$ represent the characteristic functions of
$X$, $Y$, and $(X,Y)$ respectively.

The rest follows from the proof of Proposition 1 in N. Zhang Yin (2015). For
any $\bm{\beta}\neq{\bf C}$ that satisfies
$\bm{\beta}^{T}\Sigma_{X}\bm{\beta}={\bf I}$, let $\bm{\beta}_{1}$ be the
projection of $\bm{\beta}$ to the subspace spanned by ${\bf C}$ with an inner
product induced by $\Sigma_{X}$ (that is, $\Sigma_{X}^{0.5}\bm{\beta}_{1}$
being the projection of $\Sigma_{X}^{0.5}\bm{\beta}$ to the subspace spanned
by $\Sigma_{X}^{0.5}{\bf C}$ under the Euclidean metric) and
$\bm{\beta}_{2}=\bm{\beta}-\bm{\beta}_{1}$. Then since
$\Sigma_{X}^{0.5}\bm{\beta}$ and $\Sigma_{X}^{0.5}{\bf C}$ are both orthogonal
subspaces, we have $\|{\bf
C}^{\dagger}\bm{\beta}_{1}\|=\|(\Sigma_{X}^{0.5}{\bf
C})^{\dagger}(\Sigma_{X}^{0.5}\bm{\beta}_{1})\|\leq\|(\Sigma_{X}^{0.5}{\bf
C})^{\dagger}(\Sigma_{X}^{0.5}\bm{\beta})\|\leq 1$, where ${\dagger}$
represents the pseudo inverse. Note that $\bm{\beta}_{1}$ and ${\bf C}$ have
the same column space, so for any ${\bf z}\in{\mathbb{R}}^{p}$, we have

$\|\bm{\beta}_{1}{\bf z}\|\leq\|{\bf C}{\bf z}\|.$ (7)

Then we proved that ${\bf C}$ is the solution to $\arg\min_{{\bf
C}^{T}\Sigma_{X}{\bf C}=I_{d}}\nu^{2}({\bf C}^{T}{\bf X},{\bf Y})$
asymptotically:

$\displaystyle\nu^{2}(\bm{\beta}^{T}X,Y,\alpha)=\int|\mathrm{E}e^{i\langle
t,\bm{\beta}^{T}X\rangle+i\langle s,Y\rangle}-\mathrm{E}e^{i\langle
t,\bm{\beta}^{T}X\rangle}\mathrm{E}e^{i\langle
s,Y\rangle}|^{2}\big{/}(\|t\|^{d+\alpha}\|s\|^{1+\alpha}){\mathrm{d}}t{\mathrm{d}}s$
$\displaystyle=$ $\displaystyle\int|\mathrm{E}e^{i\langle
t,\bm{\beta}_{2}^{T}X\rangle}|^{2}|\mathrm{E}e^{i\langle
t,\bm{\beta}_{1}^{T}X\rangle+i\langle s,Y\rangle}-\mathrm{E}e^{i\langle
t,\bm{\beta}_{1}^{T}X\rangle}\mathrm{E}e^{i\langle
s,Y\rangle}|^{2}\big{/}(\|t\|^{d+\alpha}\|s\|^{1+\alpha}){\mathrm{d}}t{\mathrm{d}}s$
$\displaystyle\leq$ $\displaystyle\int|\mathrm{E}e^{i\langle
t,\bm{\beta}_{1}^{T}X\rangle+i\langle s,Y\rangle}-\mathrm{E}e^{i\langle
t,\bm{\beta}_{1}^{T}X\rangle}\mathrm{E}e^{i\langle
s,Y\rangle}|^{2}\big{/}(\|t\|^{d+\alpha}\|s\|^{1+\alpha}){\mathrm{d}}t{\mathrm{d}}s$
$\displaystyle=$
$\displaystyle\nu^{2}_{n}(\bm{\beta}_{1}^{T}X,Y,\alpha)\leq\nu^{2}({\bf
C}^{T}X,Y,\alpha),$

where the last step follows from Equation (7). It is easy to verify that the
equality only holds when $\bm{\beta}={\bf C}Q$ for some rotation matrix $Q$.

It remains to shows that $\nu^{2}_{n}({\bf C}^{T}{\bf X},{\bf Y},\alpha)$ is
the empirical estimate of the random variable $\nu^{2}({\bf
C}^{T}X,Y,\alpha)$, which means that $\nu^{2}_{n}({\bf C}^{T}{\bf X},{\bf
Y},\alpha)\stackrel{{\scriptstyle a.s.}}{{\to}}\nu^{2}({\bf C}^{T}{\bf X},{\bf
Y},\alpha)$ (almost sure convergence) as $n\to\infty$. The result holds
following the proof of Lemma 2 in the supplementary material of N. Zhang Yin
(2015). ∎

###  3.1 Convergence analysis

We investigate the convergence property of the proposed algorithm in this
section. In fact, the proposed algorithm generates solutions that converge to
a stationary point of $F_{\eta}({\bf C})$ as $t\to\infty$. In addition, the
algorithm converges to the solution when well-initialized.

######  Theorem 3.2.

(a) Any accumulation point of the sequence $\left\\{\hat{{\bf
C}}^{(t)}\right\\}_{t\geq 0}$ generated by the proposed algorithm converges is
a stationary point of $F_{\eta}({\bf C})$ over the Stiefel manifold.

(b) If in addition, the global maximizer $\hat{{\bf C}}$ it is the unique
stationary point in its neighborhood $\mathcal{N}$, and $F_{\eta}({\bf
C})-F_{\eta}(\hat{{\bf C}})\leq-c\|{\bf C}-\hat{{\bf C}}\|_{F}^{2}$ for any
${\bf C}$ in $\mathcal{N}$ and some $c>0$. Then when the initialization
$\hat{{\bf C}}^{(0)}$ is sufficiently close to $\hat{{\bf C}}$, the sequence
$\left\\{\hat{{\bf C}}^{(t)}\right\\}_{t\geq 0}$ converges to $\hat{{\bf C}}$.

###### Proof.

(a) Due to the line search strategy in Algorithm 1, the objective value of the
objective function is monotonically nondecreasing and as a result,
$\nu^{2}_{n}(\hat{{\bf C}}^{(t)\,T}{\bf X},{\bf Y},\alpha)$ converges. Let
$\tilde{{\bf C}}$ be any accumulation point of the sequence $\hat{{\bf
C}}^{(t)}$, then $\nabla_{{\bf C}}\nu^{2}_{n}({\bf C}^{T}{\bf X},{\bf
Y},\alpha)|_{{\bf C}=\hat{{\bf C}}^{(t)}}=0$, since otherwise the objective
function will continue to increase.

(b) Since the gradient of $F_{\eta}({\bf C})$ is continuous, $\max_{{\bf
C}:\|{\bf C}-\hat{{\bf C}}\|_{F}\leq\epsilon}\|F_{\eta}({\bf C})\|$ converges
to zero as $\epsilon\rightarrow 0$. As a result, we may choose
$\epsilon^{\prime}>0$ such that for

$\mathcal{N}_{\epsilon^{\prime}}=\mathcal{N}\cap\\{F_{\eta}({\bf
C})-F_{\eta}(\hat{{\bf C}})>-\epsilon^{\prime}\\},$

and any $\hat{{\bf C}}^{(t)}\in\mathcal{N}_{\epsilon^{\prime}}$, $\|\hat{{\bf
C}}^{(t)}-\hat{{\bf C}}\|_{F}\leq\sqrt{\epsilon^{\prime}/c}$ and the gradient
$F_{\eta}^{\prime}(\hat{{\bf C}}^{(t)})$ is so small such that the next
iteration $\hat{{\bf C}}^{(t+1)}$ remains in $\mathcal{N}$. Since the
functional value $F_{\eta}(\hat{{\bf C}}^{(t)})$ is nonincreasing, $\hat{{\bf
C}}^{(t+1)}$ lies in $\mathcal{N}_{\epsilon^{\prime}}$ as well. As $\hat{{\bf
C}}$ is the unique stationary point in $\mathcal{N}_{\epsilon^{\prime}}$, part
(a) implies that the algorithm converges to $\hat{{\bf C}}$.

∎

##  4 Numerical Studies

In this section, we perform a comparative analysis of several algorithms
including the proposed robust SDR (rSDR), the SQP algorithm (Sheng Yin,
2013),the MMRN algorithm (Wu Chen, 2021) and the HSIC algorithm (N. Zhang Yin,
2015).

The problem in Equation (4) is nonlinear and the proposed algorithm, rSDR,
needs an good initialization. The solutions of the sliced inverse regression
(SIR, KC. Li (19912)) and the directional regression (DR, B. Li Wang (2007))
are used in the initialization of Algorithm 1. Let $\bm{\beta}_{1}$ and
$\bm{\beta}_{2}$ be two solutions of SDR obtained by SIR and DR, respectively.
We select one of $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$ with larger dCov as our
initial value of $\bm{\beta}$. Let $\hat{\bm{\Sigma}}$ be the sample
covariance of $\\{{\bm{x}}\\}_{i=1}^{n}$. The initial matrix ${\bf
C}^{(0)}=\hat{\bm{\Sigma}}^{1/2}_{X}\bm{\beta}$ is evaluated in Algorithm 1.

The proposed algorithm has an parameter $\alpha$ which governs robustness to
outliers. A smaller $\alpha$ usually enhances the robustness of Algorithm 1.
However, an excessively small $\alpha$ often results in numerous local minimum
values for the problem. Therefore, $\alpha$ is tuned through 5-fold cross-
validation. The value of $\alpha$ is fine-tuned from $\\{i/10\\}_{i=1}^{9}$ by
5-fold CV. Specifically, we partition the datasets
$\\{(y_{i},{\bm{x}}_{i})\\}_{i=1}^{n}$ into training and validation sets. For
each $\alpha$ value, we apply Algorithm 1 to the training set, yielding a
subspace $\bm{\beta}_{\alpha}$. We then assess the $0.5$-dCov of the
validation set. This process is repeated for all 5 folds, and the average
$0.5$-dCov is computed. We choose the $\alpha$ value associated with the
highest average and execute Algorithm 1 again to derive the estimated
subspace. It is important to note that if the dataset is contaminated with
outliers, the validation set will also contain outliers. Traditional dCov or
covariance calculations may be significantly impacted by these outliers.
Therefore, opting for a more robust variance statistic is crucial. In this
context, we select the $0.5$-dCov as the measure for the test set.

The SQP algorithm utilizes sequential quadratic programming to solve the dCov-
based SDR model (equivalent to Equation (4) with $\alpha=1$). While the SQP
method performs well when the dimension ($p$) and sample size ($n$) are
relatively small, it becomes computationally difficult for moderately high-
dimensional settings (Wu Chen, 2021). MMRN was later proposed as an efficient
alternative to solve the same model using Riemannian Newton’s method. Both SQP
and MMRN correspond to rSDR with $\alpha=1$, but none of them is robust
against outliers. The Hilbert-Schmidt Independence Criterion (HSIC) method (N.
Zhang Yin, 2015) addresses the single-index SDR model ($d=1$) by maximizing
the HSIC covariance between $\bm{\beta}^{T}{\bf X}$ and ${\bf Y}$.

In the first simulation, we compare rSDR with SQP and MMRN in both robust and
non-robust settings. Our results demonstrate that rSDR with a smaller $\alpha$
can effectively estimate the underlying subspace and efficiently solve the SDR
model. Additionally, even in the presence of outliers in the response, rSDR
can still estimate the subspace accurately, while SQP and MMRN fail to do so.

In the second simulation, we explore the application of rSDR in outlier
detection. By reducing the data dimension, we extend a dCor-based outlier
detection method (Wang Li, 2017) to high-dimensional cases. We compare rSDR
with PCA in dimensionality reduction and outlier detection to showcase the
applicability of robust SDR in outlier detection.

Furthermore, we present three real data examples: the New Zealand horse
mussels, cardiomyopathy microarray data, and auto MPG data. In the New Zealand
horse mussels dataset, we reduce the data dimension to 1 and compare rSDR with
HSIC. Notably, HSIC is only applicable when $d=1$, so we do not include it in
other simulations or real data examples.

###  4.1 Simulation Data

Let
$\tilde{\beta}_{1}=(1,0,0,\cdots,0)^{T},\tilde{\beta}_{2}=(0,1,0,\cdots,0)^{T},\tilde{\beta}_{3}=(1,0.5,1,\cdots,0)^{T}$
be three $p$-dimensional vectors. We further rotate the vectors
$\tilde{\beta}_{i}$ by a random rotation matrix $R_{d}\in SO(p)$ (the special
orthogonal group of dimension $p$), i.e.,
$\beta_{i}=R_{d}^{\top}\tilde{\beta}_{i}$. We consider the following three
models

  1. (A)

$Y=(\beta_{1}^{T}X)^{2}+(\beta_{2}^{T}X)+0.1\epsilon$,

  2. (B)

$Y=\mathrm{sign}(2\beta_{1}^{T}X+\epsilon_{1})\times\log|2\beta_{2}^{T}X+4+\epsilon_{2}|$,

  3. (C)

$Y=\exp(\beta_{3}^{T}X)\epsilon$,

where $X\in{\mathbb{R}}^{p}$ follows from (1) ${\mathcal{N}}(0,{\mathbf{I}})$
and (2) $U[-2,2]^{p}$ and $\epsilon,\epsilon_{1},\epsilon_{2}$ are standard
normal distributed. We analyze the principal angles between the true subspace
$\bm{\beta}$ and the estimated subspace $\hat{\bm{\beta}}$ obtained using
different SDR methods, namely rSDR, MMRN, and SQP. To further investigate the
robustness of these methods, we introduce additional noise by adding the
response with a value of $50\times\bm{1}^{T}X$ with a probability of 0.1. We
then calculate the principal angles between the true subspace and the
estimated subspaces in this robust setting. Both simulation scenarios are
conducted for two settings: $(n,p)=(100,6)$ and $(n,p)=(500,20)$. We repeat
the simulations 100 times and report the mean and standard deviation of the
principal angles for both the non-robust and robust cases in Tables 4.1 and
4.1, respectively. It is worth noting that the underlying subspace for model
(A) and (B) is represented by $\bm{\beta}=[\beta_{1},\beta_{2}]$, resulting in
a value of $d=2$. On the other hand, the underlying subspace for model (C) is
represented by $\bm{\beta}=\beta_{3}$, resulting in a value of $d=1$.

From Table 4.1 we observe that rSDR performs better than MMRN and SQP in model
(A) and (B) even without outliers. MMRN converges faster than SDR in model (A)
and (B). When $(n,p)=(500,20)$, rSDR and MMRN are faster than SQP. Table 4.1
reports the principal angles and execution time of the three estimators in the
scenario where the outliers present. Table 4.1 shows that the principal angles
between the true subspace and the estimated subspace produced by rSDR are
smaller than MMRN and SQP which implies that rSDR is more robust. Moreover,
rSDR converges faster than MMRN and SQP in most settings; particularly in
model (C).

The mean and standard deviation (in parentheses) of the principal angles and
the running times (seconds) over 100 repetitions of SQP, MMRN and rSDR in
nonrobust settings.  $(n,p)$ Model SQP MMRN rSDR Angle Time(s) Angle Time(s)
Angle Time(s) (100,6) A(1) 0.27(0.09) 0.16(0.12) 0.27(0.09) 0.19(0.09)
0.27(0.09) 0.17(0.13) A(2) 0.25(0.08) 0.13(0.13) 0.25(0.08) 0.15(0.12)
0.25(0.08) 0.18(0.13) B(1) 0.28(0.09) 0.10(0.02) 0.28(0.09) 0.19(0.12)
0.28(0.09) 0.20(0.14) B(2) 0.22(0.08) 0.11(0.04) 0.22(0.08) 0.32(0.59)
0.21(0.08) 0.23(0.18) C(1) 0.20(0.07) 0.24(0.32) 0.20(0.07) 0.25(0.11)
0.19(0.06) 0.08(0.05) C(2) 0.32(0.12) 0.14(0.24) 0.31(0.12) 0.38(0.17)
0.32(0.12) 0.08(0.05) (500,20) A(1) 0.24(0.04) 2.98(0.56) 0.24(0.04)
0.90(0.14) 0.24(0.04) 1.41(0.77) A(2) 0.23(0.04) 3.33(3.11) 0.23(0.04)
0.90(1.36) 0.23(0.04) 1.65(0.80) B(1) 0.24(0.04) 3.17(0.63) 0.24(0.04)
0.90(0.14) 0.24(0.04) 1.54(0.91) B(2) 0.19(0.03) 4.55(1.29) 0.19(0.03)
0.81(0.13) 0.18(0.03) 1.56(0.79) C(1) 0.16(0.03) 2.54(0.22) 0.16(0.03)
1.52(0.26) 0.17(0.03) 0.66(0.42) C(2) 0.25(0.04) 3.27(0.71) 0.25(0.04)
4.04(1.23) 0.28(0.05) 0.72(0.49)

The mean and standard deviation (in parentheses) of the principal angle, and
the running time (seconds) over 100 repetitions of SQP, MMRN and rSDR in
robust settings.  $(n,p)$ Model SQP MMRN rSDR Angle Time(s) Angle Time(s)
Angle Time(s) (100,6) A(1) 0.51(0.28) 0.20(0.21) 0.49(0.27) 0.28(0.20)
0.32(0.12) 0.19(0.13) A(2) 0.45(0.28) 0.13(0.12) 0.44(0.27) 0.26(0.19)
0.27(0.11) 0.17(0.12) B(1) 0.52(0.26) 0.13(0.08) 0.51(0.26) 0.32(0.26)
0.33(0.12) 0.18(0.13) B(2) 0.42(0.22) 0.12(0.07) 0.42(0.22) 0.25(0.22)
0.24(0.08) 0.21(0.17) C(1) 0.39(0.24) 0.19(0.21) 0.38(0.23) 0.32(0.16)
0.26(0.10) 0.08(0.06) C(2) 0.47(0.23) 0.17(0.27) 0.46(0.22) 0.51(0.25)
0.40(0.16) 0.09(0.06) (500,20) A(1) 0.82(0.30) 4.12(1.29) 0.82(0.30)
2.92(1.47) 0.25(0.04) 1.47(0.87) A(2) 0.93(0.42) 16.26(42.21) 0.92(0.42)
6.85(13.70) 0.24(0.04) 1.43(0.85) B(1) 0.91(0.35) 4.20(1.48) 0.90(0.35)
3.36(2.23) 0.26(0.04) 1.74(0.92) B(2) 0.60(0.36) 5.32(2.26) 0.60(0.36)
2.84(3.54) 0.19(0.03) 1.64(0.99) C(1) 0.35(0.14) 3.02(0.34) 0.35(0.14)
4.41(2.15) 0.22(0.04) 0.63(0.46) C(2) 0.89(0.28) 4.53(1.01) 0.85(0.29)
13.39(6.27) 0.35(0.07) 0.88(0.61)

###  4.2 Outlier Detection Simulation Studies

Our proposed SDR method can be effectively utilized for outlier detection.
Wang Li (2017) introduced a novel outlier detection measure based on the
distance correlation (dCor) given by

$\displaystyle\mathcal{D}_{i}({\bf X},{\bf
Y})=\frac{1}{p}\sum_{k=1}^{p}\left(\mathrm{dCor}({\bf X}_{k},{\bf
Y})-\mathrm{dCor}({\bf X}^{(i)}_{k},{\bf Y}^{(i)})\right)^{2},$ (8)

where $\mathrm{dCor}({\bf X}_{k},{\bf Y})$ represents the dCor between the
$k$-th predictor and the response ${\bf Y}$. The dCor between ${\bf X}$ and
${\bf Y}$ is defined as

$\displaystyle\mathrm{dCor}^{2}({\bf X},{\bf Y})=\frac{\mathrm{dCov}^{2}({\bf
X},{\bf Y})}{\sqrt{\mathrm{dCov}^{2}({\bf X},{\bf X})\mathrm{dCov}^{2}({\bf
Y},{\bf Y})}}.$

It is evident that if the $i$-th data point $({\bf X}^{(i)}_{k},{\bf
Y}^{(i)})$ exhibits a high value of the measure $\hat{\mathcal{D}}_{i}$, it is
more likely to be an outlier observation. The method employs a bootstrap
procedure to determine the threshold $\hat{F}_{\gamma}$. At a given
significance level $\gamma$, the $i$-th observation is identified as an
outlier if $\hat{\mathcal{D}}_{i}>\hat{F}_{\gamma}$, where $\hat{F}_{\gamma}$
represents the upper $\gamma$-th quantile of the cumulative distribution
function of $\mathcal{D}_{i}$ under the null hypothesis. Specifically, a
bootstrap sample $\mathcal{D}^{[b]}_{i}$ is formed by drawing with replacement
from ${1,\cdots,n}$, denoted as ${i^{[b]}_{(1)},\cdots,i^{[b]}_{(n)}}$, and an
estimator $\hat{\mathcal{D}}^{[b]}_{i}$ is computed for each sample. The
threshold $\hat{F}_{\gamma}$ is determined by calculating the upper
$\gamma$-th quantile of the cumulative distribution function of
$\hat{\mathcal{D}}^{[b]}_{i}$.

The algorithm proposed by Wang Li (2017), which is based on the outlier
detection measure defined in Equation (8), involves calculating the covariance
distance between ${\bf X}$ and ${\bf Y}$ in each dimension and with the
removal of each sample. As a result, its computational complexity is
$O(pn^{3})$, where the computation of dCov requires pairwise distance
calculations between the columns of ${\bf X}$ and ${\bf Y}$. A natural
approach to enhance their method is to reduce the dimensionality of the
dataset ${\bf X}$. Their method can be naturally extended to detect outlier
locations by computing

$\displaystyle\mathcal{D}_{i}(\bar{{\bf X}},{\bf
Y})=\frac{1}{p}\sum_{k=1}^{p}\left(\mathrm{dCor}(\bar{{\bf X}}_{k},{\bf
Y})-\mathrm{dCor}(\bar{{\bf X}}^{(i)}_{k},{\bf Y}^{(i)})\right)^{2},$ (9)

where $\bar{{\bf X}}\in{\mathbb{R}}^{d\times n}$ is the $d$-dimensional data
obtained by dimension reduction. Nevertheless, the conventional approach to
dimension reduction is unsuitable in the presence of outliers. Therefore, we
employ the robust SDR as a means to both reduce the data’s dimensionality and
identify outlier positions. For the sake of comparison, we also implement
principal component analysis (PCA) (Wold ., 1987) for dimension reduction.

We consider an autoregressive correlation structure with
$\bm{\Sigma}=(\rho_{j,k})_{p\times p}=0.5^{|j-k|}$ and generate the data as
follows: $X_{i}$ follows a multivariate normal distribution
${\mathcal{N}}(0,\bm{\Sigma})$, and the linear model is defined as
$Y_{i}=X_{i}\beta+\epsilon_{i}$, where $\beta=(1,1,1,1,1,0,\cdots,0)^{T}$ and
$\epsilon_{i}\sim{\mathcal{N}}(0,1)$. We have a total of $n=100$ samples, and
among them there are $10$ outliers. The outliers are generated using
$\kappa_{i}=X_{i}\gamma$, where $\gamma=(0,0,0,0,0,1,1,\cdots)$. We did four
sets of simulations for various values of $p=200,400,800,1000$. To test the
hypothesis of whether the $i$-th observation is influential or not, we employ
a bootstrap procedure and utilize a threshold rule to determine whether an
individual is an outlier. We evaluate the performance of this outlier
identification procedure by comparing the receiver operating characteristic
(ROC) curves.

The ROC curves are depicted in Figure 1. In the figure, the curve labeled as
‘PCA-2’ represents the ROC curve generated by $\bar{{\bf X}}^{PCA}$ with a
dimensionality of $d=2$, while the curve labeled as ‘rSDR-0.2-2’ corresponds
to the curve produced by $\bar{{\bf X}}^{DR}$ with $\alpha=0.2$ and $d=2$.
Similarly, the remaining labels follow similar settings. It can be observed
that the curves generated by rSDR with $d=3$ consistently surpass those
produced by rSDR with $d=2$, and both outperform the curves generated by PCA.
This suggests that the proposed rSDR method effectively captures the
underlying structure of the data, and the resulting transformed data
$\bar{{\bf X}}^{DR}$ can be utilized for outlier detection. Notably, despite
the true subspace being two-dimensional, $\bar{{\bf X}}^{DR}$ with $d=3$
outperforms its two-dimensional counterpart. We speculate that the higher
dimensionality preserves more information due to the presence of outliers.

(a) $p=200$

(b) $p=400$

(c) $p=800$

(d) $p=1000$

Figure 1: ROC curves of outlier detection. The proposed robust SDR method with
$\alpha=0.5$ with projection dimension $3$ has the highest ROC in these four
simulation sample size settings.

###  4.3 Real Data Example: New Zealand Horse Mussels

A sample of 201 horse mussels (Modiolus modiolus) was collected at 5 sites in
the Marlborough Sounds at the Northeast of New Zealand’s South Island and this
dataset was discussed by Cook (2009). The response variable is muscle mass
$M$, the edible portion of the mussel, in grams. The quantitative predictors
are all related to characteristics of the mussel shells: shell width $W$ (in
mm), shell height $H$ (in mm), shell length $H$ (in mm) and shell mass $S$ (in
grams).

To process the data, a nonlinear transformation of the predictors was
recommended by Cook (2009) as $X=(L,W^{0.36},S^{0.11})$. Each column of the
data $X$ is further standardized by
$\tilde{X}=(\frac{L-\hat{\mu}_{L}}{\hat{\sigma}(L)},\frac{W^{0.36}-\hat{\mu}_{W^{0.36}}}{\hat{\sigma}(W^{0.36})},\frac{S^{0.11}-\hat{\mu}_{S^{0.11}}}{\hat{\sigma}(S^{0.11})})$
where $\hat{\mu}_{\cdot}$ is the sample mean and $\hat{\sigma}(\cdot)$ is the
sample standard deviation, since $L$ is on a larger scale than the other
predictors. Consequently, the predictors will have mean $0$ and variance $1$.
The rSDR model with $d=1$ would be appropriate to model this dataset, as shown
in Figure 2, where we fit two second-degree polynomial regression models of
the single index $\hat{\beta}^{T}\tilde{X}$ by rSDR with $\alpha=0.2$ and
$\alpha=1$. We compare our method rSDP with $\alpha=0.2$ and $\alpha=1$, SQP
and the Hilbert-Schmidt Independence Criterion (HSIC) method, proposed by N.
Zhang Yin (2015) for solving the special case of the SDR model, namely $d=1$.
Table 4.3 provides the estimated bases $\beta$ from these four methods. The
estimates of SDR with $\alpha=1$ and SQP are similar, and this result is
expected since SQP and rSDR with $\alpha=1$ solve the same model with
different algorithms. The estimated $\hat{\beta}$ from all four methods
indicate that the standardized shell mass predictor, $\tilde{X_{3}}$, is more
significant than the other two predictors while the rSDR with $\alpha=0.2$
produces a smaller value in the coefficient of $\tilde{X_{3}}$. However, rSDR
with $\alpha=0.2$ produces a model with a slightly larger R-squared value than
the other methods, which implies a better fit of the dataset.

(a) rSDR-0.2

(b) rSDR-1

Figure 2: The second-degree polynomial fitting of the single-index model in
the New Zealand Horse Mussels data using rSDR ($\alpha=0.2$) (left) and rSDR
($\alpha=1$) (right).

Estimated bases
$\hat{\beta}=[\hat{\beta}_{1},\hat{\beta}_{2},\hat{\beta}_{3}]^{T}\in{\mathbb{R}}^{3}$
of the central subspace in the New Zealand Horse Mussels data from various
methods and their adjusted R-squared values.  Method rSDR ($\alpha=0.2$) rSDR
($\alpha=1$) HSIC SQP $\hat{\beta}_{1}$ 0.2871 0.1832 0.1897 0.1831
$\hat{\beta}_{2}$ 0.0872 -0.0270 -0.0604 -0.0269 $\hat{\beta}_{3}$ 0.6391
0.8510 0.9800 0.8509 Adjusted R-squared 0.7026 0.6979 0.6962 0.6979

###  4.4 Real Data Example: Cardiomyopathy Microarray Data

The cardiomyopathy microarray dataset consists of 30 samples and 6319
predictors, originally used by Segal . (2003) to evaluate regression-based
approaches for microarray analysis. The focus of many researchers, Zou Yuan
(2008) and R. Li . (2012), has been to investigate the relationship between
the overexpression of a G protein-coupled receptor (Ro1) in mice and the 6319
associated genes. However, due to the high dimensionality of the data compared
to the limited number of samples, the sample covariance matrix is not
invertible. To address this issue, several methods have been proposed,
including SIS (Sure Independence Screening, Fan Lv (2008)), DCSIS (Distance
Correlation SIS, R. Li . (2012)), BCSIS (Ball Correlation SIS, Pan . (2019)),
and SDRLS (Sequential Dimension Reduction for Large $p$ Small $n$ problem, Yin
Hilafu (2015)). While SIS, DCSIS, and BCSIS are feature screening methods that
rank predictors based on a utility measure, they may not be robust against
outliers. Specifically, a set of predictors $\mathcal{A}=\\{i\mid
U(X_{i},Y)>\tau,i=1,\cdots,n\\}$ is determined for some threshold $\tau$ and
pre-selected utility measure $U$. SDRLS takes a different approach. SDRLS
partitions the data set into $X=[X_{1},X_{2}]$ with $\dim(X_{1})<n$ and
applies the SDR model on $(X_{1},[X_{2},Y])$ to obtain $R(X_{1})$. The
dimension of $R(X_{1})$ is chosen some integer that is smaller than
$\dim(X_{1})<n$ and thus a new predictor $[R(X_{1}),R_{2}]$ is obtained with a
smaller dimension. SDRLS iteratively repeats this process to achieve a
dimension smaller than the number of samples.

In this experiment, we utilized the SDRLS method to reduce the dimensionality
of the cardiomyopathy microarray data and assess the rSDR against heavy-tailed
predictors. The final dimension of the dataset was reduced to $p=19$, while
the dimension of the central subspace was set to $d=2$. The central subspace
is denoted as $\bm{\beta}=[\beta_{1},\beta_{2}]$. Indexes derived from this
reduction were obtained by projecting the processed cardiomyopathy microarray
dataset ${\bf X}$ onto the subspaces: ${\bf Z}_{1}=\beta_{1}^{T}{\bf X}$ and
${\bf Z}_{2}=\beta_{2}^{T}{\bf X}$. We performed linear and nonlinear
regression to model the response variable “Ro1” using predictors ${\bf Z}_{1}$
and ${\bf Z}_{2}$. In the nonlinear model, we introduced squared terms (${\bf
Z}_{1}^{2}$, ${\bf Z}_{2}^{2}$) and an interaction term (${\bf
Z}_{1}\times{\bf Z}_{2}$) in addition to the linear model. The regression
results are presented in Table 4.4. The findings demonstrate that our proposed
method, rSDR, with a smaller value of $\alpha$, outperforms the non-robust
version ($\alpha=1$) in both linear and nonlinear models.

Adjusted R-squared and F-value of models from SQP, rSDR in Cardiomyopathy
Microarray dataset.  Adjusted R-squared rSDR ($\alpha=0.2$) rSDR
($\alpha=0.5$) rSDR ($\alpha=1$) Linear 0.826 0.817 0.804 Nonlinear 0.882
0.871 0.867 F-value rSDR ($\alpha=0.2$) rSDR ($\alpha=0.5$) rSDR ($\alpha=1$)
Linear 70.1 65.8 60.6 Nonlinear 44.4 40.3 38.9

###  4.5 Real Data Example: Auto MPG data

We also employ the auto fuel economy data to illustrate the advantage of our
rSDR method. The dataset contains city-cycle fuel consumption in miles per
gallon (MPG) and 7 predictors: cylinders, displacement, horsepower, weight,
acceleration, model year and origin. As suggested in Sheng Yin (2016), we
avoid using “origin”, because it correlates with “cylinders” closely. Missing
values are deleted, and 392 observations are left for study. In order to
investigate the city-cycle fuel consumption in miles per gallon, we assume
that this data set fits a sufficient dimension reduction model. As shown in
Figure 3, there exist outliers in “horsepower” and “acceleration”. “cylinders”
and “displacement” that are not normally distributed. Therefore, rSDR is
appropriate for this data set.

Following the suggestion of Sheng Yin (2016), we use the dimension $d=2$ of
the central subspace. Let the subspace be $\bm{\beta}=[\beta_{1},\beta_{2}]$
and the auto MPG data be denoted as ${\bf X}$ where each column is centered
and scaled to make the variance as $1$. The following procedures are similar
to those done in the cardiomyopathy microarray data. The indexes are derived
by rSDR, after the linear and nonlinear regression models are constructed to
measure the goodness of fit of the two pair of indexes to “mpg”. In Table 4.5,
the adjusted R-squared and F-value of the linear model produced by rSDR with
$\alpha=0.2$ are larger than other non-robust models however it does not show
much superior in the nonlinear regression model (Figure 4).

Figure 3: Boxplots of predictors in the auto mpg data set. Some predictors
such as “horsepower” and “acceleration” have outlying observations.

(a) SQP-index-1

(b) SQP-index-2

(c) rSDR-index-1

(d) rSDR-index-2

Figure 4: Scatter plots of “mpg” versus the indexes produced by SQP and rSDR
($\alpha=0.2$).

Adjusted R-squared and F-value of models from SQP, rSDR in MPG dataset.
Adjusted R-squared rSDR ($\alpha=0.2$) rSDR ($\alpha=0.5$) rSDR ($\alpha=1$)
Linear 0.807 0.806 0.804 Nonlinear 0.850 0.853 0.845 F-value rSDR
($\alpha=0.2$) rSDR ($\alpha=0.5$) rSDR ($\alpha=1$) Linear 817 816 807
Nonlinear 444 456 427

##  5 Discussion

In this article, the proposed rSDR using $\alpha$-dCov is robust against
outliers in both the response and predictors. Further, the proposed manifold-
learning estimation method is less sensitive to the choice of the initial
estimators. Both simulation and real-world data applications show that the
proposed method outperforms the existing methods. The proposed method does not
suffer from multicollinearity which could impact the performance of the
traditional SDR methods in high-dimensional data analysis. Simulation and
real-world data studies show its advantages in terms of computational
efficiency and robustness against outliers.

## Acknowledgments

The authors would like to thank the Editor, the Associate Editor and the
reviewers for their constructive and insightful comments that greatly improved
the manuscript. This work was partially supported by NSF grants (DMS-1924792,
DMS-2318925 and CNS-1818500).

## Author contribution

Hsin-Hsiung Huang: conceptualization, methodology, formal analysis,
investigation, writing-original draft preparation, writing-review,
supervision. Feng Yu: conceptualization, methodology, formal analysis,
investigation, writing review, programming and numerical results. Teng Zhang:
conceptualization, methodology, formal analysis, investigation, writing-
review, supervision. All authors have read and agreed to the published version
of the manuscript.

## Declarations

The authors declare no competing interests.

## References

  * P. Absil . (2009) absil2009optimizationAbsil, P., Mahony, R. Sepulchre, R.  2009\.  Optimization Algorithms on Matrix Manifolds Optimization algorithms on matrix manifolds.  Princeton University Press.  https://books.google.com/books?id=NSQGQeLN3NcC
  * PA. Absil Malick (2012) absil2012projectionAbsil, PA. Malick, J.  2012\.  Projection-like retractions on matrix manifolds Projection-like retractions on matrix manifolds.  SIAM Journal on Optimization221135–158.
  * Cook (2009) cook2009regressionCook, RD.  2009\.  Regression graphics: Ideas for studying regressions through graphics Regression graphics: Ideas for studying regressions through graphics.  John Wiley & Sons.
  * Cook Weisberg (1991) cook1991discussionCook, RD. Weisberg, S.  1991\.  Discussion of sliced inverse regression for dimension reduction Discussion of sliced inverse regression for dimension reduction.  Journal of the American Statistical Association86414328–332.
  * Dalmau-Cedeno Oviedo (2017) dalmau2017projectionDalmau-Cedeno, O. Oviedo, H.  2017\.  A projection method for optimization problems on the stiefel manifold A projection method for optimization problems on the stiefel manifold.  Mexican conference on pattern recognition Mexican conference on pattern recognition ( 84–93).
  * Fan Lv (2008) fan2008sureFan, J. Lv, J.  2008\.  Sure independence screening for ultrahigh dimensional feature space Sure independence screening for ultrahigh dimensional feature space.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)705849–911.
  * Hall Li (1993) hall1993almostHall, P. Li, KC.  1993\.  On almost linearity of low dimensional projections from high dimensional data On almost linearity of low dimensional projections from high dimensional data.  The annals of Statistics867–889.
  * B. Li Wang (2007) li2007directionalLi, B. Wang, S.  2007\.  On directional regression for dimension reduction On directional regression for dimension reduction.  Journal of the American Statistical Association102479997–1008.
  * KC. Li (19911) Li1991Li, KC.  19911\.  Sliced Inverse Regression for Dimension Reduction Sliced inverse regression for dimension reduction.  Journal of the American Statistical Association86414316-327.  https://www.tandfonline.com/doi/abs/10.1080/01621459.1991.10475035  10.1080/01621459.1991.10475035
  * KC. Li (19912) li1991slicedLi, KC.  19912\.  Sliced inverse regression for dimension reduction Sliced inverse regression for dimension reduction.  Journal of the American Statistical Association86414316–327.
  * R. Li . (2012) li2012featureLi, R., Zhong, W. Zhu, L.  2012\.  Feature screening via distance correlation learning Feature screening via distance correlation learning.  Journal of the American Statistical Association1074991129–1139.
  * Pan . (2019) pan2019genericPan, W., Wang, X., Xiao, W. Zhu, H.  2019\.  A generic sure independence screening procedure A generic sure independence screening procedure.  Journal of the American Statistical Association114526928–937.
  * Segal . (2003) segal2003regressionSegal, MR., Dahlquist, KD. Conklin, BR.  2003\.  Regression approaches for microarray data analysis Regression approaches for microarray data analysis.  Journal of Computational Biology106961–980.
  * Sheng Yin (2013) sheng2013directionSheng, W. Yin, X.  2013\.  Direction estimation in single-index models via distance covariance Direction estimation in single-index models via distance covariance.  Journal of Multivariate Analysis122148–161.
  * Sheng Yin (2016) sheng2016sufficientSheng, W. Yin, X.  2016\.  Sufficient dimension reduction via distance covariance Sufficient dimension reduction via distance covariance.  Journal of Computational and Graphical Statistics25191–104.
  * Székely Rizzo (2009) szekely2009brownianSzékely, GJ. Rizzo, ML.  2009\.  Brownian distance covariance Brownian distance covariance.  The annals of applied statistics341236–1265.
  * Székely Rizzo (2014) 10.1214/14-AOS1255Székely, GJ. Rizzo, ML.  2014\.  Partial distance correlation with methods for dissimilarities Partial distance correlation with methods for dissimilarities.  The Annals of Statistics4262382 – 2412.  https://doi.org/10.1214/14-AOS1255  10.1214/14-AOS1255
  * Székely . (2007) 10.1214/009053607000000505Székely, GJ., Rizzo, ML. Bakirov, NK.  2007\.  Measuring and testing dependence by correlation of distances Measuring and testing dependence by correlation of distances.  The Annals of Statistics3562769 – 2794.  https://doi.org/10.1214/009053607000000505  10.1214/009053607000000505
  * Wang Li (2017) wang2017outlierWang, T. Li, Z.  2017\.  Outlier detection in high-dimensional regression model Outlier detection in high-dimensional regression model.  Communications in Statistics-Theory and Methods46146947–6958.
  * Wen Yin (2013) wen2012Wen, Z. Yin, W.  2013\.  A feasible method for optimization with orthogonality constraints A feasible method for optimization with orthogonality constraints.  Mathematical Programming1421397–434.  https://doi.org/10.1007/s10107-012-0584-1  10.1007/s10107-012-0584-1
  * Wold . (1987) wold1987principalWold, S., Esbensen, K. Geladi, P.  1987\.  Principal component analysis Principal component analysis.  Chemometrics and intelligent laboratory systems21-337–52.
  * Wu Chen (2021) wu2021mmWu, R. Chen, X.  2021\.  MM algorithms for distance covariance based sufficient dimension reduction and sufficient variable selection Mm algorithms for distance covariance based sufficient dimension reduction and sufficient variable selection.  Computational Statistics & Data Analysis155107089.
  * Xia (2008) xia2008multipleXia, Y.  2008\.  A multiple-index model and dimension reduction A multiple-index model and dimension reduction.  Journal of the American Statistical Association1034841631–1640.
  * Xia . (2002) xia2002adaptiveXia, Y., Tong, H., Li, WK. Zhu, LX.  2002\.  An adaptive estimation of dimension reduction space An adaptive estimation of dimension reduction space.  Journal of the Royal Statistical Society Series B: Statistical Methodology643363–410.
  * Yin Cook (2005) yin2005directionYin, X. Cook, RD.  2005\.  Direction estimation in single-index regressions Direction estimation in single-index regressions.  Biometrika922371–384.
  * Yin Hilafu (2015) yin2015sequentialYin, X. Hilafu, H.  2015\.  Sequential sufficient dimension reduction for large p, small n problems Sequential sufficient dimension reduction for large p, small n problems.  Journal of the Royal Statistical Society: Series B: Statistical Methodology879–892.
  * Yin Li (2011) Yin11-AOS950Yin, X. Li, B.  2011\.  Sufficient dimension reduction based on an ensemble of minimum average variance estimators Sufficient dimension reduction based on an ensemble of minimum average variance estimators.  The Annals of Statistics3963392 – 3416.  https://doi.org/10.1214/11-AOS950  10.1214/11-AOS950
  * Yin . (2008) yin2008successiveYin, X., Li, B. Cook, RD.  2008\.  Successive direction extraction for estimating the central subspace in a multiple-index regression Successive direction extraction for estimating the central subspace in a multiple-index regression.  Journal of Multivariate Analysis9981733–1757.
  * J. Zhang Chen (2019) zhang2019robustZhang, J. Chen, X.  2019\.  Robust sufficient dimension reduction via ball covariance Robust sufficient dimension reduction via ball covariance.  Computational Statistics & Data Analysis140144–154.
  * N. Zhang Yin (2015) zhang2015directionZhang, N. Yin, X.  2015\.  Direction estimation in single-index regressions via Hilbert-Schmidt independence criterion Direction estimation in single-index regressions via hilbert-schmidt independence criterion.  Statistica Sinica743–758.
  * Zhu . (2019) zhu2019orthodrZhu, R., Zhang, J., Zhao, R., Xu, P., Zhou, W. Zhang, X.  2019\.  orthoDr: Semiparametric Dimension Reduction via Orthogonality Constrained Optimization orthodr: Semiparametric dimension reduction via orthogonality constrained optimization.  The R JournalVol. 11/2.
  * Zou Yuan (2008) zou2008regularizedZou, H. Yuan, M.  2008\.  Regularized simultaneous model selection in multiple quantiles regression Regularized simultaneous model selection in multiple quantiles regression.  Computational Statistics & Data Analysis52125296–5304.