# Do you know what $q$-means?

João F. Doriguello<EMAIL_ADDRESS>Centre for Quantum Technologies, National
University of Singapore, Singapore  Alessandro Luongo<EMAIL_ADDRESS>Centre
for Quantum Technologies, National University of Singapore, Singapore
Inveriant Pte. Ltd.  Ewin Tang<EMAIL_ADDRESS>Paul G. Allen School of
Computer Science & Engineering, University of Washington

###### Abstract

Clustering is one of the most important tools for analysis of large datasets,
and perhaps the most popular clustering algorithm is Lloyd’s iteration for
$k$-means. This iteration takes $N$ vectors
$v_{1},\dots,v_{N}\in\mathbb{R}^{d}$ and outputs $k$ centroids
$c_{1},\dots,c_{k}\in\mathbb{R}^{d}$; these partition the vectors into
clusters based on which centroid is closest to a particular vector. We present
an overall improved version of the “$q$-means” algorithm, the quantum
algorithm originally proposed by Kerenidis, Landman, Luongo, and Prakash which
performs $\varepsilon$-$k$-means, an approximate version of $k$-means
clustering [KLLP19]. This algorithm does not rely on the quantum linear
algebra primitives of prior work, instead only using its QRAM to prepare and
measure simple states based on the current iteration’s clusters. The time
complexity is
$O\big{(}\frac{k^{2}}{\varepsilon^{2}}(\sqrt{k}d+\log(Nd))\big{)}$ and
maintains the polylogarithmic dependence on $N$ while improving the dependence
on most of the other parameters. We also present a “dequantized” algorithm for
$\varepsilon$-$k$-means which runs in
$O\big{(}\frac{k^{2}}{\varepsilon^{2}}(kd+\log(Nd))\big{)}$ time. Notably,
this classical algorithm matches the polylogarithmic dependence on $N$
attained by the quantum algorithms.

##  1 Introduction

Among machine learning problems, data clustering and the $k$-means problem are
of particular relevance and have attracted much attention in the past [HW79,
KM99, LVV03]. Here the task is to find an assignment of each vector from a
dataset of size $N$ to one of $k$ labels (for a given $k$ assumed to be known)
such that similar vectors are assigned to the same cluster. To be more
precise, in the $k$-means problem we are given $N$ vectors
$v_{1},\dots,v_{N}\in\mathbb{R}^{d}$, given as columns in a matrix
$V\in\mathbb{R}^{N\times d}$, and a positive integer $k$, and the task is to
find $k$ centroids $c_{1},\dots,c_{k}\in\mathbb{R}^{d}$ such that the cost
function

$\displaystyle\sum_{i=1}^{N}\min_{j\in[k]}\|v_{i}-c_{j}\|^{2}$

is minimized, where $\|v_{i}-c_{j}\|$ is the Euclidean distance between
$v_{i}$ and $c_{j}$ and $[n]:=\\{1,\dots,n\\}$ for $n\in\mathbb{N}$. Since the
$k$-means problem is known to be NP-hard [Das08, Vat09, MNV12], several
polynomial-time algorithms [KMN+02, JKS14, ANFSW19, BGJK20] have been
developed to obtain _approximate_ solutions. One such algorithm is the
$k$-means algorithm, or Lloyd’s iterations, introduced by Lloyd [Llo82]. The
$k$-means algorithm is a heuristic algorithm that iteratively updates the
centroids $c_{1},\dots,c_{k}$ until some desired precision is reached. At each
time step $t$, the algorithm clusters the data in $k$ clusters denoted by the
sets $\mathcal{C}_{j}^{t}\subseteq[N]$, $j\in[k]$, each with centroid
$c_{j}^{t}$, and then updates the centroids based on such clustering. More
precisely, the $k$-means algorithm starts with initial centroids
$c_{1}^{0},\dots,c_{k}^{0}\in\mathbb{R}^{d}$, picked either randomly or
through some pre-processing routine as in the $k$-means++ algorithm [AV07],
and alternates between two steps:

  * •

Each vector $v_{i}$, for $i\in[N]$, is assigned a cluster
$\mathcal{C}_{\ell_{i}^{t}}^{t}$ where
$\ell^{t}_{i}=\arg\min_{j\in[k]}\|v_{i}-c_{j}^{t}\|$;

  * •

The new centroids $\\{c_{j}^{t+1}\\}_{j=1}^{k}$ are updated based on the new
clusters,
$c_{j}^{t+1}=\frac{1}{|\mathcal{C}_{j}^{t}|}\sum_{i\in\mathcal{C}_{j}^{t}}v_{i}$.

The above steps are repeated until
$\frac{1}{k}\sum_{j=1}^{k}\|c_{j}^{t}-c_{j}^{t+1}\|\leq\tau$ for a given
threshold $\tau>0$, in which case we say that the algorithm has converged.

On the other hand, quantum computation is a fast-paced field that promises a
new paradigm to solve certain tasks more efficiently, well-known examples
being Shor’s algorithm [Sho94, Sho99] and Grover’s search [Gro96, Gro97]. The
subfield of quantum machine learning aims to offer computational advantages in
machine learning, with many algorithms proposed to attain these advantages
[LMR14, KP16, KP20a, CGJ18, LMR13, AHKZ20]. A notable line of work in this
subfield is that of quantum versions of the $k$-means algorithm, to improve
the classical running time of $O(kNd)$ per-iteration. Aimeur, Brassard, and
Gambs [ABG13] proposed a $O(N^{3/2}/\sqrt{k})$-time quantum algorithm based on
quantum minimum finding [DH96] for executing a single Lloyd’s iteration for
the $k$-median clustering problem (which is related to the $k$-means problem).
Lloyd, Mohseni, and Rebentrost [LMR13] gave a
$O\big{(}\frac{kN\log{d}}{\varepsilon}\big{)}$-time quantum algorithm for
executing a single Lloyd’s iteration from the $k$-means problem. The authors
used an efficient subroutine for quantum distance estimation with error
$\varepsilon$ assuming quantum access to the dataset via a QRAM (Quantum
Random Access Memory) [GLM08a, GLM08b].

More recently, Kerenidis, Landman, Luongo, and Prakash [KLLP19] proposed a
quantum version of Lloyd’s iterations and $k$-means++, which they named
$q$-means, using quantum linear algebra subroutines and assuming the input
vectors be stored in a QRAM. Interestingly enough, their runtime only depends
polylogarithmically on the size $N$ of the dataset. On the other hand,
however, their quantum algorithm only performs the iteration approximately and
does not give provable approximation guarantees, since Lloyd’s iterations also
do not. Instead, the authors show that $q$-means performs a robust version of
$k$-means which they name $\varepsilon$-$k$-means, and then show through
experiments that the approximation does not worsen the quality of the
centroids. In the $\varepsilon$-$k$-means algorithm, a noise parameter
$\varepsilon>0$ is introduced both in the distance estimation and centroid
update steps. It alternates now between the steps:

  * •

Each vector $v_{i}$, for $i\in[N]$, is assigned a cluster
$\mathcal{C}_{\ell_{i}^{t}}^{t}$ such that
$\ell_{i}^{t}\in\\{l:|\|v_{i}-c_{i}^{\ast}\|-\|v_{i}-c_{l}^{t}\||\leq\varepsilon\\}$,
where $c^{\ast}_{i}$ is the closet centroid to $v_{i}$;

  * •

The new centroids $\\{c_{j}^{t+1}\\}_{j=1}^{k}$ are updated such that
$\big{\|}c_{j}^{t+1}-\frac{1}{|\mathcal{C}_{j}^{t}|}\sum_{i\in\mathcal{C}_{j}^{t}}v_{i}\big{\|}\leq\varepsilon$.

Again, it is iterated until convergence.

The overall idea behind $q$-means as proposed in [KLLP19] is to first create
the quantum state

$\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{N}|i\rangle|\ell_{i}^{t}\rangle$

using distance estimation and quantum minimum finding [DH96], followed by
measuring the label register $|\ell_{i}^{t}\rangle$ to obtain

$\displaystyle|\chi_{j}^{t}\rangle=\frac{1}{\sqrt{|\mathcal{C}_{j}^{t}|}}\sum_{i\in\mathcal{C}_{j}^{t}}|i\rangle$

for some random $j\in[k]$. Here $|\chi_{j}^{t}\rangle$ is the quantum state of
the characteristic vector $\chi_{j}^{t}\in\mathbb{R}^{N}$ defined as
$(\chi_{j}^{t})_{i}=\frac{1}{|\mathcal{C}_{j}^{t}|}$ if
$i\in\mathcal{C}^{t}_{j}$ and $0$ if $i\not\in\mathcal{C}^{t}_{j}$. The
$q$-means algorithm then proceeds to perform a matrix multiplication with
matrix $V^{\top}$ and the quantum state $|\chi_{j}^{t}\rangle$ by using
quantum linear algebra techniques. This leads to

$\displaystyle|c_{j}^{t+1}\rangle=\frac{1}{\sqrt{|\mathcal{C}_{j}^{t}|}}\sum_{i\in\mathcal{C}_{j}^{t}}|v_{i}\rangle$

by the crucial fact that $c_{j}^{t+1}=V^{\top}\chi_{j}^{t}$. Finally, quantum
tomography [KP20b] is performed on $|c^{t+1}\rangle$ in order to retrieve a
classical description. Since quantum linear algebra techniques are used to
obtain $|c_{j}^{t+1}\rangle$ from $|\chi_{j}^{t}\rangle$, the final runtime
depends on quantities like the condition number $\kappa(V)$ of the matrix $V$
and the matrix dependent parameter $\mu(V)$ defined as

$\displaystyle\mu(V)=\min_{p\in[0,1]}\left(\|V\|_{F},\sqrt{s_{2p}(V)s_{1-2p}(V^{\top})}\right),$
(1)

where $\|V\|_{F}:=\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{d}V_{ij}^{2}}$ is the
Frobenius norm of $V$ and $s_{p}(V):=\max_{i\in[N]}\sum_{j=1}^{d}V_{ij}^{p}$.
Their final runtime is the following.

######  Theorem 1 ([KLLP19, Theorem 3.1]).

For a parameter $\varepsilon>0$ and data matrix $V\in\mathbb{R}^{N\times d}$
with $1\leq\|v_{i}\|^{2}\leq\eta$ and condition number $\kappa(V)$, the
$q$-means algorithm outputs centroids consistent with the classical
$\varepsilon$-$k$-means algorithm in
$\widetilde{O}\big{(}kd\frac{\eta}{\varepsilon^{2}}\kappa(V)(\mu(V)+k\frac{\eta}{\varepsilon})+k^{2}\frac{\eta^{3/2}}{\varepsilon^{2}}\kappa(V)\mu(V)\big{)}$
time per iteration, where $\mu(V)$ is defined in Eq. (1).

###  1.1 Our results

In this work, we provide an improved $q$-means algorithm that avoids the need
for quantum linear algebra subroutines and still keeps the polylogarithmically
dependence on the size $N$ of the dataset. The main idea is to sample from the
quantum states
$\frac{1}{\|V\|_{F}}\sum_{i=1}^{N}\|v_{i}\||i\rangle|\ell_{i}^{t}\rangle$,
which can be efficiently constructed assuming quantum access to the matrix $V$
via a QRAM, and then to classically define the new centroids $c_{j}^{t+1}$
according to the sampled points. Our approach contrasts with the one from
[KLLP19] by (i) creating a slightly different state than
$\frac{1}{\sqrt{N}}\sum_{i=1}^{N}|i\rangle|\ell_{i}^{t}\rangle$ and (ii)
avoiding the quantum states $|c_{j}^{t+1}\rangle$ all together as an
intermediary step to obtaining the classical centroids $c_{j}^{t+1}$.
Furthermore, in the same spirit of previous “dequantization” results [Tan19,
Tan21, GLT18, GST22], we provide a classical algorithm for
$\varepsilon$-$k$-means based on $\ell_{2}$-sampling that has a very similar
runtime compared to our $q$-means, which includes the polylogarithmically
dependence on $N$. Here we assume access to a data structure that allows
efficient $\ell_{2}$-sampling from the rows of $V$ and from the vector of its
row norms.

######  Theorem 2.

Let $V\in\mathbb{R}^{N\times d}$ with $\min_{i\in[N]}\|v_{i}\|\geq 1$ and a
parameter $\varepsilon>0$. Then the $\varepsilon$-$k$-means algorithm can be
executed with probability at least $1-\delta$ and in time per iteration

$\displaystyle
O\left(\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}\left(kd+\log(Nd)\right)\log\frac{k}{\delta}\right).$

Moreover, there is a quantum algorithm that outputs centroids consistent with
the classical $\varepsilon$-$k$-means with probability at least $1-\delta$ and
in time per iteration

$\displaystyle
O\left(\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}(\sqrt{k}d+\log(Nd))\log\frac{k\|V\|_{F}}{\delta
N}\right).$

###  1.2 Future work

Even though the time complexity of $q$-means has been overall improved from
[KLLP19] to Theorem 2 by avoiding quantum linear algebra primitives, we
believe that it can be further improved. In particularly, we conjecture that
the time dependence on $\varepsilon$ can be decreased to $O(1/\varepsilon)$,
thus setting apart the current classical and quantum complexities more
substantially. Another open question is to consider a quantum version of
$k$-means++ using similar techniques employed here or assume data vectors with
special structure properties, e.g. well-clusterable datasets as in [KLLP19],
in order to obtain tighter runtimes.

Finally, we mention the independent and very recent work of Jaiswal [Jai23],
who quantized the highly parallel, sampling-based approximation scheme of
[BGJK20] and thus obtained a quantum algorithm for the $k$-means problem with
_provable_ guarantees, as opposed to our results, which are heuristic. More
specifically, the quantum algorithm of [Jai23] outputs a centroid set that
minimizes the $k$-means cost function up to relative error $\varepsilon$. Due
to such guarantee, though, their final runtime depends exponentially in $k$
and $1/\varepsilon$ (but maintains the polylogarithmically dependence on $N$).
It would be interesting to bridge both our results and obtain time
complexities polynomial in $k$ and $1/\varepsilon$ together with provable
guarantees. This will probably require stronger assumptions on the input.

#### Acknowledgements.

ET thanks Haotian Jiang for helping with the analysis. ET is supported by the
NSF GRFP (DGE-1762114). JFD and AL are supported by the National Research
Foundation, Singapore and A*STAR under its CQT Bridging Grant and its Quantum
Engineering Programme under grant NRF2021-QEP2-02-P05.

##  2 Preliminaries

Define $[n]:=\\{1,\dots,n\\}$. In the rest of the paper, our dataset is a
matrix $V\in\mathbb{R}^{N\times d}$ comprised of vectors
$v_{i}\in\mathbb{R}^{d}$, $i\in[N]$, such that $\min_{i\in[N]}\|v_{i}\|\geq
1$, where $N$ is the number of samples and $d$ is the number of features.
Given a vector $x\in\mathbb{R}^{d}$, let
$\|x\|:=\sqrt{\sum_{i=1}^{d}x_{i}^{2}}$ and by $\mathcal{D}_{x}$ we denote the
distribution over $[d]$ with probability density function
$\mathcal{D}_{x}(i)=x_{i}^{2}/\|x\|^{2}$. For $V\in\mathbb{R}^{N\times d}$,
let $\|V\|_{F}:=\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{d}|V_{ij}|^{2}}$ be its
Frobenius norm. In the next sections, we will make use of the quantum minimum
finding subroutine.

######  Fact 3 (Minimum finding [DH96]).

Given $\delta\in(0,1)$ and oracle access to $v\in\mathbb{R}^{d}$, there is a
quantum algorithm that outputs ${\rm argmin}_{i\in[n]}(v_{i})$ with
probability at least $1-\delta$ in $O(\sqrt{d}\log(1/\delta))$ time.

We will also need the following concentration inequality.

######  Fact 4 (McDiarmid’s Inequality (Bernstein form)).

Let $f(X_{1},\dots,X_{n})$ be a function of independent random variables
$X_{1},\dots,X_{n}$ such that, for all $i\in[n]$ and choices of
$X_{1},\dots,X_{i-1},X_{i+1},\dots,X_{n}$,

$\displaystyle\operatorname*{\max}_{X_{i},X^{\prime}_{i}}|f(X_{1},\dots,X_{i},\dots,X_{n})-f(X_{1},\dots,X^{\prime}_{i},\dots,X_{n})|\leq\delta,$
$\displaystyle\operatorname*{\mathbb{E}}_{X_{i}}[f(X_{1},\dots,X_{n})^{2}]-\operatorname*{\mathbb{E}}_{X_{i}}[f(X_{1},\dots,X_{n})]^{2}\leq\sigma^{2}.$

Then, for all $t>0$,

$\displaystyle\operatorname{Pr}[f(X_{1},\dots,X_{n})\geq
t]\leq\exp\left(-\frac{t^{2}/2}{n\sigma^{2}+\delta t/3}\right).$

#### Data structures.

In our classical algorithm, we shall assume the existence of a low-overhead
data structure storing the dataset matrix $V\in\mathbb{R}^{N\times d}$.

######  Lemma 5.

Given a matrix $V\in\mathbb{R}^{N\times d}$ with $w$ non-zero entries and the
vector $\widetilde{V}\in\mathbb{R}^{N}$ whose $i$-th entry is $\|v_{i}\|$,
there is a data structure that stores $V$ in $O(w\log(Nd))$ space and supports
the operations:

  1. 1.

Reading and updating an entry of $V$ in $O(\log(Nd))$ time;

  2. 2.

Finding $\widetilde{V}_{i}$ in $O(\log{N})$ time;

  3. 3.

Finding $\|V\|_{F}^{2}$ in $O(1)$ time;

  4. 4.

Sampling from $\mathcal{D}_{\widetilde{V}}$ and $\mathcal{D}_{v_{i}}$ in
$O(\log(Nd))$ time.

In our quantum algorithm, we shall assume quantum access via a QRAM to a KP-
tree data structure [Pra14, KP16] that stores the dataset matrix
$V\in\mathbb{R}^{N\times d}$.

######  Lemma 6 (KP-tree stored in QRAM [KP16]).

Given a matrix $V\in\mathbb{R}^{N\times d}$ with $w$ non-zero entries, there
is a data structure that stores $V$ in $O(w\log^{2}(Nd))$ space such that

  1. 1.

Reading and updating an entry of $V$ takes $O(\log^{2}(Nd))$ time;

  2. 2.

A quantum algorithm with access to the date structure can perform the
following operations in $O(\log(Nd))$ time:

     * •

$|i\rangle|0\rangle\mapsto|i\rangle|v_{i}\rangle$ for $i\in[N]$;

     * •

$|0\rangle\mapsto\frac{1}{\|V\|_{F}}\sum_{i=1}^{N}\|v_{i}\||i\rangle$.

##  3 Classical $\varepsilon$-$k$-means algorithm

Algorithm 1 $\varepsilon$-$k$-means algorithm

1:Data matrix $V\in\mathbb{R}^{N\times d}$ with $\min_{i\in[N]}\|v_{i}\|\geq
1$ stored in the data structure from Lemma 5. Precision parameters
$\delta,\varepsilon$.

2:Vectors $c_{1},\dots,c_{k}\in\mathbb{R}^{d}$ corresponding to centroids.

3:Select $k$ initial centroids $c_{1}^{0},\dots,c_{k}^{0}$ and store them in a
RAM data structure.

4:for $t=0$ until convergence do

5: Sample
$p=O\big{(}\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}\log\frac{k}{\delta}\big{)}$
rows $P=(s_{1},\dots,s_{p})\in[N]^{p}$ uniformly from $[N]$.

6: Sample
$q=O\big{(}\frac{\|V\|_{F}^{2}}{N}\frac{k}{\varepsilon^{2}}\log\frac{k}{\delta}\big{)}$
rows $Q=(r_{1},\dots,r_{q})\in[N]^{q}$ from $\mathcal{D}_{\widetilde{V}}$.

7: For $i\in[p]$, label
$\ell_{i}^{(P)}=\arg\min_{j\in[k]}\|v_{s_{i}}-c_{j}^{t}\|$. Let
$P_{j}:=\\{i\in P|\ell_{i}^{(P)}=j\\}$ for $j\in[k]$.

8: For $i\in[q]$, label
$\ell_{i}^{(Q)}=\arg\min_{j\in[k]}\|v_{r_{i}}-c_{j}^{t}\|$. Let
$Q_{j}:=\\{i\in Q|\ell_{i}^{(Q)}=j\\}$ for $j\in[k]$.

9: For $j\in[k]$, record the new centroids:
$\displaystyle\widetilde{c}_{j}^{t+1}=\frac{p}{N|P_{j}|}\sum_{i\in
Q_{j}}\frac{\|V\|_{F}^{2}}{q\|v_{r_{i}}\|^{2}}v_{r_{i}}.$

10:end for

######  Theorem 7.

Given a data matrix $V\in\mathbb{R}^{N\times d}$ with
$\min_{i\in[N]}\|v_{i}\|\geq 1$ stored in the date structure from Lemma 5 and
a parameter $\varepsilon>0$, Algorithm 1 outputs centroids consistent with the
$\varepsilon$-$k$-means algorithm with probability at least $1-\delta$ in time
per iteration

$\displaystyle
O\left(\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}\left(kd+\log(Nd)\right)\log\frac{k}{\delta}\right).$

###### Proof.

We start with the error analysis. Let $\chi_{j}^{t}\in\mathbb{R}^{N}$ be the
characteristic vector for cluster $j\in[k]$ at iteration $t$ scaled to
$\ell_{1}$-norm, i.e., $(\chi_{j}^{t})_{i}=\frac{1}{|\mathcal{C}_{j}^{t}|}$ if
$i\in\mathcal{C}^{t}_{j}$ and $0$ if $i\not\in\mathcal{C}^{t}_{j}$. Let
$X_{1},\dots,X_{q}\in\mathbb{R}^{N\times d}$ be i.i.d. random matrices formed
by sampling a row $i$ from $V$ (with probability
$\|v_{i}\|^{2}/\|V\|_{F}^{2}$), setting the $i$-th row of $X$ to
$\|V\|_{F}^{2}\frac{v_{i}}{\|v_{i}\|^{2}}$ and the rest to zero. Define
$\widehat{V}=\frac{1}{q}\sum_{j=1}^{q}X_{j}$. Then
$\mathbb{E}[\widehat{V}]=V$. Finally, let
$\lambda_{j}=\frac{p}{|P_{j}|}\frac{|\mathcal{C}_{j}^{t}|}{N}$.

We note that the standard $k$-means and Algorithm 1 can be stated,
respectively, as $c_{j}^{t+1}=V^{\top}\chi_{j}^{t}$ and
$\widetilde{c}_{j}^{t+1}=\lambda_{j}\widehat{V}^{\top}\chi_{j}^{t}$. Therefore

$\displaystyle\|\widetilde{c}_{j}^{t+1}-c_{j}^{t+1}\|=\|\lambda_{j}\widehat{V}^{\top}\chi_{j}^{t}-\lambda_{j}V^{\top}\chi_{j}^{t}+\lambda_{j}V^{\top}\chi_{j}^{t}-V^{\top}\chi_{j}^{t}\|\leq|\lambda_{j}-1|\|V^{\top}\chi_{j}^{t}\|+|\lambda_{j}|\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|,$

so we bound $|\lambda_{j}-1|\leq\frac{\varepsilon}{2\|V^{\top}\chi_{j}^{t}\|}$
and
$\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|\leq\frac{\varepsilon}{4}\leq\frac{\varepsilon}{2|\lambda_{j}|}$.
Let us start with bounding $|\lambda_{j}-1|$. Notice that $|P_{j}|$ is a
binomial random variable with mean $p|\mathcal{C}_{j}^{t}|/N$. By a Chernoff
bound,

$\displaystyle\operatorname{Pr}\left[\left|\frac{|P_{j}|}{p}\frac{N}{|\mathcal{C}_{j}^{t}|}-1\right|\geq\varepsilon_{0}\right]\leq
2e^{-\varepsilon_{0}^{2}p|\mathcal{C}_{j}^{t}|/3N}.$

It suffices to take
$p=O\big{(}\frac{\|V\|_{F}^{2}N}{\varepsilon^{2}|\mathcal{C}_{j}^{t}|^{2}}\log\frac{k}{\delta}\big{)}$
in order to estimate
$|\lambda_{j}^{-1}-1|\leq\frac{\varepsilon}{2\|V\|_{F}\|\chi_{j}^{t}\|}$, and
so $|\lambda_{j}-1|\leq\frac{\varepsilon}{\|V\|_{F}\|\chi_{j}^{t}\|}$, with
probability at least $1-\delta/2k$ (using that
$\|V^{\top}\chi_{j}^{t}\|\leq\|V\|_{F}\|\chi_{j}^{t}\|$ and
$\|\chi_{j}^{t}\|^{2}=1/|\mathcal{C}_{j}^{t}|$). By a union bound, the bound
on $\lambda_{j}$ holds for all $j\in[k]$ with probability at least
$1-\delta/2$.

Regarding the bound on $\widehat{V}$, we shall use McDiarmid’s inequality to
prove $\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|\leq\frac{\varepsilon}{4}$. For
such, let $f(X_{1},\dots,X_{q})=\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|$. We
first bound, for all $i\in[q]$,

$\displaystyle|f(X_{1},\dots,X_{i},\dots,X_{q})-f(X_{1},\dots,X^{\prime}_{i},\dots,X_{q})|$
$\displaystyle\leq\|(\widehat{V}-\widehat{V}^{\prime})^{\top}\chi_{j}^{t}\|$
$\displaystyle=\frac{1}{q}\|(X_{i}-X^{\prime}_{i})^{\top}\chi_{j}^{t}\|$
$\displaystyle\leq\frac{2}{q|\mathcal{C}_{j}^{t}|}\operatorname*{\max}_{u\in[N]}\frac{\|V\|_{F}^{2}}{\|v_{u}\|^{2}}\|v_{u}\|$
$\displaystyle\leq\frac{2\|V\|_{F}^{2}}{q|\mathcal{C}_{j}^{t}|},$

using that $\min_{i\in[N]}\|v_{i}\|\geq 1$. Second, we bound the variance: for
all $i\in[q]$,

$\displaystyle{\operatorname*{\mathbb{E}}_{X_{i}}}\bigg{[}{\operatorname*{\mathbb{E}}_{X^{\prime}_{i}}}\big{[}f(X_{1},\dots,X_{i},\dots,X_{q})-f(X_{1},\dots,X^{\prime}_{i},\dots,X_{q})\big{]}^{2}\bigg{]}$
$\displaystyle\leq\frac{1}{q^{2}}{\operatorname*{\mathbb{E}}_{X_{i}}}\bigg{[}{\operatorname*{\mathbb{E}}_{X^{\prime}_{i}}}\big{[}\|(X_{i}-X^{\prime}_{i})^{\top}\chi_{j}^{t}\|\big{]}^{2}\bigg{]}$
$\displaystyle\leq\frac{1}{q^{2}}{\operatorname*{\mathbb{E}}_{X_{i}}}\bigg{[}{\operatorname*{\mathbb{E}}_{X^{\prime}_{i}}}\big{[}\|X_{i}^{\top}\chi_{j}^{t}\|+\|(X_{i}^{\prime})^{\top}\chi_{j}^{t}\|\big{]}^{2}\bigg{]}$
$\displaystyle\leq\frac{4}{q^{2}}{\operatorname*{\mathbb{E}}_{X_{i}}}[\|X_{i}^{\top}\chi_{j}^{t}\|^{2}]$
$\displaystyle=\frac{4}{q^{2}}(\chi_{j}^{t})^{\top}\operatorname*{\mathbb{E}}_{X}[XX^{\top}]\chi_{j}^{t}$
$\displaystyle=\frac{4}{q^{2}}(\chi_{j}^{t})^{\top}(\|V\|_{F}^{2}\mathbb{I})\chi_{j}^{t}$
$\displaystyle=\frac{4\|V\|_{F}^{2}}{q^{2}|\mathcal{C}_{j}^{t}|}.$

McDiarmid’s inequality then leads to

$\displaystyle\operatorname{Pr}[\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|\geq\varepsilon/4]\leq\exp\left(-\frac{\varepsilon^{2}/32}{\frac{4\|V\|_{F}^{2}}{q|\mathcal{C}_{j}^{t}|}+\frac{\varepsilon\|V\|_{F}^{2}}{6q|\mathcal{C}_{j}^{t}|}}\right)\leq\exp\left(-\frac{q\varepsilon^{2}|\mathcal{C}_{j}^{t}|}{(128+6\varepsilon)\|V\|_{F}^{2}}\right).$

It suffices to take
$q=O\big{(}\frac{\|V\|_{F}^{2}}{|\mathcal{C}_{j}^{t}|}\frac{1}{\varepsilon^{2}}\log\frac{k}{\delta}\big{)}$
in order to approximate
$\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|\leq\varepsilon/4$ with probability at
least $1-\delta/2k$. By a union bound,
$\|(\widehat{V}-V)^{\top}\chi_{j}^{t}\|\leq\varepsilon/4$ for all $j\in[k]$
with probability at least $1-\delta/2$. All in all, we have
$\|\widetilde{c}_{j}^{t+1}-c_{j}^{t+1}\|\leq\varepsilon$ with probability at
least $1-\delta$.

Finally, the runtime of each step is:

  * •

Sampling rows $P,Q$ and querying the corresponding vectors takes
$O((p+q)\log(Nd))$ time;

  * •

Obtaining the labels $\ell_{i}^{(P)}$ and $\ell_{j}^{(Q)}$ for $i\in[p]$ and
$j\in[q]$ requires $O((p+q)kd)$ time;

  * •

Computing all $k$ new centroids by adding $q$ $d$-dimensional vectors takes
$O(qkd)$ time.

This means the overall complexity is (already using that
$|\mathcal{C}_{j}^{t}|/N=\Omega(1/k)$)

$O\big{(}(p+q)(kd+\log(Nd))\big{)}=O\left(\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}\left(kd+\log(Nd)\right)\log\frac{k}{\delta}\right).\qed$

##  4 Quantum $q$-means

Algorithm 2 $q$-means algorithm

1:Data matrix $V\in\mathbb{R}^{N\times d}$ with $\min_{i\in[N]}\|v_{i}\|\geq
1$ stored in the data structure from Lemma 6. Precision parameters
$\delta,\varepsilon$.

2:Vectors $c_{1},\dots,c_{k}\in\mathbb{R}^{d}$ corresponding to centroids.

3:Select $k$ initial centroids $c_{1}^{0},\dots,c_{k}^{0}$ and store them in a
QRAM data structure.

4:for $t=0$ until convergence do

5: Sample
$p=O\big{(}\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}\log\frac{k}{\delta}\big{)}$
rows $P=(s_{1},\dots,s_{p})$ and $R=(r_{1},\dots,r_{p})$ from $[N]$ and $[k]$,
respectively, by preparing the state
$\frac{1}{\|V\|_{F}}\sum_{i=1}^{N}\|v_{i}\||i\rangle|\ell_{i}^{t}\rangle$
(Theorem 9) and measuring both registers to obtain
$(s_{i},r_{i})\in\mathcal{C}_{r_{i}}^{t}\times[k]$. Let $P_{j}=\\{i\in
P|s_{i}=j\\}$.

6: For each $j\in[k]$, pick $Q_{j}\subseteq P_{j}$ at random such that
$|Q_{j}|/|P_{j}|=1/k$.

7: For $j\in[k]$, record the new centroids:
$\displaystyle\widetilde{c}_{j}^{t+1}=\frac{1}{N|Q_{j}|}\sum_{i\in
Q_{j}}\frac{\|V\|_{F}^{2}}{\|v_{r_{i}}\|^{2}}v_{r_{i}}.$

8:end for

######  Theorem 8.

Given a data matrix $V\in\mathbb{R}^{N\times d}$ with
$\min_{i\in[N]}\|v_{i}\|\geq 1$ stored in the data structure from Lemma 6 and
a parameter $\varepsilon>0$, Algorithm 2 outputs centroids consistent with the
$\varepsilon$-$k$-means with probability at least $1-\delta$ and in time per
iteration

$\displaystyle
O\left(\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}(\sqrt{k}d+\log(Nd))\log\frac{k\|V\|_{F}}{\delta
N}\right).$

###### Proof.

We start with the error analysis. Let $\chi_{j}^{t}\in\mathbb{R}^{N}$ be the
characteristic vector for cluster $j\in[k]$ at iteration $t$ scaled to
$\ell_{1}$-norm, i.e., $(\chi_{j}^{t})_{i}=\frac{1}{|\mathcal{C}_{j}^{t}|}$ if
$i\in\mathcal{C}^{t}_{j}$ and $0$ if $i\not\in\mathcal{C}^{t}_{j}$. Let
$X_{1},\dots,X_{q}\in\mathbb{R}^{N\times d}$ be i.i.d. random matrices formed
by sampling a row $j$ from $[k]$ with probability
$\frac{1}{\|V\|_{F}^{2}}\sum_{i\in\mathcal{C}_{j}^{t}}\|v_{i}\|^{2}$, followed
by sampling a row $i$ from $\mathcal{C}_{j}^{t}$ with probability
$\|v_{i}\|^{2}/\sum_{i\in\mathcal{C}_{j}^{t}}\|v_{i}\|^{2}$ (so a row $i$ is
then sampled from $[N]$ with probability $\|v_{i}\|^{2}/\|V\|_{F}^{2}$),
setting the $i$-th row of $X$ to $\|V\|_{F}^{2}\frac{v_{i}}{\|v_{i}\|^{2}}$
and the rest to zero. Define $\widehat{V}=\frac{1}{q}\sum_{j=1}^{q}X_{j}$
where $q:=p/k$. Then $\mathbb{E}[\widehat{V}]=V$. Finally, let
$\lambda_{j}=\frac{q}{|Q_{j}|}\frac{|\mathcal{C}_{j}^{t}|}{N}=\frac{p}{|P_{j}|}\frac{|\mathcal{C}_{j}^{t}|}{N}$.

We note that the standard $k$-means and Algorithm 2 can be stated,
respectively, as $c_{j}^{t+1}=V^{\top}\chi_{j}^{t}$ and
$\widetilde{c}_{j}^{t+1}=\lambda_{j}\widehat{V}^{\top}\chi_{j}^{t}$. The rest
of the error analysis is exactly the same as Theorem 7, i.e., the centroids
$\widetilde{c}_{j}^{t+1}$ are such that
$\|\widetilde{c}_{j}^{t+1}-c_{j}^{t+1}\|\leq\varepsilon$ with probability at
least $1-\delta/2$. The only other step prone to fail is producing the states
$\frac{1}{\|V\|_{F}}\sum_{i=1}^{N}\|v_{i}\||i\rangle|\ell_{i}^{t}\rangle$ with
success probability at least $1-\delta/2$. Finally, the runtime of each step
is:

  * •

Producing and then measuring the state
$\frac{1}{\|V\|_{F}}\sum_{i=1}^{N}\|v_{i}\||i\rangle|\ell_{i}^{t}\rangle$ a
number of $p$ times takes $O\big{(}p(\sqrt{k}d\log(p/\delta)+\log(Nd))\big{)}$
time according to Lemma 9;

  * •

Computing all $k$ new centroids by adding $q=p/k$ $d$-dimensional vectors
takes $O(pd)$ time.

This means the overall complexity is (already using that
$|\mathcal{C}_{j}^{t}|/N=\Omega(1/k)$)

$O\big{(}p(\sqrt{k}d\log(p/\delta)+\log(Nd))\big{)}=O\left(\frac{\|V\|_{F}^{2}}{N}\frac{k^{2}}{\varepsilon^{2}}(\sqrt{k}d+\log(Nd))\log\frac{k\|V\|_{F}}{\delta
N}\right).\qed$

######  Lemma 9.

Let a data matrix $V\in\mathbb{R}^{N\times d}$ and a centroid matrix
$C^{t}\in\mathbb{R}^{k\times d}$ be stored in the data structure from Lemma 6.
For any $\delta>0$, there is a quantum algorithm that prepares the state

$\displaystyle\frac{1}{\|V\|_{F}}\sum_{i=1}^{N}\|v_{i}\||i\rangle|\ell^{t}_{i}\rangle$

with probability at least $1-\delta/2$ in
$O(\sqrt{k}d\log(1/\delta)+\log(Nd))$ time.

###### Proof.

The first step is to use the data structure from Lemma 6 to create the state

$\displaystyle\frac{1}{\sqrt{k}\|V\|_{F}}\sum_{i=1}^{N}\sum_{j=1}^{k}\|v_{i}\||i,v_{i}\rangle|j,c_{j}^{t}\rangle|0\rangle$

in $O(\log(Nd))$ time, and then to compute the distance between each centroid
$c_{j}^{t}$ and vector $v_{i}$ in superposition via the mapping

$\displaystyle\frac{1}{\sqrt{k}\|V\|_{F}}\sum_{i=1}^{N}\sum_{j=1}^{k}|i,v_{i}\rangle|j,c_{j}^{t}\rangle|0\rangle\mapsto\frac{1}{\sqrt{k}\|V\|_{F}}\sum_{i=1}^{N}\sum_{j=1}^{k}|i,v_{i}\rangle|j,c_{j}^{t}\rangle|\|v_{i}-c_{j}^{t}\|\rangle$

by using a classical circuit to calculate $\|v_{i}-c_{j}^{t}\|$ in $O(d)$
time.

The second step is to assign a cluster
$\ell^{t}_{i}:=\arg\min_{j\in[k]}\|v_{i}-c_{j}^{t}\|$ to each vector $v_{i}$.
This is done by using the quantum minimum finding subroutine (Fact 3) in order
to map

$\displaystyle|v_{i}\rangle|0\rangle\mapsto|v_{i}\rangle|\ell^{t}_{i}\rangle$

with probability at least $1-\delta/2$ by having access to the mapping
$|v_{i},c_{j}^{t}\rangle|0\rangle\mapsto|v_{i},c_{j}^{t}\rangle|\|v_{i}-c_{j}^{t}\|\rangle$.
The time cost of finding the minimum is $O(\sqrt{k}T_{U}\log(1/\delta))$,
where $T_{U}=O(d\log{d})$ is the cost of the map
$|v_{i},c_{j}^{t}\rangle|0\rangle\mapsto|v_{i},c_{j}^{t}\rangle|\|v_{i}-c_{j}^{t}\|\rangle$.
By uncomputing the distances calculated in the first step, we arrive at the
desired state. ∎

## References

  * [ABG13] Esma Aïmeur, Gilles Brassard, and Sébastien Gambs.  Quantum speed-up for unsupervised learning.  Machine Learning, 90:261–287, 2013.
  * [AHKZ20] Jonathan Allcock, Chang-Yu Hsieh, Iordanis Kerenidis, and Shengyu Zhang.  Quantum algorithms for feedforward neural networks.  ACM Transactions on Quantum Computing, 1(1):1–24, 2020.
  * [ANFSW19] Sara Ahmadian, Ashkan Norouzi-Fard, Ola Svensson, and Justin Ward.  Better guarantees for $k$-means and euclidean $k$-median by primal-dual algorithms.  SIAM Journal on Computing, 49(4):FOCS17–97, 2019.
  * [AV07] David Arthur and Sergei Vassilvitskii.  $k$-means++: the advantages of careful seeding.  In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms, pages 1027–1035, 2007.
  * [BGJK20] Anup Bhattacharya, Dishant Goyal, Ragesh Jaiswal, and Amit Kumar.  On sampling based algorithms for $k$-means.  In Nitin Saxena and Sunil Simon, editors, 40th IARCS Annual Conference on Foundations of Software Technology and Theoretical Computer Science (FSTTCS 2020), volume 182 of Leibniz International Proceedings in Informatics (LIPIcs), pages 13:1–13:17, Dagstuhl, Germany, 2020. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.
  * [CGJ18] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation.  arXiv preprint arXiv:1804.01973, 2018.
  * [Das08] Sanjoy Dasgupta.  The hardness of $k$-means clustering.  Technical report, University of California, San Diego, 2008.
  * [DH96] Christoph Dürr and Peter Høyer.  A quantum algorithm for finding the minimum.  arXiv preprint quant-ph/9607014, 1996.
  * [GLM08a] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Architectures for a quantum random access memory.  Physical Review A, 78(5):052310, 2008.
  * [GLM08b] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum random access memory.  Physical review letters, 100(16):160501, 2008.
  * [GLT18] András Gilyén, Seth Lloyd, and Ewin Tang.  Quantum-inspired low-rank stochastic regression with logarithmic dependence on the dimension.  arXiv preprint arXiv:1811.04909, 2018.
  * [Gro96] Lov K. Grover.  A fast quantum mechanical algorithm for database search.  In Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, pages 212–219, 1996.
  * [Gro97] Lov K. Grover.  Quantum mechanics helps in searching for a needle in a haystack.  Physical review letters, 79(2):325, 1997.
  * [GST22] András Gilyén, Zhao Song, and Ewin Tang.  An improved quantum-inspired algorithm for linear regression.  Quantum, 6:754, 2022.
  * [HW79] John A. Hartigan and Manchek A. Wong.  Algorithm AS 136: A $k$-means clustering algorithm.  Journal of the royal statistical society. series c (applied statistics), 28(1):100–108, 1979.
  * [Jai23] Ragesh Jaiswal.  A quantum approximation scheme for $k$-means.  arXiv preprint arXiv:2308.08167, 2023.
  * [JKS14] Ragesh Jaiswal, Amit Kumar, and Sandeep Sen.  A simple ${D}^{2}$-sampling based PTAS for $k$-means and other clustering problems.  Algorithmica, 70(1):22–46, 2014.
  * [KLLP19] Iordanis Kerenidis, Jonas Landman, Alessandro Luongo, and Anupam Prakash.  q-means: A quantum algorithm for unsupervised machine learning.  Advances in Neural Information Processing Systems, 32, 2019.
  * [KM99] K. Krishna and M. Narasimha Murty.  Genetic $k$-means algorithm.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 29(3):433–439, 1999.
  * [KMN+02] Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman, and Angela Y. Wu.  A local search approximation algorithm for $k$-means clustering.  In Proceedings of the eighteenth annual symposium on Computational geometry, pages 10–18, 2002.
  * [KP16] Iordanis Kerenidis and Anupam Prakash.  Quantum recommendation systems.  arXiv preprint arXiv:1603.08675, 2016.
  * [KP20a] Iordanis Kerenidis and Anupam Prakash.  Quantum gradient descent for linear systems and least squares.  Physical Review A, 101(2):022316, 2020.
  * [KP20b] Iordanis Kerenidis and Anupam Prakash.  A quantum interior point method for LPs and SDPs.  ACM Transactions on Quantum Computing, 1(1):1–32, 2020.
  * [Llo82] Stuart Lloyd.  Least squares quantization in PCM.  IEEE transactions on information theory, 28(2):129–137, 1982.
  * [LMR13] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum algorithms for supervised and unsupervised machine learning.  arXiv preprint arXiv:1307.0411, 2013.
  * [LMR14] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum principal component analysis.  Nature Physics, 10(9):631–633, 2014.
  * [LVV03] Aristidis Likas, Nikos Vlassis, and Jakob J. Verbeek.  The global $k$-means clustering algorithm.  Pattern recognition, 36(2):451–461, 2003.
  * [MNV12] Meena Mahajan, Prajakta Nimbhorkar, and Kasturi Varadarajan.  The planar $k$-means problem is NP-hard.  Theoretical Computer Science, 442:13–21, 2012.
  * [Pra14] Anupam Prakash.  Quantum algorithms for linear algebra and machine learning.  University of California, Berkeley, 2014.
  * [Sho94] Peter W. Shor.  Algorithms for quantum computation: discrete logarithms and factoring.  In Proceedings 35th annual symposium on foundations of computer science, pages 124–134. Ieee, 1994.
  * [Sho99] Peter W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM review, 41(2):303–332, 1999.
  * [Tan19] Ewin Tang.  A quantum-inspired classical algorithm for recommendation systems.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 217–228, 2019.
  * [Tan21] Ewin Tang.  Quantum principal component analysis only achieves an exponential speedup because of its state preparation assumptions.  Physical Review Letters, 127(6):060503, 2021.
  * [Vat09] Andrea Vattani.  The hardness of $k$-means clustering in the plane.  Technical report, University of California, San Diego, 2009.