Quantum (Inspired) $D^2$-sampling with Applications

## Abstract

Abstract D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling is a fundamental component of sampling-based clustering algorithms such as k 𝑘 k italic_k -means++.
Given a dataset V ⊂ ℝ d 𝑉 superscript ℝ 𝑑 V\subset\mathbb{R}^{d} italic_V ⊂ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT with N 𝑁 N italic_N points and a center set C ⊂ ℝ d 𝐶 superscript ℝ 𝑑 C\subset\mathbb{R}^{d} italic_C ⊂ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT , D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling refers to picking a point from V 𝑉 V italic_V where the sampling probability of a point is proportional to its squared distance from the nearest center in C 𝐶 C italic_C .
Starting with C = { } 𝐶 C=\{\} italic_C = { } and iteratively D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling and updating C 𝐶 C italic_C in k 𝑘 k italic_k rounds is precisely k 𝑘 k italic_k -means++ seeding that runs in O ⁢ ( N ⁢ k ⁢ d ) 𝑂 𝑁 𝑘 𝑑 O(Nkd) italic_O ( italic_N italic_k italic_d ) time and gives O ⁢ ( log ⁡ ( k ) ) 𝑂 𝑘 O(\log{k}) italic_O ( roman_log ( start_ARG italic_k end_ARG ) ) -approximation in expectation for the k 𝑘 k italic_k -means problem.
We give a quantum algorithm for (approximate) D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling in the QRAM model that results in a quantum implementation of k 𝑘 k italic_k -means++ that runs in time O ~ ⁢ ( ζ 2 ⁢ k 2 ) ~ 𝑂 superscript 𝜁 2 superscript 𝑘 2 \tilde{O}(\zeta^{2}k^{2}) over~ start_ARG italic_O end_ARG ( italic_ζ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) .
Here ζ 𝜁 \zeta italic_ζ is the aspect ratio ( i.e., largest to smallest interpoint distance ) and O ~ ~ 𝑂 \tilde{O} over~ start_ARG italic_O end_ARG hides polylogarithmic factors in N , d , k 𝑁 𝑑 𝑘 N,d,k italic_N , italic_d , italic_k .
It can be shown through a robust approximation analysis of k 𝑘 k italic_k -means++ that the quantum version preserves its O ⁢ ( log ⁡ ( k ) ) 𝑂 𝑘 O(\log{k}) italic_O ( roman_log ( start_ARG italic_k end_ARG ) ) approximation guarantee.
Further, we show that our quantum algorithm for D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling can be dequantized using the sample-query access model of Tang [ 30 ] . This results in a fast quantum-inspired classical implementation of k 𝑘 k italic_k -means++, which we call QI- k 𝑘 k italic_k -means++ , with a running time O ⁢ ( N ⁢ d ) + O ~ ⁢ ( ζ 2 ⁢ k 2 ⁢ d ) 𝑂 𝑁 𝑑 ~ 𝑂 superscript 𝜁 2 superscript 𝑘 2 𝑑 O(Nd)+\tilde{O}(\zeta^{2}k^{2}d) italic_O ( italic_N italic_d ) + over~ start_ARG italic_O end_ARG ( italic_ζ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_k start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d ) , where the O ⁢ ( N ⁢ d ) 𝑂 𝑁 𝑑 O(Nd) italic_O ( italic_N italic_d ) term is for setting up the sample-query access data structure.
Experimental investigations show promising results for QI- k 𝑘 k italic_k -means++ on large datasets with bounded aspect ratio.
Finally, we use our quantum D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling with the known D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling-based classical approximation scheme ( i.e., ( 1 + ε ) 1 𝜀 (1+\varepsilon) ( 1 + italic_ε ) -approximation for any given ε > 0 𝜀 0 \varepsilon>0 italic_ε > 0 ) to obtain the first quantum approximation scheme for the k 𝑘 k italic_k -means problem with polylogarithmic running time dependence on N 𝑁 N italic_N .

## 1 Introduction

Data clustering and the $k$-means problem, in particular, have many applications in data processing.
The $k$-means problem is defined as: given a set of $N$ points $V=\{v_{1},...,v_{N}\}\subset\mathbb{R}^{d}$, and a positive integer $k$, find a set $C\subset\mathbb{R}^{d}$ of $k$ centers such that the cost function,

$$ $\Phi(V,C)\equiv\sum_{v\in V}\min_{c\in C}D^{2}(v,c),$ $$

is minimised. Here, $D(v,c)\equiv\norm{v-c}$ is the Euclidean distance between points $v$ and $c$.
Partitioning the points based on the closest center in the center set $C$ gives a natural clustering of the data points.
Due to its applications in data processing, a lot of work has been done in designing algorithms from theoretical and practical standpoints.
The $k$-means problem is known to be $\mathsf{NP}$-hard, so it is unlikely to have a polynomial time algorithm.
Much research has been done on designing polynomial time approximation algorithms for the $k$-means problem.
However, the algorithm used in practice to solve $k$-means instances is a heuristic, popularly known as the $k$-means algorithm (not to be confused with the $k$-means problem).
This heuristic, also known as Lloyd’s iterations [], iteratively improves the solution in several rounds.
The heuristic starts with an arbitrarily chosen set of $k$ centers. In every iteration, it (i) partitions the points based on the nearest center and (ii) updates the center set to the centroids of the $k$ partitions.
In the classical computational model, it is easy to see that every Lloyd iteration costs $O(Nkd)$ time.
This hill-climbing approach may get stuck at a local minimum or take a huge amount of time to converge, and hence, it does not give provable guarantees on the quality of the final solution or the running time.
In practice, Lloyd’s iterations are usually preceded by the $k$-means++ algorithm [], a fast sampling-based approach for picking the initial $k$ centers that also gives an approximation guarantee.
So, Lloyd’s iterations, preceded by the $k$-means++ algorithm, give the best of both worlds, theory and practice.
Hence, it is unsurprising that a lot of work has been done on these two algorithms.
The work ranges from efficiency improvements in specific settings to implementations in distributed and parallel models.
With the quantum computing revolution imminent, it is natural to talk about quantum versions of these algorithms and quantum algorithms for the $k$-means problem in general.

Early work on the $k$-means problem within the quantum setting involved efficiency gains from quantizing Lloyd’s iterations.
In particular, Aimeur, Brassard, and Gambs [] gave an $O(\frac{N^{3/2}}{\sqrt{k}})$ time algorithm for executing a single Lloyd’s iteration for the Metric $k$-median clustering problem that is similar to the $k$-means problem.
This was using the quantum minimum finding algorithm of Durr and Hoyer [].
Using quantum distance estimation techniques assuming quantum data access, Lloyd, Mohseni, and Rebentrost [] gave an $O(kN\log{d})$ time algorithm for the execution of a single Lloyd’s iteration for the $k$-means problem.
More recently, [] gave an approximate quantization of the $k$-means++ method and Lloyd’s iteration assuming QRAM data structure [] access to the data.
Interestingly, the running time has only a polylogarithmic dependence on the size $N$ of the dataset.
The algorithm uses quantum linear algebra procedures, and hence, there is dependence on certain parameters that appear in such procedures, such as the condition number $\kappa(V)$.
Since Lloyd’s iterations do not give an approximation guarantee, its quantum version is also a heuristic without a provable approximation guarantee.
A quantum implementation of $k$-means++ was also given in [].
Note that $k$-means++ gives an approximation guarantee of $O(\log{k})$ in expectation for the $k$-means problem.
However, this approximation guarantee does not immediately extend to the quantum implementation of $k$-means++ since the quantum procedure introduces $D^{2}$sampling errors, and one must show that the approximation analysis is robust against those errors. This was indeed missing from [].
$D^{2}$-sampling is a fundamental component of $k$-means++.
Given the dataset $V$ and a center set $C$, $D^{2}$-sampling refers to picking a point from $V$ such that the sampling probability is proportional to its squared distance from the nearest center in $C$. Starting with $C=\{\}$ and $D^{2}$-sampling and updating $C$ in $k$ rounds is precisely the $k$-means++ algorithm.
In a recent line of work [, ], a robust version of $k$-means++ has been analyzed where the $D^{2}$-sampling probability has a multiplicative $(1\pm\varepsilon)$ error. This is called noisy $k$-means++, and it is now known [] that this algorithm gives $O(\log{k})$ approximation. This shows that $k$-means++ is indeed robust against small relative $D^{2}$-sampling errors.
This further means that the quantum implementation of $k$-means++ in [] gives an approximation guarantee of $O(\log{k})$.

In this work, we design quantum and quantum-inspired classical algorithms for $D^{2}$-sampling-based clustering algorithms.
We start with $k$-means++.
Quantizing(^1^11Quantizing is a term that is used to refer to designing a quantum version of an algorithm with a significant efficiency advantage over the corresponding classical algorithm.) $k$-means++ reduces to quantizing the $D^{2}$-sampling step.
Most of the ideas for quantizing $D^{2}$-sampling such as coherent amplitude estimation, median estimation, distance estimation, minimum finding, etc., have already been developed in previous works on quantum unsupervised learning (see [] for examples of several such tools). In particular, [] used these tools to give a quantum implementation of $D^{2}$-sampling. We provide a similar quantum algorithm and then combine this with the known $O(\log{k})$ approximation guarantee for noisy $k$-means++ to make the following observation.

###### Theorem 1 .

There is a quantum implementation of $k$-means++ that runs in time $\tilde{O}(\zeta^{2}k^{2})$ and gives an $O(\log{k})$ factor approximate solution for the $k$-means problem with a probability of at least $0.99$. Here, $\tilde{O}$ hides $\log^{2}{(Nd)}$ and $\log^{2}{(kd)}$ terms.(^2^22The output of $k$-means++ is a subset of data points, which can be stated as a subset of indices. If the description of centers is the expected output, the running time will also include a factor of $d$.)

Dependence on aspect ratio $\zeta$: Note that the running time depends on the aspect ratio $\zeta$, which is defined to be the ratio of the maximum to the minimum interpoint distance. For the remaining discussion, it is important to understand that the $k$-means problem remains $\mathsf{NP}$-hard to approximate (beyond a fixed constant) even when restricted to instances where $\zeta$ upper bounded by a small constant. See [] for such an hardness reduction.(^3^33The hardness of approximation argument in [] goes through the vertex cover problem. For a given graph $G=(V,E)$, the $k$-means instance includes $|E|$ points in a $|V|$-dimensional space. The point corresponding to an edge $(i,j)$ has $1$ in coordinates $i,j$ and 0 0 everywhere else. It must be clear that $\zeta=\sqrt{2}$ in this case.)

The above result is another addition to the growing area of Quantum Machine Learning (QML) where the goal is to design quantum algorithms for machine learning problems with large speed-ups over their classical counterparts.
There was a flurry of activity in QML where significant speedups were shown for problems such as recommendation systems [].
However, it was not clear whether such speedups could entirely be attributed to the quantum data access model (QRAM) upon which these works were built. The works of Tang [, , ] gave much clarity on this aspect. These works showed that similar speedups (up to some polynomial factors) are achievable in the classical model using a classical counterpart of QRAM, now popularly known as the sample-query (SQ) access data structure.
In other words, most QML algorithms could be dequantized.
Interestingly, by setting up the sample-query access data structure, the dequantized algorithms give rise to quantum-inspired classical algorithms.
The interesting property of such a quantum-inspired algorithm is that apart from the linear time spent on setting up the sample-query access data structure, which can be thought of as preprocessing time, the algorithm is extremely fast (otherwise, it would not have given the quantum speedup).
Our second contribution is to dequantize $D^{2}$-sampling using Tang’s sample-query access model.
This turns out to be simple since the sample-query access model primarily involves sampling vectors based on squared $\ell_{2}$ norms, which is how $D^{2}$-sampling also works. Based on this, we give a quantum-inspired classical implementation of $k$-means++.
The following theorem formally states our main result on this.

###### Theorem 2 .

There is a classical implementation of $k$-means++ (which we call QI-$k$-means++) that runs in time $O(Nd)+\tilde{O}(\zeta^{2}k^{2}d)$. Moreover, with probability at least $0.99$, QI-$k$-means++ outputs $k$ centers that give $O(\log{k})$-approximation in expectation for the $k$-means problem.

Note that the $O(Nd)$ term in the running time is for setting up the sample-query-access data structure, which can be thought of as a preprocessing operation.
This setting of the data structure can easily be parallelized on a multi-core processor.
Once the data structure is available, the running time is sublinear in $N$. The linear dependence on the dimension $d$ can also be changed to sublinear at the cost of worsening the dependence on the aspect ratio $\zeta$. This is through implementing a small relative error $D^{2}$-sampling procedure, which gives an implementation of the noisy $k$-means++ and which we know also gives $O(\log{k})$-approximation.
We state this result formally below.

###### Theorem 3 .

There is a classical implementation of noisy $k$-means++ (which we call QI-noisy-$k$-means++) that runs in time $O(Nd)+\tilde{O}(\zeta^{6}k^{2})$. Moreover, with probability at least $0.99$, QI-noisy-$k$-means++ outputs $k$ centers that give $O(\log{k})$-approximation in expectation for the $k$-means problem.

There have been several fast classical implementations of $k$-means++.
These implementations are fast under certain assumptions on the data set or achieved through trading efficiency with an approximation guarantee or both. Bachem et al.  [] gave a Monte Carlo Markov Chain (MCMC) based implementation of $k$-means++ with a running time $O(k^{3}d\log^{2}{n}\log{k})$. However, this algorithm, called $\mathsf{K}$-$\mathsf{MC}^{2}$, preserves the $O(\log{k})$ approximation guarantee of $k$-means++ only when the dataset satisfies some strong properties that are $\mathsf{NP}$-hard to check. The MCMC-based algorithm was improved by the same authors in [] who gave an algorithm with running time $O(Nd)+O(k^{2}d\log{k})$. The approximation guarantee of this algorithm has an additive term (in addition to the multiplicative $O(\log{k})$ factor), which could possibly be large. They also give another algorithm without the additive factor and with running time $O(Nd)+O(k^{3}d\log{k})$, but that works under the assumption made in []. More recently, a fast implementation based on the multi-tree embedding of the data was given by Cohen-Addad et al.  []. Their implementation runs in time $\tilde{O}(Nd+kN^{\varepsilon}\log^{1+\varepsilon}{\zeta})$ and gives an approximation guarantee of $O\left(\frac{\log{k}}{\varepsilon^{3}}\right)$ for any $\varepsilon>0$.(^4^44The second term in the running time expression is independent of dimension $d$ since we can apply the Johnson Lindenstrauss theorem to replace $d$ with $\log{N}$.)
Here $\zeta$ is the aspect ratio, $\tilde{O}$ notation hides logarithmic dependences on $N$, and the $O(Nd)$ term can be thought of as the preprocessing step where the multi-tree embedding data structure is set up. Let us compare our sample-query access data structure based on fast implementation with these known implementations.
In the assumption-free case (i.e., no assumptions on the data set), if we leave aside the preprocessing costs, which is the same $O(Nd)$, our implementation has a better dependence on $N$ ($\log^{2}{N}$ versus $N^{\varepsilon}$) but worse dependence on the aspect ratio ($\zeta^{2}$ versus $\log^{1+\varepsilon}{\zeta}$) and $k$ ($k^{2}$ versus $k$).
So, in cases where the aspect ratio is not large (e.g., binarized data such as MNIST []) and $k$ is reasonably small, our implementation should be expected to run fast and give good solutions. Note that unlike the result of [], our results have no tradeoff between efficiency and approximation ratio.
Another advantage of the sample-and-query data structure is that it is very efficient to update, i.e., add or delete a point that costs $O(\log{N})$ time. This is useful in scenarios where new data points get added periodically, and one needs to recluster the data.
Even though our work is mainly theoretical in nature, we conduct some experimental investigations to see the running time advantage of QI-$k$-means++ over $k$-means++ (the classical implementation).

Finally, we design a quantum approximation scheme for the $k$-means problem with polylogarithmic dependence on the size $N$ of the data set.
We do this by quantizing the highly parallel, $D^{2}$-sampling-based approximation scheme of [].
An approximation scheme is an algorithm that, in addition to the dataset and $k$, takes an error parameter $\varepsilon>0$ as input and outputs a solution with a cost within $(1+\varepsilon)$ factor of the optimal.
The $k$-means problem is $\mathsf{NP}$-hard to efficiently approximate beyond an approximation factor of $1.07$ [].
The tradeoff in obtaining this fine-grained approximation of $(1+\varepsilon)$ is that the running time of our algorithm has an exponential dependence on $k$ and error parameter $\varepsilon$.
In the classical setting, such algorithms are categorized as Fixed Parameter Approximation Schemes (fpt-AS).
Such $(1+\varepsilon)$-approximation algorithms can have exponential running time dependence on the parameter (e.g., the number of clusters $k$ in our setting).
The practical motivation for studying Fixed-Parameter Tractability (FPT) for computationally hard problems is that when the parameter is small (e.g., number of clusters $k\sim 5$), the running time is not prohibitively large.
We state our main result as the following theorem, which we will prove in the remainder of the paper.

###### Theorem 4 .

Let $0<\varepsilon<1/2$ be the error parameter.
There is a quantum algorithm that, when given QRAM data structure access to a dataset $V\in\mathbb{R}^{N\times d}$, runs in time $\tilde{O}\left(2^{\tilde{O}(\frac{k}{\varepsilon})}d\zeta^{2}\right)$ and outputs a $k$ center set $C\in\mathbb{R}^{k\times d}$ such that with high probability $\Phi(V,C)\leq(1+\varepsilon)\cdot OPT$. Here, $\zeta$ is the aspect ratio, i.e., the ratio of the maximum to the minimum distance between two given points in $V$.(^5^55The $\tilde{O}$ notation hides logarithmic factors in $N$. The $\tilde{O}$ in the exponent hides logarithmic factors in $k$ and $1/\varepsilon$.)

We end on the note that the above quantum approximation scheme can also be dequantized in the sample-query-access model of Tang [].

###### Theorem 1 .

###### Theorem 2 .

###### Theorem 3 .

###### Theorem 4 .

### 1.1 Related work

We have already discussed past research works on quantum versions of the $k$-means algorithm (i.e., Lloyd’s iterations).
This includes [], [], and [].
All these have been built using various quantum tools and techniques developed for various problems in quantum unsupervised learning, such as coherent amplitude and median estimation, distance estimation, minimum finding, etc. See [] for examples of several such tools.
We have also discussed previous works on fast implementations of the $k$-means++ algorithm.
Other directions on quantum $k$-means includes adiabatic algorithms (e.g., []) and algorithms using the QAOA framework (e.g., [, ]). However, these are without provable guarantees.
A recent work of Doriguello et al. [] improves the results of [] on Lloyd’s iterations and obtains a quantum algorithm with better running time dependence and removes the dependence on the data matrix dependent parameters such as the condition number.
They also dequantize their algorithm in the sample-query access model of Tang [].
However, since Lloyd’s iterations do not give any approximation guarantee, their algorithms also do not have a provable approximation guarantee.
In another recent work [], the quantum algorithm of [] is used in a supervised setting to construct decision trees.

## 2 Quantum (Inspired) D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -Sampling

In this section, we discuss the key ideas in the design of our quantum algorithm for $D^{2}$-sampling and its dequantization in the sample-query (SQ) access model of Tang []. Let us first look at the main ideas of our quantum algorithm in the following subsection.

### 2.1 Quantum D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling

We will work under the assumption that the minimum distance between two data points is $1$, which can be achieved using scaling. This makes the aspect ratio $\zeta$ simply the maximum distance between two data points.
We will use $i$ for an index into the rows of the data matrix $V\in\mathbb{R}^{N\times d}$, and $j$ for an index into the rows of the center matrix $C\in\mathbb{R}^{m\times d}$.
We would ideally like to design a quantum algorithm that performs the transformation:

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{{D(v_{i},c_{j})}}$ (1) $$

The known quantum tools such as swap test followed by coherent amplitude estimation, and median estimation allow one to obtain a state that is an approximation of the ideal state that is shown on the right in (1). Moreover, this can be done in time that has only polylogarithmic dependence on the input size.
The details are provided in the Appendix. For the current high-level discussion, we will assume that the ideal state can be prepared.
We would like to perform sampling from the $D^{2}$ distribution over $V$ with respect to the center set $C$. This means that we would like to sample an index $i\in[N]$ with probability:

$$ $\mathbf{Pr}[i]=\frac{\min_{j\in[m]}{D(v_{i},c_{j})^{2}}}{\sum_{i^{\prime}\in[N ]}\min_{j\in[m]}{D(v_{i^{\prime}},c_{j})^{2}}}$ (2) $$

So, if we can prepare the state:

$$ $\frac{1}{\sqrt{N}}\sum_{i\in[N]}\ket{i}\ket{\min_{j\in[m]}D(v_{i},c_{j})},$ $$

then we can obtain a good estimate of the denominator of the right expression in (2), which is basically the $k$-means cost with respect to the current center set $C$. This follows from the fact that measuring the second register of the above state gives the distance of a random data point to its closest center in $C$. Repeatedly preparing the state and measuring helps obtain a close estimate of the $k$-means cost. Now that the $k$-means cost (i.e., the denominator of (2)) is known, using a few more standard quantum transformations (controlled rotations), we can pull out the state of the second register as part of the amplitude of an ancilla qubit. In other words, an approximate version of the following quantum state can be prepared:

$$ $\frac{1}{\sqrt{N}}\sum_{i\in[N]}\ket{i}\left(\beta_{i}\ket{0}+\sqrt{1-\beta_{i }}\ket{1}\right),$ $$

where $\beta_{i}=\frac{1}{Z}\min_{j\in[m]}D(v_{i},c_{j})$ and $Z$ is some appropriate normalization term. So, by measuring the above state and doing rejection sampling (ignoring the measurement when the ancilla qubit is $1$), we obtain a $D^{2}$-sample.
Some of the quantum states shown above are ’ideal’ quantum states that we have used to convey the main ideas of the quantum algorithm. The real states will be close approximations of these states, and we need to carefully account for all the errors introduced. All the details are given in the Appendix. Finally, note that this also gives a quantum implementation of the $k$-means++ algorithm since the algorithm simply performs $D^{2}$-sampling in $k$ iterations while updating the center set.

In the following subsection, we discuss how to dequantize this quantum algorithm in the sample-query access model of Tang []. This gives a quantum-inspired classical algorithm for $D^{2}$-sampling.

### 2.2 Quantum inspired D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling

Ewin Tang’s thesis [] nicely categorizes quantum machine learning (QML) algorithms into ones where the quantum advantage (i.e., polylogarithmic running time) is solely because of the quantum access to the data and the ones where it is not (e.g., the HHL algorithm []).
Quantum access means that a (normalized) data vector $v\coloneqq(v(1),...,v(d))^{T}$ can be loaded onto the quantum workspace in $O(1)$ time as $\ket{v}=v(1)\ket{1}+...+v(d)\ket{d}$.
One implication of having the quantum state $\ket{v}$ is that $i$ can be sampled with probability $\frac{v(i)^{2}}{\sum_{j}v(j)^{2}}$ by simply making a measurement in the standard basis.
One of the key insights in Tang’s Thesis [] is that QML algorithms that benefit solely from the quantum data access can be ‘dequantized’ if we work with an appropriate classical counterpart of the quantum data access that allows sampling access of the kind described above (in addition to the classical access to the data).
This is called sample-query access, or SQ access in short.
Dequantization means that a classical algorithm can simulate the steps of the quantum algorithm using SQ access with similar running time dependence.
The nice property of the SQ access data structures is that it can be constructed classically in linear time.
Note that linear time means that we lose the quantum advantage.
However, if we keep this aside and count the SQ access data structure construction as preprocessing, the remaining computation has a similar advantage as that of the quantum algorithm.
This is a level playing field with the QML algorithm that allows fair comparison, specifically given that setting up quantum access (called QRAM) also takes linear time.
Such classical algorithms obtained by dequantization are called quantum inspired algorithms.
Our quantum algorithms based on $D^{2}$-sampling fall into the category of QML algorithms that can be dequantized in Tang’s SQ access model.
This section discusses the quantum-inspired classical algorithm we obtain by dequantizing our quantum algorithm using the SQ access data structure.
Naturally, we start the discussion by looking at the definition of the SQ access data structure.

###### Definition 1 (Query access, Definition 1.1 in [ 30 ] ) .

For a vector $\vec{v}\in\mathbb{C}^{n}$, we have query access to $\vec{v}$, denoted by $Q(\vec{v})$, if for all $i\in[n]$, we can query for $\vec{v}(i)$. We use $\mathbf{q}(\vec{v})$ to denote the time (cost) of such a query.

###### Definition 2 (SQ-access to a vector, Definition 1.2 in [ 30 ] ) .

For a vector $\vec{v}\in\mathbb{C}^{n}$, we have sampling and query access to $\vec{v}$, denoted by $SQ(\vec{v})$, if we can:

- 1.
query for entries in $\vec{v}$ as in $Q(\vec{v})$. The time cost is denoted by $\mathbf{q}(\vec{v})$.
- 2.
obtain independent samples $i\in[n]$ where the probability of sampling $i$ is $\frac{\vec{v}(i)^{2}}{\norm{\vec{v}}^{2}}$. This distribution is denoted by $\mathcal{D}_{\vec{v}}$.
The time cost is denoted by $\mathbf{s}(\vec{v})$.
- 3.
query for $\norm{\vec{v}}$. The time cost is denoted by $\mathbf{n}(\vec{v})$.

Let $\mathbf{sq}(\vec{v})\equiv\max{\{\mathbf{q}(\vec{v}),\mathbf{s}(\vec{v}),
\mathbf{n}(\vec{v})\}}$ denote the time cost of SQ-access.

Figure: Figure 1: A tree data structure to enable sample-query access to an example vector of dimension $n=4$. Index $i$ can be sampled with probability $\frac{|\vec{v}_{i}|^{2}}{\sum_{j}|\vec{v_{j}|^{2}}}$ in $O(\log{n})$ time by traversing down the tree.
Refer to caption: extracted/2405.13351v1/tree.png

A simple tree-based data structure (see Figure 1) supports sample-query access for vectors and matrices. The details can be found in []. Here, we give a summary of the running time of various operations.

###### Lemma 1 (Remark 4.12 in [ 30 ] ) .

There is a data structure for storing a vector $\vec{v}\in\mathbb{R}^{n}$ supporting the following operations:

- 1.
Reading and updating an entry of $\vec{v}$ in $O(\log{n})$ time. So, $\mathbf{q}(\vec{v})=O(\log{n})$.
- 2.
Finding $\norm{\vec{v}}^{2}$ in $O(1)$ time. So, $\mathbf{n}(\vec{v})=O(1)$.
- 3.
Sampling from $\mathcal{D}_{\vec{v}}$ in $O(\log{n})$ time. So, $\mathbf{s}(\vec{v})=O(\log{n})$

It follows from the above that $\mathbf{sq}(\vec{v})=O(\log{n})$.
Moreover, the time required to construct the data structure is $O(n)$.

Let us outline the key ideas in the quantum-inspired classical algorithm in the SQ access model defined above.
Let $v_{1},...,v_{N}$ denote the data vectors, i.e., $v_{i}=V(i,.)$.
Suppose we have set up SQ access for the vectors $v_{1},...,v_{N}$, $(\norm{v_{1}}^{2},...,\norm{v_{N}}^{2})^{T}$ (together denoted by $SQ(V)$), and $m$ centers $c_{1},...,c_{m}$ (denoted by $SQ(c_{1}),...,SQ(c_{m})$).
Our goal is to set up SQ access for the vector $w$ defined as:

$$ $w\coloneqq\left(\min_{j\in[m]}D(v_{1},c_{j}),...,\min_{j\in[m]}D(v_{N},c_{j}) \right)^{T}$ (3) $$

The reason is that once we have SQ access, $SQ(w)$, a sampling query returns an index with distribution $\mathcal{D}_{w}$, which is precisely the $D^{2}$-sampling distribution.
We enable SQ access for $w$ in a sequence of steps:

$$ $(1)\ SQ(V),SQ(c_{1}),...,SQ(c_{m})\quad(2)\ SQ(u_{1}),...,SQ(u_{m})\quad(3)SQ( w).$ $$

Here $u_{j}$ is the vector defined as:

$$ $\forall i\in[N],u_{j}(i)=D(v_{i},c_{j})=\norm{v_{i}-c_{j}},$ (4) $$

the vector of distances of $N$ data points from the $j^{th}$ center.
The above steps are a gross simplification of the real steps to understand the high-level ideas and draw a parallel to the corresponding quantum steps.
In actual implementation, we may not be able to enable sample-query access for $w$ but something known as oversampling and query access, which is a generalization SQ access.
Much of the technical effort is spent designing these oversampling query accesses. These details are given in the Appendix.

Let us draw a parallel to the quantum algorithm to see that the above SQ-based algorithm is indeed quantum-inspired. The relevant quantum states involved in the $D^{2}$-sampling algorithm are given below:

$$ $\displaystyle(1)\ \forall i,\ket{v_{i}},\forall j,\ket{c_{j}},\frac{1}{\sqrt{N }}\sum_{i}\ket{i}\ket{v_{i}}\quad(2)\ \frac{1}{\sqrt{N}}\sum_{i}\ket{i}\ket{D( v_{i},c_{1})}...\ket{D(v_{i},c_{m})}$ $\displaystyle(3)\ \frac{1}{\sqrt{N}}\sum_{i}\ket{i}\ket{\min_{j\in[m]}D(v_{i}, c_{j})}\quad(4)\ \frac{1}{\sqrt{N}}\sum_{i\in[N]}\ket{i}\left(\beta_{i}\ket{0} +\sqrt{1-\beta_{i}}\ket{1}\right)$ $$

The correspondence between quantum states and SQ accesses is apparent. So, every step in the quantum algorithm can be simulated using oversampling query access, which leads to a quantum-inspired classical algorithm since the SQ access framework is completely classical.
Let us now discuss some key elements of the quantum-inspired classical implementation of $k$-means++ that results from the above quantum-inspired classical $D^{2}$-sampling.

###### Definition 1 (Query access, Definition 1.1 in [ 30 ] ) .

###### Definition 2 (SQ-access to a vector, Definition 1.2 in [ 30 ] ) .

###### Lemma 1 (Remark 4.12 in [ 30 ] ) .

### 2.3 QI- k 𝑘 k italic_k -means++

If we use the quantum-inspired classical $D^{2}$-sampling iteratively in $k$ rounds to pick $k$ centers, the resulting algorithm is a quantum-inspired classical implementation of $k$-means++, which we call the QI-$k$-means++ algorithm. This can be seen as a fast implementation of $k$-means++. We give a high-level outline of the algorithm and discuss some of its important properties.

Figure: Algorithm 1 QI-$k$-means++($V,k$)

Setting up the sample-query access for the dataset $V\in\mathbb{R}^{N\times d}$ costs $O(Nd)$ time. This can be thought of as the preprocessing time. Setting up the sample-query access is highly parallelizable, as can be seen from the tree-based data structure in Figure 1.
On a shared memory, multi-processor system with $M$ processing units, the task can be done in $O(Nd/M)$ time. After the preprocessing in step (1), the cost of the remaining $D^{2}$-sampling steps is $\tilde{O}(\zeta^{2}k^{2}d)$.
Here, again, it is possible to parallelize on a multi-processor system since much of the time in $D^{2}$-sampling is spent on waiting for rejection sampling to succeed. The rejection sampling can be performed independently on multiple processors until one of the processing units succeeds. This cuts down the time to $\tilde{O}(\zeta^{2}k^{2}d/M)$.
We implement the QI-$k$-means++ algorithm to compare the $k$-means cost and time with the classical implementation of $k$-means++ that has a running time $O(Nkd)$.

## 3 A Quantum Approximation Scheme

We convert the $D^{2}$-sampling-based approximation scheme of [] to a Quantum version. The approximation scheme is simple and highly parallel, which can be described in a few lines.

Figure: Algorithm 2 Approximation Scheme

We will discuss the quantization of the above three steps of the approximation scheme of [], thus obtaining a quantum approximation scheme. (^6^66Steps (2) and (3) in the algorithm are within a loop for probability amplification. For simplicity, this loop is skipped in this high-level description.)

(Constant approximation)
The first step requires finding a constant factor approximate solution for the $k$-means problem.
Even though several constant factor approximation algorithms are known, we need one with a quantum counterpart that runs in time that is polylogarithmic in the input size $N$.
One such algorithm is the $k$-means++ seeding algorithm [].
[] give an approximate quantum version of $D^{2}$-sampling.
The approximation guarantee of the $k$-means++ algorithm is $O(\log{k})$ instead of the constant approximation required in the approximation scheme of [].
It is known from the work of [] that if the $D^{2}$-sampling in $k$-means++ is continued for $2k$ steps instead of stopping after sampling $k$ centers, then we obtain a center set of size $2k$ that is a $(2,O(1))$-pseudo approximate solution.
This means that this $2k$-size center set has a $k$-means cost that is some constant times the optimal.
Such a pseudo-approximate solution is sufficient for the approximation scheme of [] to work.
We show that the pseudo-approximation guarantee of [] also holds when using the approximate quantum version of the $D^{2}$-sampling procedure.

($D^{2}$-sampling)
The second step of [] involves $D^{2}$-sampling, which we already discussed how to quantize.
This is no different than the $D^{2}$-sampling involved in the $k$-means++ algorithm of the previous step.
The sampling in this step is simpler since the center set $C$, with respect to which the $D^{2}$-sampling is performed, does not change (as is the case with the $k$-means++ algorithm.)

(All subsets)
Since the number of points sampled in the previous step is $poly(\frac{k}{\varepsilon})$, we need to consider a list of $\left(\frac{k}{\varepsilon}\right)^{\tilde{O}(\frac{k}{\varepsilon})}$ tuples of subsets, each giving a $k$-center set (a tuple $(S_{1},...,S_{k})$ defines $(\mu(S_{1}),...,\mu(S_{k}))$).
We need to compute the $k$-means cost for every $k$ center set in the list and then pick the one with the least cost.
We give quantization of the above steps.
(^7^77Note that when picking the center set with the least cost, we can get quadratic improvement in the search for the best $k$-center set using quantum search. Given that the search space is of size $\left(\frac{k}{\varepsilon}\right)^{\tilde{O}(\frac{k}{\varepsilon})}$, this results only in a constant factor improvement in the exponent. So, for simplicity, we leave out the quantum search from the discussion.)

Note that the quantization of the classical steps of [] will incur precision errors.
So, we first need to ensure that the approximation guarantee of [] is robust against small errors in distance estimates, $D^{2}$-sampling probabilities, and $k$-means cost estimates.
We must carefully account for errors and ensure that the quantum algorithm retains the $(1+\varepsilon)$ approximation guarantee of the robust version of []. We give the detailed analysis in the Appendix.

## 4 Experiments

We provide an experimental investigation of the QI-$k$-means++ algorithm. We implemented our code in C++ for both $k$-means++ and QI-$k$-means++ and the results provided here are averaged over $5$ runs. No pre-processing/dimensionality reductions were done on any of the datasets used. The experiments were performed on an AMD Ryzen 9 5900HX 4.68 Ghz 8 Core processor (parallelization on multicore was not used). We must set up the sample and query data structure for the dataset only once to be able to perform seeding for any value of $k$, with each run taking time only polylogarithmic in $N$. Specifically, if we wanted to perform seeding for $q$ different values of $k$, say $\{k_{1},...,k_{q}\}$, then the total runtime for QI-$k$-means++ would be $O(Nd)+\tilde{O}(\zeta^{2}d\sum_{j\in[q]}k_{j}^{2})$ as compared to $O(Nd\sum_{j\in[q]}k_{j})$ for $k$-means++. This is one reason to consider the ${O}(Nd)$ term as pre-processing.

We study the runtime (including setup time for the sample and query data structure) of QI-k-means++ on two extreme types of datasets to demonstrate its possible use cases. The first is binarized MNIST [] (70,000 data points, each being a $28\times 28$ pixel image . Each pixel has value either 0 or 1, as opposed to the original MNIST, which had values from 0 to 255) and the second is IRIS [] (150 data points with 4 feature values). The graphs below show the cumulative runtime (i.e., the total time for setting up and calculating cluster centers for all values of $k$ uptill a certain value).

**Table 1: Clustering cost for binarized MNIST (costs are scaled down by a factor of $10^{6}$)**
| k | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| QI-k-means++ | 8.43 | 8.13 | 7.81 | 7.43 | 7.29 | 7.09 | 6.79 | 6.78 | 6.66 |
| k-means++ | 8.35 | 8.10 | 7.74 | 7.38 | 7.44 | 7.17 | 7.12 | 6.91 | 6.83 |

**Table 2: Clustering cost for IRIS (costs are rounded to 2 decimal places)**
| k | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| QI-k-means++ | 538.82 | 253.41 | 112.53 | 96.04 | 83.72 | 64.92 | 58.62 | 56.25 | 49.92 |
| k-means++ | 773.26 | 227.50 | 115.33 | 93.45 | 83.21 | 60.14 | 57.98 | 53.99 | 50.35 |

Advantageous regime: For a dataset with a large number of points and small aspect ratio (for example, binary MNIST), we find that for one iteration of seeding, our algorithm performs slower but quickly catches up to be significantly faster when we require seeding to be done for multiple values of $k$ (which is generally the case in practice since in the unsupervised setting, the optimal value of $k$ is not known beforehand). This is because most of the time is spent in setting up the data structure, after which the sampling becomes significantly faster. For example, see left plot for the MNIST dataset. Notice that the cumulative runtime for QI-$k$-means++ is almost constant, while for $k$-means++, it does not scale well due to the multiplicative dependence on $N$.

Disadvantageous regime: For a dataset with a small number of points and a large aspect ratio (for example, IRIS), we find that the runtime is dominated by the sampling term and due to the quadratic dependence on $k$ and $\zeta$, QI-$k$-means++ may not scale well. For example, see right plot for the IRIS dataset.

## Appendix A Quantum D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling (proof of Theorem 1 )

We will work under the assumption that the minimum distance between two data points is $1$, which can be achieved using scaling. This makes the aspect ratio $\zeta$ simply the maximum distance between two data points.
We will use $i$ for an index into the rows of the data matrix $V\in\mathbb{R}^{N\times d}$, and $j$ for an index into the rows of the center matrix $C\in\mathbb{R}^{k\times d}$.
We would ideally like to design a quantum algorithm that performs the transformation:

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{{D(v_{i},c_{j})}}$ $$

Let us call the state on the right $\ket{\Psi_{ideal}}$.
This is an ideal quantum state for us since $\ket{\Psi_{ideal}}$ helps to perform $D^{2}$-sampling and to find the $k$-means cost of clustering, which are the main components of the approximation scheme of [] that we intend to use.
One caveat is that we will only be able to perform the following transformation (instead of the abovementioned transformation)

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{\psi_{i,j}},$ $$

where $\ket{\psi_{i,j}}$ is an approximation for $\ket{\mathcal{D}(v_{i},c_{j})}$ in a sense that we will make precise below.
We will use $\ket{\Psi_{real}}$ to denote the state $\ket{i}\ket{j}\ket{\psi_{i,j}}$.
This state is prepared using tools such as swap test followed by coherent amplitude estimation, and median estimation.
Since these tools and techniques are known from previous works [, , ], we summarise the discussion (see Section 4.1 and 4.2 in []) in the following lemma.

###### Lemma 2 ( [ 20 ] and [ 31 ] ) .

Assume for a data matrix $V\in\mathbb{R}^{N\times d}$ and a center set matrix $C\in\mathbb{R}^{t\times d}$ that the following unitaries: (i) $\ket{i}\ket{0}\rightarrow\ket{i}\ket{v_{i}}$, (ii) $\ket{j}\ket{0}\rightarrow\ket{j}\ket{c_{j}}$ can be performed in time $T$ and the norms of the vectors are known. For any $\Delta>0$, there is a quantum algorithm that in time $O\left(\frac{T\log{\frac{1}{\Delta}}}{\varepsilon}\right)$ computes:

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{\psi_{i,j}},$ $$

where $\ket{\psi_{i,j}}$ satisfies the following two conditions for every $i\in[N]$ and $j\in[t]$:

- (i)
$\norm{\ket{\psi_{i,j}}-\ket{0^{\otimes\ell}}\ket{{\tilde{D}(v_{i},c_{j})}}}
\leq\sqrt{2\Delta}$, and
- (ii)
For every $i,j$, $\tilde{D}(v_{i},c_{j})\in(1\pm\varepsilon)\cdot D(v_{i},c_{j})$.

In the subsequent discussions, we will use $T$ as the time to access the QRAM data structure [], i.e., for the transitions $\ket{i}\ket{0}\rightarrow\ket{i}\ket{v_{i}}$ and $\ket{j}\ket{0}\rightarrow\ket{j}\ket{c_{j}}$ as given in the above lemma.
This is known to be $T=O(\log^{2}{(Nd)})$. Moreover, the time to update each entry in this data structure is also $T=O(\log^{2}{(Nd)})$.
This is the logarithmic factor that is hidden in the $\tilde{O}$ notation.
In the following subsections, we discuss the utilities of $\ket{\Psi_{real}}$ for the various components of the approximation scheme of [].
During these discussions, it will be easier to see the utility first with the ideal state $\ket{\Psi_{ideal}}$ before the real state $\ket{\Psi_{real}}$ that can actually be prepared.
We will see how $\ket{\Psi_{real}}$ is sufficient within a reasonable error bound.

###### Lemma 2 ( [ 20 ] and [ 31 ] ) .

### A.1 Finding distance to closest center

Let us see how we can estimate the distance of any point to its closest center in a center set $C$ with $t\leq k$ centers.
We can use the transformation $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{{D(v_{i},c_{j})}}$ to prepare the following state for any $i$:

$$ $\ket{i}\ket{D(v_{i},c_{1})}\ket{D(v_{i},c_{2})}...\ket{D(v_{i},c_{t})}$ $$

We can then iteratively compare and swap pairs of registers to prepare the state $\ket{i}\ket{\min_{j\in[t]}D(v_{i},c_{j})}$.
If we apply the same procedure to $\ket{i}\ket{\psi_{i,1}}...\ket{\psi_{i,t}}$, then with probability at least $(1-2\Delta)^{t}$, the resulting state will be $\ket{i}\ket{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})}}$. So, the contents of the second register will be an estimate of the distance of the $i^{th}$ point to its closest center in the center set $C$. This further means that the following state can be prepared with probability at least $(1-2\Delta)^{Nt}$:(^8^88The state prepared is actually $\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\left(\alpha\ket{\min_{j\in[t]}{\tilde{ D}(v_{i},c_{j})}}+\beta\ket{G}\right)$ with $|\alpha|^{2}\geq(1-2\Delta)^{Nk}$. However, instead of working with this state, subsequent discussions become much simpler if we assume that $\ket{\Psi_{C}}$ is prepared with probability $|\alpha|^{2}$.)

$$ $\ket{\Psi_{C}}\equiv\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\ket{\min_{j\in[t]} {\tilde{D}(v_{i},c_{j})}}.$ $$

This quantum state can be used to find the approximate clustering cost of the center set $C$, which we discuss in the following subsection.
However, before we do that, let us summarise the main ideas of this subsection in the following lemma.

###### Lemma 3 .

There is a quantum algorithm that, with probability at least $(1-2\Delta)^{Nt}$, prepares the quantum state $\ket{\Psi_{C}}$ in time $O\left(\frac{Tt\log{\frac{1}{\Delta}}}{\varepsilon}\right)$.

###### Lemma 3 .

### A.2 Computing cost of clustering

Suppose we want to compute the $t$-means cost, $\Phi(V,C)\equiv\sum_{i=1}^{N}\min_{j\in[t]}{D^{2}(v_{i},c_{j})}$, of the clustering given by a $t$ center set $C\in\mathbb{R}^{t\times d}$.
We can prepare $m$ copies of the state $\ket{\Psi_{C}}$ and then estimate the cost of clustering by measuring $m$ copies of this quantum state and summing the squares of the second registers.
If $m$ is sufficiently large, we obtain a close estimate of $\Phi(V,C)$.
To show this formally, we will use the following Hoeffding tail inequality.

###### Theorem 5 (Hoeffding bound) .

Let $X_{1},...,X_{m}$ be independent, bounded random variables such that $X_{i}\in[a,b]$. Let $S_{m}=X_{1}+...+X_{m}$ . Then for any $\theta>0$, we have:

$$ $\mathbf{Pr}[|S_{m}-\mathbf{E}[S_{m}]|\geq\theta]\leq 2\cdot e^{\frac{-2\theta^ {2}}{m(b-a)^{2}}}.$ $$

Let $X_{1},...,X_{m}$ denote the square of the measured value of the second register in $\ket{\Psi_{C}}$.
These are random values in the range $[1,\zeta^{2}]$, where $\zeta=\max_{i,j}\tilde{D}(v_{i},v_{j})\in(1\pm\varepsilon)\cdot\max_{i,j}{D(v_
{i},v_{j})}$.
First, we note the expectation of these random variables equals $\frac{\tilde{\Phi}(V,C)}{N}$, where $\tilde{\Phi}(V,C)\equiv\sum_{i=1}^{N}\min_{j\in[t]}\tilde{D}(v_{i},c_{j})^{2}
\in(1\pm\varepsilon)^{2}\cdot\Phi(V,C)$.
We define the variable $S_{m}=X_{1}+X_{2}+...+X_{m}$ and apply the Hoeffding bound on these bounded random variables to get a concentration result that can then be used.

###### Lemma 4 .

Let $\alpha_{m}=S_{m}\cdot\frac{N}{m}$. If $m=\frac{\zeta^{2}\ln{(10N)}}{\varepsilon^{2}}$, then we have:

$$ $\mathbf{Pr}[\alpha_{m}\in(1\pm\varepsilon)\cdot\tilde{\Phi}(V,C)]\geq 1-\frac{ 1}{5N}.$ $$

###### Proof.

We know that $\mathbf{E}[S_{m}]=\frac{m}{N}\cdot\tilde{\Phi}(V,C)$
From the Hoeffding tail inequality, we get the following:

$$ $\displaystyle\mathbf{Pr}[|S_{m}-\mathbf{E}[S_{m}]|\geq\varepsilon\cdot\mathbf{ E}[S_{m}]]$ $\displaystyle\leq$ $\displaystyle 2\cdot e^{\frac{-2\varepsilon^{2}\mathbf{E}[S_{m}]^{2}}{m\zeta^{ 2}}}$ $\displaystyle=$ $\displaystyle 2\cdot e^{\frac{-2\varepsilon^{2}m}{\zeta^{2}}\cdot\left(\frac{ \tilde{\Phi}(V,C)}{N}\right)^{2}}$ $\displaystyle\leq$ $\displaystyle 2\cdot e^{-\ln{(10N)}}\leq\frac{1}{5N}.$ $$

This implies that:

$$ $\mathbf{Pr}[|\alpha_{m}-\tilde{\Phi}(V,C)|\geq\varepsilon\cdot\tilde{\Phi}(V,C )]\leq\frac{1}{5N}.$ $$

This completes the proof of the lemma.
∎

So, conditioned on having $m$ copies of the state $\ket{\Psi_{C}}$, we get an estimate of the clustering cost within a relative error of $(1\pm\varepsilon)^{3}$ with probability at least $(1-1/5N)$. Removing the conditioning, we get the same with probability at least $(1-2\Delta)^{Ntm}\cdot(1-1/5N)$.

###### Theorem 5 (Hoeffding bound) .

###### Lemma 4 .

###### Proof.

### A.3 D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling

$D^{2}$-sampling from the point set $V$ with respect to a center set $C\in\mathbb{R}^{t\times d}$ with $t$ centers, samples $v_{i}$ with probability proportional to $\min_{j\in[t]}{D^{2}(v_{i},c_{j})}$.
Let us see if we can use our state $\ket{\Psi_{C}}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\ket{\min_{j\in[t]}{
\tilde{D}(v_{i},c_{j})}}$ is useful to perform this sampling.
If we can pull out the value of the second register as the amplitude, then the measurement will give us close to $D^{2}$-sampling.
This is possible since we have an estimate of the clustering cost from the previous subsection.
We can use controlled rotations on an ancilla qubit to prepare the state:

$$ $\ket{\Psi_{sample}}\equiv\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\left(\beta_{i }\ket{0}+\sqrt{1-\beta_{i}}\ket{1}\right),$ $$

where $\beta_{i}=\frac{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})}}{\sqrt{2\cdot\tilde{
\Phi}(V,C)}}$.
So, the probability of measurement of $(i,0)$ is $\frac{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})^{2}}}{2\cdot\tilde{\Phi}(V,C)}$. Since we do rejection sampling (ignoring $(.,1)$’s that are sampled with probability $\frac{1}{2}$), we end up sampling with a distribution where the probability of sampling $i$ is $\frac{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})^{2}}}{\tilde{\Phi}(V,C)}\in(1\pm
\varepsilon)^{4}\cdot\frac{\min_{j\in[t]}{D(v_{i},c_{j})^{2}}}{\Phi(V,C)}$.
This means that points get sampled with a probability close to the actual $D^{2}$-sampling probability. As mentioned earlier, this is sufficient for the approximation guarantees of [] to hold.
We summarise the observations of this section in the next lemma.

###### Lemma 5 .

Let $0<\delta\leq 1$ and let $m=\frac{\zeta^{2}\ln{10N}}{\delta^{2}}$.
Given a dataset $V\in\mathbb{R}^{N\times d}$ and a center set $C\in\mathbb{R}^{t\times d}$ in the QRAM model, there is a quantum algorithm that runs in time $\tilde{O}\left(\frac{Ttm\log{\frac{1}{\Delta}}}{\delta}\right)$ and with probability at least $(1-1/5N)\cdot(1-2\Delta)^{Ntm}$ outputs a sample with
$\mathcal{D}^{2}$ distribution such that the distance function $\mathcal{D}$ is $\delta$-close to $D$.

###### Proof.

The proof follows from Lemma 3 and the preceding discussion on cost estimation in the previous subsection.
∎

The above lemma says that for $\Delta=\tilde{O}(\frac{\varepsilon^{2}}{Nt\zeta^{2}})$, we obtain the required $D^{2}$-sample with probability $\geq(1-1/N)$.
We can now give proof of Theorem 1 assembling the quantum tools for $D^{2}$-sampling in this section.
Note that we obtain samples from $\mathcal{D}^{2}$ distribution such that $\mathcal{D}$ is close to $D$.
Using our quantum sampling algorithm to pick $k$ centers in $k$ rounds is not the classical $k$-means++ algorithm that requires that the samples be picked with $D^{2}$ distribution and not something close. The resulting algorithm, however, has been analysed in the literature. If the sampling distribution is within $(1\pm\varepsilon)$ of the $D^{2}$ distribution, then the algorithm is called the noisy-$k$-means++ algorithm [, ]. We will use the following result in [].

###### Theorem 6 (Theorem 1.1 in [ 17 ] ) .

The noisy-$k$-means++ algorithm is $O(\log{k})$ approximate, in expectation.

Theorem 1, restated below, follows from the above theorem and Lemma 5.
Note that the running time of the quantum implementation is $\sum_{t}\tilde{O}(\zeta^{2}t)=\tilde{O}(\zeta^{2}k^{2})$.

###### Theorem (Restatement of Theorem 1 ) .

There is a quantum implementation of $k$-means++ that runs in time $\tilde{O}(\zeta^{2}k^{2})$ and gives an $O(\log{k})$ factor approximate solution for the $k$-means problem with a probability of at least $0.99$. Here, $\tilde{O}$ hides $\log^{2}{(Nd)}$ and $\log^{2}{(kd)}$ terms.(^9^99The output of $k$-means++ is a subset of data points, which can be stated as a subset of indices. If the description of centers is the expected output, the running time will also include a factor of $d$.)

###### Lemma 5 .

###### Proof.

###### Theorem 6 (Theorem 1.1 in [ 17 ] ) .

###### Theorem (Restatement of Theorem 1 ) .

## Appendix B Quantum Inspired D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -Sampling

Ewin Tang’s thesis [] nicely categorizes quantum machine learning (QML) algorithms into ones where the quantum advantage (i.e., polylogarithmic running time) is solely because of the quantum access to the data and the ones where it is not (e.g., the HHL algorithm []).
Quantum access means that a (normalized) data vector $\vec{v}\coloneqq(\vec{v}(1),...,\vec{v}(d))^{T}$ can be loaded onto the quantum workspace in $O(1)$ time as $\ket{v}=\vec{v}(1)\ket{1}+...+\vec{v}(d)\ket{d}$.
In this section, we will use the overhead arrow notation for vectors since we would want to distinguish between vectors and their approximations using different overhead symbols.
One implication of having the quantum state $\ket{v}$ is that $i$ can be sampled with probability $\frac{\vec{v}(i)^{2}}{\sum_{j}\vec{v}(j)^{2}}$ by simply making a measurement in the standard basis.
One of the key insights in Tang’s Thesis [] is that QML algorithms that benefit solely from the quantum data access can be ‘dequantized’ if we work with an appropriate classical counterpart of the quantum data access that allows sampling access of the kind described above (in addition to the classical access to the data).
This is called sample-query access, or SQ access in short.
Dequantization means that a classical algorithm can simulate the steps of the quantum algorithm using SQ access with similar running time dependence.
The nice property of the SQ access data structures is that it can be constructed classically in linear time.
Note that linear time means that we lose the quantum advantage.
However, if we keep this aside and count the SQ access data structure construction as preprocessing, the remaining computation has a similar advantage as that of the quantum algorithm.
This is a level playing field with the QML algorithm that allows fair comparison, specifically given that setting up quantum access (called QRAM) also takes linear time.
Such classical algorithms obtained by dequantization are called quantum inspired algorithms.
Our quantum algorithms based on $D^{2}$-sampling fall into the category of QML algorithms that can be dequantized in Tang’s SQ access model.
This section discusses the quantum-inspired classical algorithm we obtain by dequantizing our quantum algorithm using the SQ access data structure.
Naturally, we start the discussion by looking at the definition of the SQ access data structure.

###### Definition 3 (Query access, Definition 1.1 in [ 30 ] ) .

For a vector $\vec{v}\in\mathbb{C}^{n}$, we have $Q(\vec{v})$, query access to $\vec{v}$, if for all $i\in[n]$, we can query for $\vec{v}(i)$.
The time (cost) of such a query is denoted by $\mathbf{q}(\vec{v})$.

###### Definition 4 (SQ-access to a vector, Definition 1.2 in [ 30 ] ) .

For a vector $\vec{v}\in\mathbb{C}^{n}$, we have $SQ(\vec{v})$, sampling and query access to $\vec{v}$, if we can:

- 1.
query for entries in $\vec{v}$ as in $Q(\vec{v})$. The time cost is denoted by $\mathbf{q}(\vec{v})$.
- 2.
obtain independent samples $i\in[n]$ where the probability of sampling $i$ is $\frac{\vec{v}(i)^{2}}{\norm{\vec{v}}^{2}}$. This distribution is denoted by $\mathcal{D}_{\vec{v}}$.
The time cost is denoted by $\mathbf{s}(\vec{v})$.
- 3.
query for $\norm{\vec{v}}$. The time cost is denoted by $\mathbf{n}(\vec{v})$.

Let $\mathbf{sq}(\vec{v})\equiv\max{\{\mathbf{q}(\vec{v}),\mathbf{s}(\vec{v}),
\mathbf{n}(\vec{v})\}}$ denote the time cost of SQ-access.

In quantum computation an $n$-dimensional state vector $\ket{v_{0}}$ evolves using unitary transformations $\ket{v_{1}}=U\ket{v_{0}}$, where $U$ is a $2^{n}\times 2^{n}$ unitary matrix.
To be able to simulate such operations within the SQ-access model, we must make sure that SQ-access is available for intermediate states of the quantum operation. This means that we must (i) define an appropriate SQ-access notion for matrices and (ii) show closure properties for SQ-access.
Fortunately, this has been shown in Tang’s Thesis [].
A more flexible version of SQ-access is needed to show closure properties. This is called oversampling SQ-access. SQ-access for matrices and oversampling SQ-access is defined below.

###### Definition 5 (SQ-access to a matrix) .

For a matrix $A\in\mathbb{C}^{m\times n}$, we have $SQ(A)$ if we have $SQ(A(i,.))$ for all $i\in[m]$ (i.e., SQ-access to the row vectors), and $SQ(\vec{a})$, where $\vec{a}$ is the vector of row norms of $A$ (i.e., $\vec{a}=(\norm{A(1,.)},...,\norm{A(m,.)})^{T}$). The complexity of operations on a matrix is defined as: $\mathbf{q}(A)\coloneqq\max{\{\mathbf{q}(A(i,.)),\mathbf{q}(\vec{a})\}}$, $\mathbf{n}(A)\coloneqq\mathbf{n}(\vec{a})$, $\mathbf{s}(A)\coloneqq\max{\{\mathbf{s}(A(i,.)),\mathbf{s}(\vec{a})\}}$, and $\mathbf{sq}(A)\coloneqq\max{\{\mathbf{q}(A),\mathbf{n}(A),\mathbf{s}(A)\}}$

###### Definition 6 (Oversampling and query access: Definition 4.4 in [ 30 ] ) .

For $\vec{v}\in\mathbb{C}^{n}$ and $\phi\geq 1$, we have $SQ_{\phi}(\vec{v})$, $\phi$-oversampling and query access to $\vec{v}$ ($\phi$-OSQ access, in short), if we have $Q(\vec{v})$ and $SQ(\tilde{v})$ for $\tilde{v}\in\mathbb{C}^{n}$ a vector satisfying (i) $\norm{\tilde{v}}^{2}=\phi\norm{\vec{v}}^{2}$ and (ii) $|\tilde{v}(i)|^{2}\geq|\vec{v}(i)|^{2}$ for all $i\in[n]$. Denote $\mathbf{s}_{\phi}(\vec{v})\coloneqq\mathbf{s}(\tilde{v})$, $\mathbf{q}_{\phi}(\vec{v})\coloneqq\mathbf{q}(\tilde{v})$, $\mathbf{n}_{\phi}(\vec{v})\coloneqq\mathbf{n}(\tilde{v})$, and $\mathbf{sq}_{\phi}(\vec{v})\coloneqq\max{\{\mathbf{s}_{\phi}(\vec{v}),\mathbf{
q}_{\phi}(\vec{v}),\mathbf{q}(\vec{v}),\mathbf{n}_{\phi}(\vec{v})\}}$.

Note that $\mathcal{D}_{\tilde{v}}$ for $\tilde{v}$ can be seen as an oversampling distribution of $\mathcal{D}_{\vec{v}}$ by a factor of $\phi$.
This is because the sampling probability of index $i$ is given by

$$ $\mathcal{D}_{\tilde{v}}(i)=\frac{|\tilde{v}(i)|^{2}}{\norm{\tilde{v}}^{2}}= \frac{|\tilde{v}(i)|^{2}}{\phi\norm{\vec{v}}^{2}}\geq\frac{|\vec{v}(i)|^{2}}{ \phi\norm{\vec{v}}^{2}}=\frac{\mathcal{D}_{\vec{v}}(i)}{\phi}.$ $$

$\phi$-OSQ access also allows one to sample from $\mathcal{D}_{\vec{v}}$ and find $\norm{\vec{v}}$ by paying a small time penalty (as a function of $\phi$ and allowed error). This is captured in the following lemma from [].

###### Lemma 6 (Lemma 4.5 in [ 30 ] ) .

Let $0<\delta\leq 1$ and $0<\varepsilon\leq 1$.
Suppose we are given $SQ_{\phi}(\vec{v})$.
Then we can sample from $\mathcal{D}_{\vec{v}}$ with proability at least $(1-\delta)$ in time $O\left(\phi\cdot\mathbf{sq}_{\phi}(\vec{v})\cdot\log{\frac{1}{\delta}}\right)$. We can also estimate $\norm{\vec{v}}$ to within $(1\pm\varepsilon)$ factor with probability at least $(1-\delta)$ in $O\left(\frac{\phi}{\varepsilon^{2}}\cdot\mathbf{sq}_{\phi}(\vec{v})\cdot\log{
\frac{1}{\delta}}\right)$ time.

We can also define oversampling query access to a matrix.

###### Definition 7 (Oversampling and query access for a matrix) .

For a matrix $A\in\mathbb{C}^{N\times d}$ and $\phi\geq 1$, we have $SQ_{\phi}(A)$, $\phi$-oversampling and query access to $A$ ($\phi$-OSQ access, in short), if we have $Q(A)$ and $SQ(\tilde{A})$ for $\tilde{A}\in\mathbb{C}^{N\times d}$ satisfying $\norm{\tilde{A}}^{2}_{F}=\phi\norm{A}^{2}_{F}$ and $|\tilde{A}(i,j)|^{2}\geq|A(i,j)|^{2}$ for all $(i,j)\in[N]\times[d]$. The $\phi$-oversampling complexity is given by $\mathbf{sq}_{\phi}(A)=\max{\{\mathbf{q}(A),\mathbf{sq}(\tilde{A})\}}$.

There is a simple tree-based data structure that supports sample-query access for vectors and matrices. The details can be found in []. Here, we give a summary of the running time of various operations.

###### Lemma 7 (Remark 4.12 in [ 30 ] ) .

There is a data structure for storing a vector $\vec{v}\in\mathbb{R}^{n}$ with supporting the following operations:

- 1.
Reading and updating an entry of $\vec{v}$ in $O(\log{n})$ time. So, $\mathbf{q}(\vec{v})=O(\log{n})$.
- 2.
Finding $\norm{\vec{v}}^{2}$ in $O(1)$ time. So, $\mathbf{n}(\vec{v})=O(1)$.
- 3.
Sampling from $\mathcal{D}_{\vec{v}}$ in $O(\log{n})$ time. So, $\mathbf{s}(\vec{v})=O(\log{n})$

It follows from the above that $\mathbf{sq}(\vec{v})=O(\log{n})$. Moreover, the time to construct the data structure is $O(n)$.

A similar lemma holds for storing matrices.

###### Lemma 8 (Remark 4.12 in [ 30 ] ) .

For any matrix $A\in\mathbb{R}^{m\times n}$, let $\vec{A}$ denote the vector of row norms of $A$, i.e., the $i^{th}$ entry of $\vec{A}$, denoted by $\vec{A}(i)$, is $\norm{A(i,.)}$.
There is a data structure for storing a matrix $A\in\mathbb{R}^{m\times n}$ with supporting the following operations:

- 1.
Reading and updating an entry of $A$ in $O(\log{mn})$ time.
- 2.
Finding $\vec{A}(i)$ in $O(\log{m})$ time.
- 3.
Finding $\norm{A}_{F}^{2}$ in $O(1)$ time.
- 4.
Sampling from $\mathcal{D}_{\vec{A}}$ and $\mathcal{D}_{A(i,.)}$ in $O(\log{mn})$ time.

The time to construct the data structure is $O(mn)$.

Let us now see how $D^{2}$-sampling can be done using the above SQ-access data structures.
Let $V\in\mathbb{R}^{N\times d}$ be the data matrix with the row vectors $\vec{v}_{1},...,\vec{v}_{N}$ as data points.
Let $C\in\mathbb{R}^{m\times d}$ be the center set matrix with the centers being the $i$ row vectors $\vec{c}_{1},...,\vec{c}_{m}$.
Consider $D^{2}$-sampling a point from $V$ with respect to the center set $C$.
The probability of sampling the $i^{th}$ point is given by:

$$ $\mathbf{Pr}[i]=\frac{\min_{t\in[m]}\norm{\vec{v}_{i}-\vec{c}_{t}}^{2}}{\sum_{s \in[N]}\min_{t\in[m]}\norm{\vec{v}_{s}-\vec{c}_{t}}^{2}}$ $$

Consider vectors $\vec{w}_{1},...,\vec{w}_{m}\in\mathbb{R}^{N}$ defined as:

$$ $\vec{w}_{j}\coloneqq(\norm{\vec{v}_{1}-\vec{c}_{j}},\norm{\vec{v}_{2}-\vec{c}_ {j}},...,\norm{\vec{v}_{N}-\vec{c}_{j}})^{T}.$ $$

The first step towards $D^{2}$-sampling using the SQ-access data structure for matrices $V$ and $C$ is to check whether SQ-access for vectors $\vec{w}_{1},...,\vec{w}_{m}$ can be made available.
In the next lemma, we show that $\phi$-oversampling SQ-access for an appropriately chosen value of $\phi$ is possible.

###### Lemma 9 .

Let $\alpha$-oversampling SQ access to the data matrix $V\in\mathbb{R}^{N\times d}$ be available as $SQ_{\alpha}(V)$ and let $\beta$ oversampling SQ access be available for a vector $\vec{c}\in\mathbb{R}^{d}$ as $SQ_{\beta}(\vec{c})$. The we have $\gamma$ oversampling SQ access, $SQ_{\gamma}(\vec{w})$, for vector $\vec{w}\coloneqq(\norm{V(1,.)-\vec{c}},...,\norm{V(N,.)-\vec{c}})^{T}$ with $\gamma=\frac{2}{\norm{\vec{w}}^{2}}\cdot\left(\alpha\norm{V}^{2}+N\beta\norm{
\vec{c}}^{2}\right)$. Moreover, the complexity is related as $\mathbf{sq}_{\gamma}(\vec{w})\leq d\cdot(\mathbf{sq}_{\alpha}(V)+\mathbf{sq}_{
\beta}(\vec{c}))$.

###### Proof.

$SQ_{\alpha}(V)$ means that $Q(V)$ and $SQ(\tilde{V})$ are available for a matrix $\tilde{V}$ such that (i) $\norm{\tilde{V}}^{2}=\alpha\norm{V}^{2}$ and (ii) $|\tilde{V}(i,j)|\geq|V(i,j)|$ for all $(i,j)\in[N]\times[d]$.
Similarly, $SQ_{\beta}(\vec{c})$ means that $Q(\vec{c})$ and $SQ(\tilde{c})$ are available for a vector $\tilde{c}$ such that (i) $\norm{\tilde{c}}^{2}=\beta\norm{\vec{c}}^{2}$ and (ii) $|\tilde{c}(t)|\geq|\vec{c}(t)|$ for all $t\in[d]$.
To show $SQ_{\gamma}(\vec{w})$, we first need to show $Q(\vec{w})$. This is simple since the $j^{th}$ coordinate of $\vec{w}$ can be found using all the coordinates of $V(j,.)$ and $\vec{c}$ which are available since we have $Q(V)$ and $Q(\vec{c})$.
So, $\mathbf{q}(\vec{w})=d\cdot(\mathbf{q}(V)+\mathbf{q}(\vec{c}))$.
We now need to show $SQ(\tilde{w})$ for an appropriately chosen $\tilde{w}$. We consider:

$$ $\tilde{w}(j)=\sqrt{2\left(\norm{\tilde{V}(j,.)}^{2}+\norm{\tilde{c}}^{2}\right )}.$ $$

Before we see why SQ access to $\tilde{w}$ is possible, let us first observe that $SQ(\tilde{w})$ would give a $\gamma$-oversampling SQ access to $\vec{w}$. Towards this, we note that

$$ $\norm{\tilde{w}}^{2}=\sum_{j}|\tilde{w}(j)|^{2}=2\sum_{j}\left(\norm{\tilde{V} (j,.)}^{2}+\norm{\tilde{c}}^{2}\right)=\frac{2}{\norm{\vec{w}}^{2}}\left( \alpha\norm{V}^{2}+N\beta\norm{\vec{c}}^{2}\right)\cdot\norm{\vec{w}}^{2}= \gamma\norm{\vec{w}}^{2}.$ $$

Moreover, we have:

$$ $|\vec{w}(j)|=\norm{V(j,.)-\vec{c}}\leq\norm{V(j,.)}+\norm{\vec{c}}\leq\norm{ \tilde{V}(j,.)}+\norm{\tilde{c}}\leq\sqrt{2\left(\norm{\tilde{V}(j,.)}^{2}+ \norm{\tilde{c}}^{2}\right)}=|\tilde{w}(j)|.$ $$

Finally, we must show $SQ(\tilde{w})$.
Coordinates of $\tilde{w}$ can be found using $\norm{\tilde{V}(j,.)}$’s and $\norm{\tilde{c}}$ which are available through oversampling access to $V$ and $\vec{c}$. $\norm{\tilde{w}}^{2}$ is available through $\norm{\tilde{V}}^{2}$ and $\norm{\tilde{c}}^{2}$.
Finally, $\mathcal{D}_{\tilde{w}}$ can be sampled by first selecting $\tilde{V}$ with probability $\frac{\norm{\tilde{V}}^{2}}{\norm{\tilde{V}}^{2}+N\norm{\tilde{c}}^{2}}$ or $\tilde{c}$ with probability $\frac{N\norm{\tilde{c}}^{2}}{\norm{\tilde{V}}^{2}+N\norm{\tilde{c}}^{2}}$ and then sampling from $\mathcal{D}_{\tilde{v}}$ or $\mathcal{U}_{[N]}$ (uniform distribution over $[N]$), respectively. Given this, the probability of sampling $i$ works out to be:

$$ $\mathbf{Pr}[i]=\frac{\norm{\tilde{V}}^{2}}{\norm{\tilde{V}}^{2}+N\norm{\tilde{ c}}^{2}}\cdot\frac{\norm{\tilde{V}(i,.)}^{2}}{\norm{\tilde{V}}^{2}}+\frac{N \norm{\tilde{c}}^{2}}{\norm{\tilde{V}}^{2}+N\norm{\tilde{c}}^{2}}\cdot\frac{1} {N}=\frac{\tilde{w}(i)^{2}}{\norm{\tilde{w}}^{2}}=\mathcal{D}_{\tilde{w}}(i).$ $$

From the above discussion, it follows that $\mathbf{sq}_{\gamma}(\vec{w})\leq d\cdot(\mathbf{sq}_{\alpha}(V)+\mathbf{sq}_{
\beta}(\vec{c}))$.
∎

The previous lemma gives a way to obtain oversampling query access to distances from a single center.
In $D^{2}$-sampling, we have a set $C$ of $m$ centers, and we sample a point based on the squared distance from the nearest center.
To enable this operation in the SQ framework, we need sample and query access to a vector where the $i^{th}$ coordinate is the distance of $i^{th}$ data point to the nearest center in $C$.
In the following lemma, we show how to do this. The proof closely follows Lemma 4.6 of [].

###### Lemma 10 .

Let $\vec{u}_{1},...,\vec{u}_{m}\in(\mathbb{R}^{\geq 0})^{N}$ and $\phi_{1},...,\phi_{m}\in\mathbb{R}^{\geq 1}$.
Let

$$ $\vec{w}\coloneqq\left(\min_{i\in[m]}{\vec{u}_{i}(1)},\min_{i\in[m]}{\vec{u}_{i }(2)},...,\min_{i\in[m]}{\vec{u}_{i}(N)}\right)^{T}.$ $$

Suppose we have $SQ_{\phi_{i}}(\vec{u}_{i})$ for every $i\in[m]$.
Then there is a $\phi$-oversampling query access, $SQ_{\phi}(\vec{u})$ for $\phi=\frac{\sum_{i\in[m]}\phi_{i}\norm{\vec{u}_{i}}^{2}}{m\norm{\vec{w}}^{2}}$. The complexity is give by: $\mathbf{q}(\vec{w})=\sum_{i\in[m]}\mathbf{q}(\vec{u}_{i})$, $\mathbf{q}_{\phi}(\vec{w})=\sum_{i\in[m]}\mathbf{q}_{\phi_{i}}(\vec{u}_{i})$, $\mathbf{n}_{\phi}(\vec{w})=\sum_{i\in[m]}\mathbf{n}_{\phi_{i}}(\vec{u}_{i})$. The sampling complexity is $\mathbf{s}_{\phi}(\vec{w})=\max_{i\in[m]}{\mathbf{s}_{\phi_{i}}(\vec{u}_{i})}$, after a single-time pre-processing cost of $\sum_{i\in[m]}\mathbf{n}_{\phi_{i}}(\vec{u}_{i})$.

###### Proof.

We need to show $Q(\vec{w})$ and $SQ(\tilde{w})$ for an appropriate $\tilde{w}$. We can respond to each coordinate of $\vec{w}$ by querying $\vec{u}_{i}$’s and taking the minimum. This gives $\mathbf{q}(\vec{w})=\sum_{i\in[m]}\mathbf{q}(\vec{u}_{i})$.
Let $\tilde{u}_{1},...,\tilde{u}_{m}$ be the vectors corresponding to the oversampling access $SQ_{\phi_{1}}(\vec{u}_{1}),...,SQ_{\phi_{m}}(\vec{u}_{m})$.
For oversampling query access, consider:

$$ $\tilde{w}(i)=\sqrt{\frac{1}{m}\cdot\sum_{j\in[m]}|\tilde{u}_{j}(i)|^{2}}.$ $$

Each coordinate of $\tilde{w}$ can be found by querying $\tilde{u}_{i}$’s and computing the above expression. So, we have $Q(\tilde{w})$ and $\mathbf{q}_{\phi}(\vec{w})=\sum_{i\in[m]}\mathbf{q}_{\phi_{i}}(\vec{u}_{i})$. The norm of $\tilde{w}$ is given by:

$$ $\norm{\tilde{w}}^{2}=\frac{1}{m}\sum_{i\in[N]}\sum_{j\in[m]}|\tilde{u}_{j}(i)| ^{2}=\frac{1}{m}\sum_{j\in[m]}\norm{\tilde{u}_{j}}^{2}=\frac{1}{m}\sum_{j\in[m ]}\phi_{j}\norm{\vec{u}_{j}}^{2}=\left(\frac{\sum_{j\in[m]}\phi_{j}\norm{\vec{ u}_{j}}^{2}}{m\norm{\vec{w}}^{2}}\right)\norm{\vec{w}}^{2}=\phi\norm{\vec{w}}^ {2}.$ $$

So, we can query the norm of $\tilde{w}$ by querying the norms of $\tilde{u}_{j}$’s.
Hence, $\mathbf{n}_{\phi}(\vec{w})=\sum_{j\in[m]}\mathbf{n}_{\phi_{j}}(\vec{u}_{j})$.
Now, we need to show that $\tilde{w}(i)$ upper bounds $\vec{w}(i)$.
We have:

$$ $\tilde{w}(i)=\sqrt{\frac{1}{m}\sum_{j\in[m]}|\tilde{u}_{j}(i)|^{2}}\geq\sqrt{ \frac{1}{m}\sum_{j\in[m]}\min_{j^{\prime}\in[m]}|\tilde{u}_{j^{\prime}}(i)|^{2 }}=\min_{j^{\prime}\in[m]}\tilde{u}_{j^{\prime}}(i)\geq\min_{j^{\prime}\in[m]} \vec{u}_{j}(i)=\vec{w}(i).$ $$

Now we describe how to sample from $\mathcal{D}_{\tilde{w}}$. We first query the norms of $\tilde{u}_{j}$.
We pick a $j\in[m]$ with probability $\frac{\norm{\tilde{u}_{j}}^{2}}{\sum_{j^{\prime}\in[m]}\norm{\tilde{u}_{j^{
\prime}}}^{2}}$, and then sample an index $i\in[N]$ from the distribution $\mathcal{D}_{\tilde{u}_{j}}$. The probability of sampling an index $i$ is given by:

$$ $\mathbf{Pr}[i]=\sum_{j\in[m]}\frac{\norm{\tilde{u}_{j}}^{2}}{\sum_{j^{\prime} \in[m]}\norm{\tilde{u}_{j^{\prime}}}^{2}}\cdot\frac{|\tilde{u}_{j}(i)|^{2}}{ \norm{\tilde{u}_{j}}^{2}}=\frac{\frac{1}{m}\sum_{j\in[m]}|\tilde{u}_{j}(i)|^{2 }}{\frac{1}{m}\sum_{j^{\prime}\in[m]}\norm{\tilde{u}_{j^{\prime}}}^{2}}=\frac{ |\tilde{w}(i)|^{2}}{\norm{\tilde{w}}^{2}}.$ $$

This is the correct probability for sampling from $\mathcal{D}_{\tilde{w}}$. Note that from the description, the sampling complexity is $\mathbf{s}_{\phi}(\vec{w})=\max_{i\in[m]}{\mathbf{s}_{\phi_{i}}(\vec{u}_{i})}$, after a single-time pre-processing cost of $\sum_{i\in[m]}\mathbf{n}_{\phi_{i}}(\vec{u}_{i})$.
∎

We now have all the basic tools to perform $D^{2}$-sampling in the SQ access model and evaluate its time complexity.
In the remaining discussion, we will assume that the origin is centered at one of the data points so that the maximum norm of a data vector is upper bounded by the maximum interpoint distance.

###### Theorem 7 .

Let $0<\delta\leq 1$.
Let $V\in\mathbb{R}^{N\times d}$ be the dataset and let $\vec{c}_{1},...,\vec{c}_{m}\in\mathbb{R}^{d}$ be $m$ arbitrary centers from the dataset and let $SQ(V)$ and $SQ(\vec{c}_{1}),...,SQ(\vec{c}_{m})$ be available.
Let $\zeta$ denote the aspect ratio of the dataset (i.e., the ratio of the maximum to minimum interpoint distance).
Then there is a classical algorithm that with probability at least $(1-\delta)$ outputs a sample from the $D^{2}$-distribution (i.e., $D^{2}$-sample) of $V$ with respect to centers $\{\vec{c}_{1},...,\vec{c}_{m}\}$ which runs in time $O(\zeta^{2}md(\log{\frac{1}{\delta}})\log{Nd})$.

###### Proof.

Let $d_{norm}$ be the maximum norm of a data point.
Let $d_{min}$ and $d_{max}$ be the minimum and maximum distance between two data points, respectively.
Then, we have $\zeta=\frac{d_{max}}{d_{min}}\geq\frac{d_{norm}}{d_{min}}$.
Let us define vectors $\vec{u}_{1},...,\vec{u}_{m}$ where $\vec{u}_{j}(i)\coloneqq\norm{V(i,.)-\vec{c}_{j}}$.
Using Lemma 9, we can enable $SQ_{\phi_{j}}(\vec{u}_{j})$ for every $j\in[m]$, for $\phi_{j}=\frac{2}{\norm{\vec{u_{j}}}^{2}}(\norm{V}^{2}+N\norm{\vec{c}_{j}}^{2})$.
Lemma 10 gives a way to define $SQ_{\phi}(\vec{w})$, where $\vec{w}$ is defined as $\vec{w}(i)\coloneqq\min_{j\in[m]}|\vec{u}_{j}(i)|$ and $\phi=\frac{\sum_{j\in[m]}\phi_{j}\norm{\vec{u}_{j}}^{2}}{m\norm{\vec{w}}^{2}}$.
Note that for the $\vec{w}$ defined, $\mathcal{D}_{\vec{w}}$ is precisely the $D^{2}$ distribution over $V$ with respect to centers $\vec{c}_{1},...,\vec{c}_{m}$ and from Lemma 6, we know how to use $SQ_{\phi}(\vec{w})$ to sample from $\mathcal{D}_{\vec{w}}$ with probability at least $(1-\delta)$ in time $O(\phi\cdot\mathbf{sq}_{\phi}(\vec{w})\cdot\log{\frac{1}{\delta}})$. For calculating the complexity, we first get a bound on $\phi$:

$$ $\phi=\frac{\sum_{j\in[m]}\phi_{j}\norm{\vec{u}_{j}}^{2}}{m\norm{\vec{w}}^{2}}= \frac{\sum_{j\in[m]}\frac{2}{\norm{\vec{u_{j}}}^{2}}(\norm{V}^{2}+N\norm{\vec{ c}_{j}}^{2})\norm{\vec{u}_{j}}^{2}}{m\norm{\vec{w}}^{2}}=\frac{2\norm{V}^{2}}{ \norm{\vec{w}}^{2}}+\frac{2N}{\norm{\vec{w}}^{2}}\frac{1}{m}\sum_{j\in[m]} \norm{\vec{c}_{j}}^{2}.$ $$

The $i^{th}$ coordinate of $\vec{w}$ is the least distance of $V(i,.)$ from a center. This gives us $\norm{\vec{w}}^{2}\geq(N-m)\cdot d_{min}^{2}$. Also, $\norm{V}^{2}\leq Nd_{norm}^{2}$ and $\norm{\vec{c}_{j}}^{2}\leq d_{norm}^{2}$. Putting these bounds in the above equation (and that $m\leq N/2$), we get $\phi\leq 8\frac{d_{norm}^{2}}{d_{min}^{2}}\leq 8\zeta^{2}$.
From Lemmas 7 and 8, we have $\mathbf{sq}(V)=O(\log{Nd})$, $\mathbf{sq}(\vec{c}_{j})=O(\log{d})$.
From Lemma 9, we have $\mathbf{sq}_{\phi_{j}}(\vec{u}_{j})=d\cdot(\mathbf{sq}(V)+\mathbf{sq}(\vec{c}_
{j}))=O(d\log{Nd})$ for every $j\in[m]$.
From Lemma 10, we have $\mathbf{sq}_{\phi}(\vec{w})=\sum_{j\in[m]}\mathbf{sq}_{\phi_{j}}(\vec{u}_{j})=
O(md\log{Nd})$.
So, from Lemma 6, the complexity of $D^{2}$-sampling a point is $O(\phi\cdot\mathbf{sq}_{\phi}(\vec{w})\cdot\log{\frac{1}{\delta}})=O(\zeta^{2}
md(\log{\frac{1}{\delta}})\log{Nd})$.
∎

The $k$-means++ algorithm is simply $D^{2}$-sampling $k$ centers iteratively while updating the list of centers.
We summarise the discussion of this section by giving the complexity of the $k$-means++ algorithm and the proof of Theorem 2. We restate the theorem for ease of reading.

###### Theorem (Restatement of Theorem 2 ) .

There is a classical implementation of $k$-means++ (which we call QI-$k$-means++) that runs in time $O(Nd)+\tilde{O}(\zeta^{2}k^{2}d)$. Moreover, with probability at least $0.99$, QI-$k$-means++ outputs $k$ centers that give $O(\log{k})$-approximation in expectation for the $k$-means problem.

###### Proof.

The $Nd$ term in the running time expression is for setting up $SQ(V)$, the sample-query access to the data matrix $V$. This may also be regarded as the preprocessing step. This is followed by a $D^{2}$-sampling $k$ centers iteratively while updating the set of centers by including the new centers sampled. From the previous theorem, we know that when we have $m$ centers, the cost of obtaining a $D^{2}$-sample is $O(\zeta^{2}md(\log{\frac{1}{\delta}})\log{Nd})$ with probability at least $(1-\delta)$.
To obtain an overall probability of success of $0.99$, we need to have a success probability of at least $1-O(1/k)$ in each iteration. This makes the overall $D^{2}$-sampling cost $O(\zeta^{2}k^{2}d\log{k}\log{Nd})$. So, the overall time complexity of QI-$k$-means++ is $O\left(Nd\right)+\tilde{O}(\zeta^{2}k^{2}d)$. Since with probability at least $0.99$, we obtain $k$ perfectly $D^{2}$-sampled center, the approximation guarantee remains the same as that of $k$-means++, which is $O(\log{k})$.(^10^1010What happens with the remaining 0.01 probability? $D^{2}$-sampling in the SQ model succeeding with probability at least $(1-\delta)$ means that with probability at least $(1-\delta)$ a center is sampled with $D^{2}$ distribution, and with the remaining probability, no point is sampled.)
∎

###### Definition 3 (Query access, Definition 1.1 in [ 30 ] ) .

###### Definition 4 (SQ-access to a vector, Definition 1.2 in [ 30 ] ) .

###### Definition 5 (SQ-access to a matrix) .

###### Definition 6 (Oversampling and query access: Definition 4.4 in [ 30 ] ) .

###### Lemma 6 (Lemma 4.5 in [ 30 ] ) .

###### Definition 7 (Oversampling and query access for a matrix) .

###### Lemma 7 (Remark 4.12 in [ 30 ] ) .

###### Lemma 8 (Remark 4.12 in [ 30 ] ) .

###### Lemma 9 .

###### Proof.

###### Lemma 10 .

###### Proof.

###### Theorem 7 .

###### Proof.

###### Theorem (Restatement of Theorem 2 ) .

###### Proof.

### B.1 Quantum inspired approximate D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling

Sampling from a distribution that is $\varepsilon$-close to the $D^{2}$ distribution is sufficient to obtain $O(\log{k})$ approximation.
This follows from the analysis of the noisy-$k$-means++ algorithm [].
This relaxation (without any approximation penalty beyond constant factors) allows us to shave off the factor of $d$ from the time complexity of the quantum-inspired algorithm.
However, as we will see in the analysis of this section, it is at the cost of worsening the dependence on the aspect ratio.
So, the algorithm in this section should be used for high-dimensional cases where the aspect ratio is small.
The factor of $d$ in the sampling component of QI-$k$-means++ appears because for enabling the oversampling and query access for the vector $\vec{w}\coloneqq(\norm{V(1,.)-\vec{c}},...,\norm{V(N,.)-\vec{c}})^{T}$ ($\vec{c}$ is a center), we need to access the $d$ coordinates of $V(i,.)$ and $\vec{c}$ to be able to answer the query access to the $i^{th}$ index of $\vec{w}$, which is $\norm{V(i,.)-\vec{c}}$.
Note that $\norm{V(i,.)-\vec{c}}=\sqrt{\norm{V(i,.)}^{2}+\norm{\vec{c}}^{2}-2\langle V(i,
.),\vec{c}\rangle}$, where $\langle V(i,.),\vec{c}\rangle$ denotes the dot product. The terms $\norm{V(i,.)}^{2}$ and $\norm{\vec{c}}^{2}$ are available in $O(1)$ time. The dot product term is costly and adds a factor of $d$. So, the key idea in eliminating the factor of $d$ is to replace the dot product term with an appropriate approximation that requires a much smaller access time. This uses the standard sampling trick - randomly sample a subset of indices, compute the dot product restricted to these indices of the subset, and scale (see Lemma 5.17 in []). The number of indices we need to sample to ensure that the resulting approximation is within $(1\pm\varepsilon)$ of $\norm{V(i,.)-\vec{c}}$ for every $i$ with probability at least $(1-\delta)$ is $\tilde{O}\left(\frac{\log{(1/\delta)}}{\varepsilon^{2}}\cdot\frac{\norm{\vec{c
}}^{2}\cdot\norm{V(i,.)}^{2}}{\norm{V(i,.)-\vec{c}}^{2}}\right)=\tilde{O}\left
(\frac{\zeta^{4}\log{1/\delta}}{\varepsilon^{2}}\right)$. The remaining steps are the same as those in the exact version, except that $d$ is replaced with $\tilde{O}\left(\frac{\zeta^{4}\log{1/\delta}}{\varepsilon^{2}}\right)$. This results in the overall running time $O(Nd)+\tilde{O}\left(\zeta^{6}k^{2}\right)$.

## Appendix C A Quantum Approximation Scheme (proof of Theorem 4 )

We start the discussion with the $D^{2}$-sampling method.
In particular, we would like to check the robustness of the approximation guarantee provided by the $D^{2}$-sampling method against errors in estimating the distances between points.
We will show that the $D^{2}$-sampling method gives a constant pseudo-approximation even under sampling errors.

### C.1 Pseudoapproximation using D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling

Let the matrix $V\in\mathbb{R}^{N\times d}$ denote the dataset, where row $i$ contains the $i^{th}$ data point $v_{i}\in\mathbb{R}^{d}$.
Let the matrix $C\in\mathbb{R}^{t\times d}$ denote any $t$-center set, where row $i$ contains the $i^{th}$ center $c_{i}\in\mathbb{R}^{d}$ out of the $t$ centers.
Sampling a data point using the $D^{2}$ distribution w.r.t. (short for with respect to) a center set $C$ means that the datapoint $v_{i}$ gets sampled with probability proportional to the squared distance to its nearest center in the center set $C$.
This is also known as $D^{2}$ sampling w.r.t. center set $C$.
More formally, data points are sampled using the distribution
$\left(\frac{D^{2}(v_{1},C)}{\sum_{j}D^{2}(v_{j},C)},...,\frac{D^{2}(v_{N},C)}{
\sum_{j}D^{2}(v_{j},C)}\right)$, where $D^{2}(v_{j},C)\equiv\min_{c\in C}D^{2}(v_{j},c)$.
For the special case $C=\emptyset$, $D^{2}$ sampling is the same as uniform sampling.
The $k$-means++ seeding algorithm starts with an empty center set $C$ and, over $k$ iterations, adds a center to $C$ in every iteration by $D^{2}$ sampling w.r.t. the current center set $C$.
It is known from the result of [] that this $k$-means++ algorithm above gives an $O(\log{k})$ approximation in expectation.
It is also known from the result of [] that if $2k$ centers are sampled, instead of $k$ (i.e., the for-loop runs from $1$ to $2k$), the cost with respect to these $2k$ centers is at most some constant times the optimal $k$-means cost.
Such an algorithm is called a pseudo approximation algorithm.
Such a pseudo approximation algorithm is sufficient for the approximation scheme of [].
So, we will quantize the following constant factor pseudo-approximation algorithm.

Figure: Algorithm 3 A pseudo-approximation algorithm based on $D^{2}$-sampling.

In the quantum simulation of the above sampling procedure, there will be small errors in the sampling probabilities in each iteration. We need to ensure that the constant approximation guarantee of the above procedure is robust against small errors in the sampling probabilities owing to errors in distance estimation.
We will work with a relative error of $(1\pm\delta)$ for small $\delta$.
Following is a crucial lemma from [] that we will need to show the pseudo approximation property of Algorithm 1.

###### Lemma 11 (Lemma 3.2 in [ 3 ] ) .

Let $A$ be an arbitrary optimal cluster, and let $C$ be an arbitrary set of centers. Let $c$ be a center chosen from $A$ with $D^{2}$-sampling with respect to $C$.
Then $\mathbf{E}[cost(A,C\cup\{c\})]\leq 8\cdot OPT(A)$.

The above lemma is used as a black box in the analysis of Algorithm 1 in [].
The following version of the lemma holds for distance estimates with a relative error of $(1\pm\delta)$ and gives a constant factor approximation guarantee.
Since Lemma 11 is used as a black box in the analysis of Algorithm 1, replacing this lemma with Lemma 12 also gives a constant factor approximation to the $k$-means objective. We will use the following notion of the closeness of two distance functions.

###### Definition 8 .

A distance function $D_{1}$ is said to be $\delta$-close to distance function $D_{2}$, denoted by $D_{1}\sim_{\delta}D_{2}$, if for every pair of points $x,y\in\mathbb{R}^{d}$, $D_{1}(x,y)\in(1\pm\delta)\cdot D_{2}(x,y)$.(^11^1111We use the notation that for positive reals $P,Q$, $P\in(1\pm\delta)\cdot Q$ if $(1-\delta)\cdot Q\leq P\leq(1+\delta)\cdot Q$.)

###### Lemma 12 .

Let $0<\delta\leq 1/2$.
Let $A$ be an arbitrary optimal cluster and $C$ be an arbitrary set of centers. Let $c$ be a center chosen from $A$ with $\mathcal{D}^{2}$-sampling with respect to $C$, where $\mathcal{D}\sim_{\delta}D$.
Then $\mathbf{E}[cost(A,C\cup\{c\})]\leq 72\cdot OPT(A)$.

###### Proof.

Let $D(a)$ denote the distance of the point $a$ from the nearest center in $C$ and let $\mathcal{D}(a)$ denote the estimated distance. We have $\mathcal{D}(a)\in D(a)\cdot(1\pm\delta)$. The following expression gives the expectation:

$$ $\displaystyle\sum_{a_{0}\in A}\frac{\mathcal{D}^{2}(a_{0})}{\sum_{a\in A} \mathcal{D}^{2}(a)}\cdot\sum_{a^{\prime}\in A}\min{\left(D^{2}(a^{\prime}),D^{ 2}(a^{\prime},a_{0})\right)}$ $$

Note that for all $a_{0},a\in A$, $D(a_{0})\leq D(a)+D(a,a_{0})$. This gives $\mathcal{D}(a_{0})\leq\frac{1+\delta}{1-\delta}\cdot\mathcal{D}(a)+(1+\delta)
\cdot D(a_{0},a)$, which further gives
$\mathcal{D}^{2}(a_{0})\leq 2\left(\frac{1+\delta}{1-\delta}\right)^{2}\cdot
\mathcal{D}^{2}(a)+2(1+\delta)^{2}\cdot D^{2}(a_{0},a)$ and $\mathcal{D}^{2}(a_{0})\leq\frac{2}{|A|}\left(\frac{1+\delta}{1-\delta}\right)^
{2}\cdot\sum_{a\in A}\mathcal{D}^{2}(a)+\frac{2}{|A|}(1+\delta)^{2}\cdot\sum_{
a\in A}D^{2}(a_{0},a)$.
We use this to obtain the following upper bound on the expectation $\mathbf{E}[cost(A,C\cup\{c\})]$:

$$ $\displaystyle\sum_{a_{0}\in A}\frac{\mathcal{D}^{2}(a_{0})}{\sum_{a\in A} \mathcal{D}^{2}(a)}\cdot\sum_{a^{\prime}\in A}\min{\left(D^{2}(a^{\prime}),D^{ 2}(a^{\prime},a_{0})\right)}$ $\displaystyle\leq$ $\displaystyle\sum_{a_{0}\in A}\frac{\left(\frac{2}{|A|}\left(\frac{1+\delta}{1 -\delta}\right)^{2}\sum_{a\in A}\mathcal{D}^{2}(a)\right)}{\sum_{a\in A} \mathcal{D}^{2}(a)}\cdot\sum_{a^{\prime}\in A}\min{\left(D^{2}(a^{\prime}),D^{ 2}(a^{\prime},a_{0})\right)}+$ $\displaystyle\sum_{a_{0}\in A}\frac{\left(\frac{2}{|A|}(1+\delta)^{2}\sum_{a \in A}D^{2}(a_{0},a)\right)}{\sum_{a\in A}\mathcal{D}^{2}(a)}\cdot\sum_{a^{ \prime}\in A}\min{\left(D^{2}(a^{\prime}),D^{2}(a^{\prime},a_{0})\right)}$ $\displaystyle\leq$ $\displaystyle\sum_{a_{0}\in A}\sum_{a^{\prime}\in A}\frac{2}{|A|}\left(\frac{1 +\delta}{1-\delta}\right)^{2}D^{2}(a^{\prime},a_{0})+\sum_{a_{0}\in A}\sum_{a \in A}\frac{2}{|A|}\left(\frac{1+\delta}{1-\delta}\right)^{2}D(a_{0},a)^{2}$ $\displaystyle=$ $\displaystyle\frac{4}{|A|}\left(\frac{1+\delta}{1-\delta}\right)^{2}\sum_{a_{0 }\in A}\sum_{a\in A}D^{2}(a_{0},a)$ $\displaystyle=$ $\displaystyle 8\left(\frac{1+\delta}{1-\delta}\right)^{2}OPT(A)$ $\displaystyle\leq$ $\displaystyle 72\cdot OPT(A).$ $$

This completes the proof of the lemma.
∎

We will use this lemma in the approximation scheme of [].
However, this lemma may be of independent interest as this gives a quantum pseudo approximation algorithm with a constant factor approximation that runs in time that is polylogarithmic in the data size and linear in $k$ and $d$. We will discuss this quantum algorithm in the next section.

###### Lemma 11 (Lemma 3.2 in [ 3 ] ) .

###### Definition 8 .

###### Lemma 12 .

###### Proof.

### C.2 Approximation scheme of [ 8 ]

A high-level description of the approximation scheme of [] was given in the introduction.
A more detailed pseudocode is given in Algorithm 4.

Figure: Algorithm 4 Algorithm of []. For our Quantum application, we will replaces $D$ with $\delta$-close $\mathcal{D}$.

In addition to the input instance $(V,k)$ and error parameter $\varepsilon$, the algorithm is also given a constant approximate solution $C_{init}$, which is used for $D^{2}$-sampling.
A pseudoapproximate solution $C_{init}$ is sufficient for the analysis in [].
The discussion from the previous subsection gives a robust algorithm that outputs a pseudoapproximate solution even under errors in distance estimates.
So, the input requirement of Algorithm 4 can be met.
Now, the main ingredient being $D^{2}$-sampling, we need to ensure that errors in distance estimate do not seriously impact the approximation analysis of Algorithm 4.
We state the main theorem of [] before giving the analogous statement for the modified algorithm where $D$ is replaced with $\mathcal{D}$ that is $\delta$-close to $D$.

###### Theorem 8 (Theorem 1 in [ 8 ] ) .

Let $0<\varepsilon\leq 1/2$ be the error parameter, $V\in\mathbb{R}^{N\times d}$ be the dataset, $k$ be a positive integer, and let $C_{init}$ be a constant (pseudo)approximate solution for dataset $V$.
Let $\mathcal{L}$ be the list returned by Algorithm 4 on input $(V,k,\varepsilon,C_{init})$ using the Euclidean distance function $D$. Then with probability at least $3/4$, $\mathcal{L}$ contains a center set $C$ such that $\Phi(V,C)\leq(1+\varepsilon)\cdot OPT$.
Moreover, $|\mathcal{L}|=\tilde{O}\left(2^{\tilde{O}(\frac{k}{\varepsilon})}\right)$ and the running time of the algorithm is $O(Nd|\mathcal{L}|)$.

We give the analogous theorem with access to the Euclidean distance function $D$ replaced with a function $\mathcal{D}$ that is $\delta$-close to $D$.

###### Theorem 9 .

Let $0<\varepsilon\leq\frac{1}{2}$ be the error parameter, $0<\delta<1/2$ be the closeness parameter, $V\in\mathbb{R}^{N\times d}$ be the dataset, $k$ be a positive integer, and let $C_{init}$ be a constant (pseudo)approximate solution for dataset $V$.
Let $\mathcal{L}$ be the list returned by Algorithm 4 on input $(V,k,\varepsilon,C_{init})$ using the distance function $\mathcal{D}$ that is $\delta$-close to the Euclidean distance function $D$. Then with probability at least $3/4$, $\mathcal{L}$ contains a center set $C$ such that $\Phi(V,C)\leq(1+\varepsilon)\cdot OPT$.
Moreover, $|\mathcal{L}|=2^{\tilde{O}(\frac{k}{\varepsilon})}$ and the running time of the algorithm is $O(Nd|\mathcal{L}|)$.

The proof of the above theorem closely follows the proof of Theorem 8 of [].
This is similar to the proof of Theorem 12 that we saw earlier, closely following the proof of Lemma 11.
The minor changes are related to approximate distance estimates using $\mathcal{D}$ instead of real estimates using $D$.
The statement of Theorem 9 is not surprising in this light.
Instead of repeating the entire proof of [], we point out the one change in their argument caused by using $\mathcal{D}$ instead of $D$ as the distance function.
The analysis of [] works by partitioning the points in any optimal cluster $X_{j}$ into those that are close to $C_{init}$ and those that are far.
For the far points, it is shown that when doing $D^{2}$-sampling, a far point will be sampled with probability at least $\gamma$ times the uniform sampling probability (see Lemma 21 in [], which is a full version of []).
It then argues that a reasonable size set of $D^{2}$-sampled points will contain a uniform sub-sample. A combination of the uniform sub-sample along with copies of points in $C_{init}$ gives a good center for this optimal cluster $X_{j}$.
Replacing $D$ with $\mathcal{D}$ decrease the value of $\gamma$ by a multiplicative factor of $\frac{(1-\delta)^{2}}{(1+\delta)^{2}}\geq(1-\delta)^{4}$.
This means that the number of points sampled should increase by a factor of $O(\frac{1}{(1-\delta)^{4}})$.
This means that the list size increases to $\tilde{O}\left(2^{\tilde{O}(\frac{k}{\varepsilon(1-\delta)})}\right)$.
Note that when $\delta\leq\frac{1}{2}$, the list size and running time retain the same form as that in [] (i.e., $|\mathcal{L}|=2^{\tilde{O}(\frac{k}{\varepsilon})}$ and time $O(Nd|\mathcal{L}|)$).
A detailed proof of Theorem 9 is provided in the next section of the Appendix.

###### Theorem 8 (Theorem 1 in [ 8 ] ) .

###### Theorem 9 .

### C.3 The quantum algorithm

The main idea of the quantum version of the approximation scheme of [] is to $D^{2}$-sample a set $S$ of $poly(k/\varepsilon)$ points with respect to a center set $C_{init}$ that gives constant approximation to the $k$-means objective.
We then try out all various subsets of possibilities for the $k$ center set out of $S$ and $C_{init}$ and then find the one with the least cost.
$C_{init}$ is itself found quantumly by iteratively $D^{2}$-sampling $2k$ points.
So, most of the quantum ideas developed for the quantum version of $k$-means++ can be reused in the quantum approximation scheme.
An additional effort is in to calculate the $k$-means cost of a list of $L=2^{\tilde{O}(\frac{k}{\varepsilon})}$ $k$-center sets and pick the one with the least cost.
Here, instead of adding the relevant steps only, we add all the steps of the quantum algorithm for better readability.

We will work under the assumption that the minimum distance between two data points is $1$, which can be achieved using scaling. This makes the aspect ratio $\zeta$ simply the maximum distance between two data points.
We will use $i$ for an index into the rows of the data matrix $V\in\mathbb{R}^{N\times d}$, and $j$ for an index into the rows of the center matrix $C\in\mathbb{R}^{2k\times d}$.
We would ideally like to design a quantum algorithm that performs the transformation:

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{{D(v_{i},c_{j})}}$ $$

Let us call the state on the right $\ket{\Psi_{ideal}}$.
This is an ideal quantum state for us since $\ket{\Psi_{ideal}}$ helps to perform $D^{2}$-sampling and to find the $k$-means cost of clustering, which are the main components of the approximation scheme of [] that we intend to use.
One caveat is that we will only be able to perform the following transformation (instead of the abovementioned transformation)

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{\psi_{i,j}},$ $$

where $\ket{\psi_{i,j}}$ is an approximation for $\ket{\mathcal{D}(v_{i},c_{j})}$ in a sense that we will make precise below.
We will use $\ket{\Psi_{real}}$ to denote the state $\ket{i}\ket{j}\ket{\psi_{i,j}}$.
This state is prepared using tools such as swap test followed by coherent amplitude estimation, and median estimation.
Since these tools and techniques are known from previous works [, , ], we summarise the discussion (see Section 4.1 and 4.2 in []) in the following lemma.

###### Lemma 13 ( [ 20 ] and [ 31 ] ) .

Assume for a data matrix $V\in\mathbb{R}^{N\times d}$ and a center set matrix $C\in\mathbb{R}^{t\times d}$ that the following unitaries: (i) $\ket{i}\ket{0}\rightarrow\ket{i}\ket{v_{i}}$, (ii) $\ket{j}\ket{0}\rightarrow\ket{j}\ket{c_{j}}$ can be performed in time $T$ and the norms of the vectors are known. For any $\Delta>0$, there is a quantum algorithm that in time $O\left(\frac{T\log{\frac{1}{\Delta}}}{\varepsilon}\right)$ computes:

$$ $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{\psi_{i,j}},$ $$

where $\ket{\psi_{i,j}}$ satisfies the following two conditions for every $i\in[N]$ and $j\in[t]$:

- (i)
$\norm{\ket{\psi_{i,j}}-\ket{0^{\otimes\ell}}\ket{{\tilde{D}(v_{i},c_{j})}}}
\leq\sqrt{2\Delta}$, and
- (ii)
For every $i,j$, $\tilde{D}(v_{i},c_{j})\in(1\pm\varepsilon)\cdot D(v_{i},c_{j})$.

In the subsequent discussions, we will use $T$ as the time to access the QRAM data structure [], i.e., for the transitions $\ket{i}\ket{0}\rightarrow\ket{i}\ket{v_{i}}$ and $\ket{j}\ket{0}\rightarrow\ket{j}\ket{c_{j}}$ as given in the above lemma.
This is known to be $T=O(\log^{2}{(Nd)})$. Moreover, the time to update each entry in this data structure is also $T=O(\log^{2}{(Nd)})$.
This is the logarithmic factor that is hidden in the $\tilde{O}$ notation.
In the following subsections, we discuss the utilities of $\ket{\Psi_{real}}$ for the various components of the approximation scheme of [].
During these discussions, it will be easier to see the utility first with the ideal state $\ket{\Psi_{ideal}}$ before the real state $\ket{\Psi_{real}}$ that can actually be prepared.
We will see how $\ket{\Psi_{real}}$ is sufficient within a reasonable error bound.

###### Lemma 13 ( [ 20 ] and [ 31 ] ) .

### C.4 Finding distance to closest center

Let us see how we can estimate the distance of any point to its closest center in a center set $C\coloneqq C_{init}$ with $t\coloneqq 2k$ centers.
We can use the transformation $\ket{i}\ket{j}\ket{0}\rightarrow\ket{i}\ket{j}\ket{{D(v_{i},c_{j})}}$ to prepare the following state for any $i$:

$$ $\ket{i}\ket{D(v_{i},c_{1})}\ket{D(v_{i},c_{2})}...\ket{D(v_{i},c_{t})}$ $$

We can then iteratively compare and swap pairs of registers to prepare the state $\ket{i}\ket{\min_{j\in[t]}D(v_{i},c_{j})}$.
If we apply the same procedure to $\ket{i}\ket{\psi_{i,1}}...\ket{\psi_{i,t}}$, then with probability at least $(1-2\Delta)^{t}$, the resulting state will be $\ket{i}\ket{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})}}$. So, the contents of the second register will be an estimate of the distance of the $i^{th}$ point to its closest center in the center set $C$. This further means that the following state can be prepared with probability at least $(1-2\Delta)^{Nt}$:(^12^1212The state prepared is actually $\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\left(\alpha\ket{\min_{j\in[t]}{\tilde{ D}(v_{i},c_{j})}}+\beta\ket{G}\right)$ with $|\alpha|^{2}\geq(1-2\Delta)^{Nt}$. However, instead of working with this state, subsequent discussions become much simpler if we assume that $\ket{\Psi_{C}}$ is prepared with probability $|\alpha|^{2}$.)

$$ $\ket{\Psi_{C}}\equiv\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\ket{\min_{j\in[t]} {\tilde{D}(v_{i},c_{j})}}.$ $$

This quantum state can be used to find the approximate clustering cost of the center set $C$, which we discuss in the following subsection.
However, before we do that, let us summarise the main ideas of this subsection in the following lemma.

###### Lemma 14 .

There is a quantum algorithm that, with probability at least $(1-2\Delta)^{Nt}$, prepares the quantum state $\ket{\Psi_{C}}$ in time $O\left(\frac{Tt\log{\frac{1}{\Delta}}}{\varepsilon}\right)$.

###### Lemma 14 .

### C.5 Computing cost of clustering

Suppose we want to compute the $k$-means cost, $\Phi(V,C)\equiv\sum_{i=1}^{N}\min_{j\in[t]}{D^{2}(v_{i},c_{j})}$, of the clustering given by a $t$ center set $C$.
We can prepare $m$ copies of the state $\ket{\Psi_{C}}$ and then estimate the cost of clustering by measuring $m$ copies of this quantum state and summing the squares of the second registers.
If $m$ is sufficiently large, we obtain a close estimate of $\Phi(V,C)$.
To show this formally, we will use the following Hoeffding tail inequality.

###### Theorem 10 (Hoeffding bound) .

Let $X_{1},...,X_{m}$ be independent, bounded random variables such that $X_{i}\in[a,b]$. Let $S_{m}=X_{1}+...+X_{m}$ . Then for any $\theta>0$, we have:

$$ $\mathbf{Pr}[|S_{m}-\mathbf{E}[S_{m}]|\geq\theta]\leq 2\cdot e^{\frac{-2\theta^ {2}}{m(b-a)^{2}}}.$ $$

Let $X_{1},...,X_{m}$ denote the square of the measured value of the second register in $\ket{\Psi_{C}}$.
These are random values in the range $[1,\zeta^{2}]$, where $\zeta=\max_{i,j}\tilde{D}(v_{i},c_{j})\in(1\pm\varepsilon)\cdot\max_{i,j}{D(v_
{i},v_{j})}$.
First, we note the expectation of these random variables equals $\frac{\tilde{\Phi}(V,C)}{N}$, where $\tilde{\Phi}(V,C)\equiv\sum_{i=1}^{N}\min_{j\in[t]}\tilde{D}(v_{i},c_{j})^{2}
\in(1\pm\varepsilon)^{2}\cdot\Phi(V,C)$.
We define the variable $S_{m}=X_{1}+X_{2}+...+X_{m}$ and apply the Hoeffding bound on these bounded random variables to get a concentration result that can then be used.

###### Lemma 15 .

Let $\alpha_{m}=S_{m}\cdot\frac{N}{m}$ and $L>0$. If $m=\frac{\zeta^{2}\ln{(10L)}}{\varepsilon^{2}}$, then we have:

$$ $\mathbf{Pr}[\alpha_{m}\in(1\pm\varepsilon)\cdot\tilde{\Phi}(V,C)]\geq 1-\frac{ 1}{5L}.$ $$

###### Proof.

We know that $\mathbf{E}[S_{m}]=\frac{m}{N}\cdot\tilde{\Phi}(V,C)$
From the Hoeffding tail inequality, we get the following:

$$ $\displaystyle\mathbf{Pr}[|S_{m}-\mathbf{E}[S_{m}]|\geq\varepsilon\cdot\mathbf{ E}[S_{m}]]$ $\displaystyle\leq$ $\displaystyle 2\cdot e^{\frac{-2\varepsilon^{2}\mathbf{E}[S_{m}]^{2}}{m\zeta^{ 2}}}$ $\displaystyle=$ $\displaystyle 2\cdot e^{\frac{-2\varepsilon^{2}m}{\zeta^{2}}\cdot\left(\frac{ \tilde{\Phi}(V,C)}{N}\right)^{2}}$ $\displaystyle\leq$ $\displaystyle 2\cdot e^{-\ln{(10L)}}\leq\frac{1}{5L}.$ $$

This implies that:

$$ $\mathbf{Pr}[|\alpha_{m}-\tilde{\Phi}(V,C)|\geq\varepsilon\cdot\tilde{\Phi}(V,C )]\leq\frac{1}{5L}.$ $$

This completes the proof of the lemma.
∎

So, conditioned on having $m$ copies of the state $\ket{\Psi_{C}}$, we get an estimate of the clustering cost within a relative error of $(1\pm\varepsilon)^{3}$ with probability at least $(1-\frac{1}{5L})$. Removing the conditioning, we get the same with probability at least $(1-2\Delta)^{Nkm}\cdot(1-\frac{1}{5L})$.
We want to use the above cost estimation technique to calculate the cost for a list of center sets $C_{1},...,C_{L}$, and then pick the center set from the list with the least cost. We must apply the union bound appropriately to do this with high probability. We summarise these results in the following lemma.
Let us first set some of the parameters with values that we will use to implement the approximation scheme of [].

- •
$L$ denotes the size of the list of $k$-center sets we will iterate over to find the one with the least cost. This quantity is bounded as $L=\left(\frac{k}{\varepsilon}\right)^{O(\frac{k}{\varepsilon})}$.
- •
$m$ is the number of copies of the state $\ket{\Psi_{C}}$ made to estimate the cost of the center set $C$. This, as given is Lemma 15 is $m=\frac{\zeta^{2}\ln{(10L)}}{\varepsilon^{2}}$, where $\zeta=(1+\varepsilon)\cdot\max_{i,j}{D(v_{i},c_{j})}$.

###### Lemma 16 .

Let $L=\left(\frac{k}{\varepsilon}\right)^{O(\frac{k}{\varepsilon})}$, $m=\frac{\zeta^{2}\ln{(10L)}}{\varepsilon^{2}}$, and $\Delta=O\left(\frac{1}{NkmL}\right)$.
Given a point set $V$ and a list of center sets $C_{1},...,C_{L}$ in the QRAM model, there is a quantum algorithm that runs in time
$\tilde{O}\left(2^{\tilde{O}(\frac{k}{\varepsilon})}T\zeta^{2}\right)$
and outputs an index $l$ such that $\Phi(V,C_{l})\leq(1+\varepsilon)^{3}\min_{j\in L}\Phi(V,C_{j})$ with probbaility at least $\frac{3}{5}$.

###### Proof.

The algorithm estimates the cost of $C_{1},...,C_{L}$ using $m$ copies each of $\ket{\Psi_{C_{1}}},...,\ket{\Psi_{C_{L}}}$ and picks the index with the minimum value in time $O\left(\frac{TkmL\log{\frac{1}{\Delta}}}{\varepsilon}\right)$.
Plugging the values of $L,m$, and $\Delta$ we get the running time stated in the lemma.

Let us bound the error probability of this procedure.
By Lemma 3, the probability that we do not have the correct $m$ copies each of $\ket{\Psi_{C_{1}}},...,\ket{\Psi_{C_{L}}}$ is bounded by $1-(1-2\Delta)^{NkmL}$.
Conditioned on having $\ket{\Psi_{C_{1}}},...,\ket{\Psi_{C_{L}}}$, the probability that there exists an index $j\in[L]$, where the estimate is off by more than a $(1\pm\varepsilon)^{3}$ factor is upper bounded by $\frac{1}{5}$ by the union bound. So, the probability that the algorithm will find an index $l$ such that $\Phi(V,C_{l})>(1+\varepsilon)^{3}\min_{j\in[L]}{\Phi(V,C_{j})}$ is upper bounded by $1-(1-2\Delta)^{NkmL}+\frac{1}{5}$. This probability is at most $\frac{2}{5}$ since $\Delta=O(\frac{1}{NkmL})$. This completes the proof of the lemma.
∎

###### Theorem 10 (Hoeffding bound) .

###### Lemma 15 .

###### Proof.

###### Lemma 16 .

###### Proof.

### C.6 D 2 superscript 𝐷 2 D^{2} italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT -sampling

$D^{2}$-sampling from the point set $V$ with respect to a center set $C\in\mathbb{R}^{t\times d}$ with $t$ centers, samples $v_{i}$ with probability proportional to $\min_{j\in[t]}{D^{2}(v_{i},c_{j})}$.
Let us see if we can use our state $\ket{\Psi_{C}}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\ket{\min_{j\in[t]}{
\tilde{D}(v_{i},c_{j})}}$ is useful to perform this sampling.
If we can pull out the value of the second register as the amplitude, then the measurement will give us close to $D^{2}$-sampling.
This is possible since we have an estimate of the clustering cost from the previous subsection.
We can use controlled rotations on an ancilla qubit to prepare the state:

$$ $\ket{\Psi_{sample}}\equiv\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\ket{i}\left(\beta_{i }\ket{0}+\sqrt{1-\beta_{i}}\ket{1}\right),$ $$

where $\beta_{i}=\frac{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})}}{\sqrt{2\cdot\tilde{
\Phi}(V,C)}}$.
So, the probability of measurement of $(i,0)$ is $\frac{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})^{2}}}{2\cdot\tilde{\Phi}(V,C)}$. Since we do rejection sampling (ignoring $(.,1)$’s that are sampled with probability $\frac{1}{2}$), we end up sampling with a distribution where the probability of sampling $i$ is $\frac{\min_{j\in[t]}{\tilde{D}(v_{i},c_{j})^{2}}}{\tilde{\Phi}(V,C)}\in(1\pm
\varepsilon)^{4}\cdot\frac{\min_{j\in[t]}{D(v_{i},c_{j})^{2}}}{\Phi(V,C)}$.
This means that points get sampled with a probability close to the actual $D^{2}$-sampling probability. As we have mentioned earlier, this is sufficient for the approximation guarantees of [] to hold.
We summarise the observations of this section in the next lemma.

###### Lemma 17 .

Given a dataset $V\in\mathbb{R}^{N\times d}$ and a center set $C\in\mathbb{R}^{t\times d}$ in the QRAM model, there is a quantum algorithm that runs in time $O\left(\frac{TtmS\log{\frac{1}{\Delta}}}{\delta}\right)$ and with probability at least $(1-2\Delta)^{NtmS}$ outputs $S$ independent samples with
$\mathcal{D}^{2}$ distribution such that the distance function $\mathcal{D}$ is $\delta$-close to $D$.

###### Proof.

The proof follows from Lemma 14 and the preceding discussion.
∎

The above lemma says that for $\Delta=\tilde{O}(\frac{1}{NkmS})$, we obtain the required samples with high probability.
We can now give proof of Theorem 4 assembling the quantum tools of this section. We restate the theorem below for ease of reading.

###### Theorem (Restatement of Theorem 4 ) .

Let $0<\varepsilon<1/2$ be the error parameter.
There is a quantum algorithm that, when given QRAM data structure access to a dataset $V\in\mathbb{R}^{N\times d}$, runs in time $\tilde{O}\left(2^{\tilde{O}(\frac{k}{\varepsilon})}d\zeta^{2}\right)$ and outputs a $k$ center set $C\in\mathbb{R}^{k\times d}$ such that with high probability $\Phi(V,C)\leq(1+\varepsilon)\cdot OPT$. Here, $\zeta$ is the aspect ratio, i.e., the ratio of the maximum to the minimum distance between two given points in $V$.(^13^1313The $\tilde{O}$ notation hides logarithmic factors in $N$. The $\tilde{O}$ in the exponent hides logarithmic factors in $k$ and $1/\varepsilon$.)

###### Proof of Theorem 4 .

The first requirement for executing the algorithm of [] is a constant pseudo approximation algorithm using which we obtain the initial center set $C_{init}$.
By Lemma 12, we know that $2k$ points sampled using $\mathcal{D}^{2}$-sampling gives such a center set. From Lemma 17, this can be done quantumly in time $\tilde{O}(\frac{k^{2}\zeta^{2}}{\varepsilon^{2}})$.
The algorithm of [] has an outer repeat loop for probability amplification.
Within the outer loop, $poly(\frac{k}{\varepsilon})$ points are $\mathcal{D}^{2}$-sampled with respect to the center set $C_{init}$ (line 6).
This can again be done quantumly using Lemma 17 in time $\tilde{O}(\zeta^{2}(k/\varepsilon)^{O(1)})$.
We can then classically process the point set $M$ (see line 7 in Algorithm 4) and create the QRAM data structure for the list $C_{1},...,C_{L}$ of $k$-center sets that correspond to all possible disjoint subsets of $M$ (see line 8 in Algorithm 4). This takes time $\tilde{O}(Lkd)$, where $L=\left(\frac{k}{\varepsilon}\right)^{O(\frac{k}{\varepsilon})}$.
Theorem 9 shows that at least one center set in the list gives $(1+\varepsilon)$-approximation.
We use this fact in conjunction with the result of Lemma 16 to get that the underlying quantum algorithm runs in time $\tilde{O}(2^{\tilde{O}(\frac{k}{\varepsilon})}d\zeta^{2})$ and with high probability outputs a center set $C^{\prime}$ such that $\Phi(V,C^{\prime})\leq(1+\varepsilon)^{4}\cdot OPT$.(^14^1414We needed $(1+\varepsilon)$, but got $(1+\varepsilon)^{4}$ instead. However, this can be handled with $\varepsilon^{\prime}=\varepsilon/5$.)
∎

###### Lemma 17 .

###### Proof.

###### Theorem (Restatement of Theorem 4 ) .

###### Proof of Theorem 4 .

## Appendix D A Robust Approximation Scheme for k 𝑘 k italic_k -means (proof of Theorem 9 )

This section gives the proof of Theorem 9.
We restate the theorem and the algorithm for readability.
We state the algorithm of Bhattacharya et al.  [] replacing $D^{2}$-sampling with $\mathcal{D}^{2}$-sampling.

Figure: Algorithm 5 Algorithm of [] with $D$ replaced with $\delta$-close $\mathcal{D}$

###### Theorem (Restatement of Theorem 9 ) .

Let $0<\varepsilon\leq\frac{1}{2}$ be the error parameter, $0<\delta<1/2$ be the closeness parameter, $V\in\mathbb{R}^{N\times d}$ be the dataset, $k$ be a positive integer, and let $C_{init}$ be a constant (pseudo)approximate solution for dataset $V$.
Let $\mathcal{L}$ be the list returned by Algorithm 5 on input $(V,k,\varepsilon,C_{init})$ using the distance function $\mathcal{D}$ that is $\delta$-close to the Euclidean distance function $D$. Then with probability at least $3/4$, $\mathcal{L}$ contains a center set $C$ such that $\Phi(V,C)\leq(1+\varepsilon)\cdot OPT$.
Moreover, $|\mathcal{L}|=2^{\tilde{O}(\frac{k}{\varepsilon})}$ and the running time of the algorithm is $O(Nd|\mathcal{L}|)$.

Note that what this theorem essentially says is that if the error in the distance estimates is bounded within a multiplicative factor of $(1\pm\delta)$ with $\delta\leq 1/2$, then the guarantees of the algorithm of Bhattacharya et al.  [] do not change except for some constant factors (that get absorbed in the $\tilde{O}$ of the exponent).

We need to define a few quantities that will be used to prove the above theorem.
Let $\{X_{1},...,X_{k}\}$ be the optimal clusters.
It is well-known that the center that optimizes the $1$-means cost for any pointset is the centroid of the pointset.
We denote the centroid of any pointset $Y\subset\mathbb{R}^{d}$ with $\mu(Y)\equiv\frac{\sum_{y\in Y}y}{|Y|}$.
This follows from the following well-known fact:

###### Fact 1 .

For any $Y\subset\mathbb{R}^{d}$ and any $c\in\mathbb{R}^{d}$, we have $\sum_{y\in Y}\norm{y-c}^{2}=\sum_{y\in Y}\norm{y-\mu(Y)}^{2}+|Y|\cdot\norm{c-
\mu(Y)}^{2}$.

We have used the norm notation instead of $D(.,.)$ for the Euclidean distance function. This is to easily distinguish between the usage of $D$ and $\mathcal{D}$.
We define the optimal $1$-means cost for pointset $X_{i}$ as $\Delta(X_{i})\equiv\sum_{x\in X_{i}}\norm{x-\mu(X_{i})}^{2}$.
We denote the optimal $k$-means cost using $OPT\equiv\sum_{i=1}^{k}\Delta(X_{i})$.
We will use the following sampling lemma from Inaba et al.  [], which says that the centroid of a small set of uniformly sampled points from the dataset $Y$ is a good center with respect to the $1$-means cost for dataset Y.

###### Lemma 18 ( [ 19 ] ) .

Let $S$ be a set of points obtained by independently sampling $M$ points with replacement uniformly at random from a point set $Y\subset\mathbb{R}^{d}$. Then for any $\delta>0$,

$$ $\mathbf{Pr}\left[\Phi(\mu(S),Y)\leq\left(1+\frac{1}{\delta M}\right)\cdot \Delta(Y)\right]\geq(1-\delta).$ $$

We will also use the following approximate triangle inequality that holds for squared distances.

###### Fact 2 (Approximate triangle inequality) .

For any $x,y,z\in\mathbb{R}^{d}$, we have $\norm{x-z}^{2}\leq 2\cdot(\norm{x-y}^{2}+\norm{y-z}^{2})$.

The proof of Theorem 9 follows from the following theorem, which we will prove in the remaining section.

###### Theorem 11 .

Let $0<\varepsilon,\delta\leq 1/2$.
Let $\mathcal{L}$ denote the list returned by Algorithm 5 on input ($V,k,\varepsilon,C_{init}$).
Then with probability at least $3/4$, $\mathcal{L}$ contains a center set $C$ such that

$$ $\Phi\left(C,\{X_{1},...,X_{k}\}\right)\leq\left(1+\frac{\varepsilon}{2}\right) \cdot\sum_{i=1}^{k}\Delta(X_{i})+\frac{\varepsilon}{2}\cdot OPT.$ $$

Moreover, $|\mathcal{L}|=2^{\tilde{O}(\frac{k}{\varepsilon})}$ and the running time of the algorithm is $O(nd|\mathcal{L}|)$.

Let $C_{init}$ be an $(\alpha,\beta)$-approximate solution to the $k$-means problem on the dataset $V$. That is:

$$ $\Phi(V,C_{init})\leq\alpha\cdot OPT\quad\textrm{and}\quad|C_{init}|\leq\beta \cdot k.$ $$

Let us now analyze the algorithm. Note that the outer repetition (lines 5-12) operation is to amplify the probability of the list containing a good $k$ center set.
We will show that the probability of finding a good $k$ center set in one iteration is at least $(3/4)^{k}$, and the theorem will follow from standard probability analysis. So, we shall focus on only one iteration of the outer loop. Consider the multi-set $M$ obtained on line (7) of Algorithm 5.
We will show that with probability at least $(3/4)^{k}$, there are disjoint multi-subsets $T_{1},...,T_{k}$ each of size $\tau$ such that for every $i=1,...,k$:

$$ $\Phi(\mu(T_{i}),X_{i})\leq\left(1+\frac{\varepsilon}{2}\right)\cdot\Delta(X_{i })+\frac{\varepsilon}{2k}\cdot OPT.$ (5) $$

Since we try all possible subsets of $M$ in line (8), we will get the desired result.
More precisely, we will argue in the following manner:
consider the multi-set $C^{\prime}=\{\tau k\textrm{ copies of each element in $C_{init}$}\}$.
We can interpret $C^{\prime}$ as a union of multi-sets $C_{1}^{\prime},C_{2}^{\prime},...,C_{k}^{\prime}$, where $C_{i}^{\prime}=\{\tau\textrm{ copies of each element in $C_{init}$}\}$.
Also, since $M$ consists of $\rho k$ independently sampled points, we can interpret $M$ as a union of multi-sets $M_{1}^{\prime},M_{2}^{\prime},...,M_{k}^{\prime}$ where $M_{1}^{\prime}$ is the first $\rho$ points sampled, $M_{2}^{\prime}$ is the second $\rho$ points and so on.
For all $i=1,...,k$, let $M_{i}=C_{i}^{\prime}\cup(M_{i}^{\prime}\cap X_{i})$.(^15^1515$M_{i}^{\prime}\cap X_{i}$ in this case, denotes those points in the multi-set $M_{i}^{\prime}$ that belongs to $X_{i}$.)
We will show that for every $i\in\{1,...,k\}$, with probability at least $(3/4)$, $M_{i}$ contains a subset $T_{i}$ of size $\tau$ that satisfies Eqn. (5).
Note that $T_{i}$’s being disjoint follows from the definition of $M_{i}$.
It will be sufficient to prove the following lemma.

###### Lemma 19 .

Consider the sets $M_{1},...,M_{k}$ as defined above. For any $i\in\{1,...,k\}$,

$$ $\mathbf{Pr}\left[\exists T_{i}\subseteq M_{i}\textrm{ s.t. }|T_{i}|=\tau \textrm{ and }\left(\Phi(\mu(T_{i}),X_{i})\leq\left(1+\frac{\varepsilon}{2} \right)\cdot\Delta(X_{i})+\frac{\varepsilon}{2k}OPT\right)\right]\geq\frac{3}{ 4}.$ $$

We prove the above lemma in the remaining discussion.
We do a case analysis that is based on whether $\frac{\Phi(C_{init},X_{i})}{\Phi(C_{init},V)}$ is large or small for a particular $i\in\{1,...,k\}$.

- •
Case-I $\left(\Phi(C_{init},X_{i})\leq\frac{\varepsilon}{6\alpha k}\cdot\Phi(C_{init},
V)\right)$:
Here, we will show that there is a subset $T_{i}\subseteq C_{i}^{\prime}\subseteq M_{i}$ that satisfies Eqn. (5).
- •
Case-II $\left(\Phi(C_{init},X_{i})>\frac{\varepsilon}{6\alpha k}\cdot\Phi(C_{init},V)\right)$:
Here we will show that $M_{i}$ contains a subset $T_{i}$ such that $\Phi(\mu(T_{i}),X_{i})\leq\left(1+\frac{\varepsilon}{2}\right)\cdot\Delta(X_{i})$ and hence $T_{i}$ also satisfies Eqn. (5).

We will discuss these two cases next.

###### Theorem (Restatement of Theorem 9 ) .

###### Fact 1 .

###### Lemma 18 ( [ 19 ] ) .

###### Fact 2 (Approximate triangle inequality) .

###### Theorem 11 .

###### Lemma 19 .

#### Case-I: ( Φ ⁢ ( C i ⁢ n ⁢ i ⁢ t , X i ) ≤ ε 6 ⁢ α ⁢ k ⋅ Φ ⁢ ( C i ⁢ n ⁢ i ⁢ t , V ) ) Φ subscript 𝐶 𝑖 𝑛 𝑖 𝑡 subscript 𝑋 𝑖 ⋅ 𝜀 6 𝛼 𝑘 Φ subscript 𝐶 𝑖 𝑛 𝑖 𝑡 𝑉 \left(\Phi(C_{init},X_{i})\leq\frac{\varepsilon}{6\alpha k}\cdot\Phi(C_{init},%
V)\right) ( roman_Φ ( italic_C start_POSTSUBSCRIPT italic_i italic_n italic_i italic_t end_POSTSUBSCRIPT , italic_X start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ≤ divide start_ARG italic_ε end_ARG start_ARG 6 italic_α italic_k end_ARG ⋅ roman_Φ ( italic_C start_POSTSUBSCRIPT italic_i italic_n italic_i italic_t end_POSTSUBSCRIPT , italic_V ) )

First, observe that:

$$ $\Phi(C_{init},X_{i})\leq\frac{\varepsilon}{6k}\cdot OPT\qquad\textrm{(since $C _{init}$ is an $\alpha$-approximate solution)}$ (6) $$

For any point $x\in V$, let $c(x)$ denote the center in the set $C$ that is closest to $x$.
That is, $c(x)=\arg\min_{c\in C_{init}}{\norm{c-x}}$.
Given this definition, note that:

$$ $\sum_{x\in X_{i}}\norm{x-c(x)}^{2}=\Phi(C_{init},X_{i})$ (7) $$

We define the multi-set $X_{i}^{\prime}=\{c(x):x\in X_{i}\}$.
Let $m$ and $m^{\prime}$ denote the centroids of the point sets $X_{i}$ and $X_{i}^{\prime}$, respectively.
So, we have $\Delta(X_{i})=\Phi(m,X_{i})$ and $\Delta(X_{i}^{\prime})=\Phi(m^{\prime},X_{i}^{\prime})$.
We will show that $\Delta(X_{i})\approx\Delta(X_{i}^{\prime})$.
First, we bound the distance between $m$ and $m^{\prime}$.

###### Lemma 20 .

$\norm{m-m^{\prime}}^{2}\leq\frac{\Phi(C_{init},X_{i})}{|X_{i}|}$.

###### Proof.

We have:

$$ $\displaystyle\norm{m-m^{\prime}}^{2}=\frac{\norm{\sum_{x\in X_{i}}(x-c(x))}^{2 }}{|X_{i}|^{2}}\leq\frac{\sum_{x\in X_{j}}\norm{(x-c(x))}^{2}}{|X_{j}|}=\frac{ \Phi(C_{init},X_{j})}{|X_{j}|}.$ $$

where the second to last inequality follows from Cauchy-Schwartz.
∎

###### Lemma 21 .

$\Delta(X_{i}^{\prime})\leq 2\cdot\Phi(C_{init},X_{i})+2\cdot\Delta(X_{i})$.

###### Proof.

We have:

$$ $\displaystyle\Delta(X_{i}^{\prime})$ $\displaystyle=$ $\displaystyle\sum_{x\in X_{i}}\norm{c(x)-m^{\prime}}^{2}\leq\sum_{x\in X_{i}} \norm{c(x)-m}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact\ \ref{fact:2})}}}{{\leq}}$ $\displaystyle 2\cdot\sum_{x\in X_{i}}(\norm{c(x)-x}^{2}+\norm{x-m}^{2})=2\cdot \Phi(C_{init},X_{i})+2\cdot\Delta(X_{i})$ $$

This completes the proof of the lemma.
∎

We now show that a good center for $X_{i}^{\prime}$ will also be a good center for $X_{i}$.

###### Lemma 22 .

Let $m^{\prime\prime}$ be a point such that $\Phi(m^{\prime\prime},X_{i}^{\prime})\leq(1+\frac{\varepsilon}{8})\cdot\Delta(
X_{i}^{\prime})$.
Then $\Phi(m^{\prime\prime},X_{i})\leq(1+\frac{\varepsilon}{2})\cdot\Delta(X_{i})+
\frac{\varepsilon}{2k}\cdot OPT$.

###### Proof.

We have:

$$ $\displaystyle\Phi(m^{\prime\prime},X_{i})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:1})}}}{{=}}$ $\displaystyle\sum_{x\in X_{i}}\norm{x-m}^{2}+|X_{i}|\cdot\norm{m-m^{\prime \prime}}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:2})}}}{{\leq}}$ $\displaystyle\Delta(X_{i})+2|X_{i}|\cdot(\norm{m-m^{\prime}}^{2}+\norm{m^{ \prime}-m^{\prime\prime}}^{2})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma~{}\ref{lemma:case1-1})}}}{{ \leq}}$ $\displaystyle\Delta(X_{i})+2\cdot\Phi(C_{init},X_{i})+2|X_{i}|\cdot\norm{m^{ \prime}-m^{\prime\prime}}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:1})}}}{{\leq}}$ $\displaystyle\Delta(X_{i})+2\cdot\Phi(C_{init},X_{i})+2(\Phi(m^{\prime\prime}, X_{i}^{\prime})-\Delta(X_{i}^{\prime}))$ $\displaystyle\leq$ $\displaystyle\Delta(X_{i})+2\cdot\Phi(C,X_{i})+\frac{\varepsilon}{4}\cdot \Delta(X_{i}^{\prime})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma~{}\ref{lemma:case1-2})}}}{{ \leq}}$ $\displaystyle\Delta(X_{i})+2\cdot\Phi(C_{init},X_{i})+\frac{\varepsilon}{2} \cdot(\Phi(C_{init},X_{i})+\Delta(X_{i}))$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Eqn.~{}\ref{eqn:case1-1})}}}{{\leq}}$ $\displaystyle\left(1+\frac{\varepsilon}{2}\right)\cdot\Delta(X_{i})+\frac{ \varepsilon}{2k}\cdot OPT.$ $$

This completes the proof of the lemma.
∎

We know from Lemma 18 that there exists a (multi) subset of $X_{i}^{\prime}$ of size $\tau$ such that the mean of these points satisfies the condition of the lemma above.
Since $C_{i}^{\prime}$ contains at least $\tau$ copies of every element of $C_{init}$, there is guaranteed to be a subset $T_{i}\subseteq C_{i}^{\prime}$ that satisfies Eqn. (5).
So, for any index $i\in\{1,...,k\}$ such that $\frac{\Phi(C_{init},X_{i})}{\Phi(C_{init},V)}\leq\frac{\varepsilon}{6\alpha k}$, $M_{i}$ has a good subset $T_{i}$ with probability $1$.

###### Lemma 20 .

###### Proof.

###### Lemma 21 .

###### Proof.

###### Lemma 22 .

###### Proof.

#### Case-II: ( Φ ⁢ ( C i ⁢ n ⁢ i ⁢ t , X i ) > ε 6 ⁢ α ⁢ k ⋅ Φ ⁢ ( C i ⁢ n ⁢ i ⁢ t , V ) ) Φ subscript 𝐶 𝑖 𝑛 𝑖 𝑡 subscript 𝑋 𝑖 ⋅ 𝜀 6 𝛼 𝑘 Φ subscript 𝐶 𝑖 𝑛 𝑖 𝑡 𝑉 \left(\Phi(C_{init},X_{i})>\frac{\varepsilon}{6\alpha k}\cdot\Phi(C_{init},V)\right) ( roman_Φ ( italic_C start_POSTSUBSCRIPT italic_i italic_n italic_i italic_t end_POSTSUBSCRIPT , italic_X start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) > divide start_ARG italic_ε end_ARG start_ARG 6 italic_α italic_k end_ARG ⋅ roman_Φ ( italic_C start_POSTSUBSCRIPT italic_i italic_n italic_i italic_t end_POSTSUBSCRIPT , italic_V ) )

If we can show that a $\mathcal{D}^{2}$-sampled set with respect to center set $C_{init}$ has a subset $S$ that may be considered a uniform sample from $X_{i}$, then we can use Lemma 18 to argue that $M_{i}$ has a subset $T_{i}$ such that $\mu(T_{i})$ is a good center for $X_{i}$.
Note that since $\frac{\Phi(C_{init},X_{j})}{\Phi(C_{init},V)}>\frac{\varepsilon}{6\alpha k}$, and $D$ is $\delta$-close to $\mathcal{D}$, we can argue that if we $\mathcal{D}^{2}$-sample $poly(\frac{k}{\varepsilon})$ elements, then we will get a good representation from $X_{i}$.
However, some of the points from $X_{i}$ may be very close to one of the centers in $C_{init}$ and hence will have a very small chance of being $\mathcal{D}^{2}$-sampled.
In such a case, no subset $S$ of a $\mathcal{D}^{2}$-sampled set will behave like a uniform sample from $X_{i}$.
So, we need to argue more carefully taking into consideration the fact that there may be points in $X_{i}$ for which the chance of being $\mathcal{D}^{2}$-sampled may be very small.
Here is the high-level argument that we will build:

- •
Consider the set $X_{i}^{\prime}$ which is same as $X_{i}$ except that points in $X_{i}$ that are very close to $C_{init}$ have been “collapsed” to their closest center in $C_{init}$.
- •
Argue that a good center for the set $X_{i}^{\prime}$ is a good center for $X_{i}$.
- •
Show that a convex combination of copies of centers in $C_{init}$ (i.e., $C_{i}^{\prime}$) and $\mathcal{D}^{2}$-sampled points from $X_{i}$ gives a good center for the set $X_{i}^{\prime}$.

The closeness of point in $X_{i}$ to points in $C_{init}$ is quantified using radius $R$ that is defined by the equation:

$$ $R^{2}\stackrel{{\scriptstyle defn.}}{{=}}\frac{\varepsilon^{2}}{41}\cdot\frac{ \Phi(C_{init},X_{i})}{|X_{i}|}.$ (8) $$

Let $X_{i}^{near}$ be the points in $X_{i}$ that are within a distance of $R$ from a point in set $C_{init}$ and $X_{i}^{far}=X_{i}\setminus X_{i}^{near}$. That is,
$X_{i}^{near}=\{x\in X_{i}:\min_{c\in C_{init}}{\norm{x-c}}\leq R\}$ and $X_{i}^{far}=X_{i}\setminus X_{i}^{near}.$
Using these, we define the multi-set $X_{i}^{\prime}$ as:

$$ $X_{i}^{\prime}=X_{i}^{far}\cup\{c(x):x\in X_{i}^{near}\}$ $$

Note that $|X_{i}|=|X_{i}^{\prime}|$.
Let $m=\mu(X_{i})$, $m^{\prime}=\mu(X_{i}^{\prime})$.
Let $n=|X_{i}|$ and $\bar{n}=|X_{i}^{near}|$.
We first show a lower bound on $\Delta(X_{i})$ in terms of $R$.

###### Lemma 23 .

$\Delta(X_{i})\geq\frac{16\bar{n}}{\varepsilon^{2}}R^{2}$.

###### Proof.

Let $c=\arg\min_{c^{\prime}\in C_{init}}{\norm{m-c^{\prime}}}$. We do a case analysis:

- 1.
Case 1: $\norm{m-c}\geq\frac{5}{\varepsilon}\cdot R$
Consider any point $p\in X_{i}^{near}$. From triangle inequality, we have:
$\norm{p-m}\geq\norm{c(p)-m}-\norm{c(p)-p}\geq\frac{5}{\varepsilon}\cdot R-R
\geq\frac{4}{\varepsilon}\cdot R.$
This gives: $\Delta(X_{i})\geq\sum_{p\in X_{i}^{near}}\norm{p-m}^{2}\geq\frac{16\bar{n}}{
\varepsilon^{2}}\cdot R^{2}$.
- 2.
Case 2: $\norm{m-c}<\frac{5}{\varepsilon}\cdot R$
In this case, we have:
$\displaystyle\Delta(X_{i})=\Phi(c,X_{i})-n\cdot\norm{m-c}^{2}\geq\Phi(C_{init}
,X_{i})-n\cdot\norm{m-c}^{2}\geq\frac{41n}{\varepsilon^{2}}\cdot R^{2}-\frac{2
5n}{\varepsilon^{2}}\cdot R^{2}\geq\frac{16\bar{n}}{\varepsilon^{2}}\cdot R^{2}.$

This completes the proof of the lemma.
∎

We now bound the distance between $m$ and $m^{\prime}$ in terms of $R$.

###### Lemma 24 .

$\norm{m-m^{\prime}}^{2}\leq\frac{\bar{n}}{n}\cdot R^{2}$.

###### Proof.

Since $|X_{i}|=|X_{i}^{\prime}|$ and the only difference between $X_{i}$ and $X_{i}^{\prime}$ are the points corresponding to $X_{i}^{near}$, we have:

$$ $\norm{m-m^{\prime}}^{2}=\frac{1}{(n)^{2}}\left|\left|\sum_{p\in X_{i}^{near}}( p-c(p))\right|\right|^{2}\leq\frac{\bar{n}}{(n)^{2}}\sum_{p\in X_{i}^{near}} \norm{p-c(p)}^{2}\leq\frac{\bar{n}^{2}}{(n)^{2}}R^{2}\leq\frac{\bar{n}}{n}R^{2}.$ $$

The second inequality above follows from the Cauchy-Schwarz inequality.
∎

We now show that $\Delta(X_{i})$ and $\Delta(X_{i}^{\prime})$ are close.

###### Lemma 25 .

$\Delta(X_{i}^{\prime})\leq 4\bar{n}R^{2}+2\cdot\Delta(X_{i})$.

###### Proof.

The lemma follows from the following sequence of inequalities:

$$ $\displaystyle\Delta(X_{i}^{\prime})$ $\displaystyle=$ $\displaystyle\sum_{p\in X_{i}^{near}}\norm{c(p)-m^{\prime}}^{2}+\sum_{p\in X_{ i}^{far}}\norm{p-m^{\prime}}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:2})}}}{{\leq}}$ $\displaystyle\sum_{p\in X_{i}^{near}}2\cdot\left(\norm{c(p)-p}^{2}+\norm{p-m^{ \prime}}^{2}\right)+\sum_{p\in X_{i}^{far}}\norm{p-m^{\prime}}^{2}$ $\displaystyle\leq$ $\displaystyle 2\bar{n}R^{2}+2\cdot\Phi(m^{\prime},X_{i})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:1})}}}{{\leq}}$ $\displaystyle 2\bar{n}R^{2}+2\cdot(\Phi(m,X_{i})+n\cdot\norm{m-m^{\prime}}^{2})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma~{}\ref{lemma:case2-2})}}}{{ \leq}}$ $\displaystyle 4\bar{n}R^{2}+2\cdot\Delta(X_{i})$ $$

This completes the proof of the lemma.
∎

We now argue that any center that is good for $X_{i}^{\prime}$ is also good for $X_{i}$.

###### Lemma 26 .

Let $m^{\prime\prime}$ be such that $\Phi(m^{\prime\prime},X_{i}^{\prime})\leq\left(1+\frac{\varepsilon}{16}\right)
\cdot\Delta(X_{i}^{\prime})$.
Then $\Phi(m^{\prime\prime},X_{i})\leq\left(1+\frac{\varepsilon}{2}\right)\cdot
\Delta(X_{i})$.

###### Proof.

The lemma follows from the following inequalities:

$$ $\displaystyle\Phi(m^{\prime\prime},X_{i})$ $\displaystyle=$ $\displaystyle\sum_{p\in X_{i}}\norm{m^{\prime\prime}-p}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:1})}}}{{=}}$ $\displaystyle\sum_{p\in X_{i}}\norm{m-p}^{2}+n\cdot\norm{m-m^{\prime\prime}}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:2})}}}{{\leq}}$ $\displaystyle\Delta(X_{i})+2n\left(\norm{m-m^{\prime}}^{2}+\norm{m^{\prime}-m^ {\prime\prime}}^{2}\right)$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma~{}\ref{lemma:case2-2})}}}{{ \leq}}$ $\displaystyle\Delta(X_{i})+2\bar{n}R^{2}+2n\cdot\norm{m^{\prime}-m^{\prime \prime}}^{2}$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Fact~{}\ref{fact:1})}}}{{\leq}}$ $\displaystyle\Delta(X_{i})+2\bar{n}R^{2}+2\cdot\left(\Phi(m^{\prime\prime},X_{ i}^{\prime})-\Delta(X_{i}^{\prime})\right)$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma\ hypothesis)}}}{{\leq}}$ $\displaystyle\Delta(X_{i})+2\bar{n}R^{2}+\frac{\varepsilon}{8}\cdot\Delta(X_{i }^{\prime})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma~{}\ref{lemma:case2-3})}}}{{ \leq}}$ $\displaystyle\Delta(X_{i})+2\bar{n}R^{2}+\frac{\varepsilon}{2}\cdot\bar{n}R^{2 }+\frac{\varepsilon}{4}\cdot\Delta(X_{i})$ $\displaystyle\stackrel{{\scriptstyle\tiny{(Lemma~{}\ref{lemma:case2-1})}}}{{ \leq}}$ $\displaystyle\left(1+\frac{\varepsilon}{2}\right)\cdot\Delta(X_{i}).$ $$

This completes the proof of the lemma.
∎

Given the above lemma, all we need to argue is that our algorithm indeed considers a center $m^{\prime\prime}$ such that $\Phi(m^{\prime\prime},X_{i}^{\prime})\leq(1+\varepsilon/16)\cdot\Delta(X_{i}^{
\prime})$.
For this we would need about $\Omega(\frac{1}{\varepsilon})$ uniform samples from $X_{i}^{\prime}$.
However, our algorithm can only sample using $\mathcal{D}^{2}$-sampling w.r.t. $C_{init}$.
For ease of notation, let $c(X_{i}^{near})$ denote the multi-set $\{c(p):p\in X_{i}^{near}\}$.
Recall that $X_{i}^{\prime}$ consists of $X_{i}^{far}$ and $c(X_{i}^{near})$.
The first observation we make is that the probability of sampling an element from $X_{i}^{far}$ is reasonably large (proportional to $\frac{\varepsilon}{k}$).
Using this fact, we show how to sample from $X_{i}^{\prime}$ (almost uniformly).
Finally, we show how to convert this almost uniform sampling to uniform sampling (at the cost of increasing the size of sample).

###### Lemma 27 .

Let $x$ be a sample from $\mathcal{D}^{2}$-sampling w.r.t. $C_{init}$.
Then, $\mathbf{Pr}[x\in X_{i}^{far}]\geq\frac{\varepsilon}{8\alpha k}$.
Further, for any point $p\in X_{i}^{far}$, $\mathbf{Pr}[x=p]\geq\frac{\gamma}{|X_{i}|}$, where $\gamma$ denotes $\frac{(1-\delta)^{4}\varepsilon^{3}}{246\alpha k}$.

###### Proof.

Since points are $\mathcal{D}^{2}$-sampled with $\mathcal{D}$ being $\delta$-close to $D$, we
note that $\sum_{p\in X_{i}^{near}}\mathbf{Pr}[x=p]\leq\frac{(1+\delta)^{2}}{(1-\delta)^{
2}}\cdot\frac{R^{2}}{\Phi(C_{init},V)}\cdot|X_{i}|\leq\frac{(1+\delta)^{2}}{(1
-\delta)^{2}}\cdot\frac{\varepsilon^{2}}{41}\cdot\frac{\Phi(C_{init},X_{i})}{
\Phi(C_{init},V)}$.
Therefore, the fact that we are in case II implies that:

$$ $\mathbf{Pr}[x\in X_{i}^{far}]\geq\mathbf{Pr}[x\in X_{i}]-\mathbf{Pr}[x\in X_{i }^{near}]\geq\frac{\Phi(C_{init},X_{i})}{\Phi(C_{init},V)}-\frac{(1+\delta)^{2 }}{(1-\delta)^{2}}\cdot\frac{\varepsilon^{2}}{41}\cdot\frac{\Phi(C_{init},X_{i })}{\Phi(C_{init},V)}\geq\frac{\varepsilon}{8\alpha k}.$ $$

Also, if $x\in X_{i}^{far}$, then $\Phi(C_{init},\{x\})\geq R^{2}=\frac{\varepsilon^{2}}{41}\cdot\frac{\Phi(C_{
init},X_{i})}{|X_{i}|}$.
Therefore the probability that a point $p\in X_{i}^{far}$ gets $\mathcal{D}^{2}$-sampled is,

$$ $\displaystyle\mathbf{Pr}[x=p]\geq\frac{(1-\delta)^{2}}{(1+\delta)^{2}}\cdot \frac{\Phi(C_{init},\{x\})}{\Phi(C_{init},V)}$ $\displaystyle\geq$ $\displaystyle\frac{(1-\delta)^{2}}{(1+\delta)^{2}}\cdot\frac{\varepsilon}{6 \alpha k}\cdot\frac{R^{2}}{\Phi(C_{init},X_{i})}$ $\displaystyle\geq$ $\displaystyle\frac{(1-\delta)^{2}}{(1+\delta)^{2}}\cdot\frac{\varepsilon}{6 \alpha k}\cdot\frac{\varepsilon^{2}}{41}\cdot\frac{1}{|X_{i}|}$ $\displaystyle\geq$ $\displaystyle\frac{(1-\delta)^{2}}{(1+\delta)^{2}}\cdot\frac{\varepsilon^{3}}{ 246\alpha k}\cdot\frac{1}{|X_{i}|}\geq\frac{(1-\delta)^{4}\varepsilon^{3}}{246 \alpha k}\cdot\frac{1}{|X_{i}|}.$ $$

This completes the proof of the lemma.
∎

Let $O_{1},\ldots O_{\rho}$ be $\rho$ points sampled independently using $\mathcal{D}^{2}$-sampling w.r.t. $C_{init}$.
We construct a new set of random variables $Y_{1},\ldots,Y_{\rho}$.
Each variable $Y_{u}$ will depend on $O_{u}$ only and will take values either in $X_{i}^{\prime}$ or will be $\bot$.
These variables are defined as follows: if $O_{u}\notin X_{i}^{far}$, we set $Y_{u}$ to $\bot$.
Otherwise, we assign $Y_{u}$ to one of the following random variables with equal probability:
(i) $O_{u}$ or (ii) a random element of the multi-set $c(X_{i}^{near})$.
The following observation follows from Lemma 27.

###### Corollary 1 .

For a fixed index $u$, and an element $x\in X_{i}^{\prime}$, $\mathbf{Pr}[Y_{u}=x]\geq\frac{\gamma^{\prime}}{|X_{i}^{\prime}|},$ where $\gamma^{\prime}=\gamma/2$.

###### Proof.

If $x\in X_{i}^{far}$, then we know from Lemma 27 that $O_{u}$ is $x$ with probability at least $\frac{\gamma}{|X_{i}^{\prime}|}$ (note that
$X_{i}^{\prime}$ and $X_{i}$ have the same cardinality).
Conditioned on this event, $Y_{u}$ will be equal to $O_{u}$ with probability $1/2$.
Now suppose $x\in c(X_{i}^{near})$. Lemma 27 implies that $O_{u}$ is an element of $X_{i}^{far}$ with probability at least $\frac{\varepsilon}{8\alpha k}$.
Conditioned on this event, $Y_{u}$ will be equal to $x$ with probability at least $\frac{1}{2}\cdot\frac{1}{|c(X_{i}^{near})|}$.
Therefore, the probability that $O_{u}$ is equal to $x$ is at least $\frac{\varepsilon}{8\alpha k}\cdot\frac{1}{2|c(X_{i}^{near})|}\geq\frac{
\varepsilon}{16\alpha k|X_{i}^{\prime}|}\geq\frac{\gamma^{\prime}}{|X_{i}^{
\prime}|}$.
∎

Corollary 1 shows that we can obtain samples from $X_{i}^{\prime}$ which are nearly uniform (up to a constant factor).
To convert this to a set of uniform samples, we use the idea of [jks].
For an element $x\in X_{i}^{\prime}$, let $\gamma_{x}$ be such that $\frac{\gamma_{x}}{|X_{i}^{\prime}|}$ denotes the probability that the random variable $Y_{u}$ is equal to $x$ (note that this is independent of $u$).
Corollary 1 implies that $\gamma_{x}\geq\gamma^{\prime}$.
We define a new set of independent random variables $Z_{1},\ldots,Z_{\rho}$.
The random variable $Z_{u}$ will depend on $Y_{u}$ only.
If $Y_{u}$ is $\bot$, $Z_{u}$ is also $\bot$.
If $Y_{u}$ is equal to $x\in X_{i}^{\prime}$, then $Z_{u}$ takes the value $x$ with probability $\frac{\gamma^{\prime}}{\gamma_{x}}$, and $\bot$ with the remaining probability.
We can now prove the key lemma.

###### Lemma 28 .

Let $\rho$ be $\frac{256}{\gamma^{\prime}\cdot\varepsilon}$, and $m^{\prime\prime}$ denote the mean of the non-null samples from $Z_{1},\ldots,Z_{\rho}$. Then, with a probability at least $(3/4)$,
$\Phi(m^{\prime\prime},X_{i}^{\prime})\leq(1+\frac{\varepsilon}{16})\cdot\Delta
(X_{i}^{\prime})$.

###### Proof.

Note that a random variable $Z_{u}$ is equal to a specific element of $X_{i}^{\prime}$ with probability equal to $\frac{\gamma^{\prime}}{|X_{i}^{\prime}|}$.
Therefore, it takes $\bot$ value with probability $1-\gamma^{\prime}$.
Now consider a different set of iid random variables $Z_{u}^{\prime}$, $1\leq u\leq\rho$ as follows: each $Z_{u}$ tosses a coin with a probability of Heads being $\gamma^{\prime}$.
If we get Tails, it gets value $\bot$ otherwise, it is equal to a random element of $X_{i}^{\prime}$.
It is easy to check that the joint distribution of the random variables $Z_{u}^{\prime}$ is identical to that
of the random variables $Z_{u}$.
Thus, it suffices to prove the statement of the lemma for the random variables $Z_{u}^{\prime}$.

Now we condition on the coin tosses of the random variables $Z_{u}^{\prime}$.
Let $n^{\prime}$ be the number of random variables which are not $\bot$.
($n^{\prime}$ is a deterministic quantity because we have conditioned on the coin tosses).
Let $m^{\prime\prime}$ be the mean of such non-$\bot$ variables among $Z_{1}^{\prime},\ldots,Z_{\rho}^{\prime}$.
If $n^{\prime}$ happens to be larger than $\frac{128}{\varepsilon}$, Lemma 18 implies that with probability at least $(7/8)$,
$\Phi(m^{\prime\prime},X_{i}^{\prime})\leq(1+\frac{\varepsilon}{16})\cdot\Delta
(X_{i}^{\prime})$.

Finally, observe that the expected number of non-$\bot$ random variables is $\gamma^{\prime}\cdot\rho\geq\frac{256}{\varepsilon}$.
Therefore, with probability at least $\frac{7}{8}$ (using Chernoff-Hoeffding), the number of non-$\bot$ elements will be at least $\frac{128}{\varepsilon}$.
∎

Let $C^{(\rho)}$ denote the multi-set obtained by taking $\rho$ copies of each of the centers in $C_{init}$.
Now observe that all the non-$\bot$ elements among $Y_{1},\ldots,Y_{\rho}$ are elements of $\{O_{1},\ldots,O_{\rho}\}\cup C^{(\rho)}$, and so the same must hold for $Z_{1},\ldots,Z_{\rho}$.
Moreover, since we only need a uniform subset of size $\tau=\frac{128}{\varepsilon}$, $C_{i}^{\prime}$ suffices instead of $C^{(\rho)}$.
This implies that in steps 8-9 of the algorithm, we would have tried adding the point $m^{\prime\prime}$ as described in Lemma 28.
This means that $M_{i}$ contains a subset $T_{i}$ such that $\Phi(\mu(T_{i}),X_{i})\leq(1+\frac{\varepsilon}{2})\cdot\Delta(X_{i})$ with probability at least $3/4$.
This concludes the proof of Theorem 11.

###### Lemma 23 .

###### Proof.

###### Lemma 24 .

###### Proof.

###### Lemma 25 .

###### Proof.

###### Lemma 26 .

###### Proof.

###### Lemma 27 .

###### Proof.

###### Corollary 1 .

###### Proof.

###### Lemma 28 .

###### Proof.

## Appendix E Experiment Details

### E.1 Additional Evaluation

We also provide additional experiments on the DIGITS dataset from [](1797 data points, with each being a $8\times 8$ image of pixels having values between 0 and 16). As before, we plot the total cumulative runtime as described in section 4. To elaborate, the time corresponding to $k=4$ for QI-k-means++ is the sum of the time spent in setting up the sample and query data structure and then calculating the cluster centers for $k=2,3,4$. The trend observed is similar to that observed with binarized MNIST.

[Uncaptioned image]: extracted/2405.13351v1/runtime_digits.png

**Table 3: Clustering cost for DIGITS (scaled down by a factor of $10^{6}$)**
| k | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| QI-k-means++ | 3.52 | 3.20 | 2.94 | 2.93 | 2.77 | 2.43 | 2.38 | 2.35 | 2.27 |
| k-means++ | 3.47 | 3.14 | 2.96 | 2.74 | 2.67 | 2.35 | 2.37 | 2.27 | 2.17 |

### E.2 Variances of Clustering Costs

We also report the variance of the clustering costs mentioned in section 4 (recall that the experiments were conducted over 5 runs). The seeds were chosen randomly for each run.

**Table 4: Variance of Clustering cost for binarized MNIST (values are scaled down by a factor of $10^{11}$)**
| k | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| QI-k-means++ | 1.89 | 3.19 | 0.97 | 1.36 | 1.44 | 1.24 | 0.17 | 0.20 | 0.47 |
| k-means++ | 1.93 | 3.29 | 1.90 | 1.44 | 1.57 | 1.09 | 0.73 | 0.32 | 0.52 |

**Table 5: Variance of Clustering cost for IRIS (values are rounded to 2 decimal places)**
| k | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| QI-k-means++ | 246523 | 48805 | 335 | 830 | 236 | 99 | 72 | 55 | 51 |
| k-means++ | 423718 | 27544 | 420 | 561 | 765 | 54 | 93 | 35 | 47 |

**Table 6: Variance of Clustering cost for DIGITS (values are scaled down by a factor of $10^{10}$)**
| k | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| QI-k-means++ | 3.77 | 4.05 | 1.44 | 5.96 | 1.10 | 2.37 | 1.10 | 0.57 | 1.07 |
| k-means++ | 1.91 | 2.00 | 2.24 | 3.92 | 2.28 | 0.38 | 1.34 | 1.51 | 1.06 |

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: How can someone choose which algorithm to use between the ones in Theorem 2 and Theorem 3? For example, is there an efficient way to compute $\zeta$ in advance?
soundness: 3
strengths: The paper gives novel algorithms for k-means, which is a fundamental problem for machine learning. The paper is also fairly well written.
summary: This paper provides a quantum algorithm as well as a dequantized (i.e quantum inspired classical) version of it for the problem of $D^2$ sampling. This type of sampling is important for the kmeans++ algorithm and is defined as sampling points in $\mathbb{R}^d$ with probability proportional to the nearest center from a set $C$ of cluster centers in $\mathbb{R}^d$. The classical implementation, which gives a $O(log k)$ approximation guarantee, involves a preprocessing step that builds a sample-query access data structure in time O(Nd). Depending on the implementation, the running time after this data structure is set up is $\tilde{O}(\zeta^2k^2d)$ or $\tilde{O}(\zeta^6k^2)$, where $\zeta$ is the ratio between the largest and smallest pairwise distance in the pointset. The results are incomparable with prior work, since there is no algorithm that is the best in all parameter regimes. This work does asymptotically improve on all prior works for sufficiently small values of $\zeta$ and $k$ though. The authors also give a quantum fixed parameter approximation scheme to get a $(1+\varepsilon)$ approximation for the k-means problem for any $\varepsilon>0$. The k-means problem has been shown to be NP-hard to approximate below a 1.07 constant factor due to prior work and therefore, this approximation scheme has exponential dependence on $1/\varepsilon$, while it is efficient in terms of the other parameters. Finally, the authors provide experiments that verify the performance of the algorithms and make the comparison with the classical k-means++ algorithm.
weaknesses: The performance of the algorithm is only better in certain parameter regimes compared to prior work.

Minor comments:
-Line 94: I would call it “result” instead of “observation”
-In many places throughout the paper, I believe the tilde notation is a bit abused. It should only be used to hide factors logarithmic in the function displayed, not arbitrary ones as in Theorem 1 for example.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: The paper is nicely written, overall. However, the discussion through the prior work is not very clearly written. May be a table of relevant known results will be useful to understand the context of the algorithm in terms of efficiency. Also a more detailed discussion on aspect ratio and its relevance (and if there are fundamental obstacles in improving the running time in terms of this parameter) will be very useful. A table of relevant known results might include:

1. Key prior quantum and classical algorithms for k-means/k-means++
2. Their runtime complexities
3. Their approximation guarantees
4. Key parameters they depend on (e.g. aspect ratio)
soundness: 4
strengths: The paper designs new algorithms for a classical and significant computational problem: the  $k$-means problem. Any new algorithm for classical problems such as $k$-means is a contribution to the field of ML and data science.   Thus the paper has broad applicability (in particular will be of interest to ICLR audience).  This is the main strength of the paper.
summary: The present paper gives new algorithms; both quantum and classical, for  the $k$-means problem. $k$-means problem is a significant clustering problem with broad spectrum of applications. Given a set $V$ of points in ${\mathbf R}^{d}$, find a set $C$ of $k$ points, called centers, in $\mathbf{R}^d$, so that  the sum of squared distance from points in $V$ to the closest point in $C$ is minimized. Because of its foundational algorithmic nature, this computational problem has been of attention from researchers from several perspectives. 

The present paper investigates new quantum approaches to solving this problem. Given it is NP-hard, the goal is to develop approximation algorithms and an algorithm known as $k$-means++ (extending the well known $k$-means algorithm by Lloyd) by Arthur and Vassilvitskii has a $O(\log k)$ approximation guarantee in expectation. A basic ingredient of $k$-means++ algorithm is called $D^2$-$\mathit{sampling}$: a procedure that given a data set $V$ and a center set $C$, samples a point in $V$ with probability proportional to its squared distance from the nearest center in $C$. 

The main contribution of this paper is a new quantum algorithm for $D^2$-sampling and a subsequent de-quantization of it. This algorithm then leads to new implementation of  (both quantum and classical) of $k$-means++. While there are other known quantum implementation of  $D^2$-sampling, the present algorithm appears to have different efficiency parameters and hence appears to perform better in certain regimes than known ones. The main parameter that the  running time of the algorithm depends on is the ${\mathit aspect ratio}$ $\zeta$: the ratio between the largest and smallest inter-point distance in the dataset. In particular, the quantum algorithm that the author/s develop runs in time $\tilde{O}(\zeta^2k^2)$ with $\tilde{O}$ hiding the poly-log dependencies on relevant parameters. The author/s show how this can be de-quantized to get a classical algorithm with running time $O(Nd)+\tilde{O}(\zeta^2k^2d)$.   The paper also presents some initial experimental results.
weaknesses: The main weakness is  the algorithm's dependency on the aspect ratio. The algorithm has quadratic dependency on the aspect ratio. That makes one wonder whether it will be really  useful for broad applications. While the author/s briefly discuss this dependency in the paper, it is not at all clear why there needs to be a quadratic dependency. Is it possible to prove a lower bound (in a meaningful restricted model) showing that this dependency is needed? The authors cite a paper by Cohen-Addad et al that presents an algorithm that has closer to $\log$ dependency on $\zeta$ (at the expense of having a $N^{\theta(1)}$ factor in the running time) which appears to be much more attractive. The paper will be substantially stronger if it can bring more clarity on this issue. Following are concrete points that the authors may address to make the paper stronger. 

1. Provide a more detailed justification for why the algorithm presented has quadratic dependence on the aspect ratio.
2. Discuss whether the author/s  believe this dependence is fundamental or if there may be ways to improve it.
3. Consider proving a lower bound, to show whether this dependence is necessary (if the authors believe it is). 
4. Compare their aspect ratio dependence more thoroughly to other algorithms like the Cohen-Addad et al. result, analyzing the tradeoffs in more depth.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: There are some (minor) questions about the results:

1. In page 4, there is discussion about dependencies in the *assumption-free* case but seems like there are no results for the cases with assumptions. What comparisons and improvements can be made in scenarios with assumptions?

2. Considering different algorithms appear to have quite different dependencies with regards to N, d, k, aspect ratio, matrix conditions, etc., what are some impossibility results available for the k-means problem and how does the algorithm's runtime dependencies compare with it? 

3. Can you elaborate on whether there are methods to improve the dependency on the aspect ratio $\zeta$? It seems to be a bit counterintuitive. For example, if there are two clusters of points, and the two clusters are very well-separated, intuition might suggest that k-means with $k=2$ should be relatively straightforward and should be easier as the two clusters are separated out even more. But, in this case $\zeta$ grows and the bound on the runtime worsens. Additionally, $\zeta$ becomes very large if there is even a single pair of points that are close by. Do these observations suggest that the dependencies on $\zeta$ can be tightened?
soundness: 3
strengths: Starting from originality, the paper extends ideas in quantum machine learning (primarily the sample-access-query model of Tang) to the classical k-means problem, providing novel algorithms and results. Connections to existing literature are well-discussed, and discussions of the strengths and weaknesses of the new algorithm compared to classical algorithms are provided.

In terms of quality and clarity, it is quite clear that effort has been made to polish the text. The problem formulations, motivations, related work, and the design of the algorithm from a quantum-inspired perspective are all clearly stated. 

It is interesting to see that the new algorithm has better dependency in N but worse dependency in the aspect ratio and k. While this does not show a uniform improvement in all aspects, it is still a good contribution.
summary: The paper discusses a novel quantum-inspired algorithm for the k-means problem. Theoretical guarantees of the new algorithm are provided, showing it outperforms classical algorithms in certain areas. Experiments compare runtime of the quantum inspired k-means++ versus the classical k-means++.
weaknesses: One of the downsides of the algorithm's runtime dependencies is the dependence of the preprocessing cost O(Nd), which is ignored when comparing against existing results. If we add in O(Nd) to the equation, then both the new algorithm and classical algorithms has a linear dependence on $N$. Perhaps it would be meaningful to have experiments that showcase runtime *after preprocessing* so the plots in Figures 2-3 can be more reflective of the runtime dependencies after preprocessing.

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: Page 1:
I find the abstract to have too much low-level information :)

Page 4:
Can you please expand on related work?
How about the connections of your work to computational complexity?

Page 5:
Please give some extra background knowledge in quantum computing and Dirac bra-ket notation :)

Page 6:
Can you please elaborate on the concept of dequantization and provide some examples?

Page 7:
I do not see why the states $(1),\ \dots,\ (4)$ are the right ones involved in the quantum $D^2$-sampling algorithm :)

Page 8:
I do not understand the Lines 412 -- 417.
What does pseudo approximate mean here?

Page 9:
I am not an expert with experiments :)

Page 10:
Where is the Conclusion? :) Can you please add a Conclusion and Future Work section?
soundness: 3
strengths: $k$-means++ is an important tool in ML, therefore this paper is significant as it gives novel quantum or classical implementations of it.
summary: The authors give novel quantum and classical implementations of $k$-means++.
The main ingredient of their approach is a novel $D^2$-sampling method.
weaknesses: See the questions below.

Meta-review: This paper gave a quantum algorithm for D^2-sampling in clustering algorithms. The main contribution is on the theoretical side, in the sense that the proposed quantum algorithm has provable speedup over existing classical algorithms. Numerical experiments also demonstrated the benign performance of the proposed algorithm.

The reviewers in general found the results favorable for ICLR 2025. During the rebuttal period, the authors further updated the manuscript by 1) adding a new section on comparison with previous work, 2) adding a conclusion and future work section, 3) a table summarizing all results, and 4) adding a preliminary section for better accessible to a broader range of audience at ICLR. These further improved the paper, and the decision is to accept this paper at ICLR 2025.

Decision: Accept (Poster)