arXiv:1612.02099v1 [math.ST] 7 Dec 2016

Statistical and Computational Guarantees of Lloyd’s Algorithm
and Its Variants
Yu Lu1 and Harrison H. Zhou1
1

Yale University

December 8, 2016
Abstract
Clustering is a fundamental problem in statistics and machine learning. Lloyd’s algorithm, proposed in 1957, is still possibly the most widely used clustering algorithm
in practice due to its simplicity and empirical performance. However, there has been
little theoretical investigation on the statistical and computational guarantees of Lloyd’s
algorithm. This paper is an attempt to bridge this gap between practice and theory. We
investigate the performance of Lloyd’s algorithm on clustering sub-Gaussian mixtures.
Under an appropriate initialization for labels or centers, we show that Lloyd’s algorithm
converges to an exponentially small clustering error after an order of log n iterations,
where n is the sample size. The error rate is shown to be minimax optimal. For the
two-mixture case, we only require the initializer to be slightly better than random guess.
In addition, we extend the Lloyd’s algorithm and its analysis to community detection
and crowdsourcing, two problems that have received a lot of attention recently in statistics and machine learning. Two variants of Lloyd’s algorithm are proposed respectively
for community detection and crowdsourcing. On the theoretical side, we provide statistical and computational guarantees of the two algorithms, and the results improve upon
some previous signal-to-noise ratio conditions in literature for both problems. Experimental results on simulated and real data sets demonstrate competitive performance of
our algorithms to the state-of-the-art methods.

1

Introduction

Lloyd’s algorithm, proposed in 1957 by Stuart Lloyd at Bell Labs [40], is still one of the most
popular clustering algorithms used by practitioners, with a wide range of applications from
computer vision [3], to astronomy [45] and to biology [26]. Although considerable innovations
have been made on developing new provable and efficient clustering algorithms in the past
six decades, Lloyd’s algorithm has been consistently listed as one of the top ten data mining
algorithms in several recent surveys [55].
Lloyd’s algorithm is very simple and easy to implement. It starts with an initial estimate
of centers or labels and then iteratively updates the labels and the centers until convergence.
1

Despite its simplicity and a wide range of successful applications, surprisingly, there is little
theoretical analysis on explaining the effectiveness of Lloyd’s algorithm. It is well known that
there are two issues with Lloyd’s algorithm under the worst case analysis. First, as a greedy
algorithm, Lloyd’s algorithm is only guaranteed to converge to a local minimum [43]. The
k-means objective function that Lloyd’s algorithm attempts to minimize is NP-hard [14, 41].
Second, the convergence rate of Lloyd’s algorithm can be very slow. Arthur and Vassilvitskii
[6] construct a worst-case showing that Lloyd’s algorithm can require a superpolynomial
running time.
A main goal of this paper is trying to bridge this gap between theory and practice of
Lloyd’s algorithm. We analyze its performance on the Gaussian mixture model [47, 52], a
standard model for clustering, and consider the generalization to sub-Gaussian mixtures,
which includes binary observations as a special case. Specifically, we attempt to address
following questions to help understand Lloyd’s algorithm: How good does the initializer need
to be? How fast does the algorithm converge? What separation conditions do we need? What
is the clustering error rate and how it is compared with the optimal statistical accuracy?
Despite the popularity of Lloyd’s algorithm as a standard procedure for the k-means
problem, to the best of our knowledge, there is little work in statistics to understand the
algorithm. Some efforts have been made by computer scientists to develop effective initialization techniques for Lloyd’s algorithm [7, 46, 4]. Their main focus is to find polynomial-time
approximation scheme of the k-means objective function rather than to identify the cluster labels of data points, which is often the primary interest for many applications. It is
worthwhile to emphasize that when the signal-to-noise ratio is not sufficiently large, a small
k-means objective function value does not necessarily guarantee a small clustering error. Furthermore, the error rate is different from the optimal error of an exponential form. Recently,
some fascinating results by Kumar and Kannan [37] and Awasthi and Sheffet [8] show that
under certain strong separation conditions, the Lloyd’s algorithm initialized by spectral clustering correctly classifies all data points with high probability. However, they focus on the
strong consistency results and the clustering error rate of the Lloyd’s algorithm remains unclear. It is desirable to have a systematic study of Lloyd’s algorithm under various separation
conditions such that the strong consistency can be included as a special case.
Lloyd’s algorithm is an iterative procedure. Its analysis can be challenging due to dependence between iterative steps. In the statistics literature, various two-stage estimators,
or more precisely, two-step estimators, have been proposed to successfully solve some very
important non-convex problems, for example, sparse principle analysis [10, 54], community
detection [22, 23], mixture of linear regression [11, 58] and crowdsourcing [60, 21]. For all
those problems, under a strong assumption that the initial estimator is consistent, one-step
update in the second stage usually leads us to a minimax optimal estimator. However, as
observed in various simulation or real data studies in [60] and [22], the initial estimator may
perform poorly, and more iterations in the second stage keeps driving down the clustering
error. Unfortunately, due to some technical difficulties, theoretical analyses in [60] and [22]
restrict to one-step iteration.

2

1.1

Our contributions

In this paper, we give a considerably weak initialization condition under which Lloyd’s algorithm converges to the optimal label estimators of sub-Gaussian mixture model. While
previous results [37, 8] focus on exact recovery (strong consistency) of the labels, we obtain
the clustering error rates of Lloyd’s algorithm under various signal-to-noise levels. As a special case, we obtain exact recovery with high probability when the signal-to-noise level is
bigger than 4 log n. The signal-to-noise ratio condition for exact recovery is weaker than the
state-of-the-art result [8]. In contrast to previous two-stage (two-step) estimators, our analyses go beyond one-step update. We are able to show a linear convergence to the statistical
optimal error rate for Lloyd’s algorithms and its two variants for community detection and
crowdsourcing.
We illustrate our contributions here by considering the problem of clustering two-component
spherical Gaussian mixtures, with symmetric centers θ ∗ and −θ ∗ ∈ Rd and variance σ 2 . See
p
Section 2 for more details. Let n be the sample size and r = kθ ∗ k/(σ 1 + 9d/n) be the
normalized signal-to-noise ratio. We establish the following basin of attractions of Lloyd’s
algorithm.
Theorem 1.1. Assume r ≥ C and n ≥ C for a sufficiently large constant C. For symmetric, two-component spherical Gaussian mixtures, given any initial estimator of labels with
clustering error
√
1
1 2.56 + log r
− √ , w.h.p.
A0 < −
2
r
n
Lloyd’s algorithm converges linearly to an exponentially small rate after ⌈3 log n⌉ iterations,
which is the minimax rate as r → ∞ w.h.p.
The results above are extended to general number of clusters k and to (non-spherical)
sub-Gaussian distributions under an appropriate initialization condition and a signal-to-noise
ratio condition, which, to the best of our knowledge, are the weakest conditions in literature.
Now we discuss the contributions of this paper in detail. Our contributions are three folds.
First, we give statistical guarantees of Lloyd’s algorithm. Starting with constant clustering
error, by an appropriate initializer, such as spectral clustering, we show an exponentially
small clustering error rate of Lloyd’s algorithms, under a weak signal-to-noise ratio condition.
We also provide a rate-matching lower bound to show that Lloyd’s algorithm initialized by
spectral clustering is minimax optimal. When the clusters sizes are of the same order and
the distance between different centers are of the same order, our signal-to-noise condition
p
reduces to ∆ & σk 1 + kd/n, where ∆ is the minimum Euclidean distances between two
different cluster centers. Previous results on Lloyd’s algorithms focus on exact recovery of
the labels and they assume sample size n ≫ kd [37, 8]. The best known signal-to-noise ratio
condition is ∆ & σk polylog n [8]. Our condition is weaker by a polylog n factor. Moreover,
our results also hold for the high-dimensional case where d could be larger than n.
Second, we provide computational guarantees of Lloyd’s algorithm. We prove a linear
convergence rate of Lloyd’s iterations given a label initializer of constant clustering error.
3

The linear convergence rate depends on the signal-to-noise ratio. Larger signal-to-noise ratio
leads to faster convergence rate. Counterexamples are constructed in the appendix to show
that the basin of attractions we establish in Theorem 1.1 and in the general k case are almost
necessary. It is worthwhile to point out that the initialization condition in Theorem 1.1 is just
slightly stronger than random guess. As implied by the minimax lower bound (Theorem 3.3),
a necessary condition for consistently estimating the labels is ∆
σ → ∞. In addition, when
1/4
∆ & σd , we are able to prove a random initialization scheme works with high probability.
Third, we develop new proof techniques for analyzing two-stage estimators. Previous
analyses usually require the first stage’s estimator to be consistent to achieve the minimax
optimal rates [60, 22]. In contrast, we only need a constant clustering error of the first stage
estimator by introducing a new technique to analyze random indicator functions. In addition, we are able to improve previous signal-to-noise ratio conditions in community detection
and crowdsourcing by considering two variants of Lloyd’s algorithms. Simulated and real
data experiments show our algorithms are competitive as compared to the state-of-the-art
algorithms for both problems.

1.2

Related work on Gaussian mixture models

The study of Gaussian mixture model [47] has a long and rich history. We give a very brief
review here. In the original paper of Pearson [47], methods of moments estimators were first
proposed, followed by the work of [18, 39]. They all involve solving high-order moments equations, which is computationally challenging as the sample size or the dimensionality grows.
Recently, third-order moments was proved to be sufficient by using tensor decomposition
techniques [12, 5, 28]. However, their sample size requirement is in high-order polynomial
of d and k thus can not be extended to the high-dimensional case. Another line of research
focuses on spectral projections and their variants. By using different dimension reduction
technique, they keep improving the separation condition for general Gaussian mixture models [42, 1, 33, 44, 25] and spherical Gaussian mixture models [53]. These works almost all
focus on estimating the centers. Little is known about the convergence rate of estimating the
cluster labels.
Another popular algorithm for Gaussian mixture model is the EM algorithm [19]. [57]
analyzed the local convergence of EM for well-separated Gaussian mixtures. [15] showed that
a two-round variant of EM algorithm with a special initializer is consistent under a strong
√
separation condition that ∆ & σ d. Recently, a lot of attention has been gained on the
computational guarantees of EM algorithm. Two component, symmetric Gaussian mixture
is most widely studied due to its simple structure. For this model, [9] first proves the linear
convergence of EM if the center is initialized in a small neighborhood of the true parameters,
which implies an exponentially small error rate of labels after one-step label update. [36]
extends their basin of attraction to be the intersection of a half space and a ball near the
origin. More recently, [56, 16] prove the global convergence of the EM algorithm given infinite
samples. While these results are encouraging, it is unclear whether their technique can be
generalized to the k-mixtures or non-Gaussian cases. [30] provides some interesting negative
4

results of EM algorithms. When k ≥ 3, counterexamples are constructed to show there is no
general global convergence of EM algorithm by uniformly initializing the centers from data
points.

1.3

Organization and Notation

The rest of this paper is organized as follows. We introduce the sub-Gaussian mixture problem
and Lloyd’s algorithm in Section 2. Section 3 gives statistical and computational guarantees of Lloyd’s algorithm, and discusses their implications when using spectral initialization.
In Section 4, we consider two variants of Lloyd’s algorithms for community detection and
crowdsourcing and establish their theoretical properties. A number of experiments are conducted in Section 5 that confirm the sharpness of our theoretical findings and demonstrate
the effectiveness of our algorithms on real data. In Section 6, we discuss our results on random initialization and the error of estimating centers. Finally, we present the main proofs in
Section 7, with more technical part of the proofs deferred to the appendix.
Throughout the paper, let [m] = {1, 2, · · · , m} for any positive integer m. For any vector
a, kak is the ℓ2 norm of a. I{·} is the indicator function. Given two positive sequences {an }
and {bn }, an & bn means there is a universal constant C such that an ≥ Cbn for all n, and
define an . bn vice versa. We write an ≍ bn if an & bn and an . bn . Denote an = o(bn ) if
an /bn → 0 as n → ∞.

2

Model and Lloyd’s algorithm

In this section, we first introduce the mixture of sub-Gaussians, a standard model for k-means,
then briefly review Lloyd’s algorithm, followed by introducing spectral clustering algorithms
as a way of initialization.

2.1

Mixture of sub-Gaussians

Suppose we observe independent samples y1 , y2 , · · · , yn ∈ Rd from a mixture of k subGaussian distributions,
yi = θzi + wi for i ∈ [n],
(1)
where z1 , z2 , · · · , zn ∈ [k] are the underlying labels, and θ1 , · · · , θk ∈ Rd are unknown centers
of those k distributions. We assume the noise {wi , i ∈ [n]} are independent zero mean subGaussian vectors with parameter σ > 0, i.e.
Eeha,wi i ≤ e

σ 2 kak2
2

,

for all i ∈ [n] and a ∈ Rd .

(2)

A special case of model (1) is the symmetric, two-component mixture model, in which the
two centers are θ ∗ and −θ ∗ . We observe y1 , y2 , · · · , yn from the following generative model,
yi = zi θ ∗ + ξi ,

5

(3)

where zi ∈ {−1, 1} and {ξi , i ∈ [n]} are independent Gaussian noise with covariance matrix
σ 2 Id . Here with a little abuse of notation, we also use zi ∈ {−1, 1} to denote the underlying
labels. As arguably the simplest mixture model, this special model has recently gained
some attention in studying the convergence of EM algorithm [9, 36, 56]. Other examples of
model (1) are planted partition model [42] for random graphs, stochastic block model [27]
for network data analysis and Dawid-Skene model [17] for crowdsourcing (see Section 4 for
more tails).
For these clustering problems, our main goal is to recover the unknown labels zi rather
than to estimate the centers {θj }. Note that the cluster structure is invariant to the permutations of label symbols. We define the mis-clustering rate of estimated labels ẑ1 , · · · , ẑn
as
" n
#
1X
I {π(ẑi ) 6= zi } ,
(4)
L(ẑ, z) = inf
π∈Sk n
i=1

where Sk is the collection of all the mappings from [k] to [k].

2.2

Lloyd’s Algorithm

Lloyd’s algorithm was originally proposed to solve the following k-means problem. Given n
vectors y1 , y2 , · · · , yn ∈ Rd and an integer k, the goal is to find k points θ1 , θ2 , · · · , θk ∈ Rd
to minimize the following objective function
X
(5)
min kyi − θj k2 .
i∈[n]

j∈[k]

The problem above is equivalent to find θ̂1 , · · · , θ̂k ∈ Rd and ẑ1 , · · · , ẑn ∈ [k] such that
(θ̂, ẑ) = argmin
(θ,z)

X

i∈[n]

yi −

k
X

2

θj I{zi = j}

.

(6)

j=1

From a statistical point of view, (6) is the maximum likelihood estimator of model (1) with
spherical Gaussian noise. It has been proved by Pollard [50, 48] the strong consistency and
central limit theorem of using (5) to estimate the centers θ.
Lloyd’s algorithm is simply motivated by the following observation of (6). If we fix θ, ẑi
is the index of the center that yi is closest to, and if we fix z, θ̂j is the sample mean of those
yi with zi = j. If we start with an initial estimate of centers or labels, we can iteratively
update the labels and centers. Therefore, we have the following Lloyd’s algorithm.
1. Get an initial estimate of the centers or the labels.
2. Repeat the following iteration until convergence.
2a. For h = 1, 2, · · · k,

Pn
(s)
yi I{ẑi = h}
(s)
.
θ̂h = Pi=1
(s)
n
I{ẑ
=
h}
i=1
i
6

(7)

2b. For i = 1, 2, · · · n,

(s+1)

ẑi

(s)

= argmin kyi − θ̂h k2 .

(8)

h∈[k]

For model (3), since the centers are parametrized by a single vector θ ∗ , the center update
step (7) can be simplified to
n

θ̂

(s)

=

1 X (s)
ẑi yi .
n

(9)

i=1

Lloyd’s algorithm is guaranteed to converge because it is a greedy algorithm that keeps
decreasing the value of objective function (6). However, it may only converge to a local
minimum and thus a good initializer is needed. One possible way of initializing Lloyd’s
algorithm is to use spectral methods. Note that the signal matrix θ in model (1) is low rank.
Spectral methods first do a de-noising step by projecting the data onto the subspace spanned
by its top singular vectors, which approximately preserves the clustering structure. Since
solving the k-means problem is NP-hard in general [14], we run an approximated k-means
algorithm [51] on the projected data matrix, whose running time is polynomial in n, k and
d. There are many versions of spectral clustering algorithms and here is the pseudo code of
the one used in this paper.
1. Compute the SVD of the data matrix Y = [y1 , · · · , yn ] = U DV ′ . Let Uk be the first k
columns of U .
2. Project y1 , · · · , yn onto Uk , i.e. let yˆi = Uk Uk′ yi for i ∈ [n].
3. Run a O(1)-approximation algorithm [51] for k-means problem on the columns of projected matrix Ŷ = [ŷ1 , · · · , ŷn ].
We refer this algorithm as spectral clustering algorithm henceforth. A lot of progress has
been made in studying theoretical properties of this spectral clustering algorithm. For more
details, we refer to [34, 8] and references therein.

3

Clustering sub-Gaussian mixtures

In this section, we present the results of using Lloyd’s algorithm to clustering sub-Gaussian
mixtures. In Section 3.1, we give the convergence results of Lloyd’s algorithm to the symmetric, two component mixture model (3). Then we extend it to the general k mixture model
(1) in Section 3.2. A minimax lower bound is established in Section 3.3.

3.1

Two mixtures

A key quantity in determining the basin of attraction of Lloyd’s algorithm is the following
normalized signal-to-noise ratio
kθ ∗ k
,
r= √
σ 1+η
7

√
where η = 9d/n. Here we normalize the signal-to-noise ratio by 1 + η because the statistical
√
√
precision of estimating θ ∗ is at the order of σ η. If kθ ∗ k . σ η, information theoretically we
could not distinguish between two centers with positive probability, even when the labels are
P
(s)
known. Let As = n1 ni=1 I{ẑi 6= zi } be the mis-clustering rate at the iteration s of Lloyd’s
algorithm, and θ̂ (s) be the estimated centers at step s, s = 0, 1, 2, · · · . We have the following
theorem that characterizes the behavior of As .
Theorem 3.1. Assume n ≥ C and r ≥ C for a sufficiently large constant C. For any given
(data dependent) initializer satisfying
√


1
4
1 2.56 + log r
(0)
∗
−√
kθ ∗ k,
(10)
or kθ̂ − θ k ≤ 1 −
A0 ≤ −
2
r
r
n
with probability 1 − ν, we have


8
As + 2
r



2
As + 2 +
r

r

(11)

for all s ≥ 3 log n

(12)

4 log n
, for all s ≥ 0
n


∗ 2
with probability greater than 1 − ν − n−3 − 2 exp − kθ3σk2 , and
As+1 ≤

kθ ∗ k2
As ≤ exp −
16σ 2




,



∗ k2
. Moreover, if r → ∞ as n → ∞,
with probability greater than 1 − ν − 5n−1 − 8 exp − kθ
16σ2


∗ 2
the error rate in (12) can be improved to exp −(1 + o(1)) kθ2σk2 .

The proof of Theorem 3.1 is given in Section 7. We will first show that after a few
iterations, the mis-clustering ratio As is sufficiently small (for example, smaller than 1/8) by
proving (11), then prove mis-clustering rate of Lloyd’s algorithm is exponentially small with
an exponent determined by the signal-to-noise ratio. Note that the mis-clustering rate As
√
only takes discrete values in {0, n−1 , 2n−1 , · · · , 1}. If kθ ∗ k > 4σ log n, (12) implies As < 1/n,
and thus Theorem 3.1 guarantees a zero clustering error after ⌈3 log n⌉ Lloyd’s iterations with
high probability.
By Theorem 3.1, the convergence of Lloyd’s algorithm has three stages. In the first
iteration, it escapes from a small neighborhood around 1/2 and achieves an error rate As ≤
1/2 − c for a small constant c. Then it has a linear convergence rate, which depends on
the signal-to-noise ratio. Finally, similar to other two-stage estimators [60, 22], once the the
mis-clustering rate As is sufficiently small, it converges to the optimal statistical precision of
the problem after one or more iterations.
As we shall see in Section 3.3, a necessary condition for consistently estimating the labels
is kθ ∗ k/σ → ∞. Therefore, our condition on the initializer is just slightly stronger than
the random initialization error 1/2. Balakrishnan et al. [9] studied the convergence of EM
algorithm under the same model. They require an initializer θ̂ (0) ∈ B(θ ∗ , kθ ∗ k/4) under the
assumption that r is sufficiently large. Here we replace 1/4 by a factor close to 1 in Theorem
8

3.1. More specifically, the factor in the initialization condition is determined by the signalto-noise ratio. The stronger the signal-to-noise ratio is, the weaker initialization condition
we need.

3.2

k mixtures

Now we consider general number of mixtures. To better present our results, we first introduce
(s)
some notation. For all h ∈ [k], let Th∗ be the true cluster h and Th be the estimated cluster
(s)
(s)
(s)
(s)
h at iteration s. Define n∗h = |Th∗ |, nh = |Th | and nhg = |Th∗ ∩ Tg |. The mis-clustering
rate at iteration s can be written as
n

As =

1X
1
(s)
I{ẑi 6= zi } =
n
n
i=1

X

(s)

ngh .

g6=h∈[k]2

We define a cluster-wise mis-clustering rate at iteration s as

P
(s) 
(s) P

n
n
g6=h∈[k] hg
g6=h∈[k] gh
.
,
Gs = max
(s)

n∗h
h∈[k] 
nh

The first term in the maximum operator of definition of Gs can be understood as the false
positive rate of cluster h and the second term is the true negative rate of cluster h. It is easy
to see the relationship that As ≤ Gs .
(s)
Let ∆ = ming6=h∈[k] kθg − θhk be the signal strength. For h ∈ [k], let θ̂h be the estimated
center of cluster h at iteration s. Define our error rate of estimating centers at iteration s as
Λs = max
h∈[k]

1 (s)
kθ̂ − θh k.
∆ h

Besides signal-to-noise ratio, there are other two factors that determine the convergence of
Lloyd’s algorithm, the maximum signal strength and the minimum cluster size. As argued
by [30], when one cluster is very far away from other clusters, local search may fail to find
the global optimum, which indicates initialization should depend on the maximum signal
strength. Since the cluster sizes influence the accuracy to estimate centers, we a lower bound
on the size of the smallest clusters. Define λ = maxg6=h∈[k] kθg − θh k/∆ and α = minh∈k n∗h /n.
Similar to the two-cluster case, we define a normalized signal-to-noise ratio
r
α
∆
.
(13)
rk =
σ 1 + kd/n
Now we are ready to present our results for the k mixtures.
√
Theorem 3.2. Assume nα2 ≥ Ck log n and rk ≥ C k for a sufficiently large constant C.
Given any (data dependent) initializer satisfying


1
6
1
4
1
−√
or Λ0 ≤ − √ ,
(14)
G0 <
2
rk λ
2
rk
9

with probability 1 − ν, we have
Gs+1

C1
C1
≤ 2 Gs + 2 +
rk
rk

r

5k log n
α2 n

for all s ≥ 1

for some constant C1 ≤ C with probability greater than 1 − ν − n−3 , and


∆2
As ≤ exp −
for all s ≥ 4 log n
16σ 2

(15)

(16)

√
with probability greater than 1 − ν − 4/n − 2 exp (−∆/σ). Moreover, if k = o(rk ) and

∆2
k log n = o(nα2 ) as n → ∞, the exponent in (16) can be improved to exp −(1 + o(1)) 8σ
2 .

In Theorem 3.2, we establish similar results as in the two-mixture case. The initialization
condition (14) here is slightly stronger due to the asymmetry. Similarly, after initialization
that satisfies (14), the convergence has three stages. It first escapes from a small neighborhood
of 1/2, then converges linearly, and finally it achieves an exponentially small mis-clustering
rate after ⌈4 log n⌉ iterations.
From Theorem 3.2, Lloyd’s algorithm does not require to know or to estimate the covariance structure of Gaussian mixtures. Likelihood-based algorithms, such as EM algorithm
and methods of moments, need to estimate the covariance structure. Lloyd’s algorithm works
for any general sub-Gaussian mixtures satisfying (2), and thus is robust to different models.
This demonstrates a key advantage of using Lloyd’s algorithm than other likelihood-based
algorithms.
Theoretically, the initialization condition (14) in Theorem 3.2 is sharp in a sense that we
give a counterexample in Section A.5 showing that Lloyd’s algorithm may not converge when
G0 = 1/(2λ) or Λ = 1/2.
Now let us give a sufficient condition for the spectral clustering initializer to fall into the
basin of attraction (14). By Claim 1 in Section 3.2 of [34] and Lemma 7.1, for each center θh ,
√
there is a center µh returned by spectral clustering algorithm such that kµh −θh k . ( k/rk )∆
with probability greater than 1 − exp(−n/4), which implies the following corollary.
√
Corollary 3.1. Assume k = o(rk ) and k log n = o(nα2 ) as n → ∞. Let ẑ be the output of
Lloyd’s algorithm initialized by spectral clustering algorithm after ⌈4 log n⌉ iterations. Then


∆2
(17)
ℓ(ẑ, z) ≤ exp −(1 + o(1)) 2
8σ
with probability greater than 1 − 5n−1 − 2 exp (−∆/σ).
This corollary gives a sufficient separation (signal-to-noise ratio) condition for clustering
sub-Gaussian mixtures. When n ≥ kd and all the cluster sizes are of the same order, our
separation condition simplifies to ∆ & kσ. When there are finite number of clusters (k is
finite), our separation condition further simplifies to ∆ & σ. To our knowledge, this is the
first result to give an explicit exponentially small error rate for estimating the labels. Previous
10

studies mostly focus on exact recovery of the cluster labels and have no explicit statistical
convergence rate. Furthermore, our results hold for all the range of d, while previous results
all require n ≫ kd.
Previously, the best known separation condition on the Lloyd-type algorithm is rk &
√
k poly log n, under the assumption that n ≫ kd/α [8]. Corollary 3.1 improves it by a
poly log n factor. For Gaussian mixtures, among all algorithms including Lloyd’s algorithm,
p
the best known separation condition is ∆ & σ(α−1/2 + k2 + k log(nk)) under the assumption
that αn & k(d + log k) [1]. Our condition is weaker when the cluster sizes are of the same
√
order. When k is finite, our condition is weaker by a log n factor. As we shall see in Section
5, k is usually small in practice and can often be regarded as a constant.

3.3

Minimax lower bound

To show that the mis-clustering rate rate in Corollary 3.1 cannot be improved, we present a
rate matching lower bound in this section. Define a parameter space as follows,

Θ = (θ, z), θ = [θ1 , · · · , θk ] ∈ Rd×k , ∆ ≤ min kθg − θh k,
g6=h
(18)

z : [n] → [k], |{i ∈ [n], zi = u}| ≥ αn, ∀u ∈ [k] .
We have the following minimax lower bound.
i.i.d

Theorem 3.3. (Lower Bound) For model (1), assume independent gaussian noise wij ∼
∆
N (0, σ 2 ), then when σ log(k/α)
→ ∞,
∆2
inf sup Eℓ(ẑ, z) ≥ exp −(1 + o(1)) 2
ẑ (z,θ)∈Θ
8σ


If

∆
σ



.

+ log(k/α) = O(1), inf ẑ sup(z,θ)∈Θ Eℓ(ẑ, z) ≥ c for some constant c > 0.

This lower bound result shows that if the signal-to-noise ratio ∆/σ is at the constant
order, the worst case mis-clustering rate is lower bounded by a constant. In other words, a
necessary condition for consistently estimating the labels is ∆/σ → ∞. To achieve strong
√
consistency, we need ∆/σ & log n. Theorem 3.3 indicates the necessity of separation
conditions in estimating the labels. Previous results that use methods of moments to estimate
the cluster centers [28, 5], on the contrary, do not require separation conditions. This reveals
the difference between estimating the cluster centers and the labels.
Together with Corollary 3.1, Theorem 3.3 gives us the minimax rate of estimating the
underlying labels of Gaussian mixtures.

4

Applications

In this section, we generalize the results in Section 3 to community detection and crowdsourcing by considering two variants of Lloyd’s algorithm.
11

4.1

Community Detection

We have a network of k communities. There are n nodes in the network and {zi ∈ [k], i ∈ [n]}
is the community assignment. We observe a symmetric adjacency matrix A ∈ {0, 1}n×n
with zero diagonals, which is generated by the following stochastic block model [27] (SBM
thereafter),
(
Ber( na ), if zi = zj
Aij ∼
Ber( nb ), otherwise,
for 1 ≤ i < j ≤ n, where 0 < b < a < n and {zi } are unknown parameters. The SBM is a
special case of the Planted Partition Model [42] in the theoretical computer science literature.
As in [22], we use a single parameter β to control the community sizes, namely, there is a
n
β > 0 such that βk
≤ |{i ∈ [n], zi = g}| for all g ∈ [k]. The goal of community detection is to
estimate community structure z using the adjacency matrix A. Note that EAi· only takes k
different values, and the SBM can be viewed as a special case of model (1) with d = n and
σ 2 = na , ignoring the fact that A is symmetric with zero diagonals.
When a and b are of an order of n, the SBM satisfies the sub-Gaussian noise condition (2)
2
2(a−b)2
2n
×
and we have ∆2 ≥ (a−b)
=
and σ 2 ≤ a/n. Initialized by the spectral clustering
2
βk
βkn
n
algorithm, let ẑ be the labels returned by running Lloyd’s algorithm on the rows of adjacency
matrix A with ⌈4 log n⌉ iterations. By considering the fact that A is symmetric with zero
diagonals, we can slightly modify the proof of Corollary 3.1 to obtain the following result.
Corollary 4.1. Assume

(a−b)2
β 3 ak 4

→ ∞ and β 2 k3nlog n → ∞ as n → ∞, then


(a − b)2
ℓ(z, ẑ) ≤ exp −(1 + o(1))
4aβk

(19)

with probability tending to 1.
We omit the proof of Corollary 4.1. A more interesting case of the SBM is the sparse case
when a = o(n) and b = o(n), in which the sub-Gaussian assumption (2) does not hold. By
utilizing the community structure, we consider the following variants of Lloyd’s algorithm.
As pointed out by [13], under the sparse setting of a . log n, the trimmed adjacency
matrix Aτ is a better estimator of EA than A under the spectral norm. Therefore, we run
spectral clustering on the Aτ to get an initial estimator of z. Then we update the labels
node by node. For the i-th node in the network, we estimate its connecting probability to
community j (defined as Bzi j ) based on our previous estimated labels. Then we assign its
label to be the j that maximizes B̂zi j , which is expected to be close to a/n or b/n, depending
on whether zi = j. The following theorem gives a theoretical guarantee of the CommuLloyd
algorithm.
Theorem 4.1. Assume n ≥ 6k2 log n, (a − b)2 ≥ C0 aβ 2 k2 log(βk) for a sufficiently large
constant C0 . Let As be the mis-clustering rate at iteration s of the CommuLloyd algorithm.
Then for any initializer satisfies G0 ≤ 1/4 with probability 1 − ν, we have


(a − b)2
4
(22)
As+1 ≤ exp −
+ As , ∀ 1 ≤ s ≤ 3 log n
2Cβak
5
12

Algorithm 1: CommuLloyd
Input: Adjacency matrix A. Number of communities k.. Trimming threshold τ .
Output: Estimated labels ẑ1 , · · · , ẑn .
1. Trim the adjacency matrix:
P
1a. calculate the degree of each node di = nj=1 Aij , ∀i ∈ [n].
1b. trim the adjacency matrix Aτij = Aij I{di ≤ τ }, ∀i ∈ [n], j ∈ [n].

2. Run spectral clustering algorithm on the trimmed matrix Aτ .
3. Run following iterations until converge,
Pn
j=1 Aij I{ẑj = h}
P
, ∀ h ∈ [k], i ∈ [n].
B̂ih =
n
j=1 I{ẑj = h}
ẑi = argmax B̂ih ,
h∈[k]

∀ i ∈ [n].

(20)
(21)

with probability greater than 1 − n−1 − ν.
The proof of Theorem 4.1 is given in Section B of the appendix. We show the linear convergence of CommuLloyd algorithm given the first stage group-wise mis-clustering rate G0 ≤ 1/4.
p
In fact, this initialization assumption can be relaxed to G0 ≤ 1/2 − a log(βk)βk/(a − b). To
better present our results, we simplify it to 1/4. By Lemma 7 in [22] and Theorem 3.1 in [8],
we can guarantee a group-wise initialization error of 1/4 when (a − b)2 & aβ 2 k3 . Therefore,
we obtain an exponentially small error rate under the signal-to-noise ratio condition
(a − b)2 & aβ 2 k3 log β.
Theorem 4.1 provides theoretical justifications of the phenomenon observed in the numerical experiments of [22]. More iterations achieve better mis-clustering rate. The Lloyd
iterations enjoy similar theoretical performance as the likelihood based algorithm proposed in
[22]. While they require a global initialization error to be o(1/(k log k)), we require a clusterwise initialization error to be smaller than 1/4. Moreover, the CommuLloyd algorithm is
computationally much more efficient than Algorithm 1 in [22], which requires obtaining n
different initializers and hence running SVD on the network n times. Theoretically, we relax
the assumption in [22] that a ≍ b and β is a constant, and we improve the best known signalto-noise ratio condition[22] by a log k factor. When β is a constant and (a − b)2 /(ak3 ) → ∞
as n → ∞ , we are able to match the minimax rate obtained in [59].

4.2

Crowdsourcing

Crowdsourcing is an efficient and inexpensive way to collect a large amount of labels for
supervised machine learning problems. We refer to [60, 21] and references therein for more
13

details. The following Dawid-Skene model is the most popular model considered in the
crowdsourcing literature.
Suppose there are m workers to label n items. Each item belongs to one of the k categories.
Denote by z = (z1 , z2 , · · · , zn ) ∈ [k]n the true labels of n items and by Xij the label of
worker i given to item j. Our goal is to estimate the true labels z using {Xij }. DawidSkene model assumes that the workers are independent and that given zj = g, the labels
given by worker i are independently drawn from a multinomial distribution with parameter
πig∗ = (πig1 , · · · , πigk ), i.e.,
P {Xij = h|zj = g} = πigh
for all g, h ∈ [k], i ∈ [m], j ∈ [n]. Dawid-Skene model can be seen as a special case of the
mixture model (1) via the following transformation. For j ∈ [n], let

′
yj = I{X1j = 1}, · · · , I{X1j = k}, · · · , I{Xmj = 1}, · · · , I{Xmj = k} .

Then given zj = g, we have E[yj |zj = g] = θg with θg = (π1g1 , π1g2 , · · · , πmgk ). By defining
wijh = 1{Xij = h} − πizj h , our observations yj can be decomposed as the following signalplus-noise model (1) with
wj = [w1j1 , w1j2 , · · · , wmjk ]′ .
Therefore, we consider the following variants of the Lloyd’s algorithm in Algorithm 2. We
iteratively estimate workers’ accuracy and items’ labels. Iterations (24) and (25) are actually
Lloyd’s iteration (7) and (8) under the transformed Dawid-Skene model. Iterations (24) and
(25) are also closely related to the EM update for Dawid-Skene model [60]. (24) is the same
as the M-step. In (25), we first use least squares to approximate the log-likelihood function
of multinomial distributions, then we do a hard labeling step instead of soft labeling in the
E-step of EM algorithm.
As suggested by Corollary 3.1, we could use spectral clustering to initialize the labels.
But the majority voting initializer, as described in the first step of Algorithm 2, is a more
natural and commonly used initializer of crowdsourcing. By regarding all the workers have
the same accuracy, it simply estimates the labels by aggregating the results of every worker
with equal weights. The majority voting initializer is computationally more efficient than
spectral methods. We have the following result of the majority voting estimator.
Theorem 4.2. Under the Dawid-Skene model, the majority voting initializer satisfies


n
1X
mV (π)
I {ẑj 6= zj } ≤ exp −
n
2
j=1

with probability greater than 1 − exp (−mV (π)/2 + log k), where
v
v
2
um
um
X
X
u
u
1
πigg − t
πigh  .
V (π) = min t
g6=h m
i=1

14

i=1

(26)

Algorithm 2: CrowdLloyd
Input: {Xij , i ∈ [m], j ∈ [n]}. Number of possible labels k.
Output: Estimated labels ẑ1 , · · · , ẑn .
1. Initialization via Majority Voting:
ẑj = argmax
g∈[k]

m
X

I{Xij = g}

i=1

∀ j ∈ [n].

2. Run following iterations until converge.
Pn
j=1 I{Xij = h, ẑj = g}
∀ i ∈ [m], g, h ∈ [k]2 .
π̂igh = Pk Pn
I{X
=
h,
ẑ
=
g}
ij
j
h=1
j=1


m X
k
X
ẑj = argmin 
(I{Xij = g} − π̂igh )2  ∀ j ∈ [n].
h∈[k]

(23)

(24)

(25)

i=1 g=1

The proof of Theorem 4.2 is given in Section C.2. To our knowledge, this theorem is the
first theoretical characterization of the Majority Voting estimator under the general DawidSkene model. mV (π) measures the collective accuracy of m workers on estimating one of the
labels. If there are two labels g and h that are confusing for most of the workers, the value
of V (π) is small and majority voting estimator may not have good performance. Under the
1−p
for g 6= h, we have
special case of one-coin model, where πigg = p and πigh = k−1
p
√
mV (π) = m( p − (1 − p)/(k − 1))2 & log m
q
as long as p ≥ k1 + logmm . That is, for Majority Voting estimator to be consistent, the
p
accuracy of workers giving true labels only needs to be log m/m better than random guess.
With Theorem 4.2 and Theorem 3.2, we are able to to obtain an error rate upper bound
of Algorithm 2. It can be proved that the sub-Gaussian parameter of wi is 2 (see the proof
P Pk
2
of Corollary 4.2 for more details). Define ∆2 = ming6=h m
i=1
u=1 (πigu − πihu ) . Then we
have the following corollary of Theorem 3.2.
Corollary 4.2. Assume mV (π) ≥ C log k, α∆2 ≥ Ck and nα ≥ k∆2 log n for a sufficiently
large constant C. Let ẑ = (ẑ1 , · · · , ẑn ) be the estimated labels returned by running Algorithm
2. We have


n
1X
∆2
I {ẑi 6= zi } ≤ exp −
(27)
n
64
i=1

with probability greater than 1 − exp(−mV (π)/4) − 4/n − 2 exp (−∆/2).

The proof of Theorem 4.2 is given in Section C.2. We achieve an exponentially small
error rate for estimating the labels under Dawid-Skene model. By Theorem 4.2 in [21], we
15

can achieve the minimax optimal rate by an additional likelihood based Neyman-Pearson
testing step followed by Algorithm 2. Previous provable results [60, 21] assume πigh ≥ ρ > 0
for all (i, g, h) ∈ [m] × [k]2 , which is quite restrictive since a very good worker could have
πigh = 0 for some h 6= g. By considering least square update (25), we successfully get rid
of this assumption. In contrast to algorithms that using spectral methods as initializer [60],
Theorem 4.2 does not need an eigenvalue condition on the confusion matrix.

5

Numerical Analysis

5.1

Gaussian Mixture Model

5.1.1

Simulation

−1
−2
1

2

−4

k=50
k=25
k=10
k=5

−5

−5
0

−3

log(As)

−3

SNR=6
SNR=7
SNR=8
SNR=9

−4

log(As)

−2

−1

In this section, we provide simulation results that are designed to test our theoretical findings
in Theorem 3.2. As predicted by Theorem 3.2, there are three stages of the convergence of
Lloyd’s algorithm. Given an initial estimator of the labels with error smaller than 1/(2λ),
Lloyd’s algorithm escapes from a small neighborhood of 1/(2λ) in the first iteration. Then it
has a geometrically decay mis-clustering rate during the second stage. Finally, it converges to
the statistical precision of the problem and achieves an exponentially small error. In all the
results reported below, we plot the logarithm of the mis-clustering rate versus the iteration
count. Each curve plotted is an average of 10 independent trails.

3

4

5

6

0

iteration count

1

2

3

4

5

6

iteration count

(a)

(b)

Figure 1: Plots of the log mis-clustering rate As versus iteration count.
In the first simulation, we independently generate n = 1000 samples from a mixture of k =
10 spherical Gaussians. Each cluster has 100 data points. The centers of those 10 gaussians
are orthogonal unit vectors in Rd with d = 100. The four curves in Figure 1a corresponds to
16

four choices of the standard deviation σ = 2/SN R with SN R ∈ {6, 7, 8, 9}. We use the same
initializer with a group-wise mis-clustering rate of 0.454 for all the experiments. As shown
by plots in Figure 1a, there is a linear convergence of the log mis-clustering rate at the first 2
iterations and then the convergence slows down. After four iterations, the log mis-clustering
1
SN R2 . Moreover, we observe in Figure 1a that the
rate plateaus roughly at the level of − 16
convergence rate increases as SN R increases, as expected from Theorem 3.2.
The second experiment investigates the role of k in the convergence of Lloyd’s algorithm.
We first randomly generate k ∈ {5, 10, 25, 50} centers that are orthogonal unit vectors in R100 .
Then we add independent gaussian noise N (0, 0.252 ) to each coordinate of those centers to
get 1000/k samples for each cluster. In total, we have n = 1000 independent samples from
a mixture of spherical Gaussian distributions. We set the same initializer as in the first
simulation to make a group-wise mis-clustering rate of 0.37. As observed in Figure 1b,
the larger the k, the slower Lloyd’s algorithm converges. This agrees with the phenomena
predicted by Theorem 3.2. It is interesting to observe that when k is large, Lloyd’s algorithm
does not converge to the optimal mis-clustering rate in Figure 1b. When the value of k is
doubled, the convergence rate of log As slows down roughly by half, which matches the k-term
in (15).
5.1.2

Gene Microarray Data

To demonstrate the effectiveness of Lloyd’s algorithm on high dimensional clustering problem,
we use six gene microarray data sets that were preprocessed and analyzed in [20, 24]. In each
data set, we are given measurements of the expression levels of d genes for n people. There
are several underlying classes (normal or diseased) of those people. Our goal is to cluster
them based on their gene expression levels, with the number of groups k given. One common
feature of those six data sets is that the dimensionality d is much larger than the number of
samples n. As shown in Table 1, while n is usually smaller than two hundreds, the number
of genes can be as large as ten thousands.
We compare three different methods on these six data sets. RandLloyd is Lloyd’s algorithms initialized by randomly drawing k rows of the original data as k centers. SpecLloyd
runs RandLloyd on the first k eigenvectors of the data, followed by Lloyd’s iterations on the
original data. IF-HCT-PCA, proposed in [32], does a variable selection step first and then applies spectral clustering on selected variables. All three methods need random initializations
at the k-means step and thus have algorithmic randomness. We randomly initialize centers
for 30 times and use the initializer that give us the minimum k-means objective value.
In Table 1, we report the average mis-clustering rate rates over 30 independent experiments. We also report the standard deviation if it is larger than .01. As we can see,
RandLloyd is comparable to other two methods on these gene microarray data. It achieves
the best clustering error rate when the clustering task is relatively easy (signal-to-noise ratio
is strong), for example, Leukemia and Lymphoma data sets. Note that both IF-HCT-PCA
and SpecLloyd have a dimension reduction step. This agrees with the phenomena suggested
by Theorem 3.2 that Lloyd’s algorithm also works on the high dimensional setting when the
17

Date set
Brain
Leukemia
Lung(1)
Lymphoma
Prostate
SRBCT

k
5
2
2
3
2
4

n
42
72
181
62
102
63

d
5597
3571
12533
4026
6033
2308

IF-HCT-PCA
.262
.069
.033
.069
.382
.444

SpecLloyd
.119
.028
.083
.177
.422
.492

RandLloyd
.262(.04)
.028
.098
.016
.422
.492

Table 1: Clustering error rate of three clustering methods on six gene microarray datasets.
signal-to-noise ratio is strong enough.

5.2

Community Detection

In this section, we report the results of empirical studies to compare the CommuLloyd algorithm proposed in Section 4.1 with other methods. In particular, we compare to spectral
clustering (Spec) and spectral clustering with one step refinement studied in [22]. Since
the algorithm in [22] is motivated by the maximum likelihood estimator, we refer it to as
Spec+MLE.
5.2.1

Simulation

Using the same settings as in [22], we generate data from SBM under three different scenarios: (1) dense network with equal size communities; (2) sparse network with equal size
communities; and (3) dense network with unequal size communities. The detailed set-up of
three simulations is described in Table 2. We conduct each experiment with 10 independent
repetitions and report the average results of each algorithm in Figure 2.
For all the experiments, we use spectral clustering methods in [38] as our initializer and
plot the logarithm mis-clustering rate versus the Lloyd iteration or MLE iteration counts.
Iteration 0 represents the spectral clustering initializer. As shown by Figure 2, CommuLloyd significantly outperforms spectral clustering and spectral clustering with one step MLE
refinement under all three settings. Surprisingly, it is even better than multiple steps MLE
refinements. In panel (a), MLE refinements converge to a sub-optimal mis-clustering rate,
and in panel (b), it does not improve the spectral clustering initializer. This phenomenon
may be due to the stronger initialization condition for MLE refinements in [22], while in
Theorem 4.1, CommuLloyd needs only a constant-error initializer to converge to the optimal
mis-clustering rate rate. The convergence rate of CommuLloyd for the these two cases, however, is faster than the linear convergence predicted by Theorem 4.1. As for the unbalance
case in panel (c), CommuLloyd and MLE has similar performance and they both exhibit a
linear convergence rate.

18

Setting
Balanced Case
Sparse Case
Unbalanced Case

n
2000
2000
1000

k
10
4
4

Within
0.20
0.019
0.35

Between
0.11
0.005
0.22

Community Size
equal
equal
100,200,300,400

5

10

15

iteration count

(a)

−2.2
−2.6
−3.4

−2.6
0

CommuLloyd
Spec+MLE

−3.0

log(As)

−2.2

CommuLloyd
Spec+MLE

−2.4

−2.0

log(As)

−1.5

CommuLloyd
Spec+MLE

−2.5

log(As)

−1.0

Table 2: The settings of three simulations on community detection. Within is the withincommunity connecting probability a/n. Between is the between-community connecting probability b/n. Equal means the communities are of equal size.

0

5

10

iteration count

(b)

15

0

5

10

15

iteration count

(c)

Figure 2: Plots of log mis-clustering rate As versus the iteration.
5.2.2

Political Blog Data

We now examine the performance of Algorithm 2 on a political blog dataset [2]. Each node
in this data represents a web blog collected right after the 2004 US presidential election. Two
nodes are connected if one contains a hyperlink to the other, and we ignore the direction of
the hyperlink. After pre-processing [35], there are 1222 nodes and 16714 edges. Depending
on whether liberal or conservative, the network naturally has two communities. We use the
label given in [2] as our ground truth.
Table 3 reports the number of mis-labeled nodes of our algorithm on political blog
dataset. Started from spectral clustering with an error rate of 437/1222, CommuLloyd
achieves 56/1222 error rate after three iterations. After that, it fluctuates between 55 and 56
mis-labeled nodes. Our algorithm slightly outperforms the state-of-the-art algorithm SCORE
[31], which reports an error rate of 58/1222.

5.3

Crowdsourcing

In this section, we compare CrowdLloyd algorithm with three other crowdsourcing methods:
(1) Opt-D&S, which is Spectral methods with one step EM update proposed in [60], (2) MV,
the majority voting estimator and (3) MV+EM, which is majority voting estimator followed
19

Methods
No.mis-clustered

Spec
437

Spec+MLE(1st)
132

Spec+MLE(10th)
61

SCORE
58

CommuLloyd
56

Table 3: Number of mis-labeled blogs for different community detection methods. Here Spec
is the spectral clustering method described in [38]. We refer the one-step MLE refinement of
Spec as Spec+MLE(1st), and 10 steps MLE refinements as Spec+MLE(10th). SCORE is a
community detection algorithm for Degree-corrected SBM proposed in [31].
by EM updates. Among those four algorithms, Opt-D&S, MV and CrowdLloyd have provable
error rates.
When there are ties, the majority voting estimator is not clearly defined. To avoid
algorithmic randomness, we return the first label index when there are ties for CrowdLloyd
and MV estimator. Since EM algorithm takes soft labels as input, we use weights rather than
hard label of majority voting estimator in MV+EM. In practice, there are usually missing
data. We denote by 0 as the category of missing data and modify the label update step in
CrowdLloyd as


k 
2
X
X
(s)
(s+1)
(28)
I{Xij = g} − π̂igh .
ẑj
= argmin 
h∈[k]

5.3.1

i,Xij 6=0 g=1

Simulation

We follow the same setting as in the experiments of [60], where there are m = 100 workers
and n = 1000 items. Each item belongs to one of the k = 2 categories. For worker i, the
diagonal entries of his or her confusion matrix πigg are independently and uniformly drawn
from the interval [0.3, 0.9]. We set πigh = 1−πigg for all i ∈ [n] and (h, g) ∈ {(1, 2), (2, 1)}. To
simulate missing entries, we assume each entry Xij is observed independently with probability
p ∈ {0.2, 0.5, 1}.
We replicate experiments independently for 10 times and report the average results in
Table 4. Among three provable algorithms, CrowdLloyd is consistently better than Opt D&S
(1st iteration) and MV under different missing levels. CrowdLloyd, MV+EM and Opt D&S
(10th iteration) have almost identical performance. As expected, CrowdLloyd is slightly worse
because it does not use likelihood function to update the confusion matrices. As discussed
before, an additional Neyman-Pearson testing step followed by CrowdLloyd algorithm would
give the same output as EM updates.
5.3.2

Real Data

We now compare CrowdLloyd algorithm to other three algorithms on five real crowdsourcing
data sets: (1) bird dataset is a binary task of labeling two bird species, (2) recognizing textual
entailment (RTE) dataset contains a set of sentence pairs and the task is to label whether
the second sentence can be inferred from the first one, (3) TREC dataset is a binary task
20

Methods
Opt D&S (1st iteration)
Opt D&S (10th iteration)
MV
MV + EM
CrowdLloyd

p = 1.0
0.06
0.07
2.14
0.07
0.07

p = 0.5
1.32
1.10
8.56
1.10
1.14

p = 0.2
10.32
8.22
18.97
8.17
8.19

Table 4: mis-clustering rate rate (%) of different crowdsourcing algorithms on simulated data.
#
1
2
3
4
5

Date set
Bird
RTE
TREC
Dog
Web

Opt D&S
10.09
7.12
29.80
16.89
15.86

MV
24.07
8.13
34.86
19.58
26.93

MV+EM
11.11
7.34
30.02
16.66
15.74

CrowdLloyd
12.04
6.88
29.19
15.99
14.25

Table 5: mis-clustering rate rate (%) of different crowdsourcing algorithms on five real data
sets
of assessing the quality of information retrieval in TREC 2011 crowdsourcing track, (4) dog
dataset is a subset of Stanford dog dataset that requires labeling one of the four dog breeds,
and (5) web dataset is a set of query-URL pairs for workers to label on a relevance scale from
1 to 5.
The mis-clustering rates of four methods are summarized in Table 5. CrowdLloyd slightly
outperforms the other three on four of those five data sets. One possible explanation is
that the least squares update (25) works for general sub-Gaussian distributions and thus is
potentially more robust than the likelihood based EM update.

6
6.1

Discussion
Random Initialization

In this section, we explore the condition under which random initialization converges to the
global optimum. Here we consider the symmetric two-mixture model (3). Note that Theorem
3.1 holds for any initializer (could be data dependent) satisfying condition (10). If we have a
data independent initializer, condition (10) can be weakened and we have the following result
on random initialization of the labels.
Theorem 6.1. Assume
kθ ∗ k ≥ Cσ

p

log(1/δ) + (d log(1/δ)(1 + d/n))1/4
21

for a sufficiently large constant C. Suppose we independently draw 3 log(1/δ) random initializers. Then Lloyd’s algorithm initialized by one of them can correctly classifies all the data
points with probability greater than 1 − δ.
The proof of Theorem 6.1 is given in appendix. Klusowski and Brinda [36] showed similar
√
results for EM algorithm under the condition that kθ ∗ k & σ d log d. We improve their
condition to kθ ∗ k & σ(d log d)1/4 by choosing δ = 1/d. Recently, Xu et al. [56] and Daskalakis
et al. [16] proved the global convergence of EM algorithm on the population likelihood
function. However, it is not clear what signal-to-noise ratio condition they need for the finite
sample analysis. Moreover, our results can be generalized to k-mixtures and the sub-Gaussian
case, but [36], [56] and [16] all restrict their analyses to the two mixture Gaussian case.
As for the general k mixture case, one popular way to initialize Lloyd’s algorithm is to
initialize k centers by sampling without replacement from the data. However, when k ≥
√
3, Jin et al. [30] show even for well-separated (∆ & dσ) spherical Gaussian mixtures,
EM algorithm, initialized by above strategy, convergences to bad critical points with high
probability. They prove it by constructing a bad case of two cluster of centers that are very
far away with each other. We expect similar results for Lloyd’s algorithm.

6.2

Convergence of the centers

As a byproduct of Theorem 3.2, we obtain the convergence rate of estimating Gaussian
mixture’s centers using Lloyd’s algorithm.
√
Theorem 6.2. Assume nα2 > Ck log n and rk ≥ C k for a sufficiently large constant C.
Then given any initializer satisfying


1
1
6
1
4
G0 <
−√
or Λ0 ≤ − √
(29)
2
rk λ
2
rk
with probability 1 − ν, we have


∆2
C2 σ 2 (d + log n)
2
+ σ exp −
for s ≥ 4 log n,
− θh k ≤
n
24σ 2
 4
for some constant C2 with probability greater than 1 − 2 exp − ∆
σ − n − ν.
(s)
max kθ̂h
h∈[k]

2

(30)

We show a linear convergence rate of estimating the centers. When the labels are known,
p
the minimax optimal rate of estimating the centers under Euclidean distance is σ d/n.
√
Thus, we obtain the optimal rate under the condition that nα2 & k log n, rk & k, d & log n
and ∆ & σ log(n/d).

7

Proofs

This section is devoted to prove Theorem 3.1. The proofs of other theorems are collected in
the appendix.
22

7.1

Preliminaries

We rewrite the symmetric two-mixture model (3) as
zi yi = θ ∗ + wi ,

i = 1, 2, · · · , n,

(31)

where wi = zi ξi ∼ N (0, σ 2 Id ). To facilitate our proof, we present four lemmas on the
concentration behavior of wi . Their proofs are given in Section D of the appendix.
P
Lemma 7.1. For any S ⊆ [n], define WS = i∈S wi . Then,
p
(32)
kWS k ≤ σ 2(n + 9d)|S| for all S ⊆ [n]
with probability greater than 1 − exp(−0.1n).
P
∗ 2
Lemma 7.2. Let w̄ = n1 ni=1 wi . Then we have w̄′ θ ∗ ≥ − kθ√nk and kw̄k2 ≤


∗ 2
with probability greater than 1 − 2 exp − kθ3σk2 .

3dσ2
n

+

kθ ∗ k2
n

P
Lemma 7.3. The maximum eigenvalue of ni=1 wi wi′ is upper bounded by 1.62(n + 4d)σ 2
with probability greater than 1 − exp(−0.1n).
Lemma 7.4. For any fixed i ∈ [n], S ⊆ [n], t > 0 and δ > 0, we have
*

+
p

 2 

3σ 2 (t |S| + d + log(1/δ)) 
1 X
t t
P
wi ,
wj ≥
≤ exp − min
,
+ δ.


|S|
|S|
4d 4
j∈S

Let E be the intersection of high probability events in Lemma 7.1, Lemma 7.2, Lemma

7.3 and label initialization condition (10). Then P{E c } ≤ ν + n−3 + exp −kθ ∗ k2 /(3σ 2 ) . For
the center initialization condition, we refer to the proof of Theorem 3.2 in the appendix for
more details. In the following, we deterministically analyze the error of estimating centers
and labels on event E.
Error of the centers We first decompose the error of estimating θ ∗ into two parts, where
the first part comes from the error of estimating labels zi , and the second part is due to the
stochastic error wi ,
θ̂

(s)

−θ

∗

=
(s)

Using (31) and the fact that ẑi
θ̂ (s) − θ ∗ = −
(s)

where ci = I{ẑi

n

n

n

i=1

i=1

i=1

1X
1X
1 X (s)
ẑi yi −
zi yi +
zi yi − θ ∗
n
n
n
(s)

− zi = −2I{ẑi

6= zi }zi , we obtain

n

n

i=1

i=1

2X
1X
(s)
I{ẑi 6= zi }(θ ∗ + wi ) +
wi = −2As θ ∗ − 2R + w̄
n
n

6= zi }, As =

1
n

Pn

i=1 ci

and R =

23

1
n

Pn

i=1 ci wi .

(33)

(s+1)

(s+1)

Error of the labels Now we analyze the error of ẑi
. The iteration of ẑi
is equivalent
to
D
E
D
E
(s+1)
ẑi
= argmin kgyi − θ̂ (s) k22 = argmax gyi , θ̂ (s) = argmax gzi θ ∗ + wi , θ̂ (s) ,
g∈{−1,1}

(s+1)

which gives us I{ẑi
D

θ ∗ + wi , θ̂ (s)

E

g∈{−1,1}

g∈{−1,1}

6= zi } = I

E
o
nD
θ ∗ + wi , θ̂ (s) ≤ 0 . From (33),

= hθ ∗ + wi , (1 − 2As )θ ∗ − 2R + w̄i

= (1 − 2As )kθ ∗ k2 − (2R − w̄)′ θ ∗ + hwi , θ ∗ − 2As θ ∗ + 2R − w̄i .

Applying Lemma 7.1 with S = {i ∈ [n], ci = 1}, we obtain
kRk ≤

σp
kθ ∗ k p
2(n + 9d)nAs ≤
2As ,
n
r

which, together with Lemma 7.2, implies

√
2 2As ∗ 2
1
(2R − w̄) θ ≤ 2kRkkθ k − w̄ θ ≤
kθ k + √ kθ ∗ k2 .
r
n
′ ∗

∗

′ ∗

Consequently, on event E, we have
(s+1)

I{ẑi


6= zi } ≤ I β0 kθ ∗ k2 + hwi , θ ∗ i + hwi , −2As θ ∗ + 2R − w̄i ≤ 0 ,

(34)

√

where β0 = 1 − 2As − 2 r2As − √1n ≥ 1 − 2As − 2r − √1n . Now we are ready to prove Theorem
3.1. In the following, we are going to prove (11) and (12) separately based on two different
decompositions of the RHS of (34).

7.2

Proof of (11) in Theorem 3.1

Using the fact that for any a, b ∈ R and c > 0, I{a + b ≤ 0} ≤ I{a ≤ c} + I{b ≤ −c} ≤ I{a ≤
2
c} + bc2 , we obtain
(s+1)
I{ẑi

(wi′ (2R − w̄ − 2As θ ∗ ))2
6 zi } ≤ I βkθ k ≤ − hwi , θ i +
=
,
δ2 kθ ∗ k4


∗ 2

∗

with β = β0 − δ and δ = 3.12/r. Taking an average of the equation above over i ∈ [n], we
have As+1 ≤ I1 + I2 , where
n

I1 =

1X 
I hwi , θ ∗ i ≤ −βkθ ∗ k2 ,
n

(35)

i=1

I2 =

n
X
2
1
wi′ (2R − w̄ − 2As θ ∗ ) .
2
∗
4
nδ kθ k
i=1

24

(36)

Upper Bound I1 Note that nAs only takes discrete values {1, 2, · · · , ⌊ n2 ⌋}. Let us define
√1
γa = 1 − 2a − 5.12
r − n and

Ti (a) = I hwi , θ ∗ i ≤ −γa kθ ∗ k2

⌊n⌋

for a ∈ D , { n1 , n2 , · · · , n2 }. For any fixed a, {Ti (a), i = 1, 2, · · · , n} are independent
Bernoulli random variables, then Hoeffding’s inequality implies
( n
)
r
4 log(n/2)
1X
1
P
≤ 4.
(Ti (a) − ETi (a)) ≥
n
n
n
i=1

Thus by union bound, we have
)
# r
(
" n
1
4 log(n/2)
1X
≤ 3.
(Ti (a) − ETi (a)) ≥
P max
a∈D n
n
n

(37)

i=1

Now it remains to upper bound the expectation of Ti (a). Since hwi , θ ∗ i is sub-Gaussian with
parameter σkθ ∗ k, Chernoff’s bound yields
 2 ∗ 2
γ kθ k
ETi (a) ≤ exp − a 2
.
(38)
2σ
Consider the event
(
E1 =

)
 2 ∗ 2 r
n
γa kθ k
4 log(n/2)
1X
,∀ a ∈ D .
Ti (a) ≤ exp −
+
n
2σ 2
n
i=1

Equation (37) together with (38) implies P{E1c } ≤ n−3 . Therefore,
! r
2 kθ ∗ k2
γA
4 log(n/2)
s
+
I1 ≤ exp −
2σ 2
n

(39)

with probability greater than 1 − n−3 .
Upper Bound I2

Lemma 7.3 implies
I2 ≤
≤

k2R − w̄ − 2As θ ∗ k4
λmax
nδ2 kθ ∗ k4

n
X
i=1

wi wi′

!

1.62(1 + 4d/n)σ 2
k2R − w̄ − 2As θ ∗ k2 ,
δ2 kθ ∗ k4

where λmax (B) stands for the maximum eigenvalue of the symmetric matrix B. By the
Cauchy-Schwarz inequality, Lemma 7.1 and Lemma 7.2, we have

k2R − w̄ − 2As θ ∗ k2 ≤ (1 + 1 + 4) 4kRk2 + kw̄k2 + A2s kθ ∗ k2


1
1
8
2
∗ 2
As + 2 + + As .
≤ 6kθ k
r2
r
n
25

Since δ = 3.12/r, we obtain

1
1
8
As + 2 + + A2s .
r2
r
n
≤ I1 + I2 . Combining (39) and (40), we have
! r
2 kθ ∗ k2
γA
8
1
4 log n
s
+ 2 As + 2 + A2s
+
As+1 ≤ exp −
2σ 2
n
r
r
I2 ≤

Recall that As+1

(40)

√

on event E ∩ E1 . Let τ = 12 − √1n − 2.56+r log r . From above inequality and the assumption that
A0 ≤ τ, it can be proved by induction
that As ≤ τ for all s ≥ 0 when r ≥ 32 and n ≥ 1600.
√
2 log r
for all s ≥ 0. Plugging this bound of γAs into above
Consequently, we have γAs ≥
r
inequality completes the proof.

7.3

Proof of (12) in Theorem 3.1

By (11) in Theorem 3.1 and the fact that As ≤ τ when r ≥ 32,
r
2
4 log n
As+1 ≤ 0.5As + 2 +
,
r
n
p
which implies As ≤ 4r −2 + 5 (log n)/n for all s ≥ log n. Then when s ≥ log n, (34) implies
(s+1)

I{ẑi
6= zi }
√



3 1 + 10 log n
∗ 2
∗
∗
√
kθ k + hwi , θ − 2As θ + 2R − w̄i ≤ 0 .
≤ I
1− −
r
n

To get a better upper bound for mis-clustering rate, we further decompose hwi , 2R − 2As θ ∗ − w̄i
as hwi , 2R − 2As θ ∗ i − hwi , w̄i. Following similar arguments as in the proof of Theorem 3.1(a),
we obtain As+1 ≤ J1 + J2 + J3 , where J1 is I1 defined in (35) with β = 1 − 6.72/r − (3 +
√
√
10 log n)/ n, J2 and J3 are defined as follows,
J2 =
J3 =

n
X
2
r2
wi′ (2R − 2As θ ∗ ) ,
∗
4
8.1nkθ k
i=1



n
1
2
1X
∗ 2
I
+√
kθ k ≤ − hwi , w̄i .
n
2r
n

(41)
(42)

i=1

∗ 2

From (38), when r ≥ 32 and n ≥ 40000, EJ1 ≤ exp(− kθ8σk2 ). Analogous to (40), we have


8
8
1
8
+
As ≤ As
J2 ≤ 2 As + A2s ≤
r
r2 r2
r

√

 ∗
∗
∗ 2
in Lemma 7.4
d + kθσ k kθσ k and δ = exp − kθ4σk2
on event E. Choosing S = [n], t =
yields


kθ ∗ k2
.
EJ3 ≤ 2 exp −
4σ 2
26

Combine the pieces and note that As+1 ≤ 1,
EAs+1 ≤ E [(J1 + J2 + J3 )I{E}] + E [As+1 I{E c }]


1
kθ ∗ k2
≤
EAs + 3 exp −
.
r
8σ 2

(43)

Then, by induction, we have
EAs ≤

1
r s−⌈log n⌉

kθ ∗ k2
+ 4 exp −
8σ 2




,

when r ≥ 32. Setting s = ⌈3 log n⌉ and applying Markov’s inequality, we obtain


1
kθ ∗ k2
− log t + 2 .
P {As ≥ t} ≤ 4 exp −
2
8σ
n t


∗ 2
∗ k2
If kθ8σk2 ≤ 2 log n, a choice of t = exp − kθ
gives us the desired result. Otherwise, since
16σ2

As only takes discrete values {0, n1 , n2 , · · · , 1}, choosing t = n1 , we have


5
1
1
≤ 4n exp(−2 log n) + ≤ .
P {As > 0} = P As ≥
n
n
n

Therefore, the proof of (12) is complete. When r → ∞, we can slightly modify the proof to
get the improved error exponent. For more details, we refer to the proof of Theorem 3.2 in
Section A.

A
A.1

Proof of sub-Gaussian mixtures
Preliminaries

Let us begin by introducing some notaion.
For any S o⊆ [n], define WS =
n
(s)
(s)
∗
that Tg = {i ∈ [n], zi = g} and Tg = i ∈ [n], ẑi = g , let us define

(s)

Then we have nh

P

i∈S

wi . Recall

n
o
(s)
(s)
(s)
Sgh = i ∈ [n], zi = g, ẑi = h = Tg∗ ∩ Th .

=

P

(s)
g∈[k] ngh

(s)
g∈[k] nhg .
(s)
(s)
(s)
ngh , nh and Sgh

and n∗h =

P

In the rest of the proof, we will

when there is no ambiguity. Also,
sometimes drop the upper index (s) of
we suppress the dependence of k by writing rk as r.
Analogous to Lemma 7.1 - Lemma 7.4, we have following technical lemmas on the concentration behavior of sub-Gaussian random vectors {wi }.
Lemma A.1.
kWS k ≤ σ

p

3(n + d)|S| for all S ⊆ [n].

with probability greater than 1 − exp(−0.3n).
27

(44)

Lemma A.2.

n
X

λmax

wi wi′

i=1

!

≤ 6σ 2 (n + d).

(45)

with probability greater than 1 − exp(−0.5n).
Lemma A.3. For any fixed i ∈ [n], S ⊆ [n], t > 0 and δ > 0,
*

+
p

 2 

2
X
3σ (t |S| + d + log(1/δ)) 
1
t t
P
wi ,
wj ≥
,
≤ exp − min
+ δ.


|S|
|S|
4d 4
j∈S

Lemma A.4.

kWTh∗ k ≤ 3σ

q

(d + log n)|Th∗ | for all h ∈ [k]

(46)

with probability greater than 1 − n−3 .

Lemma A.5. For any fixed θ1 , · · · , θk ∈ Rd and a > 0, we have
 2 2 q
X 
a ∆
2
∗
I akθh − θg k ≤ hwi , kθh − θg ki ≤ ng exp −
+ 5n∗g log n
2
2σ
∗

(47)

i∈Tg

for all g 6= h ∈ [k]2 with probability greater than 1 − n−3 .

A.2

Two Key Lemmas

The following two lemmas give the iterative relationship between the error of estimating
centers and the error of estimating labels. Let E be the intersection of high probability
events in Lemma A.1, Lemma A.2 Lemma A.4, Lemma A.5 and the initialization condition
(14). Then we have P{E c } ≤ 3n−3 + ν. In the rest part of the proof, if not otherwise stated,
we all condition on the event E and the following analysis are deterministic.
Lemma A.6. On event E, if Gs ≤ 12 , then we have

 p
3
3
Λs ≤ + min
kGs + 2Gs Λs−1 , λGs .
r
r
1−ǫ
2

(48)

and r ≥ 36ǫ−2 , then

2 r
5k log n
2
28
Gs+1 ≤ 4 2 + 2 Λs +
.
(49)
ǫ r
ǫ r
α2 n
1 P
Proof of Lemma A.6. For any B ⊆ [n], define ȲB = |B|
i∈B yi . The error of estimated
centers at step s can be written as
1 X X
1 X
(s)
(yi − θh ) +
(yi − θh )
θ̂h − θh =
nh
nh
Lemma A.7. On event E, if Λs ≤

i∈Shh

=

a6=h i∈Sah

X nah

1 X
wi +
ȲSah − θh
nh
nh
i∈Shh

a6=h

28

(s−1)

(s−1)

k ≤ kyi − θ̂a
k for any i ∈ Sah .
According to our label update step, we have kyi − θ̂h
(s−1)
(s−1)
than θ̂a
, so is the average of {yi , i ∈ Sah }.
This means for any i ∈ Sah , yi is closer to θ̂h
Thus, we have
(s−1)
k ≤ kȲSah − θ̂a(s−1) k.
kȲSah − θ̂h
Consequently, triangle inequality gives us
(s−1)

ȲSah − θh ≤ ȲSah − θa + kθ̂a(s−1) − θa k + kθ̂h

− θh k,

which, combined with Lemma A.1 and the definition of Λs−1 , yields
p
ȲSah − θh ≤ σ 3(n + d)/nah + 2Λs−1 ∆.

Taking a weighted sum over a 6= h ∈ [k], we get
p
X nah
X nah
σ 3(n + d) X √
ȲSah − θh
≤
nah + 2Λs−1 ∆
nh
nh
nh
a6=h
a6=h
a6=h
p
σ 3(n + d) p
≤
(k − 1)Gs + 2Gs Λs−1 ∆,
√
nh
P
where the Last inequality is due to Cauchy-Schwartz and the fact that a6=h nah ≤ Gs nh .
P
Note that WShh = WTh∗ − a6=h WSha . Triangle inequality and Lemma A.5 imply
p
p
p
p
kWShh k ≤ 3σ d + log n n∗h + σ 3(n + d) n∗h − nhh .
Since Gs ≤ 21 , we have

nh ≥ nhh ≥ n∗h (1 − Gs ) ≥

1
1 ∗
nh ≥ αn.
2
2

(50)

Combining the pieces, we obtain
(s)
θ̂h

r
d + log n
k(n + d)
+ 3σ
Gs + 2Gs Λs−1 ∆
≤ 3σ
αn
αn

p
3
(1 + kGs ) + 2Gs Λs−1 ∆.
≤
r
r

− θh

(51)
(s)

Therefore, we get the first term in (48). To prove the second term, we decompose θ̂h
differently.
(s)

θ̂h

=

n
n
o
1 X
(s)
(θzi + wi ) I ẑi = h
nh

=

1
nh

=

k
X
a=1

i=1
k X
n
X
a=1 i=1

n
o
1 X
(s)
wi
θa I zi = a, ẑi = h +
nh
i∈Th

1
nah
θa +
WTh .
nh
nh
29

(52)

(s)

Then, the error of θ̂h can be upper bounded as
(s)

θ̂h − θh

k
X
nah

=

a=1

nh

(θa − θh ) +

X nah
1
1
WTh ≤
(θa − θh ) +
WTh .
nh
nh
nh
a6=h

By triangle inequality,
X nah

a6=h

nh

(θa − θh ) ≤

X nah

a6=h

nh

kθa − θh k ≤ λ∆

X nah

a6=h

nh

≤ λ∆Gs .

This, together with Lemma A.1 and (50), implies
s


3
3(n + d)
(s)
≤ λGs +
θ̂h − θh ≤ λ∆Gs + σ
∆
nh
r

(53)

(54)

for all h ∈ [k]. The proof is complete.
Proof of Lemma A.7. For any g =
6 h ∈ [k] × [k],
o
o
n
n
(s)
(s+1)
=h
≤ I kθg + wi − θ̂h k2 ≤ kθg + wi − θ̂g(s)k2
I zi = g, ẑi
Eo
D
n
(s)
(s)
= I kθg − θ̂h k2 − kθg − θ̂g(s) k2 ≤ 2 wi , θ̂h − θ̂g(s) .

(55)

Triangle inequality implies


2
(s)
(s)
kθg − θ̂h k2 ≥ kθg − θh k − kθh − θ̂h k ≥ (1 − Λs )2 kθg − θh k2 .

Using the fact that (1 − x)2 − y 2 ≥ (1 − x − y)2 when y(1 − x − y) ≥ 0, we obtain
(s)

kθg − θ̂h k2 − kθg − θ̂g(s)k2 = (1 − 2Λs )2 kθg − θh k2 ≥ ǫ2 kθg − θh k2 .

(56)

(s)

Denote by ∆h = θ̂h − θh for h ∈ [k]. Then,
n
o
(s+1)
I zi = g, ẑi
=h

≤ I ǫ2 kθg − θh k2 ≤ 2 hwi , θh − θg + ∆h − ∆g i
 2

 2

ǫ
ǫ 2
2
≤ I
kθg − θh k ≤ 2 hwi , θh − θg i + I
∆ ≤ 2 hwi , ∆h − ∆g i .
2
2
Taking a sum over i ∈ Tg∗ and using Markov’s inequality on the second term, we obtain
(s+1)

ngh

≤

 X
X  ǫ2
2
16
I
kθg − θh k2 ≤ hwi , θh − θg i +
wi′ (∆h − ∆g )
4
4
4
ǫ ∆
∗
∗

i∈Tg

i∈Tg

30

(57)

n 2
o
Note that I ǫ4 kθg − θh k2 ≤ hwi , θh − θg i are independent Bernoulli random variables. By
Lemma A.5, the first term in RHS of (57) can be upper bounded by
 4 2 q
ǫ ∆
∗
(58)
ng exp −
+ 5n∗g log n.
32σ 2
By Lemma A.2, the second term in RHS of (57) can be upper bounded by
X 16
2 96(n∗g + d)σ 2
′
≤
w
(∆
−
∆
)
k∆g − ∆h k2 .
g
h
i
ǫ4 ∆ 4
ǫ4 ∆ 4
∗

(59)

i∈Tg

Combining (57), (58) and (59) and using the fact that k∆g − ∆h k2 ≤ 4Λ2s ∆2 , we get
(s+1)
ngh

≤

n∗g exp



ǫ4 ∆ 2
−
32σ 2



+

q

5n∗g log n +

384(n∗g + d)σ 2 2
Λs .
ǫ4 ∆ 2

Consequently,
max
g∈[k]

X n(s+1)
gh
h6=g

n∗g

r
 4 2
5 log n
ǫ ∆
384
≤ k exp −
+ 4 2 Λ2s .
+k
2
32σ
αn
ǫ r

(60)

Since Λs ≤ 1/2 and r ≥ 20ǫ−2 , the RHS of (60) is smaller that 1/2 when αn ≥ 32k2 log n.
Thus,
1
1
(s+1)
(s+1)
≥ nhh ≥ n∗h ≥ αn
nh
2
2
for all h ∈ [k] and we have
max
h∈[k]

X n(s+1)
gh
(s+1)

g6=h

nh

 4 2 r
ǫ ∆
768
5k log n
2
+ 4 2 Λ2s ,
+
≤ exp −
2
2
α
32σ
α n
ǫ r

(61)

which, together with (60), implies
Gs+1

 r
ǫ4 ∆ 2
768
5k log n
≤ exp −
+ log(2/α) +
+ 4 2 Λ2s
2
2
32σ
α n
ǫ r


Under the assumptions that ǫ4 α∆2 /σ 2 ≥ r 2 ǫ4 ≥ 36, we have the desired result (49).

A.3

Proof of Theorem 3.2

Proof of (15) in Theorem 3.2. From Lemma A.6, a necessary condition for Λ0 ≤
G0 ≤

( 12

√6 ) 1 .
r λ

√7
r

1
2

−

√4
r

is

in Lemma A.7, we have G1 ≤ 0.35. Plugging it into Lemma
√
A.6 gives us Λ1 ≤ 0.4, under the assumption that r ≥ 16 k. Then it can be easily proved
by induction that Gs ≤ 0.35 and Λs ≤ 0.4 for all s ≥ 1. Consequently, Lemma A.6 yields
−

Setting ǫ =

Λs ≤

3 3p
1
+
kGs + Gs ≤ + Gs
r r
2
31

which, combined with (49), implies
Gs+1

C
C
≤ 2+ 2
r
r



1
+ 2Gs + G2s
4



+

r

2C
3C
5k log n
≤ 2 + 2 Gs +
2
α n
r
r

r

5k log n
α2 n

for some constant C. Here we have chosen ǫ = 1/5 in Lemma 49 to get the first inequality.
Proof of (16) in Theorem 3.2. From the proof of Lemma A.6, the error of estimating θh at
(s)
iteration s can be written as θ̂h − θh = n1h WTh∗ + uh , with

 p
p
3
kGs + Gs ∆ ≤ Gs ∆
kuh k ≤
r

(62)

In addition, by Lemma A.6 and Lemma A.7, there is a constant C1 such that
C1 C1
3 p
+
Λs−1 + 0.7Λs−1 +
Λs ≤ + Gs + 2Gs Λs−1 ≤
r
r
r



C1 k log n
α2 n

1/4

for all s ≥ 1. Therefore, when r is large enough, we have
Λ ≤ C2 r

−1

+ C2



k log n
α2 n

1/4

for all s ≥ log n. Then by (56), we have
n
o

(s+1)
I zi = g, ẑi
= h ≤ I β1 kθg − θh k2 ≤ 2 hwi , θh − θg + ∆h − ∆g i

where (1 − 2Λs )2 ≥ β1 := 1 − 4C2 r −1 − 4C2



k log n
α2 n

1/4

.

In order to prove that As attains convergence rates (16), we first upper bound the expectation of As and then derive the high probability bound using Markov’s inequality. Similar
to the two-mixture case, we need to upper bound the inner product hwi , ∆h − ∆g i more carefully. Note that {Th∗ , h ∈ [k]} are deterministic sets, we could use concentration equalities to
upper bound WTh∗ and uh parts separately.
n
o
(s+1)
Let vh = n1h WTh∗ for h ∈ [k] and we decompose I zi = g, ẑi
= h into three terms.
n
o

(s+1)
I zi = g, ẑi
=h
≤ I βkθg − θh k2 ≤ 2 hwi , θh − θg i

+I β2 ∆2 ≤ 2 hwi , uh − ug i

+I β4 ∆2 ≤ 2 hwi , vh − vg i ,

where β2 and β4 will be specified later and β = β1 − β2 − β4 . Taking a sum over h ∈ [k] and
i ∈ [n], we obtain
EAs+1 ≤ EJ1 + EJ2 + EJ3

32

with
J1 =

n
X 1X

I βkθzi − θh k2 ≤ 2 hwi , θh − θzi i
n

h∈[k]

J2 =

i=1

n
X 1X

I β2 ∆2 ≤ 2 hwi , uh − uzi i .
n

h∈[k]

J3 =

(63)

(64)

i=1

n
X 1X

I β4 ∆2 ≤ 2 hwi , vzi − vh i .
n

h∈[k]

(65)

i=1

Let us first consider the expectation of J1 . Using Chernoff’s bound, we have





β 2 kθh − θg k2
β 2 ∆2
2
≤ exp −
.
P βkθg − θh k ≤ 2 hwi , θh − θg i ≤ exp −
8σ 2
8σ 2

Thus,

β 2 ∆2
EJ1 ≤ k exp −
8σ 2


with γ = β 2 −

8σ2 log k
∆2



γ∆2
= exp − 2
8σ




,

≥ β 2 − 8/r 2 .

We use Markov Inequality to upper bound J2 . Markov’s inequality and Lemma A.2 give
us
n

1X 
I β2 ∆2 ≤ 2 hwi , uh − uzi i
n
i=1

≤
≤

X X
2
4
wi′ (uh − ug )
2
4
nβ2 ∆
∗
g∈[k] i∈Tg

24σ 2 X ∗
(ng + d)kuh − ug k2 .
nβ22 ∆4
g∈[k]

(62) implies
√
96σ 2 k(n + kd)
12 k
96σ 2 Gs X X ∗
(ng + d) ≤
As =
J2 ≤
As .
r
nβ22 ∆2
αnβ22 ∆2
h∈[k] g∈[k]

Here the second inequality is due to the fact that Gs ≤ As /α. And we choose β2 =
in the last equality.
Finally, we upper bound the expectation of J3 . Given zi = g, we have

P β4 ∆2 ≤ 2 hwi , vg − vh i




β4 2
β4 2
∆ ≤ hwi , vg i + P − ∆ ≥ hwi , vh i
≤ P
4
4






β4 2
β4 2
1
1
∗
≤ P
+ P − ∆ ≥ wi , ∗ WTh
∆ ≤ wi , ∗ WTg∗
8
ng
8
nh
33

p

8k/r

Choosing t = max{

√

d∆ ∆2
σ , σ2 },



∆2
in Lemma A.3, and
δ = exp − 4σ
2

!
√
3 max{ dσ∆, ∆2 } 3σ 2 d + ∆2
√
,
+
αn
αn

64
8
β4 =
≥ 2
r
∆


we obtain P β4 ∆2 ≤ 2 hwi , vg − vh i ≤ 2 exp(−∆2 /(4σ 2 )), where we have used the assumption that n∗g ≥ αn and αn ≥ 36r 2 . Thus,


∆2
EJ3 ≤ 2k exp − 2 ,
σ
Combining the pieces, we have
EAs+1 ≤ E [J1 ] + E [J2 I{E}] + E [J3 ] + P{E c }
√




∆2
12 k
γ∆2
EAs + 2k exp − 2 ,
≤ exp − 2 +
8σ
r
σ
p
with γ = (β1 − 8k/r − 64/r)2 − 8/r 2 = 1 − o(1). Here only prove the case that r → ∞. For
the finite case, all the o(1) in the following proof can be substituted by a small constant.




1
∆2
2
3
∆2
EAs ≤ s−⌈log r⌉ + 2 exp −(1 − η) 2 + 3 ≤ 2 exp −(1 − η) 2 + 3
8σ
n
8σ
n
2
when s ≥ 4 log n. By Markov’s inequality, for any t > 0,


1
2
∆2
3
P {As ≥ t} ≤ EAs ≤ exp −(1 − η) 2 + 3 .
t
t
8σ
n t


8σ ∆2
∆2
and we have
If (1 − η) 8σ
2 ≤ 2 log n, choose t = exp −(1 − η − ∆ ) 8σ 2

(66)






8σ ∆2
∆
4
P As ≥ exp −(1 − η −
)
≤ + 2 exp −
.
∆ 8σ 2
n
σ

Otherwise, since As only takes discrete values of {0, n1 , · · · , 1}, choosing t =
to


1
4
3
P {As > 0} = P As ≥
≤ 2n exp(−2 log n) + 2 ≤ .
n
n
n

1
n

in (66) leads

The proof is complete.

A.4

Proof of Lower Bounds

The key difficulty in proving the lower bound is to deal with infinitum of all label permutations. Here we adapt the proof idea from [23]. We define a subset of the parameter space,
in which a large portion of the labels in each cluster are fixed. Then any permutation other
than identity gives us bigger mis-clustering rate.

34

P
Proof of Theorem 3.3. For any z ∈ [k]n , let us define nu (z) = ni=1 I{zi = u}. Let z ∗ ∈ [k]n
satisfying n1 (z ∗ ) ≤ n2 (z ∗ ) ≤ · · · , ≤ nk (z ∗ ) with n1 (z ∗ ) = n2 (z ∗ ) = ⌊αn⌋. It is easy to check
the existence of z ∗ . For each u ∈ [k], we choose a subset of {i : z ∗ (i) = u} with cardinality
k
⌈nu (z ∗ ) − αn
4k ⌉, denoted by Tu . Define T = ∪u=1 Tu and
Z ∗ = {z ∈ Z0 , zi = zi∗ for all i ∈ T }
P
A key observation is that for any z 6= z̃ ∈ Z ∗ , we have n1 ni=1 I{zi 6= z̃i } ≤

(67)
k αn
n 4k

=

α
4

and

n

1X
1
αn
α
I{π(zi ) 6= z̃i } ≥ (αn −
)≥
n
n
4k
2

(68)

i=1

P
for all π ∈ Sk 6= Ik . Thus, ℓ(z, z̃) = n1 ni=1 I{zi 6= z̃i } for all z, z̃ ∈ Z ∗ . Then following the
same arguments in the proof of Theorem 2 in [23], we can obtain


α 1 X 1
inf sup Eℓ(ẑ, z) ≥
inf (P1 {ẑi = 2} + P2 {ẑi = 1})
(69)
ẑ z∈Z
6 |T c |
2k2 ẑi
c
i∈T

Here Pt , t ∈ {1, 2} denote the probability distribution of our data given zi = t. By NeymanPearson Lemma, the infimum of the right hand side of (69) is achieved by the likelihood ratio
test
ẑi = argmin kyi − θg k2 .
g∈{1,2}

Thus,
inf
ẑi





1
1
P1 {ẑi = 2} + P2 {ẑi = 1}
= P kθ1 + wi − θ2 k2 ≤ kwi k2
2
2

= P kθ1 − θ2 k2 ≤ 2 hwi , θ1 − θ2 i .

Since wij , j ∈ [d] are independent N (0, σ 2 ), hwi , θ1 − θ2 i ∼ N (0, σ 2 kθ1 − θ2 k2 ). Let Φ(t) be
the cumulative function of N (0, 1) random variable. By calculating the derivatives, it can be
easily proved that
Z ∞
2
t
2
1
1
e−t /2 .
e−x /2 dx ≥ √
1 − Φ(t) = √
2+1
t
2π t
2π
Then when kθ1 − θ2 k ≥ σ, we have




1
kθ1 − θ2 k2
σ
2
P hwi , θ1 − θ2 i ≥ kθ1 − θ2 k ≥ √
exp −
.
2
8σ 2
2πkθ1 − θ2 k
Consequently,
∆2
2k∆
inf sup Eℓ(ẑ, z) ≥ exp − 2 − 2 log
ẑ z∈Z
8σ
ασ


The proof is complete.
35



.

A.5

A counterexample

Now we give a counterexample to show that the initialization condition in Theorem 3.2 is
almost necessary. Consider a noiseless case with 6 equal size clusters, as showed in Figure
3. Suppose kθi − θj k = ∆ for i 6= j ∈ [3] and kθi+3 − θi k = λ∆ for i ∈ [3]. We are given
an initializer that mixes cluster i and cluster i + 3, i ∈ [3], with m/(2λ) data points from
cluster i + 3 and m − m/(2λ) data points from cluster i for some integer m. Consequently,
(0)
the estimated (initialized) center θ̂i lies in the middle of two true centers. For the next
(0)
label update step, since there are ties, we may assign half points at θ1 to cluster θ̂1 and
(0)
another half to cluster θ̂3 and the estimated centers remain the same as before. Therefore,
(0)
{θ̂i , i ∈ [3]} is a stationary point for Lloyd’s algorithm and it may not convergence to the
true centers. We would like to note that this counterexample is a worst case in theoretical
analysis, which may not happen in practice.

Figure 3: A counter example showing that Lloyd’s algorithm may not converge when G0 =
or Λ0 = 21 .

A.6

1
2λ

Proofs of Random Initialization

Proof of Theorem 6.1. For a data independent initializer, we have the following result.
Lemma A.8. For any δ < 1/4, if we have a data independent initializer satisfying
r
√ !
18 d
σ
2 log(1/δ)
1
1+
A0 ≤ − ∗
2 kθ k
n
r

(70)

the conclusion of Theorem 3.1 continues to hold with probability greater than 1 − δ.
Given Lemma A.8, let us study the condition under which completely random initialization falls into the basin of attraction (70). When we randomly assign initial labels, we have
|A0 − 1/2| ≥ n−1/2 with probability greater than 0.3. We can boost this probability to 1− 1/δ
by independently drawing 3 log(1/δ) different random initializers. Recall our loss function
36

(4). We can assume A0 < 1/2 − n−1/2 . Otherwise, we may flip all the labels in the initial
step. Combining this with condition (70), we have the following corollary.

θ̂ (1)

Now it remains to prove Lemma A.8. We focus on the first iteration. By (33), we have
P
− θ ∗ = −2A0 θ ∗ + R with R = n1 ni=1 (1 − 2ci )wi . Then
(1)

I{ẑi

6= zi } = I {hθ ∗ + wi , (1 − 2A0 )θ ∗ + Ri ≤ 0}




1
hR, θ ∗ + wi i
1
≤ I hwi , θ ∗ i ≤ − kθ ∗ k2 + I
≤ − kθ ∗ k2 .
2
1 − 2A0
2

Since the initializer is data independent, {ci } ⊆ {−1, 1}n are independent of {wi }. Chernoff’s
bound implies
p
n
σkθ ∗ k 2 log(1/δ)
1X
′ ∗
∗
√
(1 − 2cj )wi θ ≥ −
hθ , Ri =
n
n
j=1

with probability greater than 1 − δ. Consequently, we obtain


n
n
1 ∗ 2
1X
1X 
∗
I hwi , θ i ≤ − kθ k +
I hR, wi i ≤ −βkθ ∗ k2 ,
A1 ≤
n
2
n
i=1

(71)

i=1

√
σ 2 log(1/δ)
where β = 12 − A0 − √nkθ∗ k . By Chernoff’s bound and Hoeffding’s inequality, the first
term in the RHS of (71) can be upper bounded by
 r

kθ ∗ k2
log(1/δ)
+
exp −
2
8σ
n
with probability greater than 1 − δ. Markov’s inequality and proof of Lemma 7.3 give us
n
n
X
1X 
9σ 2 (1 + d/n)
1
′
2
∗ 2
(R
w
)
≤
I hR, wi i ≤ −βkθ k ≤
kRk2 .
i
n
nβ 2 kθ ∗ k4
β 2 kθ ∗ k4
i=1

i=1

with probability greater than 1 − exp(−n/2). From (80), we have kRk2 ≤ 4σ
probability greater than 1 − δ. Combining the pieces, we obtain

 r
kθ ∗ k2
log(1/δ) 36σ 4 (1 + d/n) d log(1/δ)
A1 ≤ exp −
+
+
8σ 2
n
β 2 kθ ∗ k4
n

2 d log(1/δ)

n

with

(72)

with probability greater than 1 − 3δ. Under the condition that
!
r
p
18σ d(1 + d/n)
1
σ
2 log(1/δ)
A0 ≤ − ∗
1+
,
2 kθ k
n
kθ ∗ k
we have A1 ≤ 14 with probability greater than 1 − 3δ. Then using the result of Theorem 3.1,
the proof is complete.

37

B
B.1

Proofs of Community Detection
Proof of Theorem 4.1
(s)

We use the same notation as in the Gaussian case. Let Tg∗ be the true cluster g and Tg be
the estimated cluster g at iteration s. And we will drop the dependency of s when there is
no ambiguity in the context. As is the mis-clustering rate at iteration s and Gs is the groupwise mis-clustering rate at iteration s. Similar to the analysis of Gaussian mixture
model,
√
C0 aβk log(βk)
1
the proof consists of two steps. We first prove that given G0 ≤ 2 − ǫ0 with ǫ0 ≥
a−b
for a sufficiently large constant C0 , we have Gs ≤ 13 for all 1 ≤ s ≤ ⌊3 log n⌋ with high
probability. Then we prove the mis-clustering rate As geometrically decays to the rate (22).
Given zi = g, we decompose (20) as follows.
n
1 X
(s)
Aij I{ẑj = h} =
nh
j=1

=
=

k
n
n
1 X
1 XX
(s)
(s)
Bgl I{ẑj = h, zj = l}
wij I{ẑj = h} +
nh
nh
j=1
j=1 l=1
P
X
1
a ngh
b l6=g nlh
wij +
+
nh
n nh
n
nh
j∈Th


ngh
a ngh
b
1
Wi (Th ) +
+
1−
nh
n nh
n
nh

Consequently, the error of estimating zi at iteration s + 1 can be upper bounded as
(s+1)

≤
≤
≤
≤

I{ẑi
6= g, zi = g}


n
n


X
X
(s)
(s)
Aij I{ẑj = g}, zi = g
I max
Aij I{ẑj = h} ≥

 h6=g
j=1
j=1



ngh
1
1
a − b ngg
I max Wi (Th ) − Wi (Tg ) ≥
− max
h6=g nh
h6=g nh
ng
n
ng


1
1
2ǫ(a − b)
I max Wi (Th ) − Wi (Tg ) ≥
h6=g nh
ng
n




ǫ(a − b)
1
ǫ(a − b)
1
Wi (Tg ) ≤ −
I max Wi (Th ) ≥
+I
h6=g nh
n
ng
n

(73)

where the second inequality is due to our induction assumption that ngg ≥ (1 − 1/2 + ǫ)ng
and ngh ≤ (1/2 − ǫ)nh . Union bound implies
(s+1)
I{ẑi





k
X
ǫ(a − b)n∗zi
ǫ(a − b)nh
6 zi } ≤
=
I Wi (Th ) ≥
+ I Wi (Tzi ) ≤ −
.
n
n

(74)

h=1

Now we give two upper bounds on the sum of indicator variables uniformly over {Th }.
There proofs are deferred to Section B.2.

38

Lemma B.1. Given T ∗ ⊆ [n] with cardinality greater than n/(βk) and ǫ > 0 such that
ǫ2 (a − b)2 ≥ C0 βak log(n/m) for a sufficiently large constant C0 . There is an universal
constant C such that with probability greater than 1 − n−2 ,


 2

m
X
X
ǫ (a − b)2
|T ∆T ∗ |
ǫ(a − b)|T | 
≤ m exp −
(75)
+
I
wij ≥


n
Cβak
5
i=1

j∈T

holds for all T ⊆ [n] with |T ∆T ∗ | ≤ 21 |T ∗ |.

Lemma B.2. Given T ∗ ⊆ [n] with cardinality greater than n/(βk) and ǫ > 0 such that
ǫ2 (a − b)2 ≥ C0 β 2 ak 2 log(n/m) for a sufficiently large constant C0 . There is an universal
constant C such that with probability greater than 1 − n−2 ,


 2

m
X
X
ǫ (a − b)2
|T ∆T ∗ |
ǫ(a − b)|T | 
(76)
≤ m exp −
+
I
wij ≥


n
Cβ 2 ak 2
4βk
i=1

j∈T

holds for all T ⊆ [n] with |T ∆T ∗ | ≤ 21 |T ∗ |.

P
In Lemma B.1 and Lemma B.2, we present results for the upper tail m
i=1 I {Wi (Th ) ≥ ǫ(a − b)|T |/n}.
Pm
By slightly modifying the proof, the same results hold for the lower tail i=1 I {Wi (Th ) ≤ −ǫ(a − b)|T |/n}.
Taking an sum over i ∈ Tg∗ in (74) and using Lemma B.2, we obtain
X
(s+1)
I{ẑi
6= zi }
i∈Tg∗

 X 

XX 
ǫ(a − b)nh
ǫ(a − b)nh
≤
I Wi (Th ) ≥
+
I Wi (Tg ) ≤ −
k
k
∗
∗
h6=g i∈Tg

i∈Tg


ǫ2 (a − b)2
∗
≤ kng exp −
+
Cβ 2 ak2
 2

ǫ (a − b)2
∗
≤ kng exp −
+
Cβ 2 ak2


k
1 X
|Th ∆Th∗ |
4βk
h=1

n
As .
2βk

with probability greater than 1 − kn−2 . Consequently,

 2
1 X (s+1)
1
ǫ (a − b)2
+ Gs .
ngh ≤ k exp −
∗
2
2
ng
Cβ ak
2
h6=g

with probability greater than 1 − kn−2 . Here we have used the fact that As ≤ Gs . Using
P
(s+1)
similar arguments, we obtain the same high probability upper bound for n1h g6=h ngh .
Therefore, we have
 2

ǫ (a − b)2
1
Gs+1 ≤ k exp −
+ Gs
2
2
Cβ ak
2
√

aβk
, we
with probability greater than 1−k2 n−2 . When s = 0, we choose ǫ = ǫ0 . Since ǫ0 ≥ C0a−b
have G1 ≤ 1/12 + G0 /2 ≤ 1/3 when C0 ≥ (C + 1) log 12. For s ≥ 1, we choose ǫ = 1/3 and it

39

is straight forward to prove by induction that Gs+1 ≤ 1/12+1/6 ≤ 1/3 for all 0 ≤ s ≤ 3 log n,
with probability greater than 1 − n−1 , provided n ≥ 3k2 log n.
Now we are ready to give the convergence rate of As . Since Gs ≤ 1/3 for all s ∈ [1, 3 log n],
(74) holds for all s ∈ [1, 3 log n]. Taking average over i ∈ [n] in (74) and using Lemma B.1
with ǫ = 1/3, we obtain
As+1 ≤





k
k
n
a−b
a−b
1XX
1 XX
I Wi (Th ) ≥
I Wi (Tg ) ≤ −
+
n
6k
n g=1
6k
∗
h=1 i=1

i∈Tg


k
2 X
(a − b)2
+
|Th ∆Th∗ |
≤ 2k exp −
Cβak 2
5n
h=1


2
(a − b)
4
≤ exp −
+ As
2Cβak 2
5


with probability greater than 1 − kn−2 , where C is some universal constant.

B.2

Proof of Lemma B.1 and Lemma B.2

For any fixed T with |T ∆T ∗ | ≤ γn/k, let


X
ǫ(a − b)|T | 
.
ai (T ) = I
wij ≥


n
j∈T

By Bernstein’s inequality, the success probability of ai (T ) is upper bounded by





X
ǫ(a − b)|T | 
ǫ(a − b)|T | ǫ2 (a − b)2 |T |
P
wij ≥
,
≤ exp −C1 min


n
n
an
j∈T


C2 ǫ2 (a − b)2
≤ exp −
,
βak
for some universal constant C1 and C2 . Here the last inequality is due to the fact that
|T | ≥ |T ∗ | − |T ∆T ∗ | ≥ n/(2βk). Note that {ai (T ), i ∈ [n]} are independent Bernoulli
random variables. By Bennett’s inequality (see Lemma 9 in [22]), for any t > 0,
)
(m



X
t
ai (T ) ≥ mp + t
≤ exp t − (mp + t) log 1 +
P
mp
i=1



t
,
(77)
≤ exp −t log
emp
where



n
1X
C2 ǫ2 (a − b)2
p=
Eai (T ) ≤ exp −
, exp (−R) .
n
βak
i=1

40

To prove (75) holds for all T ⊆ [n] with |T ∆T ∗ | ≤ γn
k , we need a chaining argument. For
∗
s = m exp(−R/30), define D0 = {T ⊆ [n], |T ∆T | ≤ s} and

Dr = T ⊆ [n], 2r−1 s ≤ |T ∆T ∗ | ≤ 2r s .
Then D ⊆ ∪ur=0 Dr with u = ⌈log n⌉. Union bound implies
(
)
m
X
1
ai (T ) ≥ mp + |T ∆T ∗ | + s
A , P ∃ T ∈ D, s.t.
5
i=1
(
)
m
u
X
X
1
ai (T ) ≥ mp + |T ∆T ∗ | + s
P ∃ T ∈ Dr , s.t.
≤
5
r=0
i=1
(
)
m
u
X
X
1
∗
ai (T ) ≥ mp + |T ∆T | + s
|Dr | max P
≤
T ∈Dr
5
r=0

i=1

S
Note that T = T ∗ \ (T ∗ ∩ T c ) (T ∩ (T ∗ )c ). Given T ∗ and |T ∆T ∗ | = u, there are at most

n u
u 2 possible choices of T . Thus, the cardinality of Dr is upper bounded by




2r s  
2r s 

 en 
X
2en
n u X 2en u
r
r
r
.
≤ 2 s exp 2 s log
≤
exp
2
s
log
2 ≤
u
2r s
2r−2 s
u
u=0

u=0

Combining this with (77), we obtain

(m
X

and

1
ai (T ) ≥ mp + |T ∆T ∗ | + s
|D0 | max P
T ∈D0
5
i=1





s
2en
≤ exp −s log
+ s log
emp
s


Rs
≤ exp −
2

)

(m
X

)
1
∗
ai (T ) ≥ mp + |T ∆T | + s
|Dr | max P
T ∈Dr
5
i=1
 r
 r 
 en 
2
2 s
r
≤ exp − s log
+ 2 s log r−2
10
10emp
2 s
 r 
2 Rs
≤ exp −
20
2

n
for all r ≥ 1, where we have used the assumption that (a−b)
βak ≥ C log m for a sufficiently large
constant C. Thus,

 r 
 X


u
Rs
2 Rs
Rs
A ≤ exp −
exp −
+
≤ 3 exp −
.
2
20
10
r=1


2 (a−b)2
n
≥ 43 R+30 log m
. Consequently, A ≤ 3 exp −3s log ns ≤
By our assumption, we have R = C2 ǫ βak
n−2 . Therefore, the proof Lemma B.1 is complete.
Lemma B.2 follows from almost identical

R
arguments except that we choose s = m exp − 8βk .
41

C
C.1

Proofs of Crowdsourcing
Proof of Theorem 4.2

Suppose zj = g, then Chernoff bound gives us
)
(
m
m
o
n
X
X
(0)
I{Xij = g}
I{Xij = h} ≥
= P ∃h 6= g,
P ẑj 6= zj
i=1

i=1

X

≤

P

h6=g

(m
X

m
XY

≤

h6=g i=1

i=1

)

(I{Xij = h} − I{Xij = g}) > 0

E exp (λI{Xij = h} − λI{Xij = g}) .

for all λ > 0. Since Xij is from a multinomial distribution (πig1 , πig2 , · · · , πigk ), we have
m
Y

=
(a)
≤
=

E exp (λI{Xij = h} − λI{Xij = g})

i=1
m 
Y
i=1
m
Y
i=1


πigh eλ + πigg e−λ + (1 − πigh − πigg )



exp πigh eλ + πigg e−λ − πigh − πigg



exp Aeλ + Be−λ − A − B ,

P
Pm
x
where A = m
i=1 πigh and B =
i=1 πigg . Here inequality (a) is due to the fact that 1+x ≤ e
for all x. Choosing λ = 0.5 log(B/A) yields
 v
v
2 
u
u
m
m
o
n
X
uX
uX


(0)
πigg − t
πigh  
≤
P ẑj 6= zj
exp − t
i=1

h6=g

i=1

≤ k exp (−V (π))

2
pP
pPm
m
π
−
π
where V (π) = minh6=g
igg
i=1
i=1 igh . Then using Markov’s inequality, we
have


n
1 X

(0)
P
I{ẑj 6= zj } ≥ t ≤ exp (−V (π) + log k − log t)
n

j=1

for all t > 0. A choice of t = 1/(4k) completes the proof.

42

C.2

Proof of Corollary 4.2

Let us first calculate the sub-Gaussian parameter of wi . For any λ > 0 and a = (a11 , a12 , · · · , amk ) ∈
Rmk with kak = 1,
!
k
m
X
Y
wuih auh .
E exp λ
E exp (λ ha, wi i) =
u=1

h=1

For any fixed u, define a Multinoulli random variable X such that
P{X = auh } = πuzi h , h ∈ [k].
Recall the definition of wuih that wuih = I{Xui=h } − πuzi h . Then,
!
k
X
2
2
E exp λ
wuih auh = Eeλ(X−EX) ≤ e2λ maxh∈[k] auh ,
h=1

where the last inequality is because |X| ≤ maxh∈[k] |auh | is a bounded random variable.
Consequently, we have
!
m
X

2
2
max auh ≤ E exp 2λ2 .
E exp (λ ha, wi i) ≤ E exp 2λ
u=1

h∈[k]

Therefore, the sub-Gaussian parameter of wi is upper bounded by 2.

D

Proofs of Technical Lemmas

To help readers better understand the results, we first give proofs for spherical Gaussians and
then extend it to the general sub-Gaussians. Without loss of generality, we assume σ = 1 in
this section.

D.1

Proofs for spherical Gaussians

We denote by N (0, Σ) the normal random vectors with mean zero and covariance matrix Σ.
For Z ∼ N (0, 1), we have the following two facts on its moment generating functions.
EeλZ1 = e
2
1
EeλZ1 = √
1 − 2λ

λ2
2

(78)

for all λ < 1/2.

Proof of Lemma 7.1. For any fixed S ⊆ [n], WS =
2

.

P

i∈S

(79)

wi ∼ N (0, |S|Id ). By (79), we have

EeλkWS k = (1 − 2λ|S|)−d/2
for λ < 1/(2|S|). Then Chernoff’s bound yields



d
2
P kWS k ≥ t ≤ exp −λt − log(1 − 2λ|S|) .
2
43

(80)

Choosing λ = 0.49/|S| and t = 1.62(n + 4d)|S|,



d
P kWS k ≥ 1.62(n + 4d)|S| ≤ exp −0.7938n − 3.17d + log 50 .
2
2



Since there are 2n subsets of [n], an union bound argument gives us




p
P max kWS k − 1.62(n + 4d)|S| ≥ 0 ≤ 2n exp (−0.7938n) ≤ exp(−0.1n).
S

i.i.d

Proof of Lemma 7.2. Since wi′ θ ∗ ∼ N (0, kθ ∗ k2 ), we have w̄′ θ ∗ ∼ N (0, kθ ∗ k2 /n). Consequently,




kθ ∗ k2
kθ ∗ k2
′ ∗
√
P −w̄ θ ≥
.
≤ exp −
2
n
Choosing λ = 0.4/n in (80), we have



d
2
P kw̄k ≥ t ≤ exp −0.4nt + log 5 .
2

1
A choice of t = n 3dσ 2 + kθ ∗ k2 yields the desired result.
Proof of Lemma 7.3. Let A = [w1 , w2 , · · · , wn ]. Then
sup

n
X

kak=1 i=1

(a′ wi )2 = sup kAak22 .
kak=1

Let B1 be the unit ball in Rd and C be a ǫ-net of B1 such that for any a ∈ B1 , there is a b ∈ C
satisfying ka − bk ≤ ǫ. Then we have |C| ≤ (1 + 2/ǫ)d [49, Lemma 4.1] and
kAak ≤ kAbk + kA(a − b)k

≤ max kAbk + ka − bk sup kAak
b∈C

a∈B1

≤ max kAbk + ǫ sup kAak.
b∈C

a∈B1

Taking a supreme over a ∈ B1 on both sides and rearranging, we get
1
max kAbk.
sup kAak ≤
1 − ǫ b∈C
a∈B1
i.i.d

For any fixed b ∈ B1 , b′ wi ∼ N (0, 1). Then by (79),
)
( n
n
Y
X
′
2
′
2
−0.25t
Ee0.25(b wi ) ≤ exp (−0.25t + 0.35n) ,
(b wi ) ≥ t ≤ e
P
i=1

i=1

where we have set λ = 0.25. Then (81) and union bound give us


X 
2
P sup kAak ≥ t
≤
P kAbk2 ≥ (1 − ǫ)2 t
a∈B1

b∈C


≤ exp −0.25(1 − ǫ)2 t + 0.35n + d log(1 + 2/ǫ) .

A choice of ǫ = 0.05 and t = 2(n + 9d) yields the desired result.
44

(81)

Proof of Lemma 7.4. For each i ∈ [n], define Gi = S ∩ {i}c . Then,
*
+ *
+
X
X
wi ,
wl ≤ wi ,
wl + kwi k2
l∈S

(82)

l∈Gi

P

WGi = l∈Gi wl is an isotopic Gaussian random vector with variance smaller than |S| on
every direction. By (78) and (79), for any independent Z1 ∼ N (0, 1) and Z2 ∼ N (0, 1),
EeλZ1 Z2 = Ee

2
λ2 Z1
2

=√

1
1 − λ2

for all λ ∈ [0, 1). Since wi and WGi are independent, Chernoff bound implies


o
n
p

d
2
P hwi , WGi i ≥ |S|t ≤ exp −λt − log(1 − λ ) ≤ exp −λt + λ2 d .
2
t 1
, 2 } yields
for λ ≤ 1/2. Choosing λ = min{ 2d


 2 
n
o
p
t t
P hwi , WGi i ≥ |S|t ≤ exp − min
,
.
4d 4

(83)

From (80), a choice of λ = 1/3 gives us
P



kwi k22




t
≥ t ≤ exp − + 0.55d .
3

Setting t = 3(d + 4 log n), we obtain

P kwi k22 ≥ 3σ 2 (d + log(1/δ)) ≤ δ

(84)

with probability greater than 1 − δ. Combining (D.1), (83) and (84) completes the proof.

D.2

Proofs of sub-Gaussians

Proof of Lemma A.1. For any fixed S ⊆ [n], WS is a d dimensional random vector satisfying
Eeha,WS i =

Y

i∈S

Eeha,wi i ≤ e

|S|kak2
2

for all a ∈ Rd . By [29, Theorem 2.1], we have
n
o
√
P kWS k2 ≤ |S|(d + 2 dt + 2t) ≤ exp(−t).

(85)

Note that there are at most 2n possible choices of S. Using the union bound completes the
proof.

45

Proof of Lemma A.2. For any fixed b ∈ Rd with kbk = 1, b′ wi′ s are independent sub-gaussian
random variables with parameter 1. Using [29, Theorem 2.1] again,
)
( n
X
√
(86)
(b′ wi )2 ≤ n + 2 nt + 2t ≤ exp(−t).
P
i=1

Following the same argument as in the proof of Lemma 7.3 and choosing t = 2d + n, we
obtain the desired result.
Proof of Lemma A.3. Similarly to the proof of Lemma 7.4, we upper bound hwi , WS i as
hwi , WS i ≤ hwi , WGi i + kwi k2 .
Since wi and WGi are independent sub-Gaussian random vectors, we have
λhwi ,WGi i

Ee

≤ Ee

λ2 kWG k2
i
2

≤ Ee

λ2 |S|kzk2
2

for λ2 ≤ 1/|S|, where z ∼ N (0, Id ). Here the last inequality is due to [29, Remark 2.3]. Then
83) also holds for sub-Gaussian random vectors. Applying [29, Theorem 2.1] on wi , we have
o
n
p
P kwi k2 ≥ d + 2 d log(1/δ) + 2 log(1/δ) ≤ δ
for δ > 0. Following the same arguments as in the proof of Lemma 7.4, the proof is complete.

Proof of Lemma A.4. For any h ∈ [k], (85) implies
n
o
p
P kWTh∗ k2 ≤ |Th∗ |(d + 4 d log n + 8 log n) ≤ n−4 .

Since n ≥ k, a union bound argument completes the proof.

Proof of Lemma A.5. Let ui = I akθh − θg k2 ≤ hwi , kθh − θg ki . Then ui , i ∈ Tg∗ are independent Bernoulli random variables. Hoeffding’s inequality implies
q
X
X
ui ≤
Eui + 2 n∗g log n
i∈Tg∗

i∈Tg∗

with probability greater than 1 − n−4 . Using Chernoff bound and the sub-Gaussian property
of wi , we have
Eui ≤ E exp −λakθh − θg k2 + λ hwi , kθh − θg ki

≤ exp −λakθh − θg k2 + λ2 kθh − θg k2 /2

= exp −a2 kθh − θg k2 /2 ,



where we set λ = a in the last equality. By union bound and the fact that n ≥ k, we get the
desired result.
46

References
[1] Dimitris Achlioptas and Frank McSherry. On spectral learning of mixtures of distributions. In International Conference on Computational Learning Theory, pages 458–469.
Springer, 2005.
[2] Lada A Adamic and Natalie Glance. The political blogosphere and the 2004 us election:
divided they blog. In Proceedings of the 3rd international workshop on Link discovery,
pages 36–43. ACM, 2005.
[3] Pankaj K Agarwal and Nabil H Mustafa. k-means projective clustering. In Proceedings of
the twenty-third ACM SIGMOD-SIGACT-SIGART symposium on Principles of database
systems, pages 155–165. ACM, 2004.
[4] Ankit Aggarwal, Amit Deshpande, and Ravi Kannan. Adaptive sampling for k-means
clustering. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 15–28. Springer, 2009.
[5] Animashree Anandkumar, Daniel Hsu, and Sham M Kakade. A method of moments for
mixture models and hidden markov models. In COLT, volume 1, page 4, 2012.
[6] David Arthur and Sergei Vassilvitskii. How slow is the k-means method? In Proceedings
of the twenty-second annual symposium on Computational geometry, pages 144–153.
ACM, 2006.
[7] David Arthur and Sergei Vassilvitskii. k-means++: The advantages of careful seeding.
In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms,
pages 1027–1035. Society for Industrial and Applied Mathematics, 2007.
[8] Pranjal Awasthi and Or Sheffet. Improved spectral-norm bounds for clustering. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 37–49. Springer, 2012.
[9] Sivaraman Balakrishnan, Martin J Wainwright, and Bin Yu. Statistical guarantees
for the em algorithm: From population to sample-based analysis. arXiv preprint
arXiv:1408.2156, 2014.
[10] T Tony Cai, Zongming Ma, Yihong Wu, et al. Sparse pca: Optimal rates and adaptive
estimation. The Annals of Statistics, 41(6):3074–3110, 2013.
[11] Arun Tejasvi Chaganty and Percy Liang. Spectral experts for estimating mixtures of
linear regressions. In ICML (3), pages 1040–1048, 2013.
[12] Joseph T Chang. Full reconstruction of markov models on evolutionary trees: identifiability and consistency. Mathematical biosciences, 137(1):51–73, 1996.

47

[13] Peter Chin, Anup Rao, and Van Vu. Stochastic block model and community detection
in the sparse graphs: A spectral algorithm with optimal rate of recovery. arXiv preprint
arXiv:1501.05021, 2(4), 2015.
[14] Sanjoy Dasgupta. The hardness of k-means clustering. Department of Computer Science
and Engineering, University of California, San Diego, 2008.
[15] Sanjoy Dasgupta and Leonard Schulman. A probabilistic analysis of em for mixtures of
separated, spherical gaussians. Journal of Machine Learning Research, 8(Feb):203–226,
2007.
[16] Constantinos Daskalakis, Christos Tzamos, and Manolis Zampetakis. Ten steps of em
suffice for mixtures of two gaussians. arXiv preprint arXiv:1609.00368, 2016.
[17] Alexander Philip Dawid and Allan M Skene. Maximum likelihood estimation of observer
error-rates using the em algorithm. Applied statistics, pages 20–28, 1979.
[18] Neil E Day. Estimating the components of a mixture of normal distributions. Biometrika,
56(3):463–474, 1969.
[19] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from
incomplete data via the em algorithm. Journal of the royal statistical society. Series B
(methodological), pages 1–38, 1977.
[20] Marcel Dettling. Bagboosting for tumor classification with gene expression data. Bioinformatics, 20(18):3583–3593, 2004.
[21] Chao Gao, Yu Lu, and Dengyong Zhou. Exact exponent in optimal rates for crowdsourcing. arXiv preprint arXiv:1605.07696, 2016.
[22] Chao Gao, Zongming Ma, Anderson Y Zhang, and Harrison H Zhou. Achieving optimal
misclassification proportion in stochastic block model. arXiv preprint arXiv:1505.03772,
2015.
[23] Chao Gao, Zongming Ma, Anderson Y Zhang, and Harrison H Zhou. Community detection in degree-corrected block models. arXiv preprint arXiv:1607.06993, 2016.
[24] Gavin J Gordon, Roderick V Jensen, Li-Li Hsiao, Steven R Gullans, Joshua E Blumenstock, Sridhar Ramaswamy, William G Richards, David J Sugarbaker, and Raphael
Bueno. Translation of microarray data into clinically relevant cancer diagnostic tests
using gene expression ratios in lung cancer and mesothelioma. Cancer research,
62(17):4963–4967, 2002.
[25] Moritz Hardt and Eric Price. Tight bounds for learning a mixture of two gaussians. In
Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing,
pages 753–760. ACM, 2015.

48

[26] Ralf Herwig, Albert J Poustka, Christine Müller, Christof Bull, Hans Lehrach, and
John O’Brien. Large-scale clustering of cdna-fingerprinting data. Genome research,
9(11):1093–1105, 1999.
[27] Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt. Stochastic blockmodels: First steps. Social networks, 5(2):109–137, 1983.
[28] Daniel Hsu and Sham M Kakade. Learning mixtures of spherical gaussians: moment
methods and spectral decompositions. In Proceedings of the 4th conference on Innovations in Theoretical Computer Science, pages 11–20. ACM, 2013.
[29] Daniel Hsu, Sham M Kakade, Tong Zhang, et al. A tail inequality for quadratic forms
of subgaussian random vectors. Electron. Commun. Probab, 17(52):1–6, 2012.
[30] Chi Jin, Yuchen Zhang, Sivaraman Balakrishnan, Martin J. Wainwright, and Michael
Jordan. Local maxima in the likelihood of gaussian mixture models: Structural results
and algorithmic consequences. arXiv preprint arxiv:1609.00978, 2016.
[31] Jiashun Jin et al. Fast community detection by score. The Annals of Statistics, 43(1):57–
89, 2015.
[32] Jiashun Jin and Wanjie Wang. Important feature pca for high dimensional clustering.
arXiv preprint arXiv:1407.5241, 2014.
[33] Ravindran Kannan, Hadi Salmasian, and Santosh Vempala. The spectral method for
general mixture models. In International Conference on Computational Learning Theory,
pages 444–457. Springer, 2005.
[34] Ravindran Kannan and Santosh Vempala. Spectral algorithms. Found. Trends Theor.
Comput. Sci., pages 157–288, 2009.
[35] Brian Karrer and Mark EJ Newman. Stochastic blockmodels and community structure
in networks. Physical Review E, 83(1):016107, 2011.
[36] Jason M Klusowski and WD Brinda. Statistical guarantees for estimating the centers of
a two-component gaussian mixture by em. arXiv preprint arXiv:1608.02280, 2016.
[37] Amit Kumar and Ravindran Kannan. Clustering with spectral norm and the k-means
algorithm. In Foundations of Computer Science (FOCS), 2010 51st Annual IEEE Symposium on, pages 299–308. IEEE, 2010.
[38] Jing Lei and Alessandro Rinaldo. Consistency of spectral clustering in sparse stochastic
block models. arXiv preprint arxiv:1312.2050, 2013.
[39] Bruce G Lindsay and Prasanta Basak. Multivariate normal mixtures: a fast consistent
method of moments. Journal of the American Statistical Association, 88(422):468–476,
1993.
49

[40] Stuart Lloyd. Least squares quantization in pcm. IEEE transactions on information
theory, 28(2):129–137, 1982.
[41] Meena Mahajan, Prajakta Nimbhorkar, and Kasturi Varadarajan. The planar k-means
problem is np-hard. In International Workshop on Algorithms and Computation, pages
274–285. Springer, 2009.
[42] Frank McSherry. Spectral partitioning of random graphs. In Foundations of Computer
Science, 2001. Proceedings. 42nd IEEE Symposium on, pages 529–537. IEEE, 2001.
[43] Glenn W Milligan. An examination of the effect of six types of error perturbation on
fifteen clustering algorithms. Psychometrika, 45(3):325–342, 1980.
[44] Ankur Moitra and Gregory Valiant. Settling the polynomial learnability of mixtures
of gaussians. In Foundations of Computer Science (FOCS), 2010 51st Annual IEEE
Symposium on, pages 93–102. IEEE, 2010.
[45] I Ordovás-Pascual and J Sánchez Almeida. A fast version of the k-means classification
algorithm for astronomical applications. Astronomy & Astrophysics, 565:A53, 2014.
[46] Rafail Ostrovsky, Yuval Rabani, Leonard J Schulman, and Chaitanya Swamy. The
effectiveness of lloyd-type methods for the k-means problem. In 2006 47th Annual IEEE
Symposium on Foundations of Computer Science (FOCS’06), pages 165–176. IEEE,
2006.
[47] Karl Pearson. Contributions to the mathematical theory of evolution. Philosophical
Transactions of the Royal Society of London. A, 185:71–110, 1894.
[48] David Pollard. A central limit theorem for k-means clustering. The Annals of Probability,
pages 919–926, 1982.
[49] David Pollard. Empirical processes: theory and applications. In NSF-CBMS regional
conference series in probability and statistics, pages i–86. JSTOR, 1990.
[50] David Pollard et al. Strong consistency of k-means clustering. The Annals of Statistics,
9(1):135–140, 1981.
[51] Mingjun Song and Sanguthevar Rajasekaran. Fast algorithms for constant approximation k-means clustering. Trans. MLDM, 3(2):67–79, 2010.
[52] D Michael Titterington, Adrian FM Smith, and Udi E Makov. Statistical analysis of
finite mixture distributions. Wiley,, 1985.
[53] Santosh Vempala and Grant Wang. A spectral algorithm for learning mixture models.
Journal of Computer and System Sciences, 68(4):841–860, 2004.
[54] Zhaoran Wang, Huanran Lu, and Han Liu. Nonconvex statistical optimization:
Minimax-optimal sparse pca in polynomial time. arXiv preprint arXiv:1408.5352, 2014.
50

[55] Xindong Wu, Vipin Kumar, J Ross Quinlan, Joydeep Ghosh, Qiang Yang, Hiroshi Motoda, Geoffrey J McLachlan, Angus Ng, Bing Liu, S Yu Philip, et al. Top 10 algorithms
in data mining. Knowledge and information systems, 14(1):1–37, 2008.
[56] Ji Xu, Daniel Hsu, and Arian Maleki. Global analysis of expectation maximization for
mixtures of two gaussians. arXiv preprint arXiv:1608.07630, 2016.
[57] Lei Xu and Michael I Jordan. On convergence properties of the em algorithm for gaussian
mixtures. Neural computation, 8(1):129–151, 1996.
[58] Xinyang Yi, Constantine Caramanis, and Sujay Sanghavi. Alternating minimization for
mixed linear regression. In ICML, pages 613–621, 2014.
[59] Anderson Y Zhang and Harrison H Zhou. Minimax rates of community detection in
stochastic block models. arXiv preprint arXiv:1507.05313, 2015.
[60] Yuchen Zhang, Xi Chen, Denny Zhou, and Michael I Jordan. Spectral methods meet
em: A provably optimal algorithm for crowdsourcing. In Advances in neural information
processing systems, pages 1260–1268, 2014.

51

