arXiv:1401.3424v4 [] 1 Apr 2016

Marčenko-Pastur law for Tyler’s M-estimator
Teng Zhang
Department of Mathematics,
University of Central Florida
Orlando, Florida 32816, USA
teng.zhang@ucf.edu
Xiuyuan Cheng
Applied Mathematics Program,
Yale University, New Haven, CT 06511, USA
xiuyuan.cheng@yale.edu
Amit Singer
Department of Mathematics and PACM,
Princeton University
Princeton, New Jersey 08544, USA
amits@math.princeton.edu
April 4, 2016
Abstract
This paper studies the limiting behavior of Tyler’s M-estimator for
the scatter matrix, in the regime that the number of samples n and their
dimension p both go to infinity, and p/n converges to a constant y with
0 < y < 1. We prove that when the data samples x1 , . . . , xn are identically and independently generated from the Gaussian distribution N (0, I),
¯
the operator norm
Pn of the⊤ difference between a properly scaled Tyler’s Mestimator and i=1 xi xi /n tends to zero. As a result, the spectral distribution of Tyler’s M-estimator converges weakly to the Marčenko-Pastur
distribution.

1

Introduction

Many statistical estimators and signal processing algorithms require the estimation of the covariance matrix of the data samples. When the underlying
distribution of the data samples x1 , . . . , xn ∈ Rp is assumed to have
Pnzero mean,
a commonly used estimator is the sample covariance matrix Sn = i=1 xi x⊤
i /n.
1

However, the estimator Sn is sensitive to outliers, and performs poorly in
terms of statistical efficiency (i.e., it has a large variance) for heavy-tailed distributions, e.g., when the tail decays slower than the Gaussian tail.
A popular robust covariance estimator is an M-estimator introduced by
Tyler [20], denoted by Σ̂, which is the unique solution to
n

Σ̂ =

p X xi x⊤
i
, tr(Σ̂) = 1.
−1 x
n i=1 x⊤
i
i Σ̂

(1)

Tyler’s M-estimator gives the “shape” of the covariance, but is missing its magnitude. However, for many applications the “shape” of the covariance suffices,
for example, the principal components can be obtained from the “shape”.
Compared with the sample covariance estimator, Tyler’s M-estimator is more
robust to heavy-tailed elliptical distributions. The density function of elliptical
distributions in Rp takes the form
f (x; Σ, µ) = |Σ|−1/2 g{(x − µ)⊤ Σ−1 (x − µ)},
R∞
where g is some nonnegative function such that 0 xp−1 g(x) dx is finite. This
family of distributions is a natural generalization of the Gaussian distribution
by allowing heavier or lighter tails while maintaining the elliptical geometry of
the equidensity contours. Elliptical distributions are considered important in
portfolio theory and financial data, and we refer to the work by El Karoui [11,
Section 4] for further discussion. Besides, elliptical distributions are used by
Ollila and Tyler [19] in modeling radar data, where the empirical distributions
are heavy-tailed because of outliers.
Tyler [20] showed that when a data set follows an unknown elliptical distribution (with mean zero), Tyler’s M-estimator is the most robust covariance
estimator in the sense of minimizing the maximum asymptotic variance. This
property suggests that Tyler’s M-estimator should be more accurate than the
sample covariance estimator for elliptically distributed data. Empirically, it has
been shown to outperform the sample covariance estimator in applications such
as finance in the work by Frahm and Jaekel [13], anomaly detection in wireless sensor networks by Chen et al. [4], antenna array processing by Ollila and
Koivune [18], and radar detection by Ollila and Tyler [19].

1.1

Asymptotic analysis in a high-dimensional setting

Many scientific domains customarily deal with sets of high dimensional data
samples, and therefore it is increasingly common to work with data sets where
the number of variables, p, is of the same order of magnitude as the number
of observations, n. Under this high-dimensional setting, the asymptotic spectral properties of Sn at the limit of infinite number of samples and infinite
dimensions have been well studied by Johnstone [15]. A noticeable example
is the convergence of the spectral distribution. Denoting the eigenvalues of a
matrix A by λ1 (A), . . . , λn (A), its spectral distribution is a discrete probability
2

measure

n

1X
P = P (·|A) =
δλ (A)
n i=1 i

with δs denoting Dirac measure at s ∈ R. Marčenko and Pastur [16] showed that
when the entries of {xi }ni=1 are Gaussian independent identically distributed
random variables with mean 0 and variance 1, p, n → ∞ and p/n → y, where
0 < y ≤ 1, the spectral distribution of the eigenvalues of Sn converges weakly
to the Marčenko-Pastur distribution defined by
p
1 y (y+ − x)(x − y− )
√
ρMP,y (x) =
1[y− ,y+ ] , where y± = (1 ± y)2 .
(2)
2π
x
Tyler’s M-estimator is closely related to and can be considered as a special
case of Maronna’s M-estimator, which is defined by
n

Σ̄ =

1X
−1
u(x⊤
xi )xi x⊤
i Σ̄
i
n i=1

(3)

for a nonnegative function u : [0, ∞) → [0, ∞). The properties of Maronna’s
M-estimator in the high-dimensional regime when p, n → ∞, p/n → y and
0 < y < 1 have been analyzed in recent works by Couillet et al. [7, 8], which obtained convergence results for a properly scaled Maronna’s M-estimator under
the assumptions that u(x) is nonnegative, nonincreasing and continuous; xu(x)
is nondecreasing and bounded and supx xu(x) > 1. Moreover, spiked random
matrix models were also studied by Couillet [5]. However, these results do not
apply to Tyler’s M-estimator, although Frahm and Jaekel [13] have conjectured
that the spectral distribution converges weakly to the Marčenko-Pastur distribution. Some works focused on the performance of Tyler’s M-estimator for the
case p, n → ∞ and p/n → 0: p
Dümbgen p
[10] showed that the condition number
p/n), and Frahm and Glombek [12]
of Tyler’s estimator is 1 + 4 p/n + o( p
showed that the spectral distribution of n/p(Σ̄ − I) converges weakly to a
semicircle distribution.

1.2

Main results

In this paper, we analyze Tyler’s M-estimator in the high-dimensional setting.
Our main results, Theorem 2.3 and Corollary 2.4, show that as p, n → ∞ and
p/n → y, 0 < y < 1, the spectral distribution of a properly scaled Tyler’s
M-estimator converges weakly to the Marčenko-Pastur distribution ρMP,y (x).
Based on the properties of Tyler’s M-estimator, this paper analyzes the spectral
distribution when data samples are i.i.d. drawn from other distributions, such
as elliptical distributions.
When data samples are generated from elliptical distributions, the spectral distribution of the sample covariance estimator has been studied by El
Karoui [11, Theorem 2]. Compared to Corollary 2.4, the limiting spectral distribution of Sn is much more complicated, and therefore our result might be
more applicable in practice.
3

High-dimensional analysis of Maronna’s M-estimator of the covariance are
generally obtained by showing that the operator norm of the difference between
M-estimator and a standard Wishart matrix (or sample covariance matrix)
tends to 0: Dümbgen [10] proved it by a linear expansion of the M-estimator,
and Couillet et al. [7, 8] proved it by representing Maronna’s M-estimator as a
weighted sum of xi x⊤
i and prove the uniform convergence of the weights. We
follow the same direction while giving an alternate proof for the convergence of
the weights, by considering the weights as the solution to an optimization problem, which can handle Tyler’s M-estimator that is not covered by the results
in Couillet et al. [7, 8]. We remark that this approach can also be applied to
Maronna’s M-estimator to prove some of the results in Couillet et al. [7, 8].
The rest of the paper is organized as follows. In Section 2 we introduce
the representation of Tyler’s M-estimator as a linear combination of xi x⊤
i and
present the main result that when the data set is i.i.d. sampled from the Gaussian distribution N (0, I), a properly scaled Tyler’s is asymptotically equivalent
¯
to Sn in the sense that kpΣ̂ − Sn k → 0. As a result, the spectral distribution
of Tyler’s M-estimator converges weakly to the Marčenko-Pastur distribution.
We also extend the result to elliptical distributions. The technical proofs are
given in Section 3. While some Lemmas and technical proofs are also used by
Couillet et al. [7, 8] (for example, Lemma 3.2 and the analysis in the proof of
Theorem 2.3 are similar to Couillet et al. [7, Lemma 2, Theorem 1] and Couillet
et al. [8, Lemma 6], we still include them for the completeness of the paper.
As for notations, we use c, c′ , C, C ′ to denote any fixed constants as p, n → ∞
(though they may depend on y). Depending on the context, they might denote
different values in different equations.

2

Tyler’s M-estimator in the High-dimensional
Regime

We introduce the representation of Tyler’s M-estimator as a linear combination
of xi x⊤
i in Section 2.1, and present the main result in Section 2.2 that when
the data set is i.i.d. sampled from the Gaussian distribution N (0, I), kpΣ̂ − Sn k
¯
converges to 0 almost surely. Based on this observation, we prove that the
spectral distribution of Tyler’s M-estimator converges weakly the MarčenkoPastur distribution in Section 2.3. The generalization of the results to more
general settings is also discussed in Section 2.3.

2.1

Properties of Tyler’s M-estimator

The analysis for Tyler’s M-estimator in this paper is based on the following representation, whose proof is deferred to Section 3. We remark that equation (5) in
Lemma 2.1 has appeared in the work by Wiesel [21, (27)] and Hardt and Moitra
[14, Section A] as “covariance estimation in scaled Gaussian distributions” and
“Barthe’s convex program”, but its connection to Tyler’s M-estimator has not
been rigorously justified yet.
4

Lemma 2.1. Tyler’s M-estimator can be written as
Σ̂ =

n
X
i=1

n

. X
ŵi xi x⊤
ŵi xi x⊤
tr
i ,
i

(4)

i=1

where {ŵi }ni=1 are uniquely defined by
(ŵ1 , . . . , ŵn ) =

arg
min
P

wi >0,

2.2

n
i=1

wi =1

−

n
X

ln wi +

i=1

n
X

n
wi xi x⊤
ln det
i .
p
i=1

(5)

Isotropic Gaussian Distribution

In this subsection, we assume that {xi }ni=1 ⊂ Rp are i.i.d. drawn from N (0, I).
¯
The main result, Theorem 2.3, characterizes the convergence and convergence
rate of Tyler’s M-estimator to Sn in terms of the operator norm. Its proof
applies Lemma 2.2, whose proof is rather technical and deferred to Section 3.
Tyler’s M-estimator does not exist when p > n (see the argument by Zhang [22,
Theorem III.1]) and it is not unique when p = n (one may check that when
xi = ei for all 1 ≤ i ≤ p, all diagonal matrices with trace 1 satisfy (1)). As a
result, throughout the paper we assume y < 1.
Lemma 2.2. If {xi }ni=1 are i.i.d. sampled from N (0, I), then max1≤i≤n |n ŵi −
¯
1| converges to 0 almost surely as p, n → ∞: There exist C, c, c′ > 0 such that
′
for any ε < c ,


2
Pr max |n ŵi − 1| ≤ ε ≥ 1 − Cne−cε n .
(6)
1≤i≤n

Theorem 2.3. Suppose that {xi }ni=1 are i.i.d. sampled from N (0, I), p, n → ∞
¯
and p/n = y, where 0 < y < 1, then the operator norm of the difference between
Sn and a scaled Tyler’s M-estimator converges to 0 almost surely, and there
exist C, c, c′ > 0 such that for any ε < c′ ,
!
n
2
1X
⊤
xi xi ≤ ε ≥ 1 − Cne−cε n .
(7)
Pr
p Σ̂ −
n i=1
Theorem 2.3 implies that all first order properties of the sample covariance
matrix extend to Tyler’s estimator. The strategy of the proof for Theorem 2.3
is as follows. According to Lemma 2.1, a scaled Tyler’s
is a linear
Pn M-estimator
⊤
combination of xi x⊤
,
i.e.,
it
can
be
written
as
ŵ
x
x
(up
to a scali
i=1 i i i
ing). Then Lemma 2.2 shows that nŵi converges to 1 uniformly. Based on the
following matrix analysis, Theorem 2.3 is concluded.
Proof of Theorem 2.3. We first prove that for ε < c′ ,
!
n
n
X
2
1X
⊤
⊤
Pr
ŵi xi xi −
xi xi ≤ ε ≥ 1 − Cne−cε n .
n
i=1
i=1
5

(8)

Let Bn =

Pn

i=1 (ŵi

− 1/n)xi x⊤
i =

kBn k = sup v ⊤ Bn v = sup
kvk=1

≤ sup

Pn

n
X

kvk=1 i=1

n
X

kvk=1 i=1

ŵ −

1
1
n

ŵi xi x⊤
i −

i=1

∞

(ŵi −

Pn

i=1

xi x⊤
i /n, then

1
)(v ⊤ xi )2
n
n

(v ⊤ xi )2 ≤ knŵ − 1k∞ k

1X
xi x⊤
i k.
n i=1

Since knŵ − 1k∞ → 0 with probability P
estimated in (6), and Davidson and
n
Szarek [9, Theorem II.13] showed that k i=1 xi x⊤
i /nk is bounded above by
√ 2
(1 + 2 y) with probability 1 − C exp(−cn), (8) is proved.
Second, since
n
X

ŵi xi x⊤
≤
i

i=1

Pr k

n
X
i=1

n
X
i=1

ŵi xi x⊤
i −

n
X

xi x⊤
i /n +

′
ŵi xi x⊤
i k< C

xi x⊤
i /n ,

i=1

i=1

!

n
X

> 1 − Cn exp(−cn).

(9)

Pn
Pn
⊤
Besides, tr( i=1 ŵi xi x⊤
i ) =
i=1 ŵi xi xi → p in the same rate as in (9):
applying the concentration of high-dimensional Gaussian measure on the sphere
by Barvinok [2, Corollary 2.3], we have
)#
)
( n
" ( n
X
X
⊤
⊤
ŵi xi xi > p/(1 − ε)
(10)
ŵi xi xi < p(1 − ε) , Pr
max Pr
i=1

i=1



≤ max Pr






2
min kxi k2 < p(1 − ε) , Pr max kxi k2 > p/(1 − ε)
< ne−ε p/4 .

1≤i≤n

1≤i≤n

Combining (9), (10) and (4),
n
X
i=1

ŵi xi x⊤
i − p Σ̂ =

n
X

ŵi xi x⊤
i

i=1

(

)
n
X
ŵi xi x⊤
1 − p/tr(
i )

(11)

i=1

converges in the same rate as specified in (8). (7) is then proved by combining
(8), (11) and the triangle inequality.
p
From the probabilistic estimation (7) we obtain a convergence rate of O( ln n/n).
The logarithmic factor is due to a “max” bound of {ŵi }ni=1 in Lemma 2.2, while
in
√ fact, an “average” bound is expected. As a result, we conjecture that this
ln n factor could be possibly removed by a more rigorous argument.

2.3

More general distributions and spectral distribution

We remark that Theorem 2.3 can be extended from the setting of the normal distribution N (0, I) to any elliptical distribution µp , which is characterized
¯
by the probability density function µp (x) = C(gp ) det(Tp )−1/2 gp (x⊤ Tp−1 x),
6

Tp is a positive definite matrix in Rp×p , gp : [0, ∞) → [0, ∞) satisfies
Rwhere
∞
p−1
< ∞, and C(gp ) is a normalization parameter that only depends
0 gp (x)x
Pn
on gp . Then ktr(Tp )Σ̂ − i=1 xi x⊤
i /nk → 0 almost surely as p, n → ∞. The
analysis is based on Theorem 2.3, the affine equivariance property of Tyler’s
M-estimator, and the fact that Tyler’s M-estimator is unchanged if {xi }ni=1 are
replaced by {ci xi }ni=1 .
Another direction of generalization of Theorem 2.3 is the model by Couillet et
al. [7]: The elements of {xi }ni=1 are i.i.d. sampled from an either real or circularly
symmetric complex distribution with E(xij ) = 0, E(x2ij ) = 1, and E(|xij |8+η ) <
α for some η, α > 0. Then, following the proof in this paper (while replacing
Pn
Lemma 3.2 by [7, Lemma 2]), one can show that p Σ̂ − i=1 xi x⊤
i /n → 0
almost surely as p, n → ∞.
We have the following results on the weak convergence of the spectral distribution of Tyler’s M-estimator, where the first part proves the conjecture by
Frahm and Jaekel [13].
Corollary 2.4.
• If {xi }ni=1 are i.i.d. sampled from N (0, I), then the spec¯
tral measure P (·|pΣ̂) converges weakly to the Marčenko-Pastur distribution.
• If {xi }ni=1 are i.i.d. sampled from an elliptical distribution C(gp )gp (x⊤ Tp−1 x)
such that the spectral measure of Tp converges weakly to a distribution
H on R. Then the spectral measure P (·|tr(Tp )Σ̂) converges
weakly to a
R
probabilistic measure ρ whose Stieltjes transform s(z) = 1/(x − z)ρ( dx)
(z ∈ C \ R) is given implicitly by
Z
1
s(z) =
dH(t).
t{1 − y − yz s(z)} − z
P
This corollary can be proved by combining ktr(Tp )Σ̂ − ni=1 xi x⊤
i /nk → 0,
the analysis on the perturbation
of
eigenvalues
by
Bhatia
[3,
Corollary
III.4.2],
Pn
/n
by
Marčenko
and
Pastur
[16],
Bai
and
the spectral measure of i=1 xi x⊤
i
Silverstein [16, 1, (6.1.2)] and Slutsky’s Lemma.

3
3.1

Proof of Lemmas
Proof of Lemma 2.1

We start with the definition
(ẑ1 , . . . , ẑn ) = P
arg min ln det
n
i=1

and
Σ̂z =

zi =1

n
X

n
X

eẑi xi x⊤
i .

i=1

7

i=1

e zi xi x⊤
i



(12)

(13)

The solution to (12) is unique, which follows from the convexity of the
objective function (see Wiesel [21, Lemma 4]). Besides,
noticing the equivaP
lence between (12) and (5) (by plugging wi = ezi /( ni=1 ezi ) and zi = ln wi −
Pn
( i=1 ln wi − 1)/n, there exists c1 > 0 such that Σ̂z = c1 Σ̂.
Next we will prove that Σ̂z satisfies
n
X
i=1

xi x⊤
i
−1
x⊤
i Σ̂z xi

= cΣ̂z , for some c > 0.

(14)

By checking the directional
derivative of the objective function in (12), for
Pn
any (δ1 , . . . , δn ) with i=1 δi = 0,
n
X

−1
δi eẑi x⊤
i Σ̂z xi = 0.

i=1

Therefore, there exists c2 such that
−1
eẑi x⊤
i Σ̂z xi = c2 , for all 1 ≤ i ≤ n.

(15)

Therefore, (14) is proved by applying (15) and (13):
n
X
i=1

xi x⊤
i
−1
x⊤
i Σ̂z xi

=

n
X

eẑi xi x⊤
i /c2 = Σ̂z /c2 ,

i=1

Since Σ̂z = c1 Σ̂, (14) also holds when Σ̂z is replaced by Σ̂:
n
X
i=1

xi x⊤
i
−1 x
x⊤
i
i Σ̂

= cΣ̂, for some c > 0.

(16)

At last, we will prove that Σ̂ satisfies the definition of Tyler’s M-estimator
in (1), that is, the constant c in (16) is given by c = n/p. For the objective
function
n
X
−1
ln(x⊤
xi ) + c ln det(Σ),
F (Σ) =
i Σ
i=1

its derivative with respect to Σ−1 is given by
n
X
i=1

⊤ −1
xi )−1 xi − cΣ.
x⊤
i (xi Σ

Therefore, Σ̂ is a stationary point of F (Σ). Since F (Σ) is geodesically convex
(argument follows directly from Wiesel [21] and Zhang [22]), Σ̂ is the global
minimizer of F (Σ).
the minimizer of F (Σ) exists only when c = n/p. Since F (aI) =
PnHowever,
⊤
i=1 ln(xi xi ) − n ln a + c p ln a, we have
(
as a → 0, if c > n/p
F (aI) → −∞
.
as a → ∞, if c < n/p
Therefore, the constant c in (16) is given by c = n/p, and Lemma 2.1 is
proved.
8

3.2

Proof of Lemma 2.2

We start with an outline of the proof, which consists of three parts. First, we
rewrite the constrained optimization problem (5) to the problem of finding the
root of g(w), which will be defined in (17). Since the root of g(w) is nŵ − 1, we
only need to show the convergence of the root of g(w). Second, we will show
that g(0) converges to 0, ∇g(0̄) is large and the variation of ∇g(w) is bounded.
¯
¯
Finally, we will use a perturbation analysis and the observations on g(0) and
¯
∇g(w) to show that the root of g(w) converges to 0.
¯
The proof depends on Lemma 3.2, Lemma 3.3 and Lemma 3.1, and their
proofs are postponed to subsequent sections.
Lemma 3.1. For a function f (w) : Rp → Rp , assume that ∇f (0) = I, and
¯
k∇f (w) − ∇f (0)k∞ = maxi≤i≤p k∇fi (w) − ∇fi (0)k∞ < C5 kwk∞ for kwk∞ ≤
¯
¯
1, and kf (0)k∞ < min(1/9C5 , 1/3). Then there exists w̃ such that kw̃k∞ <
¯
3kf (0)k∞ and f (w̃) = 0.
¯
P
Lemma 3.2. If xi ∼ N (0, I) for all 1 ≤ i ≤ n, and S = ni=1 xi x⊤
i /n, then
¯
there exists c, C, c′ > 0 such that for any ε < c′ ,


2
1
−1
S
x
−
1|
<
ε
≥ 1 − Cne−cε n .
Pr max | x⊤
i
i
1≤i≤n p
−1
Lemma 3.3. For the n × n matrix A defined by Aij = (x⊤
xj )2 /(n p), (a)
i S
kAk∞ < 2 with probability 1 − Cn exp(−cn).
(b) There exists c = c(p, n) > 0 and C2 = C2 (y) > 0 such that k(I − A +
c11⊤ )−1 k∞ < C2 with probability 1 − Cn exp(−cn).

We start the first part of the proof with the construction of g(w). We let
g(w) = ∇G(w + 1),

(17)

where
G(w) = −

n
X
i=1

ln wi +

n
n
X
n
c0 X
wi xi x⊤
)
+
wi − n)2 ,
ln det(
(
i
p
2
i=1
i=1

(18)

and the constant c0 will be specified later before (34).
It is easyP
to prove that the minimizer P
of G(w) and the zeros of ∇G(w)
n
n
must satisfy i=1 wi = n (otherwise nw/( i=1 wi ) is a better minimizer and
∇G(w)
is nonzero). Therefore
minimizing (18) is equivalent
P
P
Pn to minimizing
− ni=1 ln wi + n/p · ln det( ni=1 wi xi x⊤
)
with
constraint
i
i=1 wi = n, which
is the
same
as
(5)
except
for
the
constraint.
Noticing
that
a
scaling
of w increases
Pn
Pn
)
by
a
constant
only
depending
on the
− i=1 ln wi + n/p · ln det( i=1 wi xi x⊤
i
scale, the minimizer of (18) is unique and it is nŵ, where ŵ is defined in (5).
By the convexity of its equivalent problem (12), the root of g(w) is also unique
and it is nŵ − 1.

9

For the second part of the proof, we start by proving that g(0) is small. By
¯
calculation, the i-th component of function g(w) is
gi (w) = −

n
n
X
X
1
n
−1
wi .
wi xi x⊤
x i + c0
+ x⊤
i )
i (nS +
wi + 1
p
i=1
i=1

Applying Lemma 3.2,
2

Pr (kg(0)k∞ < ε) ≥ 1 − Cne−cε n .
¯

(19)

Now we will prove that ∇g(0̄) is bounded from below. By calculation, its (i, j)th entry is

∇g(w)

i,j

= I(i = j)

n
−1 o2
X
1
n n ⊤
⊤
x j + c0 .
w
x
x
nS
+
x
−
i
i
i
i
(wi + 1)2
p
i=1

Applying Lemma 3.3,
k{∇g(0̄)}−1 k∞ < C2 with probability 1 − Cne−cn .

(20)

Now we bound the variation of ∇g(w) in the region kwk∞ < 1/2. Apply
|1/(wi + 1)2 − 1| < 3|wi − 1| ≤ 3kwk∞ and coordinatewise comparison,
n
|∇i,j g(w) − ∇i,j g(0)| ≤ I(i = j) (3kwk∞ ) + 3kwk∞ · |Aij |.
¯
p
Therefore, the variation of ∇g(w) is bounded by
k∇g(w) − ∇g(0̄)k∞ < (3 + 3nkAk∞ /p)kwk∞ .

(21)

At last we finish the third part of the proof of Lemma 2.2 by applying
Lemma 3.1 to f (w) = {∇g(0̄)}−1 g(w/2). It is easy to verify that ∇f (0) = I.
¯
Due to (19) and (20), kf (0)k∞ ≤ k(∇g(0̄))−1 k∞ kg(0̄)k∞ → 0 in the same rate as
¯
−cn
in (19) and kf (0)k∞ < min(1/9C5 , 1/3) holds with probability 1−Cne
. Due
¯
to (19), (21), and the boundedness of kAk∞ (Lemma 3.3), k∇f (w) − ∇(0̄)k∞ <
C5 kwk∞ also holds with probability 1 − Cne−cn . Therefore the assumption in
Lemma 3.1 holds with probability 1 − Cne−cn and there exists w̃ such that
f (w̃) = 0 and
kw̃k∞ < 3kf (0)k∞ .
(22)
¯
When f (w̃) = 0, we have g(2w̃) = 0 and by previous discussion 2w̃ = nŵ−1.
therefore (22) gives
knŵ − 1k∞ < 6kf (0)k∞ .
¯
Since kf (0)k∞ converges to 0 in the rate as in (19), knŵ − 1k∞ converges
¯
in the same rate and Lemma 2.2 is proved.

10

3.2.1

Proof of Lemma 3.1

Proof. When kwk∞ ≤ 1,
fj (w) − fj (0) =

Z

1

t=0

ej w⊤ , ∇f (t w) dt

(23)

Z 1
ej w⊤ , ∇f (t w) − ∇f (0) + I dt = wj +
w⊤ {∇f (t w) − ∇f (0)}ej dt
¯
¯
t=0
t=0
Z 1
w⊤ {∇f (t w) − ∇f (0)}k∞ ≤ wj + C5 kwk2∞ .
≤wj + k
¯
t=0
=

Z

1

Similarly
fj (w) − fj (0) ≥ −C5 kwk2∞ + wj .

(24)

To prove it, we consider the continuous mapping h(w) = w−f (w)/(4+9C5 )
and will prove that h maps A to itself, where
A = {w : w ∈ [−3η, 3η]n } and η = kf (0)k∞ .
1. |wi | < 2η. Then apply (23) and (24) (they are applicable since for any w ∈
A, kwk∞ ≤ 1), we have |fi (w)| < |fi (0)| + C5 kwk2∞ + |wi | ≤ η + C5 (3η)2 + 3η <
(4 + 9C5 )η (η 2 < η since η < 1). Therefore, |hi (w)| ≤ |wi | + |fi (w)|/(4 + 9C5 ) ≤
3η.
2. wi > 2η, then applying (24),
fi (w) ≥ −|fi (0)| + wi − C5 kwk2∞ ≥ −η + 2η − C5 (3η)2 .
¯
Since η < 1/9C5 , we have fi (w) < 0 and therefore hi (w) ≤ wi ≤ 3η.
Similar to case 1 we can prove that hi (w) ≥ −3η. Therefore |hi (w)| < 3η.
3. Similar to case 2, when wi < −2η, |hi (w)| < 3η.
Therefore the continuous mapping h maps the convex, compact set A to
itself. By Schauder fixed point theorem, h(x) has a fixed point in A and
Lemma 3.1 is proved with w̃ being the fixed point.
3.2.2

Proof of Lemma 3.2

Assuming the SVD decomposition of X is X = U ΣV ⊤ , where U ∈ Rn×p and
U ⊤ U = I. Since xi ∼ N (0, I) for all 1 ≤ i ≤ n, U is uniformly distributed over
¯
the space of all orthogonal n × p matrices. Since
1
XS −1 X = (U ΣV ⊤ )( V Σ2 V ⊤ )−1 (U ΣV ⊤ ),
n

(25)

2
if we write the row of U by u1 , . . . , un , then n1 xi S −1 xi = u⊤
i ui = kui k .
Since U can be considered as the first p columns of a random n×n orthogonal
matrix (with haar measure over the set of all n × n orthogonal matrices), ui
can be considered as the first p entries from a random vector of length n that
is sampled from the uniform sphere in Rn .

11

Pp
Pn
Therefore, kui k2 ∼ j=1 gj2 / j=1 gj2 for i.i.d. random variables {gj }nj=1 ∼
N (0, 1). Applying the the concentration result by Barvinok [2, Corollary 2.3],
we have
)
( n
X
2
n
2
≤ e−ε n/4
(26)
gi ≥
Pr
1
−
ε
i=1

and

(

Pr

n
X
i=1

)

gi2 ≤ n(1 − ε)

2

≤ e−ε

n/4

,

(27)

therefore
)
(


p
X
p(1 − ε)2
p
p
2
2
Pr
gi ≤
≥ Pr p(1 − ε) ≤
≤ ku1 k ≤
n
n(1 − ε)2
1−ε
i=1
(
)
n
X
2
2
n
gi2 ≤
+ Pr n(1 − ε) ≤
≥ 1 − 2e−ε p/4 − 2e−ε n/4 .
1
−
ε
i=1

For ε ≤ 0.1, we have


n
p
p o
1
−1
S
x
−
1|
≤
ε
≥ 1 − n Pr |ku1 k2 − | > ε
Pr max | x⊤
i
i
1≤i≤n p
n
n



2
p
p(1 − ε/3)
≤ ku1 k2 ≤
(28)
≥1 − n 1 − Pr
n
n(1 − ε/3)2
2

≥1 − 2ne−ε

p/36

2

− 2ne−ε

n/36

,

(29)

where the second inequality follows from 1−3ε ≤ (1−ε)2 and 1/(1 − ε)2 ≤ 1+3ε.
3.2.3

Proof of Lemma 3.3
P
(a) Since kAk∞ = max1≤i≤n ( 1≤j≤n Aij ), and
X

1≤j≤n

Aij =

X

1≤j≤n

X
1 ⊤ −1
−1
−1
−1
xi S xj x⊤
xi = x⊤
(
xj x⊤
xi /np
j S
i S
j )S
np
1≤j≤n

(30)

−1
=x⊤
(nS)S −1 xi /np
i S

=

−1
x⊤
xi /p,
i S

(31)

it follows from (29) with ε = 0.1 that kAk∞ < 2 holds with probability 1 −
Cn exp(−cn).
(b) We first prove that there exists C3 = C3 (y) such that
kA − c0 11⊤ k∞ ≤ C3 < 1 with probability 1 − Cn exp(−cn).

(32)

We start with the proof of (32) with another lemma:
Lemma 3.4. There exists a c4 > 0 such that with probability 1 − C exp(−cn),
n
X

√
I(x⊤
1 xj > c4 p) > 0.75n.

j=1

12

Davidson and S. Szarek [9, Theorem II.13] showed that There exists C4 =
−1
C4 (y) such that kSk < C4 with probability 1−Cn exp(−cn). Therefore x⊤
xj ≥
i S
⊤
xi xj /C4 and Lemma 3.4 implies that for any 1 ≤ i ≤ n:
n
X
j=1

√
−1
I(x⊤
xj > c4 p/C4 ) > 0.75 with probability 1 − C exp(−cn). (33)
i S

Let c0 = (c4 /C4 )2 /n, then (33) implies
X
X
−1
|Ai,j − c| ≤
|Ai,j | − 0.25c n ≤ x⊤
xi /p − 0.25(c4 /C4 )2 , (34)
i S

1≤j≤n

1≤j≤n

where the last step follows from (31).
−1
Applying the estimation of x⊤
xi /p in (29) and a union bound argument
i S
over all 1 ≤ i ≤ n to (34), (32) is proved for C3 = 1 + η − 0.25(c4 /C4 )2 .
Lemma 3.3(b) follows from (32) with C2 = 1/(1 − C3 ), where the expansion
of (I − A + c11⊤ )−1 exists since kA + c11⊤ k ≤ kA + c11⊤ k∞ < 1. Applying
kB1 B2 k∞ ≤ kB1 k∞ kB2 k∞ , we have
k(I−A+c11⊤ )−1 k∞ = k

3.2.4

∞
X

k=0

(c11⊤ −A)k k∞ ≤

∞
X

k=0

kc11⊤ −Akk∞ ≤

∞
X

k=0

C3i =

1
.
1 − C3
(35)

Proof of Lemma 3.4

We first show that there exists c4 such that for all p,
√
E{I(|x⊤
1 x2 | > c4 p)} ≥ 0.85.

(36)

WLOG we rotate x1 such that it is nonzero only at the first coordinate, and
x2 = (g1 , g2 , ..., gp ) where gi ∼ N (0, 1). Then |x⊤
1 x2 | = |g1 | kx1 k.
Notice that kx1 k2 is the sum of p√independent χ21 distribution and Eχ21 = 1,
1 − Ce−cn . Besides,
by central√limit theorem, kx1 k ≤ 2p with probability
√
−1
Pr(|g1 | > 2 c4 ) ≥ 0.85 for c4 = Φ (1 − 0.85/2)/ 2. Therefore (36) is proved
by combining the estimations on |g1 |, x1 and |x⊤
1 x2 | = |g1 | kx1 k.
To obtain Lemma 3.4 from (36), we apply Hoeffding’s inequality to the
√
indicator function I(|x⊤
i xj | > c4 p) over all 1 ≤ j ≤ n, j 6= i.

4

Summary

We showed that Tyler’s M-estimator is asymptotically equivalent to Sn in the
sense that kpΣ̂ − Sn k → 0 as p, n → ∞ and p/n → y, where 0 < y < 1 and
data samples follow the distribution of N (0, I). We also proved the conjecture
¯
that the spectral distribution of Tyler’s M-estimator converges weakly to the
Marčenko-Pastur distribution, and extended the results to elliptical distributions.
13

There are several possible future directions of this work. First, it would be
interesting to analyze the second order statistics of Tyler’s M-estimator, considering that Couillet [6] has already investigated Maronna’s M-estimators. Second, we would like to theoretically quantify the behavior of Tyler’s M-estimator
in the spiked covariance model by Couillet [5], which includes the analysis of
the distribution of the top eigenvalue for the null cases and the analysis of the
non-null case. A recent work by Morales-Jimenez et al. [17] on the non-null case
introduced a mixture model that consists of a Gaussian distribution and some
deterministic or random outliers, and analyzed the performance of Maronna’s
M-estimator. Analyzing the performance of Tyler’s M-estimator in this model
would be another possible future direction.

Acknowledgements
A. Singer was partially supported by Award Number FA9550-12-1-0317 and
FA9550-13-1-0076 from AFOSR, by Award Number R01GM090200 from the
NIGMS, and by Award Number LTR DTD 06-05-2012 from the Simons Foundation.

References
[1] Z. Bai and J. Silverstein. Spectral Analysis of Large Dimensional Random
Matrices. Springer series in statistics. Springer, 2009.
[2] A. Barvinok. Math 710: Measure concentration. Lecture notes, Department
of Mathematics, University of Michigan, 2005.
[3] R. Bhatia. Matrix Analysis. Number 169 in Graduate Texts in Mathematics. Springer, New York, 1997.
[4] Y. Chen, A. Wiesel, and A. Hero. Robust shrinkage estimation of highdimensional covariance matrices. Signal Processing, IEEE Transactions on,
59(9):4097 –4107, sept. 2011.
[5] R. Couillet. Robust spiked random matrices and a robust g-music estimator. Journal of Multivariate Analysis, 140:139 – 161, 2015.
[6] R. Couillet, A. Kammoun, and F. Pascal. Second order statistics of robust
estimators of scatter. application to GLRT detection for elliptical signals.
Journal of Multivariate Analysis, 143:249 – 274, 2016.
[7] R. Couillet, F. Pascal, and J. Silverstein. Robust estimates of covariance
matrices in the large dimensional regime. Information Theory, IEEE Transactions on, 60(11):7269–7278, Nov 2014.
[8] R. Couillet, F. Pascal, and J. W. Silverstein. The random matrix regime
of Maronna’s M-estimator with elliptically distributed samples. Journal of
Multivariate Analysis, 139:56 – 78, 2015.
14

[9] K. Davidson and S. Szarek. Local operator theory, random matrices and
banach spaces. Handbook of the Geometry of Banach Spaces, 1:317, 2001.
[10] L. Dümbgen. On Tyler’s M-functional of scatter in high dimension. Annals
of the Institute of Statistical Mathematics, 50(3):471–491, 1998.
[11] N. El Karoui. Concentration of measure and spectra of random matrices:
Applications to correlation matrices, elliptical distributions and beyond.
The Annals of Applied Probability, 19(6):2362–2405, 12 2009.
[12] G. Frahm and K. Glombek. Semicircle law of Tyler’s M-estimator for
scatter. Statistics & Probability Letters, 82(5):959 – 964, 2012.
[13] G. Frahm and U. Jaekel. Tyler’s M-estimator, random matrix theory, and
generalized elliptical distributions with applications to finance. Technical
report, 2007.
[14] M. Hardt and A. Moitra. Algorithms and hardness for robust subspace
recovery. In S. Shalev-Shwartz and I. Steinwart, editors, COLT, volume 30
of JMLR Proceedings, pages 354–375. JMLR.org, 2013.
[15] I. M. Johnstone. High dimensional statistical inference and random matrices. In International Congress of Mathematicians. Vol. I, pages 307–333.
Eur. Math. Soc., Zürich, 2007.
[16] V. A. Marčenko and L. A. Pastur. Distribution of eigenvalues for some sets
of random matrices. Mathematics of the USSR-Sbornik, 1(4):457, 1967.
[17] D. Morales-Jimenez, R. Couillet, and M. McKay. Large dimensional analysis of robust M-estimators of covariance with outliers. Signal Processing,
IEEE Transactions on, 63(21):5784–5797, Nov 2015.
[18] E. Ollila and V. Koivunen. Robust antenna array processing using Mestimators of pseudo-covariance. In Personal, Indoor and Mobile Radio
Communications, 2003. PIMRC 2003. 14th IEEE Proceedings on, volume 3, pages 2659–2663 vol.3, 2003.
[19] E. Ollila and D. Tyler. Distribution-free detection under complex elliptically symmetric clutter distribution. In Sensor Array and Multichannel
Signal Processing Workshop (SAM), 2012 IEEE 7th, pages 413–416, June
2012.
[20] D. E. Tyler. A distribution-free M-estimator of multivariate scatter. The
Annals of Statistics, 15(1):pp. 234–251, 1987.
[21] A. Wiesel. Geodesic convexity and covariance estimation. Signal Processing, IEEE Transactions on, 60(12):6182 –6189, dec. 2012.
[22] T. Zhang. Robust subspace recovery by Tyler’s M-estimator. Information
and Inference, 2015.

15

