A Quantum Approximation Scheme for $k$-Means

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

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

We have already discussed past research works on quantum versions of the $k$-means algorithm (i.e., Lloyd’s iterations). This includes [ABG13], [LMR13], and [KLLP19]. 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 [WKS15] for examples of several such tools. Other directions on quantum $k$-means includes adiabatic algorithms (e.g., [LMR13]) and algorithms using the QAOA framework (e.g., [OMA⁺17, FGG14]). However, these are without provable guarantees. A line of work has also suggested that quantum algorithms can outperform classical ones because of the QRAM data structure access. A more level playing field is to assume that a similar sample and query data access is available in the classical setting. Under this assumption, several “dequantization” results for unsupervised machine learning algorithms have been given. This includes [Tan19, CGL⁺20, Tan21]. It will be interesting to see if similar dequantization is possible for the quantum algorithms presented in this work since the main ingredient of our algorithm and the dequantization results is length-squared 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}$ 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 [AV07] that this $k$-means++ algorithm above gives an $O(\log{k})$ approximation in expectation. It is also known from the result of [ADK09] 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 [BGJK20]. So, we will quantize the following constant factor pseudo-approximation algorithm.

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 [AV07] needed to show the pseudo-approximation property of Algorithm 1.

The above lemma is used as a black box in the analysis of Algorithm 1 in [ADK09]. 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 1 is used as a black box in the analysis of Algorithm 1, replacing this lemma with Lemma 2 also gives a constant factor approximation to the $k$-means objective. We will use the following notion of the closeness of two distance functions.

We will use this lemma in the approximation scheme of [BGJK20]. 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.

A high-level description of the approximation scheme of [BGJK20] was given in the introduction. We give a more detailed pseudocode in Algorithm 2.

In addition to the input instance $(V,k)$ and error parameter $\varepsilon$, the algorithm is also given a constant approximate solution $C$, which is used for $D^{2}$-sampling. A pseudoapproximate solution $C$ is sufficient for the analysis in [BGJK20]. 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 2 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 2. We state the main theorem of [BGJK20] before giving the analogous statement for the modified algorithm where $D$ is replaced with $\tilde{D}$ that is $\delta$-close to $D$.

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

The proof of the above theorem closely follows the proof of Theorem 2.1 of [BGJK20]. This is similar to the proof of Theorem 2 that we saw earlier, closely following the proof of Lemma 1. The minor changes are related to approximate distance estimates using $\tilde{D}$ instead of real estimates using $D$. The statement of Theorem 2.2 is not surprising in this light. Instead of repeating the entire proof of [BGJK20], we point out the one change in their argument caused by using $\tilde{D}$ instead of $D$ as the distance function. The analysis of [BGJK20] works by partitioning the points in any optimal cluster $X_{j}$ into those that are close to $C$ 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 [GJK20], which is a full version of [BGJK20]). 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$ gives a good center for this optimal cluster $X_{j}$. Replacing $D$ with $\tilde{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})$. 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 retains the same form as that in [BGJK20] (i.e., $|\mathcal{L}|=\tilde{O}\left(2^{\tilde{O}(\frac{k}{\varepsilon})}\right)$ and time $O(Nd|\mathcal{L}|)$).

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$:

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}$:⁷⁷7The 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}$.

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.

Suppose we want to compute the $k$-means cost, $\Phi(V,C)\equiv\sum_{i=1}^{N}\min_{j\in[k]}{D^{2}(v_{i},c_{j})}$, of the clustering given by a $k$ center set $C\in\mathbb{R}^{k\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.

Let $X_{1},...,X_{m}$ denotes the square of the measured value of the second register in $\ket{\Psi_{C}}$. These are random values in the range $[1,\eta^{2}]$, where $\eta=\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[k]}\tilde{D}(v_{i},v_{j})\in(1\pm\varepsilon)\cdot\Phi(V,C)$. We define the variable $S_{t}=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.

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^{2})$ 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 [BGJK20].

• •

  $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 5 is $m=\frac{\eta^{2}\ln{(10L)}}{\varepsilon^{2}}$, where $\eta=(1+\varepsilon)\cdot\max_{i,j}{D(v_{i},v_{j})}$.

$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 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:

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\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})}}{\tilde{\Phi}(V,C)}\in(1\pm\varepsilon)\cdot\frac{\min_{j\in[t]}{D(v_{i},c_{j})}}{\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 [BGJK20] to hold. We summarise the observations of this section in the next lemma. We will need the following notion of the relative similarity of two distributions.

The above lemma says that for $\Delta=O(\frac{1}{NkS})$, we obtain the required samples with high probability. We can now give proof of Theorem 1.1 assembling the quantum tools of this section.