Clustering subgaussian mixtures by semidefinite programming
Dustin G. Mixon1 , Soledad Villar2 , and Rachel Ward2

arXiv:1602.06612v2 [stat.ML] 10 May 2016

1

Department of Mathematics and Statistics, Air Force Institute of Technology
2
Department of Mathematics, University of Texas at Austin

Abstract
We introduce a model-free relax-and-round algorithm for k-means clustering based on a
semidefinite relaxation due to Peng and Wei [PW07]. The algorithm interprets the SDP output
as a denoised version of the original data and then rounds this output to a hard clustering.
We provide a generic method for proving performance guarantees for this algorithm, and we
analyze the algorithm in the context of subgaussian mixture models. We also study the fundamental limits of estimating Gaussian centers by k-means clustering in order to compare our
approximation guarantee to the theoretically optimal k-means clustering solution.

1

Introduction

Consider the following mixture model: For each t ∈ [k] := {1, . . . , k}, let Dt be an unknown
subgaussian probability distribution over Rm , with first moment γt ∈ Rm and second moment
matrix with largest eigenvalue σt2 . For each t, an unknown number nt of random points {xt,i }i∈[nt ]
is drawn independently from Dt . Given the points {xt,i }i∈[nt ],t∈[k] along with the model order k,
the goal is to approximate the centers {γt }t∈[k] . How large must ∆ := mina6=b kγa − γb k2 be relative
to σmax := maxt σt , and how large must nmin := mint nt be relative to nmax := maxt nt , in order to
have sufficient signal for successful approximation?
For the most popular instance of this problem, where the subgaussian distributions are Gaussians, theoretical guarantees date back to the work of Dasgupta [Das99]. Dasgupta introduced
an algorithm based on random projections and showed that this algorithm well-approximates cen√
ters of Gaussians in Rm that are separated by σmax m. Since Dasgupta’s seminal work, performance guarantees for several algorithmic alternatives have emerged, including expectation maximization [DS07], spectral methods [VW04, KK10, AS12], projections (random and deterministic)
[MV10, AK05], and the method of moments [MV10]. Every existing performance guarantee has
one two forms:
(a) the algorithm correctly clusters all points according to Gaussian mixture component, or
(b) the algorithm well-approximates the center of each Gaussian (a la Dasgupta [Das99]).
Results of type (a), which include [VW04, KK10, AS12, AM07], require the minimum separation
P
between the Gaussians centers to have a multiplicative factor of polylog N , where N = kt=1 nt
is the total number of points. This stems from a requirement that every point be closer to their
Gaussian’s center (in some sense) than the other centers, so that the problem of cluster recovery
is well-posed. We note that in the case of spherical Gaussians, such highly separated Gaussian
1

components may be truncated so as to match a different data model known as the stochastic ball
model, where the semidefinite program we use in this paper is already known to be tight with high
probability [ABC+ 15, IMPV15].
Results of type (b) tend to be specifically tailored to exploit unique properties of the Gaussian
distribution, and are thus not easily generalizable to other data models. For instance, if x has
distribution N (µ, σ 2 Im ), then E(kx − µk2 ) = mσ 2 , and concentration of measure implies that in
√
high dimensions, most of the points will reside in a thin shell with center µ and radius about mσ.
This sort of behavior can be exploited to cluster even concentric Gaussians as long as the covariances
are sufficiently different. However, algorithms that perform well even with no separation between
centers require a sample complexity which is exponential in k [MV10].
In this paper, we provide a performance guarantee of type (b), but our approach is model-free. In
particular, we consider the k-means clustering objective:
minimize

k X
X

1 X
xi −
xj
|At |

t=1 i∈At

subject to

j∈At

2

(1)
2

A1 ∪ · · · ∪ Ak = {1, . . . , N },

Ai ∩ Aj = ∅

∀i, j ∈ [k], i 6= j

Letting D denote the N × N matrix defined entrywise by Dij = kxi − xj k22 , then a straightforward
calculation gives the following “lifted” expression for the k-means objective:
k X
X
t=1 i∈At

1 X
xi −
xj
|At |
j∈At

2
2

1
= Tr(DX),
2


Xij =

1
|At |

0

if i, j ∈ At
otherwise

(2)

The matrix X necessarily satisfies various convex constraints, and relaxing to such constraints leads
to the following semidefinite relaxation of (1), first introduced by Peng and Wei in [PW07]:
minimize

Tr(DX)

(3)

subject to

Tr(X) = k, X1 = 1, X ≥ 0, X  0

Here, X ≥ 0 means that X is entrywise nonnegative, whereas X  0 means that X is symmetric
and positive semidefinite.
As mentioned earlier, this semidefinite relaxation is known to be tight for a particular data model
called the stochastic ball model [NW15, ABC+ 15, IMPV15]. In this paper, we study its performance
under subgaussian mixture models, which include the stochastic ball model and the Gaussian
mixture model as special cases. The SDP is not typically tight under this general model, but the
optimizer can be interpreted as a denoised version of the data and can be rounded in order to
produce a good estimate for the centers (and therefore produce a good clustering).
To see this, let P denote the m × N matrix whose columns are the points {xt,i }i∈[nt ],t∈[k] . Notice
that whenever X has the form (2), then for each t ∈ [k], P X has |At | columns equal to the centroid
of points assigned to At . In particular, if X is k-means-optimal, then P X reports the k-meansoptimal centroids (with appropriate multiplicities). Next, we note that every SDP-feasible matrix
X ≥ 0 satisfies X > 1 = X1 = 1, and so X > is a stochastic matrix, meaning each column of P X
is still a weighted average of columns from P . Intuitively, if the SDP relaxation (3) were close to
being tight, then the SDP-optimal X would make the columns of P X close to the k-means-optimal
centroids. Empirically, this appears to be the case (see Figure 1 for an illustration). Overall, we
2

(a)

(b)

(c)

(d)

Figure 1: (a) Draw 100 points at random from each of three spherical Gaussians over R2 . These points form the
columns of a 2 × 300 matrix P . (b) Compute the 300 × 300 distance-squared matrix D from the data in (a), and solve
the k-means semidefinite relaxation (3) using SDPNAL+v0.3 [YST15]. (The computation takes about 16 seconds on
a standard MacBook Air laptop.) Given the optimizer X, compute P X and plot the columns. We interpret this as a
denoised version of the original data P . (c) The points in (b) land in three particular locations with particularly high
frequency. Take these locations to be estimators of the centers of the original Gaussians. (d) Use the estimates for
the centers in (c) to partition the original data into three subsets, thereby estimating the k-means-optimal partition.

may interpret P X as a denoised version of the original data P , and we leverage this strengthened
signal to identify good estimates for the k-means-optimal centroids.
What follows is a summary of our relax-and-round procedure for (approximately) solving the kmeans problem (1):

Relax-and-round k-means clustering procedure.
Given and m × N data matrix P = [x1 · · · xN ], do:
(i) Compute distance-squared matrix D defined by Dij = kxi − xj k22 .
(ii) Solve k-means semidefinite program (3), resulting in optimizer X.
(iii) Cluster the columns of the denoised data matrix P X.

For step (iii), we find there tends to be k vectors that appear as columns in P X with particularly
high frequency, and so we are inclined to use these as estimators for the k-mean-optimal centroids
(see Figure 1, for example). Running Lloyd’s algorithm for step (iii) also works well in practice.
To obtain theoretical guarantees, we instead find the k columns of P X for which the unit balls of
a certain radius centered at these points in Rm contain the most columns of P X (see Theorem 15
for more details). An implementation of our procedure is available on GitHut [Vil16].
Our contribution. We study performance guarantees for the k-means semidefinite relaxation (3)
when the point cloud is drawn from a subgaussian mixture model. Adapting ideas from Guédon
and Vershynin [GV14], we obtain approximation guarantees comparable with the state of the art
for learning mixtures of Gaussians despite the fact that our algorithm is a generic k-means solver
and uses no model assumptions. To the best of our knowledge, no convex relaxation has been
used before to provide theoretical guarantees for clustering mixtures of Gaussians. We also provide
conditional lower bounds on how well a k-means solver can approximate the centers of Gaussians.
Organization of this paper. In Section 2, we present a summary of our results and give a
high-level explanation of our proof techniques. We also illustrate the performance of our relax3

and-round algorithm on the MNIST dataset of handwritten digits. Our theoretical results consist
of an approximation theorem for the SDP (proved in Section 3), a denoising consequence of the
approximation (explained in Section 4), and a rounding step (presented in Section 5). We also study
the fundamental limits for estimating Gaussian centers by k-means clustering (see Section 2.2).

2

Summary of results

This paper has two main results. First, we present a relax-and-round algorithm for k-means clustering that well-approximates the centers of sufficiently separated subgaussians. Second, we provide
a conditional result on the minimum separation necessary for Gaussian center approximation by
k-means clustering. The first result establishes that the k-means SDP (3) performs well with noisy
data (despite not being tight), and the second result helps to illustrate how sharp our analysis is.
This section discusses these results, and then applies our algorithm to the MNIST dataset [LC10].

2.1

Performance guarantee for the k-means SDP

Our relax-and-round performance guarantee consists of three steps.
Step 1: Approximation. We adapt an approach used by Guédon and Vershynin to provide
approximation guarantees for a certain semidefinite program under the stochastic block model for
graph clustering [GV14].
Given the points xt,1 , . . . , xt,nt drawn independently from Dt , consider the squared-distance matrix
D and the corresponding minimizer XD of the SDP (3). We first construct a “reference” matrix R
such that the SDP (3) is tight when D = R with optimizer XR . To this end, take ∆ab := kγa −γb k2 ,
let XD denote the minimizer of (3), and let XR denote the minimizer of (3) when D is replaced by
the reference matrix R defined as

(4)
(Rab )ij := ξ + ∆2ab /2 + max 0, ∆2ab /2 + 2hra,i − rb,j , γa − γb i
where rt,i := xt,i − γt , and ξ > 0 is a parameter to be chosen later. Indeed, this choice of reference
is quite technical, as an artifact of the entries in D being statistically dependent. Despite its lack
of beauty, our choice of reference enjoys the following important property:
Lemma 1. Let 1a ∈ RN denote the indicator function for the indices i P
corresponding to points xi
drawn from the ath subgaussian. If γa 6= γb whenever a 6= b, then XR = kt=1 (1/nt )1t 1>
t .
Proof. Let X be feasible for the the SDP (3). Replacing D with R in (3), we may use the SDP
constraints X1 = 1 and X ≥ 0 to obtain the bound
Tr(RX) =

N X
N
X
i=1 j=1

Rij Xij ≥

N X
N
X

ξXij =

i=1 j=1

N
N
X
X
ξ
Xij = N ξ = Tr(RXR )
i=1

j=1

Furthermore, since γa 6= γb whenever a 6= b, and since X ≥ 0, we have that equality occurs precisely
for the X such that (Xab )ij equals zero whenever a 6= b. The other constraints on X then force XR
to have the claimed form (i.e., XR is the unique minimizer).

4

Now that we know that XR is the solution we want, it remains to demonstrate regularity in the
sense that XD ≈ XR provided the subgaussian centers are sufficiently separated. For this, we use
the following scheme:
• If hR, XD i ≈ hR, XR i then kXD − XR k2F is small (Lemma 8).
• If D ≈ R (in some specific sense) then hR, XD i ≈ hR, XR i (Lemmas 9 and 10).
• If the centers are separated by O(kσmax ), then D ≈ R.
What follows is the result of this analysis:
Theorem 2. Fix , η > 0. There exist universal constants C, c1 , c2 , c3 such that if
α = nmax /nmin . k . m

and

N > max{c1 m, c2 log(2/η), log(c3 /η)},

then kXD − XR k2F ≤  with probability ≥ 1 − 2η provided
∆2min ≥

C 2 2
k ασmax


where ∆min = mina6=b kγa − γb k2 is the minimal cluster center separation.
See Section 3 for the proof. Note that if we remove the assumption α . k . m, we obtain the
2 .
result ∆2min ≥ C (min{k, m} + α)kασmax
Step 2: Denoising. Suppose we solve the SDP (3) for an instance of the subgaussian mixture
model where ∆min is sufficiently large. Then Theorem 2 ensures that the solution XD is close to
the ground truth. At this point, it remains to convert XD into an estimate for the centers {γt }t∈[k] .
Let P denote the m × N matrix whose (a, i)th column is xa,i . Then P XR is an m × N matrix
whose (a, i)th column is γ̃a , the centroid of the ath cluster, which converges to γa as N → ∞, (and
does not change when i varies, for a fixed a), and so one might expect P XD to have its columns be
close to the γt ’s. In fact, we can interpret the columns of P XD as a denoised version of the points
(see Figure 1).
To illustrate this idea, assume the points {xa,i }i∈[n] come from N (γa , σ 2 Im ) in Rm for each a ∈ [k].
Then we have
 X

k X
n
1
2
E
kxa,i − γa k2 = mσ 2 .
(5)
N
a=1 i=1

Letting ca,i denote the (a, i)th column of P XD (i.e., the ith estimate of γa ), in Section 4 we obtain
the following denoising result:
Corollary 3. If kσ . ∆min ≤ ∆max . Kσ, then
k
n
1 XX
kca,i − γ̃a k22 . K 2 σ 2
N
a=1 i=1

with high probability as n → ∞.
√
Note that Corollary 3 guarantees denoising in the regime K  m. This is a corollary of a more
technical result (Theorem 11), which guarantees denoising for certain configurations of subgaussians
(e.g., when the γt ’s are vertices of a regular simplex) in the regime k  m.
5

At this point, we comment that one might expect this level of denoising from principal component
analysis (PCA) when the mixture of subgaussians is sufficiently nice. To see this, suppose we have
spherical Gaussians of equal entrywise variance σ 2 centered at vertices of a regular simplex. Then
in the large-sample limit, we expect PCA to approach the (k − 1)-dimensional affine subspace that
contains the k centers. Projecting onto this affine subspace will not change the variance of any
Gaussian in any of the principal components, and so one expects the mean squared deviation of
the projected points from their respective Gaussian centers to be (k − 1)σ 2 .
By contrast, we find that in practice, the SDP denoises substantially more than PCA does. For
example, Figures 1 and 3 illustrate cases in which PCA would not change the data, since the data
already lies in (k − 1)-dimensional space, and yet the SDP considerably enhances the signal. In
fact, we observe empirically that the matrix XD has low rank and that P XD has repeated columns.
This doesn’t come as a complete surprise, considering SDP optimizers are known to exhibit low
rank [Sha82, Bar95, Pat98]. Still, we observe that the optimizer tends to have rank O(k) when
clustering points from the mixture model. This is not predicted by existing bounds, and we did
not leverage this feature in our analysis, though it certainly warrants further investigation.
Step 3: Rounding. In Section 5, we present a rounding scheme that provides a clustering of the
original data from the denoised results of the SDP (Theorem 15). In general, the cost of rounding
is a factor of k in the average squared deviation of our estimates. Under the same hypothesis as
Corollary 3, we have that there exists a permutation π on {1, . . . , k} such that
k
1X
kvi − γ̃π(i) k22 . kK 2 σ 2 ,
k

(6)

i=1

where vi is what our algorithm chooses as the ith center estimate. Much like the denoising portion,
we also have a more technical result that allows one to replace the right-hand side of (6) with k 2 σ 2
for sufficiently nice configurations of subgaussians. As such, we can estimate Gaussian centers with
mean squared error O(k 2 σ 2 ) provided the centers have pairwise distance Ω(kσ). This contrasts
with the seminal work of Dasgupta [Das99], which gives mean squared error O(mσ 2 ) when the
√
pairwise distances are Ω( mσ). As such, our guarantee replaces dimensionality dependence with
√
model order–dependence, which is an improvement to the state of the art in the regime k  m.
In the next section, we indicate that model order–dependence cannot be completely removed when
using k-means to estimate the centers.
Before concluding this section, we want to clarify the nature of our approximation guarantee (6).
Since centroids correspond to a partition of Euclidean space, our guarantee says something about
how “close” our k-means partition is to the “true” partition. By contrast, the usual approximation
guarantees for relax-and-round algorithms compare values of the objective function (e.g., the kmeans value of the algorithm’s output is within a factor of 2 of minimum). Also, the latter sort
of optimal value–based approximation guarantee cannot be used to produce the sort of optimizerbased guarantee we want. To illustrate this, imagine a relax-and-round algorithm for k-means that
produces a near-optimal partition with k = 2 for data coming from a single spherical Gaussian.
We expect every subspace of co-dimension 1 to separate the data into a near-optimal partition, but
the partitions are very different from each other when the dimension m ≥ 2, and so a guarantee of
the form (6) will not hold.

6

2.2

Fundamental limits of k-means clustering

In Section 5, we provide a rounding scheme that, when applied to the output of the k-means
SDP, produces estimates of the subgaussian centers. But how good is our guarantee? Observe the
following two issues: (i) The amount of tolerable noise σ and our bound on the error maxi kvi −
γ̃π(i) k2 both depend on k. (ii) Our bound on the error does not vanish with N .
In this section, we give a conditional result that these issues are actually artifacts of k-means; that
is, both of these would arise if one were to estimate the Gaussian centers with the k-means-optimal
centroids (though these centroids might be computationally difficult to obtain). The main trick in
our argument is that, in some cases, so-called “Voronoi means” appear to serve as a good a proxy
for the k-means-optimal centroids. This trick is useful because the Voronoi means are much easier
to analyze. We start by providing some motivation for the Voronoi means.
(X )

m
Given X = {xi }N
i=1 ⊆ R , let A1
objective

(X )

t · · · t Ak
k X
X

= {1, . . . , N } denote any minimizer of the k-means
1 X
xj
|At |

xi −

t=1 i∈At

j∈At

2

,
2

and define the k-means-optimal centroids by
(X )

ct

:=

1

X

(X )
|At |
(X )
j∈At

xj .

(Note that the k-means minimizer is unique for generic X .) Given Γ = {γt }kt=1 , then for each γa ,
consider the Voronoi cell
n
o
Va(Γ) := x ∈ Rm : kx − γa k2 < kx − γb k2 ∀b 6= a .
Given a probability distribution D over Rm , define the Voronoi means by

(Γ,D)
(Γ) 
µt
:= E X X ∈ Vt .
X∼D

Finally, we say Γ ⊆ Rm is a stable isogon if
(si1) |Γ| > 1,
(si2) the symmetry group G ≤ O(m) of Γ acts transitively on Γ, and
(si3) for each γ ∈ Γ, the stabilizer Gγ has the property that
{x ∈ Rm : Qx = x ∀Q ∈ Gγ } = span{γ}.
(Examples of stable isogons include regular and quasi-regular polyhedra, as well as highly symmetric
frames [BW13].) With this background, we formulate the following conjecture:
Conjecture 4 (Voronoi Means Conjecture). Draw N points X independently from a mixture D
of equally weighted spherical Gaussians of equal variance centered at the points in a stable isogon
Γ = {γt }kt=1 . Then
(X )
(Γ,D)
min
max
ct − µπ(t) 2
π : [k]→[k] t∈{1,...,k}
permutation

converges to zero in probability as N → ∞, i.e., the k-means-optimal centroids converge to the
Voronoi means.
7

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

-1

-1

-1

-1.5

-1.5

-1.5

-2
-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

(a)

Figure 2:

-2
-2

-1.5

-1

-0.5

0

0.5

1

1.5

-2
-2

2

-1.5

(b)

-1

-0.5

0

0.5

1

1.5

2

(c)

Evidence in favor of the Voronoi Means Conjecture. For each n ∈ {102 , 103 , 104 }, the following

experiment produces (a), (b) and (c), respectively. Perform 30 trials of the following: For each γ ∈ {(±1, 0), (0, ±1)},
draw n points independently from N (γ, I), and then for these N = 4n points, run MATLAB’s built-in implementation
of k-means++ with k = 4 for 10 independent initializations; of the resulting 10 choices of centroids, store the ones
of minimal k-means value (this serves as our proxy for the k-means-optimal centroids). Plot these 30 collections
of centroids in black dots, along with the Voronoi means in red circles. The Voronoi means {(±α, 0), (0, ±α)}
were computed numerically in Mathematica as α ≈ 1.39928. Importantly, the k-means-optimal centroids appear to
converge toward the Voronoi means, not the Gaussian centers, as N → ∞.

Our conditional result will only require the Voronoi Means Conjecture in the special case where
the γt ’s form an orthoplex (see Lemma 6). Figure 2 provides some numerical evidence in favor
of this conjecture (specifically in the orthoplex case). We note that the hypothesis that Γ be a
stable isogon appears somewhat necessary. For example, simulations suggest that the conjecture
does not hold when Γ is three equally spaced points on a line (in this case, the k-means-optimal
centroids appear to be more separated than the Voronoi means). The following theorem provides
some insight as to why the stable isogon hypothesis is reasonable:
Theorem 5. Let D be a mixture of equally weighted spherical Gaussians of equal variance centered
(Γ,D)
at the points of a stable isogon Γ = {γt }kt=1 . Then there exists α > 0 such that µt
= αγt for
each t ∈ {1, . . . , k}.
See Section 6 for the proof. To interpret Theorem 5, consider k-means optimization over the
distribution D instead of a large sample X drawn from D. This optimization amounts to finding k
points C = {ct }kt=1 in Rm that minimize
k
X
t=1

E

X∼D


(C) 
kX − ct k22 X ∈ Vt
P

X∼D

(C) 

X ∈ Vt

(7)

Intuitively, the optimal C is a good proxy for the k-means-optimal centroids when N is large (and
one might make this rigorous using the plug-in principle with the Glivenko–Cantelli Theorem).
What Theorem 5 provides is that, when Γ is a stable isogon, the Voronoi means have the same
Voronoi cells as do Γ. As such, if one were to initialize Lloyd’s algorithm at the Gaussian centers to
solve (7), the algorithm converges to the Voronoi means in one step. Overall, one should interpret
Theorem 5 as a statement about how the Voronoi means locally minimize (7), whereas the Voronoi
Means Conjecture is a statement about global minimization.
8

As indicated earlier, we will use the Voronoi Means Conjecture in the special case where Γ is an
orthoplex:
Lemma 6. The standard orthoplex of dimension d, given by the first d columns of I and of −I, is
a stable isogon.
Proof. First, we have that |Γ| = 2d > 1, implying (si1). Next, every γ 0 ∈ Γ can be reached from
any γ ∈ Γ with the appropriate signed transposition, which permutes Γ, and is therefore in the
symmetry group G; as such, we conclude (si2). For (si3), pick γ ∈ Γ and let i denote its nonzero
entry. Consider the matrix Q = 2ei e>
i − I, where ei denotes the ith identity basis element. Then
Q ∈ Gγ , and the eigenspace of Q with eigenvalue 1 is span{γ}, and so we are done.
What follows is the main result of this subsection:
Theorem 7. Let k ≤ 2m be even, and let Γ = {γt }kt=1 ⊆ Rm denote the standard orthoplex of
dimension k/2. Then for every σ > 0, either
σ . ∆min /

p
log k

or

min
t∈{1,...,k}

(Γ,D)

kµt

− γt k2 & σ

p
log k,

where D denotes the mixture of equally weighted spherical Gaussians of entrywise variance σ 2
centered at the members of Γ.
See Section 7 for the proof. In words, Theorem 7 establishes that one must accept k-dependence in
either the data’s noise or the estimate’s error. It would be interesting to investigate whether other
choices of stable isogons lead to stronger k-dependence.

2.3

Numerical example: Clustering the MNIST dataset

In this section, we apply our clustering algorithm to the NMIST handwritten digits dataset [LC10].
This dataset consists of 70,000 different 28×28 grayscale images, reshaped as 784×1 vectors; 55,000
of them are considered training set, 10,000 are test set, and the remaining 5,000 are validation set.
Clustering the raw data gives poor results (due to 4’s and 9’s being similar, for example), so we
first learn meaningful features, and then cluster the data in feature space. To simplify feature
extraction, we used the first example from the TensorFlow tutorial [AAB+ 15]. This consists of a
one-layer neural network y(x) = σ(W x+b), where σ is the softmax function, W is a 784×10 matrix
to learn, and b is a 10 × 1 vector to learn. As the tutorial shows, the neural network is trained for
1,000 iterations, each iteration using batches of 100 random points from the training set.
Training the neural network amounts to finding W and b that fit the training set well. After
selecting these parameters, we run the trained neural network on the first 1,000 elements of the
test set, obtaining {y(xi )}1000
i=1 , where each y(xi ) is a 10 × 1 vector representing the probabilities of
being each digit. Since y(xi ) is a probability vector, its entries sum to 1, and so the feature space
is actually 9-dimensional.
For this experiment, we cluster {y(xi )}1000
i=1 with two different algorithms: (i) MATLAB’s builtin implementation of k-means++, and (ii) our relax-and-round algorithm based on the k-means
semidefinite relaxation (3). (The results of the latter alternative are illustrated in Figure 3.)

9

(a)

(b)

(c)

(a) After applying TensorFlow [AAB+ 15] to learn a 9-dimensional feature space of MNIST digits [LC10],
determine the features of the first 1,000 images in the MNIST test set, compute the 1000 × 1000 matrix D of squared
distances in feature space, and then solve the k-means semidefinite relaxation (3) using SDPNAL+v0.3 [YST15].
(The computation takes about 6 minutes on a standard MacBook Air laptop.) Convert the SDP-optimizer X to a
grayscale image such that white pixels denote zero entries. By inspection, this matrix is not exactly of the form (2),
but it looks close, and it certainly appears to have low rank. (b) Letting P denote the 9 × 1000 matrix whose
columns are the feature vectors to cluster, compute the denoised data P X and identify the 10 most popular locations
in R9 (denoted by red circles) among the columns of P X (denoted by black dots). For the plot, we project the
9-dimensional data onto a random 2-dimensional subspace. (c) The 10 most popular locations form our estimates
for the centers of digits in feature space. We plot these locations relative to the original data, projected in the same
2-dimensional subspace as (b).

Figure 3:

Since each run of k-means++ uses a random initialization that impacts the partition, we ran this
algorithm 100 times. In fact, the k-means value of the output varied quite a bit: the all-time low
was 39.1371, the all-time high was 280.4174, and the median was 108.2358; the all-time low was
reached in 34 out of the 100 trials. Since our relax-and-round alternative has no randomness, the
outcome is deterministic, and its k-means value was 39.1371, i.e., identical to the all-time low from
k-means++. By comparison, the k-means value of the planted solution (i.e., clustering according to
the hidden digit label) was 103.5430, and the value of the SDP (which serves as a lower bound on the
optimal k-means value) was 38.5891. As such, not only did our relax-and-round alternative produce
the best clustering that k-means++ could find, it also provided a certificate that no clustering has
a k-means value that is 1.5% better.
Recalling the nature of our approximation guarantees, we also want to know well the relax-andround algorithm’s clustering captures the ground truth. To evaluate this, we determined a labeling
of the clusters for which the resulting classification exhibited a minimal misclassification rate. (This
amounts to minimizing a linear objective over all permutation matrices, which can be relaxed to
a generically tight linear program over doubly stochastic matrices.) For k-means++, the all-time
low misclassification rate was 0.0971 (again, accomplished by 34 of the 100 trials), the all-time high
was 0.4070, and the median was 0.2083. As one might expect, the relax-and-round output had a
misclassification rate of 0.0971.

3

Proof of Theorem 2

By the following lemma, it suffices to bound Tr(R(XD − XR )):

10

Lemma 8. kXD − XR k2F ≤

5
Tr(R(XD − XR )).
nmin ∆2min

Proof. First, by Lemma 1, we have kXR k2F = k. We also claim that kXD k2F ≤ k. To see this,
first note that XD 1 = 1 and XD ≥ 0, and so the ith entry of XD v can be interpreted as a convex
combination of the entries of v. Let v be an eigenvector of XD with eigenvalue µ, and let i index
the largest entry of v (this entry is positive without loss of generality). Then µvi = (XD v)i ≤ vi ,
implying that µ ≤ 1. Since the eigenvalues of XD lie in [0, 1], we may conclude that kXD k2F ≤
Tr(XD ) = k. As such,
kXD − XR k2F = kXD k2F + kXR k2F − 2 Tr(XD XR )
≤ 2k − 2 Tr(XD XR )
= 2k + 2 Tr((XR − XD )XR ) − 2kXR k2F
= 2 Tr((XR − XD )XR ).

(8)

We will bound (8) in two different ways, and a convex combination of these bounds will give the
result. For both bounds, we let Ω denote the indices in the diagonal blocks, and Ωc the indices
in the off-diagonal blocks, and Ωt ⊂ Ω denote the indices in the diagonal block for the cluster
t. In particular, AΩ denotes the matrix that equals A on the diagonal blocks and is zero on the
off-diagonal blocks. For the first bound, we use that RΩ = ξ(11> )Ω , and that (XR − XD )Ω (11T )Ω
has non-negative entries (since both XR and XD have non-negative entries, XR 1 = XD 1 = 1, and
XR = (XR )Ω ). Recalling that RΩ = ξ, we have


k
X
1
>
2 Tr((XR − XD )XR ) =
2 Tr (XR − XD )(11 )Ωt
nt
t=1


2
≥
Tr (XR − XD )(11> )Ω
nmax
2
Tr ((XD − XR )RΩ )
(9)
=−
ξnmax

P 
For the second bound, we first write nmin XR = 11> − (11> )Ωc − kt=1 1 − nnmin
(11> )Ωt . Since
t
XR 1 = 1 = XD 1, we then have
!!

k 
X
2
n
min
2 Tr((XR − XD )XR ) =
Tr (XD − XR ) (11> )Ωc +
1−
(11> )Ωt − 11>
nmin
nt
t=1
!!

k 
X
nmin
2
> c
>
Tr (XD − XR ) (11 )Ω +
1−
(11 )Ωt
=
nmin
nt
t=1

≤

2

>

Tr((XD − XR )(11 )Ωc )
nmin
2
=
Tr(XD (11> )Ωc ),
nmin
where the last and second-to-last steps use that (XR )Ωc = 0. Next, XD ≥ 0 and RΩc ≥ (ξ +
∆2min /2)(11> )Ωc , and so we may continue:
2
2 Tr((XR − XD )XR ) ≤
Tr(XD RΩc )
nmin (ξ + ∆2min /2)
2
=
Tr((XD − XR )RΩc ),
(10)
nmin (ξ + ∆2min /2)
11

where again, the last step uses the fact that (XR )Ωc = 0. At this point, we have bounds of the form
x ≥ ay1 with a < 0 and x ≤ by2 with b > 0 (explicitly, (9) and (10)), and we seek a bound of the
form x ≤ c(y1 +y2 ). As such, we take the convex combination for a, b such that a−1 /(a−1 +b−1 ) < 0
and b−1 /(a−1 + b−1 ) > 0
x≤

a−1
b−1
1
ay
+
by2 = −1
(y1 + y2 ).
1
−1
−1
−1
−1
a +b
a +b
a + b−1

Taking a = −2/(ξnmax ) and b = 2/(nmin (ξ + ∆2min /2)) and combining with (8) then gives
kXD − XR k2F

nmin 2 −1
≤ 2 Tr((XR − XD )XR ) ≤
(nmin − nmax ) +
∆min
Tr((XD − XR )(RΩ + RΩc )),
2
4
ξ

choosing ξ > 0 sufficiently small and simplifying yields the result.
We will bound Tr(R(XD − XR )) in terms of the following: For each N × N real symmetric matrix
M , let F(M ) denote the value of the following program:
F(M ) = maximum
subject to

| Tr(M X)|

(11)

Tr(X) = k, X1 = 1, X ≥ 0, X  0

e := P1⊥ DP1⊥ and R̃ := P1⊥ RP1⊥ . Then Tr(R(XD − XR )) ≤ 2F(D̃ − R̃).
Lemma 9. Put D
Proof. Since XD and XR are both feasible in (11), we have
e D ) + Tr(RX
e D ) ≤ | Tr((D
e − R)X
e D )| ≤ F(D
e − R),
e
− Tr(DX
e R ) − Tr(RX
e R ) ≤ | Tr((D
e − R)X
e R )| ≤ F(D
e − R),
e
Tr(DX
and adding followed by reverse triangle inequality gives

 

e − R)
e ≥ Tr(DX
e R ) − Tr(DX
e D ) + Tr(RX
e D ) − Tr(RX
e R) .
2F(D

(12)

Write XDe := P1⊥ XD P1⊥ . Note that XD 1 = (XD )T 1 implies XD = XDe + (1/N )11> , and so


e D ) = Tr(DX e ) = Tr D XD − (1/N )11> = Tr(DXD ) + 1 1> D1.
Tr(DX
D
N
e R ) = Tr(DXR ) +
Similarly, Tr(DX

1 >
N 1 D1,

and so

e R ) − Tr(DX
e D ) = Tr(DXR ) − Tr(DXD ) ≥ 0,
Tr(DX
e D )−Tr(RX
e R ) = Tr(R(XD −
where the last step follows from the optimality of XD . Similarly, Tr(RX
XR )), and so (12) implies the result.
e − R).
e For an n1 × n2 matrix X, consider the matrix norm
Now it suffices to bound F(D
kXk1,∞ :=

n1
X
i=1

max |Xi,j | =

1≤j≤n2

The following lemma will be useful:
12

n1
X
i=1

kXi,. k∞ .

Lemma 10. F(M ) ≤ min {kM k1,∞ , min{k, r}kM k2→2 } where r = rank(M ).
Proof. The first bound follows from the classical version of Hölder’s inequality (recalling that Xi,j ≥
0 and X1 = 1 by design):


N X
N
N
N
N
X
X
X
X
| Tr(M X)| ≤
|Mi,j Xi,j | ≤
kMi,. k∞ 
|Xi,j | =
kMi,. k∞
i=1 j=1

i=1

j=1

i=1

The second bound is a consequence of Von Neumann’s trace inequality: if the singular values of X
and M are respectively α1 ≥ . . . ≥ αN and β1 ≥ . . . ≥ βN then

| Tr(M X)| ≤

N
X

αi βi

i=1

Since X is feasible in (11) we have α1 ≤ 1 and

| Tr(M X)| ≤

k
X

PN

i=1 αi

≤ k. Using that rank(M ) = r we get

βi ≤ min{k, r}kM k2→2

i=1

Proof of Theorem 2. Write xt,i = rt,i + γt . Then
(Dab )ij = kxa,i − xb,j k22
= k(ra,i + γa ) − (rb,j + γb )k22 = kra,i − rb,j k22 + 2hra,i − rb,j , γa − γb i + kγa − γb k22 .
Furthermore,
kra,i − rb,j k22 = kra,i k22 − 2hra,i , rb,j i + krb,j k22 = ((µ1> + G> G + 1µ> )ab )ij ,
where G is the matrix whose (a, i)th column is ra,i , and µ is the column vector whose (a, i)th entry
is kra,i k22 . Recall that

(Rab )ij = ξ + ∆2ab /2 + max 0, ∆2ab /2 + 2hra,i − rb,j , γa − γb i
Then P1⊥ (D − R)P1⊥ = P1⊥ G> GP1⊥ + P1⊥ F P1⊥ where
 2
∆ab /2 + 2hra,i − rb,j , γa − γb i if 2hra,i − rb,j , γa − γb i ≤ −∆2ab /2
(Fab )ij =
0
otherwise.
Considering Lemma 10 and that rank(G> G) ≤ m we will bound
F(M ) ≤ min{k, m}kP1⊥ G> GP1⊥ k2→2 +

1
nmin

kP1⊥ F P1⊥ k1,∞ .

For the first term:
kP1⊥ G> GP1⊥ k2→2 ≤ kG> Gk2→2 = kG> k22→2 .
13

(13)

(t)

Note that if the rows Xi , i = 1, . . . nt of G> come from a distribution with second moment matrix
(t)
1/2
Σt , then Xi has the same distribution as Σt g, where g is an isotropic random vector. Then
kG> k ≤ σmax kG̃> k where the rows of G̃> are isotropic random vectors.
By Theorem 5.39 in [Ver12], we have that there exist c1 and c2 constants depending only on the
subgaussian norm of the rows of G such that with probability ≥ 1 − η:
√

p
√
kG> k2→2 ≤ σmax N + c1 m + c2 log(2/η) .
Note that by Corollary 3.35, when the rows of G> are Gaussian random vectors we have the result
for c1 = 1 and c2 = 2.
For bounding the second term in (13), the triangle inequality gives kP1⊥ F P1⊥ k1,∞ ≤ 4kF k1,∞ .
In order to get a handle on kF k1,∞ we first compute the expected value of its entries using that
|2hra,i − rb,j , γa − γb i| obeys a folded subgaussian distribution, coming from a subgaussian with
2 ∆2 :
variance at most 8σmax
ab


E|(Fab )ij | ≤ ∆2ab /2 + E|2hra,i − rb,j , γa − γb i| P 2hra,i − rb,j , γa − γb i < −∆2ab /2
 2



∆ab 4σmax ∆ab
∆2ab
√
≤
+
exp −
2
2
64σmax
π


∆2ab
2
≤ ∆2ab exp −
assuming ∆2min > 16kσmax
, k≥2
2
64σmax
4
642 σmax
1
using e−x ≤ 2 for x > 0.
≤ ∆2ab
4
x
∆ab
2
256σmax
2
using again ∆2min > 16kσmax
,
k
2
= O(σmax
/k)

≤−

k≥2

Now we can write F = 2(L − L> ) where La,i := (Lab )ij ∈ {hra,i , γa − γb i, 0} has independent rows,
2 /k). We can then bound
and E|(Lab )ij | ≤ E|(Fab )ij | = O(σmax
kF k1,∞ ≤ 4kLk1,∞ ≤ kLsmall k1,1
where Lsmall ∈ RN ×k is a submatrix of distinct columns.
Then we have a high-probability estimate:
P(kLsmall k1,1 > t) ≤ P 2k

na
k X
X

!
|La,i | > t

≤P

a=1 i=1

na
k X
X
a=1 i=1

t
2
(|La,i | − E|La,i |) >
− c3 σmax
nmax
2k

!

Using that La,i are independent subgaussian random variables, we know there exist constants
c4 , c5 ≥ 0 such that

P

na
k X
X

!
(|La,i | − E|La,i |) > u

a=1 i=1
2
So, choosing t = 2c3 knmax σmax
+

q

N
c5



u2
≤ c4 exp −c5
N

log cη4 , we get that with probability at least 1 − η
14

s
kP1⊥ F P1⊥ k1,∞ ≤

2
8c3 knmax σmax

+4

c4
N
log
c5
η

Putting everything together, we get that there exist constants C1 , C2 , C3 such that with probability
at least 1 − 2η
5
Tr(R(XD − XR ))
nmin ∆2min
10
e − R)
e
≤
F(D
nmin ∆2min
√
2
p
√
p
2
min{k, m}
N + c1 m + c2 log(2/η) σmax
2
N log c4 /η
knmax σmax
≤ C1
+ C2
+ C3
.
2
2
nmin ∆min
nmin ∆min
nmin ∆2min

kXD − XR k2F ≤

If additionally we require N > max{c1 m, c2 log(2/η), log(c4 /η)}, we get
kXD − XR k2F ≤ C

2 (α + min{k, m})
kασmax
∆2min

Rearranging gives the result.

4

Denoising

In the special case where each Gaussian is spherical with the same entrywise variance σ 2 and the
same number n of samples, the main result is says:
kXD − XR k2F .

k2 σ2
∆2min

with high probability as n → ∞.
Let P denote the m × N matrix whose (a, i)th column is xa,i . Then P XR is an m × N matrix whose
(a, i)th column is γ̃a , a good estimate of γa , and so one might expect P XD to have its columns be
close to the γ̃a ’s. This is precisely what the following theorem gives:
√
Theorem 11. Suppose σ . ∆min / k. Let P denote the m × N matrix whose (a, i)th column is
xa,i , and let ca,i denote the (a, i)th column of P XD . Then
k
n
1 XX
kΓk22→2
kca,i − γ̃a k22 .
· kσ 2
2
N
∆
min
a=1 i=1

with high probability as n → ∞. Here, the ath column of Γ is γ̃a −

1
k

Pk

b=1 γ̃b .

The proof can be found at the end of this section. For comparison,
 X

k X
n
1
2
kxa,i − γa k2 = mσ 2 ,
E
N
a=1 i=1

15

(14)

meaning the ca,i ’s serve as “denoised” versions of the xa,i ’s provided kΓk2→2 is not too large compared to ∆min . The following lemma investigates this provision:
Lemma 12. For every choice of {γ̃a }ka=1 , we have
kΓk22→2
1
≥ ,
2
2
∆min
with equality if {γ̃a }ka=1 is a simplex. More generally, if the following are satisfied simultaneously:
Pk
(i)
a=1 γ̃a = 0,
(ii) kγ̃a k2  1 for every a ∈ {1, . . . , k}, and
(iii) |hγ̃a , γ̃b i| . 1/k for every a, b ∈ {1, . . . , k} with a 6= b,
then

kΓk22→2
. 1.
∆2min

See the end of the section for the proof. Plugging these estimates for kΓk22→2 /∆2min into Theorem 11
shows that the ca,i ’s in this case exhibit denoising to an extent that the m in (14) can be replaced
with k:
k
n
1 XX
kca,i − γ̃a k22 . kσ 2 .
N
a=1 i=1

{γ̃a }ka=1 ,

For more general choices of
one may attempt to estimate kΓk2→2 in terms of ∆max , but
this comes with a bit of loss in the denoising estimate:
Corollary 13. If kσ . ∆min ≤ ∆max . Kσ, then
k
n
1 XX
kca,i − γ̃a k22 . K 2 σ 2
N
a=1 i=1

with high probability as n → ∞.
Indeed, this doesn’t guarantee denoising unless k . K ≤
following string of inequalities to Theorem 11:

√

m. To prove this corollary, apply the

kΓk22→2 ≤ kΓk2F ≤ k∆2max . kK 2 σ 2 ,
where the second inequality uses the following lemma:
P
Lemma 14. If ka=1 γ̃a = 0, then kγ̃a k2 ≤ ∆max for every a.
Proof. Fix a. Then



X

k 
k
1 X
1
γ̃a
1
γ̃a
γ̃a
min γ̃b ,
≤
γ̃b ,
=
γ̃b − γ̃a ,
=−
kγ̃a k2 .
b6=a
kγ̃a k2
k−1
kγ̃a k2
k−1
kγ̃a k2
k−1
b=1
b6=a

b=1

Let b(a) denote the minimizer. Then Cauchy–Schwarz gives


γ̃a
∆max ≥ kγ̃a − γ̃b(a) k2 ≥ γ̃a − γ̃b(a) ,
≥ kγ̃a k2 +
kγ̃a k2
16

1
k−1 kγ̃a k2

≥ kγ̃a k2 .

Proof of Theorem 11. Without loss of generality, we have
k X
n
X

Pk

a=1 γ̃a

= 0. Write

kca,i − γ̃a k22 = kP (XD − XR )k2F ≤ kP k22→2 kXD − XR k2F .

(15)

a=1 i=1

Decompose P = Γ ⊗ 1> + G, where 1 is n-dimensional and G has i.i.d. entries from N (0, σ 2 ).
Observe that
kΓ ⊗ 1> k22→2 = k(Γ ⊗ 1> )(Γ ⊗ 1> )> k2→2 = knΓΓ> k2→2 = nkΓk22→2 .

(16)

Also, Corollary 5.35 in [Ver12] gives that
√
√
√
kGk2→2 . ( N + m)σ . N σ

(17)

with probability ≥ 1 − e−Ωm (N ) . The result then follows from estimating kP k2→2 with (16) and
(17) by triangle inequality, plugging into (15), and then applying Theorem 2.
Proof of Lemma 12. Since kΓxk2 ≤ kΓk2→2 kxk2 for every x, we have that
kΓk22→2 ≥

kγ̃a − γ̃b k22
2

for every a and b, and so
1 ∆2max
1
kΓk22→2
≥
· 2 ≥ .
2
2 ∆min
2
∆min
For the second part, let {γ̃a }ka=1 be a simplex. Without loss of generality, {γ̃a }ka=1 is centered at
the origin, each point having unit 2-norm. Then hγ̃1 , γ̃2 i = −1/(k − 1), and so
∆2min = kγ̃1 − γ̃2 k22 = kγ̃1 k22 + kγ̃2 k22 − 2hγ̃1 , γ̃2 i =

2k
.
k−1

Next, we write
Γ> Γ =

k
1
I−
11> ,
k−1
k−1

and conclude that kΓk22→2 = kΓ> Γk2→2 = k/(k − 1). Combining with our expression for ∆2min then
gives the result. For the last part, pick a and b such that ∆min = kγ̃a − γ̃b k2 . Then
∆2min = kγ̃a k22 + kγ̃b k22 − 2hγ̃a , γ̃b i & 2 − 2/k.
Also, Gershgorin implies
kΓk22→2 = kΓ> Γk2→2 . 1 + (k − 1)/k,
and so combining these estimates gives the result.

5

Rounding

Theorem 15. Take  < ∆min /8, suppose
n
o n
# (a, i) : kca,i − γ̃a k2 >  < ,
2
17

and consider the graph G of vertices {ca,i }ni=1, ka=1 such that ca,i ↔ cb,j if kca,i − cb,j k2 ≤ 2. For
each i = 1, . . . , k, select the vertex vi of maximum degree (breaking ties arbitrarily) and update G by
removing every vertex w such that kw − vi k2 ≤ 4. Then there exists a permutation π on {1, . . . , k}
such that
kvi − γ̃π(i) k2 ≤ 3
for every i ∈ {1, . . . , k}.
Proof. Let B(x, r) denote the closed 2-ball of radius r centered at x. For each i, we will determine
π(i) at the conclusion of iteration i. Denote R1 := {1, . . . , k} and Ri+1 := Ri \ {π(i)} for each
i = 2, . . . , k − 1. We claim that the following hold at the beginning of each iteration i:
[
(i) < n/2 vertices lie outside
B(γ̃a , ),
a∈Ri

(ii) ≥ n/2 vertices lie inside B(vi , 2), and
(iii) there exists a unique a ∈ Ri such that kvi − γ̃a k2 ≤ 3.
First, we show that for each i, (i) and (ii) together imply (iii). Indeed, there are enough vertices in
B(vi , 2) that one of them must reside in B(γ̃a , ) for some a ∈ Ri . Furthermore, this a is unique
since  < ∆min /6. By triangle inequality, we have kvi − γ̃a k2 ≤ 3, and so we put π(i) := a.
We now prove (i) and (ii) by induction. When i = 1, we have (i) by assumption. For (ii), note that
each B(γ̃a , ) contains ≥ n/2 of the vertices, and by triangle inequality, each has degree ≥ n/2 − 1
in G. As such, the vertex v1 of maximum degree will have degree ≥ n/2 − 1, thereby implying (ii).
Now suppose (i), (ii) and (iii) all hold for iteration i < k. By triangle inequality, (iii) implies
B(γ̃π(i) , ) ⊆ B(vi , 4). As such, the ith iteration removes all vertices in B(γ̃π(i) , ) so that (i)
continues to hold for iteration i + 1. Next,  < ∆min /8 and (iii) together imply that the removal of
vertices in B(vi , 4) preserves the vertices in B(γ̃a , ) for every a ∈ Ri+1 , and their degrees are still
≥ n/2 − 1 by the same triangle argument as before. Thus, (ii) holds for iteration i + 1.
Corollary 16. Suppose k . m, and denote S := kΓk2→2 /∆min . Pick   Skσ. Perform the
rounding scheme of Theorem 15 to columns of P XD . Then with high probability, {vi }ki=1 satisfies
kvi − γ̃π(i) k2 . Skσ
for some permutation π, provided σ . ∆min /(Sk).
By Lemma 12, we have S . 1 in the best-case scenario. In this case, our rounding scheme
√ works
in the regime σ . ∆min /k. (Note that denoising is guaranteed in the regime σ . ∆min / k). In
general, the cost of rounding is a factor of k in the average squared deviation of our estimates:
n
k
1 XX
kca,i − γ̃a k22 . S 2 kσ 2 ,
N

whereas

a=1 i=1

k
1X
kvi − γ̃π(i) k22 . S 2 k 2 σ 2 .
k
i=1

On the other hand, we are not told which of the points in {ca,i }ni=1, ka=1 correspond to any given γ̃a ,
whereas in rounding, we know that each vi corresponds to a distinct γ̃a .
Proof of Corollary 16. Draw (a, i) uniformly from {1, . . . , k} × {1, . . . , n} and take X to be the

18

random variable kca,i − γ̃a k22 . Then Markov’s inequality and Theorem 11 together give
n
o
# (a, i) : kca,i − γ̃a k2 >  = N · P(X > 2 )
≤

k
n
N 1 XX
N
n
·
kca,i − γ̃a k22 . 2 · S 2 kσ 2 . .
2 N

2
a=1 i=1

For Theorem 15 to apply, it suffices to ensure  < ∆min /8, which follows from σ . ∆min /(Sk).

6

Proof of Theorem 5

Lemma 17. Let G be the symmetry group of a stable isogon Γ ⊆ Rm , and let K and H denote the
subgroups of G that fix W = span(Γ) and its orthogonal complement W ⊥ , respectively. Then
(i)
(ii)
(iii)
(iv)

G is the direct sum of H and K,
H is finite,
H acts transitively on Γ, and
K is isomorphic to the orthogonal group O(m − r), where r is the dimension of W .

Proof. Pick Q ∈ G. Then Q permutes the points in Γ, and the permutation completely determines
how Q acts on W by linearity. In particular, W is invariant under the action of G, which in turn
implies the same for W ⊥ . Let V denote an m × r matrix whose columns form an orthonormal basis
for W , and let V⊥ denote an m × (m − r) matrix whose columns form an orthonormal basis for
W ⊥ . Then Q can be expressed as




V>
A 0
Q = V V⊥
0 B
V⊥>
We see that Q ∈ H when B = I and Q ∈ K when A = I. Let QH denote the “projection” of Q
onto H, obtained by replacing B with I, and similarly for QK .
For (i), it suffices to show that H and K are normal subgroups of G, that H ∩ K = {I}, and
that G is generated by H and K. The first is obtained by observing that K is the kernel of the
homomorphism Q 7→ QH , and similarly, H is the kernel of Q 7→ QK . The second follows from the
observation that Q ∈ H ∩ K implies A = I and B = I. The last follows from the observation that
every Q ∈ G can be factored as QH QK .
For (ii), we note that A is completely determined by how Q permutes the points in Γ, of which
≤ k! possibilities are available.
For (iii), we know that by (si2), for every pair γ, γ 0 ∈ Γ, there exists Q ∈ G such that Gγ = γ 0 .
Consider the factorization Q = QH QK . Since QK γ = γ, we therefore have QH γ = γ 0 , meaning H
also acts transitively on Γ.
For (iv), we first note that B ∈ O(m − r) is necessary in order to have Q ∈ O(m). Now pick any
B ∈ O(m − r). Then




I 0
V>
(18)
Q = V V⊥
0 B
V⊥>
has the effect of fixing each point in Γ, meaning Q ∈ G (and therefore Q ∈ K). As such, K is the
set of all Q of the form (18).
19

Lemma 18. For any stable isogon Γ, we have

P

γ∈Γ γ

= 0.

Proof. By (si1), we have |Γ| > 1. In the special case where |Γ| = 2 and the points in Γ are
linearly dependent, write Γ = {γ1 , γ2 } with γ1 6= γ2 . By (si2), we know there exists Q ∈ G such
that Qγ1 = γ2 , and so kγ1 k2 = kγ2 k2 . This combined with the assumed linear dependence gives
γ1 = ±γ2 , and since γ1 6= γ2 , we conclude γ1 + γ2 = 0, as desired.
In the remaining case, Γ contains two points (say, γ1 and γ2 ) that are linearly independent. Fix
γ0 ∈ Γ. By Lemma 17(iii), the orbit Hγ0 is all of Γ. Consider the map from H onto Γ given by
f : U 7→ U γ0 for all U ∈ H. The preimage of any member of Γ is a left coset of the stabilizer Hγ0
(which is finite by Lemma 17(ii)), and so
X
X
X X
γ.
(19)
U γ0 =
U γ0 = |Hγ0 |
γ∈Γ U ∈f −1 (γ)

U ∈H

γ∈Γ

Now pick any Q ∈ G, and consider the factorization Q = QH QK with QH ∈ H and QK ∈ K
(this exists uniquely by Lemma 17(i)). Since U γ0 ∈ span(Γ) for every U ∈ H, we have that
QK U γ0 = U γ0 , and so
X
X
X
Q
U γ0 =
QH U γ 0 =
U γ0 ,
(20)
U ∈H

U ∈H

U ∈H

where the last step
P follows from the fact that multiplication by QH ∈ H permutes the members
of H. Put x = U ∈H U γ0 . Then (20) gives that Qx = x for every Q ∈ G, and therefore Qx = x
for every Q ∈ Gγ1 ∪ Gγ2 . By (si3), this then implies that x ∈ span{γ1 } ∩ span{γ2 }, i.e., x = 0.
Combining with (19) then gives the result.
Lemma 19. Let D be a mixture of equally weighted spherical Gaussians of equal variance centered
at the points of a stable isogon Γ = {γt }kt=1 . If X ∼ D and Q is any member of the symmetry group
G of Γ, then QX ∼ D.
Proof. The probability density function of D is given by
f (x) =

k
1X
1
2
2
√
e−kx−γt k2 /2σ .
k
( 2πσ)m
t=1

As such, f (Qx) = f (x) since kQx − γt k2 = kx − Q−1 γt k2 and Q−1 ∈ G permutes the γt ’s.
(Γ)

Proof of Theorem 5. Pick Q ∈ Gγt . Then x ∈ Vt
gives
(Γ,D)

µt

= E

X∼D



(Γ) 

X X ∈ Vt

= E



X∼D

(Γ)

precisely when x ∈ QVt

. As such, Lemma 19

(Γ) 

X X ∈ QVt

(Γ) 
X Q−1 X ∈ Vt
X∼D


(Γ) 
(Γ) 
(Γ,D)
= E QY Y ∈ Vt
= Q E Y Y ∈ Vt
= Qµt
.

= E



Y ∼D

(Γ,D)

By (si3), this then implies that µt

Y ∼D

∈ span{γt }.

20

(Γ,D)

(Γ,D)

At this point, we have µt
= αt γt , where αt = hµt
, γt i/kγt k22 . For any given s ∈ {1, . . . , k},
pick Q such that Qγt = γs (which exists by (si2)). Then Lemma 19 again gives
D
 E

(Γ)
X
∈
V
, γs
hµ(Γ,D)
,
γ
i
=
X
E
s
s
s
X∼D


= E hX, γs i X ∈ Vs(Γ)
X∼D

(Γ) 
= E hX, γs i X ∈ QVt
X∼D


(Γ) 
(Γ) 
(Γ,D)
= E hQY, γs i Y ∈ Vt
= E hY, γt i Y ∈ Vt
= hµt
, γt i,
Y ∼D

Y ∼D

(Γ,D)

meaning αt = α for all t. It remains to show that hµt
, γt i > 0. To this end, note that kx−γt k22 <
(Γ)
kx − γs k22 precisely when hx, γt i > hx, γs i, and so x ∈ Vt if and only if t = arg maxa hx, γa i. By
Lemma 18,
 X

k
k
1X
1
max hx, γa i ≥
hx, γa i =
x,
γa = 0,
k
k
a∈{1,...,k}
a=1

a=1

(Γ)

with equality only if the maximizer is not unique, and so we conclude that x ∈ Vt
hx, γt i > 0. As such,
D
E


(Γ,D)
(Γ) 
(Γ) 
hµt
> 0,
, γt i =
E X X ∈ Vt , γt = E hX, γt i X ∈ Vt
X∼D

only if

X∼D

as desired.

7

Proof of Theorem 7

We start with the following:
Lemma 20. Let Γ ⊆ Rd denote the standard orthoplex of dimension d, and for any given c ≥ 0,
consider the mixture Dc of equally weighted spherical Gaussians of unit entrywise variance centered
(Γ)
at the members of cΓ. Let V1 denote the Voronoi region corresponding to the first identity basis
element, and define αd : R≥0 → R by
Z ∞ Z x1
Z x1

(Γ) 
αd (c) := E X1 X ∈ V1
= 2d
···
x1 f (x; c) dxd · · · dx1 ,
(21)
X∼Dc

x1 =0

x2 =−x1

xd =−x1

where f (·; c) denotes the probability density function of Dc :
2d

f (x; c) :=

1 X
1
2
e−kx−cγt k2 /2 .
d/2
2d
(2π)
t=1

Then αd (c) ≥ αd (0) for all c ≥ 0.
See the end of this section for the proof. For context, our proof of Theorem 7 requires a bound
on αd (c) for general c ≥ 0, and so Lemma 20 allows us to pass to the easier-to-estimate quantity
αd (0). The following lemmas estimate αd (0):
√
Lemma 21. If g ∼ N (0, I) in Rd , then Ekgk∞ & log d.
21

p
Proof. When d = 1, kgk∞ has half-normal distribution, and so its expected value is 2/π. Otherwise, d ≥ 2. Since kgk∞ ≥ maxi g(i), we will estimate E maxi g(i). To this end, take z such that
P(g(1) ≥ z) = 1/d, denote j := arg maxi g(i), and condition on the event that g(j) < z, which
occurs with probability (1 − 1/d)d :
h
i
h
i 

Eg(j) = E g(j) g(j) < z · (1 − 1/d)d + E g(j) g(j) ≥ z · 1 − (1 − 1/d)d
i


1 h
≥ E g(j) g(j) < 0 · (1/4) + z · 1 − (1 − 1/d)d .
(22)
2
P
Since g(j) ≥ d1 di=1 g(i), we have
r
 X

d
h
i
2
1
E g(j) g(j) < 0 ≥ E
g(i) g(i) < 0 ∀i = −
.
d
π
i=1

With this, we may continue (22) to get
r
p
1 2
Eg(j) ≥ −
+ z · (1 − e−1 ) & z & log d,
8 π
2

where the last step follows from rearranging 1/d = P(g(1) ≥ z) ≥ e−O(z ) .
√
Lemma 22. The function αd defined by (21) satisfies αd (0) & log d.
(Γ)

Proof. It’s straightforward to verify that x ∈ Vt
j = arg max |x(i)|

precisely when

and

i∈{1,...,d}
(Γ)

and so x ∈ Vt

γt = sign(x(j)) · ej ,

implies hx, γt i = kxk∞ . Letting g ∼ N (0, I) in Rd , then
h
i
h
i
(Γ)
(Γ)
αd (0) = E hg, γ1 i g ∈ V1
= E kgk∞ g ∈ V1
= Ekgk∞ ,

where the last step follows from the fact that kΠgk∞ has the same distribution for every signed
(Γ)
permutation Π, since this implies that the random variable is independent of the event g ∈ V1 .
The result then follows from Lemma 21.
We are now ready to prove the theorem of interest:
Proof of Theorem 7. Take Γ to be the standard orthoplex of dimension d = k/2. Observe that
Theorem 5 along with a change of variables in (21) gives
(Γ,D)

µt

= αγt = σαd (1/σ)γt .

This then implies
(Γ,D)

kµt

(Γ,D)

− γt k2 = |hµt

(Γ,D)

− γt , γt i| ≥ |hµt

, γt i| − |hγt , γt i|
√
√
= σαd (1/σ) − ∆min / 2 ≥ σαd (0) − ∆min / 2,

22

(23)

where the last step follows from Lemma 20. At this point, we consider two cases. In the first case,
(23) ≥ σαd (0)/2, which implies
p
√
(Γ,D)
kµt
− γt k2 ≥ σαd (0) − ∆min / 2 ≥ σαd (0)/2 & σ log k,
where the last step follows from Lemma 22. Since this bound is independent of t, we then get
p
(Γ,D)
min kµt
− γt k2 & σ log k.
t∈{1,...,k}

In the remaining case, we have (23) < σαd (0)/2 which one may rearrange to get
p
2
∆min
σ<
· √ . ∆min / log k.
αd (0)
2
Proof of Lemma 20. We will show that αd (c) has a nonnegative
derivative, and the result will follow
P
I
(c),
where
from the mean value theorem. First, we write αd (c) = 2d
t=1 t
Z
1
2
It (c) :=
x1 e−kx−cγt k2 /2 dx.
(24)
(2π)d/2 x∈V1(Γ)
Denote γ1 = e1 and γ2 = −e1 , i.e., the first identity basis element and its negative. We claim that
Is (·) = It (·) whenever s, t ∈ {3, . . . , 2d}. This can be seen by changing variables in (24) with any
(Γ)
signed permutation that fixes γ1 (and therefore γ2 and V1 ), and that sends γs to γt . As such, we
have αd (c) = I1 (c) + I2 (c) + (2d − 2)I3 (c), where we take γ3 = e2 without loss of generality. At this
point, we factor out the x1 -dependence in the integrands of I1 (c) and I2 (c) and observe that
 Z x1
d−1
Z x1
Z x1
√
−z 2 /2
−(x22 +···+x2d )/2
dxd · · · dx2 =
e
dz
= (2π)(d−1)/2 erf d−1 (x1 / 2)
···
e
x2 =−x1

−x1

xd =−x1

to get
1
I1 (c) + I2 (c) = √
2π

Z

∞



√
2
2
x e−(x−c) /2 + e−(x+c) /2 erf d−1 (x/ 2)dx.

0

Similarly,
1
I3 (c) =
2π

Z

∞

−x21 /2

Z

x1

−(x2 −c)2 /2

x1 e
x1 =0

e

√
dx2 erf d−2 (x1 / 2)dx1 .


x2 =−x1

At this point, we apply differentiation under the integral sign to get
Z ∞ 

√
1
2
2
0
0
I1 (c) + I2 (c) = √
x (x − c)e−(x−c) /2 − (x + c)e−(x+c) /2 erf d−1 (x/ 2)dx,
2π 0
Z ∞


√
1
2
2
2
0
I3 (c) =
xe−x /2 e−(x+c) /2 − e−(x−c) /2 erf d−2 (x/ 2)dx.
2π 0
To continue, note that
√
√
d
2(d − 1) −x2 /2 d−2
erf d−1 (x/ 2) = √
e
erf
(x/ 2).
dx
2π
With this, we integrate by parts to change the expression for I30 (c):


Z ∞

√
1
1
0
0
0
−(x−c)2 /2
−(x+c)2 /2
d−1
− I1 (c) − I2 (c) + √
I3 (c) =
e
−e
erf
(x/ 2)dx .
2d − 2
2π 0
Overall, we have
Z ∞

√
1
2
2
0
0
0
0
αd (c) = I1 (c) + I2 (c) + (2d − 2)I3 (c) = √
e−(x−c) /2 − e−(x+c) /2 erf d−1 (x/ 2)dx,
2π 0
which is nonnegative since the integrand is everywhere nonnegative.
23

Acknowledgments
The authors thank Pablo Parrilo for suggesting SDPNAL+ as a solver for the k-means semidefinite
program, and David Bowie for the music that sustained our investigation of stable isogons and
the Voronoi Means Conjecture. DGM was supported by an AFOSR Young Investigator Research
Program award, NSF Grant No. DMS-1321779, and AFOSR Grant No. F4FGA05076J002. RW
was supported in part by an NSF CAREER grant and ASOFR Young Investigator Award 9550-131-0125. The views expressed in this article are those of the authors and do not reflect the official
policy or position of the United States Air Force, Department of Defense, or the U.S. Government.

References
[AAB+ 15] Martı́n Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig
Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat,
Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga,
Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit
Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin
Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large-scale machine learning on
heterogeneous systems, 2015. Software available from tensorflow.org.
[ABC+ 15] P. Awasthi, A. Bandeira, M. Charikar, R. Krishnaswamy, S. Villar, and R. Ward. Relax, no need to round: Integrality of clustering formulations. Proceedings of the 2015
Conference on Innovations in Theoretical Computer Science, pages 191–200, 2015.
[AK05]

S. Arora and R. Kannan. Learning mixtures of arbitrary gaussians. In STOC, 2005.

[AM07]

D. Achlioptas and F. McSherry. Fast computation of low-rank approximations. Journal
of ACM, 2007.

[AS12]

P. Awasthi and O. Sheffet. Improved spectral-norm bounds for clustering. In APPROXRANDOM, 2012.

[Bar95]

Alexander I. Barvinok. Problems of distance geometry and convex properties of
quadratic maps. Discrete & Computational Geometry, 13(2):189–202, 1995.

[BW13]

H. Broome and S. Waldron. On the construction of highly symmetric tight frames and
complex polytopes. Linear Algebra and its Applications, 439(12):4135 – 4151, 2013.

[Das99]

S. Dasgupta. Learning mixtures of gaussians. In FOCS, 1999.

[DS07]

S. Dasgupta and L. Schulman. A probabilistic analysis of em for mixtures of separated,
spherical gaussians. J. Mach. Learn. Res., 8:203–226, May 2007.

[GV14]

O. Guédon and R. Vershynin.
grothendieck’s inequality, 2014.

Community detection in sparse networks via

[IMPV15] T. Iguchi, D. G. Mixon, J. Peterson, and S. Villar. Probably certifiably correct k-means
clustering, 2015.
24

[KK10]

A. Kumar and R. Kannan. Clustering with spectral norm and the k-means algorithm.
In FOCS, 2010.

[LC10]

Yann LeCun and Corinna Cortes. Mnist handwritten digit database. AT&T Labs
[Online]. Available: http://yann. lecun. com/exdb/mnist, 2010.

[MV10]

A. Moitra and G. Valiant. Settling the polynomial learnability of mixtures of gaussians. In Proceedings of the 51st Annual IEEE Symposium on Foundations of Computer
Science, 2010.

[NW15]

A. Nellore and R. Ward. Recovery guarantees for exemplar-based clustering. Information and Computation, 245(C):165–180, 2015.

[Pat98]

Gábor Pataki. On the rank of extreme matrices in semidefinite programs and the
multiplicity of optimal eigenvalues. Mathematics of operations research, 23(2):339–358,
1998.

[PW07]

J. Peng and Y. Wei. Approximating k-means-type clustering via semidefinite programming. SIAM J. on Optimization, 18(1):186–205, February 2007.

[Sha82]

Alexander Shapiro. Rank-reducibility of a symmetric matrix and sampling theory of
minimum trace factor analysis. Psychometrika, 47(2):187–199, 1982.

[Ver12]

R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. Compressed Sensing, Theory and Applications, Y. Eldar, G. Kutyniok (eds.), Cambridge
University Press, 2012.

[Vil16]

Soledad Villar. Implementation of the k-means semidefinite program, 2016. Software
available from github.com.

[VW04]

S. Vempala and G. Wang. A spectral algorithm for learning mixture models. Journal
of Computer and System Sciences, 68(2):841–860, 2004.

[YST15]

Liuqin Yang, Defeng Sun, and Kim-Chuan Toh. Sdpnal+: a majorized semismooth
newton-cg augmented lagrangian method for semidefinite programming with nonnegative constraints. Mathematical Programming Computation, 7(3):331–366, 2015.

25

