A Shrinkage Principle for Heavy-Tailed Data:
High-Dimensional Robust Low-Rank Matrix Recovery∗
Jianqing Fan, Weichen Wang, Ziwei Zhu

arXiv:1603.08315v2 [math.ST] 4 May 2017

Department of Operations Research and Financial Engineering
Princeton University.

Abstract
This paper introduces a simple principle for robust high-dimensional statistical
inference via an appropriate shrinkage on the data. This widens the scope of highdimensional techniques, reducing the moment conditions from sub-exponential or subGaussian distributions to merely bounded second or fourth moment. As an illustration
of this principle, we focus on robust estimation of the low-rank matrix Θ∗ from the trace
regression model Y = Tr(Θ∗T X) + . It encompasses four popular problems: sparse
linear models, compressed sensing, matrix completion and multi-task regression. We
propose to apply penalized least-squares approach to appropriately truncated or shrunk
data. Under only bounded 2 + δ moment condition on the response, the proposed
robust methodology yields an estimator that possesses the same statistical error rates
as previous literature with sub-Gaussian errors. For sparse linear models and multitasking regression, we further allow the design to have only bounded fourth moment
and obtain the same statistical rates, again, by appropriate shrinkage of the design
matrix. As a byproduct, we give a robust covariance matrix estimator and establish its
concentration inequality in terms of the spectral norm when the random samples have
only bounded fourth moment. Extensive simulations have been carried out to support
our theories.

Keywords: Robust Statistics, Shrinkage, Heavy-Tailed Data, Trace Regression, LowRank Matrix Recovery, High-Dimensional Statistics.
∗

The research was supported by NSF grants DMS-1206464 and DMS-1406266 and NIH grant R01GM072611-12.

1

1

Introduction

Heavy-tailed distributions are ubiquitous in modern statistical analysis and machine learning
problems. They are stylized features of high-dimensional data. By chance alone, some of
observable variables in high-dimensional datasets can have heavy or moderately heavy tails
(see right panel of Figure 1). It has been widely known that financial returns and macroeconomic variables exhibit heavy tails, and large-scale imaging datasets in biological studies
are often corrupted by heavy-tailed noises due to limited measurement precisions. Figure 1
provides some empirical evidence on this which is pandemic to high-dimensional data. These
stylized features and phenomena contradict the popular assumption of sub-Gaussian or subexponential noises in the theoretical analysis of standard statistical procedures. They also
have adverse impacts on the methods that are popularly used. Simple and effective principles
are needed for dealing with moderately heavy or heavy tailed data.
Distribution of Kurtosis

a

Frequency
15

0

5

20

10

40

20

60

25

80

Histgram of Kurtosis

0

6 10 15 20 25 30
0

t5 10

20

30

40

50

36 41

48 53

60 65

71

83

96

kurtosis

Figure 1: Distributions of kurtosis of macroeconomic variables and gene expressions.
Red dashline marks variables with empirical kurtosis equals to that of t5 -distribution. Left panel:
For 131 macroeconomics variables in Stock and Watson (2002). Right panel: For logarithm of
expression profiles of 383 genes based on RNA-seq for autism data (Gupta et al., 2014), whose
kurtosis is bigger than that of t5 among 19122 genes.

Recent years have witnessed increasing literature on the robust mean estimation when
the population distribution is heavy-tailed. Catoni (2012) proposed a novel approach that
is through minimizing a robust empirical loss. Unlike the traditional `2 loss, the robust loss
function therein penalizes large deviations, thereby making the correspondent M-estimator
insensitive to extreme values. It turns out that when the population has only finite second
moment, the estimator has exponential concentration around the true mean and enjoys the
same rate of statistical consistency as the sample average for sub-Gaussian distributions.
2

Brownlees et al. (2015) pursued the Catoni’s mean estimator further by applying it to empirical risk minimization. Fan et al. (2016) utilized the Huber loss with diverging threshold,
called robust approximation to quadratic (RA-quadratic), in a sparse regression problem
and showed that the derived M-estimator can also achieve the minimax statistical error rate.
Loh (2015) studied the statistical consistency and asymptotic normality of a general robust
M -estimator and provided a set of sufficient conditions to achieve the minimax rate in the
high-dimensional regression problem.
Another effective approach to handle heavy-tailed distribution is the so-called “median
of means” approach, which can be traced back to Nemirovsky et al. (1982). The main idea
is to first divide the whole samples into several parts and take the median of the means
from all pieces of sub-samples as the final estimator. This “median of means” estimator also
enjoys exponential large deviation bound around the true mean. Hsu and Sabato (2016)
and Minsker (2015) generalized this idea to multivariate cases and applied it to robust PCA,
high-dimensional sparse regression and matrix regression, achieving minimax optimal rates
up to logarithmic factors.
In this paper, we propose a simple and effective principle: truncation of univariate data
and more generally shrinkage of multivariate data to achieve the robustness. We will illustrate our ideas through a general model called the trace regression
Y = Tr(Θ∗T X) + ,
which embraces linear regression, matrix or vector compressed sensing, matrix completion
and multi-tasking regression as specific examples. The goal is to estimate the coefficient
matrix Θ∗ ∈ Rd1 ×d2 , which is assumed to have a nearly low-rank structure in the sense
min(d
P1 ,d2 )
that its Schatten norm is constrained:
σi (Θ∗ )q ≤ ρ for 0 ≤ q < 1, where σi (Θ∗ ) is
i=1

the ith singular value of Θ∗ , i.e., the square-root of the ith eigenvalue of Θ∗ T Θ∗ . In other
words, the singular values of Θ∗ decay fast enough so that Θ∗ can be well approximated by
a low-rank matrix. We always consider the high-dimensional setting where the sample size
n  d1 d2 . As we shall see, appropriate data shrinkage allows us to recover Θ∗ with only
bounded moment conditions on noise and design.
As the most simple and important example of low-rank trace regression, sparse linear
regression and compressed sensing have become a hot topic in statistics research in the past
3

two decades. See, for example, Tibshirani (1996), Chen et al. (2001), Fan and Li (2001),
Donoho (2006), Candes and Tao (2006), Candes and Tao (2007), Candes (2008), Nowak
et al. (2007), Fan and Lv (2008), Zou and Li (2008), Bickel et al. (2009), Zhang (2010),
Negahban and Wainwright (2012), Donoho et al. (2013). These pioneering works explore the
sparsity to achieve accurate signal recovery in high dimensions.
Recently significant progresses have been made on low-rank matrix recovery under highdimensional settings. One of the most well-studied approaches is the penalized least-squares
method. Negahban and Wainwright (2011) analyzed the nuclear norm penalization in estimating nearly low-rank matrices under the trace regression model. Specifically, they derived
non-asymptotic estimation error bounds in terms of the Frobenius norm when the noise is
sub-Gaussian. Rohde and Tsybakov (2011) proposed to use a Schatten-p quasi-norm penalty
where p ≤ 1, and they derived non-asymptotic bounds on the prediction risk and Schatten-q
risk of the estimator, where q ∈ [p, 2]. Another effective method is through nuclear norm
minimization under affine fitting constraint. Other important contributions include Recht
et al. (2010), Candes and Plan (2011), Cai and Zhang (2014), Cai and Zhang (2015), etc.
When the true low-rank matrix Θ∗ satisfies certain restricted isometry property (RIP) or
similar properties, this approach can exactly recover Θ∗ under the noiseless setting and enjoy
sharp statistical error rate with sub-Gaussian and sub-exponential noise.
There has also been great amount of work on matrix completion. Candès and Recht
(2009) considered matrix completion under noiseless settings and gave conditions under
which exact recovery is possible. Candes and Plan (2010) proposed to fill in the missing
entries of the matrix by nuclear-norm minimization subject to data constraints, and showed
that rd log2 d noisy samples suffice to recover a d × d rank-r matrix with error that is proportional to the noise level. Recht (2011) improves the results of Candès and Recht (2009)
on the number of observed entries required to reconstruct an unknown low-rank matrix. Negahban and Wainwright (2012) instead used nuclear-norm penalized least squares to recover
the matrix. They derived the statistical error of the corresponding M-estimator and showed
that it matched the information-theoretic lower bound up to logarithmic factors.
Our work aims to handle the presence of heavy-tailed, asymmetrical and heteroscedastic
noises in the general trace regression. Based on the shrinkage of data, we developed a new loss
function called the robust quadratic loss, which is constructed by plugging robust covariance
b by minimizing this new
estimators in the `2 risk function. Then we obtain the estimator Θ
4

robust quadratic loss plus nuclear-norm penalty. By tailoring the analysis of Negahban et al.
(2012) to this new loss, we can establish statistical rates in estimating the matrix Θ∗ that
are the same as those in Negahban et al. (2012) for the sub-Gaussian distributions, while
allowing the noise and design to have much heavier tails. This result is very generic and
applicable to all four specific afforementioned examples.
Our robust approach is particularly simple: it truncates or shrinks appropriately the
response variables, depending on whether the responses are univariate or multivariate. Under
the setting of sub-Gaussian design, unusually large responses are very likely to be due to
the outliers of noises. This explains why we need to truncate the responses when we have
light-tailed covariates. Under the setting of heavy-tailed covariates, we need to truncate the
designs as well. It turns out that appropriate truncation does not induce significant bias or
hurt the restricted strong convexity of the loss function. With these data robustfications, we
can then apply penalized least-squares method to recover sparse vectors or low-rank matrices.
Under only bounded moment conditions for either noise or covariates, our robust estimator
achieves the same statistical error rate as that under the case of the sub-Gaussian design and
noise. The crucial component in our analysis is the sharp spectral-norm convergence rate
of robust covariance matrices based on data shrinkage. Of course, other robustifications of
estimated covariance matrices, such as the RA-covariance estimation in Fan et al. (2016), are
also possible to enjoy similar statistical error rates, but we will only focus on the shrinkage
method, as it is easier to analyze and always semi-positive definite.
It is worth emphasis that the successful application of the shrinkage sample covariance
in multi-tasking regression inspires us to also study its statistical error in covariance estimation. It turns out that as long as the random samples {xi ∈ Rd }N
i=1 have bounded fourth
moment in the sense that supv∈S d−1 E(vT xi )4 ≤ R < ∞, where S d−1 is the d-dimensional
e n achieves the statistical error rate of
unit sphere, our `4 -norm shrinkage sample covariance Σ
p
order OP ( d log d/n) in terms of the spectral norm. This rate is the same, up to a logarithmic term, as that of the standard sample covariance matrix Σn with sub-Gaussian samples
e even outperforms
under the low-dimensional regime. Under the high-dimensional regime, Σ
Σn for sub-Gaussian random samples, since now the error rate of Σn deteriorates to OP (d/n)
p
e is still OP ( d log d/n). This means even with light-tailed data,
while the error rate of Σ
standard sample covariance can be inadmissible in terms of convergence rate when dimension
is high. Therefore, shrinkage not only overcomes heavy-tailed corruption, but also mitigates
5

curse of dimensionality. In terms of the elementwise max-norm, it is not hard to show that
appropriate elementwise truncation of the data delivers a sample covariance with statistip
cal error rate of order OP ( log d/n). This estimator can further be regularized if the true
covariance has sparsity and other structure. See, for example, Meinshausen and Bühlmann
(2006), Bickel and Levina (2008), Lam and Fan (2009), Cai and Liu (2011), Cai and Zhou
(2012), Fan et al. (2013), among others.
The paper is organized as follows. In Section 2, we introduce the trace regression model
and its four well-known examples: the linear model, matrix compressed sensing, matrix completion and multi-tasking regression. Then we develop the generalized `2 loss, the truncated
and shrinkage sample covariance and corresponding M-estimators. In Section 3, we present
our main theoretical results. We first demonstrate through Theorem 1 the conditions on
the robust covariance inputs to ensure the statistical error rate of the M-estimator. Then
we apply this theorem to all the four specific aforementioned problems and derive explicitly
the statistical error rate of our M-estimators. Section 4 derives the statistical error of the
shrinkage covariance estimator in terms of the spectral norm. Finally we present simulation
studies in Section 5, which demonstrate the advantage of our robust estimator over the standard one. The associated optimization algorithms are also discussed there. All the proofs
are relegated to the Appendix A.

2

Models and methodology

We first collect the general notation before formulating the model and methodology.

2.1

Generic Notations

We follow the common convention of using boldface letters for vectors and matrices and
using regular letters for scalars. For a vector x, define kxkq to be its `q norm; specifically,
kxk1 and kxk2 denote the `1 norm and `2 norm of x respectively. We use Rd1 d2 to denote
the space of d1 d2 -dimensional real vectors, and use Rd1 ×d2 to denote the space of d1 -by-d2
real matrices. For a matrix X ∈ Rd1 ×d2 , define kXkop , kXkN , kXkF and kXkmax to be its
operator norm, nuclear norm, Frobenius norm and elementwise max norm respectively. We
use vec(X) to denote vectorized version of X, i.e., vec(X) = (XT1 , XT2 , ..., XTd2 )T , where Xj is
6

the j th column of X. Conversely, for a vector x ∈ Rd1 d2 , we use mat(x) to denote the d1 -by-d2
matrix constructed by x, where (x(j−1)d1 +1 , ..., xjd1 )T is the j th column of mat(x). For any two
matrices A, B ∈ Rd1 ×d2 , define the inner product hA, Bi := Tr(AT B) where Tr is the trace
operator. We denote diag(M1 , · · · , Mn ) to be the block diagonal matrix with the diagonal
blocks as M1 , · · · , Mn . For two Hilbert spaces A and B, we write A ⊥ B if A and B are
∞
orthogonal to each other. For two scalar series {an }∞
n=1 and {bn }n=1 , we say an  bn if there

exist constants 0 < c1 < c2 such that c1 an ≤ bn ≤ c2 an for 1 ≤ n < ∞. For a random variable
1 √
X, define its sub-Gaussian norm kXkψ2 := supp≥1 (E |X|p ) p / p and its sub-exponential norm
1

kXkψ1 := supp≥1 (E |X|p ) p /p. For a random vector x ∈ Rd , we define its sub-Gaussian norm
kxkψ2 := supv∈S d−1 kvT xkψ2 and sub-exponential norm kxkψ1 := supv∈S d−1 kvT xkψ1 . Given
x, y ∈ R, we denote max(x, y) and min(x, y) by x ∨ y and x ∧ y respectively. Let ej be the
unit vector with the j th element 1 and other elements 0.

2.2

Trace Regression

In this paper, we consider the trace regression, a general model that encompasses the linear
regression, compressed sensing, matrix completion, multi-tasking regression, etc. Suppose
N
N
we have N matrices {Xi ∈ Rd1 ×d2 }N
i=1 and responses {Yi ∈ R}i=1 . We say {(Yi , Xi )}i=1

follow the trace regression model if
Yi = hXi , Θ∗ i + i ,

(2.1)

where Θ∗ ∈ Rd1 ×d2 is the true coefficient matrix, E Xi = 0 and {i }N
i=1 are independent
noises satisfying E(i |Xi ) = 0. Note that here we do not assume {Xi }N
i=1 are independent to
each other nor assume i is independent to Xi . Model (2.1) includes the following specific
cases.
∗
• Linear regression: d1 = d2 = d, and {Xi }N
i=1 and Θ are diagonal. Let xi and

θ ∗ denote the vectors of diagonal elements of Xi and Θ∗ respectively in the context
of linear regerssion, i.e., xi = diag(Xi ) and θ ∗ = diag(Θ∗ ). Then, (2.1) reduces to
familiar linear model: Yi = xTi θ + εi . Having a low-rank Θ∗ is then equivalent to
having a sparse θ ∗ .

7

• Compressed sensing: For matrix compressed sensing, entries of Xi jointly follow the
Gaussian distribution or other ensembles. For vector compressed sensing, we can take
X and Θ∗ as diagonal matrices.
• Matrix completion: Xi is a singleton, i.e., Xi = ej(i) eTk(i) for 1 ≤ j(i) ≤ d1 and
1 ≤ k(i) ≤ d2 . In other words, a random entry of the matrix Θ is observed along with
noise for each sample.
• Multi-tasking regression: The multi-tasking (reduced-rank) regression model is
yj = Θ∗T xj + j ,

j = 1, · · · , n,

(2.2)

where xj ∈ Rd1 is the covariate vector, yj ∈ Rd2 is the response vector, Θ∗ ∈ Rd1 ×d2
is the coefficient matrix and j ∈ Rd2 is the noise with each entry independent to each
other. See, for example, Kim and Xing (2010) and Velu and Reinsel (2013). Each
sample (yj , xj ) consists of d2 responses and is equivalent to d2 data points in (2.1), i.e.,
2
. Therefore n samples in (2.2) correspond to
{(Y(j−1)d2 +i = yji , X(j−1)d2 +i = xj eTi )}di=1

N = nd2 observations in (2.1).
In this paper, we impose rank constraint on the coefficient matrix Θ∗ . Rank constraint
can be viewed as a generalized sparsity constraint for two-dimensional matrices. For linear
regression, rank constraint is equivalent to the sparsity constraint since Θ∗ is diagonal. The
rank constraint reduces the effective number of parameters in Θ∗ and arises frequently in
many applications. Consider the Netflix problem for instance, where Θ∗ij is the intrinsic score
of film j given by customer i and we would like to recover the entire Θ∗ with only partial
observations. Given that movies of similar types or qualities should receive similar scores
from viewers, columns of Θ∗ should share colinearity, thus delivering a low-rank structure
of Θ∗ . The rationale of the model can also be understood from the celebrated factor model
in finance and econometrics (Fan and Yao, 2015), which assumes that several market risk
factors drive the returns of a large panel of stocks. Consider N × T returns Y of N stocks
(like movies) over T days (like viewers). These financial returns are driven by K factors
F (K × T matrix, representing K risk factors realized on T days) with a loading matrix
B (N × K matrix), where K is much smaller than N or T . The factor model admits the

8

following form:
Y = BF + E
where E is idiosyncratic noise. Since BF has a small rank K, BF can be regarded as
the low-rank matrix Θ∗ in the matrix completion problem. If all movies were rated by all
viewers in the Netflix problem, the ratings should also be modeled as a low-rank matrix plus
noise, namely, there should be several latent factors that drive ratings of movies. The major
challenge of the matrix completion problem is that there are many missing entries.
Being exactly low-rank is still too stringent to model the real-world situations. We instead
consider Θ∗ satisfying
∗

Bq (Θ ) :=

dX
1 ∧d2

σi (Θ∗ )q ≤ ρ,

(2.3)

i=1

where 0 ≤ q ≤ 1. Note that when q = 0, the constraint (2.3) is an exact rank constraint.
Restriction on Bq (Θ∗ ) ensures that the singular values decay fast enough; it is more general
and natural than the exact low-rank assumption. In the analysis, we can allow ρ to grow
with dimensionality and sample size.
b ∈
A popular method for estimating Θ∗ is the penalized empirical loss that solves Θ
argminΘ∈S L(Θ) + λN P(Θ), where S is a convex set in Rd1 ×d2 , L(Θ) is a loss function,
λN is a tuning parameter and P(Θ) is a rank penalization function. Most of the previous
work, e.g., Koltchinskii et al. (2011) and Negahban and Wainwright (2011), chose L(Θ) =
P
∗
2
b
1≤i≤N (Yi − hΘ, Xi i) and P(Θ) = kΘkN , and derived the rate for kΘ − Θ kF under the
assumption of sub-Gaussian or sub-exponential noise. However, the `2 loss is sensitive to
outliers and is unable to handle the data with moderately heavy or heavy tails.

2.3

Robustifying `2 Loss

We aim to accomodate heavy-tailed noise and design for the nearly low-rank matrix recovery
by robustifying the traditional `2 loss. We first notice that the `2 risk can be expressed as
R(Θ) = E L(Θ) = E(Yi − hΘ, Xi i)2

= E Yi2 − 2hΘ, E Yi Xi i + vec(Θ)T E vec(Xi )vec(Xi )T vec(Θ)
≡ E Yi2 − 2hΘ, ΣY X i + vec(Θ)T ΣXX vec(Θ).

9

(2.4)

Ignoring E Yi2 , if we substitute ΣY X and ΣXX by their corresponding sample covariances,
we obtain the empirical `2 loss. This inspires us to define a generalized `2 loss via
b Y X , Θi + 1 vec(Θ)T Σ
b XX vec(Θ),
L(Θ) = −hΣ
2

(2.5)

b Y X and Σ
b XX are estimators of E Yi Xi and E vec(Xi )vec(Xi )T respectively.
where Σ
In this paper, we study the following M-estimator of Θ∗ with the generalized `2 loss:
b ∈ argmin −hΣ
b Y X , Θi + 1 vec(Θ)T Σ
b XX vec(Θ) + λN kΘkN ,
Θ
2
Θ∈S

(2.6)

where S is a convex set in Rd1 ×d2 . To handle heavy-tailed noise and design, we need to
b Y X and Σ
b XX . For ease of presentation, we always first consider
employ robust estimators Σ
the case where the design is sub-Gaussian and the response is heavy-tailed, and then further
allow the design to have heavy-tailed distribution if it is appropriate for the specific problem
setup.
We now introduce the robust covariance estimators to be used in (2.6) by the principle of
truncation, or more generally shrinkage. The intuition is that shrinkage reduces sensitivity
of the estimator to the heavy-tailed corruption. However, shrinkage induces bias. Our
theories revolve around finding appropriate shrinkage level so as to ensure the induced bias
is not too large and the final statistical error rate is sharp. Different problem setups have
b Y X and Σ
b XX , but the principle of shrinkage of data is universal. For the
different forms of Σ
linear regression, matrix compressed sensing and matrix completion, in which the response
b Y X and Σ
b XX take the following forms:
is univariate, Σ
N
X
ei
bYX = Σ
bee = 1
Yei X
Σ
YX
N i=1

N
X
b XX = Σ
bee = 1
e i )vec(X
e i )T ,
and Σ
vec(X
XX
N i=1

(2.7)

where tilde notation means truncated versions of the random variables if they have heavy tails
and equals the original random variables (truncation threshold is infinite) if they have light
tails. Note that this construction of generalized quadratic loss is equivalent to truncating
the data first and then employing the usual quadratic loss.
For the multi-tasking regression, similar idea continues to apply. However, writing (2.2)

10

in the general form of (2.1) requires adaptation of more complicated notation. We choose
d2
n
1 XX
1 b
b
ei eTj = Σ
ΣY X =
Yeij x
ey
e
x
N i=1 j=1
d2

b XX
Σ

and

d2
n
1
1 XX
b ,··· ,Σ
b xexe ) ,
=
xi eTj )T =
diag( Σ
vec(e
xi eTj )vec(e
| xexe {z
}
N i=1 j=1
d2

(2.8)

d2

where

n

X
b xeye = 1
ei y
eiT
Σ
x
n i=1

n

X
b xexe = 1
ei x
eTi
x
and Σ
n i=1

ei and x
ei are again transformed versions of yi and xi . The tilde means shrinkage for
and y
heavy-tailed variables and means identity mapping (no shrinkage) for light-tailed variables.
The factor d−1
2 is due to the fact that n independent samples under model (2.2) are treated
as nd2 samples in (2.1). As we shall see, under only bounded moment assumptions of the
design and noise, the generalized `2 loss equipped with the proposed robust covariance estib whose statistical error rates match those established
mators yields a sharp M-estimator Θ,
in Negahban and Wainwright (2011) and Negahban and Wainwright (2012) under the setting
of sub-Gaussian design and noise.

3

Main results

b defined by (2.6). In the theoretical results
Our goal is to derive the statistical error rate of Θ
of this section, we always assume d1 , d2 ≥ 2 and ρ > 1 in (2.3). We first present the following
b − Θ∗ kF .
general theorem that gives the estimation error kΘ
b =Θ
b − Θ∗ , where Θ∗ satisfies Bq (Θ∗ ) ≤ ρ. Suppose vec(∆)
b TΣ
b XX
Theorem 1. Define ∆
b ≥ κL k∆k
b 2 , where κL is a positive constant that does not depend on ∆.
b Choose
vec(∆)
F

b Y X − mat(Σ
b XX vec(Θ∗ ))kop . Then we have for some constants C1 and C2 ,
λN ≥ 2kΣ
 2−q
b 2 ≤ C1 ρ λN
k∆k
F
κL

and

b N ≤ C2 ρ
k∆k

 λ 1−q
N

κL

.

First of all, the above result is deterministic and nonasymptotic. As we can see from
b relies on the restricted eigenvalue (RE)
the theorem above, the statistical performance of Θ
11

b XX , which was first studied by Bickel et al. (2009). When the design is subproperty of Σ
b XX to be the traditional sample covariance, whose RE property has
Gaussian, we choose Σ
been well established (e.g., Rudelson and Zhou (2013), Negahban and Wainwright (2011) and
Negahban and Wainwright (2012)). We will specify these results when we need them in the
b XX = Σ
b e e,
sequel. When the design only satisfies bounded moment conditions, we choose Σ
XX

i.e., the sample covariance of shrunk X. We show that with appropriate level of shrinkage,
b e e still retains the RE property, thus satisfying the conditions of the theorem.
Σ
XX

b 2 and k∆k
b N are proportional to
Secondly, the conclusion of the theorem says that k∆k
F
b Y X − mat(Σ
b XX vec(Θ∗ )kop . This implies
λ2−q and λ1−q respectively, but we require λN ≥ kΣ
N

N

b Y X − mat(Σ
b XX vec(Θ∗ )kop is crucial to the statistical error of Θ.
b In
that the rate of kΣ
b Y X − mat(Σ
b XX vec(Θ∗ )kop for all the
the following subsections, we will derive the rate of kΣ
aforementioned specific problems with only bounded moment conditions on the response
and design. Under such weak assumptions, we show that the proposed robust M-estimator
possesses the same rates as those presented in Negahban and Wainwright (2011, 2012) with
sub-Gaussian assumptions on the design and noise.

3.1

Linear Model

For the linear regression problem, since Θ∗ and {Xi }N
i=1 are all d×d diagonal matrices, we de∗
N
note the diagonals of Θ∗ and {Xi }N
i=1 by θ and {xi }i=1 respectively for ease of presentation.

The optimization problem in (2.6) reduces to
b ∈ argmin −Σ
b xx θ + λN kθk1 ,
b T θ + 1 θT Σ
θ
Yx
2
d
θ∈R

(3.1)

bYx = Σ
b e = N −1 PN Yei x
b xx = Σ
b xexe = N −1 PN x
ei , Σ
eTi . When the design is
where Σ
i=1
i=1 ei x
e
Yx
sub-Gaussian, we only need to truncate the response. Therefore, we choose Yei = Yei (τ ) =
ei = xi , for some threshold τ . When the design is heavy-tailed, we
sgn(Yi )(|Yi | ∧ τ ) and x
choose Yei (τ ) = sgn(Yi )(|Yi | ∧ τ1 ) and x
eij = sgn(xij )(|xij | ∧ τ2 ), where τ1 and τ2 are both
predetermined threshold values. To avoid redundancy, we will not repeat stating these
choices in lemmas or theorems in this subsection.
b in (3.1), in the following lemma, we derive
To establish the statistical error rate of θ
b Y x − mat(Σ
b XX vec(Θ∗ ))kop in (2.6) for the sub-Gaussian design and boundedthe rate of kΣ

12

moment (polynomial tail) design respectively. Note here that
b Y x − mat(Σ
b XX vec(Θ∗ ))kop = kΣ
bYx − Σ
b xx θ ∗ kmax .
kΣ
Lemma 1. Uniform convergence of cross covariance.
(a) Sub-Gaussian design. Consider the following conditions:
(C1) {xi }N
i=1 are i.i.d. sub-Gaussian vectors with kxi kψ2 ≤ κ0 < ∞, E xi = 0 and
λmin (E x1 xT1 ) ≥ κL > 0;
(C2) ∀i = 1, ..., N , E |Yi |2k ≤ M < ∞ for some k > 1.
Choose τ 

p
N/ log d. For any δ > 0, there exists a constant γ1 > 0 such that as long

as log d/N < γ1 , we have
r



b Y x (τ ) − Σ
b xx θ ∗ kmax ≥ ν1
P kΣ

δ log d 
≤ 2d1−δ ,
N

(3.2)

where ν1 is a universal constant.
(b) (Bounded moment design) Consider instead the following set of conditions:
(C1’) kθ ∗ k1 ≤ R < ∞;
(C2’) E |xij1 xij2 |2 ≤ M < ∞, 1 ≤ j1 , j2 ≤ d;
(C3’) ∀i = 1, ..., N , E |Yi |4 ≤ M < ∞.
1

Choose τ1 , τ2  (N/ log d) 4 . For any δ > 0, it holds that
r



b Y x (τ1 , τ2 ) − Σ
b xx (τ2 )θ ∗ kmax > ν2
P kΣ

δ log d 
≤ 2d1−δ ,
N

where ν2 is a universal constant.
b Y x and Σ
b xx to be the sample covariance, i.e., Σ
b Y x = ΣY x =
Remark 1. If we choose Σ
PN
PN
1
1
T
b
i=1 Yi xi and Σxx = Σxx = N
i=1 xi xi , Corollary 2 of Negahban et al. (2012) showed
N
that under the sub-Gaussian noise and design,
∗

kΣY x − Σxx θ kmax
13

r log d 
= OP
.
N

This is the same rate as what we achieved under only the bounded moment conditions on
response and design.
Next we establish the restricted strong convexity of the proposed robust `2 loss.
Lemma 2. Restricted strong convexity.
(a) Sub-Gaussian design. Under Condition (C1) of Lemma 1, it holds for certain constants C1 , C2 and any η1 > 1 that


1−η1
b xx v ≥ 1 vT Σxx v − C1 η1 log d kvk2 , ∀v ∈ Rd ≥ 1 − d
P vT Σ
− 2d exp(−C2 N ).
1
2
N
3
(3.3)
(b) Bounded moment design. If xi satisfies Condition (C2’) of Lemma 1, then it holds
for some constant C3 > 0 and any η2 > 2 that


r
T

T

b xx (τ2 )v ≥ v Σxx v − C3 η2
P v Σ


log d
2
d
kvk1 , ∀v ∈ R ≤ d2−η2 ,
N

(3.4)

1

as long as τ2  (N/ log d) 4 .
Remark 2. Comparing the results we get for sub-Gaussian design and heavy-tailed design,
we can find that the coefficients before kvk21 are different. Under the sub-Gaussian design,
that coefficient is of the order log d/N , while under the heavy-tailed design, the coefficient is
p
of order log d/N . This difference will lead to different scaling requirements for N, d and ρ
in the sequel. As we shall see, the heavy-tailed design requires stronger scaling conditions to
retain the same statistical error rate as the sub-Gaussian design for the linear model.
b as defined in (3.1).
Finally we derive the statistical error rate of θ
Theorem 2. Assume

d
P

|θi∗ |q ≤ ρ, where 0 ≤ q ≤ 1.

i=1

(a) Sub-Gaussian design: Suppose Conditions (C1) and (C2) in Lemma 1 hold. For
p
p
any δ > 0, choose τ  N/ log d and λN = 2ν1 δ log d/N , where ν1 and δ are the
same as in part (a) of Lemma 1. There exist positive constants {Ci }3i=1 such that as
q

long as ρ(log d/N )1− 2 ≤ C1 , it holds that
1− 2q 


b λN ) − θ ∗ k2 > C2 ρ δ log d
P kθ(τ,
≤ 3d1−δ
2
N
14

and


 δ log d  1−q

2
∗
b
P kθ(τ, λN ) − θ k1 > C3 ρ
≤ 3d1−δ .
N
1

(b) Bounded moment design: For any δ > 0 choose τ1 , τ2  (N/ log d) 4 and λN =
p
2ν2 δ log d/N , where ν2 and δ are the same as in part (b) of Lemma 1. Under Conditions (C1’), (C2’) and (C3’), there exist constants {Ci }6i=4 such that as long as
ρ(log d/N )

1−q
2

≤ C4 , we have
 δ log d 1− 2q 
∗ 2
b
≤ 3d1−δ
P kθ(τ1 , τ2 , λN ) − θ k2 > C5 ρ
N


and



 1−q 
b 1 , τ2 , λN ) − θ ∗ k1 > C6 ρ δ log d 2 ≤ 3d1−δ .
P kθ(τ
N

b achieves the
Remark 3. Under both sub-Gaussian and heavy-tailed design, our proposed θ
minimax optimal rate of `2 norm established by Raskutti et al. (2011). However, the difference lies in the scaling requirement on N , d and ρ. For sub-Gaussian design, we require
q

ρ(log d/N )1− 2 ≤ C1 , whereas for heavy-tailed design we need ρ(log d/N )

1−q
2

≤ C4 . Under the

high-dimensional regime that d  N  log d, the former is weaker. Therefore, heavy-tailed
design requires stronger scaling than sub-Gaussian design to achieve the optimal statistical
rate.

3.2

Matrix Compressed Sensing

For the matrix compressed sensing problem, since the design is chosen by users, we consider
solely the most popular design: the Gaussian design. We thus keep the original design matrix
e i = Xi , then
and only truncate the response. In (2.7), choose Yei = sgn(Yi )(|Yi | ∧ τ ) and X
we have
N
X
bYX = Σ
b e e (τ ) = 1
Σ
sgn(Yi )(|Yi | ∧ τ )Xi
YX
N i=1

N
X
b XX = 1
and Σ
vec(Xi )vec(Xi )T . (3.5)
N i=1

b Y X − mat(Σ
b XX vec(Θ∗ ))kop . Note
The following lemma quantifies the convergence rate of kΣ
N
b Y X − mat(Σ
b XX vec(Θ∗ )) = Σ
b Y X (τ ) − 1 P hXi , Θ∗ iXi .
that here Σ
N
i=1

15

Lemma 3. Consider the following conditions:
(C1) {vec(Xi )}N
i=1 are i.i.d. sub-Gaussian vectors with kvec(Xi )kψ2 ≤ κ0 < ∞, E Xi = 0
and λmin (E vec(Xi )vec(Xi )T ) ≥ κL > 0.
(C2) ∀i = 1, ..., N , E |Yi |2k ≤ M < ∞ for some k > 1.
There exists a constant γ > 0 such that as long as (d1 + d2 )/N < γ, it holds that
r
N

X
1
d1 + d2 
∗
b
P kΣY X (τ ) −
≤ η exp(−(d1 + d2 )),
hXi , Θ iXi kop ≥ ν
N i=1
N
where τ 

(3.6)

p
N/(d1 + d2 ) and ν and η are constants.

Remark 4. For the sample covariance ΣY X =

1
N

PN

i=1

Yi Xi , Negahban and Wainwright

(2011) showed that when the noise and design are sub-Gaussian,
kΣY X −

p
1 XN
1 XN
hXi , Θ∗ iXi kop = k
i Xi kop = OP ( (d1 + d2 )/N ).
i=1
i=1
N
N

b Y X (τ ) achieves the same rate for response with just bounded moments.
Lemma 3 shows that Σ
b in (2.6).
The following theorem gives the statistical error rate of Θ
Theorem 3. Suppose Conditions (C1) and (C2) in Lemma 3 hold and Bq (Θ∗ ) ≤ ρ. We furp
p
ther assume that vec(Xi ) is Gaussian. Choose τ  N/(d1 + d2 ) and λN = 2ν (d1 + d2 )/N ,
where ν is the same as in Lemma 3. There exist constants {Ci }4i=1 such that once ρ (d1 +
1− 2q
d2 )/N
≤ C1 , we have
 d + d 1− 2q 
1
2
∗ 2
b
P kΘ(τ, λN ) − Θ kF ≥ C2 ρ
≤ η exp(−(d1 + d2 )) + 2 exp(−N/32) ,
N


and


 1−q 
b λN ) − Θ∗ kN ≥ C3 ρ d1 + d2 2 ≤ η exp(−(d1 + d2 )) + 2 exp(−N/32) ,
P kΘ(τ,
N
where η is the same constant as in Lemma 3.
Remark 5. The Frobenius norm rate here is the same rate as established under sub-Gaussian
noise in Negahban and Wainwright (2011). When q = 0, ρ is the upper bound of the rank of
16

Θ∗ and the rate of convergence depends only on ρ(d1 + d2 )  ρ(d1 ∨ d2 ), the effective number
of independent parameters in Θ∗ , rather than the ambiem number of parameters d1 ∗ d2

3.3

Matrix Completion

In this section, we consider the matrix completion problem with heavy-tailed noises. Under
√
a conventional setting, Xi is a singleton, kΘ∗ kmax = O(1) and kΘ∗ kF = O( d1 d2 ). If we
rescale the original model as
p
p
Yi = hXi , Θ∗ i + i = h d1 d2 Xi , Θ∗ / d1 d2 i + i
√
∗
d1 d2 Xi as the new design X̌i and Θ∗ / d1 d2 as the new coefficient matrix Θ̌ ,
√
∗
then kX̌i kF = O( d1 d2 ) and kΘ̌ kF = O(1). Therefore, by rescaling, we can assume

and treat

√

without loss of generality that Θ∗ satisfies kΘ∗ kF ≤ 1 and Xi is uniformly sampled from
√
{ d1 d2 · ej eTk }1≤j≤d1 ,1≤k≤d2 .
For the matrix completion problem, in order to achieve consistent estimation, we require
√
the true coefficient matrix Θ∗ not to be overly spiky, i.e., kΘ∗ kmax ≤ RkΘ∗ kF / d1 d2 ≤
√
R/ d1 d2 . We put a similar constraint in seeking the corresponding M-estimator:
b ∈
Θ

argmin
√

kΘkmax ≤R/

b Y X (τ ), Θi + 1 vec(Θ)T Σ
b XX vec(Θ) + λN kΘkN .
−hΣ
2
d1 d2

(3.7)

This spikiness condition is proposed by Negahban and Wainwright (2012) and it is required
by the matrix completion problem per se instead of our robust estimation.
To derive robust estimation in matrix completion problem, we choose Yei = sgn(Yi )(|Yi | ∧
e i = Xi in (2.7). Then, Σ
b Y X and Σ
b XX are given by (3.5). Note that the design Xi
τ ) and X
here takes the singleton form, which leads to different scaling and consistency rates from the
setting of matrix compressed sensing.
Lemma 4. Under the following conditions:
√
(C1) kΘ∗ kF ≤ 1 and kΘ∗ kmax ≤ R/ d1 d2 , where 0 < R < ∞;
√
(C2) Xi is uniformly sampled from { d1 d2 · ej eTk }1≤j≤d1 ,1≤k≤d2 ;
(C3) ∀i = 1, ..., N , E(E(2i |Xi ))k ≤ M < ∞, where k > 1;

17

there exists a constant γ > 0 such that for any δ > 0, as long as (d1 ∨ d2 ) log(d1 + d2 )/N < γ,
N
1 X
b
P kΣY X (τ ) −
hXi , Θ∗ iXi kop > ν
N i=1



r

δ(d1 ∨ d2 ) log(d1 + d2 ) 
≤ 2(d1 + d2 )1−δ , (3.8)
N

p
where τ  ( N/((d1 ∨ d2 ) log(d1 + d2 ))) and ν is a universal constant.
P
Remark 6. Again, for ΣY X = N1 N
i=1 Yi Xi , Negahban and Wainwright (2012) proved that
p
P
N
kΣY X − N1 i=1 hXi , Θ∗ iXi kop = OP ( (d1 + d2 ) log(d1 + d2 )/N ) for sub-exponential noise.
Compared with this result, Lemma 4 achieves the same rate of convergence. By Jessen’s
inequality, condition (C3) is implied by E2k
i ≤ M < ∞.
b defined in (3.7).
Now we present the following theorem on the statistical error of Θ
Theorem 4. Suppose that the conditions of Lemma 4 hold. Consider Bq (Θ∗ ) ≤ ρ with
√
kΘ∗ kmax /kΘ∗ kF ≤ R/ d1 d2 . For any δ > 0, choose
τ

p
N/((d1 ∨ d2 ) log(d1 + d2 ))) and

λN = 2ν

p
δ(d1 ∨ d2 ) log(d1 + d2 )/N

and assume (d1 ∨ d2 ) log(d1 + d2 )/N < γ, where ν and γ are the same as in Lemma 4.
There exist universal constants {Ci }4i=1 such that with probability at least 1 − 2(d1 + d2 )1−δ −
C1 exp(−C2 (d1 + d2 )) we have
n  2
1− 2q R2 o
b λN ) − Θ∗ k2 ≤ C3 max ρ δR (d1 + d2 ) log(d1 + d2 )
kΘ(τ,
,
F
N
N
and
n  δR2 (d + d ) log(d + d )  1−q
o
1
ρR2−2q  2−q
2
1
2
1
2
∗
b
,
,
kΘ(τ, λN ) − Θ kN ≤ C4 max ρ
N
N 1−q
b is defined in (3.7).
where Θ
Remark 7. Theorem 4 achieves the same statistical error rate of Frobenius norm as established in Negahban and Wainwright (2012), which also matches the information-theoretic
lower bound established in Negahban and Wainwright (2012).

18

3.4

Multi-Task Regression

Before presenting the theoretical results, we first simplify (2.6) under the setting of multi-task
regression. According to (2.8), (2.6) can be reduced to the following form:
n

X

b ∈ argmin 1 −hΣ
b xeye , Θi + 1
ei k22 + λN kΘkN .
kΘT x
Θ
d2
n i=1
Θ∈S

(3.9)

bYX −
Recall here that n is the sample size in terms of (2.2) and N = d2 n. We also have Σ

b XX vec(Θ∗ )) = Σ
b xeye − Σ
b xexe Θ∗ /d2 .
mat(Σ
Under the sub-Gaussian design, we only need to shrink the response vector yi . We choose
ei = xi and y
ei = (kyi k2 ∧ τ )yi /kyi k2 , where τ is some threshold value that depends
for (2.8) x
on n, d1 and d2 . In other words, we keep the original design, but shrink the Euclidean norm of
the response. Note that when yi is one-dimensional, the shrinkage reduces to the truncation
yi (τ ) = sgn(yi )(|yi | ∧ τ ). When the design is only of bounded moments, we need to shrink
both the design vector xi and response vector yi by their `4 norm instead, i.e., we choose
ei = (kxi k4 ∧ τ1 )xi /kxi k4 and y
ei = (kyi k4 ∧ τ2 )yi /kyi k4 , where τ1 and τ2 are two thresholds.
x
Here shrinking by the fourth order moment in fact accelerates the convergence rate of the
induced bias so that it will match the final statistical error rate. Again, we will not repeat
stating these choices in the following lemmas and theorems for less redundancy.
Lemma 5. Convergence of gradients of robustified quadratic loss.
(a) Sub-Gaussian design. Under the following conditions:
(C1) λmax (E yi yiT ) ≤ R < ∞;
(C2) {xi }ni=1 are i.i.d. sub-Gaussian vectors with kxi kψ2 ≤ κ0 < ∞, E xi = 0 and
λmin (E xi xTi ) ≥ κL > 0.
(C3) ∀i = 1, ..., n, j1 , j2 = 1, ..., d2 and j1 6= j2 , ij1 ⊥ ij2 |xi , and ∀j = 1, ..., d1 ,
k
E E(2ij |xi ) ≤ M < ∞, where k > 1;
there exists some constant γ > 0 such that if (d1 + d2 ) log(d1 + d2 )/n < γ, we have for
any δ > 0,


r
∗

b xeye (τ ) − Σ
b xexe Θ kop ≥
P kΣ

(ν1 + δ)(d1 + d2 ) log(d1 + d2 ) 
≤ 2(d1 + d2 )1−η1 δ ,
n
19

where τ 

p
n/((d1 + d2 ) log(d1 + d2 )) and ν1 and η1 are universal constants.

(b) Bounded moment design. Consider Condition (C2’) that for any v ∈ S d1 −1 ,
E(vT xi )4 ≤ M < ∞. Under Conditions (C1), (C2’) and (C3), it holds for any δ > 0

b xeye (τ1 , τ2 ) − Σ
b xexe (τ1 )Θ∗ kop ≥
P kΣ

r

(ν2 + δ)(d1 + d2 ) log(d1 + d2 ) 
≤ 2(d1 + d2 )1−η2 δ ,
n

1
where τ1 , τ2  n/((d1 + d2 ) log(d1 + d2 )) 4 and ν2 and η2 are universal constants.
Remark 8. When the noise and design are sub-Gaussian, Negahban and Wainwright (2011)
used the covering argument to show that for regular sample covariance matrices Σxy and Σxx ,
∗

kΣxy − Σxx Θ kop

n
n
p
1X
1X
∗
T
= kΣxy −
xj xj Θ kop = k
j xTj kop = OP ( (d1 + d2 )/n).
n j=1
n j=1

Lemma 5 shows that up to just a logarithmic factor, the shrinkage sample covariance achieves
nearly the same rate of convergence for noise and design with only bounded moments.
Finally we establish the statistical error rate for the low-rank multi-tasking regression.
Theorem 5. Statistical error rate for multitask regression. Assume Bq (Θ∗ ) ≤ ρ.
(a) Sub-Gaussian design. Suppose that Conditions (C1), (C2) and (C3) in Lemma 5
hold. For any δ > 0, choose
τ

p

n/((d1 + d2 ) log(d1 + d2 )) and

λN =

2p
(ν1 + δ)(d1 + d2 ) log(d1 + d2 )/n,
d2

where ν1 is the same as in Lemma 5. There exist constants γ1 , γ2 > 0 such that if
(d1 + d2 ) log(d1 + d2 )/n < γ1 and d1 + d2 ≥ γ2 , then with probability at least 1 − 3(d1 +
d2 )1−η1 δ we have

1− 2q
b λN ) − Θ∗ k2 ≤ C1 ρ (ν1 + δ)(d1 + d2 ) log(d1 + d2 )
kΘ(τ,
F
n
and

 (ν + δ)(d + d ) log(d + d )  1−q
2
1
1
2
1
2
∗
b
kΘ(τ, λN ) − Θ kN ≤ C2 ρ
,
n

where C1 and C2 are universal constants and η1 is the same as in Lemma 5.
20

(b) Bounded moment design. Suppose instead that Conditions (C1), (C2’) and (C3)
in Lemma 5 hold. For any δ > 0, choose
1
τ1 , τ2  n/((d1 +d2 ) log(d1 +d2 )) 4

and

λN =

2p
(ν2 + δ)(d1 + d2 ) log(d1 + d2 )/n,
d2

where ν2 is the same as in Lemma 5. There exist constants γ3 , γ4 > 0 such that if
(d1 + d2 ) log(d1 + d2 )/n < γ3 and d1 + d2 ≥ γ4 , then with probability at least 1 − 3(d1 +
d2 )1−η2 δ ,
 (ν + δ)(d + d ) log(d + d ) 1− 2q
2
1
2
1
2
∗ 2
b
kΘ(τ1 , τ2 , λN ) − Θ kF ≤ C3 ρ
n
and
b 1 , τ2 , λN ) − Θ∗ kN ≤ C4 ρ
kΘ(τ

 (ν + δ)(d + d ) log(d + d )  1−q
2
2
1
2
1
2
,
n

where C3 and C4 are universal constants and η2 is the same as in Lemma 5.

4

Robust Covariance Estimation

b xeye (τ1 , τ2 ) − Σ
b xexe (τ1 )Θ∗ kop in multi-tasking regression, we find that when
In derivation of kΣ
the random sample has only bounded moments, the `4 −norm shrinkage sample covariance
p
achieves the spectral norm convergence rate of order OP ( d log d/n) in estimating the true
covariance. This is nearly the optimal rate with sub-Gaussian random samples, up to a
logarithmic factor. Here we formulate the problem and the result, whose proof is relegated
to the appendix.
Suppose we have n i.i.d. d-dimensional random vectors {xi }ni=1 with E xi = 0. Our goal
is to estimate the covariance matrix Σ = E(xi xTi ) when the distribution of {xi }ni=1 has only
ei := (kxi k4 ∧ τ )xi /kxi k4 , where k · k4 is the
fourth bounded moment. For any τ ∈ R+ , let x
`4 −norm. We propose the following shrinkage sample covariance to estimate Σ.
n

X
e n (τ ) = 1
ei x
eTi .
x
Σ
n i=1

21

(4.1)

Theorem 6. Suppose E(vT xi )4 ≤ R for any v ∈ S d−1 , then it holds that for any δ > 0,


r

e n (τ ) − Σkop ≥
P kΣ

δRd log d 
≤ d1−Cδ ,
n

(4.2)

1/4
where τ  nR/(δ log d)
and C is a universal constant.
We have several comments for the above theorem. First of all, unlike the bounded
moment conditions in the previous section, R here can go to infinity with certain rates.
Hence we also put this quantity into the rate of convergence. If indeed R < ∞, we recover
p
the true covariance matrix with statistical error of OP ( d log d/n). In comparison with the
robust covariance matrix estimator given in Fan et al. (2016), our estimator here is positive
semidefinite and is very simple to implement. Our concentration inequality here is on the
operator norm, while their result is on the element-wise max norm.
If we are concerned with statistical error in terms of elementwise max norm, then we
need to apply elementwise truncation to the random samples rather than `4 norm shrinkage.
n
P
Let x̌i satisfy x̌ij = sgn(xij )(|xij | ∧ τ ) for 1 ≤ j ≤ d1 and Σ̌n =
x̌i x̌Ti . It is not hard to
i=1
p
1
derive kΣ̌n − Σkmax = OP ( log d/n) as in Fan et al. (2016) with τ  (n/ log d) 4 . Further
regularization can be applied to Σ̌n if the true covariance is sparse. See, for example,
Meinshausen and Bühlmann (2006), Bickel and Levina (2008), Lam and Fan (2009), Cai
and Liu (2011), Cai and Zhou (2012), Fan et al. (2013), among others.
Theorem 6 is non-asymptotic as all the results in the previous sections and can be applied
e n (τ ) can outperform the
to both low-dimensional and high-dimensional regimes. In fact, Σ
sample covariance Σn even with Gaussian samples if the dimension is high. The reason is
p
that according to Theorem 5.39 in Vershynin (2010), kΣn − Σkop = OP ( d/n ∨ (d/n)).
p
When d/n is large, the d/n term will dominate d/n, thus delivering statistical error of
order d/n for Gaussian sample covariance. However, our shrinkage sample covariance always
p
retains the statistical error of order d log d/n regardless of ratio between the dimension and
sample size. Therefore, shrinkage overcomes not only heavy-tailed corruption, but also curse
of dimensionality. In Section 5.4, we conduct simulations to further illustrate this point.

22

5

Simulation Study

In this section, we first compare the numerical performance of the robust procedure and
standard procedure in compressed sensing, matrix completion and multi-tasking regression.
For each setup, we investigate three noise settings: log-normal noise, truncated Cauchy
noise and Gaussian noise. They represent heavy-tailed asymmetric distributions, heavytailed symmetric distributions and light-tailed symmetric distributions. The results from
the standard and robust methods are shown in the same color for each scenario in the
following figures so that they can be compared more easily.
All the robust procedures proposed in our work are very easy to implement; we only need
to truncate or shrink the data appropriately, and then apply the standard procedure to the
transformed data. As for the parameter tuning, we refer to the rate developed in our theories
for robust procedures and theories in Negahban and Wainwright (2011) and Negahban and
Wainwright (2012) for standard procedures. The constants before the rate are tuned for
best performance. The main message is that the robust procedure outperforms the standard
procedure under the setting with bounded moment noise, and it performs equally well as the
standard procedure under the Gaussian noise. The simulations are based on 100 independent
Monte Carlo replications. There are no standard algorithms for solving penalized trace
regression problems. Hence, besides presenting the numerical performance, we also elucidate
the algorithms that we use to solve the corresponding optimization problems, which might be
of interest to readers concerned with implementations. These algorithms are not necessarily
the same as those used in the literature.
We also compare the numerical performance of the regular sample covariance and shrinkage sample covariance as proposed in (4.1) in estimating the true covariance. We choose
d/n = 0.2, 0.5, 1 and for each ratio, we let n = 100, 200, ..., 500. Simulation results show superiority of the shrinkage sample covariance over the regular sample covariance under both
Gaussian noise and t3 noise. Therefore, the shrinkage can not only overcome the heavy-tailed
corruption, but also mitigate the curse of high dimensions.

5.1

Compressed Sensing

We first specify the parameters in the true model: Y = hXi , Θ∗ i + i . We always set
√
d1 = d2 = d, kΘ∗ kF = 5 and rank(Θ∗ ) = 5. In the simulation, we construct Θ∗ to be
23

5
P

vi viT , where vi is the ith top eigevector of the sample covariance of 100 i.i.d. centered

i=1

Gaussian random vectors with covariance Id . The design matrix Xi has i.i.d. standard
Gaussian entries. The noise distributions are characterized as follows:
d

• Log-normal: i = (Z − E Z)/50, where Z ∼ lnN (0, σ 2 ) and σ 2 = 6.25;
d

• Truncated Cauchy: i = min(Z, 103 )/10, where Z follows Cauchy distribution;
• Gaussian: i ∼ N (0, σ 2 ), where σ 2 = 0.25.
The constants above are chosen to ensure appropriate signal-to-noise ratio for better presentation. We present the numerical results in Figure 2. As we can observe from the plots,
the robust estimator has much smaller statistical error than the standard estimator under
the heavy-tailed noise, i.e., the log-normal and truncated Cauchy noise. When d = 40 or 60,
robust procedures deliver sharper estimation as the sample size increases, while the standard procedure does not necessarily do so under the heavy-tailed noise. Under the setting
of Gaussian noise, the robust estimator has nearly the same statistical performance as the
standard one, which shows that it does not hurt to use the robust procedure under the
light-tailed noise setting.
As for the implementation, we exploit the contractive Peaceman-Rachford splitting method
(PRSM) to solve the compressed sensing problem. Here we briefly introduce the general
scheme of the contractive PRSM for clarity. The contractive PRSM is for minimizing the
summation of two convex functions under linear constraint:
min

x∈Rp1 ,y∈Rp2

f1 (x) + f2 (y),
(5.1)

subject to C1 x + C2 y − c = 0,
where C1 ∈ Rp3 ×p1 , C2 ∈ Rp3 ×p2 and c ∈ Rp3 . The general iteration scheme of the contractive

24

1.0
0.5

1.0

1.5

6.0

-1.0

6.5

7.0

7.5

6.5

Log N

Log-normal Noise

0.0
-0.5
-1.0

Log Frobenius Error

6.0

-1.5

d=20, Standard
d=20, Robust
d=40, Standard
d=40, Robust
d=60, Standard
d=60, Robust

-2.0

-0.5

0.0

Log Frobenius Error

0.5

1.0
0.5

Log Frobenius Error

0.0
-0.5
-1.0

d=20, Standard
d=20, Robust
d=40, Standard
d=40, Robust
d=60, Standard
d=60, Robust

7.0

7.5

Log N

Truncated Cauchy Noise

d=20, Standard
d=20, Robust
d=40, Standard
d=40, Robust
d=60, Standard
d=60, Robust
6.0

6.5

7.0

Log N

Gaussian Noise

b − Θ∗ kF v.s. logarithmic sample size ln N for different
Figure 2: Statistical errors of lnkΘ
dimensions d in matrix compressed sensing.
PRSM is


β

(k+1)

x
=
argmin
f1 (x) − (ρ(k) )T (C1 x + C2 y(k) − c) + kC1 x + C2 y(k) − ck22 ,



2
x



1

 ρ(k+ 2 ) = ρ(k) − αβ(C1 x(k+1) + C2 y(k) − c),

1
β
 (k+1)

y
= argmin f2 (y) − (ρ(k+ 2 ) )T (C1 x(k+1) + C2 y − c) + kC1 x(k+1) + C2 y − ck22 ,



2
y




 ρ(k+1) = ρ(k+ 21 ) − αβ(C x(k+1) + C y(k+1) − c),
1
2
(5.2)
where ρ ∈ Rp3 is the Lagrangian multiplier, β is the penalty parameter and α is the relaxation
factor (Eckstein and Bertsekas (1992)). Since the parameter of interest in our work is Θ∗ , now
we use θ x ∈ Rd1 d2 and θ y ∈ Rd1 d2 to replace x and y respectively in (5.2). By substituting
N
1 X
f1 (θ x ) =
(Yi − hvec(Xi ), θ x i)2 , f2 (θ y ) = λkmat(θ y )kN , C1 = I, C2 = −I and c = 0
N i=1

into (5.2), we can obtain the PRSM algorithm for the compressed sensing problem. Let X
be a N -by-d1 d2 matrix whose rows are i.i.d. random designs {vec(Xi )}N
i=1 and Y be the
25

7.5

N -dimensional response vector. Then we have the following iteration scheme specifically for
compressed sensing.

(k)

+ 2XT Y/N ),
= (2XT X/N + β · I)−1 (β · θ (k)
θ (k+1)

y +ρ
x




 ρ(k+ 21 ) = ρ(k) − αβ(θ (k+1) − θ (k) ),
y

x

1


= vec(Sλ/β (mat(θ x − ρ(k+ 2 ) ))),
θ (k+1)
y




 ρ(k+1) = ρ(k+ 12 ) − αβ(θ (k+1) − θ (k+1) ),
y
x

(5.3)

where we choose α = 0.9 and β = 1 according to Eckstein and Bertsekas (1992) and He
et al. (2014), ρ ∈ Rd1 d2 is the Lagrangian multiplier and Sτ (z) is the singular value softthresholding function for matrix version of z ∈ Rd1 d2 . To be more specific, let Z = mat(z) ∈
Rd1 ×d2 and Z = UΛVT = U diag(λ1 , ..., λr ) VT be its singular value decomposition. Then

Sτ (z) = vec U diag((λ1 − τ )+ , (λ2 − τ )+ , ..., (λr − τ )+ ) VT , where (x)+ = max(x, 0). The
algorithm stops if kθ x − θ y k2 is smaller than some predetermined threshold, and returns
mat(θ y ) as the final estimator of Θ∗ .

5.2

Matrix Completion

We again set d1 = d2 = d and construct Θ∗ to be

5
P

√
vi viT / 5, where vi is the ith top

i=1

eigevector of the sample covariance of 100 i.i.d. centered Gaussian random vectors with
covariance Id . Each design matrix Xi takes the singleton form, which is uniformly sampled
from {ej eTk }1≤j,k≤d . The noise distributions are
d

• Log-normal: i = (Z − E Z)/250, where Z ∼ ln N (0, σ 2 ) and σ 2 = 9;
d

• Truncated Cauchy: i = min(Z, 103 )/16, where Z follows Cauchy distribution;
• Gaussian: i ∼ N (0, σ 2 ), where σ 2 = 0.25.
Again, the constants above are set for an appropriate signal-to-noise ratio for better presentation. We present the numerical results in Figure 3. Analogous to the matrix compressed
sensing, we can observe from the figure that compared with the standard procedure, the
robust procedure has significantly smaller statistical error in estimating Θ∗ under the lognormal and truncated Cauchy noise. Under Gaussian noise, the robust procedure has nearly
the same statistical performance as the standard procedure.
26

7.5

8.0

9.0

9.5

0.0
-1.5

d=50, Standard
d=50, Robust
d=100, Standard
d=100, Robust
d=150, Standard
d=150, Robust
7.5

-2.0

-2.0

8.5

8.0

8.5

Log N

Log-normal Noise

-1.0

Log Frobenius Error

-0.5

0.0
-1.5

-1.0

Log Frobenius Error

-0.5

0.0
-0.5
Log Frobenius Error

-1.0
-1.5
-2.0

d=50, Standard
d=50, Robust
d=100, Standard
d=100, Robust
d=150, Standard
d=150, Robust

9.0

9.5

d=50, Standard
d=50, Robust
d=100, Standard
d=100, Robust
d=150, Standard
d=150, Robust
7.5

8.0

Log N

Truncated Cauchy Noise

8.5

9.0

Log N

Gaussian Noise

b − Θ∗ kF v.s. logarithmic sample size ln N for different
Figure 3: Statistical errors of lnkΘ
dimensions d in matrix completion.
To solve the matrix completion problem in (3.7), we adapt the ADMM method inspired
by Fang et al. (2015). They propose to recover the matrix by minimizing the square loss
plus both nuclear norm and matrix max-norm penalizations under the entrywise max-norm
constraint. By simply setting the penalization coefficient for the matrix max-norm to be
zero, we can apply their algorithm to our problem. Let L, R, W ∈ R(d1 +d2 )×(d1 +d2 ) , which
N
P
are variables in our algorithm. Define Θn , Θs ∈ Rd1 ×d2 such that Θnij =
1{Xt =ei eTj } and
t=1

Θsij =

N
P
t=1

Yt 1{Xt =ei eTj } , so Θn and Θs are constants given the data. Below we present the

iteration scheme of our algorithm to solve the matrix completion problem (3.7). Readers
who are interested in the derivation of the algorithm can refer to Fang et al. (2015) for the
technical details.

27

9.5



L(k+1) = ΠS d1 +d2 {R(k) − ρ−1 (W(k) + λN I)},


+


!


11
12

C
C


C=
= L(k+1) + W(k) /ρ,


21
22

C
C




s
n
12
R12
ij = Π[−R,R] {(ρCij + 2Θij /N )/(ρ + 2Θij /N )}, 1 ≤ i ≤ d1 , 1 ≤ j ≤ d2


!


11
12

C
R



R(k+1) =
,

12 T
22


(R
)
C




 W(k+1) = W(k) + γρ(L(k+1) − R(k+1) ).

(5.4)

In the algorithm above, ΠS d1 +d2 (·) is the projection operator onto the space of positive
+

semidefinite matrices S+d1 +d2 , ρ is the penalization parameter which we set to be 0.1 in our
simulation and γ is the step length which is typically set to be 1.618 according to Fang et al.
(2015). We omit the stopping criteria for this algorithm here since it is complex. Readers
who are interested in the implementation details of this algorithm can refer to Fang et al.
(2015) for detailed instruction. Once the stopping criteria is satisfied, the algorithm returns
R12 as the final estimator of Θ∗ .

5.3

Multi-Tasking Regression

We again set d1 = d2 = d and construct Θ∗ to be

5
P

vi viT , where vi is the ith top eigevector

i=1

of the sample covariance of 100 i.i.d centered Gaussian random vectors with covariance Id .
The design vectors {xi }ni=1 are i.i.d. Gaussian vectors with covariance matrix Id . The noise
distributions are characterized as follows:
d

• Log-normal: i = (Z − E Z)/50, where Z ∼ lnN (0, σ 2 ) and σ 2 = 4;
d

• Truncated Cauchy: i = min(Z, 104 )/10, where Z follows Cauchy distribution;
• Gaussian: i ∼ N (0, σ 2 ), where σ 2 = 0.25.
We present the numerical results in Figure 4. Similar to the two examples before, the robust
procedure has sharper statistical accuracy in estimating Θ∗ than the standard procedure
under both heavy-tailed noises. Under Gaussian noise, the robust procedure has nearly the
same statistical performance as the standard procedure.
28

1.0

0.5

1.0

6.5

7.5

8.0

6.5

-0.5
-1.5

d=50, Standard
d=50, Robust
d=100, Standard
d=100, Robust
d=150, Standard
d=150, Robust

-2.0

-1.5
-2.0

7.0

-1.0

0.0
-1.0

-0.5

Log Frobenius Error

Log Frobenius Error

0.5

0.0

0.5
0.0
-0.5

Log Frobenius Error

-1.0
-1.5
-2.0

d=50, Standard
d=50, Robust
d=100, Standard
d=100, Robust
d=150, Standard
d=150, Robust

7.0

Log N

7.5

8.0

d=50, Standard
d=50, Robust
d=100, Standard
d=100, Robust
d=150, Standard
d=150, Robust
6.5

7.0

Truncated Cauchy Noise

8.0

Log N

Log N

Log-normal Noise

7.5

Gaussian Noise

b − Θ∗ kF v.s. logarithmic sample size ln N for different
Figure 4: Statistical errors of lnkΘ
dimensions d in multi-tasking regression.
For multi-tasking regression, we exploit the contractive PRSM method again. Let X be
the n-by-d1 design matrix and Y be the n-by-d2 response matrix. By following the general
iteration scheme (5.2), we can develop the iteration steps for the multi-tasking regression.

(k)

Θ(k+1)
= (2XT X/n + β · I)−1 (β · Θ(k)
+ 2XT Y/n),

x
y +ρ




 ρ(k+ 21 ) = ρ(k) − αβ(Θ(k+1) − Θ(k) ),
x
y
1


Θ(k+1)
= Sλ/β (Θx − ρ(k+ 2 ) ),
y




 ρ(k+1) = ρ(k+ 12 ) − αβ(Θ(k+1) − Θ(k+1) ),
x
y

(5.5)

where Θx , Θy , ρ ∈ Rd1 ×d2 and Sτ (·) is the same singular value soft thresholding function
as in (5.3). Analogous to the compressed sensing, we choose α = 0.9 and β = 1. As long
as kΘx − Θy kF is smaller than some predetermined threshold, the iteration stops and the
algorithm returns Θy as the final estimator of Θ∗ .

29

5.4

Covariance Estimation

In this subsection, we investigate the statistical error of the sample covariance Σn and
e n (τ ) proposed in Section 4 when the random samples are
shrinkage sample covariance Σ
heavy-tailed. We only consider two simple distributions: Gaussian and Student’s t3 random
samples. The dimension is set to be proportional to sample size, i.e., d/n = α with α being
0.2, 0.5, 1. n will range from 100 to 500 for each case. Regardless of how large the dimension
d is, the true covariance Σ is always set to be a diagonal matrix with the first diagonal
element equal to 4 and all the other diagonal elements equal to 1. We present our results
in Figure 5. The statistical errors are measured in terms of the spectral norm gap between
the estimator and true covariance, and our simulation is based on 1, 000 independent Monte

d=0.2n, Standard
d=0.2n, Robust
d=0.5n, Standard
d=0.5n, Robust
d=n, Standard
d=n, Robust

8
4

6

Statistical Error

2.0
1.5
1.0

Statistical Error

10

2.5

12

3.0

14

Carlo replications.

100

200

300

400

2
0

0.0

0.5

d=0.2n, Standard
d=0.2n, Robust
d=0.5n, Standard
d=0.5n, Robust
d=n, Standard
d=n, Robust
500

100

n

200

300

400

500

n

Gaussian Samples

t3 Samples

b − Σkop v.s. sample size n for different dimensions d in
Figure 5: Statistical errors of kΣ
covariance estimation.
As we can see, for Gaussian samples, as long as we fix d/n, the statistical error of both
e n (τ ) does not change, which is consistent with Theorem 5.39 in Vershynin (2010)
Σn and Σ
and Theorem 6 in our paper. Also, the higher the dimension is, the more significant the
e n (τ ) is over Σn . This validates our remark after Theorem 6 that the shrinkage
superiority of Σ
ameliorates the impact of dimensionality. Even for Gaussian data, shrinkage is meaningful
30

and provides significant improvement. For t3 distribution, since it is heavy-tailed, the regular
sample covariance does not maintain constant statistical error for a fixed d/n; instead the
error increases as the sample size increases. In contrast, our shrinkage sample covariance
still retains stable statistical error and enjoys much higher accuracy than the regular sample
e n (τ ) in
covariance. This strongly supports the sharp statistical error rate we derived for Σ
Theorem 6.

References
Bickel, P. J. and Levina, E. (2008). Covariance regularization by thresholding. Annals
of Statistics 36 2577–2604.
Bickel, P. J., Ritov, Y. and Tsybakov, A. B. (2009). Simultaneous analysis of lasso
and dantzig selector. Annals of Statistics 37 1705–1732.
Boucheron, S., Lugosi, G. and Massart, P. (2013). Concentration inequalities: A
nonasymptotic theory of independence. OUP Oxford.
Brownlees, C., Joly, E. and Lugosi, G. (2015). Empirical risk minimization for heavytailed losses. Annals of Statistics 43 2507–2536.
Cai, T. and Liu, W. (2011). Adaptive thresholding for sparse covariance matrix estimation.
Journal of the American Statistical Association 106 672–684.
Cai, T. T. and Zhang, A. (2014). Sparse representation of a polytope and recovery
of sparse signals and low-rank matrices. IEEE Transactions on Information Theory 60
122–132.
Cai, T. T. and Zhang, A. (2015). Rop: Matrix recovery via rank-one projections. Annals
of Statistics 43 102–138.
Cai, T. T. and Zhou, H. H. (2012). Optimal rates of convergence for sparse covariance
matrix estimation. Annals of Statistics 40 2389–2420.
Candes, E. and Tao, T. (2007). The dantzig selector: Statistical estimation when p is
much larger than n. Annals of Statistics 35 2313–2351.
31

Candes, E. J. (2008). The restricted isometry property and its implications for compressed
sensing. Comptes Rendus Mathematique 346 589–592.
Candes, E. J. and Plan, Y. (2010). Matrix completion with noise. Proceedings of the
IEEE 98 925–936.
Candes, E. J. and Plan, Y. (2011). Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements. IEEE Transactions on
Information Theory 57 2342–2359.
Candès, E. J. and Recht, B. (2009). Exact matrix completion via convex optimization.
Foundations of Computational mathematics 9 717–772.
Candes, E. J. and Tao, T. (2006). Near-optimal signal recovery from random projections:
Universal encoding strategies? IEEE Transactions on Information Theory 52 5406–5425.
Catoni, O. (2012). Challenging the empirical mean and empirical variance: a deviation
study. Annales de l’Institut Henri Poincaré 48 1148–1185.
Chen, S. S., Donoho, D. L. and Saunders, M. A. (2001). Atomic decomposition by
basis pursuit. SIAM review 43 129–159.
Donoho, D. L. (2006). Compressed sensing. IEEE Transactions on Information Theory
52 1289–1306.
Donoho, D. L., Johnstone, I. and Montanari, A. (2013). Accurate prediction of
phase transitions in compressed sensing via a connection to minimax denoising. IEEE
Transactions on Information Theory 59 3396–3433.
Eckstein, J. and Bertsekas, D. P. (1992). On the douglas—rachford splitting method
and the proximal point algorithm for maximal monotone operators. Mathematical Programming 55 293–318.
Fan, J., Li, Q. and Wang, Y. (2016). Robust estimation of high-dimensional mean
regression. Journal of Royal Statistical Society, Series B .
Fan, J. and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its
oracle properties. Journal of the American Statistical Association 96 1348–1360.
32

Fan, J., Liao, Y. and Mincheva, M. (2013). Large covariance estimation by thresholding
principal orthogonal complements. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 75 603–680.
Fan, J. and Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature
space. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 70
849–911.
Fan, J. and Yao, Q. (2015). Elements of Financial Econometrics. Science Press, Beijing.
Fang, X. E., Liu, H., Toh, K. C. and Zhou, W.-X. (2015). Max-norm optimization for
robust matrix recovery. Tech. rep.
Gupta, S., Ellis, S. E., Ashar, F. N., Moes, A., Bader, J. S., Zhan, J., West,
A. B. and Arking, D. E. (2014). Transcriptome analysis reveals dysregulation of innate immune response genes and neuronal activity-dependent genes in autism. Nature
communications 5.
He, B., Liu, H., Wang, Z. and Yuan, X. (2014). A strictly contractive peaceman–
rachford splitting method for convex programming. SIAM Journal on Optimization 24
1011–1040.
Hsu, D. and Sabato, S. (2016). Loss minimization and parameter estimation with heavy
tails. Journal of Machine Learning Research 17 1–40.
Ibragimov, R. and Sharakhmetov, S. (1998). Short communications: On an exact
constant for the rosenthal inequality. Theory of Probability & Its Applications 42 294–
302.
Kim, S. and Xing, E. P. (2010). Tree-guided group lasso for multi-task regression with
structured sparsity. Annals of Applied Statistics 6 1095–1117.
Koltchinskii, V., Lounici, K. and Tsybakov, A. B. (2011). Nuclear-norm penalization
and optimal rates for noisy low-rank matrix completion. Annals of Statistics 39 2302–2329.
Lam, C. and Fan, J. (2009). Sparsistency and rates of convergence in large covariance
matrix estimation. Annals of Statistics 37 4254–4278.
33

Loh, P.-L. (2015). Statistical consistency and asymptotic normality for high-dimensional
robust m-estimators. arXiv preprint arXiv:1501.00312 .
Meinshausen, N. and Bühlmann, P. (2006). High-dimensional graphs and variable selection with the lasso. Annals of Statistics 34 1436–1462.
Minsker, S. (2015). Geometric median and robust estimation in banach spaces. Bernoulli
21 2308–2335.
Negahban, S. and Wainwright, M. J. (2011). Estimation of (near) low-rank matrices
with noise and high-dimensional scaling. Annals of Statistics 39 1069–1097.
Negahban, S. and Wainwright, M. J. (2012). Restricted strong convexity and weighted
matrix completion: Optimal bounds with noise. Journal of Machine Learning Research
13 1665–1697.
Negahban, S., Yu, B., Wainwright, M. J. and Ravikumar, P. K. (2012). A unified
framework for high-dimensional analysis of m-estimators with decomposable regularizers.
Statistical Science 24 538–577.
Nemirovsky, A.-S., Yudin, D.-B. and Dawson, E.-R. (1982). Problem complexity and
method efficiency in optimization .
Nowak, R. D., Wright, S. J. et al. (2007). Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems. IEEE Journal of
Selected Topics in Signal Processing 1 586–597.
Oliveira, R. I. (2016). The lower tail of random quadratic forms with applications to
ordinary least squares. Probability Theory and Related Fields 1–20.
Raskutti, G., Wainwright, M. J. and Yu, B. (2011). Minimax rates of estimation for
high-dimensional linear regression over-balls. IEEE Transactions on Information Theory
57 6976–6994.
Recht, B. (2011). A simpler approach to matrix completion. Journal of Machine Learning
Research 12 3413–3430.

34

Recht, B., Fazel, M. and Parrilo, P. A. (2010). Guaranteed minimum-rank solutions
of linear matrix equations via nuclear norm minimization. SIAM review 52 471–501.
Rohde, A. and Tsybakov, A. B. (2011). Estimation of high-dimensional low-rank matrices. Annals of Statistics 39 887–930.
Rudelson, M. and Zhou, S. (2013). Reconstruction from anisotropic random measurements. IEEE Transactions on Information Theory 59 3434–3447.
Stock, J. H. and Watson, M. W. (2002). Macroeconomic forecasting using diffusion
indexes. Journal of Business & Economic Statistics 20 147–162.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the
Royal Statistical Society: Series B (Methodological) 58 267–288.
Tropp, J. A. (2012). User-friendly tools for random matrices: An introduction. Tech. rep.,
DTIC Document.
Tropp, J. A. (2015). An introduction to matrix concentration inequalities. arXiv preprint
arXiv:1501.01571 .
Velu, R. and Reinsel, G. C. (2013). Multivariate reduced-rank regression: theory and
applications, vol. 136. Springer Science & Business Media.
Vershynin, R. (2010). Introduction to the non-asymptotic analysis of random matrices.
arXiv preprint arXiv:1011.3027 .
Zhang, C.-H. (2010). Nearly unbiased variable selection under minimax concave penalty.
Annals of statistics 38 894–942.
Zou, H. and Li, R. (2008). One-step sparse estimates in nonconcave penalized likelihood
models. Annals of Statistics 36 1509–1533.

35

A

Proofs

Proof of Theorem 1. This Lemma is just a simple application of the theoretical framework
established in Negahban et al. (2012), but for completeness and clarity, we present the whole
b
proof here. We start from the optimality of Θ:
1
∗
∗
∗
∗ T b
b Y X , Θi+
b 1 vec(Θ)
b XX vec(Θ)+λ
b
b
b
b TΣ
−hΣ
N kΘkN ≤ −hΣY X , Θ i+ vec(Θ ) ΣXX vec(Θ )+λN kΘ kN .
2
2
b =Θ
b − Θ∗ . Simple algebra delivers that
Note that ∆
1
b XX vec(∆)
b ≤ hΣ
b Y X − mat(Σ
b XX vec(Θ∗ )), ∆i
b + λN k∆k
b N
b TΣ
vec(∆)
2
b Y X − mat(Σ
b XX vec(Θ∗ ))kop · k∆k
b N + λN k∆k
b N ≤ 2λN k∆k
b N,
≤ kΣ
(A.1)
b Y X − mat(Σ
b XX vec(Θ∗ ))kop . To bound the RHS of (A.1), we need to decompose
if λN ≥ 2kΣ
b as Negahban and Wainwright (2011) did. Let Θ∗ = UDVT be the SVD of Θ∗ , where the
∆
diagonals of D are in the decreasing order. Denote the first r columns of U and V by Ur
and Vr respectively, and define
M := {Θ ∈ Rd1 ×d2 | row(Θ) ⊆ col(Vr ), col(Θ) ⊆ col(Ur )},

(A.2)

⊥

M := {Θ ∈ Rd1 ×d2 | row(Θ) ⊥ col(Vr ), col(Θ) ⊥ col(Ur )},

where col(·) and row(·) denote the column space and row space respectively. For any ∆ ∈
Rd1 ×d2 and Hilbert space W ⊆ Rd1 ×d2 , let ∆W be the projection of ∆ onto W. We first
clarify here what ∆M , ∆M and ∆M⊥ are. Write ∆ as
"
⊥

∆ = [Ur , Ur ]

Γ11 Γ12
Γ21 Γ22

#
⊥

[Vr , Vr ]T ,

then the following equalities hold:
"
∆M = Ur Γ11 (Vr )T ,

r⊥

r⊥

∆M⊥ = U Γ22 (V )T ,

r⊥

∆M = [Ur , U ]

Γ11 Γ12
Γ21

0

#
⊥

[Vr , Vr ]T .
(A.3)

36

Applying Lemma 1 in Negahban et al. (2012) to our new loss function implies that if λN ≥
b Y X − mat(Σ
b XX vec(Θ∗ ))kop , it holds that
2kΣ
b ⊥ kN ≤ 3k∆
b kN + 4
k∆
M
M

X
j≥r+1

σj (Θ∗ ).

(A.4)

b ) ≤ 2r; we thus have
Note that rank(∆
M
b kN + 4
b N ≤ k∆
b kN + k∆
b ⊥ kN ≤ 4k∆
k∆k
M
M
M

X

X
√
b F +4
σj (Θ∗ ).
σj (Θ∗ ) ≤ 4 2rk∆k
j≥r+1

j≥r+1

(A.5)
Following the proof of Corollary 2 in Negahban and Wainwright (2011), we determine the
value of r here. For a threshold τ > 0, we choose
r = #{j ∈ {1, 2, ..., (d1 ∧ d2 )}|σj (Θ∗ ) ≥ τ }.
Then it follows that
X

σj (Θ∗ ) ≤ τ

j≥r+1

X σj (Θ∗ )
X σj (Θ∗ ) q
X
≤τ
≤ τ 1−q
σj (Θ∗ )q ≤ τ 1−q ρ.
τ
τ
j≥r+1
j≥r+1
j≥r+1

On the other hand, ρ ≥

P

(A.6)

σj (Θ∗ )q ≥ rτ q , so r ≤ ρτ −q . Combining (A.1), (A.5) and

j≤r

b TΣ
b XX vec(∆)
b ≥ κL k∆k
b 2 , we have
vec(∆)
F
√
1
b 2 ≤ 2λN (4 2rk∆k
b F + 4τ 1−q ρ),
κL k∆k
F
2
which implies that
r
r
λN ρ  32λN τ −q √ 1−q 
b
k∆kF ≤ 4
+ τ
.
κL
κL
Choosing τ = λN /κL , we have for some constant C1 ,
q

 1− 2
b F ≤ C 1 √ ρ λN
k∆k
.
κL

37

Combining this result with (A.5) and (A.6), we can further derive the statistical error rate
in terms of the nuclear norm as follows.
 1−q
b N ≤ C2 ρ λN
k∆k
,
κL
where C2 is certain positive constant.
Proof of Lemma 1. (a) We first prove for the case of the sub-Gaussian design. Recall that
P e
ei = xi in this case. Let σ
we use Yei to denote sgn(Yi )(|Yi | ∧ τ ) and x
bxj Ye (τ ) = N1 N
i=1 Yi xij .
Note that
2k

k−1
1
1
Var(Yei xij ) ≤ E(Yei2 x2ij ) ≤ E(Yi2 x2ij ) ≤ (E Yi2k ) k (E xijk−1 ) k ≤ 2M k κ20 k/(k − 1) < ∞ ,
(A.7)

which is a constant that we denote by v1 . In addition, for p > 2,
E |Yei xij |p ≤ τ p−2 E(Yi2 |xij |p ) ≤ τ p−2 M

1
k

E |xij |

kp
k−1

 k−1
k

 r kp p
.
≤ τ p−2 M κ0
k−1
1
k

By the Jensen’s inequality and then Stirling approximation, it follows that for some constants
c1 and c2 ,
E |Yei xij − E Yei xij |p ≤ 2p−1 (E |Yei xij |p + | E Yei xij |p ) ≤ 2p−1 (E |Yi xij |p + (E |Yi xij |)p )
≤

c1 p!(c2 τ )p−2
.
2

Define v := c1 ∨ v1 . According to Bernstein’s Inequality (Theorem 2.10 in Boucheron et al.
(2013)), we have for j = 1, ..., d,


r

P |b
σxj Ye (τ ) − E(Yei xij )| ≥

2vt c2 τ t 
+
≤ 2 exp(−t).
N
N

Also note that by Markov’s inequality,
r
q
v E Yi2
E((Yei − Yi )xij ) ≤ E(|Yi xij | · 1|Yi |>τ ) ≤ E(Yi2 x2ij )P (|Yi | > τ ) ≤
≤
τ2
1

p

1

vM k
.
τ
1

Note that since Var(Yi ) = β ∗T Σxx β ∗ + E 2 ≤ M k and λmin (Σxx ) ≥ κL , kβ ∗ k2 ≤ M k /κL .
38

1

1

Therefore, kxTi β ∗ kψ2 ≤ M k κ0 /κL and kxij xTi β ∗ kψ1 ≤ 2M k κ20 /κL . By Proposition 5.16
(Bernstein-type inequality) in Vershynin (2010), we have for sufficiently small t,
r



∗

T

b x x ) β − E(Yi xij )| ≥ c2
P |(Σ
j
for some constant c2 . Choose τ 

p

t
≤ exp(−t)
N

N/ log d. An application of the triangle inequality and

the union bound yields that as long as log d/N < γ1 for certain γ1 , we have for some constant
ν1 > 0,
N
1 X T ∗
b
P kΣY x (τ ) −
(x θ )xi kmax ≥ ν1
N i=1 i



r

δ log d 
≤ 2d1−δ .
N

(b) Now we switch to the case where both the noise and the design have only bounded
moments. Note that
b xexe θ ∗ kmax ≤ kΣ
b e − Σ e kmax + kΣ e − ΣY x kmax + kΣY x − Σxexe θ ∗ kmax
be −Σ
kΣ
e
e
e
e
Yx
Yx
Yx
Yx
b xexe − Σxexe )θ ∗ kmax = T1 + T2 + T3 + T4 .
+ k(Σ
We bound the four terms one by one. For 1 ≤ j ≤ d, analogous to (A.7),
Var(Yei x
eij ) ≤ E(Yei x
eij )2 ≤ E(Yi xij )2 ≤

q
E Yi4 E x4ij =: v1 < ∞.

In addition, E |Yei x
eij |p ≤ (τ1 τ2 )p−2 v1 . Therefore according to Bernstein’s Inequality (Theorem
2.10 in Boucheron et al. (2013)), we have
r



P |b
σYe xej − σYe xej | ≥
where σ
bYe xej =

1
N

2v1 t cτ1 τ2 t 
+
≤ exp(−t),
N
N

PN e
eij , σYe xej = E Yei x
eij and c is certain constant. Then by the union
i=1 Yi x

bound, we have


P |T1 | >

r

2v1 t cτ1 τ2 t 
+
≤ d exp(−t).
N
N

39

Next we bound T2 . Note that for 1 ≤ j ≤ d,
E Yei x
eij − E Yi xij = E Yei x
eij − E Yei xij + E Yei xij − E Yi xij = E Yei (e
xij − xij ) + E(Yei − Yi )xij
q
q


≤ E Yi2 (e
xij − xij )2 P (|xij | ≥ τ2 ) + E (Yei − Yi )2 x2ij P (|Yi | ≥ τ1 )
1
p
1
≤ M v1 2 + 2 ,
τ2
τ1
which delivers that T2 ≤

E

xij
(e
xTi θ ∗ )e



−E

√
M v1 (1/τ12 + 1/τ22 ). Then we bound T3 . For 1 ≤ j ≤ d,

(xTi θ ∗ )xij



=

d
X

E

xik x
eij
θk∗ (e



− xik xij ) ≤

k=1

≤

d
X

d
X

xij x
eik − xij xik |
|θk∗ | E |e

k=1


|θk∗ | E |xij (e
xik − xik )| + E |(e
xij − xij )xik |

k=1

≤

d
X



 CR
|θk∗ | E |xij (e
xik − xik )|1{|xik |>τ2 } + E |(e
xij − xij )xik |1{|xij |>τ2 } ≤ 2 .
τ2
k=1

Finally we bound T4 . For 1 ≤ j, k ≤ d, we have |e
xik x
eij | ≤ τ22 and Var(e
xij x
eik ) ≤ M =: v2 .
Therefore according to Bernstein’s inequality (Theorem 2.10 in Boucheron et al. (2013)),
r



P |b
σxej xek − σxej xek | ≥
where σ
bxej xek = (1/N )

PN

i=1

(A.8)

x
eij x
eik , σxej xek = E x
eij x
eik and c is certain constant. We therefore

have


2v2 t cτ22 t 
+
≤ exp(−t),
N
N

P |T4 | ≥ R

r

2v2 t cτ22 t 
+
≤ d2 exp(−t).
N
N

Now we choose τ1 , τ2  (N/ log d)1/4 , then it follows that for some constant ν2 > 0,
r



be −Σ
b xexe θ ∗ kmax > ν2
P kΣ
e
Yx

δ log d 
≤ 2d1−δ .
N

ei = xi in
Proof of Lemma 2. (a) We first consider the sub-Gaussian design. Again, since x
b xx to be the diagonal of
this case, we do not add any tilde above xi in this proof. Define D

40

b xx . Note that given kxi kψ2 ≤ κ0 and λmin (Σxx ) ≥ κL > 0, it follows that
Σ
∀v ∈ S d−1 ,

p
E(vT xi )4 ≤ 4κ20 ≤ (4κ20 /κL )vT Σxx v.
1

In addition, since kxij kψ2 ≤ κ0 , (E x8ij ) 4 ≤ 8κ20 ≤ (8κ20 /κL ) E x2ij . According to Lemma 5.2
in Oliveira (2016), for any η1 > 1,


1−η1
b xx v ≥ 1 vT Σxx v − c1 (1 + 2η1 log d) kD
b 1/2 vk2 ≥ 1 − d
P ∀ v ∈ Rd : vT Σ
,
xx
1
2
N
3

(A.9)

where c1 is some universal constant. For any 1 ≤ j ≤ d, we know that kxij kψ1 ≤ 2kxij k2ψ2 =
2κ20 , therefore by the Bernstein-type inequality we have for sufficiently small t,
n
 1 X

P
x2ij − E x2ij ≥ t ≤ 2 exp(−cN t2 ),
N i=1

where c depends on κ0 . Note that E x2ij ≤ 2κ20 . An application of the union bound delivers
that when t is sufficiently small,
b xx kop ≥ 2κ2 + t) ≤ 2d exp(−cN t2 ).
P (kD
0
Let t < κ20 , then the inequality above yields that for a new constant C2 > 0,
b xx kop ≥ 3κ2 ) ≤ 2d exp(−C2 N ).
P (kD
0
Combining the inequality above with (A.9), it follows for a new constant C1 > 0 that


1−η1
b xx v ≥ 1 vT Σxx v − C1 (1 + 2η1 log d) kvk2 ≥ 1 − d
P ∀ v ∈ Rd : vT Σ
− 2d exp(−C2 N ).
1
2
N
3
(b) Now we switch to the case of designs with only bounded moments. We will show
that the sample covariance of the truncated design has the restricted eigenvalue property.
Recall that for i = 1, ..., N, j = 1, ..., d, x
eij = sgn(xij )(|xij | ∧ τ2 ). Note that
b xexe v = vT (Σ
b xexe − Σxexe )v + vT (Σxexe − Σxx )v + vT Σxx v.
vT Σ

41

According to (A.8), we have for some c1 > 0,
r



b xexe − Σxexe kmax ≥
P kΣ

2M t c1 τ22 t 
+
≤ d2 exp(−t).
N
N

Given τ2  (N/ log d)1/4 , we have for some constant c2 > 0 and any η2 > 2,
r



b xexe − Σxexe kmax ≥ c2
P kΣ

η2 log d 
≤ d2−η2 .
N

In addition, for any 1 ≤ j1 , j2 ≤ d, we have

E xij1 xij2 − E x
eij1 x
eij2 ≤ E |xij1 xij2 |(1{|xij1 |≥τ2 } + 1{|xij2 |≥τ2 } )
q
√ q
≤ M ( P (|xij1 | ≥ τ2 ) + P (|xij2 | ≥ τ2 ))
s
s
r
√  E x4ij1
E x4ij2 
log d
+
≤ c3
≤ M
4
4
τ2
τ2
N
p
for some c3 > 0, which implies that kΣxx − Σxexe kmax = OP ( log d/N ). Therefore we have
r



b xexe v ≥ vT Σxx v − (c2 η2 + c3 )
P ∀v ∈ Rd , vT Σ

b 2 ≥ c1 √ ρ
Proof of Theorem 2. Suppose k∆k

 q
λN 1− 2
κL


log d
kvk21 ≤ d2−η2 .
N

for some c1 > 0, then by (A.5) and

(A.6),

q
√
√ q
b 1 ≤ 4 2rk∆k
b 2 + 4τ 1−q ρ ≤ 4 2κ 2 √ρλ− 2 k∆k
b 2 + 4κq−1 ρλ1−q
k∆k
L
N
L
N
(A.10)
q
q
√
−
−1 √
b 2.
≤ (4 2 + 4c1 ) ρκL2 λN 2 k∆k
p
q
For any δ > 1, since λN = 2ν δ log d/N , it is easy to verify that as long as ρ(log d/N )1− 2 ≤

C1 for some constant C1 , we will have by (A.10) and (3.3) in Lemma 2
 T

b Σ
b 2 ≥ 1 − d1−δ .
b xx ∆
b ≥ κL k∆k
P ∆
2
4

42

An application of Theorem 1 delivers that for constants c2 , c3 > 0,
q
√  λN 1− 2
b
k∆k2 ≤ c2 ρ
κL

b 2 ≤ c1 √ ρ
When k∆k

 q
λN 1− 2
,
κL

 λ  1−q
2
N
b
and k∆k1 ≤ c3 ρ
.
κL

b 1 ≤ c4 ρ
we can still obtain the `1 norm bound that k∆k



λN
κL

 1−q
2

for some constant c4 through (A.5) and (A.6). Overall, we can achieve the conclusion for
some constants C2 and C3 that with probability at least 1 − 3d1−δ ,
b 2 ≤ C2 ρ
k∆k
2

 δ log d 1− 2q
N


 1−q
b 1 ≤ C3 ρ δ log d 2 .
and k∆k
N

When the design only satisfies bounded moment conditions, again we first assume that
b 2 ≥ c5 √ρ(λN /κL )1− 2q for some constant c5 . Analogous to the case of the sub-Gaussian
k∆k
design, it is easy to verify that for any δ > 1, as long as ρ(log d/N )

1−q
2

≤ C4 for some constant

C4 , we will have by (A.10) and (3.3) in Lemma 2
 T

b Σ
b xx ∆
b ≥ κL k∆k
b 2 ≥ 1 − d1−δ ,
P ∆
2
2
An application of Lemma 1 delivers that for some constants c6 , c7 > 0,
q

 1− 2
b 2 ≤ c7 √ ρ λ N
k∆k
κL

 
b 1 ≤ c8 √ρ λN
and k∆k
κL

1−q
2

.

b 2 ≤ c5 √ρ(λN /κL )1− 2q , we have exactly the same steps as those in the case of the
When k∆k
sub-Gaussian design.
Proof of Lemma 3. Recall that Yei = sgn(Yi )(|Yi | ∧ τ ), then
N
1 X
b
hXi , Θ∗ iXi kop
kΣY X (τ ) −
N i=1
n

(A.11)

X
b Y X (τ ) − E(Ye X)kop + kE((Ye − Y )X)kop + k 1
≤ kΣ
hXi , Θ∗ iXi − E Y Xkop .
n i=1
Next we use the covering argument to bound each term of the RHS. Let S d−1 = {u ∈ Rd :
kuk2 = 1}, N d−1 be the 1/4−net on S d−1 and Φ(A) = supu∈N d1 −1 ,v∈N d2 −1 uT Av for any

43

matrix A ∈ Rd1 ×d2 , then we claim
kAkop ≤

16
Φ(A).
7

(A.12)

To establish the claim above, note that by the definition of the 1/4−net, for any u ∈ S d1 −1
and v ∈ S d2 −1 , there exist u1 ∈ N d1 −1 and v1 ∈ N d2 −1 such that ku − u1 k2 ≤ 1/4 and
kv − v1 k2 ≤ 1/4. Then it follows that
uT Av = uT1 Av1 + (u − u1 )T Av1 + uT1 A(v − v1 ) + (u − u1 )T A(v − v1 )
1 1
1
≤ Φ(A) + ( + + )
sup
uT Av.
4 4 16 u∈S d1 −1 ,v∈S d2 −1
Taking the superlative over u ∈ S d1 −1 and v ∈ S d2 −1 on the LHS yields (A.12).
Now fix u ∈ S d1 −1 and v ∈ S d2 −1 . Later on we always write uT Xi v as Zi and uT Xv as
Z for convenience. Consider
n

b Y X (τ ) − E(Ye X))v =
uT (Σ

1Xe
Yi Zi − E(Ye Z).
n i=1

Note that




E(Yei Zi − E Yei Zi )2 ≤ E(Yei Zi )2 ≤ E(Yi Zi )2 = E(Zi2 E(Yi2 |Xi )) = E Zi2 (hXi , Θ∗ i)2 + E Zi2 E(2i |Xi )

k  k1
1
2 4
2
2
≤ 16R2 κ40 + 2kκ20 M k /(k − 1) < ∞,
≤ 16R κ0 + (2kκ0 /(k − 1)) E E(i |Xi )
(A.13)
which we denote by v1 for convenience. Also we have

E |Yei Zi |p ≤ τ p−2 E(Yi2 |Zi |p ) = τ p−2 E hXi , Θ∗ i2 |Zi |p + 2i |Zi |p



pk 1− 1
p
 k1  k−1
k
∗ 4
p−2
2
k
2p
≤τ
E xij
EhXi , Θ i E |Zi | + E E(i |Xi )
r

p p
1
pk p 
p−2
2 2
k
≤τ
4R κ0 (κ0 2p) + M κ0
k−1

44

and it holds for constants c1 and c2 that
E |Yei Zi − E Yei Zi |p ≤ 2p−1 (E |Yei Zi |p + | E Yei Zi |p ) ≤ 2p−1 (E |Yei Zi |p + (E |Yi Zi |)p )

(A.14)

≤ c1 p!(c2 τ )p−2 .

where the last inequality uses the Stirling approximation of p!. Define v := c1 ∨ v1n
. Then
oNan
application of Bernstein’s Inequality (Theorem 2.10 in Boucheron et al. (2013)) to Yei Zi
i=1

delivers

r
N
 1 X
2vt c2 τ t 
Yei Zi − E Yei Zi | >
+
≤ exp(−t).
P |
N i=1
N
N

By taking the union bound over all (u, v) ∈ N d1 −1 × N d2 −1 and (A.12) it follows that


b e (τ ) − E(Ye X)kop
P kΣ
YX

16 
≥
7

r

2vt c2 τ t 
+
≤ exp ((d1 + d2 ) log 8 − t) ,
N
N

(A.15)

where c1 is a constant.
Next we aim to bound kE((Ye − Y )X)kop . For any u ∈ S d1 −1 and v ∈ S d2 −1 , by the
Cauchy-Schwartz inequality and the Markov inequality,
r
E((Ye − Y )Z) ≤ E(|Y Z| · 1|Y |>τ ) ≤

p
E(Y 2 Z 2 )P (|Y | > τ ) ≤

vEY 2
≤
τ2

q
1
v(2R2 κ20 + M k )
τ

Note that the inequality above holds for any (u, v). Therefore

kE((Ye − Y )X)kop ≤

q
1
v(2R2 κ20 + M k )
τ

.

(A.16)

Now we give an upper bound of the third term on the RHS of (A.11). For any u ∈ S d1 −1
and v ∈ S d2 −1 , khXi , Θ∗ iZi kψ1 ≤ Rκ20 , so by Proposition 5.16 (Bernstein-type inequality) in
Vershynin (2010) it follows that for sufficiently small t,


N
 1 X

c3 N t 2
∗
P
hXi , Θ iZi − E Y Z > t ≤ 2 exp − 2 4 ,
N i=1
R κ0
where c3 is a constant. Then an combination of the union bound over all points on N d1 −1 ×

45

.

N d2 −1 and (A.12) delivers
N
 1 X

16 
c3 N t 2 
P k
hXi , Θ∗ iXi − E Y Xkop ≥ t ≤ 2 exp (d1 + d2 ) log 8 − 2 4 ,
N i=1
7
R κ0

Finally we choose τ 

p

(A.17)

N/(d1 + d2 ). Combining (A.15), (A.16) and (A.17), we can find

a constant γ > 0 such that as long as (d1 + d2 )/N < γ, it holds that
r
N

X
1
d1 + d2 
b Y X (τ ) −
P kΣ
hXi , Θ∗ iXi kop > ν
≤ η exp(−(d1 + d2 )),
N i=1
N
where ν and η are constants.

Proof of Theorem 3. We first verify the RSC property. According to Proposition 1 in Negahban and Wainwright (2011), the following inequality holds for all ∆ ∈ Rd1 ×d2 with probability
at least 1 − 2 exp(−N/32),
r 
q
r d
p
1
d2
1
+
k∆kN .
vec(∆)T ΣXX vec(∆) ≥ k ΣXX vec(∆)k2 − c
4
N
N

(A.18)

b − Θ∗ as ∆.
b Suppose
Let κL = (1/32)λmin (ΣXX ) > 0. For ease of notations, write Θ
q


1− 2
b F ≥ c1 √ρ λN
k∆k
for some c1 > 0, then by (A.5) and (A.6),
κL
q
q
p
√
b F + 4κq−1 ρλ1−q
b N ≤ 4 2rk∆k
b F + 4τ 1−q ρ ≤ 4 2ρκ 2 λ− 2 k∆k
k∆k
L N
L
N
q
√
− 2q
−1 √
2
b F.
≤ (4 2 + 4c1 ) ρκL λN k∆k

(A.19)

p
Since we choose λN = 2ν (d1 + d2 )/N , there exists a constant C1 such that as long as
1− 2q
ρ (d1 + d2 )/N
≤ C1 , we have by (A.18) and (A.19)
b TΣ
b XX vec(∆)
b ≥ κL k∆k
b 2
vec(∆)
F
with probability at least 1 − 2 exp(−N/32). An application of Theorem 1 delivers that for

46

some constants c2 , c3 > 0, it holds with high probability that
q
√  λN 1− 2
b
k∆kF ≤ c2 ρ
κL

b F ≤ c1 √ ρ
When k∆k

 q
λN 1− 2
,
κL

 λ  1−q
2
N
b
and k∆kN ≤ c3 ρ
.
κL

b 1 ≤ c4 ρ
we can still obtain the `1 norm bound that k∆k



λN
κL

 1−q
2

through (A.5) and (A.6), where c4 is some constant. Overall, we can achieve the conclusion
that,
b 2 ≤ C2 ρ
k∆k
F

 d + d 1− 2q
1
2
N

b N ≤ C3 ρ
and k∆k

 d + d  1−q
2
1
2
N

with probability at least 1 − C4 exp(−(d1 + d2 )) for constants C2 , C3 and C4 .

Proof of Lemma 4. We follow essentially the same strategies as in proof of Lemma 3. The
only difference is that we do not use the covering argument to bound the first term in (A.11).
Instead we apply the Matrix Bernstein inequality (Theorem 6.1.1 in Tropp (2015)) to take
advantage of the singleton design under the matrix completion setting.
For any fixed u ∈ S d1 −1 and v ∈ S d2 −1 , write uT Xi v as Zi . Then we have
E(Yei Zi ) =

d1 X
d2
p

1 X
e
d1 d2 E(Yi uj(i) vk(i) ) = √
E |Yei | j(i) = j0 , k(i) = k0 · |uj0 vk0 |
d1 d2 j =1 k =1
0

d1
X

0

d2
X


1
E |Yi | j(i) = j0 , k(i) = k0 · |uj0 vk0 |
≤√
d1 d2 j =1 k =1
0
0
q
q
d1 X
d2
X
1
1
1
|uj0 vk0 | ≤ R2 + M k .
≤ R2 + M k · √
d1 d2 j =1 k =1
0

0

Since the above argument holds for all u ∈ S d1 −1 and v ∈ S d2 −1 , we have kE(Yei Xi )kop ≤
p
1
R2 + M k . In addition,
kE Yei2 XTi Xi kop = d1 d2 kE Yei2 ek(i) eTj(i) ej(i) eTk(i) kop = d1 d2 kE Yei2 ek(i) eTk(i) kop


= d1 d2 E E(Yei2 |Xi )ek(i) eTk(i)
=

=
op

d2 X
d1


X
E Yei2 |k(i) = k0 , j(i) = j0 ek0 eTk0
k0 =1 j0 =1

d1


X
1
max
E Yei2 |k(i) = k0 , j(i) = j0 ≤ d1 E Yi2 ≤ d1 (R2 + M k )

k0 =1,...,d2

j0 =1

47

op

1
Similarly we can get kE Yei2 Xi XTi kop ≤ d2 (R2 + M k ). Write Yei Xi − E Yei Xi as Ai . Therefore
1

by the triangle inequality, max(kE ATi Ai kop , kE Ai ATi kop ) ≤ (d1 ∨ d2 )(R2 + M k ), which we
√
denote by v for convenience. Since kAi kop ≤ 2 d1 d2 τ , an application of the Matrix Bernstein
inequality delivers,



b e (τ ) − E(Ye X)kop ≥ t ≤ (d1 + d2 ) exp
P kΣ
YX


−N t2 /2
√
,
v + d1 d2 τ t/3

(A.20)

Next we aim to bound kE((Yei − Yi )Xi )kop . Fix u ∈ S d1 −1 and v ∈ S d2 −1 . By the
Cauchy-Schwartz inequality and the Markov inequality,
E((Yei − Yi )Zi ) =

p

=√

d1 d2 E((Yei − Yi )uj(i) vk(i) )

d1 X
d2

1 X
E |Yei − Yi | j(i) = j0 , k(i) = k0 · |uj0 vk0 |
d1 d2 j =1 k =1
0

1
≤√
d1 d2

0

d1
X

d2
X


E |Yi | · 1|Yi |>τ j(i) = j0 , k(i) = k0 · |uj0 vk0 |

j0 =1 k0 =1

1
1
d1 X
d2
1 X
R2 + M k
R2 + M k
·√
.
|uj0 vk0 | ≤
≤
τ
τ
d1 d2 j =1 k =1
0

0

Note that the inequality above holds for any (u, v). Therefore
1

R2 + M k
kE((Yei − Yi )Xi )kop ≤
.
τ

(A.21)

Now we give an upper bound of the third term on the RHS of (A.11). Denote hXi , Θ∗ iXi −
√
EhXi , Θ∗ iXi by Bi . It is not hard to verify that kBi kop ≤ 2R d1 d2 and max(kE BTi Bi kop ,
kE Bi BTi kop ) ≤ (d1 ∨ d2 )R2 , it follows that for any t ∈ R,
N
 1 X
P
hXi , Θ∗ iXi − E Yi Xi
N i=1

Finally, choose τ =

p

op



≥ t ≤ (d1 + d2 ) exp


−N t2 /2
√
.
(d1 ∨ d2 )R2 + 2R d1 d2 t/3
(A.22)

N/((d1 ∨ d2 ) log(d1 + d2 )). Combining (A.20), (A.21) and (A.22), it

is easy to verify that for any δ > 0, there exist constants ν and γ such that the conclusion
holds as long as (d1 ∨ d2 ) log(d1 + d2 )/N < γ.

48

Proof of Theorem 4. This proof essentially follows the proof of Lemma 1 in Negahban and
Wainwright (2012). Write (d1 + d2 )/2 as d. Define a constraint set
n
p
1
k∆kmax · k∆kN
C(N ; c0 ) = ∆ ∈ Rd1 ×d2 , ∆ 6= 0| d1 d2
≤
2
k∆kF
c0 L

s

o
N
.
(d1 + d2 ) log(d1 + d2 )

b ∈
According to Case 1 in the proof of Lemma 1 in Negahban and Wainwright (2012), if ∆
/
C(N ; c0 ), we have
r
b 2 ≤ 2c0 R
k∆k
F

dX
1 ∧d2
d log d √ b
{8 rk∆kF + 4
σj (Θ∗ )},
N
j=r+1

Following the same strategies in the proof of Theorem 1 in our work, we have
r

1− 2q
√
b F ≤ c1 ρ 2c0 R d log d
k∆k
N
b ∈ C(N ; c0 ), according to Case 2 in proof of Lemma 1 in Negahban
for some constant c1 . If ∆
b F ≤
and Wainwright (2012), with probability at least 1 − C1 exp(−C2 d log d), either k∆k
√
b TΣ
b XX vec(∆)
b ≥ (1/256)k∆k
b 2 , where C1 and C2 are certain constants.
512R/ N , or vec(∆)
F
√
b
For the case where k∆kF ≤ 512R/ N , combining this fact with (A.5) and (A.6) delivers
that
p
b N ≤ 4 2ρτ − 2q 512R
√ + 4τ 1−q ρ.
k∆k
N
We minimize the RHS of the inequality above by plugging in τ =
for some constant c3 > 0,

R2
ρN

1
 2−q

. Then we have


 1
2
b N ≤ c3 ρ R 1−q 2−q .
k∆k
N

b TΣ
b XX vec(∆)
b ≥ (1/256)k∆k
b 2 , it is implied by Theorem 1 in our
For the case where vec(∆)
F
work that

q
b F ≤ c4 √ρλ1− 2
k∆k
N

for some constants c4 and c5 . Since λN

1−q

b N ≤ c5 ρλ 2
and k∆k
N
p
= ν δ(d1 ∨ d2 ) log(d1 + d2 )/N and R ≥ 1, by

49

Lemma 4, it holds with probability at least 1 − C1 exp(−C2 (d1 + d2 )) − 2(d1 + d2 )1−δ that
n  δR2 (d + d ) log(d + d ) 1− 2q R2 o
1
2
1
2
2
b
k∆kF ≤ C3 max ρ
,
N
N
and

o
n  δR2 (d + d ) log(d + d )  1−q
1
ρR2−2q  2−q
2
1
2
1
2
b
k∆kN ≤ C4 max ρ
,
N
N 1−q

for constants C3 and C4 .

ei = xi , we do not
Proof of Lemma 5. (a) First we study the sub-Gaussian design. Since x
ei , then we have
add tildes above x, xi or xij in this proof. Denote yi (kyi k2 ∧ τ )/kyi k2 by y
n

n

X
X
b xey (τ ) − 1
b yex (τ ) − 1
kΣ
xj xTj Θ∗ kop = kΣ
Θ∗T xj xTj kop
n j=1
n j=1

(A.23)

b yex (τ ) − E(e
b xx Θ∗ − E(yxT )kop
≤ kΣ
yxT )kop + kE(e
yxT − yxT )kop + kΣ
Let

"
Si =

0

ei xTi
ei xTi − E y
y

eiT
eiT − E xi y
xi y

0

#
.

Now we bound kE Spi kop for p > 2. When p is even, i.e., p = 2l, we have
"
S2l
i =

#

e i − E xi y
ei ))l
ei xTi )(xi y
((e
yi xTi − E y

0

0

ei xTi ))l
e i − E xi y
ei )(e
((xi y
yi xTi − E y

.

For any v ∈ S d2 −1 ,
eiT )l )v = E((xTi xi )l (e
ei )l−1 (vT y
ei )2 ) ≤ τ 2l−2 E((xTi xi )l (vT y
ei )2 )
vT (E(e
yi xTi xi y
yiT y
≤ τ 2l−2 E((xTi xi )l (vT yi )2 ) = τ 2l−2 E E[(vT yi )2 |xi ](xTi xi )l


T
= τ 2l−2 E (vT Θ∗ xi )2 (xTi xi )l + E[(vT i )2 |xi ](xTi xi )l .



(A.24)

Also note that for any v ∈ S d2 −1 , we have
Rkvk22 ≥ vT E yi yiT v ≥ (Θ∗ v)T Σxx (Θ∗ v) ≥ κL kΘ∗T vk22 .

50

(A.25)

Therefore it follows that kΘ∗ kop ≤

p
R/κL . Combining this with the fact that kxTi xi kψ1 ≤

2d1 κ20 , it follows that


T

∗T

E (v Θ xi )

2

(xTi xi )l



q
q
p
∗T
T
4
≤ E(v Θ xi ) E(xTi xi )2l ≤ (2 R/κL κ0 )2 (4d1 κ20 l)l .

Also by Hölder Inequality, we have
T

2

E E[(v i )

|xi ](xTi xi )l



T

2

≤ E E[(v xi ) |xi ]
1

≤ M k (2d1 κ20

k

 k1 

lk
E[(xTi xi ) k−1 ]

1− k1

lk l
).
k−1

Therefore we have
eiT )l )kop ≤ c01 (c1 τ )2l−2 (d1 l)l
kE(e
yi xTi xi y

(A.26)

for l ≥ 1, where c1 and c01 are some constants. Letting l = 1 in the equation above imei xTi E xi y
eiT kop ≤ kE y
ei xTi xi y
eiT kop ≤ c01 d1 . Therefore it holds that (e
plies that kE y
yi xTi −
ei xTi )(xi y
e i − E xi y
eiT )  2e
eiT + 2 E y
ei xTi E xi y
eiT  2ỹi xTi xi ỹiT + 2c01 d1 I. In addiEy
yi xTi xi y
tion, for any commutable positive semi-definite (PSD) matrices A and B, it is true that
(A + B)l  2l−1 (Al + Bl ) for l > 1. Then it follows that
ei xTi )(xi y
eiT − E xi y
eiT ))l  (2e
eiT + 2c01 d1 I)l
((e
yi xTi − E y
yi xTi xi y
eiT )l + (c01 d1 )l I).
 22l−1 ((e
yi xTi xi y

(A.27)

ei xTi )(xi y
eiT − E xi y
eiT ))l kop ≤ c02 (c2 τ )2l−2 (d1 l)l for some
Therefore we have kE((e
yi xTi − E y
eiT − E xi y
eiT )(e
constants c2 and c02 . Using similar methods, we can derive kE((xi y
yi xTi −
ei xTi ))l kop ≤ c03 (c3 τ )2l−2 (d2 l)l for some constants c3 and c03 . So we can achieve for some
Ey
constants c4 and c04 ,
0
kE S2l
i kop ≤ c4 (c4 τ

p
d1 ∨ d2 )2l−2 ll (d1 ∨ d2 )

(A.28)

When p is an odd number, i.e., p = 2l + 1, by (A.27) we have
q
S4l+2
i
q
q
2l+1/2
T
T
T
T
2l+1
2l+1
ei )
e i xi E x i y
ei )
eiT y
ei xTi )2l+1 + (E xi y
eiT E y
e i xi y
eiT )2l+1 )
2
diag( (e
yi x i xi y
+ (E y
, (xi y

eiT )l+1/2 + (E y
ei xTi E xi y
eiT )l+1/2 , (xi y
eiT y
ei xTi )l+1/2 + (E y
ei xTi E xi y
eiT )l+1/2 .
 22l+1/2 diag (e
yi xTi xi y

S2l+1
=
i

51

1

1

eiT y
ei xTi xi y
ei xTi ) 2 = xi y
eiT y
eiT /(kxi k2 ke
ei xTi /(kxi k2 ke
eiT ) 2 = y
Note that (xi y
yi k2 ),
yi k2 ) and (e
yi xTi xi y
so for any v ∈ S d2 −1 , we have
1

1

1

eiT )l+1/2 v = E((xTi xi )l+ 2 (e
ei )2 ) ≤ τ 2l−1 E((xTi xi )l+ 2 (vT y
ei )2 ).
ei )l− 2 (vT y
vT E(e
yi xTi xi y
yiT y
eiT )l+1/2 )kop ≤
Following the same steps as in the case of l being even, we can derive kE(e
yi xTi xi y
c05 (c5 τ )2l−1 (d1 (l + 1/2))l+1/2 , kE(xi ỹiT ỹi xTi )l+1/2 )kop ≤ c06 (c6 τ )2l−1 (d2 (l + 1/2))l+1/2 and finally
kE S2l+1
kop ≤ c07 (c7 τ
i

p
d1 ∨ d2 )2l−1 (l + 1/2)l+1/2 (d1 ∨ d2 ).

(A.29)

Define σ 2 = (c7 ∨ c07 )(d1 ∨ d2 ). Then by (A.28), we have kE S2i kop ≤ σ 2 . Also by combining
(A.28) and (A.29), we can verify the moment constraints of Lemma 6; for p > 2, we have
kE Spi kop ≤ p!((c4 ∨ c7 )τ

p
d1 ∨ d2 )p−2 σ 2 .

(A.30)

Therefore by Lemma 6, we have
n


 1X

b yex (τ ) − E(e
P kΣ
yxT )kop ≥ t ≤ P k
Si kop ≥ t
n i=1


nt2
nt
√
.
≤ (d1 + d2 ) · exp − c min 2 ,
σ (c4 ∨ c7 )τ d1 ∨ d2

(A.31)

Next we aim to bound kE(e
yi xTi − yi xTi )kop . Note that for any u ∈ S d1 −1 and v ∈ S d2 −1 ,
T

2

T

2

T

2

T

2



T

∗

2

T

2

E((v yi ) (u xi ) ) = E(E((v yi ) |xi )(u xi ) ) = E E (v Θ xi ) (u xi )



T

2

T

+ E((v i ) |xi )(u xi )

2



:= v 2 < ∞.
For any u ∈ S d1 −1 and v ∈ S d2 −1 , by the Cauchy-Schwartz inequality and the Markov
inequality,
p
ei )(uT xi ) − (vT yi )(uT xi )) ≤ E(|(vT yi )(uT xi )| · 1kyk2 >τ ) ≤ E(|(vT yi )(uT xi )|2 )P (kyk2 > τ )
E((vT y
r
√ q
q
v 2 Ekyk22
v
v
d2
1
≤
=
R2 Ekxk22 + Ekk22 ≤
R2 κ20 + M k .
2
τ
τ
τ

52

Since u and v are arbitrary, we have
kE(e
yi xTi

−

yi xTi )kop

√ q
v d2
1
≤
R2 κ20 + M k .
τ

(A.32)

Now we give an upper bound of the third term on the RHS of (A.23). For this term, we
still follow the covering argument. Let N d−1 be the 1/4−net on S d−1 . For any u ∈ S d1 −1 and
p
v ∈ S d2 −1 , kuT xi xTi Θ∗ vkψ1 ≤ R/κL κ20 , so by Proposition 5.16 (Bernstein-type inequality)
in Vershynin (2010) it follows that for sufficiently small t,
P





b xx Θ∗ v − E(uT xi )(vT yi ) > t ≤ 2 exp −c8 nt2 ,
uT Σ

where c8 is some positive constant. Then an combination of the union bound over all points
on N d1 −1 × N d2 −1 and (A.12) delivers



b xx Θ∗ − E yxT k2 ≥ 16 t ≤ 2 exp (d1 + d2 ) log 8 − c8 nt2 ,
P kΣ
7

(A.33)

where c8 is a constant.
Finally we choose τ = O(

p
n/ log(d1 + d2 )), and combining (A.31), (A.32) and (A.33)

delivers that for any δ > 0, as long as (d1 + d2 ) log(d1 + d2 )/n < γ for some constant γ > 0,
n

X
b xey (τ1 ) − 1
P kΣ
Θ∗T xj xTj kop ≥
n j=1


r

(ν + δ)(d1 + d2 ) log(d1 + d2 ) 
≤ 2(d1 + d2 )1−ηδ ,
n

where ν and η are universal constants.
(b) Now we switch to the case of designs with only bounded moments. We first show
that for two constants c3 and c4 ,

b xeye − Σxy kop ≥
P kΣ

r

(c3 + δ)(d1 + d2 ) log(d1 + d2 ) 
≤ (d1 + d2 )1−c4 δ .
n

Note that

E(vT yi )4 = E(vT Θ∗T xi + vT i )4 ≤ 16 E(vT Θ∗T xi )4 + E(vT i )4 .

53

(A.34)

By Rosenthal-type inequality in Ibragimov and Sharakhmetov (1998), we have
E(vT yi )4 ≤ 16(16R4 κ40 + M ) =: V1 .
e i := x
ei y
eiT .
In addition, E(vT xi )4 ≤ 16κ40 =: V2 . Let V := max(V1 , V2 ), Zi := xi yiT and X
Following the proof strategy for Theorem 6, we have
ei − E Z
e i kop ≤ kZ
e i kop + kE Zi kop = ke
kZ
xi k2 ke
y i k2 +

√
√
1
R ≤ (d1 d2 ) 4 τ 2 + V .

Also for any v ∈ S d−1 , we have
eT Z
e
ei )2 ) ≤ E(kxi k22 (vT yi )2 )
E(vT Z
xi k22 (vT y
i i v) = E(ke
d
d q
X
X
2
T
2
=
E(xij (v yi ) ) ≤
E(x4ij ) E(vT yi )4 ≤ V d1 .
j=1

j=1

e
e eT
eT Z
Then it follows that kE Z
i i kop ≤ V d1 . Similarly we can obtain kE Zi Zi kop ≤ V d2 . Denote
e i )T E X
e i kop , k(E Z
e i) E Z
e T kop ) ≤ kE Zi k2 ≤ V ,
d1 + d2 by d. Since max(k(E X
i

op

ei − E Z
e i )T (Z
ei − E Z
e i ))kop , kE((Z
ei − E Z
e i )(Z
ei − E Z
e i )T )kop ) ≤ V (d + 1).
max(kE((Z
By Corollary 6.2.1 in Tropp (2012), we have for some constant c1 ,
P

!
n


1 Xe
nt
nt2
e i kop > t ≤ d exp −c1
√
k
∧√
.
Zi − E Z
n i=1
V (d + 1)
dτ1 τ2 + V

(A.35)

e i kop . For any v ∈ S d−1 , it holds that
Now we bound the bias EkZi − Z
e i )v) = E(|(vT xi )(vT yi ) − (vT x
ei )(vT y
ei )|1{kxi k4 ≥τ1
E(vT (Zi − Z

or kyi k4 ≥τ2 } )

≤ E(|(vT xi )(vT yi )|(1{kxi k4 >τ1 } + 1kyi k4 >τ2 ))
p
p
≤ E((vT xi )2 (vT yi )2 )P (kxi k4 > τ1 ) + E((vT xi )2 (vT yi )2 )P (kyi k4 > τ2 )
√
√
V d1 V d2
≤
+
.
τ12
τ22
e i )kop ≤ V
Therefore we have kE(Zi − Z

√

1

d(1/τ12 + 1/τ22 ). Choose τ1  (n/ log d1 ) 4 , τ2 
54

1

(n/ log d2 ) 4 and substitute t with

p
δd log d/n. Then we obtain (A.34) by combining the

bound of bias and (A.35).
Finally, note that
b xeye − Σ
b xexe Θ∗ kop ≤ kΣ
b xeye − Σxy kop + k(Σ
b xexe − Σxx )Θ∗ kop + kΣxy − Σxx Θ∗ kop .
kΣ
b xexe established in Theorem 6, Condition (C1) and (A.34), we can
Combining consistency of Σ
reach the conclusion of the lemma.
Lemma 6 (Matrix Bernstein Inequality with Moment Constraint). Consider a finite sequence {Si }ni=1 of independent, random, Hermitian matrices with dimensions d × d. Assume
that E Si = 0 and kE S2i kop ≤ σ 2 . Also the following moment conditions hold for all 1 ≤ i ≤ n
and p ≥ 2:
kE Spi kop ≤ p!Lp−2 σ 2 ,
where L is a constant. Then for every t ≥ 0 we have
n
 1X


 −nt2 
nt2 nt 
P k
≤
d
·
exp
−
c
min
,
.
Si kop ≥ t ≤ d exp
n i=1
4σ 2 + 2Lt
σ2 L

Proof of Lemma 6. Given the moment constraints, we have for 0 < θ < 1/L,

Ee

θSi

=I+

∞
X
θp E Sp
i

p=2

p!

I+

∞
X

2

σ L

θ2 σ 2
θ I=I+
I  exp
1 − θL

p−2 p

p=2



θ2 σ 2
1 − θL


I.

Let g(θ) = θ2 /(1−θL). Owing to the master tail inequality (Theorem 3.6.1 in Tropp (2015)),
we have

n
 X

X

P k
Si kop ≥ t ≤ inf e−θt tr exp
log E eθSi
θ>0

i=1

≤

inf

i
−θt

e



tr exp nσ 2 g(θ)I

0<θ<1/L

≤

inf


de−θt exp nσ 2 g(θ) .

0<θ<1/L

Choosing θ = t/(2nσ 2 + Lt), we can reach the conclusion.
Proof of Theorem 5. (a) First consider the case of the sub-Gaussian design. Denote

1
n

Pn

j=1

xj xTj

by Σxx . Since xj is a sub-Gaussian vector, an application of Theorem 5.39 in Vershynin
55

(2010) implies that with probability at least 1 − exp(−c1 d1 ), kΣxx − Σxx kop ≤ c2

p
d1 /n and

furthermore λmin (Σxx ) ≥ 12 λmin (Σxx ) > 0 as long as d1 /n is sufficiently small. Therefore
d2
N
n
1 X b
1 XX
T b
2
b
b
b xj eT i2
vec(∆) ΣXX vec(∆) =
h∆, Xi i =
h∆,
k
N i=1
N j=1 k=1
d2
d2 X
n
n
1 XX
1 X
T b
2
b k )2
=
(xT ∆e
Tr(xj ∆ek ) =
N j=1 k=1
N k=1 j=1 j

=

(A.36)

d2
1 X
b k ) ≥ 1 λmin (Σxx )k∆k
b 2 ,
b k )T Σxx (∆e
(∆e
F
d2 k=1
2d2

with probability at least 1 − exp(−c1 d1 ). In addition, by Lemma 5, as long as we choose
p
b Y X − 1 PN Yi Xi kop with probability
λN = d22 (ν1 + δ)(d1 + d2 ) log(d1 + d2 )/n, λN ≥ 2kΣ
i=1
N
at least 1 − 2(d1 + d2 )1−η1 δ . Finally, by Theorem 1 we establish the statistical error rate:
there exist γ1 , γ2 > 0 that as long as (d1 + d2 ) log(d1 + d2 )/n < γ1 and d1 + d2 > γ2 , we have
b N , τ ) − Θ∗ k 2 ≤ C 1 ρ
kΘ(λ
F
and

 (ν + δ)(d + d ) log (d + d ) 1−q/2
1
1
2
1
2
n

 (ν + δ)(d + d ) log (d + d )  1−q
2
1
1
2
1
2
∗
b
kΘ(λN , τ ) − Θ kN ≤ C2 ρ
n

for some constants C1 and C2 .
(b) Now we switch to the designs with bounded moments. According to Theorem 6, with
probability at least 1 − d1−Cδ
,
1
r
b xexe − Σxx kop ≤
kΣ

δd1 log d1
,
n

b xexe ) ≥ 1 λmin (Σxx ) > 0 as long as d1 /n is sufficiently
which furthermore implies that λmin (Σ
2
small. Analogous to (A.36), we therefore have
b TΣ
b e e vec(∆)
b ≥
vec(∆)
XX

1
b 2.
λmin (Σxx )k∆k
F
2d2

Combining Theorem 1 with the inequality above, we can reach the final conclusion.

56

e i for ease of notations. Note that
ei x
eTi by X
Proof of Theorem 6. We denote xi xTi by Xi and x
ei − EX
e i kop ≤ kX
e i kop + kE X
e i kop = ke
kX
xi k22 +

√
√
√
R ≤ dτ 2 + R

Also for any v ∈ S d−1 , we have
e i v) = E(ke
eTX
ei )2 ) ≤ E(kxi k22 (vT xi )2 )
xi k22 (vT x
E(vT X
i
d
d q
X
X
=
E(x2ij (vT xi )2 ) ≤
E(x4ij ) E(vT xi )4 ≤ Rd
j=1

j=1

ei −
eTX
e i kop ≤ Rd. Since k(E X
e i )T E X
e i kop ≤ kE Xi k2 ≤ R, kE((X
Then it follows that kE X
i
op
e i )T (X
ei − EX
e i ))kop ≤ R(d + 1). By Theorem 5.29 (Non-commutative Bernstein-type
EX
inequality) in Vershynin (2010), we have for some constant c,
n
 1X



nt
nt2
ei − EX
e i kop > t ≤ 2d exp −c
√
P k
∧√
.
X
n i=1
R(d + 1)
dτ 2 + R

(A.37)

e i kop . For any v ∈ S d−1 , it holds that
Now we bound the bias EkXi − X
e i )v) = E(((vT xi )2 − (vT x
ei )2 )1{kxi k4 >τ } )
E(vT (Xi − X
p
≤ E((vT xi )2 1{kxi k4 >τ } ) ≤ E(vT xi )4 P (kxi k4 > τ )
r
√
R d
R2 d
= 2 .
≤
τ4
τ
√
e i )kop ≤ R d/τ 2 . Choose τ  (nR/δ log d) 41 and substitute t
Therefore we have kE(Xi − X
p
with δRd log d/n. Then we reach the final conclusion by combining the bound of bias and
(A.37).

57

