# Sketch-and-solve approaches to $k$-means clustering
by semidefinite programming

Charles Clum111Department of Mathematics, The Ohio State University, Columbus,
Ohio, USA Dustin G. Mixon11footnotemark: 1 222Translational Data Analytics
Institute, The Ohio State University, Columbus, Ohio, USA Soledad
Villar333Department of Applied Mathematics & Statistics, Johns Hopkins
University, Baltimore, Maryland, USA 444Mathematical Institute for Data
Science, Johns Hopkins University, Baltimore, Maryland, USA Kaiying
Xie555Department of Electrical and Computer Engineering, The Ohio State
University, Columbus, Ohio, USA

###### Abstract

We introduce a sketch-and-solve approach to speed up the Peng-Wei semidefinite
relaxation of $k$-means clustering. When the data is appropriately separated
we identify the $k$-means optimal clustering. Otherwise, our approach provides
a high-confidence lower bound on the optimal $k$-means value. This lower bound
is _data-driven_ ; it does not make any assumption on the data nor how it is
generated. We provide code and an extensive set of numerical experiments where
we use this approach to certify approximate optimality of clustering solutions
obtained by k-means++.

##  1 Introduction

One of the most fundamental data processing tasks is clustering. Here, one is
given a collection of objects, a notion of similarity between those objects,
and a clustering objective that scores any given partition according to how
well it clusters similar objects together. The goal is then to partition the
objects in a way that optimizes this clustering objective. For example, to
partition the vertices of a simple graph into two clusters, one might take the
clustering objective to be edge cut in the graph complement. This clustering
problem is equivalent to MAX-CUT, which is known to be NP-hard [21]. To
(approximately) solve MAX-CUT, one may pass to the Goemans–Williamson
semidefinite relaxation [15] and randomly round to a nearby partition. For
certain random graph models, it is known that this semidefinite relaxation is
tight with high probability, meaning the clustering problem is exactly solved
before the rounding step [1, 7].

The above discussion suggests a workflow to solve clustering problems: relax
to a semidefinite program (SDP) and round the solution to a partition if
necessary. Since SDPs can be solved to machine precision in polynomial time
[34], this has been a worthwhile pursuit for several clustering settings [1,
4, 31, 19, 32, 2]. However, the polynomial runtime of semidefinite programming
is notoriously slow in practice, making it infeasible to cluster thousands of
objects (say). As an alternative, one may pass to a random subset of the
objects, cluster the subset by semidefinite programming, and then infer a
partition of the full data set based on proximity to the small clusters. This
sketch-and-solve approach was studied in [32, 2] in the context of clustering
vertices of a graph into two clusters. By passing to a small subset of
objects, the SDP step is no longer burdensome, and one may prove performance
guarantees in terms of planted structure in the data.

In this paper, we develop analogous sketch-and-solve approaches to cluster
data in Euclidean space. Here, the standard clustering objective is the
$k$-means objective. Much like MAX-CUT, the $k$-means problem is NP-hard [3,
5], but one may relax to the Peng–Wei semidefinite relaxation [35] and obtain
performance guarantees [4, 19, 31, 27]. By sketching to a random subset of
data points, one may solve this relaxation quickly and then cluster the
remaining data points according to which cluster mean is closest. Taking
inspiration from [32], this approach was studied in [46] in the setting of
Gaussian mixture models, where exact recovery is possible provided the planted
clusters are sufficiently separated and the size of the sketch is
appropriately large. In Section 3, we consider data that are not necessarily
drawn from a random model. As we will see, the sketch-and-solve approach
exactly recovers planted clusters provided they are sufficiently separated and
the size of the sketch is appropriately large relative to the shape of the
clusters.

In practice, it is popular to solve the $k$-means problem using Lloyd’s method
[29] with the random initialization afforded by $k$-means++ [41]. This
algorithm produces a partition that is within an $O(\log k)$ factor of the
optimal partition in expectation. In principle, the semidefinite relaxation
could deliver an a posteriori approximation certificate for this partition,
potentially establishing that the partition is even closer to optimal than
guaranteed by the $k$-means++ theory. In Section 4, we describe sketch-and-
solve approaches to this process that certify constant-factor approximations
for data drawn from Gaussian mixture models whenever the dimension of the
space is much larger than the number of clusters (even when the clusters are
not separated!). In addition, as established in Section 5, these algorithms
give improved bounds for most of the real-world data sets considered in the
original $k$-means++ paper [41].

###  1.1 Summary of contributions

What follows is a brief description of this paper’s contributions:

  * •

A sketch-and-solve algorithm for $k$-means clustering that consists of
subsampling the input dataset, solving an SDP on the samples, and
extrapolating the solution on the samples to the entire dataset. We provide
theoretical guarantees when the data points are drawn from separated clusters.
We phrase the conditions of optimality in terms of _proximity conditions_ , a
classical concept within the computer science literature [24]. These results
make use of work from [27] that connects the clustering SDP with proximity
conditions. (See Section 3.)

  * •

An algorithm that takes a dataset as input and outputs a lower bound on the
$k$-means optimal value. This algorithm leverages the sketch-and-solve scheme
and can be used to certify approximate optimality of clusters obtained with
fast, greedy methods such as $k$-means++. We provide bounds on the tightness
of this lower bound when data is sampled from spherical Gaussians. (See
Section 4.)

  * •

Open source code and an extensive set of numerical experiments showing how
tight the lower bounds are on several real-world datasets.666Our code is
available here: https://github.com/Kkylie/Sketch-and-solve_kmeans.git

###  1.2 Related work

The classical semidefinite programming relaxation for $k$-means clustering was
proposed by Peng and Wei in [35]. A few years later, [4] proved that the
Peng–Wei relaxation is tight (i.e., it recovers the solution to the original
NP-hard problem) if the clusters are sampled from the so-called stochastic
ball model [33] with sufficiently separated balls. The proof is based on the
construction of a dual certificate. These results were significantly improved
by [19] and [27], and applied to spectral methods in [28]. In particular, [27]
shows a connection between the tightness of the SDP and the proximity
conditions from theoretical computer science [24]. The actual threshold at
which the SDP becomes tight is currently unknown, an open conjecture was posed
in [27].

The Peng–Wei SDP relaxation has also been studied in the context of clustering
mixtures of Gaussians, a classical problem in theoretical computer science.
Introduced in [9], this problem has been approached with many methods
including spectral-like methods [20], methods of moments [18, 44], integer
programming [10], and it has been studied from an information-theoretic point
of view [12] under many different settings and assumptions. The performance of
the SDP relaxation for clustering Gaussian mixtures was first studied in [31],
using proof techniques from [17]. It was later shown that the SDP error decays
exponentially on the separation of the Gaussians [13, 14]. Other conic
relaxations of $k$-means clustering have been recently studied, for instance
[37, 36]. A (quite loose) linear programming relaxation of $k$-means was
proposed in [4], with mostly negative results. A significantly better LP
relaxation was recently introduced and analyzed in [11].

Another related line of work concerns efficiently certifying optimality of
solutions to data problems via convex relaxations and dual certificates. In
[6], Bandeira proposes leveraging a dual certificate to efficiently certify
optimality of solutions obtained with other (more efficient) methods. The goal
is to combine the theoretical guarantee from (slow) convex relaxations with
solutions provided by (fast) algorithms that may not have theoretical
guarantees. This idea has been used to provide a posteriori optimality
certificates in data science problems such as point cloud registration [45],
$k$-means clustering [19], and synchronization [38]. However, in order for
this method to succeed, the relaxation must be tight.

Sketch-and-solve methods provide a looser guarantee (not optimality, but
approximate optimality) that can work in broader contexts, including ones
where no convex relaxation is involved (for instance, numerical linear algebra
algorithms [43]). Here, we consider a setting where a data problem is relaxed
via convex relaxation. First, the original dataset is subsampled to a sketch,
then the convex problem is solved in the sketch, and finally a solution is
inferred for the entire dataset. The approximation guarantees from the convex
relaxation combined with regularity assumptions on the data (and how well the
sample can represent it) can be used to derive approximation guarantees for
the general approach. These ideas were used in [32, 2] in the context of graph
clustering, and in [46] in the context of clustering mixtures of Gaussians.
Here, we provide approximation guarantees for a sketch-and-solve approach for
general $k$-means clustering, but we also show that the approach provides an
efficient algorithm that computes a high-probability lower bound on the
$k$-means objective for any dataset, with no assumptions on the data nor how
it is generated. A similar observation was made by two of the authors in a
preprint [30].

###  1.3 Roadmap

The following section contains some preliminaries and introduces notation for
the remainder of the paper. Next, Section 3 introduces a sketch-and-solve
algorithm that determines the optimal clustering provided the clusters are
appropriately separated. Section 4 then introduces sketch-and-solve algorithms
to compute lower bounds on the optimal $k$-means value of a given dataset, and
we prove that these lower bounds are nearly sharp for data drawn from Gaussian
mixtures. In Section 5, we illustrate the quality of these bounds on real-
world datasets. We discuss opportunities for future work in Section 6, and our
main results are proved in Sections 7 and 8.

##  2 Preliminaries and notation

We are interested in clustering $n$ points in $\mathbb{R}^{d}$ into $k$
clusters. Denoting the index set $[n]:=\\{1,\ldots,n\\}$, let $\Pi(n,k)$ be
the set of partitions of $[n]$ into $k$ nonempty sets, i.e., if
$\Gamma\in\Pi(n,k)$, then $|\Gamma|=k$, $\bigsqcup_{S\in\Gamma}S=[n]$, and
$|S|>1$ for each $S\in\Gamma$. Given a tuple $X:=\\{x_{i}\\}_{i\in[n]}$ of
points in $\mathbb{R}^{d}$ and a nonempty set $S\subseteq[n]$ of indices, we
denote the corresponding centroid by

$c_{S}:=\frac{1}{|S|}\sum_{i\in S}x_{i}.$

Note that we suppress the dependence of $c_{S}$ on $X$ for simplicity. With
this notation, the (normalized) $k$-means problem is given by

$\text{minimize}\qquad\frac{1}{n}\sum_{S\in\Gamma}\sum_{i\in
S}\|x_{i}-c_{S}\|^{2}\qquad\text{subject to}\qquad\Gamma\in\Pi(n,k),$

and we denote the value of this program by $\operatorname{IP}(X,k)$. This
problem is trivial when $k=1$, and so we assume $k\geq 2$ in the sequel.

Lloyd’s algorithm is a popular approach to solve the $k$-means problem in
practice. This algorithm alternates between computing cluster centroids and
re-partitioning the data points according to the nearest centroid. These
iterations are inexpensive, costing only $O(kdn)$ operations each, and the
algorithm eventually converges to a (possibly sub-optimal) fixed point. The
value $V^{(0)}$ of the $k$-means++ random initialization enjoys the following
guarantee of approximate optimality (see Theorem 3.1 in [41]):

$\operatorname{IP}(X,k)\geq\mathbb{E}L,\qquad L:=\frac{V^{(0)}}{8(\log k+2)}.$
(1)

Since the $k$-means objective monotonically decreases with each iteration of
Lloyd’s algorithm, the $k$-means++ initialization ensures a $O(\log
k)$-competitive solution to the $k$-means problem on average.

As a theory-friendly alternative, one may instead solve the $k$-means problem
by relaxing to the Peng–Wei semidefinite program [35]. Let
$1_{S}\in\mathbb{R}^{n}$ denote the indicator vector of $S\subseteq[n]$, and
encode $\Gamma\in\Pi(n,k)$ with the matrix

$Z_{\Gamma}:=\sum_{S\in\Gamma}\frac{1}{|S|}1_{S}1_{S}^{\top}\in\mathbb{R}^{n\times
n}.$

Define $D_{X}\in\mathbb{R}^{n\times n}$ by
$(D_{X})_{ij}:=\|x_{i}-x_{j}\|^{2}$. A straightforward manipulation gives

$\sum_{S\in\Gamma}\sum_{i\in
S}\|x_{i}-c_{S}\|^{2}=\frac{1}{2}\sum_{S\in\Gamma}\frac{1}{|S|}\sum_{i\in
S}\sum_{j\in
S}\|x_{i}-x_{j}\|^{2}=\frac{1}{2}\operatorname{tr}(D_{X}Z_{\Gamma}).$

Thus, the $k$-means problem is equivalently given by

$\text{minimize}\qquad\frac{1}{2n}\operatorname{tr}(D_{X}Z_{\Gamma})\qquad\text{subject
to}\qquad\Gamma\in\Pi(n,k).$

Considering the containment

$\Big{\\{}Z_{\Gamma}:\Gamma\in\Pi(n,k)\Big{\\}}\subseteq\mathcal{Z}(n,k):=\Big{\\{}Z\in\mathbb{R}^{n\times
n}:Z1=1,~{}\operatorname{tr}Z=k,~{}Z\geq 0,~{}Z\succeq 0\Big{\\}},$

we obtain the (normalized) Peng–Wei semidefinite relaxation [35]:

$\text{minimize}\qquad\frac{1}{2n}\operatorname{tr}(D_{X}Z)\qquad\text{subject
to}\qquad Z\in\mathcal{Z}(n,k).$

We denote the value of this program by $\operatorname{SDP}(X,k)$. When the
relaxation is tight, the minimizer recovers the optimal $k$-means partition,
and otherwise, the value delivers a lower bound on the optimal $k$-means
value. For many real-world clustering instances, this SDP is too time
consuming to compute in practice. To resolve this issue, we introduce a few
sketch-and-solve approaches that work similarly well with far less run time.

##  3 Sketch-and-solve clustering

In this section, we introduce a sketch-and-solve approach to SDP-based
$k$-means clustering that works well provided the clusters are sufficiently
separated. Suppose we run the Peng–Wei SDP on a random subset of the data. If
the original clusters are well separated, then this random subset satisfies a
proximity condition from [27] with high probability, which in turn implies
that the Peng–Wei SDP recovers the desired clusters in this subset. Then we
can partition the full data set according to which of these cluster centroids
is closest. This approach is summarized in Algorithm 1. The main result of
this section (Theorem 2) is a theoretical guarantee for this algorithm.

We start by introducing the necessary notation to enunciate the proximity
condition from [27]. For $X\in(\mathbb{R}^{d})^{n}$ and $\Gamma\in\Pi(n,k)$,
and for each $S,T\in\Gamma$ with $S\neq T$, we define

$\alpha_{ST}:=\min_{i\in
S}\Big{\langle}x_{i}-\frac{c_{S}+c_{T}}{2},\frac{c_{S}-c_{T}}{\|c_{S}-c_{T}\|}\Big{\rangle},\quad\beta_{ST}:=\frac{1}{2}\Big{(}\Big{(}\frac{1}{|S|}+\frac{1}{|T|}\Big{)}\sum_{R\in\Gamma}\|X_{R}\|_{2\to
2}^{2}\Big{)}^{1/2}.$

Here, $X_{R}$ denotes the $d\times n$ matrix whose $i$th column equals zero
unless $i\in R$, in which case the column equals $x_{i}-c_{R}$. Notice that
$\alpha_{ST}$ captures how close $S$ is to the hyperplane that bisects the
centroids $c_{S}$ and $c_{T}$, while $\beta_{ST}$ scales some measure of
variance of the entire dataset by the sizes of $S$ and $T$. Intuitively,
$k$-means clustering is easier when clusters are well separated, e.g., when
the following quantity is positive:

$\operatorname{prox}(X,\Gamma):=\min_{\begin{subarray}{c}S,T\in\Gamma\\\ S\neq
T\end{subarray}}\Big{(}\alpha_{ST}-\beta_{ST}\Big{)}.$

######  Proposition 1 (Theorem 2 in [27]).

For each $X\in(\mathbb{R}^{d})^{n}$, there exists at most one
$\Gamma\in\Pi(n,k)$ for which $\operatorname{prox}(X,\Gamma)>0$. Furthermore,
if such $\Gamma$ exists, then $Z_{\Gamma}$ is the unique minimizer of the
Peng–Wei semidefinite relaxation for $X$, and therefore $\Gamma$ is the unique
minimizer of the $k$-means problem.

Next, we introduce notation necessary to enunciate the main result in this
section. Given $X\in(\mathbb{R}^{d})^{n}$ and $\Gamma\in\Pi(n,k)$, consider
the quantities

$\Delta:=\min_{\begin{subarray}{c}S,T\in\Gamma\\\ S\neq
T\end{subarray}}\|c_{S}-c_{T}\|,\qquad r:=\max_{S\in\Gamma}\max_{i\in
S}\|x_{i}-c_{S}\|.$

We will consider a shape parameter of $(X,\Gamma)$ that is determined by the
following quantities:

$d,\qquad\frac{\Delta}{r},\qquad\frac{\operatorname{prox}(X,\Gamma)}{r},\qquad
k,\qquad\frac{1}{n}\min_{S\in\Gamma}|S|,\qquad\frac{1}{n}\max_{S\in\Gamma}|S|.$

Notice that the first term above is completely determined by $X$, the next two
terms are determined by $X$ and $\Gamma$, and the last three terms are
completely determined by $\Gamma$. Since its dependence on $X$ factors through
$d$, $\frac{\Delta}{r}$, and $\frac{\operatorname{prox}(X,\Gamma)}{r}$, the
shape parameter is invariant to rotation, translation, and dilation.

Data: Points $X=\\{x_{i}\\}_{i\in[n]}\in(\mathbb{R}^{d})^{n}$, number of
clusters $k$, Bernoulli rate $p$

Result: Optimal $k$-means clustering $\Gamma\in\Pi(n,k)$

Draw $W\subseteq[n]$ according to a Bernoulli process with rate $p$ and put
$X^{\prime}:=\\{x_{i}\\}_{i\in W}$

Define $D_{X^{\prime}}\in\mathbb{R}^{W\times W}$ by
$(D_{X^{\prime}})_{ij}:=\|x_{i}-x_{j}\|^{2}$ for $i,j\in W$

Solve Peng–Wei semidefinite relaxation to find optimal
$\Gamma^{\prime}\in\Pi(|W|,k)$ for $X^{\prime}$

Compute centroids $\\{c_{S^{\prime}}\\}_{S^{\prime}\in\Gamma^{\prime}}$

Output $\Gamma\in\Pi(n,k)$ that partitions $X$ according to closest
$c_{S^{\prime}}$

Algorithm 1 Sketch-and-solve algorithm for Peng–Wei semidefinite relaxation

######  Theorem 2.

There exists an explicit shape parameter
$C\colon(\mathbb{R}^{d})^{n}\times\Pi(n,k)\to[0,\infty]$ for which the
following holds:

  * (a)

Suppose $X\in(\mathbb{R}^{d})^{n}$ and $\Gamma\in\Pi(n,k)$ satisfy
$\operatorname{prox}(X,\Gamma)>0$ and $r\leq\frac{\Delta}{2}$. Then
$C(X,\Gamma)<\infty$, and for $\epsilon\in(0,\frac{1}{2})$, Algorithm 1
exactly recovers $\Gamma$ from $X$ with probability $1-\epsilon$ provided

$\mathbb{E}|W|>C(X,\Gamma)\cdot\log(1/\epsilon).$

(Here, the randomness is in the Bernoulli process in Algorithm 1.)

  * (b)

$C(X,\Gamma)$ is directly related to $d$, $k$, and
$\frac{1}{n}\max_{S\in\Gamma}|S|$, and inversely related to
$\frac{\Delta}{r}$, $\frac{\operatorname{prox}(X,\Gamma)}{r}$, and
$\frac{1}{n}\min_{S\in\Gamma}|S|$.

The proof idea for Theorem 2 is simple: identify conditions under which the
sketched data satisfy the proximity condition with probability $1-\epsilon$.
Notably, the complexity of the SDP step of Algorithm 1 depends on $|W|$, which
in turn scales according to the shape of the data rather than its size.
Indeed, when the clusters are more separated, the shape parameter is smaller,
and so the sketch size can be taken to be smaller by Theorem 2.

For the sake of illustration, we test the performance of Algorithm 1 on data
drawn according to the stochastic ball model. Fix $\\{\mu_{a}\\}_{a\in[k]}$ in
$\mathbb{R}^{d}$, and for each $a$, draw $\\{g_{a,i}\\}_{i\in[m]}$
independently according to some rotation-invariant probability distribution
$\mathcal{D}$ supported on the origin-centered unit ball in $\mathbb{R}^{d}$.
Then we write $X\sim\mathsf{SBM}(\mathcal{D},\\{\mu_{a}\\}_{a\in[k]},m)$ to
denote the data $X=\\{\mu_{a}+g_{a,i}\\}_{a\in[k],i\in[m]}$. In words, we draw
$m$ points from each of $k$ unit balls centered at the $\mu_{a}$’s, meaning
$n=km$. We will focus on the case where $d=2$, $k=2$, and $\mathcal{D}$ is the
uniform distribution on the unit ball.

First, we consider the behavior of our sketch-and-solve approach as
$n\to\infty$. (Indeed, we have the luxury of considering such behavior since
the performance of our approach does not depend on $n$, but rather on the
shape of the data.) In the limit, one may show that

$\Delta=\|\mu_{1}-\mu_{2}\|,\qquad
r=1,\qquad\operatorname{prox}(X,\Gamma)=\frac{\Delta-3}{2},\qquad\frac{1}{n}\min_{S\in\Gamma}|S|=\frac{1}{n}\max_{S\in\Gamma}|S|=\frac{1}{2}.$

By Theorem 2(b), it follows that the shape parameter $C(X,\Gamma)$ is
inversely related to $\Delta$, as one might expect. Figure 1(left) illustrates
that Algorithm 1 exactly recovers the planted clustering of the entire balls
provided they are appropriately separated and the sketch size $|W|$ isn’t too
small. Overall, the behavior reported in Figure 1(left) qualitatively matches
the prediction in Theorem 2.

Figure 1:  (left) For each $\Delta\in\\{2,2.1,\ldots,4\\}$ and
$|W|\in\\{2,4,\ldots,30\\}$, perform the following experiment $500$ times:
Consider two unit balls in $\mathbb{R}^{2}$ whose centers are separated by
$\Delta$. Draw $|W|$ points uniformly from the union of the balls, solve the
Peng–Wei SDP for this sketch in CVX [16], compute the cluster centroids, and
partition the union of balls according to the nearest centroid. Plot the
proportion of these $500$ trials for which the planted clustering was exactly
recovered. (right) For each $n\in\\{2^{3},2^{4},\ldots,2^{25}\\}$, perform the
following experiment $10$ times: Draw $n$ points in $\mathbb{R}^{2}$ according
to the uniformly distributed stochastic ball model with separation $\Delta=3$.
Run Algorithm 1 with $k=2$ and $p=\min\\{\frac{10}{n},1\\}$ so that
$\mathbb{E}|W|=\min\\{10,n\\}$. For comparison, run the Peng–Wei SDP in CVX
and MATLAB’s built-in $k$-means++ algorithm (both with $k=2$), and plot the
average runtimes.

Next, we compare the runtime of our method to traditional (non-sketching)
methods in Figure 1(right). For this, we consider the same stochastic ball
model with separation $\Delta=3$. As mentioned earlier, one cannot run the
Peng–Wei SDP directly on a large dataset in a reasonable amount of time. Also,
each iteration of $k$-means++ takes linear runtime. Meanwhile, the bulk of our
sketch-and-solve approach has a runtime that scales with the shape of the data
(i.e., it’s independent of the size of the data). Of course, the final
clustering step of our approach requires a single pass of the data, which
explains the slight increase in runtime for larger datasets.

##  4 High-confidence lower bounds

The previous section focused on a sketch-and-solve approach that recovers the
optimal $k$-means clustering whenever the data exhibits a sufficiently nice
shape (in some quantifiable sense). However, the Peng–Wei semidefinite
relaxation fails to be tight when the planted clusters are not well separated,
so we should not expect exact recovery in general. Regardless, the SDP always
delivers a lower bound on the optimal $k$-means value, which can be useful in
practice (e.g., when deciding whether to re-initialize Lloyd’s algorithm again
to find an even better clustering).

In this section, we introduce sketch-and-solve approaches for computing such
bounds in a reasonable amount of time. Our approaches are inspired by the
following simple result:

######  Lemma 3.

Consider any sequence $X:=\\{x_{i}\\}_{i\in[n]}$ in $\mathbb{R}^{d}$, draw
indices $\\{i_{j}\\}_{j\in[s]}$ uniformly from $[n]$ with (or without)
replacement, and put $Y:=\\{x_{i_{j}}\\}_{j\in[s]}$. Then

$\mathbb{E}\operatorname{SDP}(Y,k)\leq\mathbb{E}\operatorname{IP}(Y,k)\leq\operatorname{IP}(X,k).$

###### Proof.

The first inequality follows from relaxation. For the second inequality,
select $\Gamma\in\arg\operatorname{IP}(X,k)$, consider the set-valued function
$\sigma\colon[n]\to\Gamma$ that satisfies $i\in\sigma(i)$ for every $i\in[n]$,
and for each $j\in[s]$, denote the random variable

$E_{j}:=\bigg{\|}x_{i_{j}}-\frac{1}{|\sigma(i_{j})|}\sum_{i\in\sigma(i_{j})}x_{i}\bigg{\|}^{2}.$

The random variables $\\{E_{j}\\}_{j\in[s]}$ have a common distribution with
expectation $\operatorname{IP}(X,k)$, though they are dependent if the indices
$\\{i_{j}\\}_{j\in[s]}$ are drawn without replacement. Considering the random
function $f\colon[s]\to[n]$ defined by $f(j):=i_{j}$, we have

$\operatorname{IP}(Y,k)\leq\frac{1}{s}\sum_{S\in\Gamma}\sum_{j\in
f^{-1}(S)}\bigg{\|}x_{i_{j}}-\frac{1}{|f^{-1}(S)|}\sum_{j^{\prime}\in
f^{-1}(S)}x_{i_{j^{\prime}}}\bigg{\|}^{2}\leq\frac{1}{s}\sum_{j\in[s]}E_{j},$

where the second inequality follows from the fact that the centroid of a tuple
of points minimizes the sum of squared distances from those points. The result
follows by taking the expectation of both sides. ∎

By Lemma 3, we can lower bound the $k$-means value $\operatorname{IP}(X,k)$ by
estimating the expected Peng–Wei value $\mathbb{E}\operatorname{SDP}(Y,k)$ of
a random sketch $Y$. To this end, given an error rate $\epsilon>0$ and
$\ell\in\mathbb{N}$ draws of the random sketch, we may leverage concentration
inequalities to compute a random variable $B$ that is smaller than
$\mathbb{E}\operatorname{SDP}(Y,k)$ with probability $\geq 1-\epsilon$. We
provide two such random variables in Algorithms 3 and 4. Notably, the random
variable $B_{H}$ computed by Algorithm 3 is consistent in the sense that
$B_{H}$ converges in probability to $\mathbb{E}\operatorname{SDP}(Y,k)$ as
$\ell\to\infty$. Meanwhile, the random variable $B_{M}$ computed by Algorithm
4 is not consistent, but as we show in Section 5, it empirically outperforms
$B_{H}$ when $\ell$ is small. Theorems 4(a) and 5(a) give that these random
variables indeed act as lower bounds with probability $\geq 1-\epsilon$.

Of course, we would like these lower bounds to be as sharp as possible. To
evaluate how close they are to the desired $k$-means value, we consider the
Gaussian mixture model. Given means $\mu_{1},\ldots,\mu_{k}\in\mathbb{R}^{d}$,
covariances $\Sigma_{1},\ldots,\Sigma_{k}\in\mathbb{R}^{d\times d}$, and a
probability distribution $p$ over the index set $[k]$, the random vector
$x\sim\mathsf{GMM}(\\{(\mu_{t},\Sigma_{t},p(t))\\}_{t\in[k]})$ is obtained by
first drawing $T$ from $[k]$ with distribution $p$, and then drawing $x$ from
the Gaussian $\mathsf{N}(\mu_{T},\Sigma_{T})$. The Gaussian mixture model can
be thought of as a “noisy” version of the stochastic ball model. By part (b)
of the following results, our random lower bounds are nearly sharp provided
$d\gg k$, even when there is no separation between the Gaussian means. See
Section 8 for the proofs of these parts.

Data: Points $X=\\{x_{i}\\}_{i\in[n]}\in(\mathbb{R}^{d})^{n}$, number of
clusters $k$

Result: Indices $\\{i_{j}\\}_{j\in[k]}\in[n]^{k}$ of well-separated points

Put $i_{1}:=1$, and iteratively select
$i_{t+1}\in\arg\max_{i\in[n]}\min_{j\in[t]}\|x_{i}-x_{i_{j}}\|$

Algorithm 2 Deterministic $k$-means$++$ initialization

Data: Points $X=\\{x_{i}\\}_{i\in[n]}\in(\mathbb{R}^{d})^{n}$, number of
clusters $k$, sketch size $s$, number of trials $\ell$, error rate $\epsilon$

Result: Random variable $B_{H}$ such that $\operatorname{IP}(X,k)\geq B_{H}$
with probability $\geq 1-\epsilon$

Run deterministic $k$-means$++$ initialization and put
$b:=\max_{i\in[n]}\min_{j\in[k]}\|x_{i}-x_{i_{j}}\|^{2}$

Draw $\\{Y_{i}\\}_{i\in[\ell]}$ independently at random, with each $Y_{i}$
denoting $s$ points drawn uniformly from $X$ with replacement

Output
$B_{H}:=\frac{1}{\ell}\sum_{i\in[\ell]}\operatorname{SDP}(Y_{i},k)-(\frac{b^{2}}{2\ell}\log(\frac{1}{\epsilon}))^{1/2}$

Algorithm 3 Hoeffding Monte Carlo $k$-means lower bound

Data: Points $X=\\{x_{i}\\}_{i\in[n]}\in(\mathbb{R}^{d})^{n}$, number of
clusters $k$, sketch size $s$, number of trials $\ell$, error rate $\epsilon$

Result: Random variable $B_{M}$ such that $\operatorname{IP}(X,k)\geq B_{M}$
with probability $\geq 1-\epsilon$

Draw $\\{Y_{i}\\}_{i\in[\ell]}$ independently at random, with each $Y_{i}$
denoting $s$ points drawn uniformly from $X$ without replacement

Output $B_{M}:=\epsilon^{1/\ell}\min_{i\in[\ell]}\operatorname{SDP}(Y_{i},k)$

Algorithm 4 Markov Monte Carlo $k$-means lower bound

######  Theorem 4 (Performance guarantee for Algorithm 3).

  * (a)

Consider any $X:=\\{x_{i}\\}_{i\in[n]}$ in $\mathbb{R}^{d}$, any
$k,s,\ell\in\mathbb{N}$, and $\epsilon>0$, and compute the random variable
$B_{H}$ in Algorithm 3. Then

$\operatorname{IP}(X,k)\geq B_{H}$

with probability $\geq 1-\epsilon$. (Here, the probability is on Algorithm 3.)

  * (b)

Consider any $\mu_{1},\ldots,\mu_{k}\in\mathbb{R}^{d}$, draw the points
$X:=\\{x_{i}\\}_{i\in[n]}$ independently with distribution
$\mathsf{GMM}(\\{(\mu_{t},I_{d},\frac{1}{k})\\}_{t\in[k]})$, take any
$s,\ell\in\mathbb{N}$, and $\epsilon>0$, and compute the random variable
$B_{H}$ in Algorithm 3. Then

$B_{H}\geq\frac{d-6k-2}{d+1}\cdot\operatorname{IP}(X,k)$

with probability $\geq 1-\frac{1}{n}-e^{-\Omega(n/(d+\log n)^{2})}-\epsilon$
provided

$s\geq 15d\log d,\qquad\ell\geq 128(d+3\log n)^{2}\log(1/\epsilon).$

(Here, the probability is on both $X$ and Algorithm 3.)

###### Proof of Theorem 4(a).

Take $\\{x_{i_{j}}\\}_{j\in[k]}$ from Algorithm 2, consider any tuple
$\\{i^{\prime}_{t}\\}_{t\in[s]}$ of indices in $[n]$, and put
$Y:=\\{x_{i^{\prime}_{t}}\\}_{t\in[s]}$. Then

$\displaystyle\operatorname{SDP}(Y,k)$
$\displaystyle\leq\operatorname{IP}(Y,k)=\min_{\mu_{1},\ldots,\mu_{k}\in\mathbb{R}^{d}}\frac{1}{s}\sum_{t\in[s]}\min_{j\in[k]}\|x_{i^{\prime}_{t}}-\mu_{j}\|^{2}$
$\displaystyle\leq\frac{1}{s}\sum_{t\in[s]}\min_{j\in[k]}\|x_{i^{\prime}_{t}}-x_{i_{j}}\|^{2}\leq\max_{t\in[s]}\min_{j\in[k]}\|x_{i^{\prime}_{t}}-x_{i_{j}}\|^{2}\leq\max_{i\in[n]}\min_{j\in[k]}\|x_{i}-x_{i_{j}}\|^{2}.$

It follows that $\operatorname{SDP}(Y_{i},k)\leq b$ almost surely for each
$i\in[\ell]$. The result then follows from Lemma 3 and Hoeffding’s inequality:

$\mathbb{P}\\{B_{H}>\operatorname{IP}(X,k)\\}\leq\mathbb{P}\bigg{\\{}\frac{1}{\ell}\sum_{i\in[\ell]}\operatorname{SDP}(Y_{i},k)-\Big{(}\tfrac{b^{2}}{2\ell}\log(\tfrac{1}{\epsilon})\Big{)}^{1/2}>\mathbb{E}\operatorname{SDP}(Y_{1},k)\bigg{\\}}\leq\epsilon.\qed$

######  Theorem 5 (Performance guarantee for Algorithm 4).

  * (a)

Consider any $X:=\\{x_{i}\\}_{i\in[n]}$ in $\mathbb{R}^{d}$, any
$k,s,\ell\in\mathbb{N}$, and $\epsilon>0$, and compute the random variable
$B_{M}$ in Algorithm 4. Then

$\operatorname{IP}(X,k)\geq B_{M}$

with probability $\geq 1-\epsilon$. (Here, the probability is on Algorithm 4.)

  * (b)

Consider any $\mu_{1},\ldots,\mu_{k}\in\mathbb{R}^{d}$, draw the points
$X:=\\{x_{i}\\}_{i\in[n]}$ independently with distribution
$\mathsf{GMM}(\\{(\mu_{t},I_{d},\frac{1}{k})\\}_{t\in[k]})$, take any
$s,\ell\in\mathbb{N}$, and $\epsilon>0$, and compute the random variable
$B_{M}$ in Algorithm 4. Then

$B_{M}\geq\frac{d-3k-2}{d+1}\cdot\operatorname{IP}(X,k)$

with probability $\geq 1-e^{-s/(8d)}-2e^{-s/54}-e^{-n/(16d)}$ provided

$s\geq 54d\log\ell,\qquad\ell\geq d\log(1/\epsilon),\qquad\ell\geq 3.$

(Here, the probability is on both $X$ and Algorithm 4.)

###### Proof of Theorem 5(a).

By Lemma 3, we have
$\operatorname{IP}(X,k)\geq\mathbb{E}\operatorname{SDP}(Y_{1},k)$. It follows
that

$\mathbb{P}\\{B_{M}>\operatorname{IP}(X,k)\\}\leq\mathbb{P}\\{B_{M}>\mathbb{E}\operatorname{SDP}(Y_{1},k)\\}=\mathbb{P}\\{\epsilon^{1/\ell}\operatorname{SDP}(Y_{1},k)>\mathbb{E}\operatorname{SDP}(Y_{1},k)\\}^{\ell},$

where the last step applies the fact that
$\\{\operatorname{SDP}(Y_{i},k)\\}_{i\in[k]}$ are independent with identical
distribution. Markov’s inequality further bounds this upper bound by
$\epsilon$. ∎

The proof of Theorem 4(b) makes use of the following approximation ratio
afforded by deterministic $k$-means++ initialization (Algorithm 2), which may
be of independent interest.

######  Lemma 6.

Given points $X:=\\{x_{i}\\}_{i\in[n]}$ in $\mathbb{R}^{d}$, consider the
problem of minimizing the function $f\colon(\mathbb{R}^{d})^{k}\to\mathbb{R}$
defined by

$f(\\{z_{j}\\}_{j\in[k]}):=\max_{i\in[n]}\min_{j\in[k]}\|x_{i}-z_{j}\|.$

The output $\\{i_{j}\\}_{j\in[k]}$ of deterministic $k$-means$++$
initialization (Algorithm 2) indexes an approximate solution to this
optimization problem with approximation ratio $2$:

$f(\\{x_{i_{j}}\\}_{j\in[k]})\leq 2\inf(f).$

###### Proof.

First, we establish that $f$ has a global minimizer. Taking $z_{j}=0$ for
every $j$ gives $f(\\{z_{j}\\}_{j\in[k]})=\max_{i\in[n]}\|x_{i}\|=:R$. Let
$w\in\mathbb{R}^{d}$ denote a point for which $\min_{i\in[n]}\|x_{i}-w\|=R$.
For any input $\\{z_{j}\\}_{j\in[k]}$ such that $f(\\{z_{j}\\}_{j\in[k]})\leq
R$, then for each $j$, either $z_{j}$ is within $R$ of the nearest $x_{i}$, or
the function value is not changed by sending $z_{j}$ to $w$. Thus, we may
restrict the input space to the set of $\\{z_{j}\\}_{j\in[k]}$ for which every
$z_{j}$ is within $R$ of the nearest $x_{i}$. Since this set is compact and
$f$ is continuous, the extreme value theorem guarantees a global minimizer.

Let $\\{z_{j}^{\star}\\}_{j\in[k]}$ denote a global minimizer of $f$, and put
$r:=\min(f)$. Then every $x_{i}$ is within $r$ of some $z_{j}^{\star}$. Run
deterministic $k$-means$++$ initialization, and consider the quantities

$r_{t}:=\max_{i\in[n]}\min_{j\in[t]}\|x_{i}-x_{i_{j}}\|$

for $t\in[k]$. Observe that $r_{1}\geq\cdots\geq
r_{k}=f(\\{x_{i_{j}}\\}_{j\in[k]})$. We wish to demonstrate $r_{k}\leq 2r$. If
$r_{k-1}\leq 2r$, then we are done. Otherwise, we have $r_{t}>2r$ for every
$t<k$, and so the pairwise distances between $\\{x_{i_{j}}\\}_{j\in[k]}$ are
all strictly greater than $2r$. Since each $x_{i_{j}}$ is within $r$ of a
point in $\\{z_{j}^{\star}\\}_{j\in[k]}$, the triangle inequality gives that
no two $x_{i_{j}}$’s are within $r$ of the same point in
$\\{z_{j}^{\star}\\}_{j\in[k]}$. Since every $x_{i}$ is within $r$ of some
$z_{j}^{\star}$, which in turn is within $r$ of some $x_{i_{j}}$, the triangle
inequality further gives that every $x_{i}$ is within $2r$ of some
$x_{i_{j}}$. This establishes $r_{k}\leq 2r$, as desired. ∎

##  5 Numerical experiments

In this section, we use both synthetic and real-world data to test the
performance of our high-confidence lower bounds $B_{H}$ (Algorithm 3) and
$B_{M}$ (Algorithm 4).

###  5.1 Implementation details

When implementing our algorithms, we discovered a few modifications that
delivered substantial improvements in practice.

Correcting the numerical dual certificate. SDPNAL+ [39] is a fast SDP solver
that was specifically designed for SDPs with entrywise nonnegativity
constraints, such as the Peng–Wei SDP. For this reason, we used SDPNAL+ to
solve our sketched SDP, but then it delivered a dual certificate that was not
exactly dual feasible. To be explicit, recall the normalized Peng–Wei SDP for
a sketch $Y$ of size $s$:

$\text{minimize}\qquad\frac{1}{2s}\operatorname{tr}(D_{Y}Z)\qquad\text{subject
to}\qquad Z\in\mathbb{R}^{s\times s},~{}\mathcal{A}(Z)=b,~{}Z\geq
0,~{}Z\succeq 0,$

where $\mathcal{A}\colon\mathbb{R}^{s\times s}\to\mathbb{R}^{s+1}$ and
$b\in\mathbb{R}^{s+1}$ are defined by

$\displaystyle\mathcal{A}(Z):=\left[\begin{matrix}\operatorname{tr}(Z)\\\
Z1\end{matrix}\right],\qquad b:=\left[\begin{matrix}k\\\
1\end{matrix}\right].$

The corresponding dual program is

$\text{minimize}\quad\frac{1}{2s}b^{\top}y\quad\text{subject to}\quad
y\in\mathbb{R}^{s+1},~{}P,S\in\mathbb{R}^{s\times
s},~{}\mathcal{A}^{*}(y)+P+S=D_{Y},~{}P\geq 0,~{}S\succeq 0.$

In practice, SDPNAL+ returns $(y_{0},P_{0},S_{0})$ such that
$\mathcal{A}^{*}(y)+P+S\neq D_{Y}$. To correct this, we fix $S=S_{0}$ in the
dual program, thereby restricting to an easy-to-solve linear program, which we
then solve using CVX. The result is a point $(y_{1},P_{1},S_{0})$ that is
feasible in the original dual program, and so weak duality implies that
$\frac{1}{2s}b^{\top}y_{1}$ is a lower bound, as desired.

Truncating the distance matrix. For some datasets (e.g., INTRUSION), the
entries of $D_{X}$ might vary widely, and we observe that SDP solvers fail to
converge in such cases. This can be resolved by truncating. Specifically, we
apply the function $t\mapsto\min\\{t,10^{8}\\}$ to each entry of $D_{X}$
before solving the SDP. Not only does this make the problem solvable in
practice, the optimal value of the truncated problem is a lower bound on the
original optimal value, and so we can use the result.

Truncating the sketched SDP value. When we run Algorithm 3 in practice, we
find that $b$ is often so large that $B_{H}$ is negative (and therefore
useless). Recall that $b$ plays the role of an almost sure upper bound on
$\operatorname{SDP}(Y,k)$. We can replicate this behavior by selecting a
threshold $u>0$ and replacing $\operatorname{SDP}(Y,k)$ with the truncated
version $\min\\{\operatorname{SDP}(Y,k),u\\}$. This leads to the following
modification:

$B_{H}=\frac{1}{\ell}\sum_{i\in[\ell]}\min\Big{\\{}\operatorname{SDP}(Y_{i},k),u\Big{\\}}-\Big{(}\tfrac{u^{2}}{2\ell}\log(\tfrac{1}{\epsilon})\Big{)}^{1/2}.$

Following the proof of Theorem 4(a), we have that

$\displaystyle\mathbb{P}\\{B_{H}>\operatorname{IP}(X,k)\\}$
$\displaystyle\quad\leq\mathbb{P}\bigg{\\{}\frac{1}{\ell}\sum_{i\in[\ell]}\min\Big{\\{}\operatorname{SDP}(Y_{i},k),u\Big{\\}}-\Big{(}\tfrac{u^{2}}{2\ell}\log(\tfrac{1}{\epsilon})\Big{)}^{1/2}>\mathbb{E}\min\Big{\\{}\operatorname{SDP}(Y_{1},k),u\Big{\\}}\bigg{\\}}\leq\epsilon,$

and so it still delivers a high-confidence lower bound. In practice, we take
$u$ to be the $k$-means value of the clustering given by $k$-means++, and we
observe that the truncation usually occurs for only a few of the $\ell$
sketches.

Other implementation details. When solving our SDPs with SDPNAL+, we warm
start with a block diagonal primal matrix $Z$ that we obtain by running
$k$-means++ over the sketched data. Also, we solve SDPs to low precision by
default, but in cases where this only takes a few iterations, we solve the SDP
again to a higher precision.

###  5.2 Datasets and parameters

We test our algorithms on five datasets. MNIST, consisting of $60000$ points
in $784$ dimensions, is the training set images of the MNIST database of
handwritten digits [26]. NORM-10, consisting of $10000$ points in $5$
dimensions, is a synthetic dataset drawn according to a mixture of $10$
Gaussians $\mathsf{GMM}(\\{(\mu_{t},I_{5},\frac{1}{10})\\}_{t\in[10]})$ with
the centers $\\{\mu_{t}\\}_{t\in[10]}$ drawn uniformly in a $5$-dimensional
hypercube of side length $500$. Similarly NORM-25 consisting of $10000$ points
in $15$ dimensions, is a synthetic dataset drawn according to a mixture of
$25$ Gaussians $\mathsf{GMM}(\\{(\mu_{t},I_{15},\frac{1}{25})\\}_{t\in[25]})$
with the centers $\\{\mu_{t}\\}_{t\in[25]}$ drawn uniformly in a
$15$-dimensional hypercube of side length $500$. CLOUD, consisting of $1024$
points in $10$ dimensions, is Philippe Collard’s cloud cover data [8] in the
UC Irvine Machine Learning Repository. INTRUSION, consisting of $494021$
points in $34$ dimensions, is a $10\%$ subset data of continuous features
(with symbolic features removed) for network intrusion detection used in KDD
Cup 1999 competition [22], also from the UC Irvine Machine Learning
Repository. Notice that NORM-25, CLOUD, and INTRUSION are the same datasets
used in [41] to show the performance of $k$-means++ algorithm.

For each dataset, we run our algorithms in various settings. First, we run our
algorithms with a small number of trials, namely, $\ell=30$. We test $k=10$
for MNIST and $k\in\\{10,25,50\\}$ for each of the other four datasets. We
expect the Hoeffding lower bound $B_{H}$ to perform better than the Markov
lower bound $B_{M}$ when $\ell$ is larger, so we also consider the setting
where $\ell=1000$. In this setting, we took $k=25$ for NORM-25 and $k=10$ for
the other four datasets. For all of these experiments, we use sketch size
$s=300$ and error rate $\epsilon=0.01$.

###  5.3 Results

Before running our algorithms, we run the $k$-means++ algorithm $\ell$ times
on the full dataset and record the smallest $k$-means value, denoted by $\min
v_{i}$. (This number of trials equals the number of random sketch trials for
simplicity.) We take this $k$-means value to be truncation level $u$ in our
SDP truncation technique for $B_{H}$. Note that $\min v_{i}$ is an upper bound
on the optimal $k$-means objective.

For comparison, we use (1) to get high-confidence lower bounds similar to
$B_{H}$ (with value truncation) and $B_{M}$:

$L_{H}:=\frac{1}{\ell}\sum_{i\in[\ell]}\min\\{L_{i},u\\}-\Big{(}\tfrac{u^{2}}{2\ell}\log(\tfrac{1}{\epsilon})\Big{)}^{1/2},\qquad
L_{M}:=\epsilon^{1/\ell}\min_{i\in[\ell]}L_{i},$

where $\ell$, $u$, and $\epsilon$ are the same as for $B_{H}$ and $B_{M}$, and
each $L_{i}$ is an independent draw of $L$. By the proofs of Theorems 4(a) and
5(a), the random variables $L_{H}$ and $L_{M}$ are indeed lower bounds on
$\operatorname{IP}(X,k)$ with probability $1-\epsilon$. In Tables 1 and 2, we
include the average of $\\{L_{i}\\}_{i\in[\ell]}$ for comparison (denoted by
$\operatorname{avg}L_{i}$), which can be seen as a lower bound without a
probability guarantee.

Our results with $\ell=30$ are in Table 1, while our results with $\ell=1000$
are in Table 2. The bold numbers give the best lower bound among $L_{H}$,
$L_{M}$, $B_{H}$, and $B_{M}$. In all datasets except INTRUSION, our lower
bound (the best between $B_{H}$ and $B_{M}$) is at least $10$ times better
than $L_{H}$ and $L_{M}$ (and is still a significant improvement over
$\operatorname{avg}L_{i}$). However, our lower bound is worse than the
$k$-means++ lower bound for the INTRUSION dataset. This is because INTRUSION
is quite unbalanced: most of the data is concentrated at the same point, and
there are a few very distant outliers. Indeed, when we uniformly sketch down
to a small subset, most of the outliers may not be selected, and so even the
optimal $k$-means value of the subset $\operatorname{IP}(Y,k)$ is much smaller
than $\operatorname{IP}(X,k)$. (This can be demonstrated by running
$k$-means++ on $Y$ for INTRUSION.) As such, the bad performance on INTRUSION
data is due to sketching rather than relaxing. As expected, $B_{M}$ is
slightly better when $\ell$ is small, while $B_{H}$ slightly better when
$\ell$ is large. In addition, we observe that $L_{H}$ is consistently negative
(as is $B_{H}$ for INTRUSION).

Tables 1 and 2 also display three types of runtime. $T_{\mathrm{init}}$ is the
time for $\ell$ repeated initializations of $k$-means++, which is used to get
$\\{L_{i}\\}_{i\in[\ell]}$. $T_{\mathrm{k++}}$ is the time to compute $\min
v_{i}$ (i.e., the threshold $u$). $T_{\mathrm{SDP}}$ is the time to compute
$\ell$ randomly sketched SDPs. Based on the definitions of these lower bounds,
the runtime for $L_{H}$ is $T_{\mathrm{init}}+T_{\mathrm{k++}}$, the runtime
for $L_{M}$ is $T_{\mathrm{init}}$, the runtime for $B_{H}$ is
$T_{\mathrm{SDP}}+T_{\mathrm{k++}}$, and the runtime for $B_{M}$ is
$T_{\mathrm{SDP}}$. While the SDP approach takes longer than $k$-means++ in
these examples, the SDP approach usually delivers a much better lower bound.

Table 1: High-confidence $k$-means lower bounds with $\ell=30$. Dataset | $k$ | $\min v_{i}$ | $\operatorname{avg}L_{i}$ | $L_{H}$ | $L_{M}$ | $B_{H}$ | $B_{M}$ | $T_{\mathrm{init}}$ | $T_{\mathrm{k++}}$ | $T_{\mathrm{SDP}}$
---|---|---|---|---|---|---|---|---|---|---
MNIST | 10 | 3.92e1 | 1.26e0 | -9.59e0 | 1.06e0 | 2.56e1 | 2.96e1 | 1.71e1 | 4.18e2 | 2.47e2
NORM-10 | 10 | 4.97e0 | 1.10e0 | -1.07e0 | 1.24e-1 | 3.42e0 | 3.72e0 | 2.24e-1 | 1.48e-1 | 2.92e1
NORM-10 | 25 | 4.05e0 | 1.02e-1 | -1.01e0 | 8.63e-2 | 1.43e0 | 2.11e0 | 4.00e-1 | 1.97e0 | 3.98e2
NORM-10 | 50 | 3.18e0 | 7.58e-2 | -8.05e-1 | 6.33e-2 | 5.24e-1 | 1.12e0 | 6.30e-1 | 3.71e0 | 1.93e2
NORM-25 | 10 | 1.18e5 | 3.94e3 | -2.90e4 | 3.03e3 | 6.99e4 | 8.08e4 | 2.58e-1 | 1.95e-1 | 1.67e2
NORM-25 | 25 | 1.50e1 | 6.37e0 | -3.31e0 | 3.08e-1 | 9.61e0 | 1.15e1 | 4.42e-1 | 3.63e-1 | 3.44e1
NORM-25 | 50 | 1.41e1 | 3.04e-1 | -3.61e0 | 2.60e-1 | 5.23e0 | 7.48e0 | 6.88e-1 | 2.50e0 | 4.64e1
CLOUD | 10 | 5.62e3 | 2.41e2 | -1.31e3 | 1.62e2 | 2.70e3 | 3.06e3 | 1.01e-1 | 1.75e-1 | 1.46e2
CLOUD | 25 | 1.94e3 | 6.31e1 | -4.75e2 | 4.91e1 | 8.24e2 | 9.43e2 | 1.54e-1 | 1.90e-1 | 3.48e2
CLOUD | 50 | 1.09e3 | 2.99e1 | -2.72e2 | 2.33e1 | 2.57e2 | 4.54e2 | 2.19e-1 | 3.67e-1 | 4.99e2
INTRUSION | 10 | 2.36e7 | 1.00e6 | -5.55e6 | 6.03e5 | -6.43e6 | 2.93e4 | 9.46e0 | 3.76e1 | 3.43e2
INTRUSION | 25 | 2.19e6 | 7.64e4 | -5.31e5 | 5.31e4 | -6.03e5 | 1.67e3 | 1.86e1 | 1.04e2 | 4.58e2
INTRUSION | 50 | 4.52e5 | 1.36e4 | -1.11e5 | 9.27e3 | -1.25e5 | 4.56e1 | 3.43e1 | 1.51e2 | 2.16e3
Table 2: High-confidence $k$-means lower bounds with $\ell=1000$. Dataset | $k$ | $\min v_{i}$ | $\operatorname{avg}L_{i}$ | $L_{H}$ | $L_{M}$ | $B_{H}$ | $B_{M}$ | $T_{\mathrm{init}}$ | $T_{\mathrm{k++}}$ | $T_{\mathrm{SDP}}$
---|---|---|---|---|---|---|---|---|---|---
MNIST | 10 | 3.92e1 | 1.26e0 | -6.11e-1 | 1.21e0 | 3.44e1 | 3.39e1 | 5.79e2 | 1.45e4 | 8.08e3
NORM-10 | 10 | 5.02e0 | 4.05e-1 | -8.05e-2 | 1.45e-1 | 4.59e0 | 4.18e0 | 8.81e0 | 5.02e0 | 1.14e3
NORM-25 | 25 | 1.49e1 | 1.73e0 | -1.84e-1 | 3.56e-1 | 1.29e1 | 1.27e1 | 1.56e1 | 1.14e1 | 1.26e3
CLOUD | 10 | 5.62e3 | 2.31e2 | -3.89e1 | 1.72e2 | 4.06e3 | 3.14e3 | 4.05e0 | 6.10e0 | 6.27e3
INTRUSION | 10 | 2.34e7 | 9.58e5 | -1.66e5 | 6.80e5 | -1.00e6 | 1.81e4 | 3.15e2 | 2.44e3 | 2.65e4

##  6 Discussion

In this paper, we introduced sketch-and-solve approaches to $k$-means
clustering with semidefinite programming. In particular, we exactly recover
the optimal clustering with a single sketch provided the data is well
separated, and we compute a high-confidence lower bound on the $k$-means value
from multiple sketches when the data is not well separated. For future work,
one might attempt to use multiple sketches to find the optimal clustering when
the data only exhibits some separation. Next, our high-confidence lower bounds
perform poorly for unbalanced data like INTRUSION, and we suspect this stems
from our uniform sampling approach. Presumably, an importance sampling–based
alternative would perform better in such settings. Finally, the main idea of
our high-confidence lower bounds is that
$\mathbb{E}\operatorname{IP}(Y,k)\leq\operatorname{IP}(X,k)$ (see Lemma 3). It
would be interesting to show a similar relationship for SDP, namely,
$\mathbb{E}\operatorname{SDP}(Y,k)\leq\operatorname{SDP}(X,k)$. While this
bound holds empirically, we do not have a proof. Such a bound might allow one
to extend our sketch-and-solve approach to more general SDPs.

##  7 Proof of Theorem 2

For each $S\in\Gamma$, denote $S^{\prime}:=S\cap W$. To prove Theorem 2, we
will first find sufficient conditions for the following approximations to
hold:

$\alpha_{S^{\prime}T^{\prime}}\gtrapprox\alpha_{ST},\qquad\beta_{S^{\prime}T^{\prime}}\approx\beta_{ST},\qquad
c_{S^{\prime}}\approx c_{S}.$

The first two approximations combined with Proposition 1 ensure that the SDP
step of Algorithm 1 produces the clustering
$\Gamma^{\prime}=\\{S^{\prime}:S\in\Gamma\\}$. Meanwhile, the approximation
$c_{S^{\prime}}\approx c_{S}$ ensures that the final step of Algorithm 1
produces the desired clustering $\Gamma$.

######  Lemma 7.

For $X\sim\mathsf{Binomial}(n,p)$, it holds that

$\mathbb{P}\\{X\leq\tfrac{1}{2}pn\\}\leq\operatorname{exp}(-\tfrac{3}{28}pn),\qquad\mathbb{P}\\{X\geq\tfrac{3}{2}pn\\}\leq\operatorname{exp}(-\tfrac{3}{28}pn).$

###### Proof.

This is an immediate consequence of Bernstein’s inequality:

$\mathbb{P}\\{X\leq\tfrac{1}{2}pn\\}=\mathbb{P}\\{pn-X\geq\tfrac{1}{2}pn\\}\leq\operatorname{exp}\bigg{(}-\frac{\frac{1}{2}(\frac{1}{2}pn)^{2}}{np(1-p)+\frac{1}{3}(\frac{1}{2}pn)}\bigg{)}\leq\operatorname{exp}(-\tfrac{3}{28}pn).$

The other bound follows from a similar proof. ∎

######  Proposition 8 (Matrix Bernstein, see Theorem 6.6.1 in [40]).

Consider independent, random, mean-zero, real symmetric $d\times d$ matrices
$\\{X_{i}\\}_{i\in[n]}$ such that

$\lambda_{\mathrm{max}}(X_{i})\leq L~{}~{}~{}\forall
i\in[n],\qquad\qquad\bigg{\|}\sum_{i\in[n]}\mathbb{E}X_{i}^{2}\bigg{\|}_{2\to
2}\leq v.$

Then for every $t\geq 0$, it holds that

$\mathbb{P}\bigg{\\{}\lambda_{\mathrm{max}}\bigg{(}\sum_{i\in[n]}X_{i}\bigg{)}\geq
t\bigg{\\}}\leq
d\cdot\operatorname{exp}\Big{(}-\tfrac{\frac{1}{2}t^{2}}{v+\frac{1}{3}Lt}\Big{)}\leq
d\cdot\operatorname{exp}(-\tfrac{1}{4}\min\\{\tfrac{t^{2}}{v},\tfrac{3t}{L}\\}).$

######  Lemma 9.

Given independent random vectors $\\{X_{i}\\}_{i\in[n]}$ satisfying
$\mathbb{E}X_{i}=0$ and $\|X_{i}\|\leq r$ almost surely for each $i\in[n]$,
put $v:=\max_{i\in[n]}\mathbb{E}\|X_{i}\|^{2}$. Then for $t\geq 0$, it holds
that

$\mathbb{P}\bigg{\\{}\bigg{\|}\sum_{i\in[n]}X_{i}\bigg{\|}\geq
t\bigg{\\}}\leq(d+1)\operatorname{exp}\Big{(}-\tfrac{\frac{1}{2}t^{2}}{nv+\frac{1}{3}rt}\Big{)}\leq(d+1)\operatorname{exp}(-\tfrac{1}{4}\min\\{\tfrac{t^{2}}{nv},\tfrac{3t}{r}\\}).$

###### Proof.

Apply matrix Bernstein to the random matrices
$\left[\begin{smallmatrix}0&X_{i}^{\top}\\\ X_{i}&0\end{smallmatrix}\right]$.
∎

######  Lemma 10.

Given a tuple $\\{x_{i}\\}_{i\in[n]}$ of points in $\mathbb{R}^{d}$, consider
their centroid and radius:

$c:=\frac{1}{n}\sum_{i\in[n]}x_{i},\qquad r:=\max_{i\in[n]}\|x_{i}-c\|.$

Let $\\{b_{i}\\}_{i\in[n]}$ denote independent Bernoulli random variables with
success probability $p$, and consider the random variable
$n^{\prime}:=\sum_{i\in[n]}b_{i}$ and the random vector

$c^{\prime}:=\left\\{\begin{array}[]{cl}\frac{1}{n^{\prime}}\sum_{i\in[n]}b_{i}x_{i}&\text{if
}n^{\prime}>0\\\ \text{undefined}&\text{if }n^{\prime}=0.\end{array}\right.$

Provided $pn\geq 16\log(\frac{d+2}{\epsilon_{\ref{lem.centroid}}})$, it holds
that
$\|c^{\prime}-c\|^{2}<\frac{16r^{2}}{pn}\log(\frac{d+2}{\epsilon_{\ref{lem.centroid}}})$
with probability $\geq 1-\epsilon_{\ref{lem.centroid}}$.

###### Proof.

In the event $\\{n^{\prime}>0\\}$, it holds that

$\|c^{\prime}-c\|=\bigg{\|}\frac{1}{n^{\prime}}\sum_{i\in[n]}b_{i}x_{i}-c\bigg{\|}=\frac{1}{n^{\prime}}\bigg{\|}\sum_{i\in[n]}b_{i}(x_{i}-c)\bigg{\|}=\frac{1}{n^{\prime}}\bigg{\|}\sum_{i\in[n]}(b_{i}-p)(x_{i}-c)\bigg{\|},$

where the last step applies the fact that $\sum_{i\in[n]}(x_{i}-c)=0$. Then

$\displaystyle\mathbb{P}\\{\|c^{\prime}-c\|\geq t\\}$
$\displaystyle\leq\mathbb{P}(\\{\|c^{\prime}-c\|\geq
t\\}\cap\\{n^{\prime}>\tfrac{1}{2}pn\\})+\mathbb{P}\\{n^{\prime}\leq\tfrac{1}{2}pn\\}$
$\displaystyle\leq\mathbb{P}\bigg{\\{}\frac{2}{pn}\bigg{\|}\sum_{i\in[n]}(b_{i}-p)(x_{i}-c)\bigg{\|}\geq
t\bigg{\\}}+\mathbb{P}\\{n^{\prime}\leq\tfrac{1}{2}pn\\}.$ (2)

We apply Lemma 9 to the first term above and Lemma 7 to the second term.
Specifically, put $X_{i}:=(b_{i}-p)(x_{i}-c)$. Then $\mathbb{E}X_{i}=0$ and
$\|X_{i}\|\leq r$ almost surely. Furthermore,

$\mathbb{E}\|X_{i}\|^{2}=\|x_{i}-c\|^{2}\cdot\mathbb{E}(b_{i}-p)^{2}\leq
pr^{2},$

and so $v\leq pr^{2}$. With this, we continue to bound (2):

$\displaystyle\mathbb{P}\\{\|c^{\prime}-c\|\geq t\\}$
$\displaystyle\leq(d+1)\operatorname{exp}\Big{(}-\tfrac{1}{4}\min\Big{\\{}\tfrac{(\frac{pnt}{2})^{2}}{npr^{2}},\tfrac{3(\frac{pnt}{2})}{r}\Big{\\}}\Big{)}+\operatorname{exp}(-\tfrac{3}{28}pn)$
$\displaystyle\leq(d+2)\operatorname{exp}(-\tfrac{pn}{4}\min\\{\tfrac{t^{2}}{4r^{2}},\tfrac{3t}{2r},\tfrac{3}{7}\\})=(d+2)\operatorname{exp}(-\tfrac{pnt^{2}}{16r^{2}}),$
(3)

where the last step holds provided $t\leq r$. The result follows by taking
$t^{2}:=\frac{16r^{2}}{pn}\log(\frac{d+2}{\epsilon_{\ref{lem.centroid}}})$,
since then our assumption $pn\geq
16\log(\frac{d+2}{\epsilon_{\ref{lem.centroid}}})$ implies $t\leq r$. ∎

Here and throughout, we denote

$n_{\mathrm{min}}:=\min_{S\in\Gamma}|S|,\qquad
n_{\mathrm{max}}:=\max_{S\in\Gamma}|S|.$

######  Lemma 11.

Fix $S,T\in\Gamma$ with $S\neq T$ and suppose $pn_{\mathrm{min}}\geq
16\log(\frac{2(d+2)}{\epsilon_{\ref{lem.lower bound alpha}}})$. Then

$\alpha_{S^{\prime}T^{\prime}}>\alpha_{ST}-(\tfrac{2r}{\Delta}+\tfrac{3}{2})\cdot
8r\cdot\sqrt{\tfrac{\log(2(d+2)/\epsilon_{\ref{lem.lower bound
alpha}})}{pn_{\mathrm{min}}}}$

with probability $\geq 1-\epsilon_{\ref{lem.lower bound alpha}}$.

###### Proof.

Denote $m_{ST}:=\tfrac{c_{S}+c_{T}}{2}$ and
$w_{ST}:=\tfrac{c_{S}-c_{T}}{\|c_{S}-c_{T}\|}$ so that $\alpha_{ST}=\min_{i\in
S}\langle x_{i}-m_{ST},w_{ST}\rangle$, and let $j$ denote any minimizer of
$\langle x_{i}-m_{S^{\prime}T^{\prime}},w_{S^{\prime}T^{\prime}}\rangle$ over
$i\in S^{\prime}$. Then

$\displaystyle\alpha_{ST}-\alpha_{S^{\prime}T^{\prime}}$
$\displaystyle\leq\langle x_{j}-m_{ST},w_{ST}\rangle-\langle
x_{j}-m_{S^{\prime}T^{\prime}},w_{S^{\prime}T^{\prime}}\rangle$
$\displaystyle=\langle
x_{j}-m_{ST},w_{ST}-w_{S^{\prime}T^{\prime}}\rangle+\langle
m_{S^{\prime}T^{\prime}}-m_{ST},w_{S^{\prime}T^{\prime}}\rangle$
$\displaystyle\leq\|x_{j}-m_{ST}\|\|w_{S^{\prime}T^{\prime}}-w_{ST}\|+\|m_{S^{\prime}T^{\prime}}-m_{ST}\|.$
(4)

To continue, we bound each of the terms above. First, the triangle inequality
gives

$\|x_{j}-m_{ST}\|=\|x_{j}-c_{S}+c_{S}-m_{ST}\|\leq\|x_{j}-c_{S}\|+\|c_{S}-m_{ST}\|\leq
r+\tfrac{1}{2}\|c_{S}-c_{T}\|,$
$\|m_{S^{\prime}T^{\prime}}-m_{ST}\|=\tfrac{1}{2}\|(c_{S^{\prime}}+c_{T^{\prime}})-(c_{S}+c_{T})\|\leq\tfrac{1}{2}\big{(}\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|\big{)}.$

Next, we apply the triangle inequality multiple times to get

$\displaystyle\|w_{S^{\prime}T^{\prime}}-w_{ST}\|$
$\displaystyle=\|\tfrac{c_{S^{\prime}}-c_{T^{\prime}}}{\|c_{S^{\prime}}-c_{T^{\prime}}\|}-\tfrac{c_{S^{\prime}}-c_{T^{\prime}}}{\|c_{S}-c_{T}\|}+\tfrac{c_{S^{\prime}}-c_{T^{\prime}}}{\|c_{S}-c_{T}\|}-\tfrac{c_{S}-c_{T}}{\|c_{S}-c_{T}\|}\|$
$\displaystyle\leq\|\tfrac{c_{S^{\prime}}-c_{T^{\prime}}}{\|c_{S^{\prime}}-c_{T^{\prime}}\|}-\tfrac{c_{S^{\prime}}-c_{T^{\prime}}}{\|c_{S}-c_{T}\|}\|+\|\tfrac{c_{S^{\prime}}-c_{T^{\prime}}}{\|c_{S}-c_{T}\|}-\tfrac{c_{S}-c_{T}}{\|c_{S}-c_{T}\|}\|$
$\displaystyle=\tfrac{|\|c_{S}-c_{T}\|-\|c_{S^{\prime}}-c_{T^{\prime}}\||}{\|c_{S}-c_{T}\|}+\tfrac{\|c_{S^{\prime}}-c_{T^{\prime}}-c_{S}+c_{T}\|}{\|c_{S}-c_{T}\|}\leq
2\cdot\tfrac{\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|}{\|c_{S}-c_{T}\|}.$

With this, we continue (4):

$\displaystyle\alpha_{ST}-\alpha_{S^{\prime}T^{\prime}}$
$\displaystyle\leq\big{(}r+\tfrac{1}{2}\|c_{S}-c_{T}\|\big{)}\cdot
2\cdot\tfrac{\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|}{\|c_{S}-c_{T}\|}+\tfrac{1}{2}\big{(}\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|\big{)}$
$\displaystyle=\big{(}\tfrac{2r}{\|c_{S}-c_{T}\|}+\tfrac{3}{2}\big{)}\big{(}\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|\big{)}$
$\displaystyle\leq\big{(}\tfrac{2r}{\Delta}+\tfrac{3}{2}\big{)}\big{(}\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|\big{)}.$

It follows that

$\displaystyle\mathbb{P}\\{\alpha_{ST}-\alpha_{S^{\prime}T^{\prime}}\geq t\\}$
$\displaystyle\leq\mathbb{P}\\{\big{(}\tfrac{2r}{\Delta}+\tfrac{3}{2}\big{)}\big{(}\|c_{S^{\prime}}-c_{S}\|+\|c_{T^{\prime}}-c_{T}\|\big{)}\geq
t\\}$
$\displaystyle\leq\mathbb{P}\\{\big{(}\tfrac{2r}{\Delta}+\tfrac{3}{2}\big{)}\|c_{S^{\prime}}-c_{S}\|\geq\tfrac{t}{2}\\}+\mathbb{P}\\{\big{(}\tfrac{2r}{\Delta}+\tfrac{3}{2}\big{)}\|c_{T^{\prime}}-c_{T}\|\geq\tfrac{t}{2}\\}$
$\displaystyle\leq
2\max_{\begin{subarray}{c}R\in\Gamma\end{subarray}}\mathbb{P}\\{\big{(}\tfrac{2r}{\Delta}+\tfrac{3}{2}\big{)}\|c_{R^{\prime}}-c_{R}\|\geq\tfrac{t}{2}\\}.$

The result then follows from Lemma 10 by taking
$\epsilon_{\ref{lem.centroid}}:=\epsilon_{\ref{lem.lower bound alpha}}/2$. ∎

######  Lemma 12.

Fix $S,T\in\Gamma$ with $S\neq T$ and suppose
$pn_{\mathrm{min}}\geq\frac{104}{3}(\frac{n_{\mathrm{max}}}{n_{\mathrm{min}}})^{2}\log(\frac{6}{\epsilon_{\ref{lem.deviation
in sum of reciprocals}}})$. Then

$|p(\tfrac{1}{|S^{\prime}|}+\tfrac{1}{|T^{\prime}|})-(\tfrac{1}{|S|}+\tfrac{1}{|T|})|<\sqrt{\tfrac{104}{3}\tfrac{\log(6/\epsilon_{\ref{lem.deviation
in sum of reciprocals}})}{pn_{\mathrm{min}}^{3}}}$

with probability $\geq 1-\epsilon_{\ref{lem.deviation in sum of
reciprocals}}$.

###### Proof.

First, the triangle inequality gives

$|p(\tfrac{1}{|S^{\prime}|}+\tfrac{1}{|T^{\prime}|})-(\tfrac{1}{|S|}+\tfrac{1}{|T|})|\leq|\tfrac{p}{|S^{\prime}|}-\tfrac{1}{|S|}|+|\tfrac{p}{|T^{\prime}|}-\tfrac{1}{|T|}|.$

As such, it suffices to bound terms of the form

$|\tfrac{p}{|S^{\prime}|}-\tfrac{1}{|S|}|=\tfrac{|p|S|-|S^{\prime}||}{|S^{\prime}||S|}\leq\tfrac{|p|S|-|S^{\prime}||}{\frac{p}{2}|S|^{2}},$

where the last step holds in the event $\\{|S^{\prime}|\geq\frac{p}{2}|S|\\}$.
For every $t\in[0,\frac{1}{|S|}]$, Bernstein’s inequality and Lemma 7 together
give

$\displaystyle\mathbb{P}\big{\\{}|\tfrac{p}{|S^{\prime}|}-\tfrac{1}{|S|}|\geq\tfrac{t}{2}\big{\\}}$
$\displaystyle\leq\mathbb{P}\big{\\{}\big{|}|S^{\prime}|-p|S|\big{|}\geq\tfrac{p}{4}|S|^{2}t\big{\\}}+\mathbb{P}\big{\\{}|S^{\prime}|\leq\tfrac{p}{2}|S|\big{\\}}$
$\displaystyle\leq
2\operatorname{exp}\Big{(}-\tfrac{\frac{1}{2}(\frac{p}{4}|S|^{2}t)^{2}}{|S|p(1-p)+\frac{1}{3}(\frac{p}{4}|S|^{2}t)}\Big{)}+\operatorname{exp}(-\tfrac{3}{28}p|S|)$
$\displaystyle\leq
2\operatorname{exp}\Big{(}-\tfrac{\frac{1}{2}(\frac{p}{4}|S|^{2}t)^{2}}{\frac{13}{12}p|S|}\Big{)}+\operatorname{exp}(-\tfrac{3}{28}p|S|)$
$\displaystyle=2\operatorname{exp}(-\tfrac{3}{104}p|S|^{3}t^{2})+\operatorname{exp}(-\tfrac{3}{28}p|S|)\leq
3\operatorname{exp}(-\tfrac{3}{104}p|S|^{3}t^{2}).$

Finally, we combine our estimates to obtain

$\displaystyle\mathbb{P}\big{\\{}|p(\tfrac{1}{|S^{\prime}|}+\tfrac{1}{|T^{\prime}|})-(\tfrac{1}{|S|}+\tfrac{1}{|T|})|\geq
t\big{\\}}$
$\displaystyle\leq\mathbb{P}\big{\\{}|\tfrac{p}{|S^{\prime}|}-\tfrac{1}{|S|}|\geq\tfrac{t}{2}\big{\\}}+\mathbb{P}\big{\\{}|\tfrac{p}{|T^{\prime}|}-\tfrac{1}{|T|}|\geq\tfrac{t}{2}\big{\\}}$
$\displaystyle\leq
3\operatorname{exp}(-\tfrac{3}{104}p|S|^{3}t^{2})+3\operatorname{exp}(-\tfrac{3}{104}p|T|^{3}t^{2})$
$\displaystyle\leq
6\operatorname{exp}(-\tfrac{3}{104}pn_{\mathrm{min}}^{3}t^{2}).$

The result follows by taking
$t:=\sqrt{\tfrac{104}{3}\tfrac{\log(6/\epsilon_{\ref{lem.deviation in sum of
reciprocals}})}{pn_{\mathrm{min}}^{3}}}$, which is at most
$\frac{1}{n_{\mathrm{max}}}$ by assumption. ∎

######  Lemma 13.

Fix $R\in\Gamma$. Then for every $t\in[0,p|R|r^{2}]$, it holds that

$\mathbb{P}\Big{\\{}\big{|}\|X_{R^{\prime}}\|_{2\to 2}^{2}-p\|X_{R}\|_{2\to
2}^{2}\big{|}\geq
t\Big{\\}}\leq(3d+3)\cdot\operatorname{exp}(-\tfrac{t^{2}}{48p|R|r^{4}}).$

###### Proof.

Fix $R\in\Gamma$. For each $i\in R$, let $b_{i}\sim\mathsf{Bernoulli}(p)$
indicate whether $i\in R^{\prime}$, and consider the matrices

$X_{R}:=\sum_{i\in R}(x_{i}-c_{R})e_{i}^{\top},\qquad
X_{R^{\prime}}:=\sum_{i\in R}b_{i}(x_{i}-c_{R^{\prime}})e_{i}^{\top},\qquad
W_{R}:=\sum_{i\in R}b_{i}(x_{i}-c_{R})e_{i}^{\top}.$

(Here and throughout, we note that $c_{R^{\prime}}$ and any quantity defined
in terms of $c_{R^{\prime}}$ is undefined in the event that $R^{\prime}$ is
empty.) We will bound the deviation of $\|X_{R^{\prime}}\|_{2\to 2}^{2}$ from
$p\|X_{R}\|_{2\to 2}^{2}$ by applying the triangle inequality through
$\|W_{R}\|_{2\to 2}^{2}$. To facilitate this analysis, put
$e_{R}:=c_{R^{\prime}}-c_{R}$. Then
$X_{R^{\prime}}=W_{R}-e_{R}1_{R^{\prime}}^{\top}$ and
$W_{R}1_{R^{\prime}}=|R^{\prime}|e_{R}$, from which it follows that

$X_{R^{\prime}}X_{R^{\prime}}^{\top}-W_{R}W_{R}^{\top}=-|R^{\prime}|e_{R}e_{R}^{\top}.$

As such, the triangle and reverse triangle inequalities together give

$\displaystyle\big{|}\|X_{R^{\prime}}\|_{2\to 2}^{2}-p\|X_{R}\|_{2\to
2}^{2}\big{|}$ $\displaystyle\leq\big{|}\|X_{R^{\prime}}\|_{2\to
2}^{2}-\|W_{R}\|_{2\to 2}^{2}\big{|}+\big{|}\|W_{R}\|_{2\to
2}^{2}-p\|X_{R}\|_{2\to 2}^{2}\big{|}$
$\displaystyle=\big{|}\|X_{R^{\prime}}X_{R^{\prime}}^{\top}\|_{2\to
2}-\|W_{R}W_{R}^{\top}\|_{2\to 2}\big{|}+\big{|}\|W_{R}W_{R}^{\top}\|_{2\to
2}-\|pX_{R}X_{R}^{\top}\|_{2\to 2}\big{|}$
$\displaystyle\leq\|X_{R^{\prime}}X_{R^{\prime}}^{\top}-W_{R}W_{R}^{\top}\|_{2\to
2}+\|W_{R}W_{R}^{\top}-pX_{R}X_{R}^{\top}\|_{2\to 2}$
$\displaystyle=|R^{\prime}|\|e_{R}\|^{2}+\|W_{R}W_{R}^{\top}-pX_{R}X_{R}^{\top}\|_{2\to
2}.$

We will use Lemma 7 to bound $|R^{\prime}|$, Lemma 10 to bound
$\|e_{R}\|^{2}$, and Proposition 8 to bound
$\|W_{R}W_{R}^{\top}-pX_{R}X_{R}^{\top}\|_{2\to 2}$. For this third bound, we
put $X_{i}:=(b_{i}-p)(x_{i}-c_{R})(x_{i}-c_{R})^{\top}$ and observe that

$\|X_{i}\|_{2\to 2}\leq\|x_{i}-c_{R}\|^{2}\leq r^{2}=:L,$ $\bigg{\|}\sum_{i\in
R}\mathbb{E}X_{i}^{2}\bigg{\|}_{2\to 2}\leq\operatorname{tr}\bigg{(}\sum_{i\in
R}\mathbb{E}X_{i}^{2}\bigg{)}=p(1-p)\sum_{i\in R}\|x_{i}-c_{R}\|^{4}\leq
p|R|r^{4}=:v.$

As such, Lemma 7, the bound (3) (which similarly holds with our current choice
of $r$), and Proposition 8 together give

$\displaystyle\mathbb{P}\Big{\\{}\big{|}\|X_{R^{\prime}}\|_{2\to
2}^{2}-p\|X_{R}\|_{2\to 2}^{2}\big{|}\geq t\Big{\\}}$
$\displaystyle\qquad\leq\mathbb{P}\Big{\\{}|R^{\prime}|\geq\tfrac{3}{2}p|R|\Big{\\}}+\mathbb{P}\Big{\\{}\tfrac{3}{2}p|R|\cdot\|e_{R}\|^{2}\geq\tfrac{t}{2}\Big{\\}}+\mathbb{P}\bigg{\\{}\Big{\|}\sum_{i\in
R}X_{i}\Big{\|}_{2\to 2}\geq\tfrac{t}{2}\bigg{\\}}$
$\displaystyle\qquad\leq\operatorname{exp}(-\tfrac{3}{28}p|R|)+(d+2)\operatorname{exp}\Big{(}-\tfrac{p|R|(\frac{t}{3p|R|})}{16r^{2}}\Big{)}+2d\cdot\operatorname{exp}\Big{(}-\tfrac{1}{4}\min\Big{\\{}\tfrac{(\frac{t}{2})^{2}}{|R|pr^{4}},\tfrac{3(\frac{t}{2})}{r^{2}}\Big{\\}}\Big{)}$
$\displaystyle\qquad\leq(3d+3)\cdot\operatorname{exp}\Big{(}-\min\Big{\\{}\tfrac{3p|R|}{28},\tfrac{t}{48r^{2}},\tfrac{t^{2}}{16p|R|r^{4}},\tfrac{3t}{8r^{2}}\Big{\\}}\Big{)}\leq(3d+3)\cdot\operatorname{exp}(-\tfrac{t^{2}}{48p|R|r^{4}}),$

where the last step holds provided $t\in[0,p|R|r^{2}]$. ∎

######  Lemma 14.

Fix $S,T\in\Gamma$ with $S\neq T$ and suppose $pn_{\mathrm{min}}\geq
48(\frac{n_{\mathrm{max}}}{n_{\mathrm{min}}})^{2}\log(\frac{9k(d+1)}{\epsilon_{\ref{lem.deviation
in beta}}})$. Then

$|\beta_{S^{\prime}T^{\prime}}-\beta_{ST}|\leq\tfrac{1}{2}\bigg{(}\Big{(}16\sqrt{3}+\sqrt{\tfrac{104}{3}}\Big{)}\cdot\tfrac{k}{\sqrt{pn_{\mathrm{min}}}}\cdot\tfrac{n_{\mathrm{max}}}{n_{\mathrm{min}}}\cdot
r^{2}\cdot\log^{1/2}(\tfrac{9k(d+1)}{\epsilon_{\ref{lem.deviation in
beta}}})\bigg{)}^{1/2}$

with probability $\geq 1-\epsilon_{\ref{lem.deviation in beta}}$.

###### Proof.

Put $U:=p(\tfrac{1}{|S^{\prime}|}+\tfrac{1}{|T^{\prime}|})$,
$V:=\frac{1}{p}\sum_{R\in\Gamma}\|X_{R^{\prime}}\|_{2\to 2}^{2}$,
$u:=\frac{1}{|S|}+\frac{1}{|T|}$, and $v:=\sum_{R\in\Gamma}\|X_{R}\|_{2\to
2}^{2}$. Then

$\displaystyle|\beta_{S^{\prime}T^{\prime}}-\beta_{ST}|$
$\displaystyle=|\tfrac{1}{2}(UV)^{1/2}-\tfrac{1}{2}(uv)^{1/2}|$
$\displaystyle\leq\tfrac{1}{2}|UV-
uv|^{1/2}=\tfrac{1}{2}|U(V-v)+v(U-u)|^{1/2}\leq\tfrac{1}{2}(U|V-v|+v|U-u|)^{1/2},$

where the first inequality follows from the fact that $x,y\geq 0$ implies

$|x-y|^{2}=x^{2}-2xy+y^{2}\leq
x^{2}-2\min\\{x^{2},y^{2}\\}+y^{2}=|x^{2}-y^{2}|,$

while the second inequality follows from the triangle inequality. Thus, it
suffices to bound $U$, $|V-v|$, $v$, and $|U-u|$ in a high-probability event.
First, Lemma 7 gives that $|S^{\prime}|>\frac{1}{2}p|S|$ with probability
$\geq 1-\exp(-\tfrac{3}{28}p|S|)$, and similarly,
$|T^{\prime}|>\frac{1}{2}p|T|$ with probability $\geq
1-\exp(-\tfrac{3}{28}p|T|)$. A union bound therefore gives

$U=p(\tfrac{1}{|S^{\prime}|}+\tfrac{1}{|T^{\prime}|})<p(\tfrac{2}{p|S|}+\tfrac{2}{p|T|})\leq\tfrac{4}{n_{\mathrm{min}}}$

with probability $\geq 1-2\exp(-\tfrac{3}{28}pn_{\mathrm{min}})$. Next, we
apply the triangle inequality, union bound, and Lemma 13 to get

$|V-v|=\bigg{|}\frac{1}{p}\sum_{R\in\Gamma}\|X_{R^{\prime}}\|_{2\to
2}^{2}-\sum_{R\in\Gamma}\|X_{R}\|_{2\to
2}^{2}\bigg{|}\leq\frac{1}{p}\sum_{R\in\Gamma}\Big{|}\|X_{R^{\prime}}\|_{2\to
2}^{2}-p\|X_{R}\|_{2\to 2}^{2}\Big{|}<\frac{kt}{p}$

with probability $\geq
1-3k(d+1)\exp(-\frac{t^{2}}{48pn_{\mathrm{max}}r^{4}})$, provided
$t\in[0,pn_{\mathrm{min}}r^{2}]$. For $v$, we pass to the Frobenius norm:

$v=\sum_{R\in\Gamma}\|X_{R}\|_{2\to
2}^{2}\leq\sum_{R\in\Gamma}\|X_{R}\|_{F}^{2}\leq nr^{2}\leq
kn_{\mathrm{max}}r^{2}.$

Finally, we apply Lemma 12 to obtain

$|U-u|=|p(\tfrac{1}{|S^{\prime}|}+\tfrac{1}{|T^{\prime}|})-(\tfrac{1}{|S|}+\tfrac{1}{|T|})|<\sqrt{\tfrac{104}{3}\tfrac{\log(6/\epsilon_{\ref{lem.deviation
in sum of reciprocals}})}{pn_{\mathrm{min}}^{3}}}$

with probability $\geq 1-\epsilon_{\ref{lem.deviation in sum of
reciprocals}}$, provided
$pn_{\mathrm{min}}\geq\frac{104}{3}(\frac{n_{\mathrm{max}}}{n_{\mathrm{min}}})^{2}\log(\frac{6}{\epsilon_{\ref{lem.deviation
in sum of reciprocals}}})$. We will combine these bounds using a union bound.
To do so, we will bound the failure probabilities corresponding to $U$,
$|V-v|$, and $|U-u|$ by $\tfrac{\epsilon_{\ref{lem.deviation in beta}}}{3}$:

$\tfrac{\epsilon_{\ref{lem.deviation in beta}}}{3}\geq
2\exp(-\tfrac{3}{28}pn_{\mathrm{min}}),\qquad\tfrac{\epsilon_{\ref{lem.deviation
in beta}}}{3}\geq
3k(d+1)\exp(-\tfrac{t^{2}}{48pn_{\mathrm{max}}r^{4}}),\qquad\tfrac{\epsilon_{\ref{lem.deviation
in beta}}}{3}\geq\epsilon_{\ref{lem.deviation in sum of reciprocals}}.$ (5)

Rearranging the first bound in (5) gives
$pn_{\mathrm{min}}\geq\tfrac{28}{3}\log(\frac{6}{\epsilon_{\ref{lem.deviation
in beta}}})$, which is implied by our hypothesis. We change the second bound
in (5) to an equality that defines $t$ as

$t:=\sqrt{48pn_{\mathrm{max}}r^{4}\log(\tfrac{9k(d+1)}{\epsilon_{\ref{lem.deviation
in beta}}})}.$

This choice satisfies the requirement that $t\in[0,pn_{\mathrm{min}}r^{2}]$
precisely when

$pn_{\mathrm{min}}\geq
48\cdot\tfrac{n_{\mathrm{max}}}{n_{\mathrm{min}}}\cdot\log(\tfrac{9k(d+1)}{\epsilon_{\ref{lem.deviation
in beta}}}),$

which is implied by our hypothesis. Finally, we change the third bound in (5)
to an equality that defines $\epsilon_{\ref{lem.deviation in sum of
reciprocals}}:=\tfrac{\epsilon_{\ref{lem.deviation in beta}}}{3}$, which
satisfies the requirement
$pn_{\mathrm{min}}\geq\frac{104}{3}(\frac{n_{\mathrm{max}}}{n_{\mathrm{min}}})^{2}\log(\frac{6}{\epsilon_{\ref{lem.deviation
in sum of reciprocals}}})$ by our hypothesis. Putting everything together, we
have

$\displaystyle|\beta_{S^{\prime}T^{\prime}}-\beta_{ST}|$
$\displaystyle\leq\tfrac{1}{2}(U|V-v|+v|U-u|)^{1/2}$
$\displaystyle<\tfrac{1}{2}\bigg{(}\tfrac{4}{n_{\mathrm{min}}}\cdot\tfrac{k}{p}\sqrt{48pn_{\mathrm{max}}r^{4}\log(\tfrac{9k(d+1)}{\epsilon_{\ref{lem.deviation
in
beta}}})}+kn_{\mathrm{max}}r^{2}\cdot\sqrt{\tfrac{104}{3}\tfrac{\log(18/\epsilon_{\ref{lem.deviation
in beta}})}{pn_{\mathrm{min}}^{3}}}\bigg{)}^{1/2}$
$\displaystyle\leq\tfrac{1}{2}\bigg{(}\Big{(}16\sqrt{3}+\sqrt{\tfrac{104}{3}}\Big{)}\cdot\tfrac{k}{\sqrt{pn_{\mathrm{min}}}}\cdot\tfrac{n_{\mathrm{max}}}{n_{\mathrm{min}}}\cdot
r^{2}\cdot\log^{1/2}(\tfrac{9k(d+1)}{\epsilon_{\ref{lem.deviation in
beta}}})\bigg{)}^{1/2},$

where $16\sqrt{3}$ comes from the first term and $\sqrt{\tfrac{104}{3}}$ comes
from the second term. (To be clear, we used the fact that
$n_{\mathrm{max}}\geq n_{\mathrm{min}}$ to bound the first term, and the fact
that $9k(d+1)\geq 18$ to bound the second term.) ∎

######  Lemma 15.

Suppose $pn_{\mathrm{min}}\geq
48(\frac{n_{\mathrm{max}}}{n_{\mathrm{min}}})^{2}\log(\frac{18k^{3}(d+1)}{\epsilon_{\ref{lem.approx
prox condition}}})$ and

$\displaystyle\operatorname{prox}(X,\Gamma)$
$\displaystyle>(\tfrac{2r}{\Delta}+\tfrac{3}{2})\cdot
8r\cdot\sqrt{\tfrac{\log(4k^{2}(d+2)/\epsilon_{\ref{lem.approx prox
condition}})}{pn_{\mathrm{min}}}}$
$\displaystyle\qquad+\tfrac{1}{2}\bigg{(}\Big{(}16\sqrt{3}+\sqrt{\tfrac{104}{3}}\Big{)}\cdot\tfrac{k}{\sqrt{pn_{\mathrm{min}}}}\cdot\tfrac{n_{\mathrm{max}}}{n_{\mathrm{min}}}\cdot
r^{2}\cdot\log^{1/2}(\tfrac{18k^{3}(d+1)}{\epsilon_{\ref{lem.approx prox
condition}}})\bigg{)}^{1/2}.$

Then $\operatorname{prox}(X^{\prime},\Gamma^{\prime})>0$ with probability
$\geq 1-\epsilon_{\ref{lem.approx prox condition}}$.

###### Proof.

We apply Lemmas 11 and 14 with $\epsilon_{\ref{lem.lower bound
alpha}}=\epsilon_{\ref{lem.deviation in
beta}}:=\frac{\epsilon_{\ref{lem.approx prox condition}}}{2k^{2}}$. By taking
union bound over $S,T\in\Gamma$ with $S\neq T$, the random variable
$\operatorname{prox}(X,\Gamma)-\operatorname{prox}(X^{\prime},\Gamma^{\prime})$
is at most the right-hand side of the displayed inequality with probability
$\geq 1-\epsilon_{\ref{lem.approx prox condition}}$, and the result follows. ∎

######  Lemma 16.

Suppose $r\leq\frac{\Delta}{2}$ and $pn_{\mathrm{min}}\geq
16\max\\{1,(\frac{r}{\Delta/2-r})^{2}\\}\log(\tfrac{k(d+2)}{\epsilon_{\ref{lem.rounding
after sketch and solve}}})$. Then with probability $\geq
1-\epsilon_{\ref{lem.rounding after sketch and solve}}$, it simultaneously
holds that

$\|x_{i}-c_{S^{\prime}}\|<\|x_{i}-c_{T^{\prime}}\|$

for every $S,T\in\Gamma$ with $S\neq T$ and every $i\in S$.

###### Proof.

We will show that a stronger condition holds, namely that

$\max_{\begin{subarray}{c}S,T\in\Gamma\\\ S\neq T\end{subarray}}\max_{i\in
S}\|x_{i}-c_{S^{\prime}}\|<\min_{\begin{subarray}{c}S,T\in\Gamma\\\ S\neq
T\end{subarray}}\min_{i\in S}\|x_{i}-c_{T^{\prime}}\|$ (6)

with probability $\geq 1-\epsilon_{\ref{lem.rounding after sketch and
solve}}$. Denote the random variable
$E:=\max_{S\in\Gamma}\|c_{S^{\prime}}-c_{S}\|$. Then

$\displaystyle\|x_{i}-c_{S^{\prime}}\|$
$\displaystyle\leq\|x_{i}-c_{S}\|+\|c_{S^{\prime}}-c_{S}\|\leq r+E,$
$\displaystyle\|x_{i}-c_{T^{\prime}}\|$
$\displaystyle\geq\|c_{S}-c_{T}\|-\|c_{T^{\prime}}-c_{T}\|-\|x_{i}-c_{S}\|\geq\Delta-
E-r.$

Thus, the desired inequality (6) holds whenever $E<\frac{\Delta}{2}-r$. The
result then follows from Lemma 10 by taking
$\epsilon_{\ref{lem.centroid}}:=\frac{\epsilon_{\ref{lem.rounding after sketch
and solve}}}{k}$ and applying a union bound over $S\in\Gamma$. ∎

###### Proof of Theorem 2.

Lemmas 15 and 16 with $\epsilon_{\ref{lem.approx prox
condition}}=\epsilon_{\ref{lem.rounding after sketch and
solve}}:=\frac{\epsilon}{2}$ together imply that Algorithm 1 exactly recovers
$\Gamma$ from $X$ with probability $1-\epsilon$ provided both of the following
hold:

$\displaystyle pn_{\mathrm{min}}$
$\displaystyle\geq\max\Big{\\{}48(\tfrac{n_{\mathrm{max}}}{n_{\mathrm{min}}})^{2}\log(\tfrac{36k^{3}(d+1)}{\epsilon}),16\max\\{1,(\tfrac{r}{\Delta/2-r})^{2}\\}\log(\tfrac{2k(d+2)}{\epsilon})\Big{\\}},$
$\displaystyle\operatorname{prox}(X,\Gamma)$
$\displaystyle>(\tfrac{2r}{\Delta}+\tfrac{3}{2})\cdot
8r\cdot\sqrt{\tfrac{\log(8k^{2}(d+2)/\epsilon)}{pn_{\mathrm{min}}}}$
$\displaystyle\qquad+\tfrac{1}{2}\bigg{(}\Big{(}16\sqrt{3}+\sqrt{\tfrac{104}{3}}\Big{)}\cdot\tfrac{k}{\sqrt{pn_{\mathrm{min}}}}\cdot\tfrac{n_{\mathrm{max}}}{n_{\mathrm{min}}}\cdot
r^{2}\cdot\log^{1/2}(\tfrac{36k^{3}(d+1)}{\epsilon})\bigg{)}^{1/2}.$

As we now discuss, there exists an explicit shape parameter $C(X,\Gamma)>0$
such that the inequality $\mathbb{E}|W|\geq C(X,\Gamma)\cdot\log(1/\epsilon)$
implies the above conditions. Denote

$\pi_{\mathrm{min}}:=\frac{1}{n}\min_{S\in\Gamma}|S|,\qquad\pi_{\mathrm{max}}:=\frac{1}{n}\max_{S\in\Gamma}|S|,$

and observe that $pn_{\mathrm{min}}=\mathbb{E}|W|\cdot\pi_{\mathrm{min}}$.
This explains the appearance of $\mathbb{E}|W|$ in our desired inequality. To
isolate $\log(1/\epsilon)$, we will use the general observation that
$\alpha,\beta\geq\gamma>1$ implies

$\log(\alpha\beta)\leq\frac{2\log(\alpha)\log(\beta)}{\log(\gamma)}.$

Indeed,
$\frac{\log(\alpha\beta)}{\log(\alpha)\log(\beta)}=\frac{1}{\log(\alpha)}+\frac{1}{\log(\beta)}\leq\frac{2}{\log(\gamma)}$.
We apply this bound several times with $\beta=1/\epsilon$ and $\gamma=2$ so
that the following conditions imply the above conditions:

$\displaystyle\mathbb{E}|W|\cdot\pi_{\mathrm{min}}$
$\displaystyle\geq\tfrac{2}{\log(2)}\cdot\log(\tfrac{1}{\epsilon})\cdot\max\Big{\\{}48(\tfrac{\pi_{\mathrm{max}}}{\pi_{\mathrm{min}}})^{2}\log(36k^{3}(d+1)),$
$\displaystyle\hskip
195.12877pt16\max\\{1,(\tfrac{\Delta}{2r}-1)^{-2}\\}\log(2k(d+2))\Big{\\}},$
$\displaystyle\operatorname{prox}(X,\Gamma)>(\tfrac{2r}{\Delta}+\tfrac{3}{2})\cdot
8r\cdot\sqrt{\tfrac{\log(8k^{2}(d+2))}{\mathbb{E}|W|\cdot\pi_{\mathrm{min}}}\cdot\tfrac{2}{\log(2)}\cdot\log(\tfrac{1}{\epsilon})}$
$\displaystyle\qquad+\tfrac{1}{2}\bigg{(}\Big{(}16\sqrt{3}+\sqrt{\tfrac{104}{3}}\Big{)}\cdot\tfrac{k}{\sqrt{\mathbb{E}|W|\cdot\pi_{\mathrm{min}}}}\cdot\tfrac{\pi_{\mathrm{max}}}{\pi_{\mathrm{min}}}\cdot
r^{2}\cdot\log^{1/2}(36k^{3}(d+1))\cdot(\tfrac{2}{\log(2)}\cdot\log(\tfrac{1}{\epsilon}))^{1/2}\bigg{)}^{1/2}.$

Notice that we may express these conditions in terms of
$c:=\mathbb{E}|W|/\log(\tfrac{1}{\epsilon})$:

$\displaystyle c\cdot\pi_{\mathrm{min}}$
$\displaystyle\geq\tfrac{2}{\log(2)}\cdot\max\Big{\\{}48(\tfrac{\pi_{\mathrm{max}}}{\pi_{\mathrm{min}}})^{2}\log(36k^{3}(d+1)),16\max\\{1,(\tfrac{\Delta}{2r}-1)^{-2}\\}\log(2k(d+2))\Big{\\}},$
$\displaystyle\tfrac{\operatorname{prox}(X,\Gamma)}{r}$
$\displaystyle>(\tfrac{2r}{\Delta}+\tfrac{3}{2})\cdot
8\cdot\sqrt{\tfrac{\log(8k^{2}(d+2))}{c\cdot\pi_{\mathrm{min}}}\cdot\tfrac{2}{\log(2)}}$
$\displaystyle\quad+\tfrac{1}{2}\bigg{(}\Big{(}16\sqrt{3}+\sqrt{\tfrac{104}{3}}\Big{)}\cdot\tfrac{k}{\sqrt{c\cdot\pi_{\mathrm{min}}}}\cdot\tfrac{\pi_{\mathrm{max}}}{\pi_{\mathrm{min}}}\cdot\log^{1/2}(36k^{3}(d+1))\cdot(\tfrac{2}{\log(2)})^{1/2}\bigg{)}^{1/2}.$

The set of $c$ for which the first inequality holds is an interval the form
$[c_{1},\infty)$, while the set of $c$ for which the second inequality holds
takes the form $(c_{2},\infty)$. Then $C(X,\Gamma):=\max\\{c_{1},c_{2}\\}$ is
an explicit function of $d$, $\frac{\Delta}{r}$,
$\tfrac{\operatorname{prox}(X,\Gamma)}{r}$, $k$, $\pi_{\mathrm{min}}$, and
$\pi_{\mathrm{max}}$, as desired. ∎

##  8 Proof of Theorems 4(b) and 5(b)

######  Lemma 17 (cf. Lemma 10 in [31]).

Given a symmetric matrix $M\in\mathbb{R}^{n\times n}$, it holds that

$\max_{Z\in\mathcal{Z}(n,k)}|\operatorname{tr}(MZ)|\leq\min\Big{\\{}\|M\|_{*},k\|M\|_{2\to
2}\Big{\\}}.$

###### Proof.

Select any $Z\in\mathcal{Z}(n,k)$, and let
$\alpha_{1}\geq\cdots\geq\alpha_{n}$ and $\beta_{1}\geq\cdots\geq\beta_{n}$
denote the singular values of $M$ and $Z$, respectively. Then Von Neumann’s
trace inequality gives that

$|\operatorname{tr}(MZ)|\leq\sum_{i\in[n]}\alpha_{i}\beta_{i}.$

Since $Z$ is stochastic, then by the Gershgorin circle theorem, every
eigenvalue of $Z$ has modulus at most $1$. Since $Z$ is symmetric, it follows
that $\beta_{1}\leq 1$, and so

$|\operatorname{tr}(MZ)|\leq\sum_{i\in[n]}\alpha_{i}\beta_{i}\leq\beta_{1}\sum_{i\in[n]}\alpha_{i}\leq\|M\|_{*}.$

We also have $\sum_{i\in[n]}\beta_{i}=\operatorname{tr}Z=k$, and so

$|\operatorname{tr}(MZ)|\leq\sum_{i\in[n]}\alpha_{i}\beta_{i}\leq\alpha_{1}\sum_{i\in[n]}\beta_{i}=k\|M\|_{2\to
2}.\qed$

######  Lemma 18.

Given any tuple $X:=\\{x_{i}\\}_{i\in[n]}$ of points in $\mathbb{R}^{d}$ and
any orthogonal projection matrix $P\in\mathbb{R}^{d\times d}$, it holds that

$\min_{Z\in\mathcal{Z}(n,k)}\frac{1}{2}\operatorname{tr}(D_{X}Z)\geq\|PX\|_{F}^{2}-k\|PX\|_{2\to
2}^{2}.$

(Here, we abuse notation by identifying $X$ with a member of
$\mathbb{R}^{d\times n}$.)

###### Proof.

Define $\nu\in\mathbb{R}^{n}$ to have $i$th coordinate $\|x_{i}\|^{2}$. Then
we have

$D_{X}=\nu 1^{\top}-2X^{\top}X+1\nu^{\top},\qquad\|X\|_{F}^{2}=1^{\top}\nu.$

Fix $Z\in\mathcal{Z}(n,k)$. Then $Z^{\top}=Z$ and $Z1=1$, and so

$\frac{1}{2}\operatorname{tr}(D_{X}Z)=\frac{1}{2}\operatorname{tr}((\nu
1^{\top}-2X^{\top}X+1\nu^{\top})Z)=\|X\|_{F}^{2}-\operatorname{tr}(X^{\top}XZ).$
(7)

We apply Lemma 17 to get

$\displaystyle\operatorname{tr}(X^{\top}XZ)$
$\displaystyle=\operatorname{tr}(X^{\top}PXZ)+\operatorname{tr}(X^{\top}(I-P)XZ)$
$\displaystyle\leq k\|X^{\top}PX\|_{2\to
2}+\|X^{\top}(I-P)X\|_{*}=k\|PX\|_{2\to 2}^{2}+\|(I-P)X\|_{F}^{2},$

and combining with (7) gives

$\frac{1}{2}\operatorname{tr}(D_{X}Z)\geq\|X\|_{F}^{2}-\Big{(}k\|PX\|_{2\to
2}^{2}+\|(I-P)X\|_{F}^{2}\Big{)}=\|PX\|_{F}^{2}-k\|PX\|_{2\to 2}^{2}.\qed$

######  Proposition 19 (Matrix Chernoff inequality, see Theorem 5.1.1 and
(5.1.8) in [40]).

Consider a finite sequence $\\{X_{k}\\}$ of independent, random, Hermitian
matrices with common dimension d, and assume that
$0\leq\lambda_{\mathrm{min}}(X_{k})\leq\lambda_{\mathrm{max}}(X_{k})\leq L$
almost surely for each $k$. Then

$\mathbb{E}\lambda_{\mathrm{max}}\bigg{(}\sum_{k}X_{k}\bigg{)}\leq
1.72\cdot\lambda_{\mathrm{max}}\bigg{(}\sum_{k}\mathbb{E}X_{k}\bigg{)}+L\log
d.$

######  Proposition 20 (Dvoretzky–Kiefer–Wolfowitz inequality, Theorem 11.5 in
[23]).

Consider a sequence $\\{X_{k}\\}_{k\in[n]}$ of real-valued independent random
variables with common cumulative distribution function
$F\colon\mathbb{R}\to[0,1]$, and let $F_{n}\colon\mathbb{R}\to[0,1]$ denote
the random empirical distribution function defined by
$F_{n}(x):=\frac{1}{n}\sum_{k\in[n]}1_{\\{X_{k}\leq x\\}}$. Then for every
$\epsilon>0$, it holds that

$\mathbb{P}\Big{\\{}\sup_{x\in\mathbb{R}}|F_{n}(x)-F(x)|>\epsilon\Big{\\}}\leq
2e^{-2n\epsilon^{2}}.$

######  Proposition 21 (Lemma 1 in [25]).

Suppose $X$ has chi-squared distribution with $d$ degrees of freedom. Then for
each $t>0$, it holds that

$\mathbb{P}\\{X\geq d+2\sqrt{dt}+2t\\}\leq e^{-t},\qquad\mathbb{P}\\{X\leq
d-2\sqrt{dt}\\}\leq e^{-t}.$

It is sometimes easier to interact with a simpler (though weaker) version of
the upper chi-squared tail estimate:

######  Corollary 22.

Suppose $X$ has chi-squared distribution with $d$ degrees of freedom. Then
$\mathbb{P}\\{X\geq x\\}\leq e^{-x/5}$ for every $x\geq 5d$.

###### Proof.

Put $t:=x/5\geq d$, and so $d+2\sqrt{dt}+2t\leq t+2t+2t=5t$, and Proposition
21 gives

$\mathbb{P}\\{X\geq x\\}=\mathbb{P}\\{X\geq 5t\\}\leq\mathbb{P}\\{X\geq
d+2\sqrt{dt}+2t\\}\leq e^{-t}=e^{-x/5}.\qed$

######  Corollary 23.

Suppose $\\{X_{k}\\}_{k\in[n]}$ are independent random variables with chi-
squared distribution with $d$ degrees of freedom. Then
$\max_{k\in[n]}X_{k}\leq 2d+6\log n$ with probability $\geq 1-\frac{1}{n}$.

###### Proof.

By the arithmetic mean–geometric mean inequality and Proposition 21, we have

$\mathbb{P}\\{X_{k}\geq 2d+3t\\}\leq\mathbb{P}\\{X_{k}\geq
d+2\sqrt{dt}+2t\\}\leq e^{-t}.$

Take $t:=2\log n$ and apply a union bound to get the result. ∎

######  Proposition 24 (Corollary 5.35 in [42]).

Let $G$ be an $m\times n$ matrix with independent standard Gaussian entries.
Then for every $t\geq 0$, it holds that $\|G\|_{2\to
2}\leq\sqrt{m}+\sqrt{n}+t$ with probability $\geq 1-2e^{-t^{2}/2}$.

Throughout, we let $\mathbb{E}_{G}$ denote expectation conditioned on $G$, and
we assume $n\geq s\geq d\geq k$ without mention. (These inequalities are
typically implied by the hypotheses of our results, but even so, it is helpful
to keep these bounds in mind when interpreting the analysis.)

######  Lemma 25.

Let $G:=\\{g_{i}\\}_{i\in[n]}$ denote a tuple of independent random vectors
$g_{i}\sim\mathsf{N}(0,I_{d})$, let $S:=\\{i_{j}\\}_{j\in[s]}$ denote a tuple
of independent random indices $i_{j}\sim\mathsf{Unif}([n])$, and let
$H=H(G,S)\in\mathbb{R}^{d\times s}$ denote the random matrix whose $j$th
column is $g_{i_{j}}$. Assuming $s\geq 15d\log d$, it holds that
$\mathbb{E}_{G}\|H\|_{2\to 2}^{2}\leq 5s$ with probability $\geq
1-\frac{1}{n}-e^{-\Omega(n/(d+\log n)^{2})}$.

###### Proof.

Fix a threshold $\tau>0$ to be selected later, and for any vector
$x\in\mathbb{R}^{d}$, denote

$x^{-}:=\left\\{\begin{array}[]{cl}x&\text{if }\|x\|^{2}\leq\tau\\\
0&\text{otherwise}\end{array}\right\\},\qquad
x^{+}:=\left\\{\begin{array}[]{cl}x&\text{if }\|x\|^{2}>\tau\\\
0&\text{otherwise}\end{array}\right\\}.$

Letting $h_{j}$ denote the $j$th column of $H$, then we are interested in the
quantity

$\displaystyle\|H\|_{2\to 2}^{2}$
$\displaystyle=\bigg{\|}\sum_{j\in[s]}h_{j}h_{j}^{\top}\bigg{\|}_{2\to
2}\leq\bigg{\|}\sum_{j\in[s]}h_{j}^{-}{h_{j}^{-}}^{\top}\bigg{\|}_{2\to
2}+\bigg{\|}\sum_{j\in[s]}h_{j}^{+}{h_{j}^{+}}^{\top}\bigg{\|}_{2\to 2}.$

We will bound the expectation of the first term using the Matrix Chernoff
inequality (Proposition 19) and the expectation of the second term using the
Dvoretzky–Kiefer–Wolfowitz inequality (Proposition 20). First,
$\mathbb{E}_{G}[h_{j}^{-}{h_{j}^{-}}^{\top}]=\frac{1}{n}\sum_{i\in[n]}g_{i}^{-}{g_{i}^{-}}^{\top}$
for each $j\in[s]$, and so

$\lambda_{\mathrm{max}}\bigg{(}\sum_{j\in[s]}\mathbb{E}_{G}\Big{[}h_{j}^{-}{h_{j}^{-}}^{\top}\Big{]}\bigg{)}=\lambda_{\mathrm{max}}\bigg{(}\frac{s}{n}\sum_{i\in[n]}g_{i}^{-}{g_{i}^{-}}^{\top}\bigg{)}\leq\lambda_{\mathrm{max}}\bigg{(}\frac{s}{n}\sum_{i\in[n]}g_{i}g_{i}^{\top}\bigg{)}=\frac{s}{n}\|G\|_{2\to
2}^{2},$

where the inequality uses the fact that
$\frac{s}{n}\sum_{i\in[n]}g_{i}^{+}{g_{i}^{+}}^{\top}\succeq 0$. Thus, Matrix
Chernoff gives

$\mathbb{E}_{G}\bigg{\|}\sum_{j\in[s]}h_{j}^{-}{h_{j}^{-}}^{\top}\bigg{\|}_{2\to
2}\leq 1.72\cdot\frac{s}{n}\|G\|_{2\to 2}^{2}+\tau\log d.$

Next, we bound the expectation of

$\bigg{\|}\sum_{j\in[s]}h_{j}^{+}{h_{j}^{+}}^{\top}\bigg{\|}_{2\to
2}\leq\sum_{j\in[s]}\|h_{j}^{+}{h_{j}^{+}}^{\top}\|_{2\to
2}=\sum_{j\in[s]}\|h_{j}^{+}\|^{2}.$

Let $F_{G}^{+}$ denote the empirical distribution function of
$\\{\|g_{i}^{+}\|^{2}\\}_{i\in[n]}$. Then

$\frac{1}{s}\mathbb{E}_{G}\bigg{\|}\sum_{j\in[s]}h_{j}^{+}{h_{j}^{+}}^{\top}\bigg{\|}_{2\to
2}\leq\mathbb{E}_{G}\|h_{1}^{+}\|^{2}=\int_{0}^{\infty}(1-F_{G}^{+}(x))dx.$
(8)

Let $F_{G}$ denote the empirical distribution function of
$\\{\|g_{i}\|^{2}\\}_{i\in[n]}$, put $\|G\|_{1\to
2}^{2}:=\max_{i\in[n]}\|g_{i}\|^{2}$, and observe that

$F_{G}^{+}(x)=\left\\{\begin{array}[]{cl}F_{G}(\tau)&\text{if }x\leq\tau\\\
F_{G}(x)&\text{if }x>\tau\end{array}\right\\}\quad\text{and}\quad
F_{G}(x)=1\quad\forall x>\|G\|_{1\to 2}^{2}.$

In particular, we may assume $\|G\|_{1\to 2}^{2}\geq\tau$ without loss of
generality, since otherwise the integral in (8) equals zero. We will estimate
this integral in pieces:

$\int_{0}^{\infty}=\int_{0}^{\tau}+\int_{\tau}^{\|G\|_{1\to
2}^{2}}+\int_{\|G\|_{1\to 2}^{2}}^{\infty}.$

We have $\int_{\|G\|_{1\to 2}^{2}}^{\infty}=0$, and we will estimate the other
two integrals in terms of the quantity

$E(G):=\sup_{x\in\mathbb{R}}|F_{G}(x)-F(x)|,$

where $F$ denotes the cumulative distribution function of the chi-squared
distribution with $d$ degrees of freedom. (Later, we will apply the
Dvoretzky–Kiefer–Wolfowitz inequality to bound $E(G)$ with high probability on
$G$.) First, we have

$\int_{0}^{\tau}(1-F_{G}^{+}(x))dx=\tau(1-F_{G}(\tau))\leq\tau(1-F(\tau)+E(G))\leq\tau(e^{-\tau/5}+E(G)),$

where the last step applies Corollary 22 under the assumption that $\tau\geq
5d$. Next, for each $x\in(\tau,\|G\|_{1\to 2}^{2})$, we have

$1-F_{G}^{+}(x)=1-F_{G}(x)\leq 1-F(x)+E(G),$

and so integrating gives

$\int_{\tau}^{\|G\|_{1\to
2}^{2}}(1-F_{G}^{+}(x))dx\leq\int_{\tau}^{\infty}(1-F(x))dx+\|G\|_{1\to
2}^{2}\cdot E(G).$

We estimate the first term using Corollary 22:

$\int_{\tau}^{\infty}(1-F(x))dx\leq\int_{\tau}^{\infty}e^{-x/5}dx=5e^{-\tau/5}.$

All together, we have

$\mathbb{E}_{G}\|H\|_{2\to 2}^{2}\leq 1.72\cdot\frac{s}{n}\|G\|_{2\to
2}^{2}+\tau\log d+s\Big{(}\tau(e^{-\tau/5}+E(G))+5e^{-\tau/5}+\|G\|_{1\to
2}^{2}\cdot E(G)\Big{)}.$

It remains to bound this random variable in a high-probability event.
Proposition 24 gives $\|G\|_{2\to 2}^{2}=O(n)$ with high probability when
$n\geq d$. This means the first term in our bound will be $O(s)$, and so we
select $\tau:=s/\log d$ so that the second term has the same order. (Note that
our assumption $\tau\geq 5d$ then requires $s\geq 5d\log d$.) The remaining
term is $s$ times

$\tau(e^{-\tau/5}+E(G))+5e^{-\tau/5}+\|G\|_{1\to 2}^{2}\cdot
E(G)\leq(\tau+5)e^{-\tau/5}+2\|G\|_{1\to 2}^{2}\cdot E(G),$ (9)

where the inequality follows from the bound $\tau\leq\|G\|_{1\to 2}^{2}$. We
want (9) to be $O(1)$. The first term in (9) is smaller than $1$ provided
$\tau\geq 15$. For the second term in (9), we recall from Corollary 23 that
$\|G\|_{1\to 2}^{2}=O(d+\log n)$ with high probability. This suggests that we
restrict to an event in which $E(G)=O(1/(d+\log n))$, which we can estimate
using the Dvoretzky–Kiefer–Wolfowitz inequality. For this bound, the failure
probability will be $2e^{-2n\epsilon^{2}}=\exp(-\Omega(n/(d+\log n)^{2}))$,
and this informs how sharply we can bound $\|G\|_{2\to 2}^{2}$. Overall, we
restrict to an event in which three things occur simultaneously:

$\|G\|_{2\to 2}^{2}\leq 1.1n,\qquad\|G\|_{1\to 2}^{2}\leq 2d+6\log n,\qquad
E(G)\leq\frac{1}{4d+12\log n}.$

A union bound gives that this event has probability $\geq
1-\frac{1}{n}-e^{-\Omega(n/(d+\log n)^{2})}$, and over this event, it holds
that $\mathbb{E}_{G}\|H\|_{2\to 2}^{2}\leq 1.72\cdot 1.1s+s+s(1+1)\leq 5s$. ∎

######  Lemma 26.

Draw $X:=\\{x_{i}\\}_{i\in[n]}$ in $\mathbb{R}^{d}$ from a mixture of $k$
gaussians with equal weights and identity covariance. Explicitly, take any
$\mu_{1},\ldots,\mu_{k}\in\mathbb{R}^{d}$, draw $\\{t_{i}\\}_{i\in[n]}$
independently with distribution $\mathsf{Unif}([k])$, draw
$\\{g_{i}\\}_{i\in[n]}$ independently with distribution $\mathsf{N}(0,I_{d})$,
and take $x_{i}:=\mu_{t_{i}}+g_{i}$. Next, draw $\\{i_{j}\\}_{j\in[s]}$
independently with distribution $\mathsf{Unif}([n])$ and define the random
tuple $Y:=\\{y_{j}\\}_{j\in[s]}$ by $y_{j}:=x_{i_{j}}$. Then provided $s\geq
15d\log d$, it holds that

$d-6k-1\leq\mathbb{E}_{X}\operatorname{SDP}(Y,k)\leq\operatorname{IP}(X,k)\leq
d+1$

with probability $\geq 1-\frac{1}{n}-e^{-\Omega(n/(d+\log n)^{2})}$.

###### Proof.

Lemma 3 gives the middle inequality in our claim. Next, we obtain an upper
bound on $\operatorname{IP}(X,k)$ by passing to the partition
$\Gamma\in\Pi(n,k)$ defined by the level sets of the planted assignment
$i\mapsto t_{i}$:

$\operatorname{IP}(X,k)\leq\frac{1}{n}\sum_{S\in\Gamma}\sum_{i\in
S}\bigg{\|}x_{i}-\frac{1}{|S|}\sum_{j\in
S}x_{j}\bigg{\|}^{2}\leq\frac{1}{n}\sum_{S\in\Gamma}\sum_{i\in
S}\|x_{i}-\mu_{t_{i}}\|^{2}=\frac{1}{n}\sum_{i\in[n]}\|g_{i}\|^{2},$

where the second inequality follows from the fact that the centroid of a tuple
of points minimizes the sum of squared distances from those points. Next,
$\sum_{i\in[n]}\|g_{i}\|^{2}$ has chi-squared distribution with $dn$ degrees
of freedom, and so an application of Proposition 21 gives that
$\operatorname{IP}(X,k)\leq d+1$ with probability $\geq 1-e^{-n/(16d)}$. It
remains find a lower bound on $\mathbb{E}_{X}\operatorname{SDP}(Y,k)$. For
this, we will apply Lemma 18 with $P$ representing the orthogonal projection
map onto $(\operatorname{span}\\{\mu_{t}\\}_{t\in[k]})^{\perp}$. Letting
$H:=\\{h_{j}\\}_{j\in[s]}$ denote the random tuple defined by
$h_{j}:=g_{i_{j}}$, which satisfies $Ph_{j}=Pg_{i_{j}}=Px_{i_{j}}=Py_{j}$, we
have

$\operatorname{SDP}(Y,k)\geq\frac{1}{s}\Big{(}\|PY\|_{F}^{2}-k\|PY\|_{2\to
2}^{2}\Big{)}\geq\frac{1}{s}\Big{(}\|PH\|_{F}^{2}-k\|H\|_{2\to 2}^{2}\Big{)}.$

Letting $G\in\mathbb{R}^{d\times n}$ denote the matrix whose $i$th column is
$g_{i}$, we take expectations of both sides to get

$\mathbb{E}_{X}\operatorname{SDP}(Y,k)\geq\frac{1}{s}\mathbb{E}_{G}\|PH\|_{F}^{2}-\frac{k}{s}\mathbb{E}_{G}\|H\|_{2\to
2}^{2}.$ (10)

The first term can be rewritten as

$\frac{1}{s}\mathbb{E}_{G}\|PH\|_{F}^{2}=\frac{1}{s}\mathbb{E}_{G}\sum_{j\in[s]}\|Ph_{j}\|^{2}=\frac{1}{s}\sum_{j\in[s]}\mathbb{E}_{G}\|Ph_{j}\|^{2}=\frac{1}{n}\sum_{i\in[n]}\|Pg_{i}\|^{2}=\frac{1}{n}\|PG\|_{F}^{2}.$

Furthermore, $\|PG\|_{F}^{2}$ has chi-squared distribution with
$(d-k^{\prime})n$ degrees of freedom, where $k^{\prime}$ denotes the dimension
of $\operatorname{span}\\{\mu_{t}\\}_{t\in[k]}$. Thus, Proposition 21 gives
$\frac{1}{n}\|PG\|_{F}^{2}\geq d-k^{\prime}-1$ with probability $\geq
1-e^{-n/(4d)}$. The second term in (10) is bounded by Lemma 25, and the result
follows from a union bound. ∎

###### Proof of Theorem 4(b).

Recall that Lemma 26 gives

$d-6k-1\leq\mathbb{E}_{X}\operatorname{SDP}(Y,k)\leq\operatorname{IP}(X,k)\leq
d+1$

with probability $\geq 1-\frac{1}{n}-e^{-\Omega(n/(d+\log n)^{2})}$. We will
use Hoeffding’s inequality to show that

$B_{H}\geq\mathbb{E}_{X}\operatorname{SDP}(Y,k)-1$

with probability $\geq 1-\epsilon$, from which the result follows by a union
bound. Our use of Hoeffding’s inequality requires a bound on $b$. To this end,
Lemma 6 gives that

$b\leq 4\operatorname{inf}(f)^{2}\leq 4f(\\{\mu_{t}\\}_{t\in[k]})^{2},$

which in turn is at most $4$ times the maximum of $n$ independent chi-squared
random variables with $d$ degrees of freedom. Corollary 23 bounds this
quantity by $4(2d+6\log n)$, and this bound holds in the event considered in
Lemma 26. Combining the bounds $b\leq 4(2d+6\log n)$ and $\ell\geq 128(d+3\log
n)^{2}\log(1/\epsilon)$ then gives
$t:=(\tfrac{b^{2}}{2\ell}\log(\tfrac{1}{\epsilon}))^{1/2}\leq\frac{1}{2}$, and
so $t\leq 1-t$. Hoeffding’s inequality then gives

$\displaystyle\mathbb{P}\\{B_{H}<\mathbb{E}_{X}\operatorname{SDP}(Y,k)-1\\}$
$\displaystyle=\mathbb{P}\bigg{\\{}\frac{1}{\ell}\sum_{i\in[\ell]}\operatorname{SDP}(Y_{i},k)-t<\mathbb{E}_{X}\operatorname{SDP}(Y,k)-1\bigg{\\}}$
$\displaystyle\leq\exp(-\tfrac{2\ell(1-t)^{2}}{b^{2}})\leq\exp(-\tfrac{2\ell
t^{2}}{b^{2}})=\epsilon,$

where the last step applies the definition of $t$. ∎

###### Proof of Theorem 5(b).

We have from the proof of the upper bound in Lemma 26 that
$\operatorname{IP}(X,k)\leq d+1$ with probability $\geq 1-e^{-n/(16d)}$. Thus,
it suffices to show that $B_{M}\geq d-3k-2$ with probability $\geq
1-e^{-s/(8d)}-2e^{-s/54}$. To this end, we first observe that

$d-3k-2\leq(1-\tfrac{1}{d})(d-3k-1)\qquad\text{and}\qquad\epsilon^{1/\ell}\geq
1-\tfrac{1}{d}.$ (11)

Indeed, the first inequality can be seen by expanding, while the second
inequality follows from the assumption $\ell\geq d\log(1/\epsilon)$. In
particular, $d\geq 1$ implies $\frac{1}{d}\leq\log(\frac{1}{1-1/d})$, and so

$\ell\geq
d\log(1/\epsilon)\geq\frac{\log(1/\epsilon)}{\log(\frac{1}{1-1/d})},$

and rearranging gives the desired inequality. We apply (11) with a union bound
to get

$\displaystyle\mathbb{P}\Big{\\{}B_{M}<d-3k-2\Big{\\}}$
$\displaystyle\leq\mathbb{P}\Big{\\{}\epsilon^{1/\ell}\min_{i\in[\ell]}\operatorname{SDP}(Y_{i},k)<(1-\tfrac{1}{d})(d-3k-1)\Big{\\}}$
$\displaystyle\leq\mathbb{P}\Big{\\{}\min_{i\in[\ell]}\operatorname{SDP}(Y_{i},k)<d-3k-1\Big{\\}}$
$\displaystyle\leq\ell\cdot\mathbb{P}\Big{\\{}\operatorname{SDP}(Y,k)<d-3k-1\Big{\\}},$

where $Y$ is a random matrix with the same distribution as each $Y_{i}$. In
particular, the columns of $Y$ are drawn uniformly without replacement from
$X$. Let $S:=\\{i_{j}\\}_{j\in[s]}$ denote the random indices such that the
$j$th column of $Y$ is $y_{j}=x_{i_{j}}$. By the law of total probability, it
suffices to bound the conditional probability

$\mathbb{P}_{S}\Big{\\{}\operatorname{SDP}(Y,k)<d-3k-1\Big{\\}}$

uniformly over all possible realizations of $S$. Recall that
$x_{i}=\mu_{t_{i}}+g_{i}$, where $t_{i}\sim\mathsf{Unif}([k])$ and
$g_{i}\sim\mathsf{N}(0,I_{d})$. Let $H\in\mathbb{R}^{d\times s}$ denote the
random matrix whose $j$th column is $h_{j}:=g_{i_{j}}$. Similar to the proof
of the lower bound in Lemma 26, we take $P$ to be the orthogonal projection
onto the $(d-k^{\prime})$-dimensional subspace
$(\operatorname{span}\\{\mu_{t}\\}_{t\in[k]})^{\perp}$ and then apply Lemma 18
to get

$\operatorname{SDP}(Y,k)\geq\frac{1}{s}\Big{(}\|PH\|_{F}^{2}-k\|H\|_{2\to
2}^{2}\Big{)}.$

This implies

$\displaystyle\mathbb{P}_{S}\Big{\\{}\operatorname{SDP}(Y,k)<d-3k-1\Big{\\}}$
$\displaystyle\leq\mathbb{P}_{S}\Big{\\{}\tfrac{1}{s}\big{(}\|PH\|_{F}^{2}-k\|H\|_{2\to
2}^{2}\big{)}<d-3k-1\Big{\\}}$
$\displaystyle\leq\mathbb{P}_{S}\Big{\\{}\tfrac{1}{s}\|PH\|_{F}^{2}<d-k-1\Big{\\}}+\mathbb{P}_{S}\Big{\\{}\tfrac{k}{s}\|H\|_{2\to
2}^{2}>2k\Big{\\}}.$ (12)

Conditioned on $S$, the entries of $H$ are independent with standard gaussian
distribution. As such, $\|PH\|_{F}^{2}$ has chi-squared distribution with
$(d-k^{\prime})s$ degrees of freedom, and so we apply the second part of
Proposition 21 with $t=s/(4(d-k^{\prime}))$ to bound the first term in (12) by
$e^{-s/(4d)}$. Next, we apply Proposition 24 with
$t=(\sqrt{2}-\sqrt{d/s}-1)\sqrt{s}$ to bound the second term in (12) by
$2e^{-t^{2}/2}$, which in turn is at most $2e^{-s/27}$ since $s\geq
54d\log\ell\geq 54d$ by assumption. Overall, we have

$\displaystyle\mathbb{P}\Big{\\{}B_{M}<d-3k-2\Big{\\}}$
$\displaystyle\leq\ell\cdot\mathbb{P}\Big{\\{}\operatorname{SDP}(Y,k)<d-3k-1\Big{\\}}$
$\displaystyle=\ell\cdot\mathbb{E}\Big{[}\mathbb{P}_{S}\Big{\\{}\operatorname{SDP}(Y,k)<d-3k-1\Big{\\}}\Big{]}$
$\displaystyle\leq e^{s/(54d)}\cdot(e^{-s/(4d)}+2e^{-s/27})\leq
e^{-s/(8d)}+2e^{-s/54},$

as desired. ∎

## Acknowledgments

Part of this research was conducted while SV was a Research Fellow at the
Simons Institute for Computing, University of California at Berkeley. DGM was
partially supported by NSF DMS 1829955, AFOSR FA9550-18-1-0107, and an AFOSR
Young Investigator Research Program award. SV is partially supported by ONR
N00014-22-1-2126, NSF CISE 2212457, an AI2AI Amazon research award, and the
NSF–Simons Research Collaboration on the Mathematical and Scientific
Foundations of Deep Learning (MoDL) (NSF DMS 2031985).

## References

  * [1] E. Abbe, A. S. Bandeira, G. Hall, Exact recovery in the stochastic block model, IEEE Trans. Inform. Theory 62 (2015) 471–487.
  * [2] P. Abdalla, A. S. Bandeira, Community detection with a subsampled semidefinite program. Sampl. Theory Signal Process. Data Anal. 20 (2022) 1–10.
  * [3] D. Aloise, A. Deshpande, P. Hansen, P. Popat, NP-hardness of euclidean sum-of-squares clustering, Mach. Learn. 75 (2009) 245–248.
  * [4] P. Awasthi, A. S. Bandeira, M. Charikar, R. Krishnaswamy, S. Villar, R. Ward, Relax, no need to round: Integrality of clustering formulations, ITCS 2015, 191–200.
  * [5] P. Awasthi, M. Charikar, R. Krishnaswamy, A. K. Sinop, The hardness of approximation of euclidean k-means, arXiv:1502.03316
  * [6] A. S. Bandeira, A note on probably certifiably correct algorithms. C. R. Math. 354 (2016) 329–333.
  * [7] A. S. Bandeira, Random Laplacian matrices and convex relaxations, Found. Comput. Math. 18 (2018) 345–379.
  * [8] P. Collard, Cloud data set, archive.ics.uci.edu/ml/datasets/cloud
  * [9] S. Dasgupta, Learning mixtures of Gaussians, FOCS 1999, 634–644.
  * [10] D. Davis, M. Diaz, K. Wang, Clustering a mixture of gaussians with unknown covariance, arXiv preprint arXiv:2110.01602
  * [11] A. De Rosa, A. Khajavirad, The ratio-cut polytope and K-means clustering, SIAM J. Optim. 32 (2022) 173–203.
  * [12] I. Diakonikolas, D. M. Kane, A. Stewart, Statistical query lower bounds for robust estimation of high-dimensional gaussians and gaussian mixtures, FOCS 2017, 73–84.
  * [13] Y. Fei, Y. Chen, Hidden integrality of SDP relaxations for sub-Gaussian mixture models, COLT 2018, 1931–1965.
  * [14] C. Giraud, N. Verzelen, 2019\. Partial recovery bounds for clustering with the relaxed $K$-means, Math. Stat. Learn. 1 (2019) 317–374.
  * [15] M. X. Goemans, D. P. Williamson, Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming, JACM 42 (1995) 1115–1145.
  * [16] M. Grant, S. Boyd, CVX: Matlab software for disciplined convex programming, cvxr.com/cvx
  * [17] O. Guédon, R. Vershynin, Community detection in sparse networks via Grothendieck’s inequality, Probab. Theory Relat. Fields 165 (2016) 1025–1049.
  * [18] D. Hsu, S. M. Kakade, Learning mixtures of spherical gaussians: moment methods and spectral decompositions, ITCS 2013, 11–20.
  * [19] T. Iguchi, D. G. Mixon, J. Peterson, S. Villar, Probably certifiably correct k-means clustering, Math. Program. 165 (2017) 605–642.
  * [20] R. Kannan, H. Salmasian, S. Vempala, The spectral method for general mixture models, COLT 2005, 444–457.
  * [21] R. M. Karp, Reducibility among combinatorial problems, In: Complexity of Computer Computations, Springer, Boston, MA, 1972, pp. 85–103.
  * [22] KDD Cup 1999 dataset, kdd.ics.uci.edu/databases/kddcup99/kddcup99.html
  * [23] M. R. Kosorok, Introduction to empirical processes and semiparametric inference, Springer, 2008.
  * [24] A. Kumar, R. Kannan, Clustering with spectral norm and the k-means algorithm, FOCS 2010, 299–308.
  * [25] B. Laurent, P. Massart, Adaptive estimation of a quadratic functional by model selection, Ann. Stat. 28 (2000) 1302–1338.
  * [26] Y. LeCun, C. Cortes, MNIST handwritten digit database, AT&T Labs, yann.lecun.com/exdb/mnist
  * [27] X. Li, Y. Li, S. Ling, T. Strohmer, K. Wei, When do birds of a feather flock together? k-means, proximity, and conic programming, Math. Program. 179 (2020) 295–341.
  * [28] S. Ling, T. Strohmer, Certifying global optimality of graph cuts via semidefinite relaxation: A performance guarantee for spectral clustering, Found. Comput. Math. 20 (2020) 367–421.
  * [29] S. Lloyd, Least squares quantization in PCM, IEEE Trans. Inform. Theory 28 (1982) 129–137.
  * [30] D. G. Mixon, S. Villar, Monte Carlo approximation certificates for k-means clustering, arXiv:1710.00956
  * [31] D. G. Mixon, S. Villar, R. Ward, Clustering subgaussian mixtures by semidefinite programming, Inform. Inference 6 (2017) 389–415.
  * [32] D. G. Mixon, K. Xie, Sketching Semidefinite Programs for Faster Clustering, IEEE Trans. Inform. Theory 67 (2021) 6832–6840.
  * [33] A. Nellore, R. Ward, Recovery guarantees for exemplar-based clustering, Inform. Comput. 245 (2015) 165–180.
  * [34] Y. Nesterov, A. Nemirovskii, Interior-point polynomial algorithms in convex programming, SIAM, 1994.
  * [35] J. Peng, Y. Wei, Approximating k-means-type clustering via semidefinite programming, SIAM J. Optim. 18 (2007) 186–205.
  * [36] V. Piccialli, A. M. Sudoso, A. Wiegele, SOS-SDP: an exact solver for minimum sum-of-squares clustering, INFORMS J. Comput., 2022.
  * [37] M. N. Prasad, G. A. Hanasusanto, Improved conic reformulations for k-means clustering, SIAM J. Optim. 28 (2018) 3105–3126.
  * [38] D. M. Rosen, L. Carlone, A. S. Bandeira, J. J. Leonard, SE-Sync: A certifiably correct algorithm for synchronization over the special Euclidean group, Int. J. Robot. Res. 38 (2019) 95–125.
  * [39] D. F. Sun, L. Q. Yang, K. C. Toh, Sdpnal+: A majorized semismooth newton-cg augmented lagrangian method for semidefinite programming with nonnegative constraints, Math. Program. Comput. (2015) 331–366.
  * [40] J. A. Tropp, An Introduction to Matrix Concentration Inequalities, Found. Trends Mach. Learn. 8 (2015) 1–230.
  * [41] S. Vassilvitskii, D. Arthur, k-means++: The advantages of careful seeding, SODA 2006, 1027–1035.
  * [42] R. Vershynin, Introduction to the non-asymptotic analysis of random matrices, Compressed Sensing, Theory and Applications, Cambridge U. Press, 2012.
  * [43] D. P. Woodruff, Sketching as a tool for numerical linear algebra, Found. Trends Theor. Comp. Sci. 10 (2014) 1–157.
  * [44] Y. Wu, P. Yang, Optimal estimation of Gaussian mixtures via denoised method of moments, Ann. Stat. 48 (2020) 1981–2007.
  * [45] H. Yang, J. Shi, L. Carlone, Teaser: Fast and certifiable point cloud registration, IEEE Trans. Robot. 37 (2020) 314–333.
  * [46] Y. Zhuang, X. Chen, Y. Yang, Sketch-and-lift: scalable subsampled semidefinite program for K-means clustering, PMLR (2022) 9214–9246.