# Computationally Efficient and Statistically Optimal Robust High-Dimensional
Linear Regression

Yinan Shen, Jingyang Li, Jian-Feng Cai111Jian-Feng Cai’s research was
partially supported by Hong Kong RGC Grant GRF 16310620 and GRF 16309219. and
Dong Xia222Dong Xia’s research was partially supported by Hong Kong RGC Grant
ECS 26302019, GRF 16303320 and GRF 16300121.
Hong Kong University of Science and Technology

(())

###### Abstract

High-dimensional linear regression under heavy-tailed noise or outlier
corruption is challenging, both computationally and statistically. Convex
approaches have been proven statistically optimal but suffer from high
computational costs, especially since the robust loss functions are usually
non-smooth. More recently, computationally fast non-convex approaches via sub-
gradient descent are proposed, which, unfortunately, fail to deliver a
statistically consistent estimator even under sub-Gaussian noise. In this
paper, we introduce a projected sub-gradient descent algorithm for both the
sparse linear regression and low-rank linear regression problems. The
algorithm is not only computationally efficient with linear convergence but
also statistically optimal, be the noise Gaussian or heavy-tailed with a
finite $1+\varepsilon$ moment. The convergence theory is established for a
general framework and its specific applications to absolute loss, Huber loss
and quantile loss are investigated. Compared with existing non-convex methods,
ours reveals a surprising phenomenon of two-phase convergence. In phase one,
the algorithm behaves as in typical non-smooth optimization that requires
gradually decaying stepsizes. However, phase one only delivers a statistically
sub-optimal estimator, which is already observed in the existing literature.
Interestingly, during phase two, the algorithm converges linearly as if
minimizing a smooth and strongly convex objective function, and thus a
constant stepsize suffices. Underlying the phase-two convergence is the
smoothing effect of random noise to the non-smooth robust losses in an area
close but not too close to the truth. Numerical simulations confirm our
theoretical discovery and showcase the superiority of our algorithm over prior
methods.

##  1 Introduction

Let $\\{({\mathbf{X}}_{i},Y_{i})\\}_{i=1}^{n}$ be a collection of independent
observations satisfying

$Y_{i}=\langle{\mathbf{X}}_{i},{\mathbf{T}}^{\ast}\rangle+\xi_{i}$ (1)

with ${\mathbf{X}}_{i}$ and $Y_{i}$ being the covariate and response,
respectively. Here, $\langle\cdot,\cdot\rangle$ denotes the inner product in
Euclidean space. The latent noise $\xi_{i}$ is usually assumed random for ease
of analysis. We assume that ${\mathbf{T}}^{\ast}$ resides in a high-
dimensional ambient space but has a low-dimensional structure. We focus on two
high-dimensional linear regression problems: the sparse linear regression
where ${\mathbf{T}}^{\ast}\in{\mathbb{R}}^{d}$ is sparse and low-rank linear
regression where ${\mathbf{T}}^{\ast}\in{\mathbb{R}}^{d\times d}$ is low-rank.
For both problems, it is of great interest to effectively, both
computationally and statistically, recover ${\mathbf{T}}^{\ast}$ using as few
observations as possible. There is a vast literature discussing the
motivations and applications, e.g., Tibshirani (1996); Needell and Tropp
(2009); Blumensath and Davies (2009, 2010); Loh and Wainwright (2011); Shen
and Li (2017); Zhu et al. (2020); Zhao et al. (2022); Gross et al. (2010); Xia
and Koltchinskii (2016); Chen et al. (2011); Negahban and Wainwright (2011);
Chiu et al. (2021); Siddiqi et al. (2010).

The methods of high-dimensional linear regression can be broadly classified
into two categories: the convex and non-convex approaches. The typical convex
approaches introduce an additional penalization such as the matrix nuclear
norm to promote low-rank solutions and the $\ell_{1}$-norm to promote
sparsity. The convexity renders well-understood algorithms for convex
programming immediately applicable. Therefore statisticians can just focus on
studying the statistical performances without paying too much attention to the
computational implementations. For instance, by assuming the i.i.d.
(sub-)Gaussian noise with a variance $\sigma^{2}$ and the so-called restricted
isometry property (RIP) or restricted strong convexity (RSC) in low-rank
linear regression for estimating a $d\times d$ matrix, Candes and Plan (2011);
Rohde and Tsybakov (2011); Negahban and Wainwright (2011); Cai and Zhang
(2015); Chen and Wainwright (2015); Davenport and Romberg (2016);
Chandrasekaran et al. (2012) proved that the nuclear-norm penalized least
square estimator attains the error rate
$\tilde{O}_{p}(\sigma^{2}rdn^{-1})$333Here $\tilde{O}_{p}(\cdot)$ stands for
the typical big-O notation up to logarithmic factors and holds with high
probability. in the squared Frobenius norm where $r$ is the rank. This rate is
shown to be optimal in the minimax sense. See, e.g., Xia (2014); Ma and Wu
(2015) and references therein. Despite the appealing theoretical performances,
convex approaches have two major drawbacks. First, convex methods operate
directly in the high-dimensional ambient space making them run slowly and are
unscalable to ultra-high dimensional problems. Secondly, while the established
error rate is theoretically minimax optimal, the implicit constant factor
seems large and the resultant error rate is, in practice, often inferior to
that by non-convex approaches. These issues can be nicely addressed by non-
convex approaches. In low-rank linear regression, the non-convex methods could
operate on the basis of matrix factorization and resort to the gradient
descent type algorithms (Chen and Wainwright, 2015; Burer and Monteiro, 2003;
Zheng and Lafferty, 2016; Zhao et al., 2015; Xia and Yuan, 2021; Wei et al.,
2016); in sparse linear regression, the typical non-convex approaches include
the least angel regression (Efron et al., 2004), iterative hard thresholding
(Blumensath and Davies, 2009), and best subset selection (Zhu et al., 2020),
etc. Under similar conditions, it has been demonstrated that these non-convex
algorithms converge fast (more exactly, linearly) and also deliver minimax
optimal estimators under (sub-)Gaussian noise. We remark that non-convex
methods often require more tuning parameters. For instance, the matrix
factorization algorithms often require knowing the true rank (Chandrasekaran
et al., 2012; Chen and Wainwright, 2015) and the iterative hard thresholding
algorithms need the information of the underlying sparsity (Blumensath and
Davies, 2009; Efron et al., 2004).

The recent boom of data technology poses new challenges to high-dimensional
linear regression, among which the heavy-tailed noise and outliers appear
routinely in numerous applications such as diffusion-weighted imaging (Chang
et al., 2005), on-line advertising (Sun et al., 2017), and gene-expression
data analysis (Sun et al., 2020). The aforementioned approaches minimizing a
square loss become vulnerable or even completely useless when the noise has a
heavy tail or when data are partially corrupted. To mitigate the influence of
heavy-tailed noise and corruptions, a natural solution is to replace the
square loss by more robust ones. Notable examples include the absolute loss
(Candès et al., 2011; Cambier and Absil, 2016), the renowned Huber loss
(Huber, 1965), and quantile loss (Koenker and Hallock, 2001), all of which are
convex but non-smooth. Elsener and van de Geer (2018) proposed a convex
approach based on nuclear-norm penalized absolute or Huber loss for low-rank
linear regression, and proved that the estimator attains the error rate
$\tilde{O}_{p}(rd_{1}n^{-1})$ in squared Frobenius norm as long as the noise
has a non-zero density in a neighbourhood of origin. Interestingly, Elsener
and van de Geer (2018) doesn’t need to estimate the noise level for choosing
the tuning parameters. A more general framework requiring only Lipschitz and
convex loss but imposing a so-called Bernstein condition on the noise is
investigated by Alquier et al. (2019). See also Klopp et al. (2017). The
foregoing approaches still rely on convex programming and suffer from the
issue of computational inefficiency explained above. In fact, the
computational issue is severer here with non-smooth objective functions (Boyd
et al., 2004).

The sub-gradient descent algorithms based on matrix factorization were studied
recently by Charisopoulos et al. (2021); Li et al. (2020); Tong et al. (2021b)
in low-rank linear regression. These algorithms converge linearly if equipped
with a good initialization and with properly chosen stepsizes, but the
statistical performances of these non-convex approaches under heavy-tailed
noise are either largely missing or, generally, sub-optimal. To be more
specific, Li et al. (2020) only proved the exact recovery under sparse
outliers without noise; the error rates in squared Frobenius norm established
in Charisopoulos et al. (2021) and Tong et al. (2021b) turn into
$\tilde{O}_{p}(\sigma^{2})$ even when the noise is Gaussian with variance
$\sigma^{2}$, which is not statistically consistent. There is also a vast
literature of sparse linear regression under heavy-tailed noise or partial
corruption, e.g., Sun et al. (2020); Zheng et al. (2015); Belloni et al.
(2011); Alquier et al. (2019); Lecué and Mendelson (2017); Lugosi and
Mendelson (2019b); Pan et al. (2021); Wang et al. (2021, 2020); Thompson
(2020). The adpative Huber loss was investigated by Sun et al. (2020) under a
finite $1+\varepsilon$ moment condition of noise. Depersin (2020) studied the
median-of-means approach and derived sharp error rates by VC-dimension,
achieving the first sub-Gaussian sparse estimators. Belloni et al. (2011) and
Wang et al. (2020) proposed the square-root LASSO and the absolute loss,
respectively. While being statistically robust to heavy tails, these prior
works didn’t provide provably fast converging algorithms. In Foucart and Lecué
(2017), an iterative hard thresholding algorithm was proposed for sparse
linear regression by minimizing the absolute loss, which converges fast but a
sub-optimal error rate was derived.

In this paper, we propose a computationally efficient non-convex algorithm
based on the projected sub-gradient descent for high-dimensional linear
regression using robust loss functions. It is essentially the absolute loss
hard thresholding algorithm (IHT-$\ell_{1}$,Foucart and Lecué (2017)) for
sparse linear regression and the Riemannian sub-gradient descent algorithm
(RsGrad, Vandereycken (2013); Cambier and Absil (2016); Wei et al. (2016); Cai
et al. (2021a)) for low-rank linear regression. In both problems, the sub-
gradient descent algorithm admits fast computation and linear convergence.
Unlike Foucart and Lecué (2017) and Tong et al. (2021b), we discover that the
sub-gradient descent algorithm exhibits an intriguing phenomenon referred to
as the two-phase convergence. In phase one, the algorithm behaves like typical
non-smooth optimization requiring gradually decaying stepsizes through
iterations. For example, in low-rank linear regression, if the noise is
Gaussian with a variance $\sigma^{2}$ and the absolute loss is equipped, the
phase-one iterations converge linearly and reach an estimator whose squared
Frobenius norm error is $O_{p}(\sigma^{2})$. This rate has been achieved by
Tong et al. (2021b); Charisopoulos et al. (2021) whereas theirs is based on a
scaled sub-gradient (ScaledSM) descent algorithm. In phase two, the algorithm
(RsGrad) behaves like smooth optimization in that a constant stepsize
guarantees a linear convergence. As a result, after phase-two iterations,
RsGrad can output an estimator which is statistically optimal under mild
conditions, e.g. it achieves the rate $O_{p}(\sigma^{2}rdn^{-1})$ for low-rank
linear regression. The two-phase convergence phenomenon also occurs for sparse
linear regression. For both problems, minimax optimal rates can be achieved
under heavy-tailed noise.

Our contributions are multi-fold. First, we propose a projected sub-gradient
descent algorithm for robust high-dimensional linear regression that is
applicable to a wide class of non-smooth loss functions for both the sparse
linear regression and low-rank linear regression problems. While the
Riemannian sub-gradient descent (RsGrad) has been introduced by Cambier and
Absil (2016) for minimizing the absolute loss, to our best knowledge, there
exists no theoretical guarantees for its convergence. Under mild conditions,
we prove that the algorithm converges linearly. Secondly, we demonstrate the
statistical optimality of the final estimator delivered by our algorithm for
both problems, be the noise Gaussian or heavy-tailed. For instance, the
specific applications of RsGrad to low-rank linear regression with the
absolute loss, Huber loss and quantile loss confirm that a rate
$O_{p}(rdn^{-1})$ is attainable as long as the noise has a fine
$1+\varepsilon$ moment and its density satisfies mild regularity conditions.
Similar conditions have appeared in Elsener and van de Geer (2018). Unlike Fan
et al. (2021) and Minsker (2018), our rate is proportional to the noise size.
Thirdly, our analysis reveals a new phenomenon of two-phase convergence that
enables us to achieve the statistically optimal error rates, which is a
significant improvement over existing works (Foucart and Lecué, 2017; Li et
al., 2020; Tong et al., 2021b; Charisopoulos et al., 2021). While the phase-
one convergence is typical for non-smooth optimization, the phase-two
convergence, interestingly, behaves like smooth optimization. It seems that
the random noise has an effect of smoothing in the phase-two convergence.

The rest of the paper is organized as follows. In Section 2, we review three
typical examples of non-smooth robust loss functions. We provide a warm-up
example of robust sparse linear regression in Section 3 and investigate the
two-phase convergence of the iterative hard thresholding algorithm. A general
framework of studying the two-phase convergence performance of the Riemannian
sub-gradient descent algorithm for robust low-rank linear regression is
presented in Section 4. Non-convex algorithms inevitably involve important
tuning parameters, which will be discussed in Section 5. In Section 6, we
further investigate our methods under the Huber’s contamination model. We
showcase the results of numerical experiments and comparison with prior
methods in Section 7. All the proofs are relegated to the Appendix.

##  2 Robust Loss Functions

Denote ${\mathbb{M}}$ the feasible set of the optimization program. It can be
the set of matrices with rank bounded by $r$ in low-rank regression (Section
4) and the set of vectors with support size bounded by $\tilde{s}$ in sparse
regression (Section 3). The goal of linear regression is to solve

$\displaystyle\hat{\mathbf{T}}:=\underset{{\mathbf{T}}\in{\mathbb{M}}}{\arg\min}\
f({\mathbf{T}})\ \textrm{ where
}f({\mathbf{T}}):=\sum_{i=1}^{n}\rho\big{(}\langle{\mathbf{T}},{\mathbf{X}}_{i}\rangle-
Y_{i}\big{)}.$ (2)

Here $\rho(\cdot):{\mathbb{R}}\mapsto{\mathbb{R}}_{+}$ is a properly chosen
loss function. The statistical property of $\hat{\mathbf{T}}$ crucially relies
on the loss function. For instance, $\hat{\mathbf{T}}$ attained by using the
$\ell_{2}$-loss (square loss), i.e., $\rho(x):=x^{2}$, has been proven
effective (Chen and Wainwright, 2015; Wei et al., 2016; Xia and Yuan, 2021;
Jain et al., 2014; Agarwal et al., 2010) in dealing with sub-Gaussian noise
but is well-recognized sensitive to outliers and heavy-tailed noise (Huber,
1965). There are the so-called robust loss functions which are relatively
immune to outliers and heavy-tailed noise and capable to deliver a more
reliable estimate $\hat{\mathbf{T}}$. The examples of robust loss functions
include:

  1. 1.

Absolute loss ($\ell_{1}$-loss): $\rho(x):=|x|$ for any $x\in{\mathbb{R}}$;

  2. 2.

Huber loss:
$\rho_{H,\delta}(x):=x^{2}\mathbbm{1}(|x|\leq\delta)+(2\delta|x|-\delta^{2})\mathbbm{1}(|x|>\delta)$
for any $x\in{\mathbb{R}}$ where $\delta>0$ is a tuning parameter;

  3. 3.

Quantile loss: $\rho_{Q,\delta}(x):=\delta x\mathbbm{1}(x\geq
0)+(\delta-1)x\mathbbm{1}(x<0)$ for any $x\in{\mathbb{R}}$ with
$\delta:={\mathbb{P}}(\xi\leq 0)$.

The absolute loss, Huber loss and quantile loss are all convex functions. They
have appeared in the literature of sparse regression including the LASSO-type
convex approaches (Sun et al., 2020; Fan et al., 2021) and the non-convex
approaches by iterative algorithms (Foucart and Lecué, 2017). These robust
losses have also been investigated for low-rank regression, e.g., the convex
approach based on nuclear norm penalization (Elsener and van de Geer, 2018;
Klopp et al., 2017; Candès et al., 2011) and the non-convex approach based on
gradient-descent-type algorithms (Li et al., 2020; Tong et al., 2021a).
Compared with the square loss, the aforementioned robust loss functions are
non-smooth, i.e., their derivatives are dis-continuous, which brings new
challenge to its computation. Indeed, the optimizing program (2) usually
exploits the sub-gradient of robust loss functions, which is written as
$\partial f({\mathbf{T}})$. See, for instance, Charisopoulos et al. (2021) and
references therein. Recently, Thompson (2020) showed that regularization using
the sorted Huber loss could significantly improve the constant in the error
rate than that using the standard Huber loss. However, since the computation
of sorted Huber loss is challenging in our framework, we lave it for future
work.

The existing statistics literature have proved that $\hat{\mathbf{T}}$
attained by using robust loss functions enjoys the statistical robustness to
heavy-tailed noise. However, the optimal solution $\hat{\mathbf{T}}$ is often
solved by heuristic algorithms without convergence guarantees (Sun et al.,
2020) or by non-smooth optimization algorithms (Elsener and van de Geer, 2018)
which may converge slowly in general and the computational cost can be
considerably high, rendering the applicability of aforementioned robust loss
functions questionable for real-world high-dimensional problems. More
recently, Charisopoulos et al. (2021); Tong et al. (2021b); Foucart and Lecué
(2017) investigated sub-gradient descent algorithms for minimizing the robust
loss functions showing that these algorithms can converge fast if the
stepsizes are properly scheduled along the iterations. Unfortunately, their
analysis frameworks all failed to demonstrate that the statistical benefits of
robust loss functions claimed by prior statistics literature can be achieved
by the sub-gradient descent algorithms. There exists a gap between the
statistical and computational efficiencies.

In the subsequent Sections 3 and 4, we will show that a revised sub-gradient
descent algorithm equipped with a novel stepsize schedule converges fast and
achieves statistically optimal rates in various settings. Besides the
algorithmic design, a novel analysis framework also plays a critical role in
developing our statistical theories.

##  3 Warm-up: Sparse Linear Regression

We begin with a warm-up example of high-dimensional sparse linear regression
using the absolute loss. Let ${\boldsymbol{\beta}}^{*}\in{\mathbb{R}}^{d}$ be
sparse with a support size $|\text{supp}({\boldsymbol{\beta}}^{*})|=s\ll d$.
We observe a collection of independent random pairs
$\\{({\mathbf{X}}_{i},Y_{i})\\}_{i=1}^{n}$ satisfying
$Y_{i}=\langle{\mathbf{X}}_{i},{\boldsymbol{\beta}}^{*}\rangle+\xi_{i}$ where
the predictor vector ${\mathbf{X}}_{i}\in{\mathbb{R}}^{d}$ is assumed randomly
sampled. Our goal is to estimate ${\boldsymbol{\beta}}^{\ast}$ when the noise
$\xi_{i}$ may have heavy tails. For ease of exposition, we only study the
absolute loss (Foucart and Lecué, 2017) but parallel results can be derived
for other robust loss functions (see Section 4). Denote

$\displaystyle\hat{{\boldsymbol{\beta}}}:=\underset{|\text{supp}({\boldsymbol{\beta}})|\leq\tilde{s}}{\arg\min}f({\boldsymbol{\beta}}),\quad\textrm{where
}f({\boldsymbol{\beta}}):=\|{\mathbf{X}}{\boldsymbol{\beta}}-{\mathbf{Y}}\|_{1}:=\sum_{i=1}^{n}|Y_{i}-\langle{\boldsymbol{\beta}},{\mathbf{X}}_{i}\rangle|.$
(3)

Here the chosen sparsity level $\tilde{s}\in[s,d]$,
${\mathbf{X}}=({\mathbf{X}}_{1},\dots,{\mathbf{X}}_{n})^{\top}$ and
${\mathbf{Y}}=(Y_{1},\dots,Y_{n})^{\top}$. Let $\|\cdot\|_{p}$ denote the
$\ell_{p}$-norm of vectors for $1\leq p\leq 2$.

###  3.1 Algorithm: Iterative Hard Thresholding

Our method for optimizing the program (3) is essentially the projected sub-
gradient descent algorithm. The projected gradient descent algorithm has been
explored in minimizing the square loss (Agarwal et al., 2010; Chen and
Wainwright, 2015), which is easy to implement and fast computable. We resort
to the sub-gradient because the absolute loss is non-smooth.

Since projecting a vector onto the set of sparse vectors can be realized by
hard thresholding, our algorithm is often referred to as the Iterative Hard
Thresholding (Blumensath and Davies, 2009, 2010). The pseudocodes are
presented in Algorithm 1. The operator
${\mathcal{H}}_{\tilde{s}}({\mathbf{V}})$ keeps ${\mathbf{V}}$’s $\tilde{s}$
entries with the largest absolute values and zeros out the others. In
practice, we choose the sub-gradient
${\mathbf{G}}_{l}=\sum_{i=1}^{n}sign(\langle{\boldsymbol{\beta}}_{l},{\mathbf{X}}_{i}\rangle-
Y_{i}){\mathbf{X}}_{i}$ and set $sign(0)=0$. Note that Algorithm 1 allows
arbitrary initializations.

Algorithm 1 Iterative Hard Thresholding ($\ell_{1}$-loss)

Input: observations $\\{({\mathbf{X}}_{i},Y_{i})\\}_{i=1}^{n}$, max iterations
${l_{\textsf{max}}}$, step sizes $\\{\eta_{l}\\}_{l=0}^{{l_{\textsf{max}}}}$
and sparsity parameter $\tilde{s}$.

Initialization: arbitrary ${\boldsymbol{\beta}}_{0}$ with a support size
$|\text{supp}({\boldsymbol{\beta}}_{0})|\leq\tilde{s}$

for $l=0,\ldots,{l_{\textsf{max}}}$ do

Choose a sub-gradient: ${\mathbf{G}}_{l}\in\partial
f({\boldsymbol{\beta}}_{l})$

Truncation:
${\boldsymbol{\beta}}_{l+1}={\mathcal{H}}_{\tilde{s}}({\boldsymbol{\beta}}_{l}-\eta_{l}{\mathbf{G}}_{l})$

end for

Output: $\hat{\boldsymbol{\beta}}={\boldsymbol{\beta}}_{{l_{\textsf{max}}}}$

###  3.2 Theory: two-phase convergence

The convergence performance of Algorithm 1 and statistical property of
$\hat{\boldsymbol{\beta}}$ rely on the following two assumptions.

######  Assumption 1.

The covariate vectors ${\mathbf{X}}_{1},\dots,{\mathbf{X}}_{n}$ are
independent Gaussian ${\mathbf{X}}_{i}\sim
N(\boldsymbol{0},{\boldsymbol{\Sigma}}_{i})$ where
${\boldsymbol{\Sigma}}_{i}$’s are symmetric and positive definite. There exist
absolute constants $C_{{\tiny l}},C_{{\tiny u}}>0$ such that

$\displaystyle C_{{\tiny
l}}\leq\lambda_{\min}({\boldsymbol{\Sigma}}_{i})\leq\lambda_{\max}({\boldsymbol{\Sigma}}_{i})\leq
C_{{\tiny u}}.$

######  Assumption 2.

(Heavy-tailed noise 1) The noise $\xi_{1},\cdots,\xi_{n}$ are i.i.d. with the
density $h_{\xi}(\cdot)$ and distribution function $H_{\xi}(\cdot)$,
respectively. There exists an $\varepsilon>0$ such that
${\mathbb{E}}|\xi|^{1+\varepsilon}<+\infty$. The noise has median zero, i.e.,
$H_{\xi}(0)=1/2$. Denote $\gamma={\mathbb{E}}|\xi|$. There exist constants
$b_{0},b_{1}>0$ (may be dependent on $\gamma$) such that 444The lower bound
condition can be slightly relaxed to $|H_{\xi}(x)-H_{\xi}(0)|\geq|x|/b_{0}$
when $|x|\leq 8(C_{{\tiny u}}/C_{{\tiny l}})^{1/2}\gamma$.

$\displaystyle h_{\xi}(x)\geq b_{0}^{-1},$ $\displaystyle\ \ \ \textrm{ for
all }|x|\leq 8(C_{{\tiny u}}/C_{{\tiny l}})^{1/2}\gamma;$ $\displaystyle
h_{\xi}(x)\leq b_{1}^{-1},$ $\displaystyle\ \ \ \forall x\in{\mathbb{R}}.$

By Assumption 2, a simple fact is $b_{0}\geq 8(C_{{\tiny u}}/C_{{\tiny
l}})^{1/2}\gamma$ and $b_{0}\geq b_{1}$. We require the noise has a finite
$1+\varepsilon$ moment, which is essential to guarantee
$\lim_{x\to+\infty}x(1-H_{\xi}(x))=0$. This moment condition is fairly weak
compared with the existing literature. For instance, Minsker (2018); Depersin
(2020) require a finite second-order moment condition; Fan et al. (2021)
imposes a $2+\varepsilon$ moment condition. The lower bound on the density
function also appeared in Elsener and van de Geer (2018) and is a special case
of the (local) Bernstein condition (Alquier et al., 2019; Chinot et al.,
2020). The upper bound condition on the density function is mild, e.g., a
Lipschitz distribution function ensures such a uniform upper bound.

######  Theorem 1.

Suppose Assumptions 1 and 2 hold. There exist absolute constants
$C_{1},\dots,C_{6}>0$, $c_{2}^{*}\in(0,1)$ such that if $n\geq C_{1}C_{{\tiny
u}}C_{{\tiny l}}^{-1}\tilde{s}\log(2d/\tilde{s})$, $\tilde{s}\geq
C_{2}(C_{{\tiny u}}/C_{{\tiny l}})^{8}(b_{0}/b_{1})^{8}\cdot s$ and the
initial stepsize $\eta_{0}\in(nC_{{\tiny u}})^{-1}C_{{\tiny
l}}^{1/2}\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{*}\|_{2}\cdot\left[1/8,\
3/8\right]$, then with probability over
$1-\exp(-C_{3}\tilde{s}\log(2d/\tilde{s}))-3\exp(-(n\log(2d/\tilde{s}))^{1/2}\log^{-1}n)$,
the sequence $\\{{\boldsymbol{\beta}}_{l}\\}_{l\leq{l_{\textsf{max}}}}$ of
Algorithm 1 has the following dynamics:

  1. (1)

in phase one, i.e., when
$\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{*}\|_{2}\geq 8C_{{\tiny
l}}^{-1/2}\gamma$, by taking the stepsize $\eta_{l}=(1-c_{1})^{l}\eta_{0}$
with any fixed constant $c_{1}\leq C_{{\tiny l}}C_{{\tiny u}}^{-1}/64$, then
$\|{\boldsymbol{\beta}}_{l+1}-{\boldsymbol{\beta}}^{*}\|_{2}\leq(1-c_{1})^{l+1}\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{*}\|_{2}$;

  2. (2)

in phase two, i.e., when $C_{4}b_{0}C_{{\tiny l}}^{-1}n^{-1/2}\big{(}C_{{\tiny
u}}\tilde{s}\log(2d/\tilde{s})\big{)}^{1/2}\leq\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{*}\|_{2}\leq
8C_{{\tiny l}}^{-1/2}\gamma$, by taking a constant stepsize $\eta_{l}=\eta\in
C_{{\tiny l}}^{1/2}b_{1}^{2}(nb_{0}C_{{\tiny
u}})^{-1}\cdot\left[C_{5},C_{6}\right]$, then
$\|{\boldsymbol{\beta}}_{l+1}-{\boldsymbol{\beta}}^{*}\|_{2}\leq(1-c_{2}^{*})\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{*}\|_{2}$

Therefore, after at most
$O\big{(}\log(\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{*}\|_{2}/\gamma)+\log(n\gamma
b_{0}^{-1}\log(2d/\tilde{s}))\big{)}$ iterations, Algorithm 1 outputs an
estimator with the error rate

$\|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}^{2}\leq\frac{C_{4}C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot\frac{b_{0}^{2}\tilde{s}\log(2d/\tilde{s})}{n}.$

Note that in Theorem 1 if a constant depends on $b_{0}$ or/and $b_{1}$, a star
sign is placed on its top right corner. The phase one convergence is typical
in non-smooth optimization which requires decaying stepsizes; see also Foucart
and Lecué (2017). Surprisingly, the phase two convergence requires a constant
stepsize as if minimizing a smooth objective function (Blumensath and Davies,
2009; Chen and Wainwright, 2015).

##### Gaussian noise

If the noise is Gaussian $N(0,\sigma^{2})$, we have $b_{0}\asymp
b_{1}\asymp\sigma$. If further $C_{{\tiny u}}\asymp C_{{\tiny l}}\asymp 1$, we
can choose a sparsity level $\tilde{s}\gtrsim s$ so that the output of
Algorithm 1 attains the rate
$\|\hat{{\boldsymbol{\beta}}}-{\boldsymbol{\beta}}^{*}\|_{2}^{2}=O_{p}\big{(}\sigma^{2}\cdot
sn^{-1}\log(d/s)\big{)}$, which is minimax optimal (Rigollet and Tsybakov,
2011; Ye and Zhang, 2010; Raskutti et al., 2011; Lounici et al., 2011).

##### Heavy-tailed noise

Under a finite second moment condition, Sun et al. (2020) proved that a Huber-
loss-based estimator can achieve the error rate $O\big{(}s/n\cdot\log
d\big{)}$ with probability exceeding $1-d^{-\Omega(1)}$. See also related
works in Pan et al. (2021); Wang et al. (2021, 2020). A median-of-means
estimator was studied in Depersin (2020) showing that, under a finite second
moment condition, the rate $O\big{(}s/n\cdot\log(d/s)\big{)}$ can be attained
with probability exceeding $1-\Omega(1)\cdot\exp(-s\log(ed/s))$. Our Theorem 1
only requires a finite $1+\varepsilon$ moment condition, but the attained
error rate is comparable to Depersin (2020). Though our estimator is not sub-
Gaussian, most existing sub-Gaussian estimators require a finite second moment
condition (Catoni, 2012; Lugosi and Mendelson, 2019c, a; Hopkins, 2020;
Depersin and Lecué, 2022). Chinot et al. (2020) relaxed the noise moment
condition but it is unclear how their estimator can be efficiently computed
with theoretical guarantees. More importantly, unlike the aforementioned
works, our estimator is attained by a provably fast converging algorithm.

###### Proof sketch of Theorem 1.

Note that not all the entries of the sub-gradient ${\mathbf{G}}_{l}$ affect
${\boldsymbol{\beta}}_{l+1}$. Denote

$\Omega_{l+1}:=\text{supp}({\boldsymbol{\beta}}_{l+1}),\quad\Omega_{l}:=\text{supp}({\boldsymbol{\beta}}_{l}),\quad\Pi_{l}:=\text{supp}\big{(}{\mathcal{H}}_{\tilde{s}}({\mathcal{P}}_{\Omega_{l}^{\rm
c}}({\mathbf{G}}_{l}))\big{)},\quad\Omega^{*}:=\text{supp}({\boldsymbol{\beta}}^{*}),$

where the operator ${\mathcal{P}}_{\Omega}({\mathbf{G}})$ keeps the entries
${\mathbf{G}}$ in $\Omega$ unchanged and zeros out the others and $\Omega^{c}$
is complement of $\Omega$. It is clear that
$\Omega_{l+1}\subseteq\Omega_{l}\cup\Pi_{l}\subseteq\Omega_{l}\cup\Pi_{l}\cup\Omega^{*}$.
We can re-write the update rule as

$\displaystyle{\boldsymbol{\beta}}_{l+1}={\mathcal{H}}_{\tilde{s}}\left({\boldsymbol{\beta}}_{l}-\eta_{l}{\mathbf{G}}_{l}\right)={\mathcal{H}}_{\tilde{s}}\left({\boldsymbol{\beta}}_{l}-\eta_{l}{\mathcal{P}}_{\Omega_{l}\cup\Pi_{l}\cup\Omega^{*}}({\mathbf{G}}_{l})\right).$

It thus suffices to focus on the sparse vector
${\mathcal{P}}_{\Omega_{l}\cup\Pi_{l}\cup\Omega^{*}}({\mathbf{G}}_{l})$. The
following lemma depicts the two-phase regularity properties of the loss
$f({\boldsymbol{\beta}})$ around the oracle ${\boldsymbol{\beta}}^{*}$.

######  Lemma 1.

Suppose Assumptions 1 and 2 hold. There exist absolute constants
$C_{1},C_{2},C_{3},c_{0}>0$ such that if $n\geq C_{1}C_{{\tiny u}}C_{{\tiny
l}}^{-1}\tilde{s}\log(2d/\tilde{s})$, then with probability over
$1-\exp\big{(}-c_{0}\tilde{s}\log\big{(}2d/\tilde{s})\big{)}-3\exp\big{(}-(n\log(2d/\tilde{s}))^{1/2}\log^{-1}n\big{)}$,
we have

  1. (1)

for all
${\boldsymbol{\beta}}\in\left\\{{\boldsymbol{\beta}}\in{\mathbb{R}}^{d}:\;\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\geq
8C_{{\tiny l}}^{-1/2}\gamma,\
|\text{supp}({\boldsymbol{\beta}})|\leq\tilde{s}\right\\}$ and for all sub-
gradient ${\mathbf{G}}\in\partial f({\boldsymbol{\beta}})$,

$\displaystyle
f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})\geq\frac{n}{4}C_{{\tiny
l}}^{1/2}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\quad{\rm
and}\quad\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}\leq
nC_{{\tiny u}}^{1/2},$

where $\Omega:=\text{supp}({\boldsymbol{\beta}})$ and
$\Pi:=\text{supp}\big{(}{\mathcal{H}}_{\tilde{s}}({\mathcal{P}}_{\Omega^{\rm
c}}({\mathbf{G}}))\big{)}$;

  2. (2)

for all
${\boldsymbol{\beta}}\in\left\\{{\boldsymbol{\beta}}\in{\mathbb{R}}^{d}:\
C_{2}C_{{\tiny u}}^{1/2}C_{{\tiny l}}^{-1}\cdot
b_{0}\big{(}\tilde{s}/n\cdot\log(2d/\tilde{s})\big{)}^{1/2}\leq\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\leq
8C_{{\tiny l}}^{-1/2}\gamma,\
|\text{supp}({\boldsymbol{\beta}})|\leq\tilde{s}\right\\}$ and for all sub-
gradient ${\mathbf{G}}\in\partial f({\boldsymbol{\beta}})$,

$\displaystyle
f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})\geq\frac{nC_{{\tiny
l}}}{12b_{0}}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}^{2}\quad{\rm
and}\quad\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}\leq
C_{3}\frac{nC_{{\tiny
u}}}{b_{1}}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2},$

where $\Omega:=\text{supp}({\boldsymbol{\beta}})$ and
$\Pi:=\text{supp}\big{(}{\mathcal{H}}_{\tilde{s}}({\mathcal{P}}_{\Omega^{\rm
c}}({\mathbf{G}}))\big{)}$.

Based on Lemma 1, we show that during phase one

$\|{\boldsymbol{\beta}}_{l}-\eta_{l}{\mathcal{P}}_{\Omega_{l}\cup\Pi_{l}\cup\Omega^{\ast}}({\mathbf{G}}_{l})-{\boldsymbol{\beta}}^{\ast}\|_{2}^{2}\leq\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{\ast}\|_{2}^{2}-\frac{n\eta_{l}C_{{\tiny
l}}^{1/2}}{2}\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{\ast}\|_{2}+n^{2}\eta_{l}^{2}C_{{\tiny
u}}.$

The selected stepsize guarantees $\eta_{l}\asymp C_{{\tiny l}}^{1/2}C_{{\tiny
u}}^{-1}\cdot
n^{-1}(1-c_{1})^{l}\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{\ast}\|_{2}$,
by which the contraction of
$\|{\boldsymbol{\beta}}_{l+1}-{\boldsymbol{\beta}}^{\ast}\|_{2}$ is derived
using the fact
${\boldsymbol{\beta}}_{l+1}={\mathcal{H}}_{\tilde{s}}\left({\boldsymbol{\beta}}_{l}-\eta_{l}{\mathcal{P}}_{\Omega_{l}\cup\Pi_{l}\cup\Omega^{*}}({\mathbf{G}}_{l})\right)$.
During phase two, by Lemma 1, we get

$\|{\boldsymbol{\beta}}_{l}-\eta_{l}{\mathcal{P}}_{\Omega_{l}\cup\Pi_{l}\cup\Omega^{\ast}}({\mathbf{G}}_{l})-{\boldsymbol{\beta}}^{\ast}\|_{2}^{2}\leq\bigg{(}1-\frac{\eta_{l}nC_{{\tiny
l}}}{6b_{0}}+\frac{\eta_{l}^{2}C_{3}^{2}n^{2}C_{{\tiny
u}}^{2}}{b_{1}^{2}}\bigg{)}\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{\ast}\|_{2}^{2}.$

Similarly, the selected constant stepsize guarantees the contraction of
$\|{\boldsymbol{\beta}}_{l+1}-{\boldsymbol{\beta}}^{\ast}\|_{2}$. ∎

######  Remark 1 (Smoothing effect of random noise).

By Theorem 1, the phase two convergence of Algorithm 1 requires a constant
stepsize as if minimizing a smooth function. Lemma 1 shows that, in the small
region around ${\boldsymbol{\beta}}^{\ast}$ where the phase two convergence
occurs, the loss function $f({\boldsymbol{\beta}})$ indeed possesses the
regularity conditions which usually appear in smooth optimization, e.g., the
square loss (Cai et al., 2021a; Zhao et al., 2015) and logistic loss (Lyu et
al., 2021). The second phase occurs when the individual random noise starts to
(stochastically) dominate
$\|{\boldsymbol{\beta}}_{l}-{\boldsymbol{\beta}}^{\ast}\|_{2}$. This suggests
that, when ${\boldsymbol{\beta}}^{\ast}$ is close enough to
${\boldsymbol{\beta}}^{\ast}$, the random noise has the effect of smoothing
the objective function. In Figure 1, we present the shape of the stochastic
function $g(t):=n^{-1}\sum_{i=1}^{n}|\xi_{i}-t|$ and also its sub-differential
where $\xi_{i}$’s are i.i.d. $N(0,\tau^{2})$. The function is clearly non-
smooth when $\tau=0$ and becomes smoother as $\tau$ increases. Randomized
smoothing has been observed in optimization and statistics literature. For
instance, Zhang and Xia (2020) proved that the random edge-wise observational
error can naturally smooth the Edgeworth expansion of the network moment
statistics, whose empirical c.d.f. can be a step function without the
smoothing effect. In order to optimize a non-smooth convex function
$h({\boldsymbol{\beta}})$, e.g. $\|{\boldsymbol{\beta}}\|_{1}$, with optimal
convergence rates, Duchi et al. (2012) proposed to optimize the randomly-
smoothed version ${\mathbb{E}}h({\boldsymbol{\beta}}+Z)$ where $Z$ is some
random variable with a known distribution.

(a) Value of $\frac{1}{n}\sum_{i=1}^{n}|\xi_{i}-t|$ against $t$

(b) Value of $\frac{1}{n}\sum_{i=1}^{n}\text{sign}(t-\xi_{i})$ against $t$

Figure 1: Sample size $n=1000$ and $\\{\xi_{i}\\}_{i=1}^{n}$ follow
$N(0,0.1^{2})$, $N(0,1)$, $N(0,10^{2})$ respectively.

##  4 Low-rank Linear Regression

Let ${\mathbf{M}}^{\ast}$ be a $d_{1}\times d_{2}$ matrix with rank $r\ll
d_{2}$. Without loss of generality, we assume $d_{1}\geq d_{2}$. Denote
$\\{({\mathbf{X}}_{i},Y_{i})\\}_{i=1}^{n}$ the collection of i.i.d.
observations satisfying
$Y_{i}=\langle{\mathbf{X}}_{i},{\mathbf{M}}^{\ast}\rangle+\xi_{i}:={\rm
tr}({\mathbf{X}}_{i}^{\top}{\mathbf{M}}^{\ast})+\xi_{i}$ with
${\mathbf{X}}_{i}\in{\mathbb{R}}^{d_{1}\times d_{2}}$ known as the measurement
matrix. Our goal is to estimate ${\mathbf{M}}^{\ast}$ by solving the
optimization program

$\displaystyle\hat{\mathbf{M}}:=\underset{{\mathbf{M}}\in{\mathbb{M}}_{r}}{\arg\min}\
f({\mathbf{M}})\quad\textrm{ where
}f({\mathbf{M}}):=\sum_{i=1}^{n}\rho\big{(}\langle{\mathbf{M}},{\mathbf{X}}_{i}\rangle-
Y_{i}\big{)},$ (4)

where ${\mathbb{M}}_{r}:=\\{{\mathbf{M}}\in{\mathbb{R}}^{d_{1}\times
d_{2}}:\;\text{rank}({\mathbf{M}})\leq r\\}$. Here $\rho(\cdot)$ is a general
robust loss function including the absolute loss (Tong et al., 2021b;
Charisopoulos et al., 2021), Huber loss (Elsener and van de Geer, 2018; Sun et
al., 2020) and quantile loss (Alquier et al., 2019; Chen et al., 2022).

###  4.1 Algorithm: Riemannian optimization

Our algorithm is the projected sub-gradient descent on the Riemannian manifold
(RsGrad). At the $l$-th iteration with a current low-rank estimate
${\mathbf{M}}_{l}$, the algorithm consists of two major steps. It begins with
computing the Riemannian sub-gradient, which is the projection of a vanilla
sub-gradient ${\mathbf{G}}_{l}\in\partial
f({\mathbf{M}}_{l})\subset{\mathbb{R}}^{d_{1}\times d_{2}}$ onto the tangent
space of ${\mathbb{M}}_{r}$ at the point ${\mathbf{M}}_{l}$, denoted by
${\mathbb{T}}_{l}$. The second step is to update the low-rank estimate along
the direction of negative Riemannian sub-gradient and then retract it back to
the manifold ${\mathbb{M}}_{r}$, for which, it suffices to take the SVD. Here
${\rm SVD}_{r}(\cdot)$ returns the best rank-$r$ approximation by SVD. The
details can be found in Algorithm 2.

We note that the Riemannian sub-gradient is used only for the benefit of
computational efficiency. The theoretical results all hold if taking the
vanilla sub-gradient. The detailed computation of Riemannian sub-gradient can
be found in the supplementary file. RsGrad has its origin in Cambier and Absil
(2016) for minimizing the absolute loss without convergence and statistical
analysis. Our framework covers general robust loss functions, and we prove its
computational efficiency and statistical optimality in several important
applications.

##### Comparison with existing algorithms

Sub-gradient descent algorithms have been investigated by Tong et al. (2021b);
Charisopoulos et al. (2021) for low-rank linear regression with the absolute
loss. Our algorithm is different in two aspects: we take Riemannian sub-
gradient and have two phases of stepsize schedule. Factorization-based
optimizations (Zheng and Lafferty, 2016; Zhao et al., 2015) reparametrize the
program (4) by ${\mathbf{M}}={\mathbf{U}}{\mathbf{V}}^{\top}$. Then it
suffices to update $({\mathbf{U}},{\mathbf{V}})$ sequentially to minimize (2)
by sub-gradient descent algorithms. These algorithms can suffer a great loss
of computational efficiency if ${\mathbf{M}}^{\ast}$ is ill-conditioned.
Though the issue can be theoretically remedied by a proper inverse scaling
(Tong et al., 2021a), it can cause a potential computational instability
especially when ${\mathbf{M}}^{\ast}$ has small non-zero singular values.
Recently, it is discovered that optimizing (2) directly on the low-rank
manifold ${\mathbb{M}}_{r}$, called Riemannian, enjoys the fast computational
speed of factorization-based approaches and, meanwhile, converges linearly
regardless of the condition number of ${\mathbf{M}}^{\ast}$. See, e.g., Cai et
al. (2021a) and Cai et al. (2021b), for the convergence of Riemannian gradient
descent (RGrad) algorithms in minimizing strongly convex and smooth functions
with tensor-related applications. We note that RGrad is similar to the
projected gradient descent (PGD, Chen and Wainwright (2015)) except that RGrad
utilizes the Riemannian gradient while PGD takes the vanilla one.

Algorithm 2 Riemannian Sub-gradient Descent (RsGrad)

Input: observations $\\{({\mathbf{X}}_{i},Y_{i})\\}_{i=1}^{n}$, max iterations
${l_{\textsf{max}}}$, step sizes $\\{\eta_{l}\\}_{l=0}^{{l_{\textsf{max}}}}$.

Initialization: ${\mathbf{M}}_{0}\in{\mathbb{M}}_{r}$

for $l=0,\ldots,{l_{\textsf{max}}}$ do

Choose a vanilla subgradient: ${\mathbf{G}}_{l}\in\partial
f({\mathbf{M}}_{l})$

Compute the Riemannian sub-gradient:
$\widetilde{\mathbf{G}}_{l}={\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})$

Retraction to ${\mathbb{M}}_{r}$:
${\mathbf{M}}_{l+1}=\text{SVD}_{r}({\mathbf{M}}_{l}-\eta_{l}\widetilde{\mathbf{G}}_{l})$

end for

Output: $\hat{\mathbf{M}}={\mathbf{M}}_{{l_{\textsf{max}}}}$

###  4.2 General Convergence Performance

We now present the general convergence performance of RsGrad Algorithm 2,
which essentially relies on the regularity conditions of the objective
function. These conditions, in spirit, largely inherit those from existing
literature (Charisopoulos et al., 2021; Elsener and van de Geer, 2018; Alquier
et al., 2019; Tong et al., 2021a). However, as explained in Section 1, these
prior works only delivered statistically sub-optimal estimates. It turns out
that more delicate characterizations of these conditions are necessary for our
purpose. More exactly, we discover that the aforementioned robust functions
exhibit strikingly different regularity conditions near and far away from
${\mathbf{M}}^{\ast}$, referred to as the two-phase regularity conditions.
Here $\|\cdot\|_{\rm F}$ represents the Frobenius norm of a matrix.

######  Condition 1.

(Two-phase regularity conditions) Define the two regions around
${\mathbf{M}}^{\ast}$ with $\tau_{\textsf{\tiny comp}}>\tau_{\textsf{\tiny
stat}}>0$ by

${\mathbb{B}}_{1}:=\left\\{{\mathbf{M}}\in{\mathbb{M}}_{r}:\|{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm
F}\geq\tau_{\textsf{\tiny comp}}\right\\}\quad{\rm
and}\quad{\mathbb{B}}_{2}:=\left\\{{\mathbf{M}}\in{\mathbb{M}}_{r}:\tau_{\textsf{\tiny
comp}}>\|{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}\geq\tau_{\textsf{\tiny
stat}}\right\\}.$

The function $f(\cdot):{\mathbb{R}}^{d_{1}\times
d_{2}}\mapsto{\mathbb{R}}_{+}$ satisfies

  1. (1)

the rank-$r$ restricted $(\tau_{\textsf{\tiny comp}},\tau_{\textsf{\tiny
stat}},\mu_{\textsf{\tiny comp}},\mu_{\textsf{\tiny stat}})$ two-phase
sharpness with respect to ${\mathbf{M}}^{\ast}$,

$f({\mathbf{M}})-f({\mathbf{M}}^{\ast})\geq\begin{cases}\mu_{\textsf{\tiny
comp}}\|{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F},&\textrm{ for }\ \
{\mathbf{M}}\in{\mathbb{B}}_{1};\\\ \mu_{\textsf{\tiny
stat}}\|{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2},&\textrm{ for }\ \
{\mathbf{M}}\in{\mathbb{B}}_{2};\end{cases}$

  2. (2)

the rank-r restricted $(\tau_{\textsf{\tiny comp}},\tau_{\textsf{\tiny
stat}},L_{\textsf{\tiny comp}},L_{\textsf{\tiny stat}})$ two-phase sub-
gradient bound with respect to ${\mathbf{M}}^{\ast}$,

$\|{\mathbf{G}}\|_{\rm F,r}\leq\begin{cases}L_{\textsf{\tiny comp}},&\textrm{
for }\ \ {\mathbf{M}}\in{\mathbb{B}}_{1};\\\ L_{\textsf{\tiny
stat}}\|{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F},&\textrm{ for }\ \
{\mathbf{M}}\in{\mathbb{B}}_{2},\end{cases}$

where ${\mathbf{G}}\in\partial f({\mathbf{M}})$ and the truncated Frobenius
norm $\|{\mathbf{G}}\|_{\rm F,r}:=\|{\rm SVD}_{r}({\mathbf{G}})\|_{\rm F}$.

Basically, Condition 1 dictates the distinct behaviors of $f(\cdot)$ in two
neighbourhoods of ${\mathbf{M}}^{\ast}$. The quantity $\tau_{\textsf{\tiny
stat}}$ reflects the statistical limit while $\tau_{\textsf{\tiny comp}}$ is
usually the rate achieved by the computational analysis in existing literature
(Tong et al., 2021a; Charisopoulos et al., 2021) without assuming noise
distributions. These prior works only reveal the first phase regularity
conditions, i.e., on ${\mathbb{B}}_{1}$. Oftentimes, we have
$\tau_{\textsf{\tiny comp}}\gg\tau_{\textsf{\tiny stat}}$, in which case these
prior works only deliver statistically sub-optimal estimates. Deriving the
second phase regularity condition is challenging where more precise
calculations are necessary. We remark that, in statistics literature, the sub-
gradient bound is related to the Lipschitz continuity (Alquier et al., 2019)
and the sharpness condition is called the one-point-margin condition (Elsener
and van de Geer, 2018). However, their estimators are built upon convex
program so that the analysis is only made in a small neighbour of
${\mathbf{M}}^{\ast}$ where the two-phase regularity conditions are not
pivotal. We summarize the two-phase regularity conditions for various
applications in Table 1. Proposition 1 characterizes the general dynamic of
Algorithm 2.

Loss function | Noise | $\tau_{\textsf{\tiny comp}}$ | $\tau_{\textsf{\tiny stat}}$ | $\mu_{\textsf{\tiny comp}}$ | $L_{\textsf{\tiny comp}}$ | $\mu_{\textsf{\tiny stat}}$ | $L_{\textsf{\tiny stat}}$
---|---|---|---|---|---|---|---
$\ell_{1}$-Loss | $N(0,\sigma^{2})$ | $\sigma$ | $O\left(\sigma\left(\frac{rd_{1}}{n}\right)^{1/2}\right)$ | $\frac{n}{12}$ | $2n$ | $\frac{n}{12\sigma}$ | $O\left(\frac{n}{\sigma}\right)$
$\ell_{1}$-Loss | Heavy-tailed | $8\gamma$ | $O\left(b_{0}\left(\frac{rd_{1}}{n}\right)^{1/2}\right)$ | $\frac{n}{4}$ | $2n$ | $\frac{n}{12b_{0}}$ | $O\left(\frac{n}{b_{1}}\right)$
Huber $\rho_{H,\delta}$ | Heavy-tailed | $8\gamma+2\delta$ | $O\left(b_{0}\left(\frac{rd_{1}}{n}\right)^{1/2}\right)$ | $\frac{\delta n}{2}$ | $4\delta n$ | $\frac{\delta n}{3b_{0}}$ | $O\left(\frac{\delta n}{b_{1}}\right)$
Quantile $\rho_{Q,\delta}$ | Heavy-tailed | $8\gamma$ | $O\left(b_{0}\left(\frac{rd_{1}}{n}\right)^{1/2}\right)$ | $\frac{n}{8}$ | $n$ | $\frac{n}{24b_{0}}$ | $O\left(\frac{n}{b_{1}}\right)$
Table 1: Summary of the two phase regularity conditions. In fact,
$\tau_{\textsf{\tiny comp}}$ is the error rate achieved after the phase one
iterations and $\tau_{\textsf{\tiny stat}}$ is the error rate achieved after
phase two iterations, i.e., the ultimate statistical accuracy. Here
$\gamma,b_{0}$ and $b_{1}$ are as described in Assumption 2. Note that
$b_{0}\asymp\gamma$ in many cases such as the Gaussian noise, zero symmetric
Pareto noise, and Student’s t noise, etc. By viewing $b_{0}/\gamma$ as a
constant, we have $\tau_{\textsf{\tiny stat}}\ll\tau_{\textsf{\tiny comp}}$ if
the sample size is large.

######  Proposition 1.

Suppose Condition 1 holds, the initialization
${\mathbf{M}}_{0}\in{\mathbb{M}}_{r}$ satisfies
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq
c_{0}\sigma_{r}({\mathbf{M}}^{\ast})\cdot\min\\{\mu_{\textsf{\tiny
stat}}^{2}L_{\textsf{\tiny stat}}^{-2},\mu_{\textsf{\tiny
comp}}^{2}L_{\textsf{\tiny comp}}^{-2}\\}$ for a small but absolute constant
$c_{0}>0$, and the initial stepsize
$\eta_{0}\in\left[0.2\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\mu_{\textsf{\tiny
comp}}L_{\textsf{\tiny comp}}^{-2},\
0.3\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\mu_{\textsf{\tiny
comp}}L_{\textsf{\tiny comp}}^{-2}\right]$. At the $l$-th iteration of
Algorithm 2,

  1. (1)

when $\|{\mathbf{M}}_{l}-{\mathbf{M}}^{*}\|_{\rm{F}}\geq\tau_{\textsf{\tiny
comp}}$, namely in phase one, take the stepsize
$\eta_{l}=\big{(}1-c\mu_{\textsf{\tiny comp}}^{2}L_{\textsf{\tiny
comp}}^{-2}\big{)}^{l}\eta_{0}$ with any fixed constant $c\leq 0.04$, then we
have

$\displaystyle\|{\mathbf{M}}_{l+1}-{\mathbf{M}}^{*}\|_{\mathrm{F}}\leq\left(1-c\frac{\mu_{\textsf{\tiny
comp}}^{2}}{L_{\textsf{\tiny
comp}}^{2}}\right)^{l+1}\cdot\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm F};$

  2. (2)

when $\tau_{\textsf{\tiny
comp}}>\|{\mathbf{M}}_{l}-{\mathbf{M}}^{*}\|_{\rm{F}}\geq\tau_{\textsf{\tiny
stat}}$, namely in phase two, take the stepsize
$\eta_{l}=\eta\in\mu_{\textsf{\tiny stat}}L_{\textsf{\tiny
stat}}^{-2}\cdot\big{[}0.125,\ 0.75\big{]}$, then we have

$\displaystyle\|{\mathbf{M}}_{l+1}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq\left(1-\frac{\mu_{\textsf{\tiny
stat}}^{2}}{32L_{\textsf{\tiny
stat}}^{2}}\right)\cdot\|{\mathbf{M}}_{l}-{\mathbf{M}}^{*}\|_{\rm{F}}.$

By Proposition 1, the RsGrad Algorithm 2 converges linearly using the proposed
stepsizes and outputs an estimate within an $O(\tau_{\textsf{\tiny stat}})$
distance in Frobenius norm from the oracle. As discussed above, the phase one
convergence is typical in non-smooth optimization, for which the decaying
stepsizes are needed. For simplicity, we apply the geometrically decaying
stepsizes (Tong et al., 2021a; Goffin, 1977). The phase two convergence
behaves like smooth optimization so that a fixed stepsize suffices to yield
the linear convergence. Note that the convergence dynamic is free of the
matrix condition number, which is a benefit of Riemannian-type algorithms.

##### Comparison with prior works

The algorithmic dynamic of RsGrad in phase one is similar to that in
Charisopoulos et al. (2021) and Tong et al. (2021b), where the geometrically
decaying stepsizes in phase one will eventually bring the stepsize into the
level $\Omega(\mu_{\textsf{\tiny stat}}L_{\textsf{\tiny stat}}^{-2})$ which is
desired by the phase-two convergence. However, Charisopoulos et al. (2021) and
Tong et al. (2021b) will continue to shrink the stepsize geometrically and
fail to deliver a statistically optimal estimator. In contrast, we set a
constant stepsize in phase two and the ultimate estimator is statistically
optimal. The difference is observed in numerical experiments in Section 7.

###  4.3 Applications

This section presents the statistical performances of RsGrad in some concrete
examples. Due to page limit, some examples are deferred to the supplementary
file. We note that, compared with the square loss (Xia et al., 2021; Zhao et
al., 2015; Chen and Wainwright, 2015), low-rank regression using the robust
loss functions requires no signal-to-noise ratio (SNR) conditions (except for
initialization).

####  4.3.1 Absolute loss with Gaussian noise

Though motivated by heavy-tailed noise, to start with, we demonstrate that the
absolute loss can yield a statistically optimal estimator under Gaussian
noise. The objective function is
$f({\mathbf{M}})=\sum_{i=1}^{n}|Y_{i}-\langle{\mathbf{M}},{\mathbf{X}}_{i}\rangle|$.

######  Theorem 2.

Assume $\xi_{1},\dots,\xi_{n}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}N(0,\sigma^{2})$ and $\\{{\rm
vec}({\mathbf{X}}_{i})\\}_{i=1}^{n}$ satisfy Assumption 1. There exist
absolute constants $c_{0},c_{1},c_{2}\in(0,1),C_{1},C_{2}>0$ such that if the
sample size $n\geq C_{1}C_{{\tiny u}}C_{{\tiny l}}^{-1}d_{1}r$, the
initialization $\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq
c_{0}\sigma_{r}$ for ${\mathbf{M}}_{0}\in{\mathbb{M}}_{r}$, the initial
stepsize $\eta_{0}\in C_{{\tiny l}}^{1/2}(48nC_{{\tiny
u}})^{-1}\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\cdot\left[0.2,0.3\right]$,
the phase one stepsize $\eta_{l}=(1-c_{3})^{l}\eta_{0}$ with any fixed
constant $c_{3}\leq c_{1}C_{{\tiny l}}C_{{\tiny u}}^{-1}$, and the phase two
stepsize $\eta_{l}=\eta\asymp C_{{\tiny l}}C_{{\tiny u}}^{-2}\sigma/n$, then
with probability at least $1-\exp(-c_{2}rd_{1})-3\exp(-\sqrt{n}/\log n)$,
after at most $O\big{(}\log(\sigma_{r}/\sigma)+\log(n/d_{1})\big{)}$
iterations, the Algorithm 2 outputs an estimator with the error rate

$\|\hat{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2}\leq C_{2}\frac{C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot\frac{\sigma^{2}d_{1}r}{n}.$

Note that the initial stepsize $\eta_{0}$ can be as large as
$O(\sigma_{r}n^{-1})$ depending on the initialization. The attained error rate
is minimax optimal (Ma and Wu, 2015; Xia, 2014). In comparison, the algorithms
proposed in Tong et al. (2021b); Charisopoulos et al. (2021) and their
analysis only achieved the sub-optimal rate $O_{p}(\sigma^{2})$. The proof of
Theorem 2 is a combination of the Proposition 1 and the two-phase regularity
conditions of $f({\mathbf{M}})$. The two phase convergence dynamics are
observed in numerical experiments in Section 7.

####  4.3.2 Absolute loss with heavy-tailed noise

We now demonstrate the effectiveness of the absolute loss in handling heavy-
tailed noise. Note that in Theorem 3 if a constant factor depends on $b_{0}$
or/and $b_{1}$, a star sign is placed on top right of it.

######  Theorem 3.

Assume $\\{{\rm vec}({\mathbf{X}}_{i})\\}_{i=1}^{n}$ and
$\\{\xi_{i}\\}_{i=1}^{n}$ satisfy Assumptions 1 and 2, respectively. There
exist constants $c_{0}^{\ast},c_{1},c_{2}\in(0,1),C_{1},C_{2}>0$ such that if
the sample size $n\geq C_{1}C_{{\tiny u}}C_{{\tiny l}}^{-1}d_{1}r$, the
initialization $\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq
c_{0}^{\ast}\sigma_{r}$ for ${\mathbf{M}}_{0}\in{\mathbb{M}}_{r}$, the initial
stepsize $\eta_{0}\in C_{{\tiny l}}^{1/2}(16nC_{{\tiny
u}})^{-1}\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\cdot\left[0.2,0.3\right]$,
the phase one stepsize $\eta_{l}=(1-c_{3})^{l}\eta_{0}$ with any fixed
constant $c_{3}\leq c_{1}$, and the phase two stepsize $\eta_{l}=\eta\asymp
C_{{\tiny l}}C_{{\tiny u}}^{-2}b_{1}^{2}(nb_{0})^{-1}$, then with probability
at least $1-\exp(-c_{2}rd_{1})-3\exp(-\sqrt{n}/\log n)$, after at most
$O\big{(}\log(\sigma_{r}/\gamma)+\log\big{(}\gamma
n/(d_{1}rb_{0})\big{)}\big{)}$ iterations, the Algorithm 2 outputs an
estimator with the error rate

$\|\hat{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2}\leq\frac{C_{2}C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot\frac{b_{0}^{2}d_{1}r}{n}.$

The proof of Theorem 3 is an immediate consequence of Proposition 1 if we have
the two-phase regularity conditions on $f({\mathbf{M}})$ under heavy-tailed
noise, for which we borrowed a powerful tool of bounding an empirical process
(Adamczak, 2008; Van de Geer and van de Geer, 2000). This allows us to derive
a stronger tail probability bound by avoiding bounding the term
$\sum_{i=1}^{n}|\xi_{i}|$. Indeed, since $\xi$ has a finite $1+\varepsilon$
moment, bounding the sum of $|\xi_{i}|$ results into a tail probability
$1-O(n^{-\varepsilon})$. This is not a desirable tail probability for a small
$\varepsilon$. Note that the Gaussian noise also satisfies Assumption 2, in
which case the attained rate in Theorem 3 has a matching minimax lower bound.
The rate is minimax optimal w.r.t. the dimension and sample size. Note that
$b_{0}\asymp\gamma={\mathbb{E}}|\xi|$ in many cases such as the Gaussian, zero
symmetric Pareto, and Student’s t noise.

The density requirement in Assumption 2 is similar to the one in Elsener and
van de Geer (2018) whose estimator is based on the nuclear-norm penalized
convex program without considering the computational efficiency. The same
error rate $O(b_{0}^{2}d_{1}r/n)$ was also derived there. But our approach
computes fast and requires only a logarithmic factor of iterations. The
analysis frameworks in Tong et al. (2021b); Charisopoulos et al. (2021) would
produce the sub-optimal rate $O(\gamma^{2})$ under the heavy-tailed noise,
which won’t decrease as the sample size $n$ gets larger. The error rates
derived by Fan et al. (2021); Minsker (2015) and Minsker (2018) were not
proportional to the noise level and were not zero even if the observations are
noiseless. In contrast, our algorithm can exactly recover
${\mathbf{M}}^{\ast}$ when there is no noise.

####  4.3.3 Huber loss with heavy-tailed noise

Huber loss is prevalent in robust statistics (Huber, 1965; Sun et al., 2020;
Elsener and van de Geer, 2018) and defined by
$\rho_{H,\delta}(x):=x^{2}\mathbbm{1}(|x|\leq\delta)+(2\delta|x|-\delta^{2})\mathbbm{1}(|x|>\delta)$
where $\delta>0$ is referred to as the robustification parameter. The function
$\rho_{H,\delta}(\cdot)$ is Lipschitz with a constant $2\delta$. The loss
function is given by
$f({\mathbf{M}}):=\sum_{i=1}^{n}\rho_{H,\delta}(Y_{i}-\langle{\mathbf{M}},{\mathbf{X}}_{i}\rangle)$.
Due to technical issues, we need a slightly different assumption of the heavy-
tailed noise.

######  Assumption 3.

(Heavy-tailed noise 2) There exists an $\varepsilon>0$ such that
${\mathbb{E}}|\xi|^{1+\varepsilon}<+\infty$. The noise distribution function
$H_{\xi}(\cdot)$ is symmetric555The symmetric condition can also be replaced
by $\int_{-\delta}^{\delta}H_{\xi}(x)\,dx=\delta$. in that
$H_{\xi}(x)=1-H_{\xi}(-x)$. Denote $\gamma={\mathbb{E}}|\xi|$. There exist
constants $b_{0},b_{1}$ (may be dependent on $\gamma$ and $\delta$) such that

$\displaystyle H_{\xi}(x+\delta)-H_{\xi}(x-\delta)\geq 2\delta b_{0}^{-1},$
$\displaystyle\ \ \ \textrm{ for all }|x|\leq 8(C_{{\tiny u}}/C_{{\tiny
l}})^{1/2}\gamma+2(C_{{\tiny u}}/C_{{\tiny l}})^{1/2}\delta;$ $\displaystyle
H_{\xi}(x+\delta)-H_{\xi}(x-\delta)\leq 2\delta b_{1}^{-1},$ $\displaystyle\ \
\ \forall x\in{\mathbb{R}},$

where $\delta$ is the Huber loss parameter.

Theorem 4 describes the convergence and statistical performance of RsGrad
Algorithm 2 for Huber loss. Similarly, the constants dependent on $b_{0}$
or/and $b_{1}$ are marked with a star.

######  Theorem 4.

Assume $\\{{\rm vec}({\mathbf{X}}_{i})\\}_{i=1}^{n}$ and
$\\{\xi_{i}\\}_{i=1}^{n}$ satisfy Assumptions 1 and 3, respectively. There
exist constants $c_{0}^{\ast},c_{1},c_{2}\in(0,1),C_{1},C_{2}>0$ such that if
the sample size $n\geq C_{1}C_{{\tiny u}}C_{{\tiny l}}^{-1}d_{1}r$, the
initialization $\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq
c_{0}^{\ast}\sigma_{r}$ for ${\mathbf{M}}_{0}\in{\mathbb{M}}_{r}$, the initial
stepsize $\eta_{0}\in C_{{\tiny l}}^{1/2}(32n\delta C_{{\tiny
u}})^{-1}\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\cdot\left[0.2,0.3\right]$,
the phase one stepsize $\eta_{l}=(1-c_{3})^{l}\eta_{0}$ with any fixed
constant $c_{3}\leq c_{1}$, and the phase two stepsize $\eta_{l}=\eta\asymp
C_{{\tiny l}}C_{{\tiny u}}^{-2}b_{1}^{2}(n\delta b_{0})^{-1}$, then with
probability at least $1-\exp(-c_{2}rd_{1})-3\exp(-\sqrt{n}/\log n)$, after at
most $O\big{(}\log(\sigma_{r}/\gamma)+\log\big{(}\gamma
n/(b_{0}d_{1}r)\big{)}\big{)}$ iterations, the Algorithm 2 outputs an
estimator with the error rate

$\|\hat{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2}\leq\frac{C_{2}C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot\frac{b_{0}^{2}d_{1}r}{n}.$

Note that Assumption 3 implies that $b_{0}$ is lower bounded by
$\Omega(\gamma+\delta)$ and thus the error rate is lower bounded by
$\Omega\big{(}(\delta^{2}+\gamma^{2})rd_{1}/n\big{)}$. This suggests an
interesting transition with respect to the robustification parameter: when
$\delta\lesssim\gamma$, increasing $\delta$ does not affect the error rate;
when $\delta\gtrsim\gamma$, the error rate appreciates if $\delta$ becomes
larger. Sun et al. (2020) derived the statistical performance of Huber loss
for linear regression without imposing regularity conditions on the noise
density resulting into an error rate that is slower than $O(n^{-1})$.

####  4.3.4 Quantile loss with heavy-tailed noise

Quantile loss was initially proposed by Koenker and Bassett Jr (1978) and has
been a popular loss function for robust statistics (Welsh, 1989; Koenker and
Hallock, 2001; Wang et al., 2012; Tan et al., 2022). Denote the quantile loss
by $\rho_{Q,\delta}:=\delta x\mathbbm{1}(x\geq 0)+(\delta-1)x\mathbbm{1}(x<0)$
for any $x\in{\mathbb{R}}$ with $\delta:={\mathbb{P}}(\xi\leq 0)$. Notice that
the function $\rho_{Q,\delta}(\cdot)$ is Lipschitz with a constant
$\max\\{\delta,1-\delta\\}$. The loss function is given by
$f({\mathbf{M}}):=\sum_{i=1}^{n}\rho_{Q,\delta}(Y_{i}-\langle{\mathbf{M}},{\mathbf{X}}_{i}\rangle)$.
Similar to the case of Huber loss, due to technical reasons, our analysis
framework needs a slightly different assumption on the heavy-tailed noise.

######  Assumption 4.

(Heavy-tailed noise 3) There exists an $\varepsilon>0$ such that
${\mathbb{E}}|\xi|^{1+\varepsilon}<+\infty$. Suppose $\delta:=H_{\xi}(0)$ lies
in $(0,1)$. Denote $\gamma={\mathbb{E}}|\xi|$. There exist constants
$b_{0},b_{1}>0$ (may be dependent on $\gamma$) such that

$\displaystyle h_{\xi}(x)\geq b_{0}^{-1},$ $\displaystyle\ \ \ \textrm{ for
all }|x|\leq 8(C_{{\tiny u}}/C_{{\tiny l}})^{1/2}\gamma;$ $\displaystyle
h_{\xi}(x)\leq b_{1}^{-1},$ $\displaystyle\ \ \ \forall x\in{\mathbb{R}}.$

When the noise has a zero median, namely $\delta=1/2$, Assumption 4 becomes
identical to the Assumption 2.

######  Theorem 5.

Assume $\\{\text{vec}({\mathbf{X}}_{i})\\}_{i=1}^{n}$ and
$\\{\xi_{i}\\}_{i=1}^{n}$ satisfy Assumptions 1 and 4, respectively. There
exist constants $c_{0}^{\ast},c_{1},c_{2}\in(0,1),C_{1},C_{2}>0$ such that if
sample size $n\geq C_{1}C_{{\tiny u}}C_{{\tiny l}}^{-1}d_{1}r$, the
initialization $\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq
c_{0}^{\ast}\sigma_{r}$, the initial stepsize $\eta_{0}\in C_{{\tiny
l}}^{1/2}(8nC_{{\tiny
u}})^{-1}\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\cdot[0.2,0.3]$, the
phase one stepsize $\eta_{l}=(1-c_{3})^{l}\eta_{0}$ with any fixed constant
$c_{3}\leq c_{1}$, and the phase two stepsize $\eta_{l}=\eta\asymp C_{{\tiny
l}}C_{{\tiny u}}^{-2}b_{1}^{2}(nb_{0})^{-1}$, then with probability over
$1-\exp(-c_{2}rd_{1})-3\exp(-\sqrt{n}/\log n)$, after at most
$O\big{(}\log(\sigma_{r}/\gamma)+\log\big{(}\max\\{\delta^{-1},(1-\delta)^{-1}\\}\cdot\gamma
n/(d_{1}rb_{0})\big{)}\big{)}$ iterations, the Algorithm 2 outputs an
estimator with the error rate

$\|\hat{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2}\leq\frac{C_{2}C_{{\tiny
u}}\cdot\max\\{\delta^{2},\ (1-\delta)^{2}\\}}{C_{{\tiny
l}}^{2}}\cdot\frac{b_{0}^{2}rd_{1}}{n}.$

##  5 Initialization and Algorithmic Parameters Selection

##### Initialization

The sparse linear regression Algorithm 1 allows arbitrary initialization. We
only discuss the initialization for low-rank linear regression. The
convergence of Algorithm 2 crucially relies on the warm initialization
${\mathbf{M}}_{0}$. One can simply apply the shrinkage low-rank approximation
(Fan et al., 2021) for initialization but a finite $2+\varepsilon$ moment
condition on noise is required. For simplicity, we investigate the performance
of the vanilla low-rank approximation. Here $vec(\cdot)$ and $mat(\cdot)$
denote the vectorization and matricization operators, respectively.

######  Theorem 6.

(Initialization Guarantees) Suppose
$\gamma_{1}:={\mathbb{E}}|\xi|^{1+\varepsilon}<+\infty$ for some
$\varepsilon\in(0,1]$ and $\\{vec({\mathbf{X}}_{i})\\}_{i=1}^{n}$ satisfy
Assumption 1 with known covariances
$\\{{\boldsymbol{\Sigma}}_{i}\\}_{i=1}^{n}$. Initialize by

${\mathbf{M}}_{0}:=\operatorname{SVD}_{r}(n^{-1}\sum_{i=1}^{n}\text{mat}({\boldsymbol{\Sigma}}_{i}^{-1}\text{vec}({\mathbf{X}}_{i}))Y_{i}).$

For any small $c_{0}>0$, there exist constants $C,c_{1},c_{2},c_{3}>0$
depending only on $c_{0}$ such that if $n>C\max\big{\\{}C_{{\tiny u}}C_{{\tiny
l}}^{-1}\kappa^{2}d_{1}r^{2}\log d_{1},\allowbreak(C_{{\tiny
l}}^{-1}d_{1}r)^{(1+\varepsilon)/(2\varepsilon)}\sigma_{r}^{-(1+\varepsilon)/\varepsilon}(\gamma_{1}\log
d_{1})^{1/\varepsilon}\big{\\}}$, then with probability over
$1-c_{1}d_{1}^{-1}-c_{2}\log^{-1}d_{1}\allowbreak-c_{3}e^{-d_{1}}$, the
initialization satisfies $\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq
c_{0}\sigma_{r}$. Here
$\kappa:=\sigma_{1}({\mathbf{M}}^{\ast})\sigma_{r}^{-1}({\mathbf{M}}^{\ast})$
is the condition number.

Compared with Fan et al. (2021), Theorem 6 only requires a finite
$1+\varepsilon$ moment condition. We note that if the noise term $\xi_{i}$ has
a finite $2+\varepsilon$ moment, the probability term $O(\log^{-1}d_{1})$ can
be replaced by $O(n^{-\min\\{1,\varepsilon/2\\}})$. Moreover, with a finite
$2+\varepsilon$ moment, one can even obtain a sub-Gaussian initial estimator
by the median-of-means (Minsker, 2015; Lugosi and Mendelson, 2019c, b) or the
Catoni’s estimator (Catoni, 2012; Minsker, 2018). Since the shrinkage-based
approaches (Fan et al., 2021; Minsker, 2018) are already effective to provide
warm initializations, we spare no additional efforts to improve Theorem 6. If
the covariances ${\boldsymbol{\Sigma}}_{i}$’s are unknown, we suggest to apply
Fan et al. (2021) or Minsker (2018) for obtaining the initialization.

##### Selection of initial stepsize

Though the main theorems in Sections 3 and 4 require the initial stepsize
$\eta_{0}\asymp
n^{-1}\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{\ast}\|_{2}$ and
$\eta_{0}\asymp n^{-1}\|{\mathbf{M}}_{0}-{\mathbf{M}}^{\ast}\|_{\rm F}$,
respectively, the following theorems show that similar results can be attained
even if we only have an upper bound of
$\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{\ast}\|_{2}$ and
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{\ast}\|_{\rm F}$.

######  Theorem 1* (alternative of Theorem 1).

Suppose Assumptions 1 and 2 hold. There exist absolute constants
$C_{1},\dots,C_{6}>0$, $c_{2}^{*}\in(0,1)$ such that if $n\geq C_{1}C_{{\tiny
u}}C_{{\tiny l}}^{-1}\tilde{s}\log(2d/\tilde{s})$, $\tilde{s}\geq
C_{2}(C_{{\tiny u}}/C_{{\tiny l}})^{8}(b_{0}/b_{1})^{8}\cdot s$, the initial
stepsize $\eta_{0}\in(nC_{{\tiny u}})^{-1}C_{{\tiny
l}}^{1/2}D_{0}\cdot\left[1/8,\ 3/8\right]$ for any
$D_{0}\geq\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{\ast}\|_{2}$, the
phase-one stepsize $\eta_{l}=(1-c_{1})^{l}\eta_{0}$ for any $c_{1}\leq
C_{{\tiny l}}C_{{\tiny u}}^{-1}/64$, and the phase-two stepsize
$\eta_{l}=\eta\in C_{{\tiny l}}^{1/2}b_{1}^{2}(nb_{0}C_{{\tiny
u}})^{-1}\cdot[C_{3},\ C_{4}]$, then with probability exceeding
$1-\exp\big{(}-C_{5}\tilde{s}\log(2d/\tilde{s})\big{)}-3\big{(}-(n\log(2d/\tilde{s})^{1/2}\log^{-1}n\big{)}$,
after at most $O\big{(}\log(D_{0}/\gamma)+\log(n\gamma
b_{0}^{-1}\log(2d/\tilde{s}))\big{)}$ iterations, the Algorithm 1 outputs an
estimator with the error rate

$\|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{\ast}\|_{2}^{2}\leq\frac{C_{6}C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot\frac{b_{0}^{2}\tilde{s}\log(2d/\tilde{s})}{n}.$

######  Theorem 3* (alternative of Theorem 3).

Assume $\\{\xi_{i}\\}_{i=1}^{n}$ and $\\{{\rm
vec}({\mathbf{X}}_{i})\\}_{i=1}^{n}$ satisfy Assumptions 1 and 2,
respectively. There exist constants
$c_{0}^{\ast},c_{1},c_{2}\in(0,1),C_{1},C_{2}>0$ such that if $n\geq
C_{1}C_{{\tiny u}}C_{{\tiny l}}^{-1}d_{1}r$, the initialization
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq D_{0}$ for any $D_{0}\leq
c_{0}^{\ast}\sigma_{r}$, the initial stepsize $\eta_{0}\in C_{{\tiny
l}}^{1/2}(16nC_{{\tiny u}})^{-1}D_{0}\cdot\left[0.2,0.3\right]$, the phase one
stepsize $\eta_{l}=(1-c_{3})^{l}\eta_{0}$ with any fixed constant $c_{3}\leq
c_{1}$, and the phase two stepsize $\eta_{l}=\eta\asymp C_{{\tiny l}}C_{{\tiny
u}}^{-2}b_{1}^{2}(nb_{0})^{-1}$, then with probability at least
$1-\exp(-c_{2}rd_{1})-3\exp(-n^{1/2}/\log n)$, after at most
$O\big{(}\log(D_{0}/\gamma)+\log\big{(}\gamma n/(d_{1}rb_{0})\big{)}\big{)}$
iterations, the Algorithm 2 outputs an estimator with the error rate

$\|\hat{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2}\leq\frac{C_{2}C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot\frac{b_{0}^{2}d_{1}r}{n}.$

By Theorem 3*, it thus suffices to obtain an estimate of
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{\ast}\|_{\rm F}$. Towards that end, one
appropriate method is to take the average
$n^{-1}\sum_{i=1}^{n}|Y_{i}-\langle{\mathbf{M}}_{0},{\mathbf{X}}_{i}\rangle|$,
which, under mild conditions, is of the same order of
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{\ast}\|_{\rm F}$ with high probability. The
same approach is also applicable in sparse vector regression. The general
quantities $\mu_{\textsf{\tiny comp}}$ and $L_{\textsf{\tiny comp}}$ involved
in Algorithm 2 are often $\asymp n$ up to a factor of loss-related parameters,
which can be easily tuned in practice.

##### Tuning Parameters in RsGrad

Our algorithm assumes knowing the rank of ${\mathbf{M}}^{\ast}$, which is
typical in non-convex methods(Tong et al., 2021a; Chen and Wainwright, 2015;
Cai et al., 2021a). One may apply the BIC or AIC for the selection of the rank
(Cai et al., 2021a). While we believe that our algorithm still works for the
over-parameterized case (Zhuo et al., 2021; Ma and Fattahi, 2022; Luo and
Zhang, 2022) where the algorithmic rank $r$ is larger than the true rank of
${\mathbf{M}}^{\ast}$, a thorough theoretical investigation is beyond the
scope of this paper. Convex methods such as Elsener and van de Geer (2018) and
Belloni et al. (2011) need no initialization or the pre-requisite of rank. But
our primary interest is to study both the computational and statistical issues
in robust estimation. Our algorithms (i.e., the phase two stepsize) also need
the noise level $\gamma$, which is not required by the convex approaches based
on the square-root LASSO or the nuclear-norm penalized absolute loss (Elsener
and van de Geer, 2018; Belloni et al., 2011; Minsker et al., 2022).
Fortunately, $\gamma$ can be reliably estimated by our methods. The phase one
iterations terminate until
$\sum_{i=1}^{n}|Y_{i}-\langle{\mathbf{M}}_{l_{1}},{\mathbf{X}}_{i}\rangle|\asymp\sum_{i=1}^{n}|\xi_{i}|\asymp
n\gamma$ where $l_{1}$ is the number of iterations in phase one. Then the
phase two stepsize can be taken as $\eta=c\cdot f({\mathbf{M}}_{l_{1}})/n^{2}$
for some tuning parameter $c>0$. In simulation experiments, we observe that
the algorithm is very tolerant with the selection of $c$.

##### Determine the phase

Due to the geometrically decay of stepsizes during phase one, after some
iterations, the stepsize will enter the level $O(\mu_{\textsf{\tiny
stat}}L_{\textsf{\tiny stat}}^{-2})$ desired by the phase two convergence. On
the other hand, if the phase one iterations continue, the stepsize will
diminish fast and the value of objective function becomes stable. This is
indeed observed in numerical experiments. See Figure 3 and Figure 4 in Section
7. Therefore, once the stepsize falls below a pre-chosen small threshold, the
phase two iterations can be initiated and the stepsizes are fixed afterwards
until convergence.

##  6 Huber’s Contamination Model

The robustness against outliers is often studied for the Huber’s
$\varepsilon$-contamination model (Huber and Ronchetti, 2009; Chen et al.,
2016, 2018; Diakonikolas and Kane, 2019), which assumes the data are i.i.d.
with a common distribution $d\,P=(1-\epsilon)\,d\,P_{i}+\epsilon\,d\,P_{o}$
where $P_{i}$ denotes the distribution of inliers, $P_{o}$ is an arbitrary
distribution of outliers, and $\epsilon\in[0,1)$. We now demonstrate the
robustness of our algorithms against outliers. Without loss of generality,
throughout this section, we assume $\epsilon$-fraction of the responses
$\\{Y_{i}\\}_{i=1}^{n}$ are corrupted with outliers. Denote the indices of
outliers as ${\mathcal{O}}$ and inliers as ${\mathcal{I}}$ where
${\mathcal{O}}\cup{\mathcal{I}}=\\{1,\dots,n\\}$ and
$|{\mathcal{O}}|=\lceil\epsilon n\rceil$. Basically,
$({\mathbf{X}}_{i},Y_{i})$ satisfies the linear model (1) for
$i\in{\mathcal{I}}$. For outliers $i\in{\mathcal{O}}$, the response $Y_{i}$ is
arbitrarily corrupted. Note that we still assume all the covariates
${\mathbf{X}}_{i}$’s or $vec({\mathbf{X}}_{i})$’s satisfy the Assumption 1.
Interested readers are suggested to refer to Depersin (2020); Lecué and
Lerasle (2019); Chinot et al. (2020); Thompson (2020); Diakonikolas et al.
(2019) for more general settings of Huber’s $\epsilon$-contamination model.
Due to page limit, we only consider the absolute loss.

###  6.1 Sparse linear regression

Using the same notations from Section 3, the following theorem describes the
statistical guarantees of the IHT-$\ell_{1}$ algorithm under Huber’s
$\epsilon$-contamination model.

######  Theorem 7.

Suppose Assumption 1 holds and $\\{\xi_{i}\\}_{i\in{\mathcal{I}}}$ satisfy
Assumption 2. There exist absolute constants $C_{1},\dots,C_{6}>0$,
$c_{2}^{*}\in(0,1)$ such that if $n\geq C_{1}C_{{\tiny u}}C_{{\tiny
l}}^{-1}\tilde{s}\log(2d/\tilde{s})$, $\tilde{s}\geq C_{2}(C_{{\tiny
u}}/C_{{\tiny l}})^{8}(b_{0}/b_{1})^{8}\cdot s$, the corruption rate
$\epsilon\leq C_{{\tiny l}}^{1/2}(4C_{{\tiny l}}^{1/2}+C_{{\tiny
u}}^{1/2})^{-1}$, the initial stepsize $\eta_{0}\in(nC_{{\tiny
u}})^{-1}C_{{\tiny
l}}^{1/2}\cdot\allowbreak\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{*}\|_{2}\cdot\left[1/8,\
3/8\right]$, the phase one stepsize $\eta_{l}=(1-c_{1})^{l}\eta_{0}$ with any
fixed constant $c_{1}\leq C_{{\tiny l}}C_{{\tiny u}}^{-1}/64$, and the phase
two stepsize $\eta_{l}=\eta\in C_{{\tiny l}}^{1/2}b_{1}^{2}(nb_{0}C_{{\tiny
u}})^{-1}\cdot\left[C_{5},C_{6}\right]$, then with probability over
$1-2\exp(-C_{3}\tilde{s}\log(2d/\tilde{s}))-3\exp(-(n\log(2d/\tilde{s}))^{1/2}\log^{-1}n)$,
after at most
$O\big{(}\log(\|{\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}^{*}\|_{2}/\gamma)+\log(n\gamma
b_{0}^{-1}\log(2d/\tilde{s}))\wedge\log(\gamma/b_{0}\epsilon)\big{)}$
iterations, Algorithm 1 outputs an estimator with the error rate

$\|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\leq\frac{C_{4}C_{{\tiny
u}}^{1/2}}{C_{{\tiny l}}}\cdot
b_{0}\cdot\max\left\\{\Big{(}\frac{\tilde{s}\log(2d/\tilde{s})}{n}\Big{)}^{1/2},\epsilon\right\\}.$

Our theorem allows $\epsilon$ to be any value smaller than $C_{{\tiny
l}}^{1/2}(4C_{{\tiny l}}^{1/2}+C_{{\tiny u}}^{1/2})^{-1}$. When the noise of
inliers is Gaussian, the error rate matches with the minimax lower bound (Chen
et al., 2016, 2018). The rates derived by Chinot et al. (2020); Depersin
(2020) have an $O(\epsilon^{1/2})$ term and the rate by Thompson (2020) has an
$O(\epsilon\log(1/\epsilon))$ term.

###  6.2 Low-rank linear regression

Using the same notations from Section 4, the following theorem describes the
statistical guarantees of the RsGrad algorithm under Huber’s
$\epsilon$-contamination model.

######  Theorem 8.

Suppose that $\\{{\rm vec}({\mathbf{X}}_{i})\\}_{i=1}^{n}$ and
$\\{\xi_{i}\\}_{i\in{\mathcal{I}}}$ satisfy Assumptions 1 and 2, respectively.
There exist constants
$c_{0}^{\ast},c_{1},c_{2}^{\ast}\in(0,1),C_{1}^{\ast},C_{2},C_{3}>0$ such that
if the sample size $n\geq C_{3}C_{{\tiny u}}C_{{\tiny l}}^{-1}d_{1}r$, the
corruption rate $\epsilon\leq C_{{\tiny l}}^{1/2}(4C_{{\tiny
l}}^{1/2}+C_{{\tiny u}}^{1/2})^{-1}$, the initialization
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\leq c_{0}^{\ast}\sigma_{r}$,
the initial stepsize $\eta_{0}\in C_{{\tiny l}}^{1/2}(16nC_{{\tiny
u}})^{-1}\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\cdot\left[0.2,0.3\right]$,
the phase one stepsize $\eta_{l}=(1-c_{3})^{l}\eta_{0}$ with any fixed
constant $c_{3}\leq c_{1}$, and the phase two stepsize $\eta_{l}=\eta\asymp
C_{{\tiny l}}C_{{\tiny u}}^{-2}b_{1}^{2}(nb_{0})^{-1}$, then with probability
at least $1-\exp(-c_{1}rd_{1})-3\exp(-\sqrt{n}/\log n)$, after at most
$O\big{(}\log(\sigma_{r}/\gamma)+C_{1}^{\ast}\log\big{(}\gamma
n/(d_{1}rb_{0})\big{)}\wedge\log(\gamma/b_{0}\epsilon)\big{)}$ iterations, the
Algorithm 2 outputs an estimator with the error rate

$\|\hat{\mathbf{M}}-{\mathbf{M}}^{\ast}\|_{\rm F}^{2}\leq\frac{C_{2}C_{{\tiny
u}}}{C_{{\tiny l}}^{2}}\cdot
b_{0}^{2}\cdot\max\left\\{\frac{d_{1}r}{n},\epsilon^{2}\right\\}.$

Exact recovery under sparse corruptions has been studied by Li et al. (2020).
Tong et al. (2021b) derived sub-optimal rates under sub-Gaussian noise and
sparse corruptions. Our method can deal with heavy-tailed noise and is
outlier-robust and minimax optimal.

##  7 Numerical Experiments

This section presents the numerical simulation results to showcase the two-
phase convergence phenomena and compare the attained accuracy with the
existing algorithms. Due to page limit, some results are relegated to the
supplement.

###  7.1 Sparse linear regression

Our algorithm is referred to as the IHT-$\ell_{1}$. We compare with with the
adaptive Huber method from Sun et al. (2020) and the truncation method from
Fan et al. (2021). The algorithmic parameters of adaptive Huber and truncation
methods are fine-tuned for their best performances. Note that our theoretical
analysis of IHT-$\ell_{1}$ can only treat the case $\tilde{s}>s$ because of
technical difficulties. The simulation results show that ITH-$\ell_{1}$ still
performs well when setting $\tilde{s}=s$. The two-phase convergence dynamic of
ITH-$\ell_{1}$ is similar to that of the low-rank linear regression and is
thus skipped.

##### Accuracy

We set ${\boldsymbol{\beta}}^{*}=(16,4,1,0,\dots,0)^{\top}$ where $d=50$ and
$s=3$. The covariate vectors $\\{{\mathbf{X}}_{i}\\}_{i=1}^{n}$ have i.i.d.
$N(0,1)$ entries. The sample sizes are chosen from $n\in\\{50,300\\}$. The
noise can be Gaussian or have a Student’s t distribution with d.f. $\nu=2$.
Figure 2 presents the box-plot of
$\|\hat{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{\ast}\|_{2}$ based on 50
simulations.

(a) Student’s t noise with d.f. $\nu=2$, $n=50$

(b) Student’s t noise with d.f. $\nu=2$, $n=300$

(c) Gaussian noise, $n=50$

(d) Gaussian noise, $n=300$

Figure 2: Accuracy comparisons of IHT-$\ell_{1}$, the adaptive Huber (Sun et
al., 2020) and truncation method (Fan et al., 2021) for sparse linear
regression.

##### Real data

We apply IHT-$\ell_{1}$ to the NCI-60 cancer cell line dataset. Note that
since the dataset has been updated recently, our experiment results may
slightly differ from Sun et al. (2020). The dataset has two parts: Protein:
Antibody Array DBT and RNA: Affy HG-U133(A-B) - GCRMA, both of which may be
found in website https://discover.nci.nih.gov/cellminer/loadDownload.do. The
first part Protein: Antibody Array DBT contains the levels of $162$ antibodies
in $60$ cell lines (Nishizuka et al., 2003). We focus on the antibody “KRT19”,
which is the most used marker for the tumor cells disseminated in lymph nodes,
peripheral blood, and bone marrow of breast cancer patients (Nakata et al.,
2004). The second part RNA: Affy HG-U133(A-B) - GCRMA reflects the gene
transcripts in $60$ cell lines (Shankavaram et al., 2007). Our goal is to
select the significant genes in generating the protein “KRT19”. By combining
overlapped genes and multiple tested genes and removing one missing data cell
line, we obtain a design matrix ${\mathbf{X}}$ with $n=59$ and $d=14,342$. See
Sun et al. (2020); Shankavaram et al. (2007) for more details about the data
set.

Method |  MAE |  Size |  Selected Genes
---|---|---|---
Adaptive Huber |  2.13 |  11 |  MALL, TM4SF4, ANXA3, ADRB2, NRN1, AUTS2, CA2, BAMBI, EPS8L2, CEMIP, KRT19
Truncation |  2.01 |  10 |  MALL, ANXA3, ADRB2, NRN1, CA2, BAMBI, COCH, CEMIP, KRT19, EMP3
IHT-$\ell_{1}$ ($\tilde{s}=7$) |  1.91 |  7 |  MALL, TGFBI, NRN1, DSP, EPS8L2, KRT19, LPXN
IHT-$\ell_{1}$ ($\tilde{s}=12$) |  2.12 |  12 |  S100A6, EPCAM, MALL, ZIC1, TGFBI, NRN1, DSP, DKK1, EPS8L2, MYO5C, KRT19, LPXN
Table 2: Prediction performance on NCI-60 data. The mean absolute error of
test data, final model size and selected genes are reported for three methods:
IHT-$\ell_{1}$, the adaptive Huber (Sun et al., 2020) and truncation method
(Fan et al., 2021).

We test the prediction performance and report the selected genes.We split the
data into a training set and a test set. The design matrix is $\log_{2}$
transformed and later centered. The algorithmic parameters of adaptive Huber,
truncation method, and IHT-$\ell_{1}$ are fine-tuned via 10-fold cross-
validations. The mean absolute error on test data, the final model size, and
the selected genes are reported in Table 2. Note that the selected genes are
obtained by applying the algorithms to the whole dataset. Here we report the
performance of IHT-$\ell_{1}$ only for the two cases $\tilde{s}=7$ and
$\tilde{s}=12$. Table 2 shows that the gene “KRT19” is significant, which has
been selected by all the algorithms. The genes “NRN1” and “EPS8L2” are
selected by all three algorithms, whose effect to breast cancer has been
verified by Kim et al. (2011); Colas et al. (2011). The gene “MALL” is also
selected by all algorithms, which is known to be related with some other
cancers so far. The genes “TGFBI”, “DSP”, “DKK1”, “S100A6”, “EPCAM”, “ZIC1”,
“MYO5C”, “LPXN” selected by IHT-$\ell_{1}$ are also important to breast cancer
(Li et al., 2012; Yang et al., 2012; Niu et al., 2019; Zhang et al., 2017;
Osta et al., 2004; Han et al., 2018; Ye et al., 2022; Abe et al., 2016).
Besides, we note that the genes “ANXA3”, “ADRB2”, “CA2”, “MAMBI”, “CEMIP”
selected by the adaptive Huber and truncation methods are significant by
existing literature (Zhou et al., 2017; Xie et al., 2019; Annan et al., 2019;
Yu and Chai, 2017; Xue et al., 2022).

###  7.2 Low-rank linear regression

We now showcase the two-phase convergence of RsGrad and compare the achieved
accuracy with the convex-huber, convex-$\ell_{1}$ (Elsener and van de Geer,
2018), convex-$\ell_{2}$ (Candes and Plan, 2011), and the (non-convex)
ScaledSM methods (Tong et al., 2021b). The accuracy comparison is deferred to
the supplement. For ScaledSM and the phase one iterations of RsGrad, the
stepsizes are set as $\eta_{0}=c_{1}\|{\mathbf{M}}_{0}\|n^{-1}$ and
$\eta_{l}=q\cdot\eta_{l-1}$ where $c_{1}>0$, $0<q<1$ are parameters to be
tuned and $\|\cdot\|$ denotes the matrix operator norm. As Tong et al. (2021b)
suggested, we choose $q=0.91$. For the phase two iterations of RsGrad, we set
the constant stepsize at $\eta_{l}=c_{2}{\mathbb{E}}|\xi|n^{-1}$ where
$c_{2}>0$ is a tuning parameter. Here we assume ${\mathbb{E}}|\xi|$ is known
for simplicity. In practice, ${\mathbb{E}}|\xi|$ can be estimated by
$n^{-1}\sum_{i=1}^{n}|Y_{i}-\langle{\mathbf{M}}_{l_{1}},{\mathbf{X}}_{i}\rangle|$
where ${\mathbf{M}}_{l_{1}}$ is the output after phase one iterations.
Following Tong et al. (2021b), denote the signal-to-noise ratio
$\text{SNR}:=20\log_{10}\big{(}\|{\mathbf{M}}^{*}\|_{\rm{F}}/{\mathbb{E}}|\xi|\big{)}$.
The relative error is calculated by
$\|{\mathbf{M}}_{l}-{\mathbf{M}}^{\ast}\|_{\rm F}\|{\mathbf{M}}^{\ast}\|_{\rm
F}^{-1}$.

##### Two-phase convergence

We compare the convergence dynamics between RsGrad and ScaledSM. We experiment
RsGrad with the absolute loss (RsGrad-$\ell_{1}$) and Huber loss (RsGrad-
Huber), but ScaledSM was proposed only for the absolute loss. The dimensions
and ranks are fixed at $d_{1}=d_{2}=80$ and rank $r=5$, respectively. The
phase one iterations are terminated when the stepsize becomes rather small,
e.g., smaller than $10^{-10}$, and then the phase two iterations are
initiated.

The first set of simulations is to showcase the convergence performance under
Gaussian noise. The SNR is set at $\\{40,80\\}$ and the sample sizes are
varied among $\\{2.5,5\\}\times rd_{1}$. The convergence dynamics are
displayed in Figure 3, which clearly show a two-phase convergence of RsGrad.
We note that ScaledSM performs poorly when sample size is small
($n=2.5rd_{1}$). It is possibly due to the instability of inverse scaling
under a small sample size. The bottom two plots in Figure 3 show that, after
the phase one iterations, RsGrad achieves a similar perfomrance as ScaledSM.
However, with a second phase convergence, RsGrad (using either Huber or
absolute loss) eventually delivers a more accurate estimate than ScaledSM.

(a) SNR=40

(b) SNR=80

(c) SNR=40

(d) SNR=80

Figure 3: Convergence dynamics of RsGrad and ScaledSM (Tong et al., 2021b)
under Gaussian noise. Top two figures: $n=2.5rd_{1}$; bottom two figures:
$n=5rd_{1}$.

The second set of simulation is to test the convergence performance under
heavy-tailed noise. For simplicity, the noise is sampled independently from a
Student’s $t$-distribution with ${\rm d.f.}$ $\nu=2$, which has a finite first
(absolute) moment and an infinite second moment. Other parameters are selected
the same as in the Gaussian case. Figure 4 presents the convergence
performance showing that both RsGrad (either absolute loss or Huber loss) and
ScaledSM are robust to heavy-tailed noise. Similarly, we can observe the two-
phase convergence of RsGrad whose second phase iterations deliver a more
accurate estimate than ScaledSM, especially when SNR is large.

(a) SNR=40

(b) SNR=80

(c) SNR=40

(d) SNR=80

Figure 4: Convergence dynamics of RsGrad and ScaledSM (Tong et al., 2021b)
under Student’s $t$-distribution with d.f. $\nu=2$. Top two figures:
$n=2.5rd_{1}$; bottom two figures: $n=5rd_{1}$.

## References

  * Abe et al. [2016] Akihiro Abe, Yukiya Yamamoto, Sachiko Iba, Tadaharu Kanie, Akinao Okamoto, Masutaka Tokuda, Yoko Inaguma, Masamitsu Yanada, Satoko Morishima, Shuichi Mizuta, et al.  Etv6-lpxn fusion transcript generated by t (11; 12)(q12. 1; p13) in a patient with relapsing acute myeloid leukemia with nup98-hoxa9.  _Genes, Chromosomes and Cancer_ , 55(3):242–250, 2016.
  * Absil et al. [2009] P-A Absil, Robert Mahony, and Rodolphe Sepulchre.  _Optimization algorithms on matrix manifolds_.  Princeton University Press, 2009.
  * Adamczak [2008] Radoslaw Adamczak.  A tail inequality for suprema of unbounded empirical processes with applications to markov chains.  _Electronic Journal of Probability_ , 13:1000–1034, 2008\.
  * Agarwal et al. [2010] Alekh Agarwal, Sahand Negahban, and Martin J Wainwright.  Fast global convergence rates of gradient methods for high-dimensional statistical recovery.  _Advances in Neural Information Processing Systems_ , 23, 2010.
  * Alquier et al. [2019] Pierre Alquier, Vincent Cottet, and Guillaume Lecué.  Estimation bounds and sharp oracle inequalities of regularized procedures with lipschitz loss functions.  _The Annals of Statistics_ , 47(4):2117–2144, 2019.
  * Annan et al. [2019] Dorcas A Annan, Nako Maishi, Tomoyoshi Soga, Randa Dawood, Cong Li, Hiroshi Kikuchi, Takayuki Hojo, Masahiro Morimoto, Tetsuya Kitamura, Mohammad Towfik Alam, et al.  Carbonic anhydrase 2 (caii) supports tumor blood endothelial cell survival under lactic acidosis in the tumor microenvironment.  _Cell Communication and Signaling_ , 17(1):1–15, 2019.
  * Bellec et al. [2018] Pierre C Bellec, Guillaume Lecué, and Alexandre B Tsybakov.  Slope meets lasso: improved oracle bounds and optimality.  _The Annals of Statistics_ , 46(6B):3603–3642, 2018.
  * Belloni et al. [2011] Alexandre Belloni, Victor Chernozhukov, and Lie Wang.  Square-root lasso: pivotal recovery of sparse signals via conic programming.  _Biometrika_ , 98(4):791–806, 2011.
  * Blumensath and Davies [2009] Thomas Blumensath and Mike E Davies.  Iterative hard thresholding for compressed sensing.  _Applied and computational harmonic analysis_ , 27(3):265–274, 2009.
  * Blumensath and Davies [2010] Thomas Blumensath and Mike E Davies.  Normalized iterative hard thresholding: Guaranteed stability and performance.  _IEEE Journal of selected topics in signal processing_ , 4(2):298–309, 2010.
  * Boyd et al. [2004] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  _Convex optimization_.  Cambridge university press, 2004.
  * Burer and Monteiro [2003] Samuel Burer and Renato DC Monteiro.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  _Mathematical Programming_ , 95(2):329–357, 2003\.
  * Cai et al. [2021a] Jian-Feng Cai, Jingyang Li, and Dong Xia.  Generalized low-rank plus sparse tensor estimation by fast riemannian optimization.  _arXiv preprint arXiv:2103.08895_ , 2021a.
  * Cai et al. [2021b] Jian-Feng Cai, Jingyang Li, and Dong Xia.  Provable tensor-train format tensor completion by riemannian optimization.  _arXiv preprint arXiv:2108.12163_ , 2021b.
  * Cai and Zhang [2015] T Tony Cai and Anru Zhang.  Rop: Matrix recovery via rank-one projections.  _The Annals of Statistics_ , 43(1):102–138, 2015\.
  * Cambier and Absil [2016] Léopold Cambier and P-A Absil.  Robust low-rank matrix completion by riemannian optimization.  _SIAM Journal on Scientific Computing_ , 38(5):S440–S460, 2016.
  * Candes and Plan [2011] Emmanuel J Candes and Yaniv Plan.  Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements.  _IEEE Transactions on Information Theory_ , 57(4):2342–2359, 2011.
  * Candès et al. [2011] Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright.  Robust principal component analysis?  _Journal of the ACM (JACM)_ , 58(3):1–37, 2011\.
  * Catoni [2012] Olivier Catoni.  Challenging the empirical mean and empirical variance: a deviation study.  In _Annales de l’IHP Probabilités et statistiques_ , volume 48, pages 1148–1185, 2012.
  * Chandrasekaran et al. [2012] Venkat Chandrasekaran, Benjamin Recht, Pablo A Parrilo, and Alan S Willsky.  The convex geometry of linear inverse problems.  _Foundations of Computational mathematics_ , 12(6):805–849, 2012.
  * Chang et al. [2005] Lin-Ching Chang, Derek K Jones, and Carlo Pierpaoli.  Restore: robust estimation of tensors by outlier rejection.  _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , 53(5):1088–1095, 2005.
  * Charisopoulos et al. [2021] Vasileios Charisopoulos, Yudong Chen, Damek Davis, Mateo D$\acute{i}$az, Lijun Ding, and Dmitriy Drusvyatskiy.  Low-rank matrix recovery with composite optimization: good conditioning and rapid convergence.  _Foundations of Computational Mathematics_ , 21(6):1505–1593, 2021.
  * Chen et al. [2011] Jianhui Chen, Jiayu Zhou, and Jieping Ye.  Integrating low-rank and group-sparse structures for robust multi-task learning.  In _Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 42–50, 2011.
  * Chen et al. [2016] Mengjie Chen, Chao Gao, and Zhao Ren.  A general decision theory for huber’s $\epsilon$-contamination model.  _Electronic Journal of Statistics_ , 10(2):3752–3774, 2016.
  * Chen et al. [2018] Mengjie Chen, Chao Gao, and Zhao Ren.  Robust covariance and scatter matrix estimation under huber’s contamination model.  _The Annals of Statistics_ , 46(5):1932–1960, 2018.
  * Chen et al. [2022] Xi Chen, Weidong Liu, and Xiaojun Mao.  Robust reduced rank regression in a distributed setting.  _Science China Mathematics_ , 65(8):1707–1730, 2022.
  * Chen and Wainwright [2015] Yudong Chen and Martin J. Wainwright.  Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees, 2015.
  * Chinot et al. [2020] Geoffrey Chinot, Guillaume Lecué, and Matthieu Lerasle.  Robust statistical learning with lipschitz and convex loss functions.  _Probability Theory and related fields_ , 176(3-4):897–940, 2020.
  * Chiu et al. [2021] Justin Chiu, Yuntian Deng, and Alexander Rush.  Low-rank constraints for fast inference in structured models.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Colas et al. [2011] Eva Colas, Cristina Perez, Silvia Cabrera, Nuria Pedrola, Marta Monge, Josep Castellvi, Fernando Eyzaguirre, Jesus Gregorio, Anna Ruiz, Marta Llaurado, et al.  Molecular markers of endometrial carcinoma detected in uterine aspirates.  _International journal of cancer_ , 129(10):2435–2444, 2011.
  * Davenport and Romberg [2016] Mark A Davenport and Justin Romberg.  An overview of low-rank matrix recovery from incomplete observations.  _IEEE Journal of Selected Topics in Signal Processing_ , 10(4):608–622, 2016.
  * Depersin [2020] Jules Depersin.  Robust subgaussian estimation with vc-dimension.  _arXiv preprint arXiv:2004.11734_ , 2020.
  * Depersin and Lecué [2022] Jules Depersin and Guillaume Lecué.  Robust sub-gaussian estimation of a mean vector in nearly linear time.  _The Annals of Statistics_ , 50(1):511–536, 2022\.
  * Diakonikolas and Kane [2019] Ilias Diakonikolas and Daniel M Kane.  Recent advances in algorithmic high-dimensional robust statistics.  _arXiv preprint arXiv:1911.05911_ , 2019.
  * Diakonikolas et al. [2019] Ilias Diakonikolas, Gautam Kamath, Daniel Kane, Jerry Li, Ankur Moitra, and Alistair Stewart.  Robust estimators in high-dimensions without the computational intractability.  _SIAM Journal on Computing_ , 48(2):742–864, 2019\.
  * Duchi et al. [2012] John C Duchi, Peter L Bartlett, and Martin J Wainwright.  Randomized smoothing for stochastic optimization.  _SIAM Journal on Optimization_ , 22(2):674–701, 2012.
  * Efron et al. [2004] Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani.  Least angle regression.  _The Annals of statistics_ , 32(2):407–499, 2004\.
  * Elsener and van de Geer [2018] Andreas Elsener and Sara van de Geer.  Robust low-rank matrix estimation.  _The Annals of Statistics_ , 46(6B):3481–3509, 2018.
  * Fan et al. [2021] Jianqing Fan, Weichen Wang, and Ziwei Zhu.  A shrinkage principle for heavy-tailed data: High-dimensional robust low-rank matrix recovery.  _Annals of statistics_ , 49(3):1239, 2021.
  * Foucart and Lecué [2017] Simon Foucart and Guillaume Lecué.  An iht algorithm for sparse recovery from subexponential measurements.  _IEEE Signal Processing Letters_ , 24(9):1280–1283, 2017.
  * Goffin [1977] Jean-Louis Goffin.  On convergence rates of subgradient optimization methods.  _Mathematical programming_ , 13(1):329–347, 1977\.
  * Grant and Boyd [2014] Michael Grant and Stephen Boyd.  Cvx: Matlab software for disciplined convex programming, version 2.1, 2014\.
  * Gross et al. [2010] David Gross, Yi-Kai Liu, Steven T Flammia, Stephen Becker, and Jens Eisert.  Quantum state tomography via compressed sensing.  _Physical review letters_ , 105(15):150401, 2010\.
  * Han et al. [2018] Wei Han, Fang Cao, Xiao-Jiao Gao, Hua-Bing Wang, Fang Chen, Si-Jian Cai, Cong Zhang, Yong-Wei Hu, Jun Ma, Xing Gu, et al.  Zic1 acts a tumor suppressor in breast cancer by targeting survivin.  _International journal of oncology_ , 53(3):937–948, 2018.
  * Hopkins [2020] Samuel B Hopkins.  Mean estimation with sub-gaussian rates in polynomial time.  2020\.
  * Huber [1965] Peter J Huber.  A robust version of the probability ratio test.  _The Annals of Mathematical Statistics_ , pages 1753–1758, 1965.
  * Huber and Ronchetti [2009] Peter J Huber and Elvezio M Ronchetti.  Robust statistics. hoboken.  _NJ: Wiley. doi_ , 10(1002):9780470434697, 2009\.
  * Jain et al. [2014] Prateek Jain, Ambuj Tewari, and Purushottam Kar.  On iterative hard thresholding methods for high-dimensional m-estimation.  _Advances in neural information processing systems_ , 27, 2014.
  * Kim et al. [2011] Ju Hee Kim, Han-Sung Kang, Tae Woo Kim, and Sun Jung Kim.  Differential methylation hybridization profiling identifies involvement of stat1-mediated pathways in breast cancer.  _International journal of oncology_ , 39(4):955–963, 2011.
  * Klopp et al. [2017] Olga Klopp, Karim Lounici, and Alexandre B Tsybakov.  Robust matrix completion.  _Probability Theory and Related Fields_ , 169(1):523–564, 2017.
  * Koenker and Bassett Jr [1978] Roger Koenker and Gilbert Bassett Jr.  Regression quantiles.  _Econometrica: journal of the Econometric Society_ , pages 33–50, 1978.
  * Koenker and Hallock [2001] Roger Koenker and Kevin F Hallock.  Quantile regression.  _Journal of economic perspectives_ , 15(4):143–156, 2001.
  * Koltchinskii et al. [2011] Vladimir Koltchinskii, Karim Lounici, and Alexandre B Tsybakov.  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  _The Annals of Statistics_ , 39(5):2302–2329, 2011.
  * Lecué and Lerasle [2019] Guillaume Lecué and Matthieu Lerasle.  Learning from mom’s principles: Le cam’s approach.  _Stochastic Processes and their applications_ , 129(11):4385–4410, 2019.
  * Lecué and Mendelson [2017] Guillaume Lecué and Shahar Mendelson.  Sparse recovery under weak moment assumptions.  _Journal of the European Mathematical Society_ , 19(3):881–904, 2017.
  * Li et al. [2012] Bingyan Li, Gengyun Wen, Yongliang Zhao, Jian Tong, and Tom K Hei.  The role of tgfbi in mesothelioma and breast cancer: association with tumor suppression.  _BMC cancer_ , 12(1):1–12, 2012.
  * Li et al. [2020] Xiao Li, Zhihui Zhu, Anthony Man-Cho So, and Rene Vidal.  Nonconvex robust low-rank matrix recovery.  _SIAM Journal on Optimization_ , 30(1):660–686, 2020.
  * Loh and Wainwright [2011] Po-Ling Loh and Martin J Wainwright.  High-dimensional regression with noisy and missing data: Provable guarantees with non-convexity.  _Advances in neural information processing systems_ , 24, 2011.
  * Lounici et al. [2011] Karim Lounici, Massimiliano Pontil, Sara Van De Geer, and Alexandre B Tsybakov.  Oracle inequalities and optimal inference under group sparsity.  2011\.
  * Ludoux and Talagrand [1991] M Ludoux and M Talagrand.  Probability in banach spaces: isoperimetry and processes, 1991.
  * Lugosi and Mendelson [2019a] Gábor Lugosi and Shahar Mendelson.  Mean estimation and regression under heavy-tailed distributions: A survey.  _Foundations of Computational Mathematics_ , 19(5):1145–1190, 2019a.
  * Lugosi and Mendelson [2019b] Gábor Lugosi and Shahar Mendelson.  Regularization, sparse recovery, and median-of-means tournaments.  _Bernoulli_ , 25(3):2075–2106, 2019b.
  * Lugosi and Mendelson [2019c] Gabor Lugosi and Shahar Mendelson.  Risk minimization by median-of-means tournaments.  _Journal of the European Mathematical Society_ , 22(3):925–965, 2019c.
  * Luo and Zhang [2022] Yuetian Luo and Anru R Zhang.  Tensor-on-tensor regression: Riemannian optimization, over-parameterization, statistical-computational gap, and their interplay.  _arXiv preprint arXiv:2206.08756_ , 2022.
  * Lyu et al. [2021] Zhongyuan Lyu, Dong Xia, and Yuan Zhang.  Latent space model for higher-order networks and generalized tensor decomposition.  _arXiv preprint arXiv:2106.16042_ , 2021.
  * Ma and Fattahi [2022] Jianhao Ma and Salar Fattahi.  Global convergence of sub-gradient method for robust matrix recovery: Small initialization, noisy measurements, and over-parameterization.  _arXiv preprint arXiv:2202.08788_ , 2022.
  * Ma and Wu [2015] Zongming Ma and Yihong Wu.  Volume ratio, sparsity, and minimaxity under unitarily invariant norms.  _IEEE Transactions on Information Theory_ , 61(12):6939–6956, 2015.
  * Minsker [2015] Stanislav Minsker.  Geometric median and robust estimation in banach spaces.  _Bernoulli_ , 21(4):2308–2335, 2015.
  * Minsker [2018] Stanislav Minsker.  Sub-gaussian estimators of the mean of a random matrix with heavy-tailed entries.  _The Annals of Statistics_ , 46(6A):2871–2903, 2018.
  * Minsker et al. [2022] Stanislav Minsker, Mohamed Ndaoud, and Lang Wang.  Robust and tuning-free sparse linear regression via square-root slope.  _arXiv preprint arXiv:2210.16808_ , 2022.
  * Mishra et al. [2014] Bamdev Mishra, Gilles Meyer, Silvere Bonnabel, and Rodolphe Sepulchre.  Fixed-rank matrix factorizations and riemannian low-rank optimization.  _Computational Statistics_ , 29(3):591–621, 2014\.
  * Nakata et al. [2004] B Nakata, T Takashima, Y Ogawa, T Ishikawa, and K Hirakawa.  Serum cyfra 21-1 (cytokeratin-19 fragments) is a useful tumour marker for detecting disease relapse and assessing treatment efficacy in breast cancer.  _British journal of cancer_ , 91(5):873–878, 2004\.
  * Needell and Tropp [2009] Deanna Needell and Joel A Tropp.  Cosamp: Iterative signal recovery from incomplete and inaccurate samples.  _Applied and computational harmonic analysis_ , 26(3):301–321, 2009.
  * Negahban and Wainwright [2011] Sahand Negahban and Martin J Wainwright.  Estimation of (near) low-rank matrices with noise and high-dimensional scaling.  _The Annals of Statistics_ , pages 1069–1097, 2011.
  * Nishizuka et al. [2003] Satoshi Nishizuka, Lu Charboneau, Lynn Young, Sylvia Major, William C Reinhold, Mark Waltham, Hosein Kouros-Mehr, Kimberly J Bussey, Jae K Lee, Virginia Espina, et al.  Proteomic profiling of the nci-60 cancer cell lines using new high-density reverse-phase lysate microarrays.  _Proceedings of the National Academy of Sciences_ , 100(24):14229–14234, 2003.
  * Niu et al. [2019] Jie Niu, Xiao-Meng Li, Xiao Wang, Chao Liang, Yi-Dan Zhang, Hai-Ying Li, Fan-Ye Liu, Hua Sun, Song-Qiang Xie, and Dong Fang.  Dkk1 inhibits breast cancer cell migration and invasion through suppression of $\beta$-catenin/mmp7 signaling pathway.  _Cancer cell international_ , 19(1):1–13, 2019\.
  * Osta et al. [2004] Walid A Osta, Yian Chen, Kaidi Mikhitarian, Michael Mitas, Mohamed Salem, Yusuf A Hannun, David J Cole, and William E Gillanders.  Epcam is overexpressed in breast cancer and is a potential target for breast cancer gene therapy.  _Cancer research_ , 64(16):5818–5824, 2004.
  * Pan et al. [2021] Xiaoou Pan, Qiang Sun, and Wen-Xin Zhou.  Iteratively reweighted $\ell_{1}$-penalized robust regression.  _Electronic Journal of Statistics_ , 15(1):3287–3348, 2021.
  * Raskutti et al. [2011] Garvesh Raskutti, Martin J Wainwright, and Bin Yu.  Minimax rates of estimation for high-dimensional linear regression over $\ell_{q}$-balls.  _IEEE transactions on information theory_ , 57(10):6976–6994, 2011.
  * Rigollet and Tsybakov [2011] Philippe Rigollet and Alexandre Tsybakov.  Exponential screening and optimal rates of sparse estimation.  2011\.
  * Rohde and Tsybakov [2011] Angelika Rohde and Alexandre B Tsybakov.  Estimation of high-dimensional low-rank matrices.  _The Annals of Statistics_ , 39(2):887–930, 2011\.
  * Shankavaram et al. [2007] Uma T Shankavaram, William C Reinhold, Satoshi Nishizuka, Sylvia Major, Daisaku Morita, Krishna K Chary, Mark A Reimers, Uwe Scherf, Ari Kahn, Douglas Dolginow, et al.  Transcript and protein expression profiles of the nci-60 cancer cell panel: an integromic microarray study.  _Molecular cancer therapeutics_ , 6(3):820–832, 2007.
  * Shen and Li [2017] Jie Shen and Ping Li.  A tight bound of hard thresholding.  _The Journal of Machine Learning Research_ , 18(1):7650–7691, 2017.
  * Siddiqi et al. [2010] Sajid Siddiqi, Byron Boots, and Geoffrey Gordon.  Reduced-rank hidden markov models.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , pages 741–748. JMLR Workshop and Conference Proceedings, 2010.
  * Sun et al. [2020] Qiang Sun, Wen-Xin Zhou, and Jianqing Fan.  Adaptive huber regression.  _Journal of the American Statistical Association_ , 115(529):254–265, 2020.
  * Sun et al. [2017] Will Wei Sun, Junwei Lu, Han Liu, and Guang Cheng.  Provable sparse tensor decomposition.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(3):899–916, 2017.
  * Tan et al. [2022] Kean Ming Tan, Lan Wang, and Wen-Xin Zhou.  High-dimensional quantile regression: Convolution smoothing and concave regularization.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(1):205–233, 2022.
  * Thompson [2020] Philip Thompson.  Outlier-robust sparse/low-rank least-squares regression and robust matrix completion.  _arXiv preprint arXiv:2012.06750_ , 2020.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tong et al. [2021a] Tian Tong, Cong Ma, and Yuejie Chi.  Accelerating ill-conditioned low-rank matrix estimation via scaled gradient descent, 2021a.
  * Tong et al. [2021b] Tian Tong, Cong Ma, and Yuejie Chi.  Low-rank matrix recovery with scaled subgradient methods: Fast and robust convergence without the condition number.  _IEEE Transactions on Signal Processing_ , 69:2396–2409, 2021b.
  * Van de Geer and van de Geer [2000] Sara A Van de Geer and Sara van de Geer.  _Empirical Processes in M-estimation_ , volume 6.  Cambridge university press, 2000.
  * Van Der Vaart et al. [1996] Aad W Van Der Vaart, Aad van der Vaart, Adrianus Willem van der Vaart, and Jon Wellner.  _Weak convergence and empirical processes: with applications to statistics_.  Springer Science & Business Media, 1996.
  * Vandereycken [2013] Bart Vandereycken.  Low-rank matrix completion by riemannian optimization.  _SIAM Journal on Optimization_ , 23(2):1214–1236, 2013.
  * Vershynin [2018] Roman Vershynin.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * Wang et al. [2012] Lan Wang, Yichao Wu, and Runze Li.  Quantile regression for analyzing heterogeneity in ultra-high dimension.  _Journal of the American Statistical Association_ , 107(497):214–222, 2012.
  * Wang et al. [2020] Lan Wang, Bo Peng, Jelena Bradic, Runze Li, and Yunan Wu.  A tuning-free robust and efficient approach to high-dimensional regression.  _Journal of the American Statistical Association_ , 115(532):1700–1714, 2020.
  * Wang et al. [2021] Lili Wang, Chao Zheng, Wen Zhou, and Wen-Xin Zhou.  A new principle for tuning-free huber regression.  _Statistica Sinica_ , 31(4):2153–2177, 2021.
  * Wei et al. [2016] Ke Wei, Jian-Feng Cai, Tony F Chan, and Shingyu Leung.  Guarantees of riemannian optimization for low rank matrix recovery.  _SIAM Journal on Matrix Analysis and Applications_ , 37(3):1198–1222, 2016.
  * Welsh [1989] AH Welsh.  On m-processes and m-estimation.  _The Annals of Statistics_ , pages 337–361, 1989.
  * Xia [2014] Dong Xia.  Optimal schatten-q and ky-fan-k norm rate of low rank matrix estimation.  _arXiv preprint arXiv:1403.6499_ , 2014.
  * Xia [2021] Dong Xia.  Normal approximation and confidence region of singular subspaces.  _Electronic Journal of Statistics_ , 15(2):3798–3851, 2021.
  * Xia and Koltchinskii [2016] Dong Xia and Vladimir Koltchinskii.  Estimation of low rank density matrices: bounds in schatten norms and other distances.  _Electronic Journal of Statistics_ , 10(2):2717–2745, 2016.
  * Xia and Yuan [2021] Dong Xia and Ming Yuan.  Statistical inferences of linear forms for noisy matrix completion.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 83(1):58–77, 2021.
  * Xia et al. [2021] Dong Xia, Ming Yuan, and Cun-Hui Zhang.  Statistically optimal and computationally efficient low rank tensor completion from noisy entries.  _The Annals of Statistics_ , 49(1):76–99, 2021\.
  * Xie et al. [2019] Wan-Ying Xie, Ruo-Hui He, Jun Zhang, Yi-Jing He, Zan Wan, Cheng-Fang Zhou, Yong-Jun Tang, Zhi Li, Howard L Mcleod, and Jie Liu.  $\beta$-blockers inhibit the viability of breast cancer cells by regulating the erk/cox-2 signaling pathway and the drug response is affected by adrb2 single-nucleotide polymorphisms.  _Oncology Reports_ , 41(1):341–350, 2019.
  * Xue et al. [2022] Jinqi Xue, Xudong Zhu, Xinbo Qiao, Yulun Wang, Jiawen Bu, Xiaoying Zhang, Qingtian Ma, Lu Liang, Lisha Sun, and Caigang Liu.  Cemip as a potential biomarker and therapeutic target for breast cancer patients.  _International Journal of Medical Sciences_ , 19(3):434, 2022.
  * Yang et al. [2012] Linlin Yang, Yuan Chen, Tiantian Cui, Thomas Knösel, Qing Zhang, Kai Frederik Albring, Otmar Huber, and Iver Petersen.  Desmoplakin acts as a tumor suppressor by inhibition of the wnt/$\beta$-catenin signaling pathway in human lung cancer.  _Carcinogenesis_ , 33(10):1863–1870, 2012.
  * Ye and Zhang [2010] Fei Ye and Cun-Hui Zhang.  Rate minimaxity of the lasso and dantzig selector for the lq loss in lr balls.  _The Journal of Machine Learning Research_ , 11:3519–3540, 2010.
  * Ye et al. [2022] Xinxing Ye, Jie Zhou, Dandan Tong, Dandan Wang, Hui Wang, Jixue Guo, and Xinmei Kang.  E2f1 affects the therapeutic response to neoadjuvant therapy in breast cancer.  _Disease Markers_ , 2022, 2022.
  * Yu and Chai [2017] Wenzhen Yu and Hong Chai.  Inhibition of bambi reduces the viability and motility of colon cancer via activating tgf-$\beta$/smad pathway in vitro and in vivo retraction in/10.3892/ol. 2021.12505.  _Oncology Letters_ , 14(4):4793–4799, 2017.
  * Zhang et al. [2017] Shizhen Zhang, Zhen Wang, Weiwei Liu, Rui Lei, Jinlan Shan, Ling Li, and Xiaochen Wang.  Distinct prognostic values of s100 mrna expression in breast cancer.  _Scientific reports_ , 7(1):1–11, 2017.
  * Zhang and Xia [2020] Yuan Zhang and Dong Xia.  Edgeworth expansions for network moments.  _arXiv preprint arXiv:2004.06615_ , 2020.
  * Zhao et al. [2022] Peng Zhao, Yun Yang, and Qiao-Chu He.  High-dimensional linear regression via implicit regularization.  _Biometrika_ , 2022.
  * Zhao et al. [2015] Tuo Zhao, Zhaoran Wang, and Han Liu.  A nonconvex optimization framework for low rank matrix estimation.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Zheng et al. [2015] Qi Zheng, Limin Peng, and Xuming He.  Globally adaptive quantile regression with ultra-high dimensional data.  _Annals of statistics_ , 43(5):2225, 2015.
  * Zheng and Lafferty [2016] Qinqing Zheng and John Lafferty.  A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements, 2016.
  * Zhou et al. [2017] Tao Zhou, Yong Li, Li Yang, Liang Liu, Yingchao Ju, and Chunxiao Li.  Silencing of anxa3 expression by rna interference inhibits the proliferation and invasion of breast cancer cells.  _Oncology Reports_ , 37(1):388–398, 2017.
  * Zhu et al. [2020] Junxian Zhu, Canhong Wen, Jin Zhu, Heping Zhang, and Xueqin Wang.  A polynomial algorithm for best-subset selection problem.  _Proceedings of the National Academy of Sciences_ , 117(52):33117–33123, 2020.
  * Zhuo et al. [2021] Jiacheng Zhuo, Jeongyeol Kwon, Nhat Ho, and Constantine Caramanis.  On the computational and statistical complexity of over-parameterized matrix sensing.  _arXiv preprint arXiv:2102.02756_ , 2021.

Appendix of “Computationally Efficient and Statistically Optimal Robust High-
Dimensional Linear Regression”

Yinan Shen, Jingyang Li, Jian-Feng Cai and Dong Xia

Hong Kong University of Science and Technology

Section A presents materials that are not shown in the main text due to
limited length, including Riemannian optimization details and simulation
results. In Section B, we present proofs of sparse vector recovery. Section C
proves Proposition 1, convergence dynamics of low-rank regression under two-
phase regularity conditions. Section D contains four sub-sections, D.1-D.4,
which verify the regularity conditions in certain different cases, absolute
loss, Huber loss and quantile loss. The initialization analysis is provided in
Section E. Then Section F exhibits the technical lemmas. Finally, lengthy
proofs of two significant lemmas 14 matrix perturbation theories and
Proposition 3 empirical processes for low-rank regression are put in Section G
and Section H respectively. In the last section I, we analyze outlier
robustness.

The following key ingredients act as a guidance for the proof of low-rank
linear regression.

  1. 1.

Tangent Space over ${\mathbb{M}}_{r}$. In the algorithm, it updates in the way
of
${\mathbf{T}}_{l+1}=\text{SVD}_{r}({\mathbf{T}}_{l}-\eta_{l}{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l}))$.
Then in the convergence analyses, it involves the upper bound of
$\|{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})\|_{\rm{F}}^{2}$ and
$\left|\langle{\mathcal{P}}_{{\mathbb{T}}_{l}}^{\perp}{\mathbf{G}}_{l},{\mathbf{T}}_{l}-{\mathbf{T}}^{*}\rangle\right|$.
As for the first term, note that

$\displaystyle{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})={\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}+{\mathbf{G}}_{l}{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}={\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}+({\mathbf{I}}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}){\mathbf{G}}_{l}{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}$

could be written into sum of two rank at most $r$ matrices. Then we use
‘Partial Frobenius Norm’ (see Tong et al. [2021a] or Lemma 12) to bound
$\|{\mathcal{P}}_{{\mathbb{T}}_{l}}{\mathbf{G}}_{l}\|_{\rm{F}}$; see Lemma 12
and Lemma 13 for more details. As for the second term, we use

$\displaystyle\left|\langle{\mathcal{P}}_{{\mathbb{T}}_{l}}^{\perp}{\mathbf{G}}_{l},{\mathbf{T}}_{l}-{\mathbf{T}}^{*}\rangle\right|=\left|\langle{\mathbf{G}}_{l},{\mathcal{P}}_{{\mathbb{T}}_{l}}^{\perp}({\mathbf{T}}_{l}-{\mathbf{T}}^{*})\rangle\right|,$

and use
${\mathcal{P}}_{{\mathbb{T}}_{l}}^{\perp}({\mathbf{T}}_{l}-{\mathbf{T}}^{*})=-({\mathbf{I}}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}){\mathbf{T}}^{*}({\mathbf{V}}{\mathbf{V}}^{\top}-{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top})=({\mathbf{I}}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top})({\mathbf{T}}_{l}-{\mathbf{T}}^{*})({\mathbf{V}}{\mathbf{V}}^{\top}-{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top})$.
Hence, with definition of partial Frobenius norm Lemma 12, one has

$\displaystyle\left|\langle{\mathcal{P}}_{{\mathbb{T}}_{l}}^{\perp}{\mathbf{G}}_{l},{\mathbf{T}}_{l}-{\mathbf{T}}^{*}\rangle\right|\leq\|({\mathbf{I}}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top})({\mathbf{T}}_{l}-{\mathbf{T}}^{*})({\mathbf{V}}{\mathbf{V}}^{\top}-{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top})\|_{\rm{F}}\|{\mathbf{G}}_{l}\|_{\rm{F},r}.$

  2. 2.

Perturbation lemma. Based on Theorem 1 in work Xia [2021], we develop powerful
analysis scheme for matrix perturbations, where we have
$\|{\mathbf{T}}_{l+1}-{\mathbf{T}}^{*}\|_{\rm{F}}\leq\|{\mathbf{T}}_{l}-{\mathbf{T}}^{*}-\eta_{l}{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})\|_{\rm{F}}+O(\|{\mathbf{T}}_{l}-{\mathbf{T}}^{*}-\eta_{l}{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})\|_{\rm{F}}^{2})$.
For clarity, we present Theorem 1 in work Xia [2021] as Theorem 11. Lemma 18
presents symmetric matrix perturbations, based on which we have general matrix
perturbations Lemma 14.

  3. 3.

Empirical Processes. We use emprical process theories [Ludoux and Talagrand,
1991, Van Der Vaart et al., 1996, Adamczak, 2008] to prove the regularity
properties. One of the merits is that it enables strong exponential
probability.

  4. 4.

Initialization Guarantees. Note that the perturbation bound and convergence
dynamics Equation (9) require
$\|{\mathbf{M}}_{0}-{\mathbf{M}}^{*}\|_{\rm{F}}\lesssim\sigma_{r}$. Thus we
start with a warm initialization. A clean and simple spectral initialization
scheme is provided. Section E shows the proofs.

  5. 5.

Proposition 1. Proposition 1 analyzes the general convergence under two-phase
regularity property assumption and depicts the stepsize selections. Lemma 4,
6, 8 and 10 prove two-phase regularity properties under several settings.
Namely, the loss function behaves in the typical non-smooth way when it is far
from the oracle matrix ${\mathbf{M}}^{*}$, while when it gets sufficiently
close to ${\mathbf{M}}^{*}$, the loss function performs alike smooth
$\ell_{2}$ loss.

##  Appendix A Additional Results

This section collects some results which are not presented in the main context
due to the page limit.

###  A.1 Computational benefit of Riemannian sub-gradient algorithm

The Riemannian sub-gradient is fast computable. Let
${\mathbf{M}}_{l}={\mathbf{U}}_{l}{\boldsymbol{\Sigma}}_{l}{\mathbf{V}}_{l}^{\top}$
be the thin SVD of ${\mathbf{M}}_{l}$. It is well-known [Absil et al., 2009,
Vandereycken, 2013] that the tangent space ${\mathbb{T}}_{l}$ can be
characterized by
${\mathbb{T}}_{l}:=\\{{\mathbf{Z}}\in{\mathbb{R}}^{d_{1}\times
d_{2}}:{\mathbf{Z}}={\mathbf{U}}_{l}{\mathbf{R}}^{\top}+{\mathbf{L}}{\mathbf{V}}_{l}^{\top},{\mathbf{R}}\in{\mathbb{R}}^{d_{2}\times
r},{\mathbf{L}}\in{\mathbb{R}}^{d_{1}\times r}\\}$. Then projection of
${\mathbf{G}}_{l}$ onto ${\mathbb{T}}_{l}$ is

${\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})={\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}+{\mathbf{G}}_{l}{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top},$

which is of rank at most $2r$. The first step in Algorithm 2 could be written
as

$\displaystyle{~{}~{}~{}~{}}{\mathbf{M}}_{l}-\eta_{l}{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})$
$\displaystyle={\mathbf{U}}_{l}\left({\boldsymbol{\Sigma}}_{l}-\eta_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}{\mathbf{V}}_{l}\right){\mathbf{V}}_{l}^{\top}-\eta_{l}{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}\left({\mathbf{I}}_{d_{2}}-{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}\right)-\eta_{l}\left({\mathbf{I}}_{d_{1}}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}\right){\mathbf{G}}_{l}{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}.$

Let
$\left({\mathbf{I}}_{d_{2}}-{\mathbf{V}}_{l}{\mathbf{V}}_{l}^{\top}\right){\mathbf{G}}_{l}^{\top}{\mathbf{U}}_{l}={\mathbf{Q}}_{1}{\mathbf{R}}_{1}$
and
$\left({\mathbf{I}}_{d_{1}}-{\mathbf{U}}_{l}{\mathbf{U}}_{l}^{\top}\right){\mathbf{G}}_{l}{\mathbf{V}}_{l}={\mathbf{Q}}_{2}{\mathbf{R}}_{2}$
be the QR decomposition. Then the update is

$\displaystyle{\mathbf{M}}_{l}-\eta_{l}{\mathcal{P}}_{{\mathbb{T}}_{l}}({\mathbf{G}}_{l})=\left[{\mathbf{U}}_{l}\;{\mathbf{Q}}_{2}\right]\left[\begin{matrix}{\boldsymbol{\Sigma}}_{l}-\eta_{l}{\mathbf{U}}_{l}^{\top}{\mathbf{G}}_{l}{\mathbf{V}}_{l}&-\eta_{l}{\mathbf{R}}_{1}^{\top}\\\
-\eta_{l}{\mathbf{R}}_{2}&\boldsymbol{0}\end{matrix}\right]\left[{\mathbf{V}}_{l}\;{\mathbf{Q}}_{1}\right]^{\top},$

with $\left[{\mathbf{U}}_{l}\;{\mathbf{Q}}_{2}\right]$,
$\left[{\mathbf{V}}_{l}\;{\mathbf{Q}}_{1}\right]$ being orthogonal matrices.
Consequently, the final step of retraction only requires the SVD of the medium
term, a $2r\times 2r$ matrix, which can be computed using $O(r^{3})$ floating
point operations. See, e.g., Vandereycken [2013] and Mishra et al. [2014] for
more details.

###  A.2 More simulation results of low-rank linear regression

##### Accuracy

We now compare the statistical accuracy of the final estimator output by
RsGrad, ScaledSM and those by convex approaches including the nuclear norm
penalized [Elsener and van de Geer, 2018] absolute loss (convex $\ell_{1}$),
Huber loss (convex Huber) and the square loss ($\ell_{2}$-loss) by seminal
work Candes and Plan [2011]. Functions from Matlab library cvx [Grant and
Boyd, 2014] are borrowed to implement the convex methods for convex
$\ell_{1}$, convex Huber and $\ell_{2}$-loss. For each setting, every method
is repeated for $10$ times and the box-plots of error rates are presented.

The parameter settings are similar as above where both the Gaussian noise and
Student’s $t$-distribution are experimented. The comparison of statistical
accuracy under Gaussian noise is displayed in Figure 5. The top two plots
suggest that all the convex methods and ScaledSM performs poorly when sample
size is small ($n=2.5rd_{1}$). On the other hand, ScaledSM becomes much better
when sample size is large but still under-performs RsGrad, as predicted by our
theory. Similar results are also observed when noise has a Students’
$t$-distribution with d.f. $\nu=2$. See Figure 6.

(a) SNR=40

(b) SNR=80

(c) SNR=40

(d) SNR=80

Figure 5: Accuracy comparisons of RsGrad, ScaledSM, robust convex models
[Elsener and van de Geer, 2018] and convex $\ell_{2}$ loss [Candes and Plan,
2011] under Gaussian noise. Top two figures: $n=2.5rd_{1}$; bottom two
figures: $n=5rd_{1}$.

(a) SNR=40

(b) SNR=80

(c) SNR=40

(d) SNR=80

Figure 6: Accuracy comparisons of RsGrad, ScaledSM, robust convex models
[Elsener and van de Geer, 2018] and convex $\ell_{2}$ loss [Candes and Plan,
2011] under Students’ $t$-distribution noise with d.f. $\nu=2$. Top two
figures: $n=2.5rd_{1}$; bottom two figures: $n=5rd_{1}$.

From Figure 5 and Figure 6, we observe that non-convex models achieve better
accuracy, which may be attributed to more finely tuning parameters than convex
approaches. Convex ones we compare with are also proven statistically optimal
in theory and are computationally more demanding. For instance, in the case
$n=5rd_{1}$, convex $\ell_{1}$, convex Huber, convex $\ell_{2}$ need more than
1.5h, 1.5h, 3h to complete the ten repetitions, respectively. In sharp
contrast, RsGrad and ScaledSM only take around ten minutes. When it comes to
an even larger sample size like $n=10d_{1}r$, the Matlab package on our
computing platform needs 3h and 6h to complete the simulations (10
repetitions) for convex $\ell_{1}$ and convex Huber, respectively, whereas
RsGrad and ScaledSM each takes about 20 minutes to do the job.

##  Appendix B Proofs of Sparse Linear Regression

In general, there are three key technical points, Proposition 2, Theorem 9 and
Lemma 2; see below.

######  Proposition 2.

Suppose function $f(\cdot)$ is given in Equation (3). Then there exist
absolute constants $C,C_{1}>0$ such that

$\displaystyle\left|f({\boldsymbol{\beta}}+\Delta{\boldsymbol{\beta}})-f({\boldsymbol{\beta}})-{\mathbb{E}}\left[f({\boldsymbol{\beta}}+\Delta{\boldsymbol{\beta}})-f({\boldsymbol{\beta}})\right]\right|\leq
C\sqrt{nsC_{{\tiny u}}\log(2d/s)}\|\Delta{\boldsymbol{\beta}}\|_{2}$

holds for all ${\boldsymbol{\beta}}\in{\mathbb{R}}^{d}$ and
$\Delta{\boldsymbol{\beta}}\in\left\\{{\mathbf{b}}\in{\mathbb{R}}^{d}:\;\left|\text{supp}({\mathbf{b}})\right|\leq
s\right\\}$ with probability exceeding
$1-\exp(-\frac{C_{1}^{2}s\log(2d/s)}{3})-3\exp\left(-\frac{\sqrt{n\log(2d/s)}}{\log
n}\right)$.

Proposition 2 bounds deviation of
$f({\boldsymbol{\beta}}+\Delta{\boldsymbol{\beta}})-f({\boldsymbol{\beta}})$.
Its proof is based on empirical process theories and is presented later in
this section. Besides, we also need perturbation bound. Rough perturbation
bound $\|{\mathcal{H}}_{K}({\mathbf{b}})-{\mathbf{x}}\|_{2}\leq
2\|{\mathbf{b}}-{\mathbf{x}}\|_{2}$ is usually used in $\ell_{2}$ loss
[Blumensath and Davies, 2009, 2010] and is good enough to derive statiscal
optimality in $\ell_{2}$ loss. However, the coefficient $2$ is too large for
phase two convergence of $\ell_{1}$ loss. Thanksfully, work Shen and Li [2017]
provides a tight bound for hard thresholding perturbations.

######  Theorem 9 (Tight Bound for Hard Thresholding Shen and Li [2017]).

Let ${\mathbf{b}}\in{\mathbb{R}}^{d}$ be an arbitrary vector and
${\mathbf{x}}\in{\mathbb{R}}^{d}$ be any $K$-sparse signal. For any $k\geq K$,
we have the following bound:

$\left\|\mathcal{H}_{k}(\boldsymbol{b})-\boldsymbol{x}\right\|_{2}\leq\sqrt{\nu}\|\boldsymbol{b}-\boldsymbol{x}\|_{2},\quad\nu=1+\frac{\rho+\sqrt{(4+\rho)\rho}}{2},\quad\rho=\frac{\min\\{K,d-k\\}}{k-K+\min\\{K,d-k\\}}.$

Note that to some extent larger hard thresholding parameter $k$ infers smaller
upper bound coefficient $\nu$. Besides, Theorem 9 works for arbitrary
${\mathbf{b}}$ and puts no assumption on ${\mathbf{b}}$. It leads to we do not
require good initializations for sparse vector recovery but need the practical
thresholding parameter $\tilde{s}\geq\min\\{cs,d\\}$ to have convergence
guarantees. We conjecture there are tighter bounds if with small enough
$\|{\mathbf{b}}-{\mathbf{x}}\|_{2}$ based on which we could choose
$\tilde{s}=s$. We leave this for future study.

######  Lemma 2.

Suppose function $f(\cdot)$ is given in Equation (3). Suppose Assumptions 1
and 2 hold. Then for any vector
${\boldsymbol{\beta}}\in\left\\{{\mathbf{b}}\in{\mathbb{R}}^{d}:\,\|{\mathbf{b}}-{\boldsymbol{\beta}}^{*}\|_{2}\leq
8\sqrt{C_{{\tiny u}}C_{{\tiny l}}^{-1}}\gamma\right\\}$, we have

$\displaystyle{\mathbb{E}}f({\boldsymbol{\beta}})-{\mathbb{E}}f({\boldsymbol{\beta}}^{*})\geq\frac{n}{6b_{0}}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}^{2}.$

For arbitrary
${\boldsymbol{\beta}}_{1},{\boldsymbol{\beta}}_{2}\in{\mathbb{R}}^{d}$, we
have

$\displaystyle{\mathbb{E}}f({\boldsymbol{\beta}}_{1})-{\mathbb{E}}f({\boldsymbol{\beta}}_{2})\leq
nb_{1}^{-1}C_{{\tiny
u}}\cdot\left[\|{\boldsymbol{\beta}}_{1}-{\boldsymbol{\beta}}_{2}\|_{2}^{2}+2\|{\boldsymbol{\beta}}_{1}-{\boldsymbol{\beta}}_{2}\|_{2}\|{\boldsymbol{\beta}}_{2}-{\boldsymbol{\beta}}^{*}\|_{2}\right].$

Lemma 2 provides upper bound and lower bound of expectation of
$f({\boldsymbol{\beta}}_{1})-f({\boldsymbol{\beta}}_{2})$ and
$f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})$, respectively. Then we
are ready to prove regularity proporties Lemma 1.

###### Proof of Lemma 1.

We shall assume the event

${\boldsymbol{\mathcal{E}}}:=\left\\{\sup_{{\boldsymbol{\beta}}\in{\mathbb{R}}^{d},\;|\text{supp}(\Delta{\boldsymbol{\beta}})|\leq
3\tilde{s}}\left|f({\boldsymbol{\beta}}+\Delta{\boldsymbol{\beta}})-f({\boldsymbol{\beta}})-{\mathbb{E}}\left[f({\boldsymbol{\beta}}+\Delta{\boldsymbol{\beta}})-f({\boldsymbol{\beta}})\right]\right|\cdot\|\Delta{\boldsymbol{\beta}}\|_{2}^{-1}\leq
C\sqrt{n\tilde{s}C_{{\tiny u}}\log(2d/\tilde{s})}\right\\}$

holds. Proposition 2 proves ${\mathbb{P}}({\boldsymbol{\mathcal{E}}})\geq
1-\exp(-\frac{C^{2}\tilde{s}\log(2d/\tilde{s})}{3})-3\exp\left(-\frac{\sqrt{n\log(2d/\tilde{s})}}{\log
n}\right)$. Based on event ${\boldsymbol{\mathcal{E}}}$, we shall obtain the
regularity properties.

First, analyze phase one properties, namely, when
$\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\geq 8\sqrt{C_{{\tiny
l}}^{-1}}\gamma$. As for
$f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})$, under the event
${\boldsymbol{\mathcal{E}}}$, it has

$\displaystyle f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})$
$\displaystyle\geq{\mathbb{E}}f({\boldsymbol{\beta}})-{\mathbb{E}}f({\boldsymbol{\beta}}^{*})-C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}$
$\displaystyle={\mathbb{E}}\left[\|{\mathbf{X}}({\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*})-{\boldsymbol{\xi}}\|_{1}-\|{\boldsymbol{\xi}}\|_{1}\right]-C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}$
$\displaystyle\geq{\mathbb{E}}\|{\mathbf{X}}({\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*})\|_{1}-2{\mathbb{E}}\|{\boldsymbol{\xi}}\|_{1}-C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}.$

Notice that
$\left|{\mathbf{X}}_{i}^{\top}({\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*})\right|$
follows folded Gaussian distributions and then we have
${\mathbb{E}}\|{\mathbf{X}}({\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*})\|_{1}\geq
n\sqrt{\frac{2}{\pi}}\sqrt{C_{{\tiny
l}}}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}$. Thus, together
with ${\mathbb{E}}|\xi|=\gamma$,
$\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\geq 8\sqrt{C_{{\tiny
l}}^{-1}}\gamma$ and $n\geq C_{1}C_{{\tiny u}}C_{{\tiny
l}}^{-1}\tilde{s}\log(2d/\tilde{s})$, it has

$\displaystyle
f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})\geq\frac{n}{4}\sqrt{C_{{\tiny
l}}}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}.$

Then consider
$\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}$. By
definition of sub-gradient, it has

$\displaystyle\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}^{2}=\langle{\mathbf{G}},{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\rangle\leq
f({\boldsymbol{\beta}}+{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}}))-f({\boldsymbol{\beta}}).$
(5)

On the other hand, under the event ${\boldsymbol{\mathcal{E}}}$, it has

$\displaystyle{~{}~{}~{}~{}}f({\boldsymbol{\beta}}+{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}}))-f({\boldsymbol{\beta}})$
$\displaystyle\leq{\mathbb{E}}f({\boldsymbol{\beta}}+{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}}))-{\mathbb{E}}f({\boldsymbol{\beta}})+C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}$
$\displaystyle\leq{\mathbb{E}}\|{\mathbf{X}}{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{1}+C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}.$

Then together with folded Gaussian distributions, it has
${\mathbb{E}}\|{\mathbf{X}}{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{1}\leq
n\sqrt{\frac{2}{\pi}}\sqrt{C_{{\tiny
u}}}\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}$. More
details of calculations can be found in Corollary 2. Furthermore, with $n\geq
C_{1}\tilde{s}\log(2d/\tilde{s})$, it has

$\displaystyle
f({\boldsymbol{\beta}}+{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}}))-f({\boldsymbol{\beta}})\leq
n\sqrt{C_{{\tiny
u}}}\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}.$

Then the above equation and Equation (5) lead to

$\displaystyle\|{\mathcal{P}}_{\Omega\cup\Pi\cup\Omega^{*}}({\mathbf{G}})\|_{2}\leq
n\sqrt{C_{{\tiny u}}}.$

Secondly, analyze phase two properties, when $C_{2}b_{0}\sqrt{\frac{C_{{\tiny
u}}}{C_{{\tiny
l}}^{2}}\cdot\frac{\tilde{s}\log(2d/\tilde{s})}{n}}\leq\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}\leq
8\sqrt{C_{{\tiny l}}^{-1}}\gamma$. With Lemma 2 and event
${\boldsymbol{\mathcal{E}}}$, we have

$\displaystyle f({\boldsymbol{\beta}})-f({\boldsymbol{\beta}}^{*})$
$\displaystyle\geq{\mathbb{E}}f({\boldsymbol{\beta}})-{\mathbb{E}}f({\boldsymbol{\beta}}^{*})-C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}$
$\displaystyle\geq\frac{n}{6b_{0}}C_{{\tiny
l}}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}^{2}-C\sqrt{n\tilde{s}C_{{\tiny
u}}\log(2d/\tilde{s})}\|{\boldsymbol{\beta}}-{\boldsymbol{\beta}}^{*}\|_{2}$