arXiv:1605.07129v4 [] 15 Aug 2017

Sub-Gaussian estimators of the mean of a random
matrix with heavy-tailed entries
Stanislav Minsker∗
e-mail: minsker@usc.edu
Abstract: Estimation of the covariance matrix has attracted a lot of attention of the statistical research community over the years, partially due to important applications such as
Principal Component Analysis. However, frequently used empirical covariance estimator, and
its modifications, is very sensitive to the presence of outliers in the data. As P. Huber wrote
[Hub64], “...This raises a question which could have been asked already by Gauss, but which
was, as far as I know, only raised a few years ago (notably by Tukey): what happens if the
true distribution deviates slightly from the assumed normal one? As is now well known, the
sample mean then may have a catastrophically bad performance...” Motivated by Tukey’s
question, we develop a new estimator of the (element-wise) mean of a random matrix, which
includes covariance estimation problem as a special case. Assuming that the entries of a
matrix possess only finite second moment, this new estimator admits sub-Gaussian or subexponential concentration around the unknown mean in the operator norm. We explain the
key ideas behind our construction, and discuss applications to covariance estimation and
matrix completion problems.
MSC 2010 subject classifications: Primary 60B20, 62G35; secondary 62H12.
Keywords and phrases: random matrix, heavy tails, concentration inequality, covariance
estimation, matrix completion..

1. Introduction
Let Y1 , . . . , Yn ∈ Cd1 ×d2 be a sequence of independent random matrices such that all their entries
have finite second moments: E |(Yj )k,l |2 < ∞ for all 1 ≤ j ≤ n, 1 ≤ k ≤ d1 , 1 ≤ l ≤ d2 .
Let EY1 , . . . , EYn ∈ Cd1 ×d2 be the expectations evaluated element-wise, meaning that (EYj )k,l =
h P
i
E (Yj )k,l . The goal of this paper is to construct and study estimators of EȲ := E n1 nj=1 Yj
under minimal assumptions on the distributions of Y1 , . . . , Yn . In particular, we are interested
in the estimators that admit tight non-asymptotic bounds and exponential deviation inequalities
without imposing any additional assumptions (besides finite second moments) on Y1 , . . . , Yn . For
example, if Yj = Zj ZjT , where Z1 , . . . , Zn ∈ Rd are i.i.d. copies of a random vector Z such that


EZ = 0, E ZZ T = Σ and EkZk42 < ∞, formulated problem is reduced to covariance estimation
(here and in what follows, k · k2 and h·, ·i stand for the usual Euclidean norm and Euclidean dot
product respectively).
Techniques developed in this paper have direct connection to several problems in high-dimensional
statistics and statistical learning theory. In the past decade, these fields have seen numerous breakthroughs in structural estimation, concerned with a task of recovering a high-dimensional parameter that belongs to a set with “simple” structure from a small number of measurements. Examples
include sparse linear regression, low-rank matrix recovery and structured covariance estimation.
∗

Department of Mathematics, University of Southern California
1

S. Minsker/Estimators of the mean of a random matrix

2

However, theoretical recovery guarantees for popular techniques (e.g., `1 and nuclear norm minimization) usually require strong assumptions on the underlying probability distribution, such
as sub-Gaussian or bounded noise. What happens with the performance of the algorithms when
these conditions are violated, which is the case for many real data sets modeled by heavy-tailed
distributions? Can the assumptions be weakened without sacrificing the quality of theoretical
guarantees? We look at examples where the answer is positive, and describe modifications of
existing techniques that allow to achieve the improvements.
1.1. Overview of the previous work
Let us begin by briefly discussing a scalar version of the problem investigated in this paper. Assume
that X1 , . . . , Xn ∈ R are i.i.d. copies of X, where EX 2 < ∞. One of the fundamental problems
in statistics is to construct the confidence interval for the unknown mean EX based on a given
sample. A surprising fact (dating back to [NY83] where the “median of means” estimator was
introduced, along with [AMS96] and [JVV86]) is that it’s possible to construct
 a non-asymptotic

ˆ
confidence intervals In (δ) with coverage probability 1 − δ (meaning that Pr EX ∈ Iˆn (δ) ≥ 1 − δ
q
p
ˆ
for given n and δ) and “nearly optimal” length In (δ) ≤ L Var(X) log(e/δ)
, where L > 0 is an
n
absolute constant. An in-depth study of this and closely related questions was performed in [Cat12,
DLLO15] based on two different approaches. Note that the
interval
 center of any such confidence
 is
q
p
a point estimator µ̂ := µ̂(X1 , . . . , Xn , δ) that satisfies Pr |µ̂ − EX| ≥ L Var(X) log(e/δ)
≤ δ.
n
Because the only assumption on X is the existence of a second moment, it is natural to call such
an estimator “robust” 1 : it admits strong deviation bounds even for the heavy-tailed distributions
that can be used to model outliers in the data. Ideas behind these results have also been extended
to empirical risk minimization methods [LO11, BJL15] which cover a wide range of statistical
applications. Let us emphasize that the aforementioned estimators do not require any assumptions
on the “shape” of the distribution, such as unimodality or elliptical symmetry.
Generalizations of univariate results to the case of random vectors and random matrices are
not straightforward since element-wise deviation inequalities do not always translate into desired
bounds. In some cases, element-wise bounds yield inequalities for the “wrong” norm: for example,
estimating each entry of the covariance matrix results in a deviation inequality for the Frobenius
norm, while we are frequently interested in the bounds for the operator norm that can be much
smaller. An approach which often yields “dimension - free” bounds was proposed in [HS16] and
[Min15] (using generalizations of the median in higher dimensions); however, to the best of our
knowledge, results of these papers are still not sufficient to obtain deviation guarantees in the operator norm that we are mainly interested in. Under more restrictive assumptions on the sequence of
random matrices Y1 , . . . , Yn (such as kYj k ≤ M almost surely for some fixed M P
> 0, j = 1, . . . , n,
1
where k · k stands for the operator norm), behavior of the sample mean Ȳ = n nj=1 Yj has been
analyzed with the help of matrix concentration inequalities [AW02, Oli09, Tro12a].
A closely related covariance matrix estimation problem has been extensively studied in the
1
For the classical treatment of robust estimators based on the notion of a breakdown point, we refer the reader
to [HR09].

S. Minsker/Estimators of the mean of a random matrix

3

past decades. A comprehensive review is beyond the scope of this introduction, so we will just
mention few classical results and more recent work related to the current line of research. Statistical properties of the sample covariance matrix for Gaussian and sub-Gaussian observations have
been investigated in detail, see [KL17, KL16, Ver10, CZZ10, CRZ16] and references therein; under
weaker moment assumptions, sample covariance estimator has been studied in [SV13]. Some popular robust estimators of scatter are discussed in [HRVA08], including the Minimum Covariance
Determinant (MCD) estimator and the Minimum Volume Ellipsoid estimator (MVE). However,
rigorous results for these estimators are available only for elliptically symmetric distributions; see
[BDJ93] for results on MCD and [Dav92] for results on MVE. Popular Maronna’s [Mar76] and
Tyler’s [Tyl87, ZCS16] M-estimators of scatter also admit theoretical guarantees for the family
of elliptically symmetric distributions, but we are unaware of any results extending beyond this
case.
Recent papers of O. Catoni [Cat16] and I. Guilini [Giu15], Fan et al. [FWZ16b] are closest in
spirit to our work. For instance, in [Cat16] author constructs a robust estimator of the Gram
matrix of a random vector Z ∈ Rd (as well as its covariance matrix) via estimating the quadratic
form E hZ, ui2 uniformly over kuk2 = 1, and obtains error bounds for the operator norm. The
latter (univariate) estimators for the quadratic form are based on the fruitful ideas originating in
[Cat12]. However, results of these works can not be straightforwardly extended beyond covariance
estimation, and are obtained under more stringent (compared to the present paper) assumptions
on the underlying distribution (such as a known upper bound on the kurtosis of hZ, ui2 for any u
of norm 1). In [FWZ16b], authors obtain error bounds for norms other than the operator norm
which is the main focus of the present paper.
Finally, let us mention that the problem of robust matrix recovery (that is discussed as an
example below) has also received attention recently: for instance, the work [CLMW11, KLT14]
investigates robust matrix completion under the “low rank + sparse” model. In [FWZ16a], authors
study low-rank matrix recovery under the assumption that the additive noise has only (2 + ε)
moments, and obtain strong results via truncation argument. We propose a different approach
based on general techniques developed in this paper and achieve similar results for the matrix
completion problem while requiring only the finite variance of the noise.
1.2. Organization of the paper
Section 2 contains definitions, notation and background material. Our main results are introduced
in section 3. After presenting core results, we discuss applications to covariance estimation and
low-rank matrix completion in section 4, and illustrate the role of various quantities involved
in the general bounds through these examples. Sections 5 and 6 discuss adaptation to unknown
parameters that appear in our construction, and contain longer proofs.
Appendix contains proofs omitted in the main text, as wells as several technical lemmas.
2. Preliminaries
In this section, we introduce main notation and recall several useful facts from linear algebra,
matrix analysis and probability theory that we rely on in the subsequent exposition.

S. Minsker/Estimators of the mean of a random matrix

4

2.1. Definitions and notation
Given A ∈ Cd1 ×d2 , let A∗ ∈ Cd2 ×d1 be the Hermitian adjoint of A. If A is self-adjoint, we will
write λmax (A) and λmin (A) for the largest and smallest eigenvalues of A. Next, we will introduce
the matrix norms used in the paper.
p
Everywhere below, k · k stands for the operator norm kAk := λmax (A∗ A). If d1 = d2 = d,
d1 ×d2 , the nuclear norm k · k is defined as
we denote √
by tr A the trace
1
√ of A. Next, for A ∈ C
√
∗
∗
kAk1 = tr ( A A), where A A is a nonnegativep
definite matrix such that ( A∗ A)2 = A∗ A. The
Frobenius (or Hilbert-Schmidt) norm is kAkF = tr (A∗ A), and the associated inner product is
hA1 , A2 i = tr (A∗1 A2 ). Finally, set kAkmax := supi,j |ai,j |. For Y ∈ Cd , kY k2 stands for the usual
Euclidean norm of Y .
Given two self-adjoint matrices A and B, we will write A  B (or A  B) iff A − B is
nonnegative (or positive) definite.
Given a sequence Y1 , . . . , Yn of random matrices, Ej [ · ] will stand for the conditional expectation
E[ · |Y1 , . . . , Yj ].
Finally, for a, b ∈ R, set a ∨ b := max(a, b) and a ∧ b := min(a, b).
2.2. Tools from linear algebra
In this section, we collect several facts from linear algebra, matrix analysis and probability theory
that are frequently used in our arguments.
Definition 2.1. Given a real-valued function f defined on an interval T ⊆ R and a self-adjoint
A ∈ Cd×d with the eigenvalue decomposition A = U ΛU ∗ such that λj (A) ∈ T, j = 1, . . . , d, define
f (A) as f (A) = U f (Λ)U ∗ , where

 

λ1
f (λ1 )

 

..
..
f (Λ) = f 
 = 
.
.
.
λd

f (λd )

Additionally, we will often use the following facts:
Fact 2.1. Let A ∈ Cd×d be a self-adjoint matrix, and f1 , f2 be two real-valued functions such
that f1 (λj ) ≥ f2 (λj ) for j = 1, . . . , d. Then f1 (A)  f2 (A).
Fact 2.2. Let A, B ∈ Cd×d be two self-adjoint matrices such that A  B. Then λj (A) ≥
λj (B), j = 1, . . . , d, where λj (·) stands for the j-th largest eigenvalue. Moreover, tr eA ≥ tr eB .
Fact 2.3. Matrix logarithm is operator monotone: if A  0, B  0 and A  B, then log(A) 
log(B).
Proof. See [Bha97].
Fact 2.4. Let A ∈ Cd×d be a self-adjoint matrix. Then I + A +

A2
2

 0. Moreover,





A2
A2
− log I + A +
 log I − A +
.
2
2

S. Minsker/Estimators of the mean of a random matrix

5

Proof. In view of the definition of a matrix function, the first claim follows from scalar inequality
1 + t + t2 /2 > 0 for t ∈ R. Similarly, the second relation follows from the inequality − log(1 + t +
t2 /2) ≤ log(1 − t + t2 /2) for t ∈ R.
Fact 2.5 (Lieb’s concavity theorem). Given a fixed self-adjoint matrix H, the function
A 7→ tr exp (H + log(A))
is concave on the cone of positive definite matrices.
Proof. See [Lie73] and [Tro15].

2

Fact 2.6. Let f : R 7→ R be a convex function. Then A 7→ tr f (A) is convex on the set of
self-adjoint matrices. In particular, for any self-adjoint matrices A, B,


A+B
1
1
tr f
≤ tr f (A) + tr f (B).
2
2
2
Proof. This is a consequence of Peierls inequality, see Theorem 2.9 in [Car10] and the comments
following it.
Finally, we introduce the Hermitian dilation which allows to reduce many problems involving
general rectangular matrices to the case of Hermitian operators. Given the rectangular matrix
A ∈ Cd1 ×d2 , the Hermitian dilation H : Cd1 ×d2 7→ C(d1 +d2 )×(d1 +d2 ) is defined as


0 A
H(A) =
.
(2.1)
A∗ 0


AA∗
0
, it is easy to see that kH(A)k = kAk. Another tool useful in dealing
Since H(A)2 =
0
A∗ A
with rectangular matrices is the following lemma:
Lemma 2.1. Let S ∈ Cd1 ×d1 , T ∈ Cd2 ×d2 be self-adjoint matrices, and A ∈ Cd1 ×d2 . Then




S A
0 A
≥
.
A∗ T
A∗ 0
Proof. See section A.1 in the appendix.
3. Main results
Our construction has its roots in the technique proposed by O. Catoni [Cat12] for estimating
the univariate mean. Let us briefly recall the main ideas of Catoni’s approach. Assume that
ξ, ξ1 , . . . , ξn is a sequence of i.i.d. random variables such that Eξ = µ and Var(ξ) ≤ v 2 . Catoni’s
estimator is defined as follows: let ψ(x) : R 7→ R be a non-decreasing function such that for all
x ∈ R,
− log(1 − x + x2 /2) ≤ ψ(x) ≤ log(1 + x + x2 /2).
2

(3.1)

Let us mention that Lieb’s theorem is one of the key tools for proving matrix concentration inequalities, and
its power in this context was first demonstrated by J. Tropp [Tro12b].

S. Minsker/Estimators of the mean of a random matrix

6

See remark 1 below for examples of such functions. Given θ > 0, let µ̂θ be such that
n
X

ψ (θ(ξj − µ̂θ )) = 0

(3.2)

j=1

(clearly, µ̂θ always exists due to monotonicity). Set η = v

q

2t
n(1−2t/n)

and θ∗ =

q

2t
.
n(v 2 +η 2 )
−t
− 2e .

Assum-

ing that n > 2t, it is shown in [Cat12] that |µ̂θ∗ − µ| ≤ η with probability ≥ 1
We proceed by presenting a multivariate extension of the estimator µ̂θ . We will first formulate
main results for the self-adjoint matrices, and will later deduce the general case of rectangular
matrices as a corollary. Let Y1 , . . . , Yn ∈ Cd×d be a sequence of independent self-adjoint random
n
P
matrices such that σn2 :=
EYj2 < ∞. Let Ψ be such that Ψ0 (x) = ψ(x) for all x ∈ R, and set
j=1


Tbθ∗ = argminS∈Cd×d ,S=S ∗ tr

n
X


Ψ (θ(Yj − S)) ,

(3.3)

j=1

where θ > 0 is an appropriate constant. It follows from Fact 2.6 that Tbθ∗ exists, moreover, it is
unique if ψ(x) is strictly increasing. It is also not hard to see that (3.3) is equivalent to
n
X



ψ θ(Yj − Tbθ∗ ) = 0d×d .

(3.4)

j=1

P
Indeed, if Fψ (S) := tr nj=1 Ψ (θ(Yj − S)), then (3.4) simply states that the gradient of Fψ evaluated at Tbθ∗ is equal to zero; see Lemma A.1 in the appendix for more details.
To understand the properties of the estimator defined via (3.3) and (3.4), we will first consider
(0)
another estimator Tbθ that shares many important properties with Tbθ∗ but is easier to analyze.
(0)
The “preliminary estimator” Tbθ is constructed as follows: given θ > 0 and a function ψ
satisfying (3.1), set Xj := ψ (θYj ) , j = 1, . . . , n and
n
1 X
(0)
Tbθ :=
Xj .
nθ

(3.5)

j=1

(0)
In other words, Tbθ is an average of “ψ-truncated” observations. Since Xj ' θYj for small θ and
P
(0)
a smooth function ψ, we expect that Tbθ is close to n1 nj=1 EYj . In the following sections, we will
make this intuition more precise. In particular, we will establish the following (so far informally
stated) results:

Theorem.
1. Assume that the observations Y1 , . . . , Yn are i.i.d. copies of Y ∈ Cd×d and the
parameter θ is chosen properly. Then
r !
 
t
t
(0)
b
Pr
Tθ − EY ≥ σ
≤ 2d exp −
,
n
2
where σ 2 := σn2 /n = EY 2 .

S. Minsker/Estimators of the mean of a random matrix

7

2. Assume that n is large enough and θ is chosen properly. Then the estimator Tbθ∗ defined via
(3.4) satisfies the inequality
r !
 
t
t
∗
b
≤ C2 d exp −
,
Pr
Tθ − EY ≥ C1 σ0
n
2
where C1 , C2 > 0 are absolute constants and σ02 := E(Y − EY )2 .
Note that the “variance term” EY 2 appearing in the first part of the bound above is akin
to the second moment, while in the second bound it is replaced by σ02 = E(Y − EY )2 ; presence
(0)
of the term EY 2 can be explained by the fact that the estimator Tbθ is obtained via biasproducing truncation. We remark that in some applications, such as matrix completion discussed
(0)
in section 4, even the estimator Tbθ with “suboptimal” variance term suffices to obtain good
bounds.
Remark 1. Most of our results do not depend on the concrete choice of the function ψ. One
possibility is
 

log 1 + x + x2 ,
x ≥ 0,
2


(3.6)
ψ1 (x) =
2
− log 1 − x + x , x < 0.
2
Another example is
(
2
x − x2 ,
ψ2 (x) =
1/2,

x ∈ [−1, 1],
|x| > 1.

(3.7)

Since the latter function is bounded, it can provide additional advantages (such as robustness) in
applications. However, note that ψ2 (x) does not satisfy (3.1); instead, it satisfies a slightly weaker
inequality


− log 1 − x + x2 ≤ ψ2 (x) ≤ log 1 + x + x2 ,
hence all subsequent results hold for ψ2 as well, albeit with slightly worse constant factors. We
also note that both ψ1 and ψ2 are operator Lipschitz functions; see Lemma A.3 for details.
3.1. Bounds for the moment generating function
(0)
1 Pn
In this section, we will establish deviation inequalities for the estimator Tbθ = nθ
j=1 ψ(θXj ).
The lemma below is the cornerstone of our results. As before, given θ > 0, let Xj = ψ(θYj ).

Lemma 3.1. The following inequalities hold:




n
n
2
X
X
θ
Etr exp 
(Xj − θEYj ) ≤ tr exp 
EYj2  ,
2
j=1
j=1




n
n
2
X
X
θ
(θEYj − Xj ) ≤ tr exp 
EYj2  .
Etr exp 
2
j=1

j=1

(3.8)

(3.9)

S. Minsker/Estimators of the mean of a random matrix

8

Proof. Note that

E tr exp 

n
X


(Xj − θEYj ) =

j=1



n−1
X

EEn−1 tr exp 





(Xj − θEYj ) − θEYn  + ψ (θYn ) ≤

j=1



n−1
X

EEn−1 tr exp 






(Xj − θEYj ) − θEYn  + log I + θYn + θ2 Yn2 /2  ≤

j=1


E tr exp 

n−1
X


(Xj − θEYj ) + log I + θEYn + θ2 EYn2 /2 − θEYn  ,


j=1



2
where the first inequality follows from the semidefinite relation ψ(θYn )  log I + θYn + θ2 Yn2
and
Fact 2.2, and the second inequality follows from Lieb’s concavity theorem (Fact 2.5) with H =
Pn−1
j=1 (Xj − θEYj ) − θEYn and Jensen’s inequality for conditional expectation. We also note that

I +θEYn +θ2 EYn2 /2  0 since I +θYn +θ2 Yn2 /2  0 almost surely, hence log I + θEYn + θ2 EYn2 /2
is well-defined. Repeating the steps for Xn−1 , . . . , X1 , we obtain the inequality




n
n
X
X


E tr exp 
(Xj − θEYj ) ≤ tr exp 
log I + θEYj + θ2 EYj2 /2 − θEYj 
(3.10)
j=1

j=1

It remains to note that by Fact 2.1 and the inequality log(1 + x) ≤ x (that holds ∀ x > −1), for
all j = 1, . . . , n

θ2
log I + θEYj + θ2 EYj2 /2  θEYj + EYj2 ,
2


2
θ
or log I + θEYj + θ2 EYj2 /2 − θEYj  2 EYj2 . The first inequality (3.8) now follows from (3.10)
and Fact 2.2.
To
 establish the second inequality of the lemma, we use the relation −Xj = −ψ(θYj ) 

log I − θYj +
that

θ2 2
2 Yj

(which follows from (3.1) and Fact 2.1) together with the Fact 2.2 to deduce



n
X
E tr exp 
(θEYj − Xj ) ≤
j=1


E tr exp 

n
X



log I + θ(−Yj ) + θ2 Yj2 /2 − θE(−Yj )  ,


j=1

and apply inequality (3.8) to the sequence −Y1 , . . . , −Yn with

Xj = log I + θ(−Yj ) + θ2 (−Yj )2 /2 , j = 1, . . . , n.

S. Minsker/Estimators of the mean of a random matrix

9

We are ready to state and prove the main result of this section.
Theorem 3.1. Let Y1 , . . . , Yn ∈ Cd×d be a sequence of independent self-adjoint random matrices,
Pn
2
and σn2 ≥
j=1 EYj . Then for all θ > 0





n 
X
√
√
1
θ2 σn2


Pr
.
ψ (θYj ) − EYj
≥ t n ≤ 2d exp −θt n +
θ
2
j=1

In particular, setting θ =

√
t n
2 ,
σn



2
we get the “sub-Gaussian” tail bound 2d exp − 2σt2 /n , for a given
n

√

t > 0. Alternatively, setting θ = σ2n (independent of t), we obtain sub-exponential concentration
n


2t−1
with tail 2d exp − 2σ
for all t > 1/2.
2 /n
n

Remark 2. In the important special case when Yj , j = 1, . . . , n are i.i.d. copies of Y , we will
2
2
often use the following equivalent
p s 1 form of of the bound: assume that σ ≥ kEY k, then replacing
√
t by σ s and setting θ := n σ implies that
r 

s
(0)
≤ 2d exp (−s/2) ,
(3.11)
Pr Tbθ − EY ≥ σ
n
(0)
where Tbθ was defined in (3.5).

Proof. As before, set Xj := ψ (θYj ) , j = 1, . . . , n. Then


!
n
X
1
Pr λmax 
(Xj − θEYj ) ≥ s
θ
j=1





n
X
= Pr exp λmax 
(Xj − θEYj ) ≥ eθs 
j=1





n
n
2
X
X
θ
≤ e−θs E tr exp 
(Xj − θEYj ) ≤ e−θs tr exp 
EYj2 
2
j=1
j=1


n
θ2 X
EYj2  ,
≤ d exp −θs +
2
j=1

where we used Chebyshev’s inequality, the fact that eλmax (A) = λmax (eA ) and the inequality
λmax (eA ) ≤ tr eA on the second step, the first inequality of Lemma 3.1 on the third step, and
the bound tr eA ≤ d ekAk on the last step (here and below, A ∈ Cd×d is an arbitrary self-adjoint
matrix). Similarly, since −λmin (A) = λmax (−A), we have


!
n
X
1
Pr λmin 
(Xj − θEYj ) ≤ −s
θ
j=1

S. Minsker/Estimators of the mean of a random matrix







n
1X

10



(θEYj − Xj ) ≥ s
θ
j=1




n
n
2
X
X
θ
EYj2 
≤ e−θs Etr exp 
(θEYj − Xj ) ≤ e−θs tr exp 
2
j=1
j=1


n
θ2 X
≤ d exp −θs +
EYj2  ,
2
= Pr λmax 

j=1

√
where we used the second inequality of Lemma 3.1 instead. The result follows by taking s := t n
since for a self-adjoint matrix A, kAk = max ( λmax (A), −λmin (A) ).
The main weakness of the estimator Tbθ0 discussed above is the fact that the “variance term”
Pn
2 appearing in the bound is akin to the second moment (the price we pay for applyj=1 EYj
Pn
2
ing bias-producing truncation) while we would like to replace it by
j=1 E(Yj − EYj ) . This
problem will be addressed in detail in section 6. In particular, we will
q show the following: assume
that Y1 , . . . , Yn are i.i.d. copies of Y , σ02 ≥ E(Y − EY )2 , θ0 =

2t 1
n σ0 ,

and n is large enough
(n & d2 ). Then, with exponentially high probability with respect to s, the solution Tbθ∗0 of equation
p
(3.4) satisfies Tbθ∗0 − EY ≤ Cσ0 ns for an absolute constant C > 0. Another problem is the fact
Pn
2 (or its tight upper bound) a priori to choose
that one needs to know the value of
j=1 EYj
the “optimal” value of parameter θ. This issue and its resolution based on adaptive estimators is
discussed in section 5. We conclude this discussion with few additional comments.

Remark 3.
1. Sub-Gaussian guarantees provided by Theorem 3.1 hold for a given confidence
parameter t > 0 that has to be fixed a priori: in particular, the optimal value of θ depends it.
However, as it was noted in [DLLO15], this is sufficient to construct (via Lepski’s method [Lep92])
estimators that admit sub-Gaussian tails uniformly over t in a certain range. We discuss the details
in section 6 of the supplementary material.
2. Let Y1 , . . . , Yn ∈ Cd×d be i.i.d. copies of Y , and σ02 = kE(Y − EY )2 k. It is interesting
toP
compare our estimator (in particular, bound (3.11)) to the guarantees for the sample mean
n
1
j=1 Yj . Under an additional restrictive boundedness assumption requiring that kY k ≤ M
n
almost surely, the noncommutative
q Bernstein’s inequality (see Theorem 1.4 in [Tro12b]) implies
1 Pn
that n j=1 Yj − EY ≤ 2σ0 nt ∨ 43 Mn t with probability ≥ 1 − 2de−t/2 . Hence, even under
additional strong assumptions our technique allows to obtain guarantees that compare favorably
to the sample mean.
3.2. Bounds depending on the effective dimension
The bound obtained in Theorem 3.1 explicitly depends on the dimension d of random matrices.
Example is subsection 3.2.1 below shows that the dimensional factor in the right-hand side of the

S. Minsker/Estimators of the mean of a random matrix

11

inequality is unavoidable in general. However, it is possible to prove a similar inequality which
only includes the “effective dimension” defined as
P

n
2
EY
tr
j
j=1
d¯ :=
(3.12)
Pn
2
j=1 EYj
which can be much smaller than d if
following result holds:

Pn

2
j=1 EYj

has many eigenvalues that are close to 0. The

Theorem 3.2. Let Y1 , . . . , Yn ∈ Cd×d be a sequence of independent self-adjoint random matrices,
Pn
2
and σn2 ≥
j=1 EYj . Then

Pr 

n 
X
1
j=1

θ


ψ (θYj ) − EYj


¯

≥ t n ≤ 2d 1 +

Remark 4. As before, we can set θ =

Pr 

n 
X
1

θ

j=1

√
t n
2
σn


ψ (θYj ) − EYj



√

1
√
θt n





√
θ2 σn2
exp −θt n +
.
2

to get






2 /n
2
σ
t
n
exp − 2
.
≥ t n ≤ 2d¯ 1 + 2
t
2σn /n
√

p
For the values of t ≥ σn2 /n (when the bound becomes useful), it further simplifies to





n 
X
√
t2
1
ψ (θYj ) − EYj
≥ t n ≤ 4d¯exp − 2
.
Pr 
θ
2σn /n
j=1

√

For the “sub-exponential regime” with θ =

Pr 

n 
X
1
j=1

θ

n
2 ,
σn


ψ (θYj ) − EYj

we get that for all t ≥
√

1
2

∨ σn2 /n simultaneously,





2t − 1
¯

≥ t n ≤ 4d exp − 2
.
2σn /n

Proof. The argument is similar in spirit to the proof of Theorem 3.1. Details are included in
section B.2.
3.2.1. Dimensional factor in Theorem 3.1
Example below shows that the dimensional factor in Theorem 3.1 is unavoidable in general.
Assume that ψ(x) = ψ1 (x) as defined in (3.6), θ = 1, n = d, and let Yj , j ≤ d be independent and
such that ψ1 (Yj ) = γj ej eTj , where γj , j ≤ d are i.i.d. random variables with density p(x) = e−2|x| ,

S. Minsker/Estimators of the mean of a random matrix

12

and {e1 , . . . , ed } is the standard Euclidean basis. Recalling that Yj = ψ1−1 (γj )ej eTj , it is easy to

Pd
2 = E ψ −1 (γ ) 2 < ∞. Theorem 3.1 implies that
check that EYj = 0d×d , and that
1
1
j=1 EYj

Pr 

d
X


γj ej eTj ≥ s ≤ f (d)e−s

j=1

Pd

T
j=1 γj ej ej

with f (d) ≤ Cd for some absolute constant C. Since

= max (|γ1 |, . . . , |γd |), it fol

T ≥ 1 − τ log d ≥
γ
e
e
lows from Lemma 7.2 of the supplementary material that Pr
j=1 j j j
2
c(τ ) for any 0 < τ < 1/2 and some constant c(τ ) > 0. This shows that the dimensional factor
f (d) can not grow slower than d1/2−τ for any τ > 0.
 P
d

3.3. Bounds for arbitrary rectangular matrices
In this section, we will deduce results for arbitrary matrices from the bounds for self-adjoint
operators. Let Y1 , . . . , Yn ∈ Cd1 ×d2 be independent, and assume that


n
n
X
X
σn2 ≥ max 
EYj Yj∗ ,
EYj∗ Yj  .
j=1

j=1

Given θ > 0, set Xj := ψ (θH(Yj )) (where H(·) is the self-adjoint dilation, see equation (2.1)) and
define T̂ ∈ C(d1 +d2 )×(d1 +d2 ) as
Tb := Tb(θ) =

n
X
1
j=1

θ

Xj .



T̂11 T̂12
Let T̂11 ∈
T̂22 ∈
T̂12 ∈
. Since Tb is “close”
∗
T̂12
T̂22
P
P
to nj=1 H (EYj ) for the proper choice of θ, it is natural to expect that T̂12 is close to nj=1 EYj .
Cd1 ×d1 ,

Cd2 ×d2 ,

Cd1 ×d2

be such that Tb =

Corollary 3.1. Under the assumptions stated above,




n
X
√
√
θ2 σn2


T̂12 −
EYj ≥ t n ≤ 2(d1 + d2 ) exp −θt n +
Pr
2
j=1

and


Pr  T̂12 −

n
X



√
≥ t n ≤ 2d¯ 1 +

EYj

j=1

where d¯ = 2

tr (

∗
j=1 EYj Yj
Pn
∗
j=1
j=1 EYj Yj ∨

Pn

k

Pn

k k

)
EYj∗ Yj k

.

1
√
θt n





√
θ2 σn2
exp −θt n +
,
2

S. Minsker/Estimators of the mean of a random matrix

13

Proof. Note that
n
X


E H(Yj )2 = max 

j=1

n
X

n
X

EYj Yj∗ ,

j=1


EYj∗ Yj  ≤ σn2 .

j=1

Theorem 3.1 applied to self-adjoint random matrices H(Yj ) ∈ C(d1 +d2 )×(d1 +d2 ) , j = 1, . . . ,n
2 2
P
√
√
implies that Tb − n H(EYj ) ≤ t n with probability ≥ 1 − 2(d1 + d2 ) exp −θt n + θ σn .
j=1

2

It remains to apply Lemma 2.1:
Tb −

n
X

H(EYj ) =

j=1

≥

∗ −
T̂12

∗ −
T̂12

T̂
P11
n

T̂12 −

0
Pn

T̂12 −

∗
j=1 EYj

Pn

j=1 EYj

!

T̂22

∗
j=1 EYj

Pn

j=1 EYj

0

!
= T̂12 −

n
X

EYj ,

j=1

and the first inequality follows. To obtain the second inequality, it is enough to use Theorem 3.2
instead of Theorem 3.1 and note that








n
n
n
n
X
X
X
X
tr 
EH(Yj )2  = tr 
EYj Yj∗  + tr 
EYj∗ Yj  = 2tr 
EYj∗ Yj 
j=1

j=1

j=1

j=1



since for any 1 ≤ j ≤ n, tr EYj Yj∗ = Etr (Yj Yj∗ ) = Etr (Yj∗ Yj ).
T = Σ and Y , . . . , Y are
In a particular case
Y ∈ Rd is a random vector
1
n
 when
 such that EY Y P
Pn
P
n
n
∗ ,
∗Y
∗Y
its i.i.d. copies, max
EY
Y
EY
=
n
tr
Σ
and
tr
EY
=
n
tr Σ,
j j
j j
j j
j=1
j=1
j=1
√
√
hence d¯p
= 2 and the estimator T̂12 admits the following bound: if we replace t by s tr Σ and
set θ = ns √ 1 in the second bound of Corollary 3.1, then
tr Σ
r !
√
T̂12
s
Pr
− EY ≥ tr Σ
≤ 4 (1 + 1/s) e−s/2 .
n
n
2

3.4. Bounds under weaker moment assumptions
In this section, we discuss the mean estimation problem under weaker moment conditions. Namely,
assume that Y1 , . . . , Yn are independent self-adjoint random matrices such that kE|Yj |α k < ∞ for
some α ∈ (1, 2] and all 1 ≤ j ≤ n. Let ψα satisfy
− log(1 − x + cα |x|α ) ≤ ψα (x) ≤ log(1 + x + cα |x|α )
q
α−1
for all x ∈ R, where cα = α ∨ 2−α
α . The fact that such ψα exists follows from Lemma A.2 in
the appendix. For example, one can take ψα (x) = log(1 + x + cα |x|α ). The following result holds:

S. Minsker/Estimators of the mean of a random matrix

14

Pn
α . Then for any positive t and θ,
Theorem 3.3. Assume that vnα ≥
j=1 E|Yj |



n 
X
1
≥ t ≤ 2d exp (−θt + cα θα vnα ) .
ψα (θYj ) − EYj
Pr 
θ
j=1

Proof. The argument repeats the steps of Lemma 3.1 and Theorem 3.1, the only difference being
that application of Fact 2.4 is replaced by Lemma A.2.
Remark 5. In the special case when Y1 , . . . , Yn are i.i.d. copies of Y with v = kE|Y |α k1/α , setting

1/(α−1) 
α−1
s 1/α 1
t = vn1/α s α and θ = αc1α
n
v gives the inequality



1/(α−1) !
n
α−1


X
1
1
s
α
−
1
α
 ≤ 2d exp −
ψα (θYj ) − EY ≥ v
s .
Pr 
nθ
n
α
αcα
j=1

Note that for α = 2, we recover (3.11).
Before we proceed with discussion or further improvements and adaptation issues, let us demonstrate applications of developed techniques to popular problems in statistics and highlight the
advantages over existing results.
4. Examples
We present two examples which highlight the potential improvements obtained via our technique
in popular scenarios: estimation of the covariance matrix in Frobenius and operator norms, and
low-rank matrix completion problem.
4.1. Estimation of the covariance matrix in operator norm


Let Z ∈ Rd be a random vector with EZ = µ, EkZ − µk42 < ∞, Σ = E (Z − µ)(Z − µ)T , and
let Z1 , . . . , Z2n be i.i.d. copies of Z. Let us first assume that µ = 0, and define
2n
X

e 2n (θ) = 1
Σ
ψ θZj ZjT ,
2nθ
j=1

where ψ(·) satisfies (3.1). Let σ 2 ≥ EkZk22 ZZ T
from Theorem 3.1 that with probability ≥ 1 −

and θ̃ =

t 1
n σ.

It is straightforward to deduce

2de−t ,
r

e 2n (θ̃) − Σ ≤ σ
Σ
Remark 6.

q

t
.
n

1. Note that for any matrix X = λU U T of rank 1 (where kU k2 = 1),
ψ(X) = ψ(λ)U U T (since ψ(0) = 0),

T
P
2 Zj Zj
hence Σ̃2n (θ̃) = 2n1 θ̃ 2n
j=1 ψ(θ̃kZj k2 ) kZj k22 . In particular, this expression is easy to evaluate numerically; in general, computation of the estimator (3.5) requires n singular value decompositions.

S. Minsker/Estimators of the mean of a random matrix

15

(Σ)
[Ver10]; clearly,
2. Parameter σ is closely related to the effective rank defined as r(Σ) = trkΣk
p
it always true that r(Σ) ≤ d. The quantity r(Σ)kΣk has been shown to control the expected
error of the sample covariance estimator in the Gaussian setting [KL17]. Under the additional
assumption that the kurtosis of the coordinates Z (j) := hZ, ej i of Z is uniformly bounded by K, it
is possible to show that (see Lemma 2.3 in [MW17]) that σ 2 ≤ K r(Σ) kΣk2 . On the other hand,
fluctuations of the error around its expected value √
in the Gaussian case [KL17] are controlled by
the “weak variance” supv∈Rd :kvk2 =1 E1/2 hZ, vi4 ≤ KkΣk, while in our bounds fluctuations are
controlled by the “strong variance” σ 2 ; this fact leaves room for improvement of our construction
and proof techiques.

Of course, the initial assumption that µ is known is often unrealistic, hence we modify the
estimator as follows. Given θ > 0, set
1
(Z2j−1 − Z2j ) (Z2j−1 − Z2j )T ,
2
n
1 X
b
Σ2n (θ) =
ψ(θYj ).
nθ
Yj =

j=1

q
2
E (Z − µ)(Z − µ)T + tr (Σ)Σ + 2Σ2 , and θ̂ = nt σ̂1 . Our covariance estimator is
b 2n := Σ
b 2n (θ̂). The following result can be deduced from Theorem 3.1:
then defined as Σ

Let σ̂ 2 ≥

1
2

Corollary 4.1. With probability ≥ 1 − 2d e−t ,
b 2n − Σ ≤
Σ

√

r
2σ̂

t
.
n

Before presenting the proof, let us make several additional remarks.
Remark 7.
1. It is not hard to show that (see Corollary 7.1 of the supplementary material)
2
that E (Z − µ)(Z − µ)T
≥ tr (Σ) kΣk, hence it is enough to choose σ̂ 2 ≥ kΣk2 + σ02 =
2
kΣk2 + E (Z − µ)(Z − µ)T
. In view of remark 6, this expression can be further simplified
under the bounded kurtosis assumption, and one can choose σ̂ 2 ≥ kΣk2 (1 + Kr(Σ)), where K is
the uniform bound on the kurtosis of the coordinates of Z, and r(Σ) is the effective rank.
b 2n (θ) essentially halves the effective sample size. While the loss of a con2. Construction of Σ
stant factor can be deemed insignificant in non-asymptotic theoretical bounds, it is undesirable in
applications. A more natural version of the estimator based on a sample of size 2n is the U-statistic

X 1 θ
1
T
Σ̄2n (θ) = 2n
ψ
(Zi − Zj )(Zi − Zj )
.
θ
2
2
1≤i<j≤2n

Another possibility to avoid “halving” the sample size is to center the data using a robust estimator
of location, such as the spatial median or the median-of-means estimator [JLO+ 17, LM17, Min15].
Analysis of the estimators of these types is not covered in the present paper, and requires a slightly
different set of technical tools to deal with dependent summands; see [MW17] for results in this
direction.

S. Minsker/Estimators of the mean of a random matrix

16

Proof of Corollary 4.1. Note that for all j = 1, . . . , n, EYj = Σ. Since Y1 , . . . , Yn are i.i.d. random
matrices, Theorem 3.1 applies (see remark 2), giving that
r !
2t
Pr
≤ 2de−t ,
Σ̂(θ̂) − Σ ≥ σ̂
n
where σ̂ 2 ≥ EY12 . It is easy to check that
EY12 =

2
1
E (Z − µ)(Z − µ)T + tr (Σ)Σ + 2Σ2 ,
2

and result follows.
4.2. Estimation of the covariance matrix in Frobenius norm
Next, we present an estimator which achieves strong deviation guarantees in the Frobenius norm.
Estimation of the covariance matrix with respect to this norm has been previously investigated
in the literature, for instance, see [LF09], [CZZ10] and references therein; Frobenius norm is a
natural choice when one wants to understand the effect of the rank of an unknown covariance
matrix on the estimation error [Lou14]. Let Ŝ2n be the sample covariance estimator based on
Z1 , . . . , Z2n :
X (Zi − Zj )(Zi − Zj )T
1
Ŝ2n = 2n
.
2
2
1≤i<j≤2n

The following “soft thresholding” estimator has been studied in [Lou14]; here, τ > 0 is a fixed
threshold parameter:


2
τ
Ŝ2n = argminA∈Rd×d A − Ŝ2n + τ kAk1 .
(4.1)
F

b 2n , and consider
We propose to replace the sample covariance Ŝ2n by Σ


2
τ
b
b
Σ2n = argminA∈Rd×d A − Σ2n + τ kAk1 .
F

(4.2)

b τ can be written explicitly
It is not hard to see (e.g., see the proof of Theorem 1 in [Lou14]) that Σ
2n
as
d
 


X
τ
b
b 2n − τ /2, 0 vj (Σ
b 2n )vj (Σ
b 2n )T ,
Σ2n =
max λj Σ
j=1

b 2n ) and vj (Σ
b 2n ) are the eigenvalues and corresponding eigenvectors of Σ
b 2n . The folwhere λj (Σ
lowing result holds:
Theorem 4.1. For any
r
τ ≥ 4σ̂

t + log(2d)
2n

S. Minsker/Estimators of the mean of a random matrix

bτ − Σ
Σ
2n

2
F

"
≤

inf

A∈Rd×d

kA −

Σk2F

17

#
√
(1 + 2)2 2
+
τ rank(A) .
8

(4.3)

with probability ≥ 1 − e−t .
Result stated above mimics the (almost) optimal rates obtained in [Lou14] (in the situation
when no data is missing) under significantly weaker assumptions on the underlying distribution.
Proof of Theorem 4.1. The proof is based on the following lemma:
n
o
b 2n − Σ .
Lemma 4.1. Inequality (4.3) holds on the event E = τ ≥ 2 Σ
To verify this statement, it is enough to repeat the steps of the proof of Theorem 1 in [Lou14],
bτ
replacing each occurrence of the sample covariance Ŝ2n by its robust counterpart
q Σ2n .
Result then follows from corollary 4.1 that Pr(E) ≥ 1 − e−t whenever τ ≥ 4σ̂

t+log(2d)
.
2n

4.3. Matrix completion
Let A0 ∈ Rd1 ×d2 be an unknown matrix, and assume that we observe a random subset of its
entries contaminated by noise. The goal is to estimate A0 from a small number of such noisy
measurements under an additional assumption that A0 is likely to be of low rank (or can be well
approximated by a low rank matrix). More specifically, let

X = ej (d1 )eTk (d2 ), 1 ≤ j ≤ d1 , 1 ≤ k ≤ d2 ,
where ej (d1 ) and ek (d2 ) are the elements of the canonical bases of Rd1 and Rd2 respectively. Let
X have uniform distribution Π := Unif(X ) on X , and assume that the noisy linear measurement
Y has the form
Y = tr (X T A0 ) + ξ,
where E(ξ|X) = 0. Finally, assume that (X1 , Y1 ), . . . , (Xn , Yn ) are i.i.d. copies of (X, Y ).
It is easy to check that E(Y X) = d11d2 A0 , hence the natural unbiased estimator of A0 is
n
X
b = d1 d2
Yj Xj .
A
n
j=1

To incorporate the structural (low-rank) assumption on A0 , the following estimator has been
considered in the literature: let τ > 0, and define


1
τ
2
b = argmin
b + τ kAk1
kA − Ak
A
F
A∈Rd1 ×d2
d1 d2


* n
+
X
1
2
= argminA∈Rd1 ×d2 
kAk2F −
Yj Xj , A + τ kAk1  .
d1 d2
n
j=1

S. Minsker/Estimators of the mean of a random matrix

18

bs ∈ R(d1 +d2 )×(d1 +d2 ) of A
b instead, defined as
Note that one can use the symmetric version A
n
X
bs = d1 d2
A
Yj H(Xj ),
n
j=1

bs = H(A0 ), and consider the equivalent convex minimization problem
so that EA


1
2
τ
b
b
kH(A) − H(As )kF + 2τ kAk1
A = argminA∈Rd1 ×d2
d1 d2


* n
+
X
2
1
kH(A)k2F −
= argminA∈Rd1 ×d2 
Yj H(Xj ), H(A) + 2τ kAk1  .
d1 d2
n
j=1

However, strong theoretical guarantees for this estimator exist only when the “noise term” ξj is
bs with a
either bounded with probability 1, or has sub-exponential tails. We propose to replace A
robust estimator
n
d1 d2 X
b
ψ (θYj H(Xj )) ,
R=
nθ
j=1

where ψ(·) satisfies (3.1) and
θ := θ(t, n, A0 ) =

1
p
kA0 kmax ∨ Var(ξ)

r

(t + log(2(d1 + d2 )))(d1 ∧ d2 )
.
n

The reasoning behind this choice of θ is explained below. Consider




1
2 b
τ
2
b
R = argminA∈Rd1 ×d2
R, H(A) + 2τ kAk1 .
kH(A)kF −
d1 d2
d1 d2
Finally, set
b − E (Y H(X)) .
M =R
The following result holds:
Theorem 4.2. Assume that ξj is independent of Xj , j = 1, . . . , n, and that Var(ξ) < ∞. For
any
s

 t + log(2(d + d ))
p
1
2
τ ≥ 4 kA0 kmax ∨ Var(ξ)
,
n(d1 ∧ d2 )


√ !2
2
1
2
1
1
+
b τ − A0 ≤
d1 d2 τ 2 rank(A) .
inf 
R
kA − A0 k2F +
d1 d2
2
F
A∈Rd1 ×d2 d1 d2
with probability ≥ 1 − e−t .
Note that we only assume that Var(ξ) < ∞, while in [FWZ16a], a similar result is obtained under
a slightly stronger assumption requiring that E|ξ|2+ε < ∞ for some ε > 0.

S. Minsker/Estimators of the mean of a random matrix

19

Proof. Define A ⊆ R(d1 +d2 )×(d1 +d2 ) to be the image of Rd1 ×d2 under H(·):
n
o
A = B ∈ R(d1 +d2 )×(d1 +d2 ) : B = H(A) for some A ∈ Rd1 ×d2 .
We begin with the following inequality:
Lemma 4.2. Assume that τ ≥ 2kM k. Then

2
1
bτ ) − H(A0 ) ≤ inf  1 kB − H(A0 )k2 +
H(R
F
B∈A d1 d2
d1 d2
F


√ !2
1+ 2
d1 d2 τ 2 rank(B) .
2

bτ , we see that
Proof. By the definition of R




2 b
1
2
τ
b
kBkF −
R, B + τ kBk1 .
H(R ) = argminB∈A
d1 d2
d1 d2
1 Pn
b by 1 A
b
If we replace d11d2 R
j=1 Yj H(Xj ), the result follows from Theorem 1 in [KLT11]
d1 d2 s = n
immediately. To obtain the current statement, it is enough to repeat the argument of Theorem 1
bs by 1 R.
b
in [KLT11], replacing each occurrence of the matrix d11d2 A
d1 d2
To complete the proof, we will estimate each side of the inequality of Lemma 4.2. First, it is
obvious from the definition of the Frobenius norm that
1
bτ ) − H(A0 )
H(R
d1 d2

2
F

=

2
b τ − A0
R
d1 d2

2
F

.

Next, since rank(H(A)) = 2 rank(A),


√ !2
1+ 2
1
inf 
kB − H(A0 )k2F +
d1 d2 τ 2 rank(B)
B∈A d1 d2
2


√ !2
1
1+ 2
= 2 inf 
kA − A0 k2F +
d1 d2 τ 2 rank(A) .
d
×d
d
d
2
1
2
A∈R
1 2

(4.4)

(4.5)

It remains to estimate the probability of the event E = {τ ≥ 2kM k}. Let



 
σ 2 := max E Y 2 XX T , E Y 2 X T X
.
Lemma 4.3. Assume that ξj is independent of Xj , j = 1, . . . , n. Then
2
.
d1 ∧ d2
h
i
2




Proof. Note that E Y 2 XX T = E ξ 2 XX T + E tr (X T A0 ) XX T . Moreover, |tr (X T A0 )| ≤
σ 2 ≤ Var(ξ) ∨ kA0 k2max

maxi,j |(A0 )i,j | = kA0 kmax , and EXX T =

1
d1



, hence



1
1
E Y 2 XX T ≤ Var(ξ) + kA0 k2max .
d1
d1

S. Minsker/Estimators of the mean of a random matrix

20

Similarly,


1
1
E Y 2 X T X ≤ Var(ξ) + kA0 k2max .
d2
d2
Applying Theorem 3.1 (see remark 2) with
r
2(t + log(2(d1 + d2 )))
θ=

n
=

r

1
kA0 kmax ∨

p

Var(ξ)

1
2
(Var(ξ) ∨ kA0 k2max ) d1 ∧d
2

1/2

(t + log(2(d1 + d2 )))(d1 ∧ d2 )
,
n

we see that

s



p
kM k ≤ 2 kA0 kmax ∨ Var(ξ)

t + log(2(d1 + d2 ))
n(d1 ∧ d2 )

with probability ≥ 1 − e−t . Final result now follows from the combination of this inequality with
(4.4), (4.5) and Lemma 4.2.

5. Optimal choice of θ and adaptation to the unknown second moment
To make results of Theorem 3.1 useful, one has to set the value for the parameter θ which in turn
Pn
2
depends on the (usually unknown) norm σn2 =
j=1 EYj . To address this problem, we develop
a simple adaptive solution based on Lepski’s method.
Lepski’s method [Lep92] is a powerful general technique that allows to adapt to the unknown
structure of the problem - for example, bandwidth selection in nonparametric estimation, or
unknown second moment in our case. Let Y1 , . . . , Yn ∈ Cd×d be independent self-adjoint random
Pn
2
matrices with σn2 =
j=1 EYj , and assume that σmin , σmax are such that
σn
σmin ≤ √ ≤ σmax .
n
√
Parameters σmin and σmax are “crude” preliminary bounds that can differ from σn / n by several
orders of magnitude. Let σj = σmin 2j and
J = {j ∈ Z : σmin ≤ σj < 2σmax }
be a set of cardinality |J | ≤ 1 + log2 (σmax /σmin ), and for each j ∈ J set θj = θ(j, t) =
Define
n
1 X
Tn,j =
ψ(θj Yi ),
nθj
i=1

q

2t 1
n σj .

S. Minsker/Estimators of the mean of a random matrix

where ψ(·) satisfies (3.1). Finally, set
(

21

r

j∗ := min j ∈ J : ∀k > j s.t. k ∈ J , kTn,k − Tn,j k ≤ 2σk

2t
n

)
(5.1)

and Tn∗ := Tn,j∗ .
Next result shows that adaptation is possible at the cost of an additional multiplicative constant
factor 6 in the deviation bound.
Theorem 5.1. The following inequality holds for any t > 0:
r !


√
2t
2σmax
∗
Pr kTn − EY k ≥ 6(σn / n)
≤ 2d log2
e−t .
n
σmin
n
o
σn
σn
Proof. Let j̄ = min j ∈ J : σj ≥ √
(hence σj̄ ≤ 2 √
). First, we will show that j∗ ≤ j̄ with
n
n
high probability. Indeed,

(
r )
[
2t 
Pr (j∗ > j̄) ≤ Pr 
Tn,k − Tn,j̄ > 2σk
n
k∈J :k>j̄
!
r !
r
X
2t
2t
+
≤ Pr
Tn,j̄ − EY > σj̄
Pr kTn,k − EY k > σk
n
n
k∈J : k>j̄


σmax
≤ 2de−t + 2d log2
e−t ,
σmin
where we used Theorem 3.1 to bound each of the probabilities in the sum. The display above
implies that the event
(
r )
\
2t
B=
kTn,k − EY k ≤ σk
n
k∈J :k≥j̄

of probability ≥ 1 − 2d log2



2σmax
σmin



e−t is contained in E = {j∗ ≤ j̄}. Hence, on B we have that
r

kTn∗

− EY k ≤

kTn∗

− Tn,j̄ k + kTn,j̄ − EY k ≤ 2σj̄
r
r
r
σn 2t
σn 2t
σn 2t
≤ 4√
+ 2√
= 6√
,
n n
n n
n n

2t
+ σj̄
n

r

2t
n

and result follows.
Remark 8. It follows from the proof that constant factor 6 in Theorem
to 3 + ε
 5.1 can be reduced
jσ
for any ε > 0 by considering the “finer grid”, that 
is, replacing
J
by
j
∈
Z
:
σ
≤
κ
min < κσmax


 min
for some 1 < κ < 2, at the cost of replacing log2

2σmax
σmin

by log2

κσmax
σmin

/ log2 κ.

S. Minsker/Estimators of the mean of a random matrix

22

6. From bounds depending on kEY 2 k to bounds depending on kE(Y − EY )2 k
Assume that Y1 , . . . , Yn are i.i.d. copies of Y ∈ Cd×d . In this section, we build upon previously
established bounds to provide performance guarantees for the estimator defined via (3.3), (3.4).
To this end, we study a version of the steepest descent scheme for the problem (3.3) initialized at
(0)
(0)
the point Tb , namely, T̂0 := Tb and
θ

θ0

n

1 X 
ψ θk (Yj − T̂k−1 ) , k ≥ 1
T̂k = T̂k−1 +
nθk
j=1

for anpappropriate choice of θk , k ≥ 0. Note that for any non-random self-adjoint matrix S and
1
θS = ns kE(Y −S)
2 k1/2 , Theorem 3.1 implies that
r 

s
2 1/2
≤ 2d exp (−s/2) ,
Pr kTn (S) − EY k ≥ kE(Y − S) k
n
P
where Tn (S) = S + nθ1S nj=1 ψ (θS (Yj − S)). Hence, if we use random S which is “not too far”
from EY with high probability, we expect that the deviation guarantees will still hold with the
“variance parameter” close to E(Y − EY )2 .
Everywhere in this section, we will assume that one has access to some known (possibly very
crude) bounds for σ 2 = kEY 2 k and σ02 = kE(Y − EY )2 k:
Assumption 1. Let σmin , σ0,min and σmax , σ0,max be known constants such that
σmin ≤ σ ≤ σmax and σ0,min ≤ σ0 ≤ σ0,max .
6.1. Two-step estimation based on sample splitting
We will first discuss the simplest (but not the most efficient) approach based on splitting the
sample Y1 , . . . , Yn into two disjoint subsets G1 and G2 of cardinality ≥ bn/2c each, and performing
one step of the steepest descent. The main advantage of this approach is the fact that it requires
very mild assumptions. The idea is to apply Lepski’s method (as discussed in section 5) twice:
on the first step, we obtain an estimator T̂0 based
n on subsample G1 , and on the
o second step we
apply Lepski’s method again to the subsample Yj − T̂0 : 1 ≤ j ≤ n, Yj ∈ G2 .
Here is the more detailed description: set σj = 2j σmin ,
J1 = {j ∈ Z : σmin ≤ σj < 2σmax }
and σ0,j = 2j σ0,min ,
(
J2 =

j ∈ Z : σ0,min ≤ σ0,j

r !)
t
< 2 σ0,max + 12σmax
,
n

and let T̂0 be the “Lepski-type” adaptive estimator based on the subsample G1 defined as
T̂0 = T|G1 |,j1∗ (0; G1 ),

S. Minsker/Estimators of the mean of a random matrix

where
T|G1 |,j (S; G1 ) =

23

|G1 |
1 X
ψ(θj (Yi − S)),
|G1 |θj
i=1

θj =

q

2t 1
n/2 σj ,

ψ(·) satisfies (3.1) and
(
j1∗

:= min j ∈ J1 :∀k ∈ J1 s.t. k > j,
s
T|G1 |,k (0; G1 ) − T|G1 |,j (0; G1 ) ≤ 2σk

2t
|G1 |

)

T̂1 is then defined as follows:
T̂1 = T̂0 + T|G2 |,j2∗ (T̂0 ; G2 ),
where
1
T|G2 |,j (S; G2 ) =
|G2 |θ0,j

n
X

s
ψ(θ0,j (Yi − S)),

θ0,j =

i=|G1 |+1

2t 1
n/2 σ0,j

and
(
j2∗

:= min j ∈ J2 :∀k ∈ J2 s.t. k > j,
s
T|G2 |,k (T̂0 ; G2 ) − T|G2 |,j (T̂0 ; G2 ) ≤ 2σ0,k

)
2t
.
|G2 |

Theorem 6.1. With probability at least

1 − 2d 2 + log2

σmax
σmin


+ log2

σ0,max + 12σmax
σ0,min

!!
p
t/n

e−t ,

the following inequality holds:
T̂1 − EY

r !r
t
t
≤ 12 σ0 + 12σ
.
n
n

Proof. See section B.3 in the appendix.
The main feature of this result is the variance term σ0 + 12σ
compared to σ as long as t  n.

q

t
n

that can be much smaller

S. Minsker/Estimators of the mean of a random matrix

24

6.2. Results for the estimator Tbθ∗ defined via equation (3.4)
We will next show how to design an estimator with deviations controlled by “correct” variance
term without sample splitting (however, subject to the condition that the sample size is sufficiently
large). In what follows, we will make an additional assumption about the function ψ:
Assumption 2. Function ψ(·) satisfies (3.1) and is operator Lipschitz, meaning that kψ(A) − ψ(B)k ≤
LkA−Bk for all self-adjoint A, B ∈ Cd×d , with Lipschitz constant L independent of the dimension
d.
For example, we may take ψ = ψ1 or ψ = ψ2 (see Lemma A.3 for details). As before, let t > 0 be
fixed, set σ0,j = 2j σ0,min ,

J = j ∈ Z : σ0,min ≤ σ0,j < 2σ0,max ,
r

(0)

For all j ∈ J , define δj

2t 1
θ=
and θj =
n σmax
q
= σmax 2t
n and
12
= σ0,j
5

(k)
δj

r

r

2t 1
for j ∈ J .
n σ0,j

r

2t
+ 6−k
n

σmax

2t 12
− σ0,j
n
5

r

2t
n

!
(6.1)

for k ≥ 1. Next, for each j ∈ J , we define
(0)

Tn,j := Tn(0) =

n
1 X
ψ (θYi ) ,
nθ

(6.2)

i=1

(independent of j)3 , and
(k)

(k−1)

Tn,j := Tn,j

+

n

1 X  
(k−1)
ψ θj Yi − Tn,j
nθj
i=1

for k ≥ 1. Finally, we apply Lepski’s method to the collection of estimators

n
o
(k)
Tn,j : j ∈ J . To

(k)

this end, define T̂k := Tn,j ∗ , where
k

n
jk∗ = min j ∈ J : ∀l ∈ J s.t. l > j,

(k)

(k)

Tn,l − Tn,j

(k)

≤ 2δl

o

.

Note that the estimator T̂k is completely data-dependent. We are ready to state the main result
of this section:
3

(0)

Particular choice of Tn

(0)

does not matter as long as kTn − EY k is small with high probability.

S. Minsker/Estimators of the mean of a random matrix

Theorem 6.2. Let

r
τ = 1.1K

d2 + Lt
+
n

r

25

2t 1
,
n 2σ0

where K > 0 is an absolute constant, and assume that τ ≤ 1/6. Moreover, assume that
r

24
2t
σ0,max ∨ σmax
≤ 1.
5
n

(6.3)

Then for all k ≥ 0 simultaneously,
"

T̂k − EY

24
≤ 3 (1 − 6−k ) σ0
5

r

2t
+ 6−k σmax
n

r

2t
n

#






2σ
12σmax
with probability ≥ 1 − 8d 1 + 2 log2 5σ
log2 σ 0,max e−t .
0,min

0,min

Proof. See section B.4 in Appendix.
The next corollary easily follows from the preceding result. Let A be the event of probability





2σ0,max
12σmax
Pr(A) ≥ 1 − 8d 1 + 2 log2
log2
e−t
5σ0,min
σ0,min
(k)

defined in Theorem 6.2. Since by the properties of the steepest descent scheme Tn,j converges to
the solution (denoted Tbθ∗j ) of the problem (3.3), we can easily deduce the following inequality.
n o
satisfy the equations
Corollary 6.1. Let Tbθ∗j
j∈J

n

1 X 
ψ θj (Yi − Tbθ∗j ) = 0d×d ,
nθj

j ∈ J.

i=1

Then on event A, Tbθ∗j − EY

(k)

≤ limk→∞ δj

=

12
5 σ0,j

q

2t
n.

One can further apply Lepski’s method (see section 5) to the collection
completely data-dependent estimator Tb∗ that satisfies
Tb∗ − EY

72
≤ σ0
5

with high probability (in particular, on event A).

r

2t
n

n o
Tbθ∗j

j∈J

to obtain a

S. Minsker/Estimators of the mean of a random matrix

26

7. Numerical simulation results
Numerical simulation was performed for covariance estimation problem. Data was simulated as
T
d
follows: let U = U (1) , . . . , U (100) ∈ R100 be a vector with i.i.d. coordinates such that U (j) =
1
√
(ξj,1 − ξj,2 ), where ξj,1 and ξj,2 , j = 1, . . . , 100, are independent random variables with
2c(q)

probability density function
pξ (t; q) =

q
I{t ≥ 0}
(1 + t)1+q

(which belongs to the Pareto family), c(q) = Var(ξ) = (q−1)2q (q−2) and q = 4.01; in particular,
√
Var(U (j) ) = 1. Finally, let Z = ΣU , where Σ is a diagonal matrix with Σ11 = 10, Σ22 =
1
5, Σ33 = 1, and Σjj = 97
, j ≥ 4. In particular, EZ = 0 and EZZ T = Σ.
The goal of numerical experiment was to evaluate the quality of estimation of the covariance
matrix Σ as well as its first eigenvector e1 corresponding to λ1 = 10. We tested two scenarios with
sample sizes equal n to 100 and 1000. In both cases, we generated Z1 , . . . , Zn , i.i.d. copies of Z,
and centered the data via the spatial (or geometric) median defined as
cn = argminy∈R100
M

100
X

ky − Zj k2 .

j=1

b n constructed as follows: set Z 0 := Zj − M
cn for brevity,
We compared two estimators, Sbn and Σ
j
and
n
1 X 0 0T
Sbn =
Zj Zj ,
n
j=1

which is the analogue of sample covariance with “robust centering”.
b n was constructed using a version of Lepski’s method described in section 5. We provide
Next, Σ
details for completeness: set
v
u
n
u
σmax
u 1 X 0 2 0 0T
σmax := 2t
kZj k2 Zj Zj , σmin =
,
n
100
j=1


J = j ∈ Z : σmin < 1.3j ≤ σmax ,
and let ψ(·) be the function defined in (3.6). Let t = log 10, and for j ∈ J , set θj =


P
T
Σ̂n,j = nθ1 j ni=1 ψ θj Zi0 Zi0 . Finally, define

q

2t 1
n 1.3j

and

(

r )
t
j∗ := min j ∈ J : ∀k > j, kΣ̂n,k − Σ̂n,j k ≤ 1.3k
n
(note that we modified some constants compared to the “theoretical” version), and finally set
b n := Σ̂n,j∗ .
Σ

S. Minsker/Estimators of the mean of a random matrix

27

n −Σk
n −Σk
Quality of covariance estimation was evaluated via comparing kSkΣk
with kΣkΣk
over 500
runs of simulations. We also compared errors of estimation of projectors onto the first principal
component,

b

u1 (Sbn )u1 (Sbn )T − u1 (Σ)u1 (Σ)T

and

b

b n )u1 (Σ
b n )T − u1 (Σ)u1 (Σ)T ,
u1 (Σ

where u1 (·) denotes the eigenvector corresponding to the largest eigenvalue of a matrix. Histograms
illustrating performance of both estimators are presented in figures 1a and 1b (for the sample size
n = 100), and in figures 2a and 2b (for the sample size n = 1000). It is clear from the graphs that
b n performs significantly better than Sbn .
in all scenarios, Σ
Acknowledgements. I want to thank L. Goldstein, A. Juditsky, and A. Nemirovski for making
many helpful and insightful suggestions that helped to improve the quality of presentation, as well
as X. Wei for suggesting the simple proof of Corollary D.1.
References
[AMS96] N. Alon, Y. Matias, and M. Szegedy. The space complexity of approximating the
frequency moments. In Proceedings of the twenty-eighth annual ACM symposium on
Theory of computing, pages 20–29. ACM, 1996.
[AP16] Alexei Borisovich Aleksandrov and Vladimir Vsevolodovich Peller. Operator lipschitz
functions. Russian Mathematical Surveys, 71(4):605, 2016.
[AW02] R. Ahlswede and A. Winter. Strong converse for identification via quantum channels.
IEEE Trans. Inform. Theory, 48(3):569–579, 2002.
[BDJ93] R. W. Butler, P. L. Davies, and M. Jhun. Asymptotics for the minimum covariance
determinant estimator. The Annals of Statistics, pages 1385–1400, 1993.
[Bha97] R. Bhatia. Matrix analysis. Springer, 1997.
[BJL15] C. Brownlees, E. Joly, and G. Lugosi. Empirical risk minimization for heavy-tailed
losses. The Annals of Statistics, 43(6):2507–2536, 2015.
[BLM13] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration inequalities:
A nonasymptotic theory of independence. Oxford university press, 2013.
[Car10] E. Carlen. Trace inequalities and quantum entropy: an introductory course. Available
at http://www.mathphys.org/AZschool/material/AZ09-carlen.pdf, 2010.
[Cat12] O. Catoni. Challenging the empirical mean and empirical variance: a deviation study.
In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, volume 48, pages
1148–1185, 2012.
[Cat16] O. Catoni. PAC-Bayesian bounds for the Gram matrix and least squares regression
with a random design. arXiv preprint arXiv:1603.05229, 2016.
[CLMW11] E. J. Candès, X. Li, Y. Ma, and J. Wright. Robust principal component analysis?
Journal of the ACM (JACM), 58(3):11, 2011.
[CRZ16] T. T. Cai, Z. Ren, and H. H. Zhou. Estimating structured high-dimensional covariance
and precision matrices: optimal rates and adaptive estimation. Electron. J. Statist.,
10(1):1–59, 2016.
[CZZ10] T. T. Cai, C.-H. Zhang, and H. H. Zhou. Optimal rates of convergence for covariance
matrix estimation. The Annals of Statistics, 38(4):2118–2144, 2010.

S. Minsker/Estimators of the mean of a random matrix

28

[Dav92] L. Davies. The asymptotics of Rousseeuw’s minimum volume ellipsoid estimator. The
Annals of Statistics, pages 1828–1843, 1992.
[Dir13] S. Dirksen. Tail bounds via generic chaining. arXiv preprint arXiv:1309.3522, 2013.
[DLLO15] L. Devroye, M. Lerasle, G. Lugosi, and R. I. Oliveira. Sub-Gaussian mean estimators.
arXiv preprint arXiv:1509.05845, 2015.
[FWZ16a] J. Fan, W. Wang, and Z. Zhu. Robust low-rank matrix recovery. arXiv preprint
arXiv:1603.08315, 2016.
[FWZ16b] Jianqing Fan, Weichen Wang, and Yiqiao Zhong. An `∞ eigenvector perturbation bound and its application to robust covariance estimation. arXiv preprint
arXiv:1603.03516, 2016.
[Giu15] I. Giulini. PAC-Bayesian bounds for Principal Component Analysis in Hilbert spaces.
arXiv preprint arXiv:1511.06263, 2015.
[HR09] P. J. Huber and E. M. Ronchetti. Robust statistics. Wiley Series in Probability and
Statistics. John Wiley & Sons Inc., Hoboken, NJ, second edition, 2009.
[HRVA08] M. Hubert, P. J. Rousseeuw, and S. Van Aelst. High-breakdown robust multivariate
methods. Statistical Science, pages 92–119, 2008.
[HS16] D. Hsu and S. Sabato. Loss minimization and parameter estimation with heavy tails.
Journal of Machine Learning Research, 17(18):1–40, 2016.
[Hub64] P. J. Huber. Robust estimation of a location parameter. The Annals of Mathematical
Statistics, 35(1):73–101, 1964.
[JLO+ 17] Emilien Joly, Gábor Lugosi, Roberto Imbuzeiro Oliveira, et al. On the estimation of
the mean of a random vector. Electronic Journal of Statistics, 11(1):440–451, 2017.
[JVV86] Mark R Jerrum, Leslie G Valiant, and Vijay V Vazirani. Random generation of
combinatorial structures from a uniform distribution. Theoretical Computer Science,
43:169–188, 1986.
[KL16] V. Koltchinskii and K. Lounici. New asymptotic results in principal component
analysis. arXiv preprint arXiv:1601.01457, 2016.
[KL17] Vladimir Koltchinskii and Karim Lounici. Concentration inequalities and moment
bounds for sample covariance operators. Bernoulli, 23(1):110–133, 2017.
[KLT11] V. Koltchinskii, K. Lounici, and A. B. Tsybakov. Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion. The Annals of Statistics, 39(5):2302–
2329, 2011.
[KLT14] O. Klopp, K. Lounici, and A. B. Tsybakov. Robust matrix completion. arXiv preprint
arXiv:1412.8132, 2014.
[Lep92] O. Lepski. Asymptotically minimax adaptive estimation. I: Upper bounds. Optimally
adaptive estimates. Theory of Probability & Its Applications, 36(4):682–697, 1992.
[LF09] Clifford Lam and Jianqing Fan. Sparsistency and rates of convergence in large covariance matrix estimation. Annals of statistics, 37(6B):4254, 2009.
[Lie73] E. H. Lieb. Convex trace functions and the Wigner-Yanase-Dyson conjecture. Advances in Math., 11:267–288, 1973.
[LM17] Gábor Lugosi and Shahar Mendelson. Sub-Gaussian estimators of the mean of a
random vector. arXiv preprint arXiv:1702.00482, 2017.
[LO11] Matthieu Lerasle and Roberto I Oliveira. Robust empirical mean estimators. arXiv

S. Minsker/Estimators of the mean of a random matrix

29

preprint arXiv:1112.3914, 2011.
[Lou14] K. Lounici. High-dimensional covariance matrix estimation with missing observations.
Bernoulli, 20(3):1029–1058, 2014.
[Mar76] R. A. Maronna. Robust M-estimators of multivariate location and scatter. Ann.
Statist., 4(1):51–67, 01 1976.
[Min15] S. Minsker. Geometric median and robust estimation in Banach spaces. Bernoulli,
21(4):2308–2335, 2015.
[MW17] S. Minsker and X. Wei. Estimation of the covariance structure of heavy-tailed distributions. arXiv preprint arXiv:1708.00502, 2017.
[NY83] A. Nemirovski and D. Yudin. Problem complexity and method efficiency in optimization. John Wiley & Sons Inc., 1983.
[Oli09] R. I. Oliveira. Concentration of the adjacency matrix and of the Laplacian in random
graphs with independent edges. arXiv preprint arXiv:0911.0600, 2009.
[SV13] N. Srivastava and R. Vershynin. Covariance estimation for distributions with 2 + ε
moments. The Annals of Probability, 41(5):3081–3111, 2013.
[Tal14] M. Talagrand. Upper and lower bounds for stochastic processes: modern methods and
classical problems, volume 60. Springer Science & Business Media, 2014.
[Tro12a] J. A. Tropp. User-friendly tail bounds for sums of random matrices. Found. Comput.
Math., 12(4):389–434, 2012.
[Tro12b] J. A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of
computational mathematics, 12(4):389–434, 2012.
[Tro15] J. A. Tropp. An introduction to matrix concentration inequalities. arXiv preprint
arXiv:1501.01571, 2015.
[Tyl87] D. E. Tyler. A distribution-free M-estimator of multivariate scatter. The Annals of
Statistics, pages 234–251, 1987.
[Ver07] R. Vershynin. Non-asymptotic theory of random matrices: lecture notes. Available
at http://www-personal.umich.edu/~romanv/teaching/2006-07/280/lec6.pdf,
2007.
[Ver10] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv
preprint arXiv:1011.3027, 2010.
[ZCS16] T. Zhang, X. Cheng, and A. Singer. Marčenko-Pastur law for Tylers M-estimator.
Journal of Multivariate Analysis, 149:114–123, 2016.
Appendix A: Technical results
Lemma A.1. Let F : R 7→ R be a continuously differentiable function, and S ∈ Cd×d be a
self-adjoint matrix. Then the gradient of G(S) := tr F (S) is
∇G(S) = F 0 (S),
where F 0 is the derivative of F and F 0 (S) : Cd×d 7→ Cd×d is the matrix function in the sense of
definition 2.1.
Proof. We will first check the claim assuming that F is a polynomial of the form F (x) = xk , k ∈ N.
Let H = H ∗ be a self-adjoint operator, and consider the directional derivative dG(S; H) of G in

S. Minsker/Estimators of the mean of a random matrix

30

direction H:
k
 X


1 
k
k
dG(S; H) = lim tr (S + tH) − S =
tr S j−1 HS k−j
t→0 t
j=1


= tr kS k−1 H = F 0 (S), H ,

hence the claim holds for monomials. By linearity, it also holds for arbitrary polynomials. It
remains to extend the claim to arbitrary continuously differentiable function via a standard approximation argument (for instance, see [Bha97, chapter 5, section 3]).
Lemma A.2. Let 1 < α ≤ 2 and cα =

α−1
α

∨

q

2−α
α .

Then 1 + y + cα |y|α > 0 and

− log(1 + y + cα |y|α ) ≤ log(1 − y + cα |y|α )

for all y ∈ R.

Proof. To check the first claim, it is enough to note that f (y) = 1 + y + cα |y|α is convex and its
1/(α−1)

. It is easy to check that f (ym ) = 1 − ym + yαm ,
minimum is attained for ym = − αc1α
which implies that f (ym ) > 0 ⇐⇒ cα > α−1
which always holds since cα ≥ α−1
α and α > 1.
α2
For the second part, it is enough to show that (1 + cα |y|α + y)(1 + cα |y|α − y) ≥ 1 for all y ∈ R,
which is equivalent to claiming that c2α y 2α + 2cα y α ≥ y 2 , y ≥ 0. Note that for any τ ∈ (−1, 1),
p, q > 0 such that 1/p + 1/q = 1, and y ≥ 0,
y 2 = y 1−τ y 1+τ ≤
Choosing p :=

α
2(α−1) ,

q :=

by 2cα y α + c2α y 2α for cα =

α
2−α ,
α−1
α

we get y 2 ≤
q
∨ 2−α
α .

y p(1−τ ) y q(1+τ )
+
.
p
q

2(α−1) α
y
α

+

2−α 2α
α y

which is further bounded above

Lemma A.3. Functions ψ1 (x) and ψ2 (x) defined in Remark 1 are operator Lipschitz, with Lipschitz constants independent of the dimension.
Proof. Lipshitz property of ψ1 (x) follows from Theorem 1.6.1 in [AP16]. Result for ψ2 (x) follows
from Theorem 1.1.1 in the same paper.
A.1. Proof of Lemma 2.1
For a self-adjoint matrices R, Q, kRk ≥ kQk iff kR2 k ≥ kQ2 k. Clearly,


S A
A∗ T

2


=

S 2 + AA∗ SA + AT
A∗ S + T A∗ T 2 + A∗ A



S. Minsker/Estimators of the mean of a random matrix

2
S A
It implies that
≥ S 2 + AA∗ ≥ kAA∗ k and
A∗ T

2 

0 A
AA∗
0
∗
kA Ak. Since
=
, we obtain
A∗ 0
0
A∗ A




S A
A∗ T

2



0 A
A∗ 0

≥



S A
A∗ T

31

2
≥

T 2 + A∗ A

≥

2
,

and result follows.
Appendix B: Proofs omitted in the main text
B.1. Tools from probability theory and linear algebra
We recall several useful results that we need in the proofs below.
Lemma B.1 (Matrix Hoeffding inequality). Let Z1 , . . . , Zn ∈ Cd×d be a sequence of independent
self-adjoint random matrices such that for all 1 ≤ k ≤ n,

Then

Pn

j=1 Zj

EZk = 0 and kZk k ≤ Mk almost surely.


2
t
≤ t with probability ≥ 1 − 2d exp − 8 Pn M 2 .
j=1

j

Proof. See Theorem 1.3 in [Tro12b].
We conclude this section by recalling the notion of Talagrand’s generic chaining complexity (see
[Tal14]) and several related results. Given a metric space (T, ρ), let {∆n } be a nested sequence
n
of partitions of T such that card ∆0 = 1 and card ∆n ≤ 22 . For s ∈ T , let ∆n (s) be the unique
subset of ∆n containing s. The generic chaining complexity γ2 (T, ρ) is defined as
X n
γ2 (T, ρ) := inf sup
2 2 D(∆n (s))
{∆n } s∈T

n≥0

where the infimum is taken over all admissible sequences of partitions and D(A) := D(A, ρ) stands
for the diameter of a set A. The covering number N (T, ρ, ε) is defined as the smallest N ∈ N such
that there exists a subset F ⊆ T of cardinality N with the property that for all z ∈ T , ρ(z, F ) ≤ ε.
Dudley’s entropy integral bound (see [Tal14]) states that
1
γ2 (T, ρ) ≤ √
2 2−1

D(T
Z )

p
log N (T, ρ, ε/4)dε.

(B.1)

0

Cd×d -valued

We will say that
stochastic process {X(t), t ∈ T } has sub-Gaussian increments with
respect to the metric ρ if for all t1 , t2 ∈ T,
Pr (kXt1 − Xt2 k ≥ sρ(t1 , t2 )) ≤ 2de−s
where k · k is the operator norm.

2 /2

,

S. Minsker/Estimators of the mean of a random matrix

32

Lemma B.2. Let (T, ρ) be a metric space and let Cd×d -valued stochastic process {X(t), t ∈ T }
have sub-Gaussian increments with respect to ρ. There exists an absolute constant C > 0 such
that for any t0 ∈ T and any s ≥ 1,

√
sup kXt − Xt0 k ≤ C γ2 (T, ρ) + sD(T )
t∈T

with probability ≥ 1 − 2de−s .
Proof. See Theorem 3.2 in [Dir13] for a more general statement.
B.2. Proof of Theorem 3.2
Define φ(x) = max(ex − 1, 0) and Xj = ψ(θYj ). Proceeding as in the proof of Theorem 3.1, we
get that for s ≥ 0,

 




!
n
n
X
X
1
(Xj − θEYj ) ≥ s = Pr φ λmax 
(Xj − θEYj ) ≥ φ(θs)
Pr λmax 
θ
j=1
j=1






n
n
X
X
1
1
Etr exp 
Etr φ 
(Xj − θEYj ) =
(Xj − θEYj ) − I  .
≤
φ(θs)
φ(θs)
j=1

It follows from Lemma 3.1 that

Etr exp 

n
X
j=1

Set Bn2 :=


Pn

2
j=1 EYj



tr exp 


≤ tr 

n
θ2 X

2

j=1

2


Bn2 1



j=1

 0, and note that




+

θ2
2
2 kBn k

2!




n
2 X
θ
(Xj − θEYj ) ≤ tr exp 
EYj2  .
2



 θ2 p 2 
Bn I +
EYj2  − I  = tr 
2


θ2

j=1


+ ... +

θ2
2
2 kBn k

k!

k−1

θ2 2
2 Bn

2!


+ ... +

θ2 2
2 Bn

k−1

k!





p 
+ . . . Bn2 


 tr Bn2
+ . . .  =
kBn2 k


 2


θ
2
exp
kBn k − 1 .
2

p
Here we have used the fact that A  B implies SAS ∗  SBS ∗ for S = S ∗ := Bn2 , and the
P
x
xj−1
equality e x−1 = ∞
j=1 j! . We have shown that
 2



!
θ
2k − 1
n
exp
kB
2
X
n
2
1
tr Bn
Pr λmax 
(Xj − θEYj ) ≥ s ≤
2
θs
θ
kBn k
e −1
j=1
 2

 2


tr Bn2
θ
eθs
tr Bn2
θ
1
2
2
≤
exp
kB
k
−
θs
≤
exp
kB
k
−
θs
1
+
,
n
n
kBn2 k
2
kBn2 k
2
θs
eθs − 1

S. Minsker/Estimators of the mean of a random matrix

33

θs

e
1
where we used an elementary inequality eθs
≤ 1 + θs
on the last step.
−1
Combining the same steps with Fact 2.4 and the equality −λmin (A) = λmax (−A), we get


!

 2

n
X
θ
1
1
tr Bn2
2
exp
kB
k
−
θs
1
+
.
Pr λmin 
(Xj − θEYj ) ≤ −s ≤
n
θ
kBn2 k
2
θs
j=1

√
Finally, replace s by t n to get the bound in the required form.
B.3. Proof of Theorem 6.1
Let E1 be the event defined by
s

(
E1 =

T̂0 − EY

≤ 6σ

2t
n/2

)
.

By Theorem 5.1,

Pr(E1 ) ≥ 1 − 2d log2

2σmax
σmin



e−t .

(B.2)

Note that on this event,
h

2

E (Y − T̂0 ) |T̂0

i

2

2

= kE(Y − EY ) + (EY − T̂0 ) k ≤

σ02

+

r !2
t
.
12σ
n

(B.3)

In particular, on event E1 ,
h
i
E (Y − T̂0 )2 |T̂0

1/2

r
≤ σ0 + 12σ

t
≤ σ0,max + 12σmax
n

r

t
.
n

(B.4)

Next,
Pr

T̂1 − EY

!
r !s
t
2t
≥ 6 σ0 + 12σ
n
n/2

(B.5)

!
r !s
t
2t
= Pr
T|G2 |,j2∗ (T̂0 ; G2 ) − (EY − T̂0 ) ≥ 6 σ0 + 12σ
n
n/2
!
r !s
t
2t
c
≤ Pr(E1 ) + Pr
T|G2 |,j2∗ (T̂0 ; G2 ) − (EY − T̂0 ) ≥ 6 σ0 + 12σ
E1 .
n
n/2
f
Define the new probability measure by Pr(A)
= Pr(A|E1 ). Clearly, under this new measure,
subsample G2 is still independent of G1 since E1 ∈ σ(G1 ) - the sigma-algebra generated by G1 ,

S. Minsker/Estimators of the mean of a random matrix

34

f
f
and for any B ∈ σ(G2 ), Pr(B)
= Pr(B). Let Ẽ be the expectation with respect to measure Pr(·).
Then
r !s


t
2t
Pr kT|G2 |,j2∗ (T̂0 ; G2 ) − (EY − T̂0 )k ≥ 6 σ0 + 12σ
E1
n
n/2
!
r !s
t
2t
f kT|G |,j ∗ (T̂0 ; G2 ) − (EY − T̂0 )k ≥ 6 σ0 + 12σ
= Pr
2 2
n
n/2
s
!
h
i 1/2
2t
f kT|G |,j ∗ (T̂0 ; G2 ) − (EY − T̂0 )k ≥ 6 E (Y − T̂0 )2 |T̂0
≤ Pr
2 2
n/2
s
!
h
i 1/2
2t
2
f
T̂0
= ẼPr
T|G2 |,j2∗ (T̂0 ; G2 ) − (EY − T̂0 ) ≥ 6 E (Y − T̂0 ) |T̂0
n/2
!
p
2(σ0,max + 12σmax t/n)
≤ 1 − 2d log2
e−t .
σ0,min
f on the first step and (B.4) on the second step. The last inequalHere, we use the definition of Pr(·)
f
ity follows from independence of G2 from T̂0 (under Pr(·))
and Theorem 5.1 applied conditionally
on T̂0 : indeed, this can be done since (B.4) holds on E1 . It remains to combine the last bound
with (B.5) and (B.2).
B.4. Proof of Theorem 6.2
We will first state several technical results that are required in the proof. Let j ∈ J be such that
σ0,j = 2j σ0,min ≥ σ0 , and define
Xj,i (S) := ψ (θj (Yi − S)) , i = 1, . . . , n.
Moreover, set
Ln (δ, j) :=

sup
S:kS−EY k≤δ

n
n
1 X
1 X
(Xj,i (S) − EXj,i (S)) −
(Xj,i (EY ) − EXj,i (EY )) ,
nθj
nθj
i=1

i=1

and define the event

(
Ω(δ, j) =

Ln (δ, j) ≤ Kδ

r

d2 + Lt
n

where K > 0 is an absolute constant.
Lemma B.3. For K large enough,
Pr (Ω(δ, j)) ≥ 1 − 2de−t .
Proof. See section B.4.1.

)
,

S. Minsker/Estimators of the mean of a random matrix

35

Lemma B.4. For any Hermitian S,
S − EY +

θj
1
EXj,1 (S) ≤
E(Y − S)2 .
θj
2

Proof. Note that
θj2
Xj,1 (S) = ψ (θj (Y1 − S))  log I + θj (Y1 − S) + (Y1 − S)2
2

!

θj2
(Y1 − S)2 ,
2
which is a consequence of scalar inequality log(1 + x) ≤ x, x > −1 and Fact 2.1, hence we can
deduce from Fact 2.2 that


θj
1
λmax S − EY + EXj,1 (S) ≤
E(Y − S)2 .
θj
2
 θj (Y1 − S) +

At the same time,




1
1
−λmin S − EY + EXj,1 (S) = λmax − EXj,1 (S) − (S − EY ) .
θj
θj
Since
θj2
−Xj,1 (S) = −ψ (θj (Y1 − S))  log I − θj (Y1 − S) + (Y1 − S)2
2

!

by the definition of ψ(·), we conclude that


θj
1
E(Y − S)2 ,
−λmin S − EY − EXj,1 (S) ≤
θj
2
hence S − EY +

1
θj EXj,1 (S)

≤

θj
2

E(Y − S)2 .

Lemma B.5. With probability ≥ 1 − 2de−t ,
n


1 X
Xj,i (EY ) − E Xj,i (EY )
≤ σ0,j
nθj
i=1

r

2t θj 2
+ σ0 .
n
2

Proof. Result follows from Theorem 3.1 and the inequality
θj
1
EXj,1 (EY ) ≤ σ02 ,
θj
2
which is a consequence of Lemma B.4. Indeed,
n
n


1 X
1 X
1
≤
Xj,i (EY ) − E Xj,i (EY )
Xj,i (EY ) + kEXj,1 (EY )k
nθj
nθj
θj
i=1
i=1
r
2t θj 2
≤ σ0,j
+ σ0
n
2

with probability ≥ 1 − 2de−t .

S. Minsker/Estimators of the mean of a random matrix

We are ready to proceed with the proof of the theorem. Let
(
r
Tn(0)

E0 =
(0)

(0)

where Tn was defined in (6.2) as Tn
Theorem 3.1. Let

=

− EY

1
nθ

≤ σmax

Pn

i=1 ψ (θYi ),

)
,

and note that Pr(E0 ) ≥ 1 − 2de−t by


12
kmax = 1 + max k ≥ 0 : σ0,min 1.1 ≤ σmax ,
5
r
2t
γl = 1.1l σ0,min
, l ≥ 0,
n

 
 log 12σmax 


 2 5σ0,min 
 ≤ 1 + 8 log2 12σmax . Define
≤1+
log 1.1
5σ


and note that kmax

2t
n

36

k

0,min

2

(
Ωj :=

(B.6)

n
1 X
(Xi,EY − EXi,EY ) ≤ σ0,j
nθj

r

i=1

2t θj 2
+ σ0
n
2

)
∩

k\
max

Ω (γl , j) .

l=0

By Lemma B.3, Lemma B.5 and the union bound, Pr (Ωj ) ≥ 1 − 2d(2 + kmax )e−t . We will now
(k)

(k)

show by induction that on the event E0 ∩ Ωj , Tn,j − EY ≤ δj for all k ≥ 0. For k = 0, result
followsnfrom the definition of E0o
. In remains to complete the induction step k − 1 7→ k. Note that
(k−1)
(k−1)
when
Tn,j − EY ≤ δj
occurs, we have
(k)
Tn,j

− EY

=

(k−1)
Tn,j

n

1 X 
(k−1)
− EY +
ψ θj (Yi − Tn,j )
nθj
i=1

≤

sup
(k−1)

S:kS−EY k≤δj

n
1 X
S − EY +
ψ (θj (Yi − S)) .
nθj

(B.7)

i=1

Expression under the supremum in (B.7) can be decomposed as follows:
n
1 X
ψ (θj (Yi − S)) =
S − EY +
nθj
i=1

S − EY +

n
1 X
1
EXj,1 (S) +
(Xj,i (S) − EXj,i (S))
θj
nθj
i=1

n
n
1 X
1 X
−
(Xj,i (EY ) − EXj,i (EY )) +
(Xj,i (EY ) − EXj,i (EY )) .
nθj
nθj
i=1

i=1

We will treat 3 terms separately: first, it follows from Lemma B.4 that on Ωj


 
θj
1
(k−1) 2
2
sup
S − EY + EX1,S ≤
σ0 + δ j
.
θj
2
(k−1)
S:kS−EY k≤δ
j

(B.8)

S. Minsker/Estimators of the mean of a random matrix

37

Next,
n
1 X
(Xi,EY − EXi,EY ) ≤ σ0,j
nθj
i=1

r

2t θj 2
+ σ0 ,
n
2

(B.9)

n
o
(k−1)
once again by the definition of Ωj . Let ˜l = min l ≥ 0 : γl ≥ δj
(where γl was defined in
(k−1)
(B.6)), and note that ˜l ≤ kmax and γl̃ ≤ 1.1δj
. We bound the third term as
n
n
1 X
1 X
sup
(Xj,i (S) − EXj,i (S)) −
(Xj,i (EY ) − EXj,i (EY ))
nθj
nθj
(k−1)
S:kS−EY k≤δj
i=1
i=1
r
r



d2 + Lt
d2 + Lt
(k−1)
(k−1)
= Ln δ j
, j ≤ Ln γl̃ , j ≤ Kγl̃
≤ 1.1Kδj
.
n
n

(B.10)

Putting the bounds (B.8),(B.9),(B.10) together, we can estimate the supremum in (B.7) as
n
1 X
ψ (θj (Yi − S))
nθj
(k−1)
S:kS−EY k≤δj
i=1
r
r


  θ
θj
2t
d2 + Lt
j 2
(k−1) 2
(k−1)
2
σ0 + δj
+ δj
· 1.1K
≤
+ σ0 + σ0,j
2
2
n
n
!
r
r
r
2t
d2 + Lt
2t 1
(k−1)
≤ (σ0 + σ0,j )
+ δj
+
1.1K
.
n
n
n 2σ0,j

sup

S − EY +

Note that we have used bounds

θj σ02

≤ σ0

q

2t
n

(B.11)



(k−1) 2
(k−1)
and θj δj
≤ θj δ j
(indeed, inequality

(m)

≤ 1 for all j and m) to get the second inequality above. Since j was chosen
q
q
2
2t 1
1
+
≥ σ0 and τ = 1.1K d +Lt
n
n 2σ0 ≤ 6 by assumption, we have shown that

(6.3) implies that δj
such that σ0,j

(k)
Tn,j

r
− EY

≤ 2σ0,j

2t 1 (k−1)
(k)
+ δj
= δj ,
n
6
(k)

where the last equality follows from the fact that the sequence δj defined in (6.1) satisfies the
recursive relation
r
r
2t
2t 1 (k−1)
(0)
(k)
δj = σmax
, δj = 2σ0,j
+ δj
.
n
n
6
To complete the proof, it
n is enough toofollow the steps of the proof Theorem 5.1 applied to the
(k)
collection of estimators Tn,j : j ∈ J : first, let j̄ = min {j ∈ J : σ0,j ≥ σ0 }, and note that the
event
\
E0 ∩
Ωj
j≥j̄, j∈J

S. Minsker/Estimators of the mean of a random matrix

38






2σ0,max
12σmax
has probability ≥ 1−8d 1 + 2 log2 5σ
log2 σ
e−t . Moreover, on this event jk∗ ≤
0,min

0,min

j̄, hence
T̂k − EY

(k)

(k)

= Tn,j ∗ − EY
k
"
(k)

≤ 3δj̄

(k)

(k)

≤ Tn,j ∗ − Tn,j̄ + Tn,j̄ − EY
k
r
r #
24
2t
2t
≤ 3 (1 − 6−k ) σ0
,
+ 6−k σmax
5
n
n

where we used the fact that σ0,j̄ ≤ 2σ0 in the last inequality.
B.4.1. Proof of Lemma B.3
To this end, we will use a chaining argument. Recall that the function ψ(·) is operator Lipschitz
with Lipschitz constant L by assumption. Recall that Xj,i (S) := ψ (θj (Yi − S)) , i = 1, . . . , n. It
follows from Assumption 2 (see also Lemma A.3) that for any Hermitian S1 , S2 and 1 ≤ i ≤ n,
kXi,S1 − Xi,S2 k = kψ(θj (Yi − S1 )) − ψ(θj (Yi − S2 )k ≤ Lθj kS1 − S2 k.
Matrix Hoeffding’s inequality (Lemma B.1) applies with
Zi =
and Mi =

2L
n kS1

1
((Xi,S1 − EXi,S1 ) − (Xi,S2 − EXi,S2 )) , i = 1, . . . , n,
nθj

− S2 k, and yields that
n
X
i=1

Zi

√

r
s
≤ L 32kS1 − S2 k
n

with probability ≥ 1 − 2de−s .
Lemma B.6 (Covering number in the operator norm). Let B(r) be the ball of radius r > 0 in
2
Rd with respect to the operator norm k · k, centered at 0. Then the covering number N (B(r), ε) :=
N (B(r), k · k, ε) satisfies

d2
2r
N (B(r), ε) ≤
+1
.
ε
Proof. It is well known [Ver07] that
N (A, ε) ≤

|A + B(ε/2)|
,
|B(ε/2)|

where |C| denotes the Lebesgue measure of a set C, and A + C stands for the Minkowski sum of
the sets A and C. For A = B(r), we get N (B(r), ε) ≤ |B(r+ε/2)|
|B(ε/2)| . The volume of the unit ball is
given by
Z
Y
|B(r)| = cd
x2i − x2j dx1 . . . dxd ,
[−r,r]d

1≤i<j≤d

S. Minsker/Estimators of the mean of a random matrix

39

Q
2
d
2
where cd = d!4−d
and vj is the volume of the Euclidean unit ball in Rj . From here, it
j=1 vj
is easy to see that

d2
2r
|B(r + ε/2)|
=
+1
.
|B(ε/2)|
ε

Let T (δk−1 ) := S ∈ Cd×d : kS − EY k ≤ δk−1 , and define the metric
ρd (S1 , S2 ) := LkS1 − S2 k, S1 , S2 ∈ Cd×d .
P
Viewing S 7→ nθ1 j ni=1 (Xi,S −EXi,S ) as a Cd×d -valued stochastic process indexed by the elements
of the metric space (T (δk−1 ), ρd ), we can apply Lemma B.2 which implies that there exists an
absolute constant C > 0 such that for any t ≥ 1,

√
C 
Ln (δk−1 ) ≤ √ γ2 (T (δk−1 ), ρd ) + tD(T (δk−1 ), ρd )
(B.12)
n
with probability ≥ 1 − 2de−t . Recall the Dudley’s entropy integral bound (B.1):
1
γ2 (T (δk−1 ), ρd ) ≤ √
2 2−1

D(T (δZk−1 ),ρd )

p
log N (T (δk−1 ), ρd , ε/4)dε.

0

Noting that D(T (δk−1 ), ρd ) = 2L δk−1 and combining Dudley’s bound with the estimate of Lemma
B.6, we get
γ2 (T (δk−1 ), ρd ) ≤ C1 L δk−1 d,
R1
1/2
2
where C1 = 2√2−1 0 log (1 + 4/ε)dε. Bound (B.12) implies that with probability ≥ 1 − 2de−t ,
r
√
C 
d2 + Lt
Ln (δk−1 ) ≤ √ C1 Lδk−1 d + 2Lδk−1 t ≤ δk−1 · K
(B.13)
n
n
for some absolute constant K > 0.
Appendix C: Sub-Gaussian bounds that hold simultaneously for multiple
confidence levels
Sub-Gaussian tail guarantees of Theorem 3.1 do not hold for multiple confidence levels simultaneously. Here, we show how to improve on this aspect. Let 0 < tmin < tmax be fixed, and define
tj = 2j tmin ,
J = {j ∈ Z : tmin < tj < 2tmax } .
Assume that Y1 , . . . , Yn ∈ Cd×d is a sequence
q of independent self-adjoint random matrices, and
P
2tj
n
1√
2
let σn2 ≥
j=1 EYj . For j ∈ J , set θj =
n σn / n , and
Tn,j =

n
1 X
ψ(θj Yi ).
nθj
i=1

S. Minsker/Estimators of the mean of a random matrix

Define

(

σn
j∗ := min j ∈ J : ∀k ∈ J s.t. k > j, kTn,k − Tn,j k ≤ 2 √
n

40

r

2tk
n

)
,

and Tn∗ := Tn,j∗ .
Theorem C.1. For any 0 < tmin ≤ t ≤ tmax ,
σn
Pr kTn∗ − EY k ≥ 6 √
n

r !
e−t
t
≤ 2d
.
n
1 − e−t

Proof. Let tmin ≤ t ≤ tmax , and define j(t) = min {j ∈ J : tj ≥ t}; note that tj(t) ≤ 2t. Next,
we have that
(
)
r
[
σn 2tk
{j∗ > j(t)} =
Tn,k − Tn,j(t) > 2 √
n
n
k∈J :k>j(t)
) (
)
(
r
r
[
σn 2tj(t)
σn 2tk [
Tn,j(t) − EY > √
:= E c (t)
⊆
kTn,k − EY k > √
n
n
n
n
k∈J :k>j(t)

By union bound, probability of the latter event E c (t) is bounded by
2de−tj(t) +

X

2de−tk ≤ 2d

k∈J :k>j(t)

X

l

e−tj(t) ·2 ≤ 2d

l≥0

e−t
.
1 − e−t

Next, on the event E(t) := (E c (t))c ⊆ {j∗ ≤ j(t)}, we have
kTn,j∗

σn
− EY k ≤ kTn,j∗ − Tn,j(t) k + kTn,j(t) − EY k ≤ 2 √
n
r
r
√ σn 2t √ σn 2t
≤ 2 2√
+ 2√
,
n n
n n

r

2tj(t)
σn
+√
n
n

r

2tj(t)
n

and conclusion follows.
Appendix D: Supplementary results
The following lemma is a generalization of Chebyshev’s association inequality.
Lemma D.1 (FKG inequality). Let f, g : Rd → R be two functions that are non-decreasing with
respect to each coordinate. Moreover, let V = (V1 , V2 , . . . , Vd ) be a random vector taking values in
Rd . Then
Ef (V )g(V ) ≥ Ef (V )Eg(V ).
Proof. See Theorem 2.15 in [BLM13].
The following corollary is immediate.

S. Minsker/Estimators of the mean of a random matrix

41

Corollary D.1. Let Z ∈ Rd be a centered random vector with covariance matrix Σ. Then
σ02 := EkZk22 ZZ T ≥ EkZk22 EZZ T = tr Σ kΣk.


Proof. Consider any unit vector v ∈ Rd . It is enough to show E (vT Z)2 kZk22 ≥ E(vT Z)2 EkZk22 .
We make the change the coordinates by considering an orthonormal basis {v1 , · · · , vd } with
v1 = v. Letting Vi = viT Z, i = 1, 2, · · · , d, we obtain
E(vT Z)2 kZk22 = EV12 kV k22 ≥ EV12 EkV k22 ,
where the last inequality follows from Lemma D.1 inequality by setting f (V1 , . . . , Vd ) := V12 and
g (V1 , . . . , Vd ) := kV k22 .
Lemma D.2. Let γ1 , . . . , γn be independent random variables with density p(x) = e−2|x| , x ∈ R.
Then for all n > 1
1. E max (|γ1 |, . . . , |γn |) ≥ 21 log n.
2. E max (γ1 , . . . , γn ) ≥ 14 log n.


3. Pr max (|γ1 |, . . . , |γn |) ≥ 12 − τ log n ≥ c(τ ) > 0 for every 0 < τ < 1/2.
Proof. Note that, since the distribution of γj ’s is symmetric and E max (γ1 , . . . , γn ) is positive,
E max (|γ1 |, . . . , |γn |) = E max (max (γ1 , . . . , γn ) , max (−γ1 , . . . , −γn ))
≤ 2E max (γ1 , . . . , γn ) .
Next, for 1 ≤ j ≤ n, |γj | has exponential distribution with P
density p̃(x) = 2e−2x I{x ≥ 0}. It
1
follows from a well-known fact that E max (|γ1 |, . . . , |γn |) = 2 nj=1 1j ≥ 12 log n, and the first and
second inequalities follow.
A standard computation shows that E max2 (|γ1 |, . . . , |γn |) ≤ c1 log2 n for some numerical constant c1 > 0, hence Paley-Zygmund inequality implies the last claim.

S. Minsker/Estimators of the mean of a random matrix

42

Frequency

0.35

Sample covariance estimator
0.3

Robust covariance estimator
Robus t es timator er r or
b n − Σk/kΣk
kΣ

0.25

0.2

0.15

0.1

Sample covar iance er r or
k Sbn − Σk/kΣk

0.05

0

1

2

3

4

5

6

7

8

9

10

11

Error

(a) Covariance matrix estimation error

Frequency

0.9

Sample covariance estimator
0.8

Robust covariance estimator
0.7

b n)u( Σ
b n) T − u 1 (Σ)u 1 (Σ) T k
ku 1 ( Σ

0.6

0.5

0.4

0.3

ku 1 ( Sbn)u 1 ( Sbn) T − u 1 (Σ)u 1 (Σ) T k
0.2

0.1

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

(b) First principal component estimation error

Fig 1: Sample size n = 100, dimension d = 100.

0.9

Error

1

S. Minsker/Estimators of the mean of a random matrix

43

Frequency

0.5

Sample covariance estimator
0.45

Robust covariance estimator
0.4
0.35
0.3

Robus t es timator er r or
b n − Σk/kΣk
kΣ

0.25
0.2
0.15

Sample covar iance er r or
k Sbn − Σk/kΣk

0.1
0.05
0
0

1

10

20

30

40

50

60

Error

(a) Covariance matrix estimation error

Frequency

1

Sample covariance estimator
0.9

Robust covariance estimator
0.8
0.7
0.6
0.5

b n)u( Σ
b n) T − u 1 (Σ)u 1 (Σ) T k
ku 1 ( Σ

0.4
0.3

ku 1 ( Sbn)u 1 ( Sbn) T − u 1 (Σ)u 1 (Σ) T k
0.2
0.1
0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

(b) First principal component estimation error

Fig 2: Sample size n = 1000, dimension d = 100.

0.9

Error

1

