Adaptive Huber Regression: Optimality and
Phase Transition

arXiv:1706.06991v1 [] 21 Jun 2017

Qiang Sun∗, Wen-Xin Zhou† and Jianqing Fan‡

Abstract
Big data are often contaminated by outliers and heavy-tailed errors, which
makes many conventional methods inadequate. To address this challenge, we
propose the adaptive Huber regression for robust estimation and inference. The
key observation is that the robustification parameter should adapt to the sample
size, dimension and moments for optimal tradeoff between bias and robustness.
Our framework is able to handle heavy-tailed data with bounded (1+δ)-th moment for any δ > 0. We establish a sharp phase transition for robust estimation
of regression parameters in both low and high dimensions: when δ ≥ 1, the
estimator admits a sub-Gaussian-type deviation bound without sub-Gaussian
assumptions on the data, while only a slower rate is available in the regime
0 < δ < 1 and the transition is smooth and optimal. Moreover, a nonasymptotic
Bahadur representation for finite-sample inference is derived when the variance
is finite. Numerical experiments lend further support to our obtained theory.

Keywords: Bias and robustness tradeoff, adaptive Huber regression, finite-sample inference, heavy-tailed data, minimax optimality, phase transition.

1

Introduction

Modern data acquisitions have facilitated the collection of massive and high dimensional
data with complex structures. Along with holding great promises for discovering subtle
population patterns that is less achievable with small-scale data, big data have introduced
a series of new challenges to data analysis both computationally and statistically [Loh and
Wainwright (2015); Fan et al. (2017)]. During the last two decades, extensive progress has
been made towards extracting useful information from massive data with high dimensional
features and sub-Gaussian tails [Tibshirani (1996); Fan and Li (2001); Efron et al. (2004);
Bickel, Ritov and Tsybakov (2009)]. We refer to the monographs, Bühlmann and van de
∗

Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ
08544, USA. E-mail: qsun.ustc@gmail.com.
†
Department of Mathematics, University of California, San Diego, La Jolla, CA 92093, USA. E-mail:
stevewenxin@hotmail.com. Supported by NIH Grant R01-GM072611-12.
‡
Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ
08544, USA. E-mail: jqfan@princeton.edu. Supported by NSF Grants DMS-1206464 and DMS-1406266 and
NIH Grant R01-GM072611-12.

1

Geer (2011) and Hastie, Tibshirani and Wainwright (2015), for a systematic coverage of
contemporary statistical methods and theory for high dimensional data. The sub-Gaussian
condition, although is quite convenient for theoretical analysis, is not really realistic in
many practical applications since modern data are often collected with low quality. In fact,
with an enormous amount of automatically measured features, by chance alone, some of
them may exhibit tails that are much heavier than those of sub-Gaussian or sub-exponential
distributions. A recent study on functional magnetic resonance imaging (fMRI) [Eklund,
Nichols and Knutsson (2016)] shows that the principal cause of invalid fMRI inferences is
that data do not follow the assumed Gaussian shape, which speaks to the need of validating
the statistical methods being used in the field of neuroimaging. In a microarray data
example considered in Wang, Peng and Li (2015), it is observed that some gene expression
levels have heavy tails as their kurtosises are much larger than 3, the kurtosis of a normal
distribution. Outliers frequently arise in microarray experiments due to the array chip
artifacts such as uneven spray of reagents within arrays. In addition, it is well known that
financial returns and microeconomic variables often exhibit heavy tails. Therefore, it is
imperative to develop new statistical inference procedures that are robust against heavytailed errors and other potential forms of contamination.
Standard statistical procedures, such as the least squares method, often behave poorly in
the presence of heavy-tailed data [Catoni (2012)]. The term “heavy-tailed random variable
X” used in this paper is to describe the variable X for which the measure of its tail,
P(|X| > t), decays to zero polynomially in 1/t as t → ∞. An important step towards
estimation under heavy-tailedness has been made by Catoni (2012), whose focus is on
estimating the mean of a univariate probability distribution. Let X be a real-valued, nondegenerate random variable with mean µ = E(X) and variance σ 2 = var(X), and assume
that X1 , . . . , Xn are independent and identically distributed (i.i.d.) copies of X. For any
prespecified exception probability t > 0, Catoni constructs a robust mean estimator µ
bC (t)
that deviates from the true mean µ logarithmically in 1/t, that is,
"

 #
log(1/t) 1/2
P µ
bC (t) − µ ≤ Cσ
≥ 1 − 2t,
(1.1)
n

while the empirical mean deviates from the true mean only polynomially in 1/t. Further,
Devroye et al. (2016) developed adaptive sub-Gaussian estimators that are independent of
the prespecified exception probability. Beyond the problem of mean estimation, Brownlees,
Joly and Lugosi (2015) extended Catoni’s idea to study empirical risk minimization problems when the losses are not bounded. Generalizations of the univariate results to those for
matrices, such as the covariance matrices, can be found in Giulini (2016), Catoni (2016),
Minsker (2016) and Fan, Li and Wang (2017). Fan, Li and Wang (2017) modified Huber’s
procedure [Huber (1973)] to obtain a robust estimator, which is concentrated around the
true mean with exponentially high probability in the sense of (1.1), and also proposed a
robust procedure for high dimensional linear regression under asymmetric and heavy-tailed
errors with bounded second moments.
In this paper, we revisit and investigate further the robust regression problem, which
was initiated by Peter Huber in his seminal work [Huber (1973)]. From a nonasymptotic
perspective, we study the robustness of the Huber regression estimator against heavy-tailed
2

errors, which is different from the contamination model setup in Huber (1973) and is of particular interest in high dimensional data analysis [Delaigle, Hall and Jin (2011)]. Asymptotic
properties of the Huber estimator with a fixed robustification parameter (to be defined in
Section 2) have been well studied in the literature. See, for example, Huber (1973), Yohai
and Maronna (1979), Portnoy (1985), Mammen (1989) and He and Shao (1996, 2000). Fan,
Li and Wang (2017) are among the first to study robust regression with only bounded
second moments and possibly asymmetric error distributions in high dimensions. Yet few
nonasymptotic results on robust regression exist, including the concentration of measure
phenomenon, minimax optimality and finite-sample inference. To close this gap, we propose and study the Huber regression with an adaptive robustification parameter, which is
referred to as adaptive Huber regression, for robust estimation and inference. In particular,
we do not assume the symmetry and homoscedasticity for error distributions, so that the
problem of current interest is intrinsically different from median and quantile regressions.
We made four major contributions towards robust modeling in this paper. First and
foremost, we establish nonasymptotic deviation bounds for robust regression estimators
when the error variables have only finite (1 + δ)-th moments. By providing the matching
lower bound, we observe a sharp phase transition phenomenon, which is in line with that
discovered by Devroye et al. (2016) in the univariate mean estimation problem. Second,
we derive a nonasymptotic Bahadur representation for the adaptive Huber estimator when
δ ≥ 1. This provides a theoretical foundation for robust finite-sample inference. Third,
a similar phase transition for sparse adaptive Huber regression is established in high dimensions. By defining the effective dimension and effective sample size, we present the
nonasymptotic results from the two different regimes in a unified form. Last, by exploiting
the localized analysis developed in Fan et al. (2017), we remove the artificial bounded parameter constraint imposed in previous works. See, for example, Spokoiny (2012), Loh and
Wainwright (2015) and Fan, Li and Wang (2017).
The rest of the paper proceeds as follows. In Section 2, we define the Huber loss
and robustification parameter, followed by the proposal of adaptive Huber regression in
both low and high dimensions. We sharply characterize the nonasymptotic performance
of the proposed estimator in Section 3. We describe the algorithm and implementation in
Section 4. Section 5 is devoted to numerical studies. We conclude the paper with further
discussions in Section 6. Proofs of the main theorems are given in the Appendix. Due to the
page limitation, proofs of the lemmas and propositions, along with some additional results,
are collected in the supplemental material [Sun, Zhou and Fan (2017)].
Notation: We fix some notations that will be used throughout this paper. For any vector
P
u = (u1 , . . . , ud )T ∈ Rd and q ≥ 1, kukq = ( dj=1 |uj |q )1/q is the `q norm. For any vectors
P
u, v ∈ Rd , we write hu, vi = uT v. Moreover, we let kuk0 = dj=1 1(uj 6= 0) denote the
number of nonzero entries of u, and set kuk∞ = max1≤j≤d |uj |. For two sequences of
real numbers {an }n≥1 and {bn }n≥1 , an . bn denotes an ≤ Cbn for some constant C > 0
independent of n, an & bn if bn . an , and an  bn signifies the case that an . bn and bn . an .
For two scalars, we use a ∧ b = min{a, b} to denote the minimum of a and b. If A is an
m×n matrix, we use kAk to denote its spectral norm, defined by kAk = maxu∈Sn−1 kAuk2 ,
where Sn−1 = {u ∈ Rn : kuk2 = 1} is the unit sphere in Rn . For an n × n matrix A, we use
3

λmax (A) and λmin (A) to denote the maximum and minimum eigenvalues of A, respectively.
For two n × n matrices A and B, we write A  B if B − A is positive semi-definite. For a
function f : Rd → R, we use ∇f to denote its gradient vector in Rd as long as it exists.

2

Methodology

Consider a heteroscedastic linear regression model y = xT β ∗ + ε, from which we observe n
samples, {(yi , xi ) : i = 1, . . . , n}, satisfying

∗
1+δ
yi = xT
< ∞.
i β + εi , with E(εi |xi ) = 0 and vi,δ = E |εi |

(2.1)

Assuming only the finiteness of second moments, the standard ordinary least squares (OLS)
estimator, denoted by βbols , only admits a suboptimal polynomial-type deviation bound, and
thus does not concentrate around β ∗ enough for some applications (e.g., large-scale simultaneous inference). The key observation that underpins this suboptimality of the OLS
estimator is the sensitivity of quadratic loss to outliers [Huber (1973); Catoni (2012)]. To
overcome this drawback, we propose to use the Huber loss with an adaptive robustification
parameter to robustify the estimation procedure against heavy-tailed errors and other potential forms of contamination. We begin with the definitions of the Huber loss and the
related robustification parameter.
Definition 2.1 (Huber Loss and Robustification Parameter). The Huber loss `τ (·) [Huber
(1964)] is defined as
(
1 2
if |x| ≤ τ,
2x ,
`τ (x) =
1 2
τ |x| − 2 τ , if |x| > τ,
where τ > 0 is referred to as the robustification parameter that trades bias for robustness
in mean estimation problems.
The Huber loss, `τ (x), is quadratic for small values of x, and becomes linear when x
gets larger. The parameter τ controls the blending of quadratic and `1 losses, which can be
regarded as two extremes of the Huber loss with τ = ∞ and τ → 0, respectively. Comparing
with the least squares, outliers are down weighted or removed in the Huber loss. What is
different from the classical setting is that we allow the parameter τ to adapt to the sample
size, dimension and moments for a better tradeoff between bias and robustness. As τ → ∞
is needed to reduce the bias when the error distribution is asymmetric, this loss is also called
the RA-quadratic (robust approximation to quadratic) loss in Fan, Li and Wang (2017).
P
d
Define the loss function Lτ (β) = n−1 ni=1 `τ (yi − xT
i β) for β ∈ R . The Huber estimator is defined through the following convex optimization problem:
βbτ = argmin Lτ (β).

(2.2)

β∈Rd

P
In low dimensions, under the condition that vδ = n−1 ni=1 E(|εi |1+δ ) < ∞ for some
1/(1+δ) 1/2
δ > 0, we will prove that βbτ with τ  min{vδ
, v1 } nmax{1/(1+δ),1/2} (the first factor
is kept in order to show its explicit dependence on the moment) achieves the minimax
√
√
optimal rate of convergence d τ −(δ∧1)  d n− min{δ/(1+δ),1/2} , where d is the dimension
4

of β ∗ . The phase transition at δ = 1 can be easily observed. Proposition 1 in Section 3.1.1
indicates that the bias introduced by the Huber loss scales at the rate τ −δ . This suggests
that when 0 < δ ≤ 1, the approximation bias matches the optimal rate of convergence,
which is of order n−δ/(1+δ) , when holding d fixed. Additionally, robustification leads to a
sub-Gaussian-type deviation inequality in the sense of (1.1). When higher moments exist,
that is δ > 1, the approximation bias decays at a faster rate and is negligible compared to
the estimation error, which admits convergence rate in the order of n−1/2 .
n

, d

n

(n/ log d)

, d

n

impossible region

1/2

=

=

1+

0

1
2

1

Figure 1: Phase transition in terms of `2 -rate of convergence for the adaptive Huber estimator.

−δ/(1+δ)
−1/2
With fixed effective dimension, kβbτ − β ∗ k2  neff
, when 0 < δ < 1; kβbτ − β ∗ k2  neff , when
δ ≥ 1. Here neff is the effective sample size: neff = n in low dimensions while neff = n/ log d in high
dimensions.

In the high dimensional regime, we propose the following sparse adaptive Huber regression with a slightly different choice of the robustification parameter:

βbτ,λ = argmin Lτ (β)+λkβk1 , where τ  νδ
β∈Rd

λ  νδ



log d
n



n
log d

min{δ/(1+δ),1/2}

,

max{1/(1+δ),1/2}

and
(2.3)

1/(1+δ) 1/2
where νδ = min{vδ
, v1 }. When it is clear from the context, the dependence of βbτ,λ on
λ is sometimes suppressed. Let s be the size of the true support S ≡ supp(β ∗ ). Then, the
proposed sparse Huber estimator achieves the minimax optimal rate of convergence, which
min{δ/(1+δ),1/2}
√
, for estimating β ∗ in `2 -error with high probability. For
is of order s logn d
δ/(1+δ)
√
√
0 < δ ≤ 1, the obtained estimator has the convergence rate of order s λ  s logn d
,
−δ
and the approximation bias scales at the rate τ , matching the order of shrinkage bias λ.
To unify the rates of convergence in the above two different regimes, we define the
effective dimension, deff , to be d in low dimensions and s in high dimensions. In other words,
deff denotes the number of nonzero (unknown) parameters of the problem. The effective
sample size, neff , is defined as neff = n and neff = n/ log d in low and high dimensions,

5

respectively. We will establish a sharp phase transition: when δ ≥ 1, the proposed estimator
enjoys a sub-Gaussian rate of convergence, while it only achieves a slower rate when 0 <
δ < 1. This is summarized in Figure 1. Specifically, we will show that, for any δ ∈ (0, ∞),
1/(1+δ) 1/2
max{1/(1+δ),1/2}
the proposed estimators with τ  min{vδ
, v1 } neff
achieve the following
optimal rate of convergence:

3

βbτ − β ∗

1/2

2

− min{δ/(1+δ),1/2}

 deff neff

,

with large probability.

(2.4)

Nonasymptotic Theory

In this section, we study the concentration of measure phenomenon and minimax optimality
for the proposed estimators. In addition, a nonasymptotic Bahadur representation is derived
for finite-sample inference when second moments of the error variables exists.

3.1

Adaptive Huber Regression

We start with adaptive Huber regression in the low dimensional setting. First, we develop
a bias approximation bound for the Huber regression coefficient to the underlying true
parameters. We then establish the phase transition by providing both the upper and lower
bounds. Based on these results, we further develop a nonasymptotic Bahadur representation
for finite-sample inference. The results are carried out under both fixed and random designs.
3.1.1

Fix Design

We investigate the robust estimator βbτ given in (2.2) in the low dimensional regime under
fix designs. We make the following regularity condition.

Condition 3.1. The covariate vectors {xi }ni=1 are standardized such that max1≤i≤n ke
xi k∞ ≤
P
−1/2
ei ≡ Sn xi and Sn ≡ n−1 ni=1 xi xT
is
nonsingular.
1, where x
i
Note that, for any τ > 0, βbτ given in (2.2) is a sample version of the following M estimator
n

βτ∗

1X
≡ argmin E{Lτ (β)} = argmin
E{`τ (yi − xT
i β)},
n
β∈Rd
β∈Rd

(3.1)

i=1

where the expectation is taken over the error variables ε1 , . . . , εn . We call βτ∗ the Huber
regression coefficient, which is in general different from the vector of true parameters β ∗ .
The presence of the approximation error kβτ∗ −β ∗ k2 is a direct consequence of robustification
and asymmetric error distributions. Heuristically, choosing a sufficiently large τ reduces this
approximation error at the cost of losing robustness (the extreme case of τ = ∞ corresponds
to the least squares estimator). Our first result reveals that how the size of τ affects the
approximation error kβτ∗ − β ∗ k2 . Its proof is given in the supplemental material [Sun, Zhou
P
and Fan (2017)]. Recall that vδ = n−1 ni=1 vi,δ , where vi,δ = E(|εi |1+δ ) is defined in (2.1).

6

Proposition 1. Assume that vδ < ∞ for some δ > 0 and the matrix Sn = n−1
is nonsingular. Then βτ∗ must satisfy
(
2vδ τ −δ , if 0 < δ < 1, τ ≥ (4vδ m2n )1/(1+δ) ,
1/2
∗
∗
Sn (βτ − β ) 2 ≤
2vδ τ −δ , if δ ≥ 1, τ ≥ (2v1 )1/2 mn ,
−1/2

where mn ≡ max1≤i≤n kSn

Pn

T
i=1 xi xi

(3.2)

xi k2 .

The following result gives a nonasymptotic upper bound on the `2 -error with an exponential1/(1+δ) 1/2
type exception probability. Recall that νδ = min{vδ
, v1 } for any δ > 0.
Theorem 3.1 (Upper Bound). Assume that Condition 3.1 holds and that vδ is finite for
some 0 < δ < ∞. Then, for any t > 0 and τ0 ≥ νδ , with probability at least 1 − (2d + 1)e−t ,
the estimator βbτ with τ = τ0 (n/t)max{1/(1+δ),1/2} must satisfy
b
S1/2
n (βτ

∗

−β )

2

≤ 4τ0

√

 min{δ/(1+δ),1/2}
t
d
,
n

(3.3)

as long as n ≥ 32d2 t. Further assume that c` ≡ λmin (Sn ) > 0, then we must have
βbτ − β ∗

2

 
√ t min{δ/(1+δ),1/2}
.
≤ 4τ0 c−1
d
`
n

Remark 1. It is worth mentioning that the proposed robust estimator depends on both
1/(1+δ)
the prespecified confidence level and the unknown parameter vδ
(standard deviation
when δ = 1). The choice of t is user-specific and flexible, which may depend on the
particular problem of interest. See, for example, Zhou et al. (2016) for applications of the
adaptive Huber regression to large-scale multiple testing under dependence. Adaptation
to the unknown (1 + δ)-th moment is indeed a much more complex issue. In Section 5,
we suggest a simple data-driven choice of τ with nice numerical performance. A general
adaptive construction of τ can be obtained via Lepski’s method [Lepski (1991)], which
1/(1+δ)
replies on “crude” preliminary lower and upper bounds on vδ
. We leave the theoretical
guarantee and numerical implementation of this procedure for future research.
Remark 2. Note that we do not assume E(|εi |1+δ |xi ) to be a constant, and hence the
proposed method accommodates heteroscedastic regression models, where εi may possibly
depend on xi . For example, εi can take the form of σ(xi )ιi , where σ(·) : Rd → (0, ∞) is a
positive function, and ιi ’s are random variables satisfying E(ιi ) = 0 and E(|ιi |1+δ ) < ∞.
Theorem 3.1 indicates that, with only bounded (1+δ)-th moment, the adaptive Huber
estimator achieves n− min{δ/(1+δ),1/2} rate of convergence with exponentially high probability.
A natural question that arises is whether the upper bound in (3.3) is optimal. To address
this question, we provide a matching lower bound up to a logarithmic factor. Let Pδvδ be the
class of all distributions on R whose (1 + δ)-th absolute central moment equals vδ . Let X =
(x1 , . . . , xn )T = (x1 , . . . , xd ) ∈ Rn×d be the design matrix and Un = {u : u ∈ {−1, 1}n }.
Theorem 3.2 (Lower Bound). Let ε1 , . . . , εn be i.i.d. copies of the random variable ε
with distribution P. Assume that there exists u ∈ Un such that kn−1 XT ukmin ≥ α for
7

some α > 0. Suppose that δ ∈ (0, ∞). Then, for any t ∈ [0, n/2] and any estimator
b 1 , . . . , yn , t) possibly depending on t, we have
βb = β(y
"
√  min{δ/(1+δ),1/2} #
α
d t
e−2t
∗
sup
P βb − β 2 ≥ νδ
≥
,
v
cu
n
2
P∈P δ , β ∗ ∈Rd
δ

where cu = λmax (Sn ).

Theorem 3.2 reveals the nonexistence of sub-Gaussian regression estimators when δ <
1. It widens the phenomenon observed in Theorem 3.1 in Devroye et al. (2016) on the
nonexistence of sub-Gaussian mean estimators. In addition to the eigenvalue assumption,
we need to assume that there exists a u ∈ Un ⊂ Rn such that the minimum angle between
n−1 u and xj ’s is non-diminishing. This assumption comes from the intuition that the linear
subspace spanned by xj ’s is at most rank d and thus can not span the whole space Rn .
We note that this assumption naturally holds in the univariate mean estimation problem
where X = (1, . . . , 1)T and we can take u = (1, . . . , 1)T and α = 1. This lower bound
result, together with the upper bound, shows that the adaptive Huber estimator achieves
near-optimal (up to a logarithmic factor) deviations. Moreover, it indicates that the Huber
estimator with an adaptive robustification parameter τ has a sharp phase transition: when
√
δ ≥ 1, βbτ converges to the true parameter at the standard rate d n−1/2 , while only a
√
slower rate of order d n−δ/(1+δ) is available when the second moment does not exist.
In addition to estimation, construction of confidence sets is another important problem
in statistics. Polynomial-type deviation inequalities, such as those for the least square
estimators with heavy-tailed data, typically result in too wide confidence intervals. This is
because the length of the confidence interval does not only scale with n−1/2 but also scales
with the exception probability polynomially; see the discussion below (1.1).
Faster concentration of the adaptive Huber estimator enables us to construct sharper
confidence sets. To achieve this, we prove a nonasymptotic version of the Bahadur representation. In particular, we derive an exponential-type tail probability bound for the remainder
of the Bahadur representation for βbτ when the error variable has a finite second moment. In
addition to Condition 3.1, we also need the following assumption on the covariate vectors.
P
Recall that Sn ≡ n−1 ni=1 xi xT
i .
Condition 3.2. There exist constants κ, M > 0 such that
n

sup
u∈Sd−1
−1/2

e i = Sn
where x

1X T 4
ei ) exp(κ|uT x
ei |2 ) ≤ M,
(u x
n
i=1

xi ’s are as in Condition 3.1.

Condition 3.2 is needed to bound the tail probability of quadratic forms of the Hessian
matrix ∇2 Lτ (β) in a neighborhood of β ∗ by using tools from empirical processes; see the
proof of Theorem 3.3 for details. It is a non-stochastic analog of the sub-Gaussian assumption on the covariate vectors x1 , . . . , xn . To see this, note that if xi ’s are i.i.d. random
P
T
2
ei )4 eκ|u xei | →
vectors drawn from x, the law of large numbers reveals that n−1 ni=1 (uT x
T 2
e)4 eκ|u xe| almost surely as n → ∞, where x
e = Σ−1/2 x with Σ = cov(x) and
E(uT x
T
2
E(e
xT u)4 eκ|u xe| is finite provided that x is sub-Gaussian; see Condition 3.3 in Section 3.1.2
for random design analysis.
8

Theorem 3.3 (Nonasymptotic Bahadur Representation). Assume that Conditions 3.1 and
1/2
3.2 hold, and that v1 is finite. Then, for any t > 0 and τ0 ≥ v1 , the estimator βbτ given in
(2.2) with τ = τ0 (n/t)1/2 satisfies that, with probability greater than 1 − 2(d + 1)e−t ,
√
n
1X 0
1/2 dt
1/2 b
∗
`τ (εi )e
xi ≤ Aτ0 (d + t)
(3.4)
Sn ( β τ − β ) −
n
n
2
i=1

whenever n ≥ max{32d2 t, 2κ−2 (2d + t)}, where ∆n ≡ n−1
constant depending only on M and τ0−2 k∆n k.

Pn

ei x
eT
i
i=1 vi,1 x

and A > 0 is a

Theorem 3.3 provides a finite-sample approximation of βbτ by a sum of independent
random vectors. The remainder of such an approximation exhibits sub-exponential tails as
long as the error variables have finite second moments. Unlike the least squares estimator,
the adaptive Huber estimator does not admit an explicit closed-form representation, which
causes the main difficulty for analyzing its nonasymptotic properties. Theorem 3.3 reveals
that, up to a higher-order remainder, the property of βbτ mainly depends on its linear
approximation term, which is much easier to deal with.
Note that the OLS estimator βbols admits a closed-form expression that βbols − β ∗ =
Pn
−1
a heavy-tailed distribution, the worst case deviations
S−1
n (n
i=1 εi xi ). When εi follows
P
of the empirical mean, say n−1 ni=1 εi xi , can be large and suboptimal [Catoni (2012)].
Therefore, the deviations of the OLS estimator to β ∗ may increase remarkably when the
error distribution is far from Gaussian, while as indicated by Theorem 3.1, the new estimator
βbτ defined in (2.2) with a properly chosen τ will not. In addition, Theorem 3.3 shows
that although βbτ does not have a closed-form expression, βbτ − β ∗ can be approximated
Pn 0
0
−1
by S−1
n {n
i=1 `τ (εi )xi } with a higher-order remainder. Note that `τ (ε) = εi 1(|ε| ≤
τ ) + τ 1(ε > τ ) − τ 1(ε < −τ ) is a truncated version of ε. Hence, the corresponding empirical
P
average n−1 ni=1 `0τ (εi )xi is less sensitive to outliers.
Regarding the truncated random variable `0τ (εi ), the following result shows that the
differences between the first two moments of `0τ (εi ) and εi depend on τ and the moments
of εi . The higher moment εi has, the faster these differences decay as a function of τ . We
summarize this observation in the following proposition. We drop i for ease of presentation.
Proposition 2. Assume that E(ε) = 0, σ 2 = E(ε2 ) > 0 and E(|ε|2+ι ) < ∞ from some
ι ≥ 0. Then, we have

Moreover, if ι > 0,



|E`0τ (ε)| ≤ min τ −1 σ 2 , τ −1−ι E |ε|2+ι .


σ 2 − 2ι−1 τ −ι E |ε|2+ι ≤ E{`0τ (ε)}2 ≤ σ 2 .

Proposition 2, along with Theorem 3.3, reveals that the adaptive Huber estimator can
achieve nonasymptotic robustness against heavy-tailed errors, while enjoying high efficiency
when τ diverges to ∞. The proof of this result is given in the supplementary material [Sun,
Zhou and Fan (2017)].

9

3.1.2

Random Design

In this section, we derive counterparts of Theorems 3.2 and 3.3 under random designs.
Condition 3.3. The covariate vector x ∈ Rd is sub-Gaussian, i.e.,
Σ−1/2 x

ψ2

≡ sup

uT Σ−1/2 x

u∈Sd−1

ψ2

≤ A0 < ∞

for some positive constant A0 , where Σ = E(xxT ) is invertible, and k · kψ2 is the subGaussian norm defined by kXkψ2 = sup`≥1 `−1/2 (E|X|` )1/` for any random variable X
[Vershynin (2012)].
Throughout this section, for simplicity, we assume the independent error variables
εi ’s in model (2.1) are independent of xi ’s and have mean zero. The conditional heteroscedastic model can be allowed with slightly modifications as before. With this setup,
P
1/(1+δ) 1/2
vδ = n−1 ni=1 vi,δ with vi,δ = E(|εi |1+δ ). Recall that νδ ≡ min{vδ
, v1 }.
Theorem 3.4. Suppose Condition 3.3 holds.

(I) Assume that vδ < ∞ for some δ > 0. Then, for any t > 0 and τ0 ≥ νδ , the estimator
βbτ with τ = τ0 (n/t)max{δ/(1+δ),1/2} satisfies that, with probability greater than 1 −
(2d + 3)e−t ,
1/2

Σ

(βbτ − β )
∗

1/2

2

≤ C1 τ0 (log n + t)

√

 min{δ/(1+δ),1/2}
t
d
n

(3.5)

as long as n ≥ C2 (log n + t)dt, where C1 , C2 > 0 are constants depending only on A0 .
1/2
(II) Assume that v1 < ∞. Then, for any t > 0 and τ0 ≥ v1 , the estimator βbτ with
p
√
τ = τ0 n/(d + t) satisfies that as long as n ≥ C3 d (d + t),
r


d+t
1/2
1/2 b
∗
P Σ (βτ − β ) 2 ≥ C4 v1
≤ 9e−t
(3.6)
n

and

P



n

Σ1/2 (βbτ − β ∗ ) −

1X 0
`τ (εi )e
xi
n
i=1

1/2 d

2

≥ C5 v1

+t
n



≤ 10e−t ,

(3.7)

where C3 –C5 are positive constants depending only on (A0 , τ0 ).
Remark 3. Observe that the deviation bound (3.5) with δ = 1 and t  log n reveals the
p
convergence rate (log n)2 d/n under the scaling that n & (log n)2 d; while the bound (3.6)
p
gives a faster rate of order (log n + d)/n under a stronger scaling condition n & d3/2 .
The main difference comes from the approach we use to bound the quadratic forms of a
d-dimensional random vector, say ξ, as described below.
(1) Under weaker moment conditions that E(|ε|1+δ ) < ∞ for some 0 < δ ≤ 1, we first
√
use the elementary inequality kξk2 ≤ d kξk∞ , and then bound the tail probability
of kξk∞ using the union bound and a bound on the moment generating function.
10

(2) Under the finite second moment condition, we directly bound the tail probability of
kξk2 using techniques from empirical process theory [Spokoiny (2013)].
The first part of Theorem 3.4 gives the rate of convergence, whereas the second part
gives a Bahadur representation for the estimator under random designs. In the second part
of Theorem 3.4, by taking t = A log n for some constant A > 0, we see that as long as
p
n & d3/2 , the robust estimator βbτ with τ  n/(d + A log n) satisfies


n
1X 0
d + log n
βbτ = β ∗ +
`τ (εi )Σ−1 xi + O
n
n
i=1

with probability at least 1 − O(n−A ). From an asymptotic point of view, this implies that
if the dimension d grows with n in such a way that
d3/2 = o(n) as n → ∞,

(3.8)

then for any deterministic vector a ∈ Rd , the distribution of aT (βbτ −β ∗ ) coincides with that
P
of n−1 ni=1 `0τ (εi )aT Σ−1 xi asymptotically as long as kΣ−1/2 ak22 converges. If ε1 , . . . , εn
are independent copies of ε with variance σ 2 and E(|ε|2+δ ) < ∞ for some δ > 0, taking
p
P
τ  n/(d + log n) in Proposition 2 implies that n−1/2 ni=1 `0τ (εi )aT Σ−1 xi follows a normal distribution with mean zero and variance σ 2 kΣ−1/2 ak22 asymptotically. It is worth
mentioning that the scaling condition (3.8) improves that in Portnoy (1985) by getting rid
of an additional log n term.

3.2

Adaptive Huber Regression in High Dimensions

In this section, we sharply characterize the statistical performance of the sparse adaptive
Huber estimator in high dimensions when the dimension d is allowed to grow with sample
size n exponentially fast. The analysis is carried out under fix designs, and can be extended
to random designs based on a conditional argument, in a similar spirit to the previous
section. We start with a modified version of the localized restricted eigenvalue introduced
by Fan et al. (2017), which is needed for high dimensional nonasymptotic analysis under
general loss functions. Let Hτ (β) ≡ ∇2 Lτ (β) denote the Hessian matrix. Recall that
S ≡ supp(β ∗ ) ⊆ {1, . . . , d}, with |S| = s, is the true support set.
Definition 3.1 (Localized Restricted Eigenvalue, LRE). The localized restricted eigenvalue
of Hτ (β) is defined as
n
o
κ+ (m, γ, r) = sup uT Hτ (β)u : (u, β) ∈ C(m, γ, r) ,
n
o
κ− (m, γ, r) = inf uT Hτ (β)u : (u, β) ∈ C(m, γ, r) ,

where C(m, γ, r) ≡ {(u, β) ∈ Sd−1 ×Rd : ∀J ⊆ {1, . . . , d} satisfying S ⊆ J, |J| ≤ m, kuJ c k1 ≤
γkuJ k1 , kβ − β ∗ k1 ≤ r} is a local `1 cone.
The LRE is defined in a local neighborhood of β ∗ in the sense of `1 -distance. This
facilitates our proof, while Fan et al. (2017) use the `2 -distance to define the local region.
11

Condition 3.4. For Hτ (β), we say that the localized restricted eigenvalue condition holds
with the triplet (k, γ, r), denoted as LRE(k, γ, r), if 0 < κl ≤ κ− (k, γ, r) ≤ κ+ (k, γ, r) ≤
κu < ∞.
The condition above is referred to as the LRE condition [Fan et al. (2017)]. It is a unified
condition for studying generalized loss functions, whose Hessians may possibly depend on β.
For Huber loss, Condition 3.4 also depends on the error variables. The following definition
concerns the restricted eigenvalues of Sn instead of Hτ (β).
Definition 3.2 (Restricted Eigenvalue, RE). The restricted maximum and minimum eigenvalues of Sn are defined respectively as

ρ+ (m, γ) = sup uT Sn u : u ∈ C(m, γ) ,
u

ρ− (m, γ) = inf uT Sn u : u ∈ C(m, γ) ,
u

where C(m, γ) ≡ {u ∈ Sd−1 : ∀J ⊆ {1, . . . , d} satisfying S ⊆ J, |J| ≤ m, kuJ c k1 ≤ γkuJ k1 }.

Condition 3.5. For Sn , we say that the restricted eigenvalue condition holds with the
doublet (k, γ), denoted as RE(k, γ), if 0 < κl ≤ ρ− (k, γ) ≤ ρ+ (k, γ) ≤ κu < ∞.
To make Condition 3.4 practically useful, in what follows, we show that Condition 3.5
P
implies Condition 3.4 with high probability. Recall that vδ = n−1 ni=1 vi,δ .

Lemma 3.1. Assume that the covariates are standardized such that max1≤i≤n kxi k∞ ≤ 1.
For any prespecified integer k ≥ 1 and r, t > 0, suppose that τ & (vδ k)1/(1+δ) ∨r and n & k 2 t.
Then, Condition 3.5 implies Condition 3.4 with high probability, i.e., if 0 < κl ≤ ρ− (k, γ) ≤
ρ+ (k, γ) ≤ κu < ∞, we must have 0 < κl /2 ≤ κ− (k, γ, r) ≤ κ+ (k, γ, r) ≤ κu < ∞, with
probability at least 1 − e−t .
With the above preparations in place, we are now ready to present the main results on
1/(1+δ) 1/2
the adaptive Huber estimator in high dimensions. Recall that νδ ≡ min{vδ
, v1 }.
Theorem 3.5 (Upper Bound in High Dimensions). Assume Condition 3.5 holds with (k, γ) =
(2s, 3) and the covariate vectors x1 , . . . , xn are standardized such that max1≤i≤n kxi k∞ ≤ 1.
For any t > 0, let τ = τ0 (n/t)1/(1+δ) for some τ0 ≥ νδ and λ ≥ 4τ n−1 t. Then, with
probability at least 1 − (2s + 1)e−t , the penalized adaptive Huber estimator βbτ,λ defined in
(2.3) satisfies
√
βbτ,λ − β ∗ 2 ≤ 3κ−1
(3.9)
l λ s,
as long as n & s2 t. In particular, for any A > 0, by taking t = (A + 1) log d such that
τ  νδ



n
log d

max{1/(1+δ),1/2}

and λ  νδ



log d
n

min{δ/(1+δ),1/2}

then, with probability at least 1 − d−A , we have


(1 + A) log d min{δ/(1+δ),1/2}
−1 √
∗
b
βτ,λ − β 2 . κ` νδ s
.
n
12

,

(3.10)

The result above suggests that the sparse adaptive Huber estimator with an adaptive
√
robustification parameter converges at the rate s {(log d)/n}max{δ/(1+δ),1/2} with overwhelming probability. When only the second moment of the error variable exists, the
proposed estimator performs as well as the penalized OLS regression as if sub-Gaussian
errors were assumed. We advocate the use of Huber’s robust regression since sub-Gaussian
condition often fails in practice [Wang, Peng and Li (2015); Eklund, Nichols and Knutsson
(2016)]. Finally, we prove a matching lower bound for estimating β ∗ in the s-sparse set
B0 (s) ≡ {β ∈ Rd : kβk0 = s}. Recall the definition of Un in Theorem 3.2.
Theorem 3.6 (Lower Bound in High Dimensions). Let ε1 , . . . , εn be i.i.d. copies of the
random variable ε with distribution P. Suppose that Condition 3.5 holds with k = 2s
and γ = 0. Further assume that there exist a set A with |A| = s and u ∈ Un such
that kXT
A u/nkmin ≥ α for some α > 0. Then, for any A > 0 and s-sparse estimator
b 1 , . . . , yn , A) possibly depending on A, we have
βb = β(y
"
min{δ/(1+δ),1/2} #
√ 
α
s
A
log
d
≥ 2−1 d−A
sup
P βb − β ∗ 2 ≥ νδ
vδ
κ
2n
∗
u
P∈P , β ∈B0 (s)
δ

as long as n ≥ 2(A log d + log 2).

Theorem 3.6, together with the upper bound established in Theorem 3.5, indicates that
the sparse adaptive Huber estimator achieves the optimal rate of convergence in `2 -error.
The proof, which is given in Section A.6, involves constructing a sub-class of binomial
distributions for the error variables. Combining with the results in Section 3.1.1, we have
proved the unified version in (2.4) and thus the phase transition in Figure 1.

4

Algorithm and Implementation

This section is devoted to computational algorithm and numerical implementation. We
focus on the sparse adaptive Huber regression in (2.3), as (2.2) can be solved by setting λ to
be zero. The optimization problem in (2.3) is convex, and thus standard convex optimization
tools such as cutting-plane or interior point method can be exploited. However, they are
not scalable to large-scale problems.
In what follows, we suggest a fast and easily implementable method using the local adaptive majorize-minimization (LAMM) principle [Fan et al. (2017)]. We say that a function
g(β|β (k) ) majorizes f (β) at the point β (k) if
g(β|β (k) ) ≥ f (β)

and g(β (k) |β (k) ) = f (β (k) ).

To minimize a general function f (β), a majorize-minimization (MM) algorithm initializes
at β (0) , and then iteratively computes β (k+1) = argminβ∈Rd g(β|β (k) ). The objective value
of such an algorithm decreases in each step, since
major.

min.

init.

f (β (k+1) ) ≤ g(β (k+1) | β (k) ) ≤ g(β (k) | β (k) ) = f (β (k) ).

13

(4.1)

Algorithm 1 The LAMM algorithm for sparse adaptive Huber regression.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:

(0)
Algorithm: {β (k) , φk }∞
k=1 ← LAMM(λ, β , φ0 ,  )
Input: λ, β (0) , φ0 , 
Initialize: φ(`,k) ← max{φ0 , γu−1 φ(`,k−1) }
for k = 0, 1, . . . until kβ (k+1) − β (k) k2 ≤  do
Repeat
β (k+1) ← Tλ,φk (β (k) )
If gk (β (k+1) |β (k) ) < Lτ (β (k+1) ) then φk ← γu φk
Until gk (β (k+1) |β (k) ) ≥ Lτ (β (k+1) )
Return {β (k+1) , φk }
end for
Output: βb = β (k+1)

As pointed out by Fan et al. (2017), the majorization requirement only needs to hold locally
at β (k+1) when starting from β (k) . We therefore locally majorize Lτ (β) in (2.3) at β (k) by
an isotropic quadratic function
gk (β|β (k) ) ≡ Lτ (β (k) ) + ∇Lτ (β (k) ), β − β (k) +

φk
β − β (k)
2

2
,
2

where φk is a quadratic parameter such that gk (β (k+1) |β (k) ) ≥ Lτ (β (k+1) ). The isotropic
form also allows a simple analytic solution to the subsequent majorized optimization problem:


φk
(k) 2
(k)
(k)
(4.2)
β−β
+λ β 1 .
argmin ∇Lτ (β ), β − β
+
2
2
β∈Rd
It can be shown that (4.2) is minimized at


−1
(k)
β (k+1) = Tλ,φk (β (k) ) ≡ S β (k) − φ−1
∇L
(β
),
φ
λ
,
τ
k
k

where S(x, λ) is the soft-thresholding operator, defined by S(x, λ) ≡ (sign(xj ) · max{|xj | −
λ, 0}). The simplicity of this updating rule is due to the fact that (4.2) is an unconstrained
optimization problem.
To find the smallest φk such that gk (β (k+1) |β (k) ) ≥ Lτ (β (k+1) ), the basic idea of LAMM
is to start from a very small isotropic parameter φk = φ0k and then successfully inflate φk
by a factor of γu > 1 (say, 2). If the solution satisfies gk (β (k+1) |β (k) ) ≥ Lτ (β (k+1) ), we stop
this part of the algorithm and obtain β (k+1) , which makes the target value non-increasing.
We then continue with the iteration to produce next solution until the solution sequence
(k+1) − β (k) k ≤  for a sufficiently
{β (k) }∞
2
k=1 converges. A simple stopping criterion is kβ
small , say 10−5 . We refer to Fan et al. (2017) for a detailed complexity analysis of the
LAMM algorithm.

14

●

2.5

df = 1.5

2.5

df = 1.5
d = 100
d = 500
d = 5000

●

1.5

1.5

2

2

2.0

●

2.0

●

d = 100
d = 500
d = 5000
Theoretical

●

●

●

●

1.0

1.0

●

●

●

●
●

●

0.5

●

0.5

●

0

200

400

600

800

1000

1200

0

n

50

100

150

n/ log d

Figure 2: The scaling behavior of the `2 -error kβbτ − βk2 in high dimensions. Left panels
shows the `2 -error versus sample size, while the right panel shows the `2 -error versus the
effective sample size neff = n/ log d.

5

Numerical Examples

In this section, numerical experiments are carried out to support the theoretical findings.
We generate continuous responses according to the following linear model
∗
yi = x T
i β + εi ,

i = 1, . . . , n,

(5.1)

where εi ∼ tdf , a student’s t-distribution with df being the degree of freedom, and β ∗ is a
d-dimensional vector of coefficients such that
β ∗ = (5, −2, 0, 0, 3, 0, . . . , 0)T .
| {z }
d−5

In model (5.1), xi ’s are generated from the multivariate normal distribution N (0, Id ), which
are independent of εi ’s. We note that tdf has only bounded (1 + δ)-th moments, for any
δ < df − 1. Thus we take δ = df − 1 − 0.05 through out this section.
In our first example, we examine the scaling behavior of the `2 -error kβbτ − βk2 . In order
to predict the scaling behavior, data are generated according to model (5.1), with degree of
the freedom of the t-distributed error taking to be df = 1.5. As suggested by Theorem 3.5,
we choose the robustification parameter to be

max{1/(1+δ),1/2}
n
τ = 0.5 ×
.
log d
It turns out the performance of the adaptive Huber estimator is less sensitive to the choice
of τ and thus 0.5 is chosen mindlessly to demonstrate the phase transition phenomenon.
15

●

●

●

●

●

(n, d) = (500, 1000)

●

●

●

●

●

●

1.2

2.0

(n, d) = (500, 5)

●

●

●

●

●

●

●

●

●

●

●

●

1.0

●

●
●

2)

●

2)

●

●

0.8

●

1.5

●

●

0.6

●

log(

●

0.2

1.0

log(

●

0.4

●

●

0.0

●

●

●

0.5

1.0

1.5

2.0

0.5

1.0

1.5

2.0

Figure 3: The minimax optimality of Huber estimator in the `2 -error. Left and right panel
show the negative log `2 -error versus δ in low and high dimensions, respectively.
Similarly, the regularization parameter is chosen as


log d min{δ/(1+δ),1/2}
.
λ = cλ ×
n
Here we choose cλ = 0.5 mainly for simplicity. In practice, cλ and the constant in robustification parameter can be chosen by cross validation. The left panel of Figure 2 plots
the error kβb − β ∗ k2 versus the sample size n, for problem dimensions d ∈ {100, 500, 5000}
over 200 repetitions. For all three choices, the `2 -error decreases to zero as the sample size
grows. The curves shift to the right when the dimension increases, reflecting the natural
intuition that higher dimensional problems are harder. Theorem 3.5 makes a specific prediction about this scaling behavior: if we plot the `2 -error versus the effective sample size
(n/ log d), the curves should align roughly with the theoretical curve

−δ/(1+δ)
n
∗
b
kβ − β k2 
log d

for different values of d. This is verified empirically by the right panel of Figure 2, as all
the curves now lie on top of one another. This alignment phenomenon in Figure 2 is also
observed by Wainwright (2009) for sub-Gaussian errors.
Our second example is designed to demonstrate the phase transition phenomenon of the
Huber estimator in both low and high dimensions. For the low dimensional case, we take
(n, d) = (500, 5) and a set of choices of degree of freedoms (df’s): df ∈ {1.1, 1.2, . . . , 3.0}.
For the high dimensional case, we take (n, d) = (500, 1000), with the same set of choices of
df’s. Tuning parameters (τ, λ) are set similarly as before. Predicted by Theorems 3.1, 3.2,
3.5 and 3.6, we have
16

(n = 500, d = 5)
●

Huber
LSE

●

2

1.5

0.5
0.4

2

2.0

(n = 500, d = 1000)

●

1.0

0.3

●
●

●

●

●

0.5

0.2

●

●
●

0.5

●

●

●
●
●

●

●

●

●

●

1.0

●

●

●

1.5

●

●

●

2.0

0.5

●

●

●

●

●

●

●

●

1.0

●

●

1.5

●

●

●

●

2.0

Figure 4: Comparison between the (sparse) adaptive Huber estimator and the (sparse)
OLS estimator in `2 -error: result in low dimensions is shown in left panel while right panel
concerns the high dimensional setting.
1. (Low dimensional setting):

− log kβbτ − β ∗ k2 

δ
1
log(n) −
log(vδ ), 0 < δ ≤ 1,
1+δ
1+δ

2. (High dimensional setting):


− log kβbτ,λ − β ∗ k2 

 n 
1
δ
log
−
log(vδ ), 0 < δ ≤ 1,
1+δ
log d
1+δ

which are approximately log(n) × δ/(1 + δ) and log(n/ log d) × δ/(1 + δ), respectively,
when n is large enough. Figure 3 shows the negative log `2 -error versus δ, in both low
and high dimensions over 200 repetitions for each (n, d) combination, roughly resembles
the shape of the curve in Figure 1. These simulations, once again, show good agreement
with the theoretical predictions and validate the phase transition of the Huber estimator.
Moreover, we compare the `2 -error of the adaptive Huber estimator with that of the OLS
estimator under t-distributed errors with varying degrees of freedoms. Figure 4 shows
that the performance of the adaptive Huber estimator always dominates that of the OLS
estimator especially when δ is small. When δ increases, the performances of the two methods
get closer.

6

Discussion

In this paper, we propose to use the adaptive Huber estimator for robust estimation and
inference. We sharply characterize the performance of the proposed method in terms of
17

nonasymptotic deviation inequalities in both low and high dimensions. A phase transition
phenomenon is then established by providing the matching lower bounds. Furthermore, a
nonasymptotic Bahadur representation for finite-sample inference is derived under bounded
variances.
The main focus of this paper is on the conditional mean regression in the presence of
heavy-tailed noise, which automatically distinguishes our method from quantile-based robust regressions [Belloni and Chernozhukov (2011); Zheng, Peng and He (2015)]. If recovery
of the conditional mean structure is the goal, quantile regression has to make the strong
assumption that the error-distributions are symmetric around zero. This is not necessary
in our work. Moreover, our work implies many of the previous results as special cases. For
instance, the mean estimation of a univariate probability distribution considered in Catoni
(2012) can be formulated as a regression problem X = µ + ε, where ε is the error variable
with finite (1 + δ)-th moment. Fan, Li and Wang (2017) studied the robust mean regression
in high dimensions under finite variance. Our result recovers theirs by taking δ = 1 in
Theorem 3.5. We highlight that we do not need the bounded parameter constraint in their
paper by exploiting the localized analysis developed in Fan et al. (2017).
Our techniques developed here can be potentially applied and generalized to other problems, such as robust reduced rank regression, robust low-rank matrix recovery and generic
robust M -estimation, among others. In all these problems, we conjecture that the adaptive
Huber-type estimators will achieve optimal rate of convergence in the presence of heavytailed errors. Another interesting question to ask is how the folded concave penalized
regression improves upon Lasso in terms of rate of convergence, when the signals are sparse
and large. We leave these to future work.

Supplementary Material
This supplementary material contains two sections. Section 1 contains proofs for Lemmas A.1–A.3 and Propositions 1 and 2 in the main text. In Section 2, we extend the
methodology to allow possibly heavy-tailed predictors.

A
A.1

Proofs
Auxiliary Lemmas

In this section, we collect several auxiliary lemmas. Our first lemma concerns the localized
analysis that can be utilized to remove the parameter constraint in previous works. It is
established in Fan et al. (2017) and we reproduce it here for completeness.
s (β , β ) =
Lemma A.1. Let DL (β1 , β2 ) = L(β1 ) − L(β2 ) − h∇L(β2 ), β1 − β2 i and DL
1
2
∗
∗
DL (β1 , β2 ) + DL (β2 , β1 ). For βη = β + η(β − β ) with η ∈ (0, 1] and any convex loss
functions L, we have
s
s
DL
(βη , β ∗ ) ≤ ηDL
(β, β ∗ ).

The following two lemmas provide some regularity properties of the Hessian matrix in
a local vicinity of the true parameter under both fix and random designs.
18

P
Lemma A.2. Assume that Condition 3.1 holds and that vδ = n−1 ni=1 E(|εi |1+δ ) < ∞
−1/2
−1/2
for some 0 < δ ≤ 1. Then, for any t, r > 0, the Hessian matrix Sn ∇2 Lτ (β)Sn
satisfies
−t
that, with probability greater than 1 − e ,

min
λmin Sn−1/2 ∇2 Lτ (β)Sn−1/2
1/2

β∈Rd :kSn (β−β ∗ )k2 ≤r

p

≥ 1 − 4τ −2 r2 d − (2/τ )1+δ vδ + t/(2n) d.

(A.1)

Lemma A.3. Assume that Condition 3.3 holds and vδ < ∞ for some 0 < δ ≤ 1. Then,
for any t, r > 0, the Hessian matrix ∇2 Lτ (β) satisfies that, with probability greater than
1 − 2e−t ,

min
λmin Σ−1/2 ∇2 Lτ (β)Σ−1/2
β∈Rd :kΣ1/2 (β−β ∗ )k2 ≤r

r
 √

 1+δ
r d
d+t
2
− C1
≥ 1 − vδ
+
τ
τ
n

(A.2)

as long as n ≥ C2 (d + t), where C1 , C2 > 0 are constants depending only on A0 .
The proofs of Lemmas A.1–A.3 can be found in the supplemental material.

A.2

Proof of Theorem 3.1

We start by constructing an intermediate estimator, denoted by βbη = β ∗ + η(βb − β ∗ ), such
1/2
1/2 b
∗ )k ≤
that kSn (βbη −β ∗ )k2 ≤ r for some r > 0 to be specified. We take η = 1 if kSn (β−β
2
1/2 b
∗
r; otherwise, we choose η ∈ (0, 1) so that kSn (βη − β )k2 = r. Applying Lemma A.1 gives
b − ∇Lτ (β ∗ ), βb − β ∗ ,
∇Lτ (βbη ) − ∇Lτ (β ∗ ), βbη − β ∗ ≤ η ∇Lτ (β)

b = 0 according to the Karush-Kuhn-Tucker condition. By the mean value
where ∇Lτ (β)
theorem, there exists a βeη , a convex combination of β ∗ and βbη , satisfying ∇Lτ (βbη ) −
∇Lτ (β ∗ ) = ∇2 Lτ (βeη )(βbη − β ∗ ). If, there exists some constant amin > 0 such that

λmin S−1/2
∇2 Lτ (βeη )Sn−1/2 ≥ amin ,
(A.3)
n

1/2
−1/2
1/2
then we have amin kSn (βbη − β ∗ )k22 ≤ kSn ∇Lτ (β ∗ )k2 kSn (βbη − β ∗ )k2 . Cancelling the
common factor on both sides yields

b
amin S1/2
n (βη − β)
−1/2

Define the vector ξ ∗ = Sn

2

≤ S−1/2
∇Lτ (β ∗ ) 2 .
n

(A.4)

∇Lτ (β ∗ ), which can be written as

n

1X 0
ei = (e
ξ =−
`τ (εi )e
xi with x
xi1 , . . . , x
eip )T = Sn−1/2 xi .
n
∗

i=1

√
P
Moreover, write Ψj = n−1 ni=1 x
eij ψ(εi /τ ) for j = 1, . . . , d such that kξ ∗ k2 ≤ d kξ ∗ k∞ =
√
d τ max1≤j≤d |Ψj |, where ψ(u) ≡ τ −1 `0τ (τ u), u ∈ R. It is easy to see that the function
ψ(·) satisfies
− log(1 − u + |u|1+δ ) ≤ ψ(u) ≤ log(1 + u + |u|1+δ )
19

(A.5)

for all u ∈ R. It follows that
x
eij ψ(εi /τ ) ≤ x
eij 1(e
xij ≥ 0) log(1 + εi /τ + |εi /τ |1+δ )

−x
eij 1(e
xij < 0) log(1 − εi /τ + |εi /τ |1+δ ).

This, together with the inequality (1 + u)v ≤ 1 + uv for u ≥ −1 and 0 < v ≤ 1, implies
exp{e
xij ψ(εi /τ )}
≤ (1 + εi /τ + |εi /τ |1+δ )xeij 1(exij ≥0) + (1 − εi /τ + |εi /τ |1+δ )−exij 1(exij <0)
≤ 1 + (εi /τ )e
xij + |εi /τ |1+δ .

Consequently, we have
E{exp(nΨj )} =

n
Y
i=1

n
Y
E exp{e
xij ψ(εi /τ )} ≤
(1 + vi,δ /τ 1+δ ) ≤ exp(vδ n/τ 1+δ ),
i=1

where we used the inequality 1 + u ≤ eu in the last step. For any z ≥ 0, using Markov’s
inequality gives
P(Ψj ≥ vδ z) ≤ exp(−vδ nz)E{exp(nΨj )} ≤ exp{vδ n(τ −1−δ − z)}.
As long as τ ≥ (2/z)1/(1+δ) , we have P(Ψj ≥ vδ z) ≤ e−vδ nz/2 . On the other hand, it can be
similarly shown that P(−Ψj ≥ vδ z) ≤ e−vδ nz/2 . For any t > 0, taking z = 2t/(vδ n) in these
two inequalities yields that as long as τ ≥ (vδ n/t)1/(1+δ) ,
√

P kξ ∗ k2 ≥ 2 d τ n−1 t
∗

≤ P kξ k∞ ≥ 2τ n

−1



t ≤

d
X
j=1


P |Ψj | ≥ 2n−1 t ≤ 2d exp(−t).

(A.6)

√
Taking r = 2δ · 4 d τ n−1 t. Then, it follows from Lemma A.2 and the definition of
τ that with probability greater than 1 − e−t , (A.3) holds with amin = 1/2 as long as
n ≥ 16 max(2d2 t, 2δ dt) = 32d2 t. Hence, combining (A.4) and (A.6) implies that, with
probability at least 1 − (2d + 1)e−t ,
√
bη − β ∗ ) ≤ 4 d τ n−1 t < r.
S1/2
(
β
n
2

By the definition of βbη in the beginning of the proof, we have η = 1 and thus βb = βη .

A.3

Proof of Theorem 3.2

We start by defining a simple class of distributions for the response variable y as Pc,γ =

Pc+ , Pc− , where
+
−
−
P+
c ({0}) = 1 − γ, Pc ({c}) = γ, and Pc ({0}) = 1 − γ, Pc ({−c}) = γ.

−
Here, we suppress the dependence of P+
c and Pc on γ for convenience. It follows that, for
−
any 0 < δ ≤ 1, the (1 + δ)-th absolute central moment vδ of y with law either P+
c or Pc is

vδ = |c|1+δ γ(1 − γ) γ δ + (1 − γ)δ .
(A.7)

20

For i = 1, . . . , n, let (y1i , y2i ) be independent pairs of real-valued random variables satisfying
−
P(y1i = y2i = 0) = 1 − γ, P(y1i = ci , y2i = −ci ) = γ, and y1i ∼ P+
ci , y2i ∼ Pci .

Let yk = (yk1 , . . . , ykn )T for k = 1, 2, and ξ ∈ (0, 1/2]. Taking γ = log{1/(2ξ)}/(2n) with
ξ ≥ e−n /2, we obtain 1 − γ ≥ 1/2 and
n
 −γ on

P y1 = y2 = 0 = (1 − γ)n ≥ exp
≥ 2ξ.
1−γ

By assumption, we know that there is an n-dimensional vector u ∈ {−1, +1}n with each
coordinate taking −1 or 1 such that n1 kXT ukmin ≥ α. Note that this assumption naturally
holds for the mean model, where X = (1, . . . , 1)T and α can be taken as 1. Now we take c,
β1∗ and β2∗ such that c = cu for a c > 0, Xβ1∗ = cγ and β2∗ = −β1∗ , which indicates that

β1∗

2

1
−1 1
β1∗ =
XT X
XT u, and
n
n
1
−1 1
√
√
≥ cγ
XT X
XT u ≥ cγ d kXT u/nkmin /cu ≥ cγα d/cu .
n
n
2

Let βbk (yk ) be any estimator possibly depending on ξ, then the above calculation yields
n
√
√

o
max P βb1 − β1∗ 2 ≥ cγα d/cu , P βb2 − β2∗ 2 ≥ cγα d/cu

√
√
1 
≥ P βb1 − β1∗ 2 ≥ cγα d/cu or βb2 − β2∗ 2 ≥ cγα d/cu
2
 1
 1
1 b
(A.8)
≥ P β1 = βb2 ≥ P y1 = y2 ≥ (1 − γ)n ≥ ξ,
2
2
2

where we suppress the dependence of βbk on yk for simplicity. Using the fact that cγ ≥
1/(1+δ)
vδ
(γ/2)δ/(1+δ) further implies
"
√ 
δ/(1+δ) #
d
log(1/(2ξ))
α
1/(1+δ)
P βb1 − β1∗ 2 ≥ vδ
cu
2n
"
√ 
δ/(1+δ) #
_
1/(1+δ) α d log(1/(2ξ))
∗
b
≥ ξ.
P β2 − β2 2 ≥ v δ
cu
2n

Now since Pc,γ ⊆ Pδvδ , taking log{1/(2ξ)} = 2t implies the result for the case where δ ∈ (0, 1].
When δ > 1, the second moment exists. Thus using v1 < ∞ completes the proof.

A.4

Proof of Theorem 3.3

To begin with, we introduce the following notations. Define the pseudo-likelihood process
L(β) = −Lτ (β) such that βb = argmaxβ∈Rd L(β), and write D2 (β) = −∇2 EL(β). It follows
from elementary calculus that
n

D2 (β) = Sn −


1X
T
P |yi − xT
i β| > τ xi xi .
n
i=1

21

1/2

Moreover, define the parameter set Θ0 (r) = {β ∈ Rd : Sn (β − β ∗ ) 2 ≤ r}, r > 0. Taking
r1 = 4τ0 (dt/n)1/2 , it follows from Theorem 3.1 with δ = 1 that as long as n ≥ 32d2 t,

P βb ∈ Θ0 (r1 ) ≥ 1 − (2d + 1) exp(−t).
(A.9)
Throughout the following, we restrict our attention to the local vicinity Θ0 (r1 ) of β ∗ , and
use Proposition 3.1 in Spokoiny (2013) to prove (3.4). To this end, we need to verify
Conditions (L0 ) and (ED2 ) therein.
1/2

Condition (L0 ): For every β ∈ Θ0 (r) and u ∈ Sd−1 , write δ = Sn (β − β ∗ ) and note that
n
 T 2
1X
u
− Id u ≤
P |yi − xT
xi u)
i β| > τ (e
n
i=1


n
 T 2
1 X
1
T
2
2
T 2
≤
xi k2 .
vi,1 + |e
xi δ| (e
xi u) ≤ 2 u ∆n u + r max ke
1≤i≤n
nτ 2
τ
T



S−1/2
D2 (β)S−1/2
n
n

i=1

Under Condition 3.1, this implies that
max

β∈Θ0 (r)


S−1/2
D2 (β)Sn−1/2 − Id ≤ δ(r) ≡ τ −2 k∆n k + r2 d ,
n

(A.10)

which verifies Condition (L0 ).

Condition (ED2 ): Define ζ(β) = L(β) − EL(β) and note that
n

∇2 ζ(β) = −



1 X
T
1 |yi − xT
xi xT
i β| ≤ τ − P |yi − xi β| ≤ τ
i .
n
i=1

1/2

1/2

1/2

1/2

For every u, v ∈ Rd \ {0}, write u1 = Sn u/kSn uk2 and v1 = Sn v/kSn vk2 . Under
√
Condition 3.2, we have for all λ ∈ R satisfying |λ| ≤ κ n that


λ
uT ∇2 ζ(β)v
E exp
1/2
n−1/2 kS1/2
n uk2 kSn vk2

n 
Y
|λ|
λ2 T
eT
xT
2 T
2 √n |e
i u1 x
i v1 |
(e
x u1 ) (e
xi v1 ) e
≤
1+
2n i
i=1

n 
Y
λ
λ
2
2
λ2 T
√
xT
(e
xT
2 T
2 2√n (e
i u1 )
i v1 )
2
n
(e
x u1 ) (e
xi v1 ) e
≤
1+
e
2n i
i=1
 2

n
λ
1 X T 4 κ(exT u)2
≤ exp
max
(e
xi u) e i
≤ exp(M λ2 /2).
2 u∈Sd−1 n
i=1

√
Therefore, Condition (ED2 ) holds with ω = n−1/2 , ν0 = M 1/2 and g(r) = κ n.
1/2
With the above preparations, using Proposition 3.1 in Spokoiny (2013) with D0 = Sn
we have, for any 0 < t ≤ 12 κ2 n − 2d,


∆(r) ≡ sup S1/2
β − β ∗ + Sn−1/2 ∇L(β) − ∇L(β ∗ ) 2
n
β∈Θ0 (r)

≤ δ(r)r + 6M 1/2 (4d + 2t)1/2 n−1/2 r

with probability greater than 1 − e−t , where δ(r) is given in (A.10). Taking r = r1 in the
b = 0 prove (3.4).
above inequality, using (A.9) and the fact that ∇L(β)
22

A.5

Proof of Theorem 3.4

Proof of (3.5). In view of (A.3) and (A.4), lying in the heart of the arguments is to derive
deviation inequalities for kΣ−1/2 ∇Lτ (β ∗ )k2 under the moment condition that vδ < ∞ for
some 0 < δ ≤ 1.
Case 1. First, we deal with the case δ ∈ (0, 1) using arguments similar to those employed
ei = (e
in the proof of Theorem 3.1. Define vectors x
xi1 , . . . , x
eid )T = Σ−1/2 xi for i = 1, . . . , n,
and zj = (e
x1j , . . . , x
enj )T , zej = (e
zj1 , . . . , zejn ) = zj /kzj k∞ for j = 1, . . . , d. For any t > 0, it
follows that


√
−1/2
∗
−1
P Σ
∇Lτ (β ) 2 ≥ 2 d max kzj k∞ τ n t x1 , . . . , xn
1≤j≤d


−1/2
∗
−1
≤P Σ
∇Lτ (β ) ∞ ≥ 2 max kzj k∞ τ n t x1 , . . . , xn
1≤j≤d

 X

d
n
X
0
≤
P
zeij `τ (εi ) ≥ 2τ t x1 , . . . , xn .
j=1

(A.11)

i=1

Next, we define the function ψ(u) = τ −1 `0τ (τ u), u ∈ R. Again, by (A.5) we have
zeij ψ(εi /τ ) ≤ zeij 1(e
zij ≥ 0) log(1 + εi /τ + |εi /τ |1+δ )

− zeij 1(e
zij < 0) log(1 − εi /τ + |εi /τ |1+δ ),

which, together with the inequality (1 + u)v ≤ 1 + uv for all u ≥ −1 and 0 < v ≤ 1, yields
exp{e
zij ψ(εi /τ )}
≤ (1 + εi /τ + |εi /τ |1+δ )zeij 1(ezij ≥0) + (1 − εi /τ + |εi /τ |1+δ )−ezij 1(ezij <0)

≤ 1 + (εi /τ )e
zij + |εi /τ |1+δ .

Taking conditional expectations gives
 Pn
E e i=1 zeij ψ(εi /τ ) |x1 , . . . , xn
n
n
Y
Y


=
E ezeij ψ(εi /τ ) |x1 , . . . , xn ≤
1 + vi,δ τ −1−δ ≤ exp(vδ nτ −1−δ ).
i=1

i=1

By (A.11) and Markov’s inequality, it follows that for each j,
X

n
0
P
zeij `τ (εi ) ≥ 2τ t x1 , . . . , xn ≤ exp(−t),
i=1

P
almost surely. Similarly, it can be shown that P{− ni=1 zeij `0τ (εi ) ≥ 2τ t x1 , . . . , xn } ≤ e−t .
Therefore, we obtain that given x1 , . . . , xn , with conditional probability at least 1 − 2e−t ,
√
xi k∞ τ n−1 t
Σ−1/2 ∇Lτ (β ∗ ) 2 ≤ 2 d max ke
1≤i≤n

almost surely. Under Condition 3.3, we have ke
xi kψ2 ≤ A0 . For any z > 0, by the union
bound we get



P max ke
xi k∞ ≥ z ≤ nd
max
P |e
xij | ≥ z ≤ exp{log(2nd) − cz 2 }.
1≤i≤n

1≤i≤n,1≤j≤d

23

p
By taking z = c−1/2 log(2nd) + t, the last two displays jointly imply that with probability
greater than 1 − (2d + 1)e−t ,
√
Σ−1/2 ∇Lτ (β ∗ ) 2 ≤ 2c−1/2 {log(2nd) + t}1/2 d τ n−1 t.
(A.12)
Back to (A.3), recall that kΣ1/2 (βeη − β ∗ )k2 ≤ r. Then, it follows from Lemma A.3 by
√
taking r = 4.1c−1/2 {log(2nd) + t}1/2 d τ n−1 t that, with probability greater than 1 − 2e−t ,
 1
λmin Σ−1/2 ∇2 Lτ (β)Σ−1/2 ≥
2
β∈Rd :kΣ1/2 (β−β ∗ )k2 ≤r
min

as long as n ≥ C1 max{d + t, (log n + t)dt}, where C1 > 0 is a constant depending only on
A0 . This, together with (A.3) and (A.4) with Sn replaced by Σ, and (A.12), implies that
kΣ1/2 (βbη − β ∗ )k2 ≤ r with probability at least 1 − (2d + 3)e−t . By the definition of βbη in
b
the beginning of the proof, this indicates that η = 1, and thus the bound (3.5) holds for β.
Proof of (3.6). Next we focus on the finite variance case. Define the random vector
n
√
1 X 0
−1/2
∗
√
ξ = − nΣ
∇Lτ (β ) =
`τ (εi )e
xi .
n
∗

i=1

For any λ ∈ R and u ∈ Sd−1 , by the inequalities eu ≤ 1 + u + u2 et∨0 /2 and 1 + u ≤ eu that
hold for any u ∈ R, we have

n
n 
Y
Y
λ2 2  T 2 √|λ|n τ |uT xei |
T ∗
−1/2 0
T
ei ) e
ei } ≤
v E (u x
E exp(λu ξ ) =
E exp{λn
`τ (εi )u x
1+
.
2n i,1
i=1

Recall that τ = τ0

i=1

p
p
1/2
n/(d + t) with τ0 ≥ v1 . Hence, if |λ| ≤ 2(d + t), we have
√
|λ|
√ τ ≤ 2τ0 and E exp(λuT ξ ∗ ) ≤ exp(ν02 v1 λ2 /2),
n

where ν0 ≥ 1 is a constant depending only on τ0 and A0 . This verifies condition (1.18) in
p
the supplement of Spokoiny (2012) with V0 = Id and g = 2(d + t). By Corollary 1.13
√
there, we obtain that for any 2d/18 < x ≤ xc ,

P kξ ∗ k22 > ν02 v1 (d + 6x) ≤ 2 exp(−x) + 8.4 exp(−xc ),
√
where xc = (1 − log 3 )d + 1.5t. Taking x = d/3 + t in the last display yields that with
probability at least 1 − 7e−t ,
r
3d + 6t
1/2
−1/2
∗
Σ
∇Lτ (β ) 2 ≤ r0 ≡ ν0 v1
.
n
This, together with Lemma A.3 with r = 2.1r0 , implies that with probability at least
√
1 − 9e−t , kΣ1/2 (βbη − β ∗ )k2 ≤ 2r0 < r whenever n ≥ C2 d (d + t), where C2 > 0 is a
constant depending only on (τ0 , A0 ). Following the same arguments as we used in the proof
of Theorem 3.1, this proves (3.6).
Proof of (3.7). On the basis of (3.6), throughout the following we restrict our attention
to the local neighborhood of β ∗ , so that a variant of Proposition 3.1 in Spokoiny (2013)
24

can be used. To this end, we need to verify Conditions (L0 ) and (ED2 ) there. For β ∈ Rd ,
define L(β) = −Lτ (β) such that βb = argmaxβ∈Rd L(β), and write
n

1X 
T
E 1 |yi − xT
D (β) = −∇ EL(β) = Σ −
i β| > τ xi xi .
n
2

2

i=1

Moreover, define the parameter set Θ0 (r) = {β ∈ Rd : kΣ1/2 (β − β ∗ )k2 ≤ r}, r > 0.
p
For any t > 0, it follows from (3.6) that with τ = τ0 n/(d + t), βb = βbτ ∈ Θ0 (r1 ) holds
1/2 p
with probability at least 1 − 9e−t for some r1 > 0 that is proportional to v1
(d + t)/n.
Therefore, we only need to focus on the local vicinity Θ0 (r1 ) of β ∗ .
Condition (L0 ): For every β ∈ Θ0 (r) and u ∈ Sd−1 , we write δ = Σ1/2 (β − β ∗ ) such that
kδk2 ≤ r. Then, it follows that
n

 T 2
1X 
ei )
uT Σ−1/2 D2 (β)Σ−1/2 − Id u =
E 1 |yi − xT
i β| > τ (u x
n
i=1

n
1 X
ei )2 (uT x
ei )2 ≤ v1 τ −2 + 16A40 τ −2 kδk22 ≤ v1 τ −2 + 16A40 τ −2 r2 .
≤
vi,1 + E(δ T x
nτ 2
i=1

This verifies Condition (L0 ) in Spokoiny (2013) with δ(r) = v1 τ −2 + 16A40 τ −2 r2 and D0 =
1/2
Σ1/2 . We let (r, τ ) satisfy τ ≥ max(2v1 , 8A20 r), so that δ(r) ≤ 1/2 in this range.
Condition (ED2 ): Define ζ(β) = L(β) − EL(β), β ∈ Rd . For every u, v ∈ Rd , write
u1 = Σ1/2 u/kΣ1/2 uk2 and v1 = Σ1/2 v/kΣ1/2 vk2 . For any random variable X, we use
(Id − E)X to denote its centered version X − EX. Then, for any λ ∈ R, it follows from the
inequality eu ≤ 1 + u + u2 eu∨0 /2 that


√
uT ∇2 ζ(β)v
E exp λ n 1/2
kΣ uk2 kΣ1/2 vk2
n
Y

 T
e i u1 x
eT
=
E exp − λn−1/2 (Id − E)1 |yi − xT
i β| ≤ τ x
i v1
i=1

≤

≤

n 
Y

1+

n 
Y

1+

i=1

i=1

|λ|
λ2  T
eT
eT
xT
xT
2 √n (|e
i u1 x
i v1 |+E|e
i u1 x
i v1 |)
E (e
xi u1 )2 (e
xT
i v1 ) e
n

λ2  √|λ| (|exT u xeT v |+E|exTi u1 xeTi v1 |)
+ E e n i 1 i 1
n

|λ|
λ2 √|λ|n  T
eT
xT
2 √n |e
i u1 x
i v1 |
e E (e
xi u1 )2 (e
xT
i v1 ) e
n




λ2 √|λ|n  √|λ|n |exTi u1 xeTi v1 |
.
+ e E e
n
√
In particular, for λ ∈ R satisfying |λ| ≤ c0 n with c0 > 0 being a constant depending only
on A0 , using the inequality 1 + u ≤ eu we get


√
uT ∇2 ζ(β)v
sup E exp λ n 1/2
≤ exp(ν02 λ2 /2),
1/2 vk
kΣ
uk
kΣ
d
2
2
u,v∈R
25

where ν0 > 0 is a constant depending only on A0 . Therefore, Condition (ED2 ) holds with
√
ω = n−1/2 and g(r) = c0 n.
Consequently, using Proposition 3.1 in Spokoiny (2013) with D0 replaced by Σ1/2 we
derive that, for any 0 < t ≤ (c20 /2)n − 2d,
sup
β∈Θ0 (r)

Σ1/2 (β − β ∗ ) + Σ−1/2 {∇L(β) − ∇L(β ∗ )}

2


≤ δ(r)r + 6ν0 (4d + 2t)1/2 n−1/2 r

with probability greater than 1 − e−t . Recall that βb ∈ Θ0 (r1 ) with probability greater than
b = 0. Therefore, taking r = r1 in the last display proves (3.7).
1 − 9e−t , and ∇L(β)

A.6

Proof of Theorem 3.5

We start with the proof of Lemma 3.1.
Proof of Lemma 3.1. Let Hτ = ∇2 Lτ (β), where we suppress the dependence on β. Then
for any (u, β) ∈ C(k, γ, r), we have

 X
n

1
T
T
xi xi 1 |yi − xi β| ≤ τ u
u Hτ u = u
n
T

T

i=1

n
n
X

1X T 2
T
2
T
∗
(xi u) 1 |xi (β−β )| ≥ τ /2)−
(xi u) 1 |εi | > τ /2
n

≥

2 1
S1/2
n u 2−
n

≥

2 2r
S1/2
max kxi k∞ kSn1/2 uk22 − max
n u 2−
1≤i≤n
τ 1≤i≤n

i=1

i=1

2
xT
i u

n

1X
1 |εi | > τ /2 .
n

(A.13)

i=1

As kxi k∞ ≤ 1 for any 1 ≤ i ≤ n, we have

√
|xT
i u| ≤ kxi k∞ kuk1 ≤ (1 + γ)kxi k∞ kuJ k1 ≤ (1 + γ) k.

Moreover, for any t ≥ 0, applying Hoeffding’s inequality implies that with probability at
least 1 − e−t ,
r
r
 1+δ X
 1+δ
n
n

1X
2
1
t
2
t
1 |εi | > τ /2 ≤
vi,δ +
=
vδ +
.
n
τ
n
2n
τ
2n
i=1

i=1

Putting the above calculations together, we obtain
uT Hτ u ≥ S1/2
n u

2
2

− 2τ −1 r Sn1/2 u

2
2

− k(1 + γ)2 21+δ vδ τ −1−δ +

Consequently, as long as τ ≥ 8r, the following inequality

p

t/2 n−1/2 .

p
 1
3
uT Hτ u ≥ κl − k(1 + γ)2 21+δ vδ τ −1−δ + t/2 n−1/2 ≥ κl ,
4
2

(A.14)

holds uniformly over (u, β) ∈ C(k, γ, r) with probability at least 1 − e−t , where the last
inequality in (A.14) holds whenever τ & (vδ k)1/(1+δ) and n & k 2 t. On the other side, it can
be easily shown that uT Hτ u ≤ κu . This completes the proof of the lemma.
26

The following lemma is taken from Fan et al. (2017), showing that the solution βb falls
in a `1 -cone. It is critical for our analysis.
Lemma A.4 (`1 -cone Property). For any E such that S ⊆ E, if k∇Lτ (β ∗ )k∞ ≤ λ/2, we
must have: k(βb − β ∗ )E c k1 ≤ 3k(βb − β ∗ )E k1 .
We are ready to prove the main theorem of this section.

Proof of Theorem 3.5. It suffices to prove the statement for δ ∈ (0, 1]. We follow the proof
scheme of Lemma B.4 in Fan et al. (2017). We start by constructing an intermediate
estimator βbη = β ∗ + η(βb − β ∗ ) such that kβbη − β ∗ k1 ≤ r for some r > 0 to be specified. We
take η = 1 if kβb − β ∗ k1 ≤ r, and choose η ∈ (0, 1) so that kβbη − β ∗ k1 = r otherwise. Lemma
A.4 implies that βbη falls in a `1 -cone, and thus βbη is in the corresponding local `1 -cone:
k(βbη − β ∗ )S c k1 ≤ 3k(βbη − β ∗ )S k1 and kβbη − β ∗ k1 ≤ r.

(A.15)

Under Condition 3.5, it follows from Lemma 3.1 that with probability at least 1 − e−t ,
κl b
kβη − β ∗ k22 ≤ ∇Lτ (βbη ) − ∇Lτ (β ∗ ), βbη − β ∗
2

as long as τ & (vδ k)1/(1+δ) ∨ r and n & k 2 t. Applying Lemma A.1 and following the same
calculations as in Lemma B.4 of Fan et al. (2017), we obtain
 √
κl b
kβη − β ∗ k22 ≤ λ s + k∇Lτ (β ∗ )S k2 k(βbη − β ∗ )S k2 ,
2

which, combined with k∇Lτ (β ∗ )S k∞ ≤ λ/2, implies that

√
kβbη − β ∗ k2 ≤ 3κ−1
l λ s.

(A.16)

√
Inequalities in (A.15) imply that kβbη − β ∗ k1 ≤ 4k(βbη − β ∗ )S k1 ≤ 4 skβbη − β ∗ k2 ≤
b
b
b
12κ−1
l λs < r. By the construction of βη , we obtain that βη = β, and thus the desired result
holds. It remains to bound the probability of the event {k∇Lτ (β ∗ )S k∞ ≤ λ/2}. Recall
P
that the gradient of Lτ at β ∗ can be written as ∇Lτ (β ∗ ) = n−1 ni=1 `0τ (εi )xi . Following
the same arguments as in the proof of Theorem 3.1 and taking τ = τ0 (n/t)1/(1+δ) for some
τ0 ≥ νδ , we obtain that P{k(∇Lτ (β ∗ ))S k∞ ≥ 2τ n−1 t} ≤ 2se−t . This, together with (A.16),
proves (3.9). Finally, taking t = (1 + A) log d for any A > 0 yields P[k∇Lτ (β ∗ )k∞ ≥
2τ0 {(1 + A)(log d)/n}δ/(1+δ) ] ≤ (2s + 1)d−1−A . As implied by Condition 3.5 with k = 2s,
we have 2s + 1 ≤ d and thus (3.10) follows immediately.

A.7

Proof of Theorem 3.6

The proof of this theorem follows the similar argument to that of Theorem 3.2. It suffices
to prove the result for δ ∈ (0, 1]. Similar to the proof of Theorem 3.2, We start by defining
a simple class of distributions for the response variable y as Pc,γ = {Pc+ , Pc− }, where
+
−
−
P+
c ({0}) = 1 − γ, Pc ({c}) = γ, and Pc ({0}) = 1 − γ, Pc ({−c}) = γ.

27

−
Here, we suppress the dependence of P+
c and Pc on γ for convenience. It follows that, for
−
any 0 < δ ≤ 1, the (1 + δ)-th absolute central moment vδ of y with law either P+
c or Pc is

vδ = |c|1+δ γ(1 − γ){γ δ + (1 − γ)δ }.

(A.17)

We define the following s-sparse sign-ball Un as


Un = u : u ∈ {−1, 1}n .

By assumption, there exist u ∈ Un and A with |A| = s such that kXT
A ukmin /n ≥ α. Take
∗
∗
n
∗
β1 , β2 supported on A and c ∈ R such that c = cu for a c > 0, Xβ1 = cγ and β2∗ = −β1∗ .
Let P+ be the distribution of y1 = Xβ1∗ + ε and P− that of y2 = Xβ2∗ + ε. Clearly, we have
E(εi ) = 0 and E(|εi |1+δ ) = c1+δ γ(1 − γ){γ δ + (1 − γ)δ }.
Let A be the support of β1∗ . Then, we have

−1 1
1
XT u, and
(β1∗ )A = cγ XT
A XA
n
n A
√
√
kβ1∗ k2 ≥ cγ skXT
A u/nkmin /κu ≥ cγ sα/κu .

Let βbk (yk ) be any s-sparse estimator. With the above setup, we have
n

o
√
√
max P βb1 − β1∗ 2 ≥ cγ sα/κu , P βb2 − β2∗ 2 ≥ cγ sα/κu

√
√
1 
2
≥ P βb1 − β1∗ 2 ≥ cγ sα/κu or βb2 − β2∗ 2 ≥ cγ sα/κu
2

 1
1 b
≥ P β1 = βb2 ≥ P y1 = y2 = 0
2
2

(A.18)

where we suppress the dependence of βbk on yk for simplicity. For the last quantity in

the displayed inequality above, taking γ = log{1/(2t)} (2n) with t ≥ e−n /2, we obtain
1 − γ ≥ 1/2 and
n
 −γ on

P y1 = y2 = 0 = (1 − γ)n ≥ exp
≥ 2t.
1−γ
1/(1+δ)

Using the fact that cγ ≥ vδ
(γ/2)δ/(1+δ) , this further implies
"

δ/(1+δ) #
√
log{1/(2t)}
1/(1+δ)
P βb1 − β1∗ 2 ≥ vδ
α d/κu
2n
"

δ/(1+δ) #
_
√
log{1/(2t)}
1/(1+δ)
α d/κu
≥ t.
P βb2 − β2∗ 2 ≥ vδ
2n

Now since Pc,γ ⊆ Pδvδ , taking t = 2−1 d−A implies the result for the case where δ ∈ (0, 1].
When δ > 1, the second moment exists. Thus using v1 < ∞ completes the proof.

28

References
Belloni, A. and Chernozhukov, V. (2011). `1 -penalized quantile regression in highdimensional sparse models. Ann. Statist. 39 82–130.
Bickel, P. J., Ritov, Y. and Tsybakov, A. B. (2009). Simultaneous analysis of Lasso
and Dantzig selector. Ann. Statist. 37 1705–1732.
Brownlees, C., Joly, E. and Lugosi, G. (2015). Empirical risk minimization for
heavy-tailed losses. Ann. Statist. 43 2507–2536.
Bühlmann, P. and van de Geer, S. (2011). Statistics for High-Dimensional Data:
Methods, Theory and Applications. Springer, Heidelberg.
Catoni, O. (2012). Challenging the empirical mean and empirical variance: A deviation
study. Ann. Inst. Henri Poincaré Probab. Stat. 48 1148–1185.
Catoni, O. (2016). PAC-Bayesian bounds for the Gram matrix and least squares
regression with a random design. Preprint. Available at arXiv:1603.05229.
Delaigle, A., Hall, P. and Jin, J. (2011). Robustness and accuracy of methods for
high dimensional data analysis based on Student’s t-statistic. J. R. Stat. Soc. Ser. B.
Stat. Methodol. 73 283–301.
Devroye, L., Lerasle, M., Lugosi, G. and Oliveira, R. I. (2016). Sub-Gaussian
mean estimators. Ann. Statist. 44 2695–2725.
Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004).
regression. Ann. Statist. 32 407–499.

Least angle

Eklund, A., Nichols, T. and Knutsson, H. (2016). Cluster failure: Why fMRI
inferences for spatial extent have inflated false-positive rates. Proc. Natn. Acad. Sci.
USA 113 7900–7905.
Fan, J., Li, Q. and Wang, Y. (2017). Estimation of high dimensional mean regression
in the absence of symmetry and light tail assumptions. J. R. Stat. Soc. Ser. B. Stat.
Methodol. 79 247–265.
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its
oracle properties. J. Amer. Statist. Assoc. 96 1348–1360.
Fan, J., Liu, H., Sun, Q. and Zhang, T. (2017). I-LAMM for sparse learning:
Simultaneous control of algorithmic complexity and statistical error. Ann. Statist. To
appear. Available at arXiv:1507.01037.
Giulini, I. (2016). Robust principal component analysis in Hilbert spaces. Preprint.
Available at arXiv:1606.00187.
Hastie, T., Tibshirani, R. and Wainwright, M. J. (2015). Statistical Learning with
Sparsity: The Lasso and Generalizations. CRC Press.
He, X. and Shao, Q.-M. (1996). A general Bahadur representation of M -estimators and
its application to linear regression with nonstochastic designs. Ann. Statist. 24 2608–2630.
He, X. and Shao, Q.-M. (2000). On parameters of increasing dimensions. J. Multivariate
Anal. 73 120–135.
29

Huber, P. J. (1964). Robust estimation of a location parameter. Ann. Math. Statist. 35
73–101.
Huber, P. J. (1973). Robust regression: Asymptotics, conjectures and Monte Carlo. Ann.
Statist. 1 799–821.
Loh, P. and Wainwright, M. J. (2015). Regularized M -estimators with nonconvexity: Statistical and algorithmic theory for local optima. J. Mach. Learn. Res. 16 559–616.
Mammen, E. (1989). Asymptotics with increasing dimension for robust regression with
applications to the bootstrap. Ann. Statist. 17 382–400.
Minsker, S. (2016). Sub-Gaussian estimators of the mean of a random matrix with
heavy-tailed entries. Preprint. Available at arXiv:1605.07129.
Lepski, O. V. (1991). Asymptotically minimax adaptive estimation. I. Upper bounds.
Optimally adaptive estimates. Theory Probab. Appl. 36 682–697.
Portnoy, S. (1985). Asymptotic behavior of M estimators of p regression parameters
when p2 /n is large; II. Normal approximation. Ann. Statist. 13 1403–1417.
Spokoiny, V. (2012). Parametric estimation. Finite sample theory. Ann. Statist. 40
2877–2909.
Spokoiny, V. (2013). Bernstein–von Mises theorem for growing parameter dimension.
Preprint. Available at arXiv:1302.3430.
Sun, Q., Zhou, W.-X. and Fan, J. (2017). Supplement to “Adaptive Huber regression:
Optimality and phase transition.”
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. R. Stat. Soc.
Ser. B. Stat. Methodol. 58 267–288.
Vershynin, R. (2012). Introduction to the nonasymptotic analysis of random matrices.
In Compressed Sensing 210–268. Cambridge University Press, Cambridge.
Wainwright, M. J. (2009). Sharp thresholds for high-dimensional and noisy sparsity
recovery using `1 -constrained quadratic programming (lasso). IEEE Trans. Inform.
Theory 55 2183–2202.
Wang, L., Peng, B. and Li, R. (2015). A high-dimensional nonparametric multivariate
test for mean vector. J. Amer. Statist. Assoc. 110 1658–1669.
Yohai, V. J. and Maronna, R. A. (1979). Asymptotic behavior of M -estimators for the
linear model. Ann. Statist. 7 258–268.
Zheng, Q., Peng, L. and He, X. (2015). Globally adaptive quantile regression with
ultra-high dimensional data. Ann. Statist. 43 2225–2258.
Zhou, W.-X., Bose, K., Fan, J. and Liu, H. (2016). A new perspective on robust
M -estimation: Finite sample theory and applications to dependence-adjusted multiple
testing. Preprint.

30

