Probably certifiably correct k-means clustering
Takayuki Iguchi1 , Dustin G. Mixon1 , Jesse Peterson1 , and Soledad Villar2

arXiv:1509.07983v2 [] 23 Apr 2016

1

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

Abstract
Recently, Bandeira [5] introduced a new type of algorithm (the so-called probably certifiably
correct algorithm) that combines fast solvers with the optimality certificates provided by convex
relaxations. In this paper, we devise such an algorithm for the problem of k-means clustering.
First, we prove that Peng and Wei’s semidefinite relaxation of k-means [20] is tight with high
probability under a distribution of planted clusters called the stochastic ball model. Our proof
follows from a new dual certificate for integral solutions of this semidefinite program. Next, we
show how to test the optimality of a proposed k-means solution using this dual certificate in
quasilinear time. Finally, we analyze a version of spectral clustering from Peng and Wei [20]
that is designed to solve k-means in the case of two clusters. In particular, we show that this
quasilinear-time method typically recovers planted clusters under the stochastic ball model.

1

Introduction

Clustering is a central problem in unsupervised machine learning. It consists of partitioning a given
m into k subsequences such that some dissimilarity function
finite sequence {xi }N
i=1 of points in R
is minimized. Usually, this function is chosen ad hoc with an application in mind. A particularly
common choice is the k-means objective:
minimize

k X
X
t=1 i∈At

xi −

1 X
xj
|At |
j∈At

2

(1)
2

subject to A1 t · · · t Ak = {1, . . . , N }
Problem (1) is NP-hard in general [13]. A popular heuristic for solving k-means is Lloyd’s algorithm,
also known as the k-means algorithm [15]. This algorithm alternates between calculating centroids
of proto-clusters and reassigning points according to the nearest centroid. In general, Lloyd’s
algorithm (and its variants [3, 19]) may converge to local minima of the k-means objective (e.g.,
see section 5 of [4]). Furthermore, the output of Lloyd’s algorithm does not indicate how far it is
from optimal. Instead, we seek a new sort of algorithm, recently introduced by Bandeira [5]:
Definition 1. Let P be an optimization problem that depends on some input, and let D denote a
probability distribution over possible inputs. Then a probably certifiably correct (PCC) algorithm for (P, D) is an algorithm that on input D ∼ D produces a global optimizer of P with high
probability, and furthermore produces a certificate of having done so.
Most non-convex optimization methods fail to produce a certificate of global optimality. However, if a non-convex problem enjoys a convex relaxation, then solving the dual of this relaxation
1

will produce a certificate of (approximate) optimality. Along these lines, the k-means problem enjoys a semidefinite relaxation. To see this, let 1A denote the indicator function of A ⊆ {1, . . . , N },
and define the N × N matrix D by Dij := kxi − xj k22 . Then taking
k
X
1
X :=
1A 1> ,
|At | t At

(2)

t=1

the k-means objective (1) may be expressed as 12 Tr(DX). Since X satisfies several convex constraints, we may relax the region of optimization to produce a convex program, namely, the Peng–
Wei semidefinite relaxation [20] (cf. [6]):
minimize

Tr(DX)

(3)

subject to

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

Here, X ≥ 0 means that each entry of X is nonnegative, whereas X  0 means that X is symmetric
and positive semidefinite.
Recently, it was shown that under a certain random data model, this convex relaxation is
tight with high probability [4], that is, every solution to the relaxed problem (3) has the form (2),
thereby identifying an optimal clustering. As such, in this high-probability event, one may solve
the dual program to produce a certificate of optimality. However, semidefinite programming (SDP)
solvers are notoriously slow. For example, running MATLAB’s built-in implementation of Lloyd’s
algorithm on 64 points in R6 will take about 0.001 seconds, whereas a CVX implementation [11]
of the dual of (3) for the same data takes about 20 seconds. Also, Lloyd’s algorithm scales much
better than SDP solvers, and so one should expect this runtime disparity to only increase with larger
datasets. Overall, while the SDP relaxation theoretically produces a certificate in polynomial time
(e.g., by an interior-point method [18]), it is far too slow to wait for in practice.
As a fast alternative, Bandeira [5] recently devised a general technique to certify global optimality. This technique leverages several components simultaneously:
(i) A fast non-convex solver that produces the optimal solution with high probability (under
some reasonable probability distribution of problem instances).
(ii) A convex relaxation that is tight with high probability (under the same distribution).
(iii) A fast method of computing a certificate of global optimality for the output of the non-convex
solver in (i) by exploiting convex duality with the relaxation in (ii).
In the context of k-means, one might expect Lloyd’s algorithm and the Peng–Wei SDP to be suitable
choices for (i) and (ii), respectively. For (iii), one might adapt Bandeira’s original method in [5]
based on complementary slackness (see Figure 1 for an illustration). In this paper, we provide a
theoretical basis for each of these components in the context of k-means.

1.1

Technical background and overview

The first two components of a probably certifiably correct algorithm require non-convex and convex
solvers that perform well under some “reasonable” distribution of problem instances. In the context
of geometric clustering, it has become popular recently to consider a particular model of data called
the stochastic ball model, introduced in [17]:

2

y

y

x

x

Figure 1: (left) Depiction of complementary slackness. The horizontal axis represents a vector space in which
we consider a cone program (e.g., a linear or semidefinite program), and the feasibility region of this program is
highlighted in green. The dual program concerns another vector space, which we represent with the vertical axis
and feasibility region highlighted in red. The downward-sloping line represents all pairs of points (x, y) that satisfy
complementary slackness. Recall that when strong duality is satisfied, we have that x is primal-optimal and y is
dual-optimal if and only if x is primal feasible, y is dual feasible, and (x, y) satisfy complementary slackness. As
such, the intersection between the blue Cartesian product and the complementary slackness line represents all pairs
of optimizers. (right) Bandeira’s probably certifiably correct technique [5]. Given a purported primal-optimizer x,
we first check that x is primal-feasible. Next, we select y such that (x, y) satisfies complementary slackness. Finally,
we check that y is dual-feasible. By complementary slackness, y is then a dual certificate of x’s optimality in the
primal program, which can be verified by comparing their values (a la strong duality).

Definition 2 ((D, γ, n)-stochastic ball model). Let {γa }ka=1 be ball centers in Rm . For each a,
draw i.i.d. vectors {ra,i }ni=1 from some rotation-invariant distribution D supported on the unit ball.
The points from cluster a are then taken to be xa,i := ra,i + γa . We denote ∆ := mina6=b kγa − γb k2 .
Table 1 summarizes the state of the art for recovery guarantees under the stochastic ball model.
In [17], it was shown that an LP relaxation of k-medians will, with high probability, recover clusters
drawn from the stochastic ball model provided the smallest distance between ball centers is ∆ ≥
3.75. Note that exact recovery only makes sense for ∆ > 2 (i.e., when the balls are disjoint). Once
∆ > 4, any two points within a particular cluster are closer to each other than any two points from
different clusters, and so in this regime, cluster recovery follows from a simple distance thresholding.
For the k-means problem, Awasthi et al. [4] studies
the Peng–Wei semidefinite relaxation and
√
√
demonstrates exact recovery in the regime ∆ > 2 2(1 + 1/ m), where m is the dimension of the
Euclidean space.
As indicated in Table 1, we also study the Peng–Wei SDP, but our guarantee is different from [4].
In particular, we demonstrate tightness in the regime ∆ > 2 + k 2 /m, which is near-optimal for
large m. The source of this improvement is a different choice of dual certificate, which leads to the
following result (see Section 2 for details):
Theorem 3. Take X of the form (2), and let PΛ⊥ denote the orthogonal projection onto the
orthogonal complement of the span of {1At }kt=1 . Then there exists an explicit matrix Z = Z(D, X)
and scalar z = z(D, X) such that X is a solution to the semidefinite relaxation (3) if
PΛ⊥ ZPΛ⊥  zPΛ⊥ .

(4)

To prove that ∆ > 2 + k 2 /m suffices for the SDP to recover the planted clustering under the
stochastic ball model, we estimate the left- and right-hand sides of (4) with the help of standard
techniques from random matrix theory and concentration of measure; see Appendix B for the
3

Table 1:

Method

Sufficient Condition

Optimal?

Reference

Thresholding

∆>4

Yes

(simple exercise)

k-medians LP

∆≥4
∆ ≥ 3.75
∆>2

No
No
Yes

Theorem 2 in [9]
Theorem 1 in [17]
Theorem 1 in [4]

k-means LP

∆>4

Yes

Theorem 9 in [4]

k-means SDP

√
∆ > 2 2(1 + 1/ m)
∆ > 2 + k 2 /m

No
No

Theorem 3 in [4]
Theorem 9

Spectral k-means

∆ > ∆? , k = 2

Yes

Theorem 14

√

Summary of cluster recovery guarantees under the stochastic ball model. The second column reports

sufficient separation between ball centers in order for the corresponding method to provably give exact recovery
with high probability. The third column reports whether the sufficient condition on ∆ cannot be improved. Here,
∆? = ∆? (D, k) denotes the smallest value for which ∆ > ∆? implies that minimizing the k-means objective recovers
planted clusters under the (D, γ, n)-stochastic ball model with probability 1 − e−ΩD,γ (n) .

(rather technical) details. While this is an improvement over the condition from [4] in the large-m
regime, there are other regimes (e.g., k = m) for which their condition is much better, leaving open
the question of what the optimal bound is. Conjecture 4 in [4] suggests that ∆ > 2 suffices for the
k-means SDP to recover planted clusters under the stochastic ball model, but as we illustrate in
Section 2.3, this conjecture is surprisingly false.
Having accomplished component (ii) in Bandeira’s PCC technique, we tackle component (iii)
next. For this, consider the matrix
A :=

z >
11 + PΛ⊥ ZPΛ⊥ ,
N

(5)

where z and Z come from Theorem 3. Since the all-ones vector 1 lies in the span of {1At }kt=1 , we
have that 1 spans the unique leading eigenspace of A precisely when PΛ⊥ ZPΛ⊥ ≺ zPΛ⊥ , which in
turn implies that X is a k-means optimal clustering by Theorem 3. As such, component (iii) can
be accomplished by solving the following fundamental problem from linear algebra:
Problem 4. Given a symmetric matrix A ∈ Rn×n and an eigenvector v of A, determine whether the
span of v is the unique leading eigenspace, that is, the corresponding eigenvalue λ has multiplicity 1
and satisfies |λ| > |λ0 | for every other eigenvalue λ0 of A.
Interestingly, this same problem appeared in Bandeira’s original PCC theory [5], but it was left
unresolved. In this paper, we fill this gap by developing a so-called power iteration detector, which
applies the power iteration to a random initialization on the unit sphere. Due to the randomness,
the power iteration produces a test statistic that allows us to infer whether (A, v) satisfies the
desired leading eigenspace condition. In Section 3, we pose this as a hypothesis test, and we
estimate the associated error probabilities. In addition, we show how to leverage the structure of
A defined by (5) and Theorem 4 to compute the matrix–vector multiplication Ax for any given
x in only O(kmN ) operations, thereby allowing the test statistic to be computed in linear time
(up to the spectral gap of A and the desired confidence for the hypothesis test). See Figure 2 for
an illustration of the runtime of our method. Overall, the power iteration detector will deliver a
highly confident inference on whether (A, v) satisfies the leading eigenspace condition, which in
4

2

10

1

10

k-means SDP

spectral
k-means

0

10

O(N )

power
iteration
detector

−1

10

−2

10

k-means++

−3

10

−4

10

0

10

Figure 2:

1

10

2

10

3

10

4

10

5

10

6

10

Take two unit balls in R6 at distance 2.3 apart. For each N ∈ {23 , 24 , . . . , 216 }, we perform 300 trials

of the following experiment: Draw N/2 points uniformly at random from each ball, and then compute four different
functions: (a) MATLAB’s built-in implementation of k-means++, (b) a CVX implementation [11] of the k-means
SDP (3), (c) the power iteration detector (Algorithm 1) with A given by (5), and (d) spectral k-means clustering
(Algorithm 2). For each trial, we recorded the runtime in seconds. Above, we plot the average runtime along with
error bars for standard deviation. For the record, the power iteration detector failed to certify optimality (i.e., reject
H0 in (14)) in at most 3% of the trials with N ≤ 27 , but rejected H0 in every trial otherwise; similarly, spectral
k-means failed to recover the planted clusters in two of the trials with N = 23 . Our implementation of the k-means
SDP was too slow to perform trials with N ≥ 27 in a reasonable amount of time, so we only recorded runtimes for
N ∈ {23 , 24 , 25 , 26 }. As the plot illustrates, the other algorithms ran in quasilinear time, as expected.

turn certifies the optimality of X up to the prescribed confidence level. Of course, one may remove
the need for a confidence level by opting for deterministic spectral methods, but we have no idea
how to accomplish this in linear or even near-linear time.
Now that we have discussed components (ii) and (iii) in Bandeira’s PCC technique, we conclude
by discussing component (i). While we presume that there exists a fast initialization of Lloyd’s
algorithm that performs well under the stochastic ball model, we leave this investigation for future
research. Instead, Section 4 considers a spectral method introduced by Peng and Wei in [20]. We
show that when k = 2, this method performs as well as the optimizer of the original k-means problem
under the stochastic ball model. Figure 2 illustrates the quasilinear runtime of this approach.

5

1.2

Outline

In this paper, we provide a theoretical analysis of probably certifiably correct k-means clustering,
and we do so by developing components (i), (ii) and (iii) of Bandeira’s general technique. First,
we investigate (ii) in Section 2 by analyzing the tightness of the Peng–Wei SDP. In particular, we
choose a different dual certificate from the one used in [4], and our choice demonstrates tightness
in the SDP for clusters that are near-optimally close. Section 3 then addresses (iii) by providing
a fast method of computing this dual certificate given the optimal k-means partition. In fact,
a subroutine of our method (the so-called power iteration detector) resolves a gap in Bandeira’s
original PCC theory [5], and as such, we expect this to be leveraged in future PCC algorithms. We
conclude in Section 4 with some theoretical guarantees for (i). Here, we focus on the case k = 2,
and we show that a slight modification of the spectral clustering–based method in [20] manages to
recover the optimal k-means partition with high probability under the stochastic ball model. We
conclude in Section 5 with a discussion of various open problems.

2

A typically tight relaxation of k-means

This section establishes that the Peng–Wei semidefinite relaxation (3) of the k-means problem (1)
is typically tight under the stochastic ball model. First, we find a deterministic condition on the set
of points under which the relaxation finds the k-means-optimal solution. Later, we discuss when
this deterministic condition is satisfied with high probability under the stochastic ball model.

2.1

The dual program

The following is the dual program of (3):
minimize

kz +

N
X

αi

(6)

i=1

subject to

Q := zI +

N
X
i=1

N

N

XX
1
1
>
αi · (ei 1> + 1e>
βij · (ei e>
i )−
j + ej ei ) + D  0
2
2
i=1 j=i

β≥0
For notational simplicity, from this point forward, we organize indices according to clusters.
For example, 1a shall denote the indicator function of the ath cluster. Also, we shuffle the rows
and columns of X and D into blocks that correspond to clusters; for example, the (i, j)th entry of
(a,b)
the (a, b)th block of D is given by Dij . We also index α in terms of clusters; for example, the ith
entry of the ath block of α is denoted αa,i . For β, we identify
β :=

N X
N
X
i=1 j=i

1
>
βij · (ei e>
j + ej ei ).
2

Indeed, when i ≤ j, the (i, j)th entry of β is βij . We also consider β as having its rows and columns
(a,b)
shuffled according to clusters, so that the (i, j)th entry of the (a, b)th block is βij .
With this notation, the following proposition characterizes all possible dual certificates of (3):
P
Proposition 5 (Theorem 4 in [12], cf. [4]). Take X := ka=1 n1a 1a 1>
a , where na denotes the number
of points in cluster t. The following are equivalent:
6

(a) X is a solution to the semidefinite relaxation (3).
(b) Every solution to the dual program (6) satisfies
Q(a,a) 1 = 0,

β (a,a) = 0

∀a ∈ {1, . . . , k}.

(c) Every solution to the dual program (6) satisfies
αa,r = −

1
1
2
z + 2 1> D(a,a) 1 − e>
D(a,a) 1
na
na
na r

∀a ∈ {1, . . . , k}, r ∈ a.

The following subsection will leverage this result to identify a condition on D that implies that
the SDP (3) relaxation is tight.

2.2

Selecting a dual certificate

The goal is to certify when the SDP relaxation is tight. In this event, Proposition 5 characterizes
acceptable dual certificates (z, α, β), but this information fails to uniquely determine a certificate.
In this subsection, we will motivate the application of additional constraints on dual certificates so
as to identify certifiable instances.
We start by reviewing the characterization of dual certificates (z, α, β) provided in Proposition 5.
In particular, α is completely determined by z, and so z and β are the only remaining free variables.
Indeed, for every a, b ∈ {1, . . . , k}, we have
X
k X

(a,b)
1
>
>
αt,i · (et,i 1 + 1et,i )
2
t=1 i∈t
X
X
1
1
=
αa,i · ei 1> +
αb,j · 1e>
2
2 j
i∈a
j∈b



X 1
1
2 > (a,a) 1 >
1 1
> (a,a)
+
1 D
1 − ei D
1 ei 1
=−
z+
2 na nb
n2a
na
2
i∈a


X 1
1
2
(b,b)
+
1 1e>
,
1> D(b,b) 1 − e>
j D
2
nb
2 j
nb
j∈b

and so since
Q = zI +

k X
X
t=1 i∈t

1
1
αt,i · (et,i 1> + 1e>
t,i ) − β + D,
2
2

we may write Q = z(I − E) + M − B, where


1 1
1
(a,b)
E
:=
+
11>
2 na nb

X 1
2 > (a,a) 1 >
> (a,a)
M (a,b) := D(a,b) +
1
D
1
−
e
D
1
ei 1
n2a
na i
2
i∈a

X 1
2 > (b,b) 1 >
> (b,b)
1 1ej
+
1 D
1 − ej D
nb
2
n2b
j∈b

1
B (a,b) = β (a,b)
2
7

(7)

(8)

for every a, b ∈ {1, . . . , k}. The following is one way to formulate our task: Given D and a clustering
X (which in turn determines E and M ), determine whether there exist feasible z and B such that
Q  0; here, feasibility only requires B to be symmetric with nonnegative entries and B (a,a) = 0 for
every a ∈ {1, . . . , k}. We opt for a slightly more modest goal: Find z = z(D, X) and B = B(D, X)
such that Q  0 for a large family of D’s.
Before determining z and B, we first analyze E:
Lemma 6. Let E be the matrix defined by (7). Then rank(E) ∈ {1, 2}. The eigenvalue of largest
magnitude is λ ≥ k, and when rank(E) = 2, the other nonzero eigenvalue of E is negative. The
eigenvectors corresponding to nonzero eigenvalues lie in the span of {1a }ka=1 .
Proof. Writing



 k
X
>
k X
k
k
X
1 1
1
1
1 X 1
1
>
>
E=
+
1a 1b =
1a 1 + 1
1b ,
2 na nb
2
na
2
nb
a=1 b=1

a=1

b=1

we see that rank(E) ∈ {1, 2}, and it is easy to calculate 1> E1 = N k and Tr(E) = k. Observe that
λ = sup x> Ex ≥
x∈RN
kxk2 =1

1 >
1 E1 = k,
N

and combining with rank(E) ≤ 2 and Tr(E) = k then implies that the other nonzero eigenvalue (if
there is one) is negative. Finally, any eigenvector of E with a nonzero eigenvalue necessarily lies in
the column space of E, which is a subspace of span{1a }ka=1 by the definition of E.
When finding z and B such that Q = z(I − E) + M − B  0, it will be useful that I − E has
only one negative eigenvalue to correct. Let v denote the corresponding eigenvector. Then we will
pick B so that v is also an eigenvector of M − B. Since we want Q  0 for as many instances
of D as possible, we will then pick z as large as possible, thereby sending v to the nullspace of
Q. Unfortunately, the authors found that this constraint fails to uniquely determine B in general.
Instead, we impose a stronger constraint:
Q1a = 0

∀a ∈ {1, . . . , k}.

(This constraint implies Qv = 0 by Lemma 6.) To see the implications of this constraint, note that
we already necessarily have




1 >
(a,a)
(a,a)
(a,a)
(Q1a )a = (z(I − E) + M − B)1a = z(I − E
)1 + M
1−B
1 = z 1 − 11 1 = 0,
na
a
and so it remains to impose


0 = (Q1b )a = (z(I − E) + M − B)1b
= −zE

(a,b)

1+M

(a,b)

1−B

a

(a,b)

1 = −z

na + nb
1 + M (a,b) 1 − B (a,b) 1.
2na

In order for there to exist a vector B (a,b) 1 ≥ 0 that satisfies (9), z must satisfy
z

na + nb
≤ min(M (a,b) 1),
2na
8

(9)

and since z is independent of (a, b), we conclude that
z≤

2na
min(M (a,b) 1).
a,b∈{1,...,k} na + nb

(10)

min
a6=b

Again, in order to ensure z(I − E) + M − B  0 for as many instances of D as possible, we intend
to choose z as large as possible. Luckily, there is a choice of B which satisfies (9) for every (a, b),
even when z satisfies equality in (10). Indeed, we define
u(a,b) := M (a,b) 1 − z

na + nb
1,
2na

ρ(a,b) := u>
(a,b) 1,

B (a,b) :=

1
ρ(b,a)

u(a,b) u>
(b,a)

(11)

for every a, b ∈ {1, . . . , k} with a 6= b. Then by design, B immediately satisfies (9). Also, note
that ρ(a,b) = ρ(b,a) , and so B (b,a) = (B (a,b) )> , meaning B is symmetric. Finally, we necessarily
have u(a,b) ≥ 0 (and thus ρ(a,b) ≥ 0) by (10), and we implicitly require ρ(a,b) > 0 for division to be
permissible. As such, we also have B (a,b) ≥ 0, as desired.
Now that we have selected z and B, it remains to check that Q  0. By construction, we
already have Λ := span{1a }ka=1 in the nullspace of Q, and so it suffices to ensure


0  PΛ⊥ QPΛ⊥ = PΛ⊥ z(I − E) + M − B PΛ⊥ = zPΛ⊥ + PΛ⊥ (M − B)PΛ⊥ .
Here, PΛ⊥ denotes the orthogonal projection onto the orthogonal complement of Λ. Rearranging
then gives the following result:
P
Theorem 7. Take X := kt=1 n1t 1t 1>
t , where nt denotes the number of points in cluster t. Consider
M defined by (8), pick z so as to satisfy equality in (10), take B defined by (11), and let Λ denote
the span of {1t }kt=1 . Then X is a solution to the semidefinite relaxation (3) if
PΛ⊥ (B − M )PΛ⊥  zPΛ⊥ .

(12)

The next subsection leverages this sufficient condition to establish that the Peng–Wei SDP (3)
is typically tight under the stochastic ball model.

2.3

Integrality of the relaxation under the stochastic ball model

We first note that our sufficient condition (12) is implied by
kPΛ⊥ M PΛ⊥ k2→2 + kPΛ⊥ BPΛ⊥ k2→2 ≤ z.
By further analyzing the left-hand side above (see Appendix A), we arrive at the following corollary:
P
Corollary 8. Take X := kt=1 n1t 1t 1>
t , where nt denotes the number of points in cluster t. Let Ψ
denote the m × N matrix whose (a, i)th column is xa,i − ca , where
ca :=

1 X
xa,i
na
i∈a

denotes the empirical center of cluster a. Consider M defined by (8), pick z so as to satisfy equality
in (10), and take ρ(a,b) defined by (11). Then X is a solution to the semidefinite relaxation (3) if
2kΨk22→2

k
k
X
X
kP1⊥ M (a,b) 1k2 kP1⊥ M (b,a) 1k2
+
≤ z.
ρ(a,b)
a=1 b=a+1

9

In Appendix B, we leverage the stochastic ball model to bound each term in Corollary 8, and
in doing so, we identify a regime in which the data points typically satisfy the sufficient condition
given in Corollary 8:
Theorem 9. The k-means semidefinite relaxation (3) recovers the planted clusters in the (D, γ, n)stochastic ball model with probability 1 − e−ΩD,γ (n) provided ∆ > 2 + k 2 /m.
We note that Theorem 9 is an improvement to the main result of the authors’ preprint [12].
When k = o(m1/2 ), Theorem 9 is near-optimal, and in this sense, it’s a significant improvement
over the sufficient condition


√
1
∆>2 2 1+ √
(13)
m
given in [4]. However, there are regimes (e.g., k = m) for which (13) is much better, leaving open
the question of what the optimal bound is. Conjecture 4 in [4] suggests that ∆ > 2 suffices for
the k-means SDP to recover planted clusters under the stochastic ball model, but as we illustrate
below, this conjecture is surprisingly false.
Consider the special case where m = 1, D is uniform on {±1}, and k = 2. Centering the
two balls on ±∆/2, then all of the points land in {±∆/2 ± 1}. The k-means-optimal clustering
will partition the real line into two semi-infinite intervals, and so there are three possible ways of
clustering these points. Suppose exactly N/4 of the points land in each of the 4 positions. Then by
symmetry, there are only two ways to cluster: either we select the planted clusters, or we make the
left-most location its own cluster. Interestingly, a little algebra
reveals that this second alternative
√
is strictly better in the k-means sense provided ∆ < 1 + 3 ≈ 2.7320. Also, in this regime, then as
N gets large, the proportion of points in each position will be so close to 1/4 (with high probability)
that this clustering will beat the planted clusters.
√
Overall, when m = 1 and k = 2, then ∆ ≥ 1 + 3 is necessary for minimizing the k-means
objective to recover planted clusters for an arbitrary D. As a relaxation, the k-means SDP recovers
planted clusters only if minimizing the k-means objective does so as well, and so it inherits this necessary condition, thereby disproving Conjecture 4 in [4]. Furthermore, as Figure 3(left) illustrates,
a similar counterexample is available in higher dimensions.
To study when the SDP recovers the clusters, let’s continue with the case where m = 1 and
k = 2. We
√ know that minimizing k-means will recover the clusters with high probability provided
∆ > 1 + 3. However, Theorem 9 only guarantees that the SDP recovers the clusters when ∆ > 6;
in fact, (13) is slightly better here, giving that ∆ ≥ 5.6569 suffices. To shed light on the disparity,
Figure 3(center) illustrates the performance of the SDP for different values of ∆. Observe that the
SDP is often tight when ∆ is close to 2, but it doesn’t reliably recover the planted clusters until
∆ > 4. We suspect that ∆ = 4 is a phase transition for cluster recovery in this case.
Qualitatively, the biggest difference between Theorem 9 and (13) is the dependence on k that
Theorem 9 exhibits. Figure 3(right) illustrates that this comes from (12), meaning that one would
need to use a completely different dual certificate in order to remove this dependence.

3

A fast test for k-means optimality

In this section, we leverage the certificate (12) to test the optimality of a candidate k-means
solution. We first show how to solve a more general problem from linear algebra, and then we
apply our solution to devise a fast test for k-means optimality (as well as fast test for a related
PCC algorithm).

10

1.0020

1

20

1.0015

0.8

16

1.0010

0.6

12

1.0005

0.4

1.0000

0.2

0.9995

0
0

Figure 3:

0.05

0.1

0.15

0.2

8
4

2

3

4

5

2

2.2

2.4 2.6

2.8 3.0

(left) Take two unit disks in R2 with centers on the x-axis at distance 2.08 apart. Let x0 denote

the smallest possible x-coordinate in the disk on the right. For each disk, draw N/2 = 50, 000 points uniformly at
random from the perimeter. Given θ, cluster the points according to whether the x-coordinate is smaller than x0 + θ.
When θ = 0, this clustering gives the planted clusters, and the k-means objective (divided by N ) is 1. We plot this
normalized k-means objective for θ ∈ [0, 0.2]. Since N is large, this curve is very close to its expected shape, and we
see that there are clusters whose k-means value is smaller than that of the planted clustering. (center) Take two
intervals of width 2 in R, and let ∆ denote the distance between the midpoints of these intervals. For each interval,
draw 10 points at random from its endpoints, and then run the k-means SDP. For each ∆ = 2 : 0.1 : 5, after running
2, 000 trials of this experiment, we plot the proportion of trials for which the SDP relaxation was tight (dashed line)
and the proportion of trials for which the SDP recovered the planted clusters (solid line). In this case, cluster recovery
appears to exhibit a phase transition at ∆ = 4. (right) For each ∆ = 2 : 0.1 : 3 and k = 2 : 2 : 20, consider the
unit balls in R20 centered at { √∆2 ei }ki=1 , where ei denotes the ith identity basis element. Draw 100 points uniformly
from each ball, and test if the resulting data points satisfy (12). After performing 10 trials of this experiment for
each (∆, k), we shade the corresponding pixel according to the proportion of successful trials (white means every trial
satisfied (12)). This plot indicates that our certificate (12) is to blame for Theorem 9’s dependence on k.

3.1

Leading eigenvector hypothesis test

This subsection concerns Problem 4. To solve this problem, one might be inclined to apply the
power method:
Proposition 10 (Theorem 8.2.1 in [10]). Let A ∈ Rn×n be a symmetric matrix with eigenvalues
{λi }ni=1 (counting multiplicities) satisfying
|λ1 | > |λ2 | ≥ · · · ≥ |λn |,
and with corresponding orthonormal eigenvectors {vi }ni=1 . Pick a unit-norm vector q0 ∈ Rn and
consider the power iteration qj+1 := Aqj /kAqj k2 . If q0 is not orthogonal to v1 , then
(v1> qj )2

≥1−



(v1> q0 )−2

−1

 λ 2j
2

λ1

.

Notice that the above convergence guarantee depends on the quality of the initialization q0 . To
use this guarantee, draw q0 at random from the unit sphere so that q0 is not orthogonal to v1 almost
surely; one might then analyze the statistics of v1> q0 to produce statistics on the time required for
convergence. The power method is typically used to find a leading eigenvector, but for our problem,
we already have access to an eigenvector v, and we are tasked with determining whether v is the
11

Algorithm 1: Power iteration detector
Input: Symmetric matrix A ∈ Rn×n , unit eigenvector v ∈ Rn , tolerance  > 0
Output: Decision of whether to accept H0 or to reject H0 and accept H1 as given in (14)
λ ← v > Av
Draw q uniformly at random from the unit sphere in Rn
while no decision has been made do
if |q > Aq| > |λ| then
Print accept H0
else if (v > q)2 ≥ 1 −  then
Print reject H0 and accept H1
end
q ← Aq/kAqk2
end

unique leading eigenvector. Intuitively, if you run the power method from a random initialization
and it happens to converge to v, then this would have been a remarkable coincidence if v were not
the unique leading eigenvector. Since we will only run finitely many iterations, how do we decide
when we are sufficiently confident? The remainder of this subsection answers this question.
Given a symmetric matrix A ∈ Rn×n and a unit eigenvector v of A, consider the hypotheses
H0 :
H1 :

span(v) is not the unique leading eigenspace of A,
span(v) is the unique leading eigenspace of A.

(14)

To test these hypotheses, pick a tolerance  > 0 and run the power iteration detector (see Algorithm 1). This detector terminates either by accepting H0 or by rejecting H0 and accepting H1 . We
say the detector fails to reject H0 if it either accepts H0 or fails to terminate. Before analyzing
this detector, we consider the following definition:
Definition 11. Given a symmetric matrix A ∈ Rn×n and unit eigenvector v of A, put λ = v > Av,
and let λ1 denote a leading eigenvalue of A (i.e., |λ1 | = kAk2→2 ). We say (A, v) is degenerate if
(a) the eigenvalue λ of A has multiplicity ≥ 2,
(b) −λ is an eigenvalue of A, or
(c) −λ1 is an eigenvalue of A.
Theorem 12. Consider the power iteration detector (Algorithm 1), let qj denote q at the jth
2
iteration (with q0 being the initialization), and let π denote the probability that (e>
1 q0 ) < .
(i) (A, v) is degenerate only if H0 holds. If (A, v) is non-degenerate, then the power iteration
detector terminates in finite time with probability 1.
(ii) The power iteration detector incurs the following error rates:



Pr reject H0 and accept H1 H0 ≤ π ,
Pr fail to reject H0
(iii) If H1 holds, then
n
o
3 log(1/)
min j : (v > qj )2 > 1 −  ≤
+1
2 log(λ1 /λ2 )
with probability ≥ 1 − π .
12

H1



= 0.

Proof. Denote the eigenvalues of A by {λi }ni=1 (counting multiplicities), ordered in such a way that
|λ1 | ≥ · · · ≥ |λn |, and consider the corresponding orthonormal eigenvectors {vi }ni=1 , where v = vp
for some p.
For (i), first note that H1 implies that (A, v) is non-degenerate, and so the contrapositive
gives the first claim. Next, suppose (A, v) is non-degenerate. If H1 holds, then (v > qj )2 → 1 by
Proposition 10 provided q0 is not orthogonal to v, and so the power iteration detector terminates
with probability 1. Otherwise, H0 holds, and so the non-degeneracy of (A, v) implies that the
eigenspace corresponding to λ1 is the unique leading eigenspace of A, and furthermore, |λ1 | > |λ|.
Following the proof of Theorem 8.2.1 in [10], we also have
Pn
(vi> qj )2 λ2j+1
q0> A2j+1 q0
>
i
= Pi=1
qj Aqj = > 2j
.
n
> q )2 λ2j
q0 A q0
(v
i=1 i j
i
Putting r := min{i : |λi | < |λ1 |}, then
Pn
>
2 2j
>
i=1 (vi qj ) λi (λi − λ1 )
|qj Aqj − λ1 | =
Pn
>
2 2j
i=1 (vi qj ) λi
 
 2j

n
|λ1 − λn | X > 2 λi 2j
λr
1 − kPλ1 q0 k22
≤
(vi qj )
≤ |λ1 − λn |
,
2
2
λ1
λ1
kPλ1 q0 k2
kPλ1 q0 k2
i=r

where Pλ1 denotes the orthogonal projection onto the eigenspace corresponding to λ1 . As such,
|qj> Aqj | → |λ1 | > |λ| provided Pλ1 q0 6= 0, and so the power iteration detector terminates with
probability 1.
For (ii), we first consider the case of a false positive. Taking v = vp for p 6= 1, note that
(v > qj )2 > 1 −  implies
 > 1 − (v > qj )2 = kqj k22 − (vp> qj )2 =

n
X
(vi> qj )2 ≥ (v1> qj )2 .
i=1
i6=p

Also, since kAxk2 ≤ |λ1 |kxk2 for all x ∈ Rn , we have that (v1> qj )2 monotonically increases with j:

2
(v1> qj )2
(λ1 v1> qj )2
>
2
> Aqj
≥
= (v1> qj )2 .
(v1 qj+1 ) = v1
=
kAqj k2
kqj k2
kAqj k22
As such,  > (v1> qj )2 ≥ (v1> q0 )2 . Overall, when H0 holds, the power iteration detector rejects
H0 only if q0 is initialized poorly, i.e., (v1> q0 )2 < , which occurs with probability π (since q0
has a rotation-invariant probability distribution). For the false negative error rate, note that
Proposition 10 gives that H1 implies convergence (v > qj )2 → 1 provided q0 is not orthogonal to v,
i.e., with probability 1.
For (iii), we want j such that (v > qj )2 > 1 − . By Proposition 10, it suffices to have

 λ 2j
2
>
−2
< .
(v1 q0 ) − 1
λ1
In the event that (v1> q0 )2 ≥  (which has probability 1 − π ), it further suffices to have
 2j
−2 λ2

< .
λ1
Taking logs and rearranging then gives the result.
13

To estimate  and π , first note that q0 has a rotation-invariant probability distribution, and so
linearity of expectation gives
n
 > 2 1 X

 1
1
2
E (e1 q0 ) =
E (e>
= Ekq0 k22 = .
i q0 )
n
n
n
i=1

Thus, in order to make π small, we should expect to have   1/n. The following lemma gives
that such choices of  suffice for π to be small:
√
Lemma 13. If  ≥ n−1 e−2n , then π ≤ 3 n.
2
2
Proof. First, observe that (e>
1 q0 ) is equal in distribution to Z /Q, where Z has standard normal
distribution and Q has chi-squared distribution with n degrees of freedom
√ (Z and Q are independent). The probability density function of Z has a maximal value of 1/ 2π at zero, and so

 r 2a
Pr Z 2 < a ≤
.
π

Also, Lemma 1 in [14] gives


√
Pr Q ≥ n + 2 nx + 2x ≤ e−x

∀x > 0.

Therefore, picking a = 5n and x = n, the union bound gives
Pr



2
(e>
1 q0 )

 2





 r 10n
√
Z
2
<  = Pr
<  ≤ Pr Z < 5n + Pr Q > 5n ≤
+ e−n ≤ 3 n.
Q
π

Overall, if we take  = n−(2c+1) for c > 0, then if H0 is true, our detector will produce a false
positive with probability O(n−c ). On the other hand, if H1 is true, then with probability 1−O(n−c ),
our detector will reject H0 after Oδ (c log n) power iterations, provided |λ2 | ≤ (1 − δ)|λ1 |.

3.2

Testing optimality with the power iteration detector

In this subsection, we leverage the power iteration detector to test k-means optimality. Note that
the sufficient condition (12) holds if and only if v := √1N 1 is a leading eigenvector of the matrix
A :=

z >
z
11 + PΛ⊥ (B − M )PΛ⊥ = 11> + PΛ⊥ (B − D)PΛ⊥ .
N
N

(15)

(The second equality follows from distributing the PΛ⊥ ’s and recalling the definition of M in (8).)
m
As such, it suffices that (A, v) satisfy H1 in (14). Overall, given a collection of points {xi }N
i=1 ⊆ R
and a proposed partition A1 t · · · t Ak = {1, . . . , N }, we can produce the corresponding matrix A
(defined above) and then run the power iteration detector of the previous subsection to test (12). In
particular, a positive test with tolerance  will yield ≥ 1−π confidence that the proposed partition is
optimal under the k-means objective. Furthermore, as we detail below, the matrix–vector products
computed in the power iteration detector have a computationally cheap implementation.
Given an m × na matrix Φa = [xa,1 · · · xa,na ] for each a ∈ {1, . . . , k}, we follow the following
procedure to implement the corresponding function x 7→ Ax as defined in (15):
1. Compute νa ∈ Rna such that (νa )i = kxa,i k22 for every a ∈ {1, . . . , k} in O(mN ) operations.
Let ν ∈ RN denote the vector whose ath block is νa .
14

>
2. Define the function (a, b, x) 7→ D(a,b) x such that D(a,b) = νa 1> − 2Φ>
a Φb + 1νb .
Running this function costs O(m(na + nb )) operations.

3. Define the function x 7→ Dx such that D = ν1> − 2Φ> Φ + 1ν > , where Φ = [Φ1 · · · Φk ].
Running this function costs O(mN ) operations.
4. Compute µa = 12 ( n12 11> −
a

2
(a,a) 1
na I)D

for every a ∈ {1, . . . , k} in O(mN ) operations.

5. Define the function (a, b, x) 7→ M (a,b) x such that M (a,b) = D(a,b) + µa 1> + 1µ>
b .
Running this function costs O(m(na + nb )) operations.
6. Compute z = mina6=b

2na
na +nb

min(M (a,b) 1) in O(kmN ) operations.

b
1 for every a, b ∈ {1, . . . , k}, a 6= b in O(kmN ) operations.
7. Compute u(a,b) = M (a,b) 1 − z na2n+n
a

8. Compute ρ(a,b) = u>
(a,b) 1 for every a, b ∈ {1, . . . , k}, a 6= b in O(kN ) operations.
9. Define the function x 7→ Bx such that the ath block of the output is given by
(Bx)a =

>
k u
X
(a,b) u(b,a) xb
b=1
b6=a

ρ(b,a)

.

Running this function costs O(kmN ) operations.
10. Define the function x 7→ PΛ⊥ x such that PΛ⊥ = I −
Running this function costs O(N ) operations.

Pk

1
>
a=1 na 1a 1a .

11. Define the function x 7→ Ax such that A = Nz 11> + PΛ⊥ (B − D)PΛ⊥ .
Running this function costs O(kmN ) operations.
Overall, after O(kmN ) operations of preprocessing, one may compute the function x 7→ Ax for
any given x in O(kmN ) operations. (Observe that this is the same complexity as each iteration of
Lloyd’s algorithm, and as we illustrate in Figure 2, the runtimes are comparable.)
At this point, we take a short aside to illustrate the utility of the power iteration detector
beyond k-means clustering. The original problem for which a PCC algorithm was developed was
community recovery under the stochastic block model [5]. For this random graph, there are
two communities of vertices, each of size n/2, and edges are drawn independently at random with
probability p if the pair of vertices belong to the same community, and with probability q < p if
they come from different communities. Given the random edges, the maximum likelihood estimator
for the communities is given by the vertex partition of two sets of size n/2 with the minimum cut.
Given a partition of the vertices, let X denote the corresponding n × n matrix of ±1s such that
Xij = 1 precisely when i and j belong to the same community. Given the adjacency matrix A of
the random graph, one may express the cut of a partition X in terms of Tr(AX). Furthermore,
X satisfies the convex constraints Xii = 1 and X  0, and so one may relax to these constraints
to obtain a semidefinite program and hope that the relaxation is typically tight over a large region
of (p, q). Amazingly, this relaxation is typically tight precisely over the region of (p, q) for which
community recovery is information-theoretically possible [1].
Given A, put B := 2A − 11> + I,
given a vector x ∈ Rn , define the corresponding n × n
Pand
n
diagonal matrix Dx by (Dx )ii := xi j=1 Bij xj . In [5], Bandeira observes that, given a partition

15

matrix X by some means (such as the fast algorithm provided in [2]), then X = xx> is SDPoptimal if both x> 1 = 0 and the second smallest eigenvalue of Dx − B is strictly positive, meaning
the partition gives the maximum likelihood estimator for the communities. However, as Bandeira
notes, the computational bottleneck here is estimating the second smallest eigenvalue of Dx − B,
and he suggests that a randomized power method–like algorithm might suffice, but leaves the
investigation for future research.
Here, we show how the power iteration detector fills this void in the theory. First, we note
that in the interesting regime of (p, q), the number of nonzero entries in A is O(n log n) with high
probability [1]. As such, the function x 7→ Bx can exploit this sparsity to take only O(n log n)
operations. This in turn allows for the computation of the diagonal of Dx to cost O(n log n)
operations. Next, note that
kDx − Bk2→2 ≤ kDx k2→2 + k2A − 11> k2→2 + kIk2→2
≤ kDx k2→2 + k2A − 11> kF + 1 = max |(Dx )ii | + n + 1 =: λ,
i

and that λ can be computed in O(n) operations after computing the diagonal of Dx . Also, it takes
O(n) operations to verify x> 1 = 0. Assuming x> 1 = 0, then the second smallest eigenvalue of
Dx − B is strictly positive if and only if x spans the unique leading eigenspace of λI − Dx + B.
Thus, one may test this condition using the power iteration detector, and furthermore, each iteration
will take only O(n log n) operations, thanks to the sparsity of A.

4

A fast k-means solver for two clusters

The previous section illustrated how to quickly test whether a proposed solution to the k-means
problem is optimal. In particular, this test will be successful with high probability if the data
follows the stochastic ball model with ∆ > 2 + k 2 /m. It remains to find a fast k-means solver which
also performs in this regime.
In doing so, we maintain the philosophy that our algorithm should not “see” the stochastic ball
model. Indeed, we view the stochastic ball model as a method of evaluating clustering algorithms
rather than a realistic data model. For example, Lloyd’s algorithm can be viewed as an alternating
minimization of the lifted objective function:
f (A1 , . . . , Ak , c1 , . . . , ck ) :=

k X
X

kxi − ct k2 ,

A1 t · · · t Ak = {1, . . . , N }, c1 , . . . , ck ∈ Rm ,

t=1 i∈At

and since this function is minimized at the k-means optimizer (regardless of how the data is distributed), such an algorithm is acceptable. On the other hand, one might consider matching the
stochastic ball model to the data by maximizing the following function:
g(c1 , . . . , ck ) :=

N X
k
X

pD (xi − ct ),

c1 , . . . , ck ∈ Rm ,

i=1 t=1

where pD (·) denotes the density function of D, which is supported on the unit ball centered at the
origin. One could certainly devise a fast greedy method such as matching pursuit [16] to optimize
this objective function (especially if pD is smooth), but doing so violates our philosophy.
In [20], Peng and Wei showed that k-means is equivalent to the following program:
minimize
subject to

Tr(DX)
X

>

(16)
2

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

One may quickly observe that the SDP (3) we analyzed in Section 2 is a relaxation of this program.
In this section, we follow Peng and Wei [20] by considering another relaxation of (16), obtained by
discarding the X ≥ 0 constraint (this is known as the spectral clustering relaxation [7, 8]). We
first denote the m × N matrix Φ = [x1 · · · xN ]. Without loss of generality, the data set is centered
at the origin so that Φ1 = 0. Letting ν denote the N × 1 vector with νi = kxi k22 , then
2
>
>
>
Dij = kxi − xj k22 = kxi k22 − 2x>
i xj + kxj k2 = (ν1 − 2Φ Φ + 1ν )ij .

As such, D = ν1> − 2Φ> Φ + 1ν > , and so the constraints X = X > and X1 = 1 together imply an
alternative expression for the objective function:
Tr(DX) = Tr(ν1> X − 2Φ> ΦX + 1ν > X)
= Tr(ν1> X > ) − 2 Tr(Φ> ΦX) + Tr(X1ν > )
= 2ν > 1 − 2 Tr(Φ> ΦX).
We conclude that minimizing Tr(DX) is equivalent to maximizing Tr(Φ> ΦX).
Next, we observe that the feasible X in our relaxation are precisely the rank-k N ×N orthogonal
projection matrices satisfying X1 = 1. This in turn is equivalent to X having the form X =
1
>
N 11 + Y , where Y is a rank-(k − 1) N × N orthogonal projection matrix satisfying Y 1 = 0.
Discarding the Y 1 = 0 constraint produces the following relaxation of (16):
maximize
subject to

Tr(Φ> ΦY )
Y

>

(17)
2

= Y, Y = Y, Tr(Y ) = k − 1

For general values of k, this program amounts to finding k − 1 principal components of the data.
Recalling our initial clustering goal, after finding the optimal Y , it remains to take X = N1 11> + Y
and then round to a nearby member of the feasibility region in (16). In [20], Peng and Wei focus
on the k = 2 case; they reduce the rounding step to a 2-means problem on the real line, and they
establish an approximation ratio of 2 for this relax-and-round procedure. Here, we are concerned
with exact recovery under the stochastic ball model, and as such, we slightly modify the rounding
step.
When k = 2, the solution to (17) has the form Y = yy > , where y is a leading unit eigenvector
1
1
>
of Φ> Φ. Our task is to find a matrix of the form |A|
1A 1>
A + |B| 1B 1B with A t B = {1, . . . , N } that
is close to

1
>
N 11

+ yy > . To this end, it seems natural to consider
Aθ := {i : yi < θ},

Bθ := Acθ

for some threshold θ. Since the data is centered (Φ1 = 0), one may be inclined to take θ = 0, but
this will be a poor choice if the true clusters have significantly different numbers of points. Instead,
we select the θ which minimizes the k-means objective of (Aθ , Bθ ). Since we only need to consider
N − 1 choices of θ, this is plausibly tractable, although computing the k-means objective once costs
O(mN ) operations, and so some care is necessary to keep the algorithm fast.
We will show how to find the optimal (Aθ , Bθ ) in O((m + log N )N ) operations using a simple
dynamic program. Order the indices so that y1 ≤ · · · ≤ yN . Then the function to minimize is
f (i) :=

i
i
N
N
X
X
1
1XX
kxj − xj 0 k22 +
kxj − xj 0 k22 .
i
N
−
i
j=1 j 0 =1
j=i+1 j 0 =i+1
|
{z
}
|
{z
}
vic

vi

17

Algorithm 2: Spectral k-means clustering (for two clusters)
Input: m × N matrix Φ = [x1 · · · xN ] of points to be clustered
Output: Clusters AP
t B = {1, . . . , N }
1
Subtract centroid N N
i=1 xi from each column of Φ to produce Φ0
Compute leading eigenvector y of Φ>
0 Φ0
Find θ that minimizes the k-means objective of ({i : yi < θ}, {i : yi ≥ θ})
(A, B) ← ({i : yi < θ}, {i : yi ≥ θ})

Expanding the square and distributing sums gives
vi+1 = vi + 2

i
X

kxj k22

−

4x>
i+1

j=1

i
X

xj + 2ikxi+1 k22 ,

j=1

and the vic ’s satisfy a similar recursion rule. As such, one may iteratively compute the vi ’s and
vic ’s before computing the f (i)’s and then minimizing. Overall, the following procedure finds the
optimal (Aθ , Bθ ) in O((m + log N )N ) operations:
1. Sort the entries y1 ≤ · · · ≤ yN in O(N log N ) operations.
2. Iteratively compute
s1 (i) :=

i
X
j=1

xj ,

sc1 (i) :=

N
X

xj ,

s2 (i) :=

j=i+1

i
X
j=1

kxj k22 ,

sc2 (i) :=

N
X

kxj k22

j=i+1

for every i ∈ {1, . . . , N − 1} in O(mN ) operations.
2
3. Compute v1 = 0 and vi+1 = vi + 2s2 (i) − 4x>
i+1 s1 (i) + 2ikxi+1 k2 for every i ∈ {1, . . . , N − 2}
in O(mN ) operations.
c
c
c
c
> c
2
4. Compute vN
−1 = 0 and vi−1 = vi +2s2 (i)−4xi s1 (i)+2(N −i)kxi k2 for every i ∈ {N −1, . . . , 2}
in O(mN ) operations.

5. Compute f (i) = vi /i + vic /(N − i) for every i ∈ {1, . . . , N − 1} in O(N ) operations.
6. Find i that minimizes f (i) and output {1, . . . , i} and {i + 1, . . . , N } in O(N ) operations.
Note that in the special case where m = 1, the above method exactly solves the k-means problem
when k = 2 in only O(N log N ) operations, recovering the rounding step of Peng and Wei [20]. For
comparison, [23] leverages more sophisticated dynamic programming for the m = 1 case, but k is
arbitrary and the algorithm costs O(kN 2 ) operations.
See Algorithm 2 for a summary of our relax-and-round procedure. As a spectral method,
this algorithm enjoys quasilinear computational complexity; see Figure 2 for an illustration. In
particular, when computing the leading eigenvector of Φ>
0 Φ0 , each matrix–vector multiply in the
power method costs only O(mN ) operations. Furthermore, as the following result guarantees, this
algorithm performs well under the stochastic ball model:
Theorem 14. Let ∆? = ∆? (D, k) denote the smallest value for which ∆ > ∆? implies that minimizing the k-means objective recovers planted clusters under the (D, γ, n)-stochastic ball model with
probability 1 − e−ΩD,γ (n) . When k = 2, spectral k-means clustering (Algorithm 2) recovers planted
clusters under the stochastic ball model with probability 1 − e−ΩD,γ (n) provided ∆ > ∆? .
18

See Appendix C for the proof. The main idea is that the leading eigenvector of Φ0 Φ>
0 is biased
towards the difference between the ball centers, and as the following lemma establishes, this bias
encourages spectral k-means clustering to separate the planted clusters:
Lemma 15. Take two clusters contained in unit balls centered at γ and −γ with kγk2 > 1. If minimizing the k-means objective recovers these clusters, then spectral k-means clustering (Algorithm 2)
>
also recovers them, provided the leading eigenvector z of Φ0 Φ>
0 satisfies |γ z| > kzk2 .
Proof. Write Φ0 = Φ − µ1> , put θ := −µ> z, and observe that y = Φ>
0 z is a leading eigenvector of
Φ>
Φ
.
Then
0 0
yi = (xi − µ)> z = x>
(18)
i z+θ
for every i. Next, if |γ > z| > kzk2 , then a simple trigonometric argument gives that the balls (and
therefore the planted clusters) are separated by the hyperplane orthogonal to z. Combined with
(18), we then have that the clusters can be identified according to whether yi < θ or yi > θ. It
therefore suffices to minimize the k-means objective subject to partitions of this form (for arbitrary
thresholds θ), as so spectral k-means clustering succeeds.

5

Discussion

This paper discussed various facets of probably certifiably correct algorithms for k-means clustering.
There are still many questions that have yet to be answered:
• Let ∆? (D, k) denote the smallest value for which ∆ > ∆? implies that minimizing the k-means
objective recovers planted clusters under the (D, γ, n)-stochastic ball model with probability
1 − e−ΩD,γ (n) . What is ∆? ? It was conjectured in [4] that ∆? = 2, but as we demonstrated in
Subsection 2.3, this is not the case.
• Let ∆?SDP (D, k) denote the smallest value for which ∆ > ∆?SDP implies that solving the kmeans SDP recovers planted clusters under the (D, γ, n)-stochastic ball model with probability
1 − e−ΩD,γ (n) . What is ∆?SDP ? Considering Subsection 2.3 and Figure 3(center), we suspect
the SDP exhibits a performance gap: ∆?SDP > ∆? .
• Is there a single dual certificate for the k-means SDP that typically certifies planted clusters under the stochastic ball model whenever ∆ > ∆?SDP ? Does this certification have a
quasilinear-time implementation similar to Subsection 3.2?
• Is there a quasilinear-time k-means solver that typically solves k-means under the stochastic
ball model whenever ∆ > ∆? ? In particular, is there a quasilinear-time initialization of
Lloyd’s algorithm that meets this specification? Following the philosophy of Section 4, such
algorithms should be designed so as to not “see” the stochastic ball model.

Acknowledgments
The authors thank the anonymous referees, whose suggestions significantly improved this paper’s
presentation and literature review. The authors also thank Afonso S. Bandeira and Nicolas Boumal
for interesting discussions and valuable comments on an earlier version of this manuscript. DGM
was supported by an AFOSR Young Investigator Research Program award, NSF Grant No. DMS1321779, and AFOSR Grant No. F4FGA05076J002. SV was supported by Rachel Ward’s NSF
CAREER award and AFOSR Young Investigator Research Program award. The views expressed
19

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
[1] E. Abbe, A. S. Bandeira, and G. Hall. Exact recovery in the stochastic block model. arXiv
preprint arXiv:1405.3267, 2014.
[2] E. Abbe and C. Sandon. Community detection in general stochastic block models: fundamental
limits and efficient recovery algorithms. arXiv preprint arXiv:1503.00609, 2015.
[3] D. Arthur and S. Vassilvitskii. k-means++: the advantages of careful seeding. In Proceedings
of the eighteenth annual ACM-SIAM symposium on Discrete algorithms, 2007.
[4] P. Awasthi, A. Bandeira, M. Charikar, K. Ravishankar, S. Villar, and R. Ward. Relax, no
need to round: Integrality of clustering formulations. http://arxiv.org/abs/1408.4045, 2014.
[5] A. S. Bandeira. A note on probably certifiably correct algorithms.
arXiv:1509.00824, 2015.

arXiv preprint

[6] H. Chen and J. Peng. 0–1 semidefinite programming for graph-cut clustering: modelling and
approximation. Data Mining and Mathematical Programming. CRM Proceedings and Lecture
Notes of the American Mathematical Society, pages 15–40, 2008.
[7] I. S. Dhillon, Y. Guan, and B. Kulis. Kernel k-means: spectral clustering and normalized cuts.
In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery
and data mining, pages 551–556. ACM, 2004.
[8] I. S. Dhillon, Y. Guan, and B. Kulis. Weighted graph cuts without eigenvectors a multilevel
approach. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 29(11):1944–
1957, 2007.
[9] E. Elhamifar, G. Sapiro, and R. Vidal. Finding exemplars from pairwise dissimilarities via
simultaneous sparse recovery. Advances in Neural Information Processing Systems, pages 19–
27, 2012.
[10] G. H. Golub and C. F. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
[11] M. Grant, S. Boyd, and Y. Ye. Cvx: Matlab software for disciplined convex programming,
2008.
[12] T. Iguchi, D. G. Mixon, J. Peterson, and S. Villar. On the tightness of an sdp relaxation of
k-means. arXiv preprint arXiv:1505.04778, 2015.
[13] K. Jain, M. Mahdian, and A. Saberi. A new greedy approach for facility location problems.
In Proceedings of the 34th Annual ACM Symposium on Theory of Computing, 2002.
[14] B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model selection.
Annals of Statistics, pages 1302–1338, 2000.
[15] S. Lloyd. Least squares quantization in pcm. Information Theory, IEEE Transactions on,
28(2):129–137, 1982.
20

[16] S. G. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries. Signal Processing, IEEE Transactions on, 41(12):3397–3415, 1993.
[17] A. Nellore and R. Ward. Recovery guarantees for exemplar-based clustering. arXiv:1309.3256,
2013.
[18] Y. Nesterov, A. Nemirovskii, and Y. Ye. Interior-point polynomial algorithms in convex programming, volume 13. SIAM, 1994.
[19] R. Ostrovsky, Y. Rabani, L. Schulman, and C. Swamy. The effectiveness of lloyd-type methods
for the k-means problem. In Proceedings of the 47th Annual IEEE Symposium on Foundations
of Computer Science, 2006.
[20] J. Peng and Y. Wei. Approximating k-means-type clustering via semidefinite programming.
SIAM Journal on Optimization, 18(1):186–205, 2007.
[21] J. A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computational Mathematics, 12(4):389–434, 2012.
[22] R. Vershynin.
Introduction to the non-asymptotic analysis of random matrices.
arXiv:1011.3027v7, 2011.
[23] H. Wang and M. Song. Ckmeans.1d.dp: optimal k-means clustering in one dimension by
dynamic programming. The R Journal, 3(2):29–33, 2011.

A

Proof of Corollary 8

It suffices to have
kPΛ⊥ M PΛ⊥ k2→2 + kPΛ⊥ BPΛ⊥ k2→2 ≤ z.

(19)

We will bound the terms in (19) separately and then combine the bounds to derive a sufficient
condition for Theorem 7. To bound the first term in (19), let ν be the N × 1 vector whose (a, i)th
entry is kxa,i k22 , and let Φ be the m × N matrix whose (a, i)th column is xa,i . Then
2
>
>
>
D(a,i),(b,j) = kxa,i − xb,j k22 = kxa,i k22 − 2x>
a,i xb,j + kxb,j k2 = (ν1 − 2Φ Φ + 1ν )(a,i),(b,j) ,

meaning D = ν1> − 2Φ> Φ + 1ν > . With this, we appeal to the blockwise definition of M (8):
kPΛ⊥ M PΛ⊥ k2→2 = kPΛ⊥ DPΛ⊥ k2→2 = kPΛ⊥ (ν1> − 2Φ> Φ + 1ν > )PΛ⊥ k2→2
= 2kPΛ⊥ Φ> ΦPΛ⊥ k2→2 = 2kΦPΛ⊥ k22→2 = 2kΨk22→2 .
For the second term in (19), we first write the decomposition
B=

k
k

X
X


H(a,b) (B (a,b) ) + H(b,a) (B (b,a) ) ,

a=1 b=a+1

where H(a,b) : Rna ×nb → RN ×N produces a matrix whose (a, b)th block is the input matrix, and is
otherwise zero. Then
k
k


X
X
PΛ⊥ BPΛ⊥ =
PΛ⊥ H(a,b) (B (a,b) ) + H(b,a) (B (b,a) ) PΛ⊥
a=1 b=a+1

=

k
k

X
X


H(a,b) (P1⊥ B (a,b) P1⊥ ) + H(b,a) (P1⊥ B (b,a) P1⊥ ) ,

a=1 b=a+1

21

and so the triangle inequality gives
kPΛ⊥ BPΛ⊥ k2→2 ≤

k
k
X
X

kH(a,b) (P1⊥ B (a,b) P1⊥ ) + H(b,a) (P1⊥ B (b,a) P1⊥ )k2→2

a=1 b=a+1

=

k
k
X
X

kP1⊥ B (a,b) P1⊥ k2→2 ,

a=1 b=a+1

where the last equality can be verified by considering the spectrum of the square:

2
H(a,b) (P1⊥ B (a,b) P1⊥ ) + H(b,a) (P1⊥ B (b,a) P1⊥ )




= H(a,a) (P1⊥ B (a,b) P1⊥ )(P1⊥ B (a,b) P1⊥ )> + H(b,b) (P1⊥ B (a,b) P1⊥ )> (P1⊥ B (a,b) P1⊥ ) .
At this point, we use the definition of B (11) to get
kP1⊥ B (a,b) P1⊥ k2→2 =

kP1⊥ u(a,b) k2 kP1⊥ u(b,a) k2
.
ρ(a,b)

Recalling the definition of u(a,b) (11) and combining these estimates then produces the result.

B

Proof Theorem 9

In this section, we apply the certificate from Corollary 8 to the (D, γ, n)-stochastic ball model (see
Definition 2) to prove our main result. We will prove Theorem 9 with the help of several lemmas.
Lemma 16. Denote
n

ca :=

1X
xa,i ,
n

∆ab := kγa − γb k2 ,

Oab :=

i=1

γa + γb
.
2

Then the (D, γ, n)-stochastic ball model satisfies the following estimates:
kca − γa k2 < 

w.p.

1 − e−Ωm, (n)

(20)

kra,i k22 − Ekrk22 < 

w.p.

1 − e−Ω (n)

(21)

1X
kxa,i − Oab k22 − Ekr + γa − Oab k22 < 
n

w.p.

1 − e−Ω∆ab , (n)

(22)

1
n

n
X
i=1

n

i=1

Proof. Since Er = 0 and krk22 ≤ 1 almost surely, one may lift


>
0 ra,i
Xa,i :=
ra,i 0
and apply the Matrix Hoeffding inequality [21] to conclude that
 X

n
2
Pr
ra,i ≥ t ≤ me−t /8n .
i=1

2

Taking t := n then gives (20). For (21) and (22), notice that the random variables in each
sum are iid and confined to an interval almost surely, and so the result follows from Hoeffding’s
inequality.
22

Lemma 17. Under the (D, γ, n)-stochastic ball model, we have D(a,b) 1 − D(a,a) 1 = 4np + q, where
∆2
>
pi := ra,i
(γa − Oab ) + ab
4

 X

n
n
X
>
2
2
qi := 2n(xa,i − Oab ) (ca − cb ) − (γa − γb ) +
kxb,j − Oab k2 −
kxa,j − Oab k2
j=1

j=1

and |qi | ≤ (6 + 2∆ab )n with probability 1 − e−Ωm,∆ab , (n) .
Proof. Add and subtract Oab and then expand the squares to get
(a,b)
e>
1
i (D

−D

(a,a)

1) =

n
X

kxa,i −

j=1

xb,j k22

−

n
X

kxa,i − xa,j k22

j=1


n
1X
2
kxb,j − Oab k2
= n − 2(xa,i − Oab ) (cb − Oab ) +
n
j=1


n
1X
2
>
kxa,j − Oab k2
− n − 2(xa,i − Oab ) (ca − Oab ) +
n
j=1
X

n
n
X
>
2
2
= 2n(xa,i − Oab ) (ca − cb ) +
kxb,j − Oab k2 −
kxa,j − Oab k2 .


>

j=1

j=1

Add and subtract γa − γb to ca − cb and distribute over the resulting sum to obtain
(a,b)
e>
1 − D(a,a) 1) = 2n(xa,i − Oab )> (γa − γb ) + q
i (D

>
= 4n ra,i + (γa − Oab ) (γa − Oab ) + q.

Distributing and identifying kγa − Oab k22 = ∆2ab /4 explains the definition of p. To show |qi | ≤
(6 + 2∆ab )n, apply triangle and Cauchy–Schwarz to obtain


n
n
X
X
2
|qi | ≤ 2n(xa,i − Oab ) (ca − cb ) − (γa − γb ) +
kxb,j − Oab k2 −
kxa,j − Oab k22
>

j=1



j=1



≤ 2n kra,i k2 + kγa − Oa,b k2


kca − γa k2 + kcb − γb k2

+

n
X

kxb,j − Oab k22 −

j=1

n
X

kxa,j − Oab k22

j=1



n
n
X
X
∆ab
kca − γa k2 + kcb − γb k2 +
kxb,j − Oab k22 −
kxa,j − Oab k22 .
≤ 2n 1 +
2


j=1

j=1

To finish the argument, apply (20) to the first term while adding and subtracting
Ekr + γa − Oab k22 = Ekr + γb − Oab k22 ,
from the second and apply (22).
Lemma 18. Under the (D, γ, n)-stochastic ball model, we have
1 > (a,a)
1 D
1 − 2nEkrk22 ≤ 4n
n
23

w.p.

1 − e−Ω∆ab , (n) .

Proof. Add and subtract γa and expand the square to get
n

n

j=1

j=1

1 > (a,a)
1X
1X
>
ei D
1=
kxa,i − xa,j k22 = kra,i k22 − 2ra,i
(ca − γa ) +
kra,j k22 .
n
n
n
The triangle and Cauchy–Schwarz inequalities then give
1 > (a,a)
1 D
1 − 2nEkrk22
n

n
n 
X
1X
2
>
2
=
kra,i k2 − 2ra,i (ca − γa ) +
kra,j k2 − 2nEkrk22
n
j=1

i=1

≤n

≤n

1
n
1
n

n
X
i=1
n
X

kra,i k22 − Ekrk22 + 2
kra,i k22 − Ekrk22 + 2

i=1

n
X
i=1
n
X

n

1X
kra,j k22 − Ekrk22
n

>
|ra,i
(ca − γa )| + n

kca − γa k2 + n

i=1

1
n

j=1
n
X

kra,j k22 − Ekrk22

j=1

≤ 4n,
where the last step occurs with probability 1 − e−Ω∆ab , (n) by a union bound over (21) and (20).
Lemma 19. Under the (D, γ, n)-stochastic ball model, we have
1> D(a,b) 1 − 1> D(a,a) 1 ≥ n2 ∆2ab − (6 + 4∆ab )n2 

w.p.

1 − e−Ωm,∆ab , (n) .

Proof. Lemma 17 gives
1> D(a,b) 1 − 1> D(a,a) 1 = 1> (4np + q)

n 
X
∆2ab
>
≥ 4n
ra,i (γa − Oab ) +
− (6 + 2∆ab )n2 
4
i=1


n∆2ab
>
≥ 4n n(ca − γa ) (γa − Oab ) +
− (6 + 2∆ab )n2 .
4
Cauchy–Schwarz along with (20) then gives the result.
Lemma 20. Under the (D, γ, n)-stochastic ball model, there exists C = C(γ) such that
min(M (a,b) 1) ≥ n∆(∆ − 2) + Cn

min

w.p.

1 − e−Ωm,γ, (n) ,

a,b∈{1,...,k}
a6=b

where ∆ :=

min
a,b∈{1,...,k}
a6=b

∆ab .

Proof. Fix a and b. Then by Lemma 17, the following holds with probability 1 − e−Ωm,∆ab , (n) :




∆2
>
min D(a,b) 1 − D(a,a) 1 ≥ 4n min
ra,i
(γa − Oab ) + ab − (6 + 2∆ab )n
4
i∈{1,...,n}
≥ n∆2ab − 2n∆ab − (6 + 2∆ab )n,
24

where the last step is by Cauchy–Schwarz. Taking a union bound with Lemma 18 then gives
min(M (a,b) 1)


 1 1
1 > (b,b)
> (a,a)
(a,b)
(a,a)
1 D
1− 1 D
1
= min D
1−D
1 +
2 n
n


 1 1
1 > (b,b)
(a,b)
(a,a)
> (a,a)
2
2
≥ min D
1−D
1 −
1 D
1 − 2nEkrk2 + 1 D
1 − 2nEkrk2
2 n
n
≥ n∆ab (∆ab − 2) − (10 + 2∆ab )n
with probability 1 − e−Ω∆ab , (n) . The result then follows from a union bound over a and b.
Lemma 21. Suppose  ≤ 1. Then there exists C = C(∆ab , m) such that under the (D, γ, n)stochastic ball model, we have
kP1⊥ M (a,b) 1k22 ≤

4n3 ∆2ab
+ Cn3 
m

with probability 1 − e−Ωm,∆ab , (n) .
Proof. First, a quick calculation reveals


1 > (b,b)
1 1 > (a,a)
1 D
1− 1 D
1 ,
=
−
+
2 n
n


1 > (a,b)
1
1 1 > (a,a)
1
1 M
1 = 1> D(a,b) 1 −
1 D
1 + 1> D(b,b) 1 ,
n
n
2 n
n
(a,b)
e>
1
i M

(a,b)
e>
1
i D

(a,a)
e>
1
i D

from which it follows that
1
(a,b)
(a,b)
e>
1 = e>
1 − 1> M (a,b) 1
i P1⊥ M
i M
n

 

1 > (a,b)
1 > (a,a)
> (a,b)
> (a,a)
= ei D
1− 1 D
1 − ei D
1− 1 D
1
n
n
(a,b)
= e>
1 − D(a,a) 1).
i P1⊥ (D

As such, we have
kP1⊥ M (a,b) 1k22 = kP1⊥ (D(a,b) 1 − D(a,a) 1)k22
= kD(a,b) 1 − D(a,a) 1k22 − kP1 (D(a,b) 1 − D(a,a) 1)k22 .

(23)

To bound the first term, we apply the triangle inequality over Lemma 17:
kD(a,b) 1 − D(a,a) 1k2 ≤ 4nkpk2 + kqk2 ≤ 4nkpk2 + (6 + 2∆ab )n3/2 .

(24)

We proceed by bounding kpk2 . To this end, note that the pi ’s are iid random variables whose outcomes lie in a finite interval (of width determined by ∆ab ) with probability 1. As such, Hoeffding’s
inequality gives
n
1X 2
pi − Ep21 ≤ 
w.p.
1 − e−Ω∆ab , (n) .
n
i=1

25

With this, we then have
kpk22

 X

n
1
2
2
2
=n
pi − Ep1 + Ep1 ≤ nEp21 + n
n

(25)

i=1

in the same event. To determine Ep21 , first take r1 := e>
1 r. Then since the distribution of r is
rotation invariant, we may write
>
p1 = ra,1
(γa − Oab ) + kγa − Oab k22 =

∆2
∆ab
r1 + ab ,
2
4

where the second equality above is equality in distribution. We then have


∆2ab 2 ∆2ab 2 ∆4ab
∆ab
2
=
Ep1 = E
r1 +
Er1 +
.
2
4
4
16

(26)

We also note that 1 ≥ Ekrk22 = mEr12 by linearity of expectation, and so
Er12 ≤

1
.
m

(27)

Combining (24), (25), (26) and (27) then gives
kD

(a,b)

1−D

(a,a)


1k2 ≤

1/2
4n3 ∆2ab
3 4
3
+ n ∆ab + 16n 
+ (6 + 2∆ab )n3/2 .
m

(28)

To bound the second term of (23), first note that
1
kP1 (D(a,b) 1 − D(a,a) 1)k2 = √ 1> D(a,b) 1 − 1> D(a,a) 1 .
n

(29)

1> D(a,b) 1 − 1> D(a,a) 1 ≥ 1> D(a,b) 1 − 1> D(a,a) 1 ≥ n2 ∆2ab − (6 + 4∆ab )n2 

(30)

Lemma 19 then gives

with probability 1 − e−Ωm,∆ab , (n) . Using (23) to combine (28) with (29) and (30) then gives the
result.
Lemma 22. There exists C = C(γ) such that under the (D, γ, n)-stochastic ball model, we have

ρ(a,b) ≥ n2 ∆2ab − ∆(∆ − 2) − Cn2 
w.p.
1 − e−ΩD,γ, (n) .
Proof. Recall from (11) that
>
(a,b)
ρ(a,b) = u>
1 − nz = 1> M (a,b) 1 − n
(a,b) 1 = 1 M

min

min(M (a,b) 1).

a,b∈{1,...,k}
a6=b

To bound the first term, we leverage Lemma 19:
1
1> M (a,b) 1 = 1> D(a,b) 1 − (1> D(a,a) 1 + 1> D(b,b) 1)
2
 1

1  > (a,b)
=
1 D
1 − 1> D(a,a) 1 +
1> D(b,a) 1 − 1> D(b,b) 1
2
2
≥ n2 ∆2ab − (6 + 4∆ab )n2 
26

(31)

with probability 1 − e−Ωm,∆ab , (n) . To bound the second term in (31), note from Lemma 18 that
min(M (a,b) 1)


 1 1
1 > (b,b)
> (a,a)
(a,b)
(a,a)
1 D
1− 1 D
1
= min D
1−D
1 +
2 n
n


 1 1
1 > (b,b)
(a,b)
(a,a)
> (a,a)
2
2
≤ min D
1−D
1 +
1 D
1 − 2nEkrk2 + 1 D
1 − 2nEkrk2
2 n
n


≤ min D(a,b) 1 − D(a,a) 1 + 4n
with probability 1 − e−Ω∆ab , (n) . Next, Lemma 17 gives


min D(a,b) 1 − D(a,a) 1 ≤ n∆2ab + (6 + 2∆ab )n + 4n

min
i∈{1,...,n}

>
ra,i
(γa − Oab ).

By assumption, we know krk2 ≥ 1 −  with positive probability regardless of  > 0. It then follows
that
∆ab
+
r> (γa − Oab ) ≤ −
2
with some (-dependent) positive probability. As such, we may conclude that
min
i∈{1,...,n}

>
(γa − Oab ) ≤ −
ra,i

∆ab
+
2

1 − e−ΩD, (n) .

w.p.

Combining these estimates then gives
min(M (a,b) 1) ≤ n∆2ab − 2n∆ab + (10 + 2∆ab )n

w.p.

1 − e−ΩD,∆ab , (n) .

Performing a union bound over a and b then gives
min

min(M (a,b) 1) ≤ n∆2 − 2n∆ + (10 + 2∆)n

w.p.

1 − e−ΩD,γ, (n) .

a,b∈{1,...,k}
a6=b

Combining these estimates then gives the result.
Lemma 23. Under the (D, γ, n)-stochastic ball model, we have


√
(1 + )σ
√
N
w.p.
1 − e−Ωm,k,σ, (n) ,
kΨk2→2 ≤
+
m
where σ 2 := Ekrk22 for r ∼ D.
Proof. Let R denote the matrix whose (a, i)th column is ra,i . Then
h
i
>
>
Ψ = R − (c1 − γ1 )1 · · · (ck − γk )1 ,
and so the triangle inequality gives
kΨk2→2 ≤ kRk2→2 +

h

>

(c1 − γ1 )1

>

· · · (ck − γk )1

i
2→2

 X
1/2
k
2
≤ kRk2→2 + n
kca − γa k2
,
a=1

where the last estimate passes to the Frobenius norm. For the first term, since D is rotation
invariant, we may apply Theorem 5.41 in [22]:
r
N
kRk2→2 ≤ (1 + )σ
w.p.
1 − e−Ωm,σ, (n) .
m
For the second term, apply (20). The union bound then gives the result.
27

Proof of Theorem 9. First, we combine Lemmas 21, 22 and 23: For every δ > 0, there exists an
 > 0 such that
2kΨk22→2

k
k
X
X
kP1⊥ M (a,b) 1k2 kP1⊥ M (b,a) 1k2
+
ρ(a,b)
a=1 b=a+1

2

k
k
X
X
4n3 ∆2ab /m + Cn3 
1+
≤ 2 √ +  nk +
m
n2 (∆2ab − ∆(∆ − 2)) − Cn2 
a=1
b=a+1



2k
4
≤n
+
m
m

k
X

k
X

a=1 b=a+1

∆2ab
+δ
∆2ab − ∆(∆ − 2)


(32)

with probability 1 − e−ΩD,γ, (n) . Next, the uniform bound ∆ab ≥ ∆ implies
∆2ab
1
1
∆
=
≤
= .
2
2
2
1 − ∆(∆ − 2)/∆
2
∆ab − ∆(∆ − 2)
1 − ∆(∆ − 2)/∆ab
Combining this with (32) and considering Lemma 20, it then suffices to have
 
2k
4
k
∆
+
·
·
< ∆(∆ − 2).
m
m
2
2
Rearranging then gives
k(k − 1)
2k
+
,
m∆
m
which is implied by the hypothesis since ∆ ≥ 2.
∆>2+

C

Proof of Theorem 14

>
Put g = γ/kγk2 and let z have unit 2-norm. Since kΦ>
0 zk2 ≥ kΦ0 gk2 , then considering Lemma 15,
it suffices to show that the containment

 

2
m−1
>
m−1
>
>
S1 := v ∈ S
: |hg vi| ≤
⊆ v∈S
: kΦ0 vk2 < kΦ0 gk2 =: S2
∆

holds with probability 1 − e−Ωm,∆ (N ) . To this end, we will first show that each v ∈ S1 is also a
member of S2 with high probability, and then we will perform a union bound over an -net of S1 .
We start by considering kΦ> vk2 and kΦ> gk2 . Decompose xi as either γ +ri or −γ +ri depending
on whether xi belongs to the ball centered at γ or −γ. Let w with kwk2 = 1 be arbitrary. Then
2
>
2
>
>
2
>
2
>
>
>
2
(x>
i w) = ((±γ + ri ) w) = (±γ w + ri w) = (γ w) ± 2(γ w)(ri w) + (ri w) ,
2
>
2
> 2
and so E(x>
i w) = (γ w) + E(e1 r) . Linearity of expectation then gives



4
2
> 2
E (x>
= (γ > g)2 − (γ > v)2 = kγk2 (1 − (g > v)2 ) ≥ 1 − 2 .
i g) − (xi v)
∆
2
> 2
2
Since |(x>
i g) − (xi v) | ≤ 2(1 + ∆/2) almost surely, we may apply Hoeffding’s inequality to get

kΦ

>

gk22

− kΦ

>

vk22

=

N 
X
i=1

2
(x>
i g)

−

2
(x>
i v)





4
≥ N 1 − 2 − s w.p.
∆
28

1 − e−Ω∆ (s

2 /N )

.

(33)

For a properly chosen t, rearranging gives that kΦ> vk2 < kΦ> gk2 . Instead, we will use (33) to prove
>
the closely related inequality kΦ>
0 vk2 < kΦ0 gk2 . Letting µ denote the centroid of the columns of
Φ, we know by (20) that kµk2 ≤ δ with probability 1 − e−Ωm,δ (N ) . In this event, every w with
kwk2 = 1 satisfies
> >
>
>
kΦ>
0 wk2 − kΦ wk2 = k(Φ + µ1 ) wk2 − kΦ wk2

= kΦ> w + 1µ> wk2 − kΦ> wk2 ≤ k1µ> wk2 ≤

√
N δ.

(34)

Furthermore,
kΦ>
0 wk2

> >

>

= k(Φ − µ1 ) wk2 ≤ kΦwk2 + k1µ wk2 ≤

√


∆
N
+ 1 + kµk2 ,
2


where the last inequality follows from Cauchy–Schwarz along with the fact that kxi k2 ≤ ∆/2 + 1
for every i. Taking a supremum over w then gives




√
√
∆
∆
kΦ>
k
≤
N
+
1
+
kµk
N
+
1
+
δ
≤
w.p. 1 − e−Ωm,δ (N ) .
(35)
2
0 2→2
2
2
In (33), pick s = (N/2)(1 − 4/∆2 ) =: c1 (∆)N . Then taking a union bound with (34) gives
√ 2
√ 2
≤ kΦ> vk22 ≤ kΦ> gk22 c1 (∆)N ≤ kΦ>
− c1 (∆)N
kΦ>
0 gk2 + N δ
0 vk2 − N δ
with probability 1 − e−Ωm,∆,δ (N ) . Expanding both sides and rearranging then gives
√

2
> 2
>
>
kΦ>
0 vk2 ≤ kΦ0 gk2 + 2 N δ kΦ0 vk2 + kΦ0 gk2 − c1 (∆)N



∆
> 2
+1+δ
≤ kΦ0 gk2 − c1 (∆) − 4δ
N,
2
|
{z
}
c2 (∆)

where the last step follows from (35). Thus, picking δ = δ(∆) sufficiently small ensures c2 (∆) > 0.
2
> 2
>
>
>
>
Since c2 (∆)N ≤ kΦ>
0 gk2 − kΦ0 vk2 = (kΦ0 gk2 + kΦ0 vk2 )(kΦ0 gk2 − kΦ0 vk2 ), we further have
>
kΦ>
0 gk2 − kΦ0 vk2 ≥

c2 (∆)N
>
kΦ>
gk
2 + kΦ0 vk2
0

√
≥ c3 (∆) N ,

where the last inequality takes c3 (∆) := c2 (∆)/(∆/2 + 1 + δ), following (35).
At this point, we know that if v ∈ S1 , then v ∈ S2 with probability 1 − e−Ωm,∆ (N ) . It remains
to perform a union bound over an -net of S1 to conclude that S1 ⊆ S2 with high probability. To
this end, pick  < c3 (∆)/(∆/2 + 1 + δ), consider an -net N of S1 , and suppose
√
>
kΦ>
∀v ∈ N .
(36)
0 vk2 ≤ kΦ0 gk2 − c3 (∆) N
Then for every x ∈ S1 , there exists v ∈ N such that kx − vk2 ≤ , and so (35) gives


√
√
∆
>
>
>
>
kΦ0 xk2 ≤ kΦ0 k2→2 kx − vk2 + kΦ0 vk2 ≤ N
+ 1 + δ  + kΦ>
0 gk2 − c3 (∆) N < kΦ0 gk2 ,
2
as desired. To measure the probability of the success event (36), a standard volume comparison
argument establishes the existence of an -net of size |N | ≤ (1 + 2/)m ; see Lemma 5.2 in [22]. As
such, the union bound gives that (36) occurs with probability 1 − e−Ωm,∆ (N ) .
29

