arXiv:1803.06510v1 [stat.ML] 17 Mar 2018

Hidden Integrality of SDP Relaxation for Sub-Gaussian
Mixture Models
Yingjie Fei and Yudong Chen
School of Operations Research and Information Engineering
Cornell University
{yf275,yudong.chen}@cornell.edu

Abstract
We consider the problem of estimating the discrete clustering structures under Sub-Gaussian
Mixture Models. Our main results establish a hidden integrality property of a semidefinite
programming (SDP) relaxation for this problem: while the optimal solutions to the SDP are
not integer-valued in general, their estimation errors can be upper bounded in terms of the
error of an idealized integer program. The error of the integer program, and hence that of
the SDP, are further shown to decay exponentially in the signal-to-noise ratio. To the best
of our knowledge, this is the first exponentially decaying error bound for convex relaxations
of mixture models, and our results reveal the “global-to-local” mechanism that drives the
performance of the SDP relaxation.
A corollary of our results shows that in certain regimes the SDP solutions are in fact
integral and exact, improving on existing exact recovery results for convex relaxations. More
generally, our results establish sufficient conditions for the SDP to correctly recover the cluster
memberships of (1−δ) fraction of the points for any δ ∈ (0, 1). As a special case, we show that
under the d-dimensional Stochastic Ball Model, SDP
√ achieves non-trivial (sometimes exact)
recovery when the center separation is as small as 1/d, which complements previous exact
recovery results that require constant separation.

1

Introduction

We consider the Sub-Gaussian Mixture Models (SGMMs), where one is given n random points
drawn from a mixture of k sub-Gaussian distributions with different means/centers. SGMMs,
particularly its special case Gaussian Mixture Models (GMMs), are widely used in a broad range of
applications including speaker identification, background modeling and online recommendations
systems. In these applications, one is typically interested in two types of inference problems under
SGMMs:
• Clustering: (approximately) identify the cluster membership of each point, that is, which
of the k mixture components generates a given point;
• Center estimation: estimate the k centers of a mixture, that is, the means of the k components.
Standard approaches to these problems, such as k-means clustering, typically lead to integer programming problems that are non-convex and NP-hard to optimize [Aloise et al., 2009, Jain et al.,
2002, Mahajan et al., 2009]. Consequently, much work has been done in developing computationally tractable algorithms for SGMMs, including expectation maximization [Dempster et al.,
1

1977], Lloyd’s algorithm [Lloyd, 1982], spectral methods [Vempala and Wang, 2004], the method
of moments [Pearson, 1936], and many more. Among them, convex relaxations, including those
based on linear programming (LP) and semidefinite programming (SDP), have emerged as an important approach for clustering SGMMs. This approach has several attractive properties: (a) it is
solvable in polynomial time, and does not require a good initial solution to be provided; (b) it has
the flexibility to incorporate different quality metrics and additional constraints; (c) it is not restricted to specific forms of SGMMs (such as Gaussian distributions), and is robust against model
misspecification [Peng and Xia, 2005, Peng and Wei, 2007, Nellore and Ward, 2015]; (d) it can
provide a certificate for optimality.
Theoretical performance guarantees for convex relaxation methods have been studied in a
body of classical and recent work. As will be discussed in the related work section (Section 2),
these existing results often have one of the two forms:
1. How well the (rounded) solution of a relaxation optimizes a particular objective function
(e.g., the k-means or k-medians objective) compared to the original integer program, as
captured by an approximation factor [Charikar et al., 1999, Kanungo et al., 2004, Peng and
Wei, 2007, Li and Svensson, 2016];
2. When the solution of a relaxation corresponds exactly to the ground-truth clustering, a phenomenon known as exact recovery and studied in a more recent line of work [Nellore and
Ward, 2015, Awasthi et al., 2015, Mixon et al., 2017, Iguchi et al., 2017, Li et al., 2017].
In many practical scenarios, optimizing a particular objective function, and designing approximation algorithms for doing so, are often only a means to solving the two inference problems above,
namely learning the true underlying model that generates the observed data. Results on exact recovery are more directly relevant to this goal; however, such results often require very stringent
conditions on the separation or signal-to-noise ratio (SNR) of the model. In practice, convex relaxation solutions are rarely exact, even when the data are generated from the assumed model. On
the other hand, researchers have observed that the solutions, while not exact or integer-valued, are
often a good approximation to the desired solution that represents the ground truth [Mixon et al.,
2017]. Such a phenomenon is not captured by the results on exact recovery.
In this paper, we aim to strengthen our understanding of the convex relaxation approach to
SGMMs. In particular, we study the regime where solutions of convex relaxations are not integral
in general, and seek to directly characterize the estimation errors of the solutions—namely, their
distance to desired integer solution corresponding to the true underlying model.

1.1

Our contributions

For a class of SDP relaxations for SGMMs, our results reveal a perhaps surprising property of
them: while the SDP solutions are not integer-valued in general, their errors can be controlled by
that of the solutions of an idealized integer program (IP), in which one tries to estimate cluster
memberships when an oracle reveals the true centers of the SGMM. We refer to the latter program
as the Oracle Integer Program. In particular, we show that, in a precise sense to be formalized
later, the estimation errors of the SDP and Oracle IP satisfy the relationship (Theorem 1):
error(SDP) . error(IP)
under certain conditions. We refer to this property as hidden integrality of the SDP relaxation; its
proof in fact involves showing that the optimal solutions of certain intermediate linear optimization

2

problems are integral. We then further upper bound the error of the Oracle IP and show that it
decays exponentially in terms of the SNR (Theorem 2):
h

i
error(IP) . exp −Ω SNR2 ,
where the SNR is defined as the ratio of the center separation and the standard deviation of the
sub-Gaussian components. Combining these two results immediately leads to explicit bounds on
the error of the SDP solutions (Corollary 1).
1.1.1

Consequences

When the SNR is sufficiently large, the above results imply that the SDP solutions are integral and
exact up to numerical errors, hence recovering (sometimes improving) existing results on exact
recovery as a special case. Moreover, if the SNR is low and the SDP solutions are fractional, one
may obtain an explicit clustering from the SDP solutions via a simple, optimization-free rounding
procedure. We show that the error of this explicit clustering (in terms of the fraction of points
misclassified) is also bounded by the error of the Oracle IP and hence also decays exponentially
in the SNR (Theorem 3). As a consequence, we obtain sufficient conditions for misclassifying at
most δ fraction of the points for any given δ ∈ [0, 1]. Finally, we show that the SDP solutions
also lead to an efficient estimator of the cluster centers, for which estimation error bounds are
established (Theorem 4).
Significantly, our results often match and sometimes improve upon state-of-the-art performance guarantees in settings for which known results exist, and lead to new guarantees in other
less studied settings of SGMMs. For instance, a corollary of our results shows that under the
Stochastic Ball Model, SDP achieves
√ meaningful (sometimes exact) recovery even when the center separation ∆ is as small as O( 1/d), where d is the dimension (Section 4.5). For high dimensional settings, this bound generalizes existing results that focus on exact recovery and require
∆ = Ω(1). Detailed discussions of the implications of our results and comparison with existing
ones will be provided after we state our main theorems.
1.1.2

The “global-to-local” phenomenon

Our results above are obtained in two steps: (a) relating the SDP to the Oracle IP, and (b) bounding
the Oracle IP errors. Conceptually, this two-step approach allows us to decouple two types of
mechanisms that determine the performance of the SDP relaxation approach.
• On the one hand, step (a) is done by leveraging the structures of the entire dataset with n
points. In particular, certain global spectral properties of the data ensure that the error of
the SDP is non-trivial and bounded (in terms of the Oracle IP error). This step is relatively
insensitive to the specific structure of the SGMM.
• On the other hand, as shall become clear in the sequel, the Oracle IP essentially reduces to n
independent clustering problems, one for each data point. Knowing the true cluster centers,
the Oracle IP is optimal in terms of the clustering errors. No other algorithms (including
SDP relaxations) can achieve a strictly better error, due to the inherent randomness of individual data points. Step (b) above hence captures the local mechanism that determines
fine-grained error rates as a function of the SNR.
Our two-step analysis establishes the hidden integrality property as the interface between these two
types of mechanisms. As a clustering algorithm, the SDP approach is powerful enough to capture
3

these two mechanisms simultaneously, without requiring a good initial solution or sophisticated
pre-processing/post-processing steps.
We note that recent work on clustering under Stochastic Block Models reveals a related “localto-global amplification” phenomenon, which connects the original clustering problem with that of
recovering a single point when the memberships of the others are known [Abbe and Sandon, 2015,
Abbe et al., 2016, Abbe, 2017]. Our results share similar spirits with this line of work, though our
models, algorithms and proof techniques are quite different.

1.2

Paper Organization

The remainder of the paper is organized as follows. In Section 2, we discuss related work on
SGMMs and its special cases. In Section 3, we describe the problem setup for SGMMs and
provide a summary of our clustering algorithms. In Section 4, we present our main results, discuss
some of their consequences and compare them with existing results. The paper is concluded with
a discussion of future directions in Section 5. The proofs of our main theorems are deferred to the
Appendix.

2

Related work

The study of SGMMs has a long history and is still an active area of research. Here we review the
most relevant results with theoretical guarantees.
Dasgupta [1999] is among the first to obtain performance guarantees for GMMs. Subsequent
work has obtained improved guarantees, achieved by various algorithms including spectral methods. These results often establish sufficient conditions, in terms of the separation between the
cluster centers (or equivalently the SNR), for achieving (near)-exact recovery of the cluster memberships. Vempala and Wang [2004] obtain one of the best results and require SNR & (k ln n)1/4 ,
which is later generalized and extended in a long line of work including Achlioptas and McSherry
[2005], Kumar and Kannan [2010], Awasthi and Sheffet [2012]. We compare these results with
ours in Section 4.
Expectation-Maximization (EM) and Lloyd’s algorithm are among the most popular methods
for GMMs. Despite their empirical effectiveness, non-asymptotic statistical guarantees are established only recently. In particular, convergence and center estimation error bounds for EM under
GMMs with two components are derived in Balakrishnan et al. [2017], Klusowski and Brinda
[2016], with extension to multiple components given in Yan et al. [2017]. The work of Lu and
Zhou [2016] provides a general convergence analysis for Lloyd’s algorithm, which implies clustering and center estimation guarantees for random models including SGMMs. All these results
assume that one has access to a sufficiently good initial solution, typically obtained by spectral
methods. A recent breakthrough was made by Daskalakis et al. [2016], Xu et al. [2016], who
establish global convergence of randomly-initialized EM for GMMs with two symmetric components. Complementarily, Jin et al. [2016] show that EM may fail to converge under GMMs with
k ≥ 3 components due to the existence of bad local minima. Robustness of the Lloyd’s algorithm
under a semi-random GMM is studied in Awasthi and Vijayaraghavan [2017].
Most relevant to us are work on convex relaxation methods for GMMs and k-means/median
problems. A class of SDP relaxations are developed in the seminal work by Peng and Xia [2005],
Peng and Wei [2007]. Thanks to convexity, these methods do not suffer from the issues of bad
local minima faced by EM and Lloyd’s, though it is far from trivial to round their (typically fractional) solutions into valid clustering solutions with provable quality guarantees. In this direction,

4

Awasthi et al. [2015], Iguchi et al. [2017], Li et al. [2017] establish conditions for LP/SDP relaxations to achieve exact recovery. The work of Mixon et al. [2017] considers SDP relaxations as
a denoising method, and proves error bounds for a form of approximate recovery. Most of these
results are directly comparable to ours, and we discuss them in more details in Section 4 after
presenting our main theorems.
Clustering problems under Stochastic Block Models (SBMs) have also witnessed fruitful
progress on understanding convex relaxation methods; see Abbe [2017] for a survey and further references. Much work has been done on exact recovery guarantees for SDP relaxations of
SBMs [Krivelevich and Vilenchik, 2006, Oymak and Hassibi, 2011, Ames and Vavasis, 2014,
Chen et al., 2014, Amini and Levina, 2018]. A more recent line of work establishes approximate
recovery guarantees of the SDPs [Guédon and Vershynin, 2016, Montanari and Sen, 2016] in the
low SNR regime. Particularly relevant to us is the work by Fei and Chen [2017], who also establish exponentially decaying error bounds. Despite the apparent similarity in the forms of the error
bounds, our results require very different analytical techniques, due to the fundamental difference
between the geometric and probabilistic structures of SBMs and SGMMs; moreover, our results
reveal the more subtle hidden integrality property of SDP relaxations, which we believe holds
more broadly beyond specific models like SBMs and SGMMs.

3

Models and algorithms

In this section, we formally set up the clustering problem under SGMMs and describe our SDP
relaxation approach.

3.1

Notations

We first introduce some notations. Vectors and matrices are denoted by bold letters such as u and
M. For a vector u, we denote by ui its i-th entry. For a matrix M, Tr(M) denotes its trace, Mi j its
P
(i, j)-th entry, diag (M) the vector of its diagonal entries, kMk1 B i, j Mi j its entry-wise `1 norm,
Mi• its i-th row and M• j its j-th column. We write M  0 if M is symmetric positive semidefinite
(psd). The trace inner product between two matrices M and Q of the same dimension is denoted
by hM, Qi B Tr(M> Q). For a number a, M ≥ a means Mi j ≥ a, ∀i, j. We denote by 1m the
all-one column vector of dimension m. For a positive integer i, let [i] B {1, 2, . . . , i}. For two
non-negative sequences {ai } and {bi }, we write ai . bi if there exists a universal constant C > 0
such that ai ≤ Cbi for all i, and write ai  bi if ai . bi and bi . ai .
We recall that the sub-Gaussian norm of a random variable X is defined as
n
o


kXkψ2 B inf t > 0 : E exp X 2 /t2 ≤ 2 ,
and X is called sub-Gaussian if kXkψ2 < ∞. Note that Normal and bounded random variables are
sub-Gaussian. Denote by Sm−1 the m-dimensional unit `2 sphere. A random vector x ∈ Rm is
sub-Gaussian if the one dimensional marginals hx, ui are sub-Gaussian random variables for all
u ∈ Rm . The sub-Gaussian norm of x is defined as kxkψ2 B supu∈Sm−1 k hx, ui kψ2 .

3.2

Sub-Gaussian Mixture Models

We focus on Sub-Gaussian Mixture Models (SGMMs) with balanced clusters.

5

Model 1 (Sub-Gaussian Mixture Models). Let µ1 , . . . , µk ∈ Rd be k unknown cluster centers. We
observe n random points in Rd of the form
hi B µσ∗ (i) + gi ,

i ∈ [n]

where σ∗ (i) ∈ [k] is the unknown cluster label of the i-th point, and {gi } are independent subGaussian random vectors with sub-Gaussian norms kgi kψ2 ≤ τ.1 We assume that the ground-truth
clusters have equal sizes, that is, |{i ∈ [n] : σ∗ (i) = a}| = nk for each a ∈ [k].
Note that we do not require {gi } to be identically distributed or isotropic. Model 1 includes
several important mixture models as special cases:
• The spherical GMM, where {gi } are Gaussian with the covariance matrix τ2 I.
• More general GMMs with non-identical and non-diagonal covariance matrices {Σa }a∈[k] for
the k clusters.
• The Stochastic Ball Model [Nellore and Ward, 2015], where the distributions of {gi } are
supported on the unit `2 ball in Rd ; we discuss this model in details in Section 4.5
Throughout the paper we assume n ≥ 4 and k ≥ 2 to avoid degeneracy. Let σ∗ ∈ [k]n be the vector
of the true cluster labels; that is, its i-th coordinate is σ∗i ≡ σ∗ (i) (we use them interchangeably
throughout the paper.) The task is to estimate the underlying clustering σ∗ given the observed data
{hi : i ∈ [n]}. The separation of the centers of clusters a and b is denoted by ∆ab B kµa − µb k2 ,
and ∆ B mina,b∈[k] kµa − µb k2 is the minimum separation of the centers. Playing a crucial role in
our results is the quantity
∆
sB ,
(1)
τ
which is a measure of the SNR of the SGMM.

3.3

Semidefinite programming relaxation

We now describe our SDP relaxation for clustering SGMMs. To begin, note that each candidate
clustering of n points into k clusters can be represented using an assignment matrix F ∈ {0, 1}n×k
where



1 if point i is assigned to cluster a
Fia = 

0 otherwuse.
n
o
Let F B F ∈ {0, 1}n×k : F1k = 1n be the set of all possible assignment matrices. Given the
points {hi } to be clustered, a natural approach is to find an assignment F that minimizes the total
within-cluster pairwise distance. Arranging the pairwise squared distance as a matrix A ∈ Rn×n
with
Ai j = khi − h j k22 ,
for each (i, j) ∈ [n] × [n],
we can express the above objective as
X
X
D
E
Ai j 1{points i and j are assigned to the same cluster} =
Ai j (FF> )i j = FF> , A .
i, j

i, j

Therefore, the approach described above is equivalent to solving the integer program (2) below:
1

More explicitly, the sub-Gaussian assumption is equivalent to
Eehgi ,vi ≤ e

τ2 kvk2
2
2

,

6

∀v ∈ Rd .

D
E
min FF> , A

min hY, Ai
Y

F

s.t. F ∈ F

n
s.t. Y1n = 1n
k
Y is psd

(2)

n >
1>
n F = 1k
k

(3)

diag (M) = 1n
Y ∈ {0, 1}n×n ; rank(Y) = k.

n >
In program (2) the additional constraint 1>
n F = k 1k enforces that all k clusters have the same
n
size k , as we are working with an SGMM whose true clusters are balanced. Under this balanced
model, it is not hard to see that the program (2) is equivalent to the classical k-means formulation.
With a change of variable Y = FF> , we may lift the program (2) to the space of n × n matrices and
obtain the equivalent formulation (3). Both programs (2) and (3) involve non-convex combinatorial constraints and are computationally hard to solve. To obtain a tractable formulation, we drop
the non-convex rank constraint in (3) and replace the integer constraint with the box constraint
0 ≤ Y ≤ 1 (the constraint Y ≤ 1 is in fact redundant). These lead to the following SDP relaxation:

b ∈ arg min hY, Ai
Y
Y∈Rn×n

n
s.t. Y1n = 1n
k
Y is psd

(4)

diag (M) = 1n
Y ≥ 0.
The performance of the SDP is measured against the true cluster matrix Y∗ ∈ {0, 1}n×n , where for
each i, j ∈ [n],



1 if σ∗ (i) = σ∗ ( j), i.e., points i and j are in the same cluster,
Yi∗j = 

0 if σ∗ (i) = σ∗ ( j), i.e., points i and j are in different clusters,
with the convention Yii∗ = 1, ∀i ∈ [n]. The true cluster matrix Y∗ encodes the ground-truth clusb to (4) as an estimate of
tering σ∗ , and is feasible to program (4) . We view an optimal solution Y
∗
b − Y∗ k1
the true clustering Y . Our goal is to characterize the cluster recovery/estimation error kY
in terms of the number of points n, the number of clusters k, the data dimension d and the SNR
b in `1 metric; as we shall see later, this
s defined in (1). Note that here we measure the error of Y
metric is directly related to the clustering error (i.e., the fraction of misclassified points).
We remark that the SDP (4) is somewhat different from the more classical and well-known
SDP relaxation of k-means proposed by Peng and Wei [2007]. The SDP (4) is closely related to
the one considered by Amini and Levina [2018] in the context of Stochastic Block Models, though
it seems to be less studied under SGMMs, with the notable exception of Li et al. [2017].

3.4

Explicit clustering

b which is not integral
Our main results in the next section directly concern the SDP solution Y,
in general and hence does not directly correspond to an explicit clustering. In case an explicit
b using a
clustering is desired, we may easily extract cluster memberships from the solution Y
simple procedure.

7

The procedure consists of two steps given as Algorithms 1 and 2, respectively. In the first step,
b as points in Rn , and consider the `1 balls centered at each row with a certain
we treat the rows of Y
radius. The ball that contains the most rows is identified, and the indices of the rows in this ball
b This
are output and removed. The process continues iteratively with the remaining rows of Y.
n
step outputs a number of sets whose sizes are no larger than k but may not equal to each other.
Algorithm 1 First step
b ∈ Rn×n , number of points n, number of clusters to extract k.
Input: data matrix Y
1. B0 ← ∅, t ← 0, V ← [n]
S
2. While V\ ti=0 Bi , ∅:
(a) t ← t + 1
S
(b) Vt ← V\ t−1
i=0 Bi
n
bu• − Y
bw• k1 ≤
(c) For each u ∈ Vt : B(u) ← w ∈ Vt : kY

n
4k

o
.

(d) Bt ← arg maxB(u):u∈Vt |B(u)|
(e) If |Bt | > nk , then remove arbitrary elements in Bt so that |Bt | =

n
k

Output: sets {Bt }t≥1 .
In the second step, we convert the sets output by Algorithm 1 into k equal-size clusters. This
is done by identifying the k largest sets among them, and distributing points in the remaining sets
across the chosen k sets so that each of these sets contains exactly nk points.
Algorithm 2 Second step
Input: approximate clustering sets {Bt }t≥1 , number of points n, number of clusters to extract k.
1. k0 ← {Bt }t≥1
2. Choose k largest sets among {Bt }t≥1 and rename the chosen sets as {Ut }t∈[k]
3. Arbitrarily distribute elements of {Bt }t≥1 \ {Ut }t∈[k] among {Ut }t∈[k] so that each Ut has exactly nk elements
4. For each i ∈ [n]: b
σi ← t, where t is the unique index in [k] such that i ∈ Ut
b ∈ [k]n .
Output: clustering assignment vector σ
Combining the above two algorithms gives our final algorithm, cluster, for extracting an
b This procedure is given as Algorithm 3.
explicit clustering from the SDP solution Y.
Algorithm 3 cluster
b ∈ Rn×n , number of points n, number of clusters to extract k.
Input: data matrix Y
b n and k as input and get {Bt }t≥1
1. Run Algorithm 1 with Y,
b
2. Run Algorithm 2 with {Bt }t≥1 , n and k as input and get σ
b ∈ [k]n .
Output: clustering assignment σ
8

The output of the above procedure
b n, k)
b B cluster(Y,
σ
is a vector in [k]n such that point i is assigned to the b
σi -th cluster. We are interested in controlling
b relative to the ground-truth clustering σ∗ . Let Sk denote the symmetric
the clustering error of σ
group consisting of all permutations of [k]. The clustering error is defined by
err(b
σ, σ∗ ) B min
π∈Sk

o
1 n
i ∈ [n] : b
σi , π(σ∗i ) ,
n

(5)

which is the proportion of points that are misclassified, modulo permutations of the cluster labels.
Variants of the above cluster procedure have been considered before by Makarychev et al.
[2016], Mixon et al. [2017]. Our results in the next section establish that the clustering error
b − Y∗ k1 of the SDP solution Y.
b
err(b
σ, σ∗ ) is always upper bounded by the `1 error kY

4

Main results

In this section, we establish the connection between the estimation error of the SDP relaxation (4)
and that of what we call the Oracle Integer Program. Using this connection, we derive explicit
bounds on the error of the SDP, and explore their implications for clustering and center estimation.

4.1

Oracle Integer Program


Consider an idealized setting where an oracle reveals the true cluster centers µa a∈[k] . Moreover,
n o
we are given the data points h̄i
, where h̄i B µσ∗ (i) + (2c)−1 gi for c = 18 and {gi } are the same
i∈[n]
n o
realizations of the random variables in the original SGMM. In other words, h̄i are the same as the
original data points {hi }, except that the standard deviation (or more generally, the sub-Gaussian
−1
norm) of the noise
n o {gi } is scaled by (2c) = 4.
To cluster h̄i in this idealized setting, a natural approach is to simply assign each point to the
closest cluster center, so that the total distance of the points to their assigned centers is minimized.
We may formulate this procedure as an integer program, by representing each candidate clustering
assignment using an assignment matrix F ∈ F as before. Then, for each assignment matrix F, the
quantity
XX
η(F) B
kh̄ j − µa k22 F ja
j

a

is exactly the sum of the distances of each point to its assigned cluster center. The clustering
procedure above thus amounts to solving the following Oracle Integer Program (IP):
min η(F),

s.t. F ∈ F .

F

(6)

As this program is separable across the rows of F, it can be reduced to n independent optimization
problems, one for each data point h̄i .
Let F∗ ∈ F be the assignment matrix associated with the true underlying clustering of the
SGMM; that is, F ∗ja = 1 {σ∗ ( j) = a} for each j ∈ [n], a ∈ [k]. For each feasible solution F ∈ F to
the Oracle IP, it is easy to see that the quantity 12 kF − F∗ k1 is exactly the number of points that are
assigned differently in F and F∗ , and hence measures the clustering error of F with respect to the
ground truth F∗ .
9

A priori, there is no obvious connection between the estimation error of a solution F to the
above Oracle IP and that of a solution to the SDP. In particular, the latter involves a continuous
relaxation, for which the solutions are fractional in general and the true centers are unknown.
Surprisingly, we are able to establish a formal relationship between the two, and in particular
show that the error of the SDP is bounded by the error of the IP in an appropriate sense.

4.2

Errors of SDP relaxation and Oracle IP

To establish the connection between the SDP and Oracle IP, we begin with the following observation: for a solution F ∈ F to potentially be an optimal solution of the Oracle IP (6), it must satisfy
η(F) ≤ η(F∗ ) since F∗ is feasible to (6). Consequently, the quantity
(
)
1
∗
∗
max kF − F k1 : F ∈ F , η(F) ≤ η(F )
(7)
2
represents the worst-case error of a potentially optimal solution to the Oracle IP. This quantity
b to the SDP relaxation, as is
turns out to be an upper bound of the error of the optimal solution Y
shown in the theorem below.
Theorem 1 (IP bounds SDP). Under Model 1, there exist some universal constants C s > 0, C ≥ 1
for which the following holds. If the SNR satisfies
r

 kd

kd
2
s ≥ C s 
log n +
+ k ,
(8)
n
n
then we have
(
)
b − Y∗ k1
kY
kF − F∗ k1
∗
≤ 2 · max
: η(F) ≤ η(F ), F ∈ F
kY∗ k1
kF∗ k1
with probability at least 1 − n−C − 2e−n .
The proof is given in Section B, and consists of two main steps: (i) showing that with high
probability the SDP error is upper bounded by the objective value of a linear program (LP), and
(ii) showing that the LP admits an integral optimal solution and relating this solution to the quantity (7). We note that the key step (ii), which involves establishing certain hidden integrality
properties, is completely deterministic. The SNR condition (8) is required only in the probabilistic step (i). As we elaborate in Sections 4.3–4.5, our SNR condition holds even in the regime
where exact recovery is impossible, and is often milder than existing results on convex relaxations.
Sharper analysis in step (i) above will lead to potentially more relaxed conditions on the SNR.
To obtain an explicit bound on the SDP error, it suffices to upper bound the error of the Oracle
IP. This turns out to be a relatively easy task compared to directly controlling the error of the
SDP. The reason is that the Oracle IP has only finitely many feasible solutions, allowing one to
use a union-bound-like argument. Our analysis establishes that the error of the Oracle IP decays
exponentially in the SNR, as summarized in the theorem below.
Theorem 2 (Exponential rate of IP). Under Model 1, there exist universal constants C s , Cg , Ce > 0
for which the following holds. If s2 ≥ C s k, then we have
(
)
" 2#
kF − F∗ k1
s
∗
max
: η(F) ≤ η(F ), F ∈ F ≤ Cg exp −
∗
kF k1
Ce
with probability at least 1 − 32 n−1 .
10

The proof is given in Section C. An immediate consequence of Theorems 1 and 2 is that the
SDP (4) also achieves an exponentially decaying error rate.
Corollary 1 (Exponential rate of SDP). Under Model 1 and the SNR condition (8), there exist
universal constants Cm , Ce > 0 such that
" 2#
b − Y∗ k1
s
kY
≤ Cm exp −
kY∗ k1
Ce
with probability at least 1 − 2n−1 .
b using the cluster prob extracted from Y
Our last theorem concerns the explicit clustering σ
cedure described in Section 3.4. In particular, we show that the number of misclassified points
b and in turn by the error of the Oracle IP; consequently, σ
b
is upper bounded by the error in Y
misclassifies an exponentially-decaying number of points.
b
b is always upper bounded by the error in Y:
Theorem 3 (Clustering error). The error rate in σ
err(b
σ, σ∗ ) .

b − Y∗ k1
kY
.
kY∗ k1

Consequently, under Model 1 and the SNR condition (8), there exist universal constants Cm , Ce > 0
such that
" 2#
s
∗
b, σ ≤ Cm exp −
err σ
Ce
with probability at least 1 − 2n−1 .
The proof is given in Section E. Note that the above bound, in terms of the clustering error, is
optimal (up to a constant in the exponent) in view of the minimax results in Lu and Zhou [2016].

4.3

Consequences

We explore the consequences of our error bounds in Corollary 1 and Theorem 3.
• Exact recovery: If the SNR s satisfies the condition (8) and moreover s2 & log n, then


b, σ∗ < n1 , which means that err σ
b, σ∗ = 0 and the true
Theorem 3 guarantees that err σ
underlying clustering is recovered exactly. Note that these conditions can be simplified to
s2 & k + log n when n & d. In fact, in this case Corollary 1 guarantees the SDP solution
b − Y∗ k1 < 1 , so simply rounding Y
b element-wise produces the groundsatisfies the bound kY
4
∗
truth cluster matrix Y . Therefore, the SDP relaxation is able to achieve exact recovery
(sometimes called strong consistency in the literature on Stochastic Block Models [Abbe,
2017]) of the underlying clusters when the SNR is sufficiently large.
In fact, our results are applicable even in regimes with a lower SNR, for which exact recovery
of the clusters is impossible due to the overlap between points from different clusters. In such
regimes, Corollary 1 and Theorem 3 imply approximate recovery guarantees for the SDP relaxation:
• Almost exact recovery: If s satisfies the condition (8) and s2 = ω (1), then Theorem 3

b, σ∗ = o (1). That is, the SDP recovers the cluster memberships of
implies that err σ
almost all points asymptotically, which is sometimes called weak consistency.
11

• Recovery with δ-error: More generally, for any number δ ∈ (0, 1), Theorem 3 implies the
following non-asymptotic recovery guarantee: If s satisfies the condition (8) and s2 & log 1δ ,

b, σ∗ ≤ δ. That is, the SDP correctly recovers the cluster memberships of at
then err σ
least (1 − δ) fraction of the points.
We compare the above results with existing ones in Section 4.4 to follow.

4.4

Comparison with existing results

Table 1 summarizes several most representative results in the literature on clustering SGMM/GMM.
Most of them are in terms of SNR conditions required to achieve exact recovery of the underlying
clusters. Note that our results imply sufficient conditions for both exact and approximate recovery.
Most relevant to us is the work of Li et al. [2017], which considers similar SDP relaxation
formulations. They show that exact recovery is achieved when s2 & k + log n and n  d2 k3 log k.
In comparison, a special case of our Corollary 1 guarantees exact recovery whenever s2 & k +log n
and n & d, which is milder than the condition in Li et al. [2017].
The work in Lu and Zhou [2016] also proves an exponentially decaying clustering error rate,
but for a different algorithm (Lloyd’s algorithm). To achieve non-trivial approximate recovery of
the clusters, they require s2  k2 + k3 dn and k3  logn n as n → ∞. Our SNR condition in (8) has
milder dependency on k, though dependency on n and d are a bit more subtle. We do note that
under their more restricted SNR condition, Lu and Zhou [2016] are able to obtain tight constants
in the exponent of the error rate.
Finally, the work of Mixon et al. [2017] considers the SDP relaxation introduced by Peng
and Wei [2007] and provides bounds on center estimation when s2 & k2 . An intermediate result
of theirs concerns errors of the SDP solutions; under the setting of balanced clusters, their error
bound can be compared with ours after appropriate rescaling. In particular, their result implies the
b − Y∗ k2 . n22 when n is sufficiently large. This bound is non-trivial when s2 & k
error bound kY
F
s
since kY∗ k2F =
bound

n2
k .

Under the same conditions on s2 and n, our results imply the exponential error
2
b − Y∗ k2 ≤ kY
b − Y∗ k1 . n e−s2 ,
kY
F
k

which is strictly better.
To sum up, corollaries of our results provide more relaxed conditions for exact or approximate
√
recovery compared to most of the existing results listed in Table 1. Our results are weaker by a k
factor than the one in Vempala and Wang [2004], which considers spectral clustering methods and
focuses on exact recovery under spherical Gaussian mixtures. In comparison, our results apply
to the more general sub-Gaussian setting, and imply exponential error bounds for approximate
recovery guarantees under more general SNR conditions.

4.5

Stochastic Ball Model

We illustrate the power of our main results for the Stochastic Ball Model introduced in Nellore
and Ward [2015].
Model 2 (Stochastic Ball Model). Under Model 1, we assume in addition that each gi is sampled
from a rotationally invariant distribution supported on the unit `2 ball in Rd .
Under the Stochastic Ball Model, each data point hi = µσ∗ (i) + gi is sampled from the unit ball
around its cluster center in a rotationally invariant fashion. This model is a special case of SGMM
with its sub-Gaussian norm given below:
12

Paper
Vempala and Wang [2004]
Achlioptas and McSherry
[2005]
Kumar and Kannan [2010]
Awasthi and Sheffet [2012]
Lu and Zhou [2016]
Mixon et al. [2017]
Li et al. [2017]
Ours

Condition on SNR2
p

Ω k log n


Ω k log n + k2


Ω k2 · polylog (n)

Ω k · polylog
  (n)
Ω k2


Ω k2 + log n
 
Ω k2

Ω k + log n
Ω (k)

Ω k + log n

Recovery Type

Algorithm

Exact

Spectral

Exact

Spectral

Exact
Exact
Approximate
Exact
For center estimation
Exact
Approximate
Exact

Spectral
Spectral
Spectral
+Lloyd’s
SDP
SDP
SDP

Table 1: Summary of existing results on cluster recovery for GMM. Here “approximate” means
correct recovery of the memberships of at least (1 − δ) fraction of the points for a fixed constant
δ ∈ (0, 1). Some of the results listed assume that n  poly(k, d); see Section 4.4 for details.
Fact 1. Under Model 2, each gi has sub-Gaussian norm kgi kψ2 ≤ τ = C
constant C > 0.

q

1
d

for some universal

For completeness we prove this claim in Section F. Specializing Corollary 1 and Theorem 3
to the Stochastic Ball Model, we obtain the following sufficient conditions on the minimum center
separation ∆ for various types of recovery:
r
∆2 &

 log n


for exact recovery

d ,

k
k k 
 ω(1)
for almost exact recovery
log n + + + 

d ,
nd
n d 


 log δ−1 , for recovery of (1 − δ) fraction of the notes.
d

The state-of-the-art results for Stochastic Ball Models are given in Awasthi et al. [2015], Iguchi
et al. [2017], Li et al. [2017], which establish that SDP achieves exact recovery when n is sufficiently large and ∆2 ≥ 4+(k, d) for some non-negative function (·). Regardless of the values of k
and d, these results all require the separation to satisfy ∆2 = Ω(1) and thus the balls to be disjoint.
In contrast, our results above are applicable to a small-separation regime that is not covered by
these results. In particular, when n is large and k = O(1), our results guarantee that SDP achieves
approximate recovery when ∆2 & d1 , which can be arbitrarily smaller than Ω(1) when the dimension d grows. Moreover, the recovery is exact if ∆2 & logd n , which can again be arbitrarily small
as long as n does not grow exponentially fast (i.e., n = eo(d) ). Therefore, in the high dimensional
setting, our results guarantee strong performance of the SDP even when the centers are very close
and the balls overlap with each other.
It may appear a bit counter-intuitive that exact/approximate recovery is achievable when the
separation is so small. Such a result is a manifestation of the geometry in high dimensions: the
relative volume of the intersection of two balls vanishes as the dimension grows. As q
a passing
2
note, our exact recovery result above does not contradict the necessary condition ∆ ≥ 1+ 1 + d+2
given in Li et al. [2017, Corollary 4.3], as they allow n to grow arbitrarily fast, in which case with
high probability some points will land in the intersection.

13

4.6

Cluster center estimation

To conclude this section, we show that our main theorems also imply a guarantee for center esb produced by the SDP relaxation, we
timation. In particular, given the estimated cluster labels σ
may obtain an estimate of the cluster centers {µa } by
µ̂a B

k X
hi ,
n

a ∈ [k].

i:b
σi =a

That is, we simply compute the empirical means of the points within each estimated clusters. As a
corollary of our bounds on clustering errors, we obtain the following center estimation guarantee.
Theorem 4 (Cluster center estimation error). Suppose that maxa,b∈[k] ∆ab ≤ Cq ∆ for some universal constant Cq > 0. Under Model 1 and the SNR condition (8), there exist universal constants
Cm , Ce > 0 such that
r
" 2 #

 k d + log n  p
s 
+ d + log n · exp −
max min kµ̂a − µπ(a) k2 ≤ Cm τ 

a∈[k] π∈Sk
n
Ce 
with probability at least 1 − 4n−1 .
The proof is given in Section G. Note that the error is measured again up to permutation
q of the
k(d+log n)
cluster labels. Our error bound in Theorem 4 consists of two terms. The first term, τ
,
n
n
corresponds to the error of estimating a d-dimensional cluster center vector using the k data points
from that cluster. This error is unavoidable even when the true cluster labels are known. On the
other hand, the second term captures the error due to incorrect cluster labels
q for some of the points.
When s2 & log n and d & log n, we achieve the minimax optimal rate τ

5

d
n/k

for center estimation.

Conclusion

In this paper, we consider clustering problems under SGMMs using an SDP relaxation. Our analysis consists of two steps: (a) we establish the main result of this paper that the clustering error
of the SDP can be controlled by that of the idealized Oracle IP provided that an SNR condition
is satisfied; (b) we show that the error of the Oracle IP decays exponentially in the SNR. As immediate corollaries, we obtain sufficient conditions for the SDP to achieve exact and approximate
recovery, as well as error bounds for estimating the mixture centers.
As mentioned, this two-step approach allows for a certain decoupling of the computational
and statistical mechanism that determines the performance of the SDP approach. We expect that
further progress in understanding of the SDP relaxations are likely to come from improvements in
step (a). On the other hand, by modifying and sharpening step (b), one may generalize our results
to other variants of SGMMs.
Our work points to several interesting future directions. An immediate problem is extending
our results to the case of imbalanced clusters. It is also of interest to study the robustness of SDP
relaxations for SGMMs by considering various semi-random models that allow for adversarial
attacks, arbitrary outliers and model misspecification [Awasthi and Vijayaraghavan, 2017]. Other
directions worth exploring include obtaining better constants in error bounds, identifying sharp
thresholds for different types of recovery, and obtaining tight localized proximity conditions in the
lines of Li et al. [2017].
14

Acknowledgement
Y. Fei and Y. Chen were partially supported by the National Science Foundation CRII award
1657420 and grant 1704828.

Appendices
A

Additional notations

b − Y∗ k1 . For a matrix M, we write kMk∞ B maxi, j Mi j as its
We define the shorthand γ B kY
entry-wise `∞ norm, and kMkop as its spectral norm (maximum singular value). We let I and J
be the n × n identity matrix and all-one matrix, respectively. For a real number x, dxe denotes its
ceiling. We denote by Ca∗ B {i ∈ [n] : σ∗ (i) = a} the set of indices of points in cluster a, and we
define ` B Ca∗ = nk .

B

Proof of Theorem 1

In this section, we prove Theorem 1, which relates the errors of the SDP and Oracle IP formulations.

B.1

Preliminaries

We may decompose the input matrix of pairwise squared distances as
A = C + C> − 2HH> ,
where H ∈ Rn×d is the matrix whose i-th row is the point hi , and C ∈ Rn×n is the matrix
where the entries in the i-th row are identical and equal to khi k22 . The row-sum constraint in
∗
b
Dthe program
E (4)
D ensures that
E the matrix Y − Y
D has zero row Esum. Consequently, we have
b − Y∗ , C = Y
b − Y∗ , C> = 0 which implies Y
b − Y∗ , C + C> = 0.
Y
Let G B H − EH be the centered version of H. We can compute
HH> = (G + EH) (G + EH) >
= GG> + G (EH) > + (EH) G> + (EH) (EH) >
and
EHH> = EGG> + (EH) (EH) > .
Therefore, we have


HH> − EHH> = GG> − EGG> + G (EH) > + (EH) G> .
Let U ∈ Rn×k be the matrix of the left singular vectors of Y∗ ; note that U is simply a scaled version
of the true assignment matrix F∗ and takes the form

1 ∗
1
Uia = √ Fia
= √ · 1 σ∗ (i) = a .
`
`

15

For each M ∈ Rn×n , define the projection PT (M) B UU> M + MUU> − UU> MUU> and its
b is optimal and Y∗ is feasible to
orthogonal complement PT ⊥ (M) B M − PT (M). The fact that Y
the program (4), implies that
E
1 Db
0≤− Y
− Y∗ , A
D2
E D
E
b − Y∗ , HH> − EHH> + Y
b − Y∗ , EHH>
= Y
E D
E
D
b − Y∗ , EHH>
b − Y∗ , GG> − EGG> + G (EH) > + (EH) G> + Y
= Y
E
E D
E D
(i) D
b − Y∗ , EHH>
b − Y∗ , G (EH) > + (EH) G> + Y
b − Y∗ , GG> + Y
= Y

E D

E
D
b − Y∗ , PT ⊥ GG>
b − Y∗ , PT GG> + Y
= Y
D
E D
E
b − Y∗ , G (EH) > + Y
b − Y∗ , EHH>
+2 Y
C S 1 + S 2 + 2S 3 + S 4 ,
∗ has zero diagonal. The following
b
where step (i) holds since EGG> is a diagonal matrix and Y−Y
propositions control the terms S 1 , S 2 and S 4 .
q
q

kd
k
kd
Proposition 1. If s2 ≥ C
log
n
+
log
n
+
n
n
n for some universal constant C > 0, then
1
2
100 ∆ γ

with probability at least 1 − n−8 .


1
Proposition 2. If s2 ≥ Ck dn + 1 for some universal constant C > 0, then S 2 ≤ 100
∆2 γ with
probability at least 1 − e−n/2 .


P
P
b − Y∗ for
Proposition 3. We have S 4 = − 12 a,b T ab ∆2ab ≤ − 41 ∆2 γ, where T ab B i∈Ca∗ , j∈Cb∗ Y
ij
each a, b ∈ [k].
S1 ≤

The proofs are given in Sections B.4, B.5 and B.6, respectively. Combining the above propositions, we have S 1 + S 2 ≤ − 12 S 4 and therefore
1
0 ≤ S3 + S4
4

(9)

0

with probability at least 1 − (2n)−C − 2e−n for some universal constant C 0 > 0.
b − Y∗ . We have
Let us take a closer look at the quantity S 3 + 41 S 4 . Let B B Y
XXX
D
E
S3 =
B ji µa , g j
j

a i∈Ca



XXD
E  1 X 
=`
µa , g j 
B ji 
`
a
j
i∈Ca∗


X X D
E  1 X 
=`
µa − µσ∗ ( j) , g j 
B ji  ,
`
j a,σ∗ ( j)
i∈C ∗
a

P



where the last step holds since a,σ∗ ( j) i∈Ca∗ B ji = − i∈Ca∗ :a=σ∗ ( j) B ji for each j ∈ [n], which
follows from the row-sum constraint of program (4). By Proposition 3, we have


 1 X 
X X 1
S 4 = −`
∆2σ∗ ( j),a 
B ji  .
2
`
∗
j a,σ∗ ( j)
i∈C
P

P

a

16

Therefore, together with the inequality (9), we obtain that
!
1
1
0≤
S3 + S4
`
4


X X D
E
  1 X 
=
B ji 
µa − µσ∗ ( j) , g j − c∆2σ∗ ( j),a 
`
∗
j a,σ∗ ( j)
i∈Ca


 1 X 
X X
B ji  ,
=
β ja 
`
∗
j a,σ∗ ( j)
i∈C

(10)

a

where c =

1
8

and β ja

D
E
B µa − µσ∗ ( j) , g j − c∆2σ∗ ( j),a .

b − Y∗ k1 = kBk1 and observe that the
To control the RHS of (10), we recall that γ B kY
constraints of the SDP (4) implies that


X X  1 X 
γ

B ji  =
∈ (0, n].
 `
2`
j∈[n] a,σ∗ ( j)
i∈C ∗
a

On the other hand, consider the linear program
X X
max
β ja X ja
X

j a,σ∗ ( j)

∀a , σ∗ ( j), j ∈ [n]

s.t. 0 ≤ X ja ≤ 1,
X
X ja ≤ 1,

∀ j ∈ [n]

(11)

a,σ∗ ( j)

X X

X ja = R.

j a,σ∗ ( j)

The above program is parameterized by the number R ∈ (0, n]. Let us denote by V(R) the optimal value of the above program (with the convention V(R) = −∞ if the program is infeasible).
Inspecting
 γ  the equation (10) and the program (11), we find that the RHS of (10) is upper bounded
by V 2`
. We therefore conclude that
!
γ
1
1
0≤
S3 + S4 ≤ V
.
(12)
`
4
2`

B.2

Controlling γ by a linear program

We next convert the inequality (12) into an upper bound on γ in terms of the objective value of an
LP that is related to the above program (11).
If γ = 0, then the conclusion of Theorem 1 holds trivially. For γ > 0, we consider the following
two cases:
1. If

γ
2`

∈ (0, 1], it follows from equation (12) that the error γ must satisfy
γ
γ
 γ 
γ
0≤V
= β∗ ≤ β∗
=V
,
2`
2`
2`
2`

where β∗ B max j∈[n],a,σ∗ ( j) β ja . This inequality implies that
γ
γ
≤
≤ max {R ∈ {0, 1, . . . .} : V(R) ≥ 0} .
2`
2`
17

γ
2`

> 1, it follows from equation (12) that the error γ must satisfy
γ
  γ   γ 
  γ   γ 

0≤V
≤ max V
,V
= max V
,V
−1 .
2`
2`
2`
2`
2`
This inequality implies that
γ
γ
≤
≤ max {R ∈ {0, 1, . . . .} : V(R) ∨ V(R − 1) ≥ 0}
2`
2`
= 1 + max {R ∈ {0, 1, . . . .} : V(R) ≥ 0} .
lγm
Note that 2` ≥ 2, and therefore we must have 1 ≤ max {R ∈ {0, 1, . . . .} : V(R) ≥ 0}. It
follows that
γ
≤ 2 max {R ∈ {0, 1, . . . .} : V(R) ≥ 0} .
2`

2. If

Combining the two cases, we obtain the
γ
≤ 2 max {R ∈ {0, 1, . . . .} : V(R) ≥ 0} .
(13)
2`
Note that the maximization on the RHS is in fact an LP, whose constraint is expressed in terms of
the optimal value V(R) of the previous LP (11).

B.3

Converting LP to IP

Equipped with the bound (13), we are now ready to formally establish a connection between
the error of the SDP (4) and that of the Oracle IP (6). We do so by relating the LP value
max {R ∈ {0, 1, . . . .} : V(R) ≥ 0} in (13), to the quantity (7) associated with the Oracle IP. n o
Concretely, note that for each integer R ∈ {0, 1, . . .}, there exists an optimal solution w ja of
program (11) such that w ja ∈ {0, 1} for all j ∈ [n], a ∈ [k]. In this case, we have
X X




max
β
X


ja
ja




X


∗


j a,σ ( j)










∗




[n]
s.t.
X
∈
{0,
1},
∀a
,
σ
(
j),
j
∈


ja




X


(14)
V(R) = IP1 (R) B 



[n]
X
≤
1,
∀
j
∈
ja








∗


a,σ ( j)




X X










X
=
R
ja






j a,σ∗ ( j)

for each integer R. Combining the last two equations (13) and (14), we obtain that
γ
≤ 2 max {R ∈ {0, 1, . . . .} : IP1 (R) ≥ 0}
2`




max R




R,X














s.t.
R
∈
{0,
1,
.
.
.}






X X








β
X
≥
0


ja
ja








∗
j


a,σ ( j)




=2·
∗



[n]
X
∈
{0,
1},
∀a
,
σ
(
j),
j
∈


ja




X










[n]
X
≤
1,
∀
j
∈
ja








∗


a,σ ( j)




X
X










X
=
R,


ja




∗
j a,σ ( j)

18



max



X










s.t.



= 2 · IP2 B 2 · 

















X X








j a,σ∗ ( j)



X X





β ja X ja ≥ 0



j a,σ∗ ( j)
.




∗

X ja ∈ {0, 1},
∀a , σ ( j), j ∈ [n]



X





X ja ≤ 1,
∀ j ∈ [n]



X ja

(15)

a,σ∗ ( j)

Let us reparameterize the integer program IP2 by a change of variable. Recall that
n
o
F B F ∈ {0, 1}n×k : F1k = 1n
is the set of all possible assignment matrices and F∗ ∈ F is the true assignment matrix; that is,
F ∗ja = 1 {a = σ∗ ( j)} for all j ∈ [n], a ∈ [k]. Consider any feasible solution X of IP2 ; here for
P
each j ∈ [n], we may fix X j,σ∗ ( j) = − a,σ∗ ( j) X ja — doing so does not affect the feasibility and
objective value of X with respect to IP2 . Define the new variable F B F∗ + X ∈ F . The objective
value and constraints of the old variable X can be mapped to those of F. In particular, we have
X X
j a,σ∗ ( j)

X ja =

1
(F ja − F ∗ja ) = kF − F∗ k1
2
a,σ∗ ( j)

X X
j

and

moreover,
X X
j a,σ∗ ( j)

(i)

β ja X ja =

X ja ∈ {0, 1}, ∀a , σ∗ ( j), j ∈ [n]
P
X ja ≤ 1, ∀ j ∈ [n]
a,σ∗ ( j)P
X j,σ∗ ( j) = − a,σ∗ ( j) X ja , ∀ j ∈ [n]






⇐⇒ F ∈ F ;





XX

XX

j

a

β ja X ja =

XX
j

β ja F ja −

a

j

(ii)

β ja F ∗ja =

a

XX
j

β ja F ja ,

a

where steps (i) and (ii) both follow from the fact that β j,σ∗ ( j) = 0, ∀ j. It follows that IP2 has the
same optimal value as another integer program IP3 with variable F; in particular, we have


1


∗




max
kF
−
F
k
1






F 2




X
X




s.t.
β
F
≥
0
.
IP2 = IP3 B 

ja
ja








a
j










F∈F
Plugging back to equation (15), we see that the error γ satisfies
γ
≤ 2 · IP3 .
2`

(16)

We further simplify the first constraint in IP3 . Recall that h̄i B µσ∗ (i) + (2c)−1 gi for each
n o
i ∈ [n]. Note that h̄i can be viewed as data points generated from the Sub-Gaussian Mixture
Model but with (2c)−1 times the standard deviation. By definition of β ja , we have
D
E
β ja = µa − µσ∗ ( j) , g j − c∆2σ∗ ( j),a
19

 D
E

= c 2 µa − µσ∗ ( j) , (2c)−1 g j − ∆2σ∗ ( j),a
 D
E

= c 2 µa − µσ∗ ( j) , (2c)−1 g j − kµa − µσ∗ ( j) k22
 D
E

= c 2 µa − µσ∗ ( j) , (2c)−1 g j − kµa − µσ∗ ( j) k22 − k(2c)−1 g j k22 + k(2c)−1 g j k22


= c −kµσ∗ ( j) − µa + (2c)−1 g j k22 + k(2c)−1 g j k22


= c −kh̄ j − µa k22 + k(2c)−1 g j k22 .
For any F ∈ F , we then have
XX
XX

β ja F ja = c
−kh̄ j − µa k22 + k(2c)−1 g j k22 F ja
j

a

a

j



 X X

X
X

2
−1
2

= c −
kh̄ j − µa k2 F ja +
k(2c) g j k2
F ja 
a
a
j
j




X
X
X
X

(i) 
= c −
kh̄ j − µa k22 F ja +
k(2c)−1 g j k22
F ∗ja 
a
a
j
j



 X X
XX
kh̄ j − µa k22 F ja +
k(2c)−1 g j k22 F ∗ja 
= c −
a
a
j
j



 X X
XX
kh̄ j − µa k22 F ja +
kh̄ j − µσ∗ ( j) k22 F ∗ja 
= c −
a
a
j
j


 X X

XX
(ii) 
kh̄ j − µa k22 F ja +
= c −
kh̄ j − µa k22 F ∗ja  ,
a

j

j

a

where step (i) holds because a F ja = 1 = a F ∗ja , ∀ j, and step (ii) holds because F ∗ja = 1 only if
a = σ∗ ( j). Recalling the shorthand
XX
kh̄ j − µa k22 F ja
η(F) B
P

P

j

a

introduced in Section 4.1, we have the more compact expression
XX

β ja F ja = c η(F∗ ) − η(F)
j

(17)

a

It follows that for any F ∈ F , the first constraint in IP3 is satisfied if and only if
η(F) ≤ η(F∗ ).
Combining with the bound (16), we obtain that


1




∗

max
kF
−
F
k


1 




F
2


γ


.
≤ 2 · IP3 = 2 · 
∗



s.t.
η(F)
≤
η(F
)


2`










F∈F
Rearranging terms, we have the bound

γ ≤ 2` · max kF − F∗ k1 : η(F) ≤ η(F∗ ), F ∈ F .
Theorem 1 then follows from the fact that kY∗ k1 = n` and kF∗ k1 = n.
20

(18)

B.4

Proof of Proposition 1

In this section we control S 1 . We can further decompose S 1 as
D

E D


E
b − Y∗ , UU> GG> + Y
b − Y∗ , GG> UU>
S1 = Y


E
D
b − Y∗ , UU> GG> UU>
− Y


E
D

E
D
b − Y∗ , UU> GG> UU>
b − Y∗ , UU> GG> + Y
≤2 Y
C 2T 1 + T 2
By the generalized Holder’s inequality, we have


T 1 ≤ γ · kUU> GG> k∞
and


E
D
b − Y∗ , UU> GG> UU>
T2 = Y
D


E
b − Y∗ UU> , UU> GG>
= Y


≤ γ · kUU> GG> k∞
where the last inequality holds since


b − Y∗ k1 = γ.
b − Y∗ UU> k1 ≤ kY
k Y
Combining the above, we have


S 1 ≤ 3γ · kUU> GG> k∞ .


Note that there are m = nk distinct random variables in UU> GG> and let us call them
D
E
P
X1 , . . . , Xm . For each i, we can see that Xi takes the form gu(i) , |M1 i | j∈Mi g j for some index
u(i) ∈ [n] and some set Mi ⊂ [n] such that |Mi | = `. Applying Lemma 7 with M = Mi , we have
for some universal constant C > 0
r

 d log n log n d 
2
Xi ≤ Cτ 
+ √ + 
`
`
`
for each i ∈ [n] with probability at least 1 − n−10 . The union bound therefore implies
r
r

 kd log n

k
kd
2
S 1 ≤ 3Cγτ 
+
log n + 
n
n
n
with probability at least 1 − n−8 . The result follows from the condition of the proposition.

B.5

Proof of Proposition 2

In this section we control S 2 . We have
D


E
b − Y∗ , GG>
S 2 = PT ⊥ Y
h

i
b − Y∗ · kGG> kop
≤ Tr PT ⊥ Y
21

γ
· kGG> kop
`
γ
= · kG> Gkop ,
`
≤

where the last step holds since
kGG> kop = max u> GG> u
u∈Sn−1

= max kG> uk22
=

u∈Sn−1
kG> kop

= kGkop
= max kGwk22
w∈Sd−1

= max wG> Gw
w∈Sd−1
>

= kG Gkop .
Observing that G> G =

>
i∈[n] gi gi

P

and applying Lemma 6, we have
!
γ
d
2
2
S 2 ≤ · 6τ (n + d) = 6γτ k 1 +
`
n

with probability at least 1 − e−n/2 . The result follows from the condition of the proposition.

B.6

Proof of Proposition 3

We can compute


EHH>


ij




if i , j and σ∗ (i) = σ∗ ( j)
Dkµσ∗ (i) k22
E
=

 µσ∗ (i) , µσ∗ ( j) if i , j and σ∗ (i) , σ∗ ( j).



b − Y∗ into k2 of ` × ` blocks, and note that T ab B Pi∈C ∗ , j∈C ∗ Y
b − Y∗
We partition the matrix Y
a
b
ij
denotes the sum of entries within the (a, b)-th block. The constraints of program (4) implies that
1. T aa ≤ 0 for each a ∈ [k] and T ab ≥ 0 for each a , b ∈ [k];
2. T ab = T ba for each a, b ∈ [k];
P
3. −T aa = b∈[k]:b,a T ab for each a ∈ [k];
P
P
P
P
4. − a∈[k] T aa + a,b∈[k]:a,b T ab = γ and thus − a∈[k] T aa = a,b∈[k]:a,b T ab = γ2 .
b − Y∗ has zero diagonal, we can write
Since Y
X
S4 =
T aa kµa k22 + 2
a∈[k]

=−

X
a,b∈[k]:a<b

X

T ab ∆2ab

a,b∈[k]:a<b

1 X
=−
T ab ∆2ab
2 a,b∈[k]:a,b
22

T ab µa , µb

≤−

1 X
T ab ∆2
2 a,b∈[k]:a,b

1
= − ∆2 γ.
4

C

Proof of Theorem 2

We define the shorthand
)
1
∗
∗
B max kF − F k1 : η(F) ≤ η(F ), F ∈ F .
2
(

γIP

It is not hard to see that γIP takes integer values in [0, n]. If γIP = 0 then we are done. We therefore
focus on the case γIP ∈ [n].
2
2
Suppose γIP > 3nke−s /C0 for a fixed C0 > D/c. Note that
2 /C 2
0

3nke−s

(i)

≤ nk ·

2
2
1 −s2 /(2C 2 )
0 ≤ ne−s /(2C 0 ) < n
·e
k

where step (i) holds since we have assumed s2 ≥ C s k for some universal constant C s > 0. We
record an important result for our proof.
Lemma 1. Let m ≥ 4 and g ≥ 1 be integers. Let X ∈ Rm×g be a matrix such that each X ja is a
sub-Gaussian random variable with its mean equal to λ ja and its sub-Gaussian norm no larger
 
than ρ ja , and each pair X ja and Xib are independent for j , i and a, b ∈ g . Then for some
universal constant D > 0 and for any β ∈ (0, m], we have
v
u

u
t X

X
X

λ ja M ja ,
X ja M ja ≤ D dβe  ρ2ja M ja  log (3mg/β) +
j,a

j,a

j,a

m×g

∀M ∈ {0, 1}
with probability at least 1 −

: M1g ≤ 1m , kMk1 = dβe ,

1.5
m .

The proof is given in Section D.2. Define the set
n
o
M B M ∈ {0, 1}n×k : M1k ≤ 1n , kMk1 = γIP , M j,σ∗ ( j) = 0 ∀ j ∈ [n] .
For any F feasible to IP3 , we have

1
η(F∗ ) − η(F)
c
X X
(i)
=
β ja F ja

0≤

j∈[n] a∈[k]

X

=

β ja

( j,a):F ja =1,a,σ∗ ( j)

≤ max

M∈M

X X

β ja M ja

j a,σ∗ ( j)

 v

u


t
 u



X
X
X
X






2
2
2
≤ max D γIP τ 
∆σ∗ ( j),a M ja  log (3n (k − 1) /γIP ) − c
∆σ∗ ( j),a M ja 
M∈M 

j a,σ∗ ( j)
j a,σ∗ ( j)

(ii)

23


 v
u


t

 u


X
X
X
X
2





 s
2
2
2
∆σ∗ ( j),a M ja 
≤ max D γIP τ 
∆σ∗ ( j),a M ja  2 − c
M∈M 

C0
j a,σ∗ ( j)
j a,σ∗ ( j)
!
X X
D
− c · max
∆2σ∗ ( j),a M ja
≤
C0
M∈M
j a,σ∗ ( j)
where step (i) holds
n o by equation (17), step (ii) holds by Lemma 1 with g = k − 1 since only
k − 1 entries of β ja are considered for each j in the sum above (ii), and the last step holds since
P P
P P
γIP ∆2 ≤ j a,σ∗ ( j) ∆2σ∗ ( j),a M ja . Since C0 > D/c and j a,σ∗ ( j) ∆2σ∗ ( j),a M ja > 0, the RHS above
2
2
2
2
is negative, which is a contradiction. Hence, we must have γIP ≤ 3nke−s /C0 ≤ ne−s /(2C0 ) and the
result follows from the fact that kF∗ k1 = n.

D

Proof of technical lemmas

In this section we prove the technical lemmas used in the proofs of our main theorems.

D.1

Proof of Lemma 7

We record the following lemma (Lemma A.3 in Lu and Zhou [2016]).
Lemma 2. Let {xi }i∈[n] be independent sub-Gaussian random vectors such that kxi kψ2 ≤ ρ for each
i ∈ [n]. For any fixed i ∈ [n], M ⊂ [n], t > 0 and δ > 0, we have
 √

*
( 2 )!

X + 3ρ2 t |M| + d + log (1/δ) 
1


 ≤ exp − min t , t + δ.
xj ≥
P  xi ,

4d 4
|M|
|M|
j∈M


p
Taking t = 4C d log n + log n and δ = n−20 completes the proof (where C > 0 is a sufficiently large universal constant).

D.2

Proof of Lemma 1

We define
LM B

X


X ja − λ ja M ja ,

j,a

Rβ,M

v
u

u
t X


B D dβe  ρ2ja M ja  log (3mg/β),
j,a

n
o
Mβ B M ∈ {0, 1}m×g : M1g ≤ 1m , kMk1 = dβe .
To establish a uniform bound in β, we apply a discretization argument to the possible values of β.
Define the shorthand E B (0, m]. We can cover E by the sub-intervals Et B (t − 1, t] for t ∈ [m].
For each t ∈ [m] we define the probability
n
o
αt B P ∃β ∈ Et , ∃M ∈ Mβ : LM > Rβ,M .
We bound each of these probabilities as follows:
24

(i) 
αt ≤ P ∃M ∈ Mt : LM > Rt,M








 [ 

≤ P
L
>
R

M
t,M




M∈M

t
X 
≤
P LM > Rt,M ,

(19)

M∈Mt

where step (i) holds since β ∈ Et implies β ≤ dβe = t.
Note that each X ja − λ ja is an independent zero-mean sub-Gaussian random variable and
P
the squared sub-Gaussian norm of LM is at most Cψ2 j,a ρ2ja M ja where Cψ2 > 0 is a universal
constant. We apply Hoeffding inequality (Lemma 5) to bound the probability on the RHS of (19):




2 t P ρ2 M




log(3mg/t)
cD


ja
j,a

ja


P LM > Rt,M ≤ exp 
−
P



2 M


C
ρ


ψ2
j,a ja ja

≤ exp −4t log(3mg/t) ,
where c > 0 is a universal constant. Plugging this back to (19), we have for each t ∈ [m],
X

αt ≤
exp −4t log(3mg/t)
M∈Mt

!

m t
=
g exp −4t log(3mg/t)
t
 me t

gt exp −4t log(3mg/t)
≤
t

≤ exp t log(3mg/t) + t − 4t log(3mg/t)
!t

t
≤ exp −t log(3mg/t) =
,
3mg
where the last inequality follows from t ≤ t log(3mg/t) for t ∈ [m]. It follows that
n
o
P ∃β ∈ E, ∃M ∈ Mβ : LM > Rβ,M
m


o


[ n

∃β ∈ Et , ∃M ∈ Mβ : LM > Rβ,M 
≤ P




t=1

≤
≤

m
X
t=1
m
X
t=1

αt
t
3mg

It remains to show that P1 (m) ≤

!t
C P1 (m).

1.5
m .

P1 (m) ≤

Since
m 
X
t t
3m
t=1
m

≤

X  t t
1
+
3m t=2 3m
25

(20)

 t t
1
+ m · max
,
t=2,3,...,m 3m
3m

≤

the proof is completed if for each integer t = 2, 3, . . . , m, we can show the bound
equivalently f (t) B t(log 3m − log t) ≥ 2 log m. Since t ≤ m, f (t) has derivative
!
3m
0
f (t) = log 3m − log t − 1 ≥ log 3m − log
− 1 = log 3 − 1 ≥ 0.
3




t t
3m

≤

1
,
m2

or

Therefore, f (t) is non-decreasing for 2 ≤ t ≤ m and therefore f (t) ≥ f (2) = 2 log 3m − 2 log 2 ≥
2 log m. Hence, P1 (m) ≤ 1.5
m .

E

Proof of Theorem 3

We only need to prove the first part of the theorem. The second part follows immediately from the
first part and Theorem 1.
The proof follows similar lines as those of Theorem 17 and Lemma 18 in Makarychev et al.
[2016]. In the rest of the section, we work under the context of Algorithms 1 and 2. Recall that
b − Y∗ k1 /kY∗ k1 . We have the following lemma.
k0 = {Bt }t≥1 and we let  B kY
Lemma 3. There exists a partial matching π0 between [k] and [k0 ] and a universal constant C > 0
such that
[
Ca∗ ∩ Bt ≥ (1 − C) n.
t=π0 (a)

The proof is given in Section E.1. The next lemma concerns the quality of clustering by
Algorithm 3.
Lemma 4. There exists a permutation π on [k] and a universal constant C > 0 such that
[

Ca∗ ∩ Ut ≥ (1 − C) n.

t=π(a)

The proof is given in Section E.2. The result follows from combining the above lemmas and
the fact that
[
1
err(b
σ, σ∗ ) = 1 − max
C ∗ ∩ Ut .
n π∈S k t=π(a) a

E.1

Proof of Lemma 3

We define ya to be an arbitrary row of Y∗ whose index is in Ca∗ .
(
)
`
∗
b
Ga B i ∈ Ca : kYi• − ya k1 ≤
,
∀a ∈ [k]
8
[
GB
Ga ,
a∈[k]

H B V\G.
We construct a partial matching π0 between sets Ca∗ and Bt by matching every cluster Ca∗ with
the first Bt that intersects Ga , and we let π0 (a) = t. Since each i ∈ [n] belongs to some Bt , we are
able to match every Ca∗ with some Bt . The fact that we cannot match two distinct clusters Ca∗ and
Cb∗ with the same Bt as well as the rest of the proof are given by the following fact.
26

Fact 2. We have
1. For each a ∈ [k] and t ∈ [k0 ] such that t = π0 (a), we have Bt ∩ Gb = ∅ for any b ∈ [k] \ {a}
and Bt ⊂ Ga ∪ H;
2. For each a ∈ [k] and t ∈ [k0 ] such that t = π0 (a), we have
Bt ∩ Ca∗ ≥ |Ga | − |Bt ∩ H| .
3. We have

X

Bt ∩ Ca∗ ≥ |V| − 2 |H| .

t=π0 (a)

4. There exists a universal constant C > 0 such that |H| ≤ Cn.
The proof is given below.
E.1.1

Proof of Fact 2

1. Suppose that there exist Bt and b ∈ [k] such that b , a and Bt ∩ Gb , ∅. Let u ∈ Bt ∩ Ga and
v ∈ Bt ∩ Gb . Since Ga and Gb are disjoint, we know that u , v. Let w ∈ Bt . Then we have
`
4
bv• − Y
bw• k1 ≤ ` .
kY
4

bu• − Y
bw• k1 ≤
kY

Therefore

bu• − Y
bv• k1 ≤ kY
bu• − Y
bw• k1 + kY
bv• − Y
bw• k1 ≤ ` .
kY
2

This implies
bu• k1 + kY
bu• − Y
bv• k1 + kyb − Y
bv• k1
kya − yb k1 ≤ kya − Y
` ` `
≤ + + < `,
8 2 8
which is a contradiction to the fact that kya − yb k1 = 2`. To complete the proof, we note that
for any i ∈ Bt we have either i ∈ Ga or i ∈ H.
2. Fix i ∈ Ga for some a ∈ [k]. For any j ∈ Ga we have j ∈ B(i) since
bi• − Y
b j• k1 ≤ kya − Y
bi• k1 + kya − Y
b j• k1 ≤ ` .
kY
4
Therefore, by definition
|Bt | ≥ |B(i)| ≥ |Ga | .
We have
(i)

Bt ∩ Ca∗ ≥ |Bt ∩ Ga |
= |Bt | − |Bt \Ga |
(ii)

= |Bt | − |Bt ∩ H|

≥ |Ga | − |Bt ∩ H| ,
where step (i) holds since Ga ⊂ Ca∗ and step (ii) holds since Bt ⊂ Ga ∪ H.
27

3. Summing the LHS of the above equation over t = π0 (a) gives
X
X
X
Bt ∩ Ca∗ =
|Bt ∩ H|
|Ga | −
t=π0 (a)

t=π0 (a)

a∈[k]

≥

X

|Ga | −

X

|Bt ∩ H|

t≥1

a∈[k]
(i)

= |G| − |V ∩ H|

= |V| − 2 |H| ,
S
where step (i) holds since Bt ∩ H are disjoint and t≥1 Bt = V.
4. We have
|H| ·

` X b
b − Y∗ k1 ≤ kY∗ k1 =  · n`
≤
kYi• − yσ∗ (i) k1 ≤ kY
8 i∈H

where the last step follows from the fact that kY∗ k1 = n`. The result follows.

E.2

Proof of Lemma 4

Let π0 be the partial matching between Ca∗ and Bt from Lemma 3. Define π(a) = π0 (a) for π0 (a) ≤ k.
If the resulting π is a partial permutation, we extend π to a permutation defined on [k] in an
arbitrary way. We may assume that {Ut }t∈[k] are {Bt }t∈[k] WLOG, and that Ut consists of Bt and
some elements from sets Bu with u > k. We have
[

Ca∗ ∩ Ut ≥

[

Ca∗ ∩ Bt

t=π0 (a)≤k

t=π(a)

=

[

[

Ca∗ ∩ Bt −

t=π0 (a)

Ca∗ ∩ Bt

t=π0 (a)>k

[


≥ 1 − C0 n −

Ca∗ ∩ Bt

t=π0 (a)>k

where C 0 > 0 is a universal constant. Define

T 1 B t > k : t = π0 (a) for some a ∈ [k] ,

T 2 B t ∈ [k] : t , π0 (a) for any a ∈ [k] .
Note that |T 1 | = |T 2 | and for any t1 ∈ T 1 and t2 ∈ T 2 we have Bt1 ≤ Bt2 . Therefore,
[

Ca∗ ∩ Bt ≤

t=π0 (a)>k

[

Bt

t∈T 1

≤

[

Bt

t∈T 2

≤ |V| −

[
t=π0 (a)

= C n.
0

The result follows by setting C B 2C 0 .
28

Ca∗ ∩ Bt

F

Proof of Fact 1

Let g = (g1 , . . . , gd )> be a random vector drawn from a rotationally invariant distribution supported
on the `2 ball in Rd . By Bilodeau and Brenner [2008, Proposition 4.10], g can be represented as
d
d
d
g = ru, where = means equality in distribution, r = kgk2 , u is uniformly distributed on the unit
d
sphere in R
q , and r and u are independent. The random vector u is sub-Gaussian with norm
kukψ2 ≤ C

1
d

[Vershynin, 2017, Theorem 3.4.5], so its one-dimensional margin satisfies
!
t2
P {|u1 | > t} ≤ 2 exp − 0
.
C /d
d

On the other hand, we have r = kgk2 ∈ [0, 1] since g is supported on the unit ball. Putting together
the above facts gives
!
t2
P {|g1 | > t} = P {r |u1 | > t} ≤ P {|u1 | > t} ≤ 2 exp − 0
,
C /d
q
d
whence kg1 kψ2 ≤ C d1 . By rotation invariance, we know that ha, gi = g1 for all unit vector a
[Bilodeau and Brenner, 2008, Proposition
q 4.8]. Therefore, all one-dimensional margin ha, gi of g
1
d,

is sub-Gaussian with norm at most C

G

which completes the proof.

Proof of Theorem 4

ba B i ∈ [n] : b
b, and recall
For a ∈ [k], we let C
σi = a denote the estimated clusters encoded in σ
P
that our cluster center estimators are defined by µ̂a B `−1 i∈Cba hi . For each a ∈ [k], we have
kµ̂a − µa k2 ≤ k

1X
1X
1X
hi −
h j k2 + k
h j − µa k2
`
` j∈C ∗
` j∈C ∗
ba
i∈C

a

a

C Q1 + Q2 .

(21)

We bound Q1 and Q2 separately below.

G.1

Controlling Q1

We work on the event that the result Theorem 3 is true. We have
X
1 X
hi −
h j k2
Q1 = k
`
∗
∗
ba \Ca
i∈C

ba
j∈Ca \C

ba \Ca∗ = Ca∗ \C
ba so we can pair each point in C
ba \Ca∗ with a point in Ca∗ \C
ba . Let us pair
Note that C
ba \Ca∗ with the j(i)-th point in Ca∗ \C
ba , and define M B {(i, j(i))}. We then have
the i-th point in C
∗
σ, σ ). Moreover, we can write
|M| ≤ n err(b
X 

1
Q1 = k
hi − h j(i) k2
`
(i, j(i))∈M
1 X
≤
khi − h j(i) k2
`
(i, j(i))∈M

29

≤
≤

1
`
1
`

X



∆σ∗ (i),σ∗ ( j(i)) + kgi − g j(i) k2

(i, j(i))∈M

X



Cq ∆ + kgi − g j(i) k2 ,

(i, j(i))∈M

where the last step holds for some universal constant Cq > 0 given that maxa,b∈[k] ∆ab ≤ Cq ∆.
To control the RHS of the last displayed equation, we first note that for each pair i , j ∈ [n],
the vector gi − g j is sub-Gaussian with sub-Gaussian norm at most Cψ2 τ where Cψ2 > 0 is a
universal constant. By Lemma 7, there exists a universal constant C > 0 such that
q
p

kgi − g j k2 ≤ Cτ2 d log n + log n + d
q

≤ 2Cτ2 log n + d
with probability at least 1 − n−11 . By the union bound over all pairs i , j, we get
q

max kgi − g j k2 ≤ 2Cτ2 log n + d
i, j

with probability at least 1 − n−9 . On the other hand, Theorem 3 guarantees that
" 2#
s
∗
σ, σ ) ≤ Cm exp −
|M| ≤ n · err(b
Ce
with probability at least 1−2n−1 . Combining the last two bounds, we obtain that for some universal
constants C0 > 0,
!
1
Q1 ≤ · Cq ∆ + max kgi − g j k2 · |M|
i, j
`
" 2#


p
s
≤ C0 ∆ + τ log n + d · k exp −
Ce
"
#


p
s2
≤ C0 ∆ + τ log n + d · exp −
2Ce
with probability at least 1 − 3n−1 , where the last step holds because s2 & k under the SNR condition (8). Note that the inequality e x ≥ x implies that
#
"
4Ce
τ 4Ce
τ
s2
≤ 2 = ·
≤ 4Ce ,
exp −
4Ce
∆
s
∆
s
where the last step holds since we have s ≥ 1 under the SNR condition (8). Hence, we have
"
#


p
s2
Q1 ≤ C0 τ 4Ce + log n + d · exp −
4Ce
"
#
2


p
s
≤ C1 τ 1 + log n + d · exp −
4Ce
"
#
p
s2
≤ 2C1 τ log n + d · exp −
4Ce
where C1 > 0 is a universal constant.
30

G.2

Controlling Q2

We have the expression
Q2 = k

1X
g j k2 .
` j∈C ∗
a

Applying Lemma 8, we obtain that
1
Q2 ≤ · 3τ
`

q

r

d + log n ` = 3τ

k d + log n
n



with probability at least 1 − n−3 .
Theorem 4 follows from plugging the above bounds on Q1 and Q2 into the inequality (21).

H

Technical lemmas

In this section, we state several standard tail bounds that are used in the proofs of our main theorems.
The following lemma is Theorem 2.6.2 in Vershynin [2017].
Lemma 5 (Hoeffding’s inequality for Sub-Gaussians). Let X1 , . . . , XN be independent, mean zero,
sub-Gaussian random variables. Then, for every t ≥ 0 we have
 N



2
 X



ct
 ,
P 
Xi ≥ t ≤ 2 exp − PN
2 
kX
k
i
i=1
ψ2
i=1
where c > 0 is a universal constant.
The following lemma is Lemma A.2 in Lu and Zhou [2016].
Lemma 6. We have
k

n
X

2
gi g>
i kop ≤ 6τ (n + d)

i=1

with probability at least 1 −

e−0.5n .

The following lemma is extended from Lemma A.3 in Lu and Zhou [2016].
Lemma 7. Let {xi }i∈[n] be independent sub-Gaussian random vectors such that kxi kψ2 ≤ ρ for each
i ∈ [n]. For any fixed i ∈ [n], M ⊂ [n], t > 0 and δ > 0, there exists a universal constant C > 0
such that
*
+ 3ρ2 5C √|M|  pd log n + log n + d
1 X
xi ,
xj ≤
|M|
|M|
j∈M

with probability at least 1 − n−10 .
The proof is given in Section D.1. The following lemma is Lemma A.4 in Lu and Zhou [2016].
Lemma 8. For any fixed a ⊂ [k], we have
q
X

k
gi k2 ≤ 3τ d + log n |Ca∗ |
i∈Ca∗

with probability at least 1 − n−3 .
31

References
Emmanuel Abbe. Community detection and the stochastic block model: recent developments.
Journal of Machine Learning Research, to appear, 2017. URL http://www.princeton.
edu/~eabbe/publications/sbm_jmlr_4.pdf.
Emmanuel Abbe and Colin Sandon. Community detection in general stochastic block models:
Fundamental limits and efficient algorithms for recovery. In IEEE 56th Annual Symposium on
Foundations of Computer Science (FOCS), pages 670–688. IEEE, 2015.
Emmanuel Abbe, Afonso S. Bandeira, and Georgina Hall. Exact recovery in the stochastic block
model. IEEE Transactions on Information Theory, 62(1):471–487, 2016.
Dimitris Achlioptas and Frank McSherry. On spectral learning of mixtures of distributions. In
International Conference on Computational Learning Theory, pages 458–469. Springer, 2005.
Daniel Aloise, Amit Deshpande, Pierre Hansen, and Preyas Popat. NP-hardness of euclidean
sum-of-squares clustering. Machine Learning, 75(2):245–248, 2009.
Brendan P. W. Ames and Stephen A. Vavasis. Convex optimization for the planted k-disjointclique problem. Mathematical Programming, 143(1-2):299–337, 2014.
Arash A. Amini and Elizaveta Levina. On semidefinite relaxations for the block model. The
Annals of Statistics, 46(1):149–179, 2018.
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.
Pranjal Awasthi and Aravindan Vijayaraghavan. Clustering semi-random mixtures of gaussians.
arXiv preprint arXiv:1711.08841, 2017.
Pranjal Awasthi, Afonso S. Bandeira, Moses Charikar, Ravishankar Krishnaswamy, Soledad Villar, and Rachel Ward. Relax, no need to round: Integrality of clustering formulations. In
Proceedings of the 2015 Conference on Innovations in Theoretical Computer Science, pages
191–200. ACM, 2015.
Sivaraman Balakrishnan, Martin J. Wainwright, and Bin Yu. Statistical guarantees for the em
algorithm: From population to sample-based analysis. The Annals of Statistics, 45(1):77–120,
2017.
Martin Bilodeau and David Brenner. Theory of multivariate statistics. Springer Science & Business Media, 2008.
Moses Charikar, Sudipto Guha, Éva Tardos, and David B. Shmoys. A constant-factor approximation algorithm for the k-median problem. In Proceedings of the 31st Annual ACM Symposium
on Theory of Computing, pages 1–10. ACM, 1999.
Yudong Chen, Sujay Sanghavi, and Huan Xu. Improved graph clustering. IEEE Transactions on
Information Theory, 60(10):6440–6455, 2014.
Sanjoy Dasgupta. Learning mixtures of gaussians. In 40th Annual Symposium on Foundations of
Computer Science, pages 634–644. IEEE, 1999.
32

Constantinos Daskalakis, Christos Tzamos, and Manolis Zampetakis. Ten steps of em suffice for
mixtures of two gaussians. arXiv preprint arXiv:1609.00368, 2016.
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.
Yingjie Fei and Yudong Chen. Exponential error rates of SDP for block models: Beyond
Grothendieck’s inequality. arXiv preprint arXiv:1705.08391, 2017.
Olivier Guédon and Roman Vershynin.
Community detection in sparse networks via
Grothendieck’s inequality. Probability Theory and Related Fields, 165(3-4):1025–1049, 2016.
Takayuki Iguchi, Dustin G. Mixon, Jesse Peterson, and Soledad Villar. Probably certifiably correct
k-means clustering. Mathematical Programming, 165(2):605–642, 2017.
Kamal Jain, Mohammad Mahdian, and Amin Saberi. A new greedy approach for facility location
problems. In Proceedings of the Thiry-Fourth Annual ACM Symposium on Theory of Computing, pages 731–740. ACM, 2002.
Chi Jin, Yuchen Zhang, Sivaraman Balakrishnan, Martin J. Wainwright, and Michael I. Jordan.
Local maxima in the likelihood of gaussian mixture models: Structural results and algorithmic
consequences. In Advances in Neural Information Processing Systems, pages 4116–4124, 2016.
Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman, and
Angela Y. Wu. A local search approximation algorithm for k-means clustering. Computational
Geometry, 28(2-3):89–112, 2004.
Jason M. Klusowski and W. D. Brinda. Statistical guarantees for estimating the centers of a twocomponent gaussian mixture by em. arXiv preprint arXiv:1608.02280, 2016.
Michael Krivelevich and Dan Vilenchik. Semirandom models as benchmarks for coloring algorithms. In Third Workshop on Analytic Algorithmics and Combinatorics (ANALCO), pages
211–221, 2006.
Amit Kumar and Ravindran Kannan. Clustering with spectral norm and the k-means algorithm.
In Proceedings of the 2010 IEEE 51st Annual Symposium on Foundations of Computer Science,
pages 299–308. IEEE Computer Society, 2010.
Shi Li and Ola Svensson. Approximating k-median via pseudo-approximation. SIAM Journal on
Computing, 45(2):530–547, 2016.
Xiaodong Li, Yang Li, Shuyang Ling, Thomas Strohmer, and Ke Wei. When do birds of a feather
flock together? K-means, proximity, and conic programming. arXiv preprint arXiv:1710.06008,
2017.
Stuart Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28
(2):129–137, 1982.
Yu Lu and Harrison H. Zhou. Statistical and computational guarantees of Lloyd’s algorithm and
its variants. arXiv preprint arXiv:1612.02099, 2016.

33

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.
Konstantin Makarychev, Yury Makarychev, and Aravindan Vijayaraghavan. Learning communities in the presence of errors. In 29th Annual Conference on Learning Theory, pages 1258–1291,
2016.
Dustin G. Mixon, Soledad Villar, and Rachel Ward. Clustering subgaussian mixtures by semidefinite programming. Information and Inference: A Journal of the IMA, 6(4):389–415, 2017.
Andrea Montanari and Subhabrata Sen. Semidefinite programs on sparse random graphs and
their application to community detection. In Proceedings of the 48th Annual ACM SIGACT
Symposium on Theory of Computing (STOC), pages 814–827, 2016.
Abhinav Nellore and Rachel Ward. Recovery guarantees for exemplar-based clustering. Information and Computation, 245:165–180, 2015.
Samet Oymak and Babak Hassibi. Finding dense clusters via "low rank+ sparse" decomposition.
arXiv preprint arXiv:1104.5186, 2011.
Karl Pearson. Method of moments and method of maximum likelihood. Biometrika, 28(1/2):
34–59, 1936.
Jiming Peng and Yu Wei. Approximating k-means-type clustering via semidefinite programming.
SIAM Journal on Optimization, 18(1):186–205, 2007.
Jiming Peng and Yu Xia. A new theoretical framework for k-means-type clustering. In Foundations and Advances in Data Mining, pages 79–96. Springer, 2005.
Santosh Vempala and Grant Wang. A spectral algorithm for learning mixture models. Journal of
Computer and System Sciences, 68(4):841–860, 2004.
Roman Vershynin. High-dimensional probability. 2017.
Ji Xu, Daniel J. Hsu, and Arian Maleki. Global analysis of expectation maximization for mixtures
of two gaussians. In Advances in Neural Information Processing Systems, pages 2676–2684,
2016.
Bowei Yan, Mingzhang Yin, and Purnamrita Sarkar. Convergence analysis of gradient EM for
multi-component gaussian mixture. arXiv preprint arXiv:1705.08530, 2017.

34

