Quantum and Quantum-Inspired Stereographic K Nearest-Neighbour Clustering

1. Introduction

Quantum Machine Learning (QML), using quantum algorithms to learn quantum or classical systems, has attracted much research in recent years, with some algorithms possibly gaining an exponential speedup [1,2,3]. Since machine learning routines often push real-world limits of computing power, an exponential improvement to algorithm speed would allow for such systems with vastly greater capabilities [4]. Google’s `Quantum Supremacy’ experiment [5] showed that quantum computers can naturally solve specific problems with complex correlations between inputs that can be incredibly hard for traditional (“classical”) computers. Such a result suggests that machine learning models executed on quantum computers could be more effective for specific applications. It seems quite possible that quantum computing could lead to faster computation, better generalisation on less data, or both, for an appropriately designed learning model. Hence, it is of great interest to discover and model the scenarios in which such a “quantum advantage” could be achieved. A number of such “Quantum Machine Learning” algorithms are detailed in papers such as [2,6,7,8,9]. Many of these methods claim to offer exponential speedups over analogous classical algorithms. However, some significant gaps exist between theoretical prediction and implementation on the path from theory to technology. These gaps result in unforeseen technological hurdles and sometimes misconceptions, necessitating more careful case-by-case studies such as [10].

It is known from the literature that the k-nearest-neighbour clustering algorithm (kNN) can be applied to solve the problem of phase estimation in optical fibres [11,12]. A quantum version of this kNN has been developed in [2], promising an exponential speedup. However, the practical usefulness of this algorithm is under debate [13]. The encoding of classical data into quantum states has been proven to be a complex task which significantly reduces the advantage of known quantum machine learning algorithms [13]. There are claims that the speedup is reduced to only polynomial once the quantum version of the algorithm takes into account the time taken to prepare the necessary quantum states. Furthermore, for near-intermediate scale quantum (NISQ) [3] applications, we should not expect the availability of QRAM, as this assumes reliable memories and operations which are still several milestones out of reach [14]. For this reason, it is not currently possible to use the fully quantum clustering algorithm and thus we resort to using hybrid quantum-classical kNN algorithms. Any classical implementation of kNN clustering involves, among other steps, repeated evaluations of a dissimilarity and and a loss function; changing the dissimilarity leads to a different clustering. A hybrid quantum classical kNN clustering algorithm utilizes quantum methods only to estimate the dissimilarity, eliminating the need for long-lasting quantum memories. However, reproducing the dissimilarity of a classical kNN algorithm using quantum methods can be prohibitively restrictive. The quantum dissimilarity also depends on the embedding (how the classical data is encoded in quantum states) and might only approximates the classical one, introducing fundamental deviations from the classical kNN algorithm. In [10], we applied a hybrid quantum-classical algorithm with modified angle embedding to the problem of k-means clustering for 64-QAM (Quadrature Amplitude Modulation) optical-fibre data (a well-known technical problem in signal processing through optical-fibre communication links) provided by Huawei [15], and show that this currently does not yield an advantage due to both the embedding and the current speed and noise of quantum devices.

In this work, we use the same problem and datasets to bring two main but independent contributions using the generalised inverse stereographic projection. First, we embed classical 2-dimensional data by computing the ISP onto the 3-dimensional sphere, and use the resulting normalised vector as Bloch vector to produce a pure quantum state of one qubit, which we call stereographic embedding. The resulting quantum dissimilarity directly translates into the cosine dissimilarity, thus making the quantum algorithm mathematically closer to the classical k-means algorithm. This means that no inherent limitation is introduced by the embedding and any loss in performance of this hybrid algorithm can be compensated for by improving the noise level and the speed of the quantum device. We thus propose stereographic embedding as an improved quantum embedding that may lead to improvement in several quantum machine learning algorithms (although there might still not be a practical quantum time advantage).

The second contribution comes from the benchmarking of the hybrid stereographic quantum mentioned above. Since, as already mentioned, the resulting hybrid clustering algorithm is mathematically equivalent to a classical `quantum-inspired’ kNN algorithm, in order to assess its performance in the absence of noise we simply test the equivalent classical quantum-inspired kNN algorithm. This algorithm is the result of first computing the ISP of the data and then performing clustering using a novel `quantum’ centroid update. We observe an increase in accuracy and convergence performance over k-means clustering on the 2-dimensional optical-fibre data. This suggests, as a purely classical second main contribution, that an advantage in decoding 64-QAM optical-fibre data is achieved by performing clustering in the inverse stereographically projected sphere and by using the spherical centroid.

The paper is structured as follows. In the remainder of this introduction, we discuss related works and our contribution to it. In Section 2, we introduce the experimental setup generating the 64-QAM optical-fibre transmission data and define clustering, the stereographic projection and the necessary quantum concepts for the hybrid protocols. Next, Section 3 introduces the developed Stereographic Quantum kNN (SQ-kNN), while Section 4 defines the developed quantum-inspired 2D Stereographic Classical kNN (2DSC-kNN) algorithm and proves its equivalence to the SQ-kNN quantum algorithm. In Section 5, we describe the various experiments for testing the algorithms, present the obtained results, and discuss their conclusions. We end the main text in Section 6 proposing some directions for future research, some of which are partially discussed in the appendix.

1.1. Related Work

A unifying overview of several quantum algorithms is presented in [16] in a tutorial style. An overview targeting data scientists is given in [17]. The idea of using quantum information processing methods to obtain speedups for the k-means algorithm was proposed in [18]. In general, neither the best nor even the fastest method for a given problem and problem size can be uniquely ascribed to either the class of quantum or classical algorithms, as seen in the detailed discussion presented in [9]. The advantages of using local (classical) processing units alongside quantum processing units in a distributed fashion are discussed in [19]. The accuracy of (quantum) K-means has been demonstrated experimentally in [20] and in [21], while quantum circuits for loading classical data into a quantum computer are described in [22].

An algorithm is proposed in [2] that solves the problem of clustering N-dimensional vectors to M clusters in $\mathcal{O}(log\left(MN\right))$ time on a quantum computer, which is exponentially faster than $\mathcal{O}\left(\mathrm{poly}\right(MN\left)\right)$ time for the (then) best known classical algorithm. The approach detailed in [2] requires querying the QRAM [23] for preparing a `mean state’, which is then used to find the inner product between the centroid (by default the mean point) using the SWAP test [24,25,26]. However, there exist some significant caveats to this approach. Firstly, this algorithm achieves an exponential speedup only when comparing the bit-to-bit processing time with the qubit-to-qubit processing time. If one compares the bit-to-bit execution times of both algorithms, the exponential speedup disappears, as shown in [13,27]. Secondly, since stable enough quantum memories do not exist, a hybrid quantum-classical approach must be used in real-world applications. Namely, all the information must be stored in classical memories, and the states to be used in the algorithm are prepared in real time. The process of preparing quantum states from classical data is known as `Data Embedding’ since we are embedding the classical data into quantum states. This, as mentioned before [4,27], slows down the algorithm to only a polynomial advantage over classical k-means. However, we propose an approach whereby this step of embedding can be treated as a data pre-processing step, allowing us to achieve some advantages in accuracy and convergence rate, taking a step towards making the quantum approach more viable. Instead of using a quantum algorithm, classical alternatives mimicking their behaviour, collectively known as quantum-inspired algorithms, have shown much promise in achieving classically some types of advantage that are demonstrated by quantum algorithms [4,27,28,29], but as [9] remarks, the massive increase in runtime with rank, condition number, Frobenius norm, and error threshold make the algorithms proposed in [13,27] impractical for matrices arising from real-world applications. This observation is supported by [30].

Recent works such as [27] suggest that even the best QML algorithms, without state preparation assumptions, fail to achieve exponential speedups over their classical counterparts. In [4], it is pointed out that most QML algorithms are incomparable to classical algorithms since they take quantum states as input and output quantum states, and that there is no analogous classical model of computation where one could search for similar classical algorithms. In [4], the idea of matching state preparation assumptions with ${\ell }^2$-norm sampling assumptions (first proposed in [27]) is implemented by introducing a new input model, sample and query access (SQ access). In [4], the Quantum K-Means algorithm described in [2] is `de-quantised’ using the `toolkit’ developed in [27], i.e. a classical quantum-inspired algorithm is given that, with classical SQ access assumptions replacing quantum state preparation assumptions, matches the bounds and runtime of the corresponding quantum algorithm up to polynomial slowdown. From the works [4,27,31], we can conclude that the exponential speedups of many quantum machine learning algorithms that are under consideration arise not from the `quantumness’ of the algorithms but instead from strong input assumptions, since the exponential part of the speedups vanish when classical algorithms are given analogous assumptions. In other words, in a wide array of settings, these algorithms do not give exponential speedups but rather yield polynomial speedups on classical data.

The fundamental aspect that allowed for the exponential speedup in [27] is exemplified by the problem of recommendation systems. The philosophy of classical recommendation algorithms before this breakthrough was to estimate all the possible preferences of a user and then suggest one or more of the most preferred objects. A quantum algorithm in [8] promised an exponential speedup but provided a recommendation without estimating all the preferences; namely, it only provided a sample of the most preferred objects. This process of sampling, along with state preparation assumptions, was, in fact, what gave the quantum algorithm an exponential advantage. The new classical algorithm obtains comparable speedups also by only providing samples rather than solving the whole preference problem. In [4], it is argued that the time taken to create the quantum state should be included for comparison since the time taken is not insignificant; it is also claimed that for every such linear algebraic quantum machine learning algorithm, a polynomially slower classical algorithm can be constructed by using the binary tree data structure described in [27]. Since then, more sampling algorithms have shown that multiple quantum exponential speedups are not due to the quantum algorithms themselves but due to the way data is provided to the algorithms and how the quantum algorithm provides the solutions [4,31,32,33]. Notably, in [33] it is argued that there exist competing classical algorithms for all linear algebraic subroutines and thus for many quantum machine learning algorithms. However, as pointed out in [9] and proven in [30], significant caveats exist to these aforementioned results of quantum-inspired algorithms. The polynomial factor in these algorithms often contains a very high power of the rank and condition number, making them suitable only for sparse low-rank matrices. Matrices of real-world data are often relatively high in rank and hence unfavourable for such sampling-based quantum-inspired approaches. Whether such sampling algorithms can be used also highly depends on the specific application and whether or not samples of the solution instead of the complete data are suitable. It should be pointed out that quantum machine learning algorithms generally do not provide an advantage if such complete data is needed.

The method of encoding classical data into quantum states contributes to the complexity and performance of the algorithm. An extensive analysis and testing of the hybrid quantum-classical implementation of the quantum k-means algorithm using angle embedding can be found in [10]. In this work, the use of the ISP is proposed. Others have explored this procedure [34,35,36] as well; however, the motivation, implementation, and use vary significantly, as well as the procedure for embedding data points into quantum states. There has also been no extensive testing of the proposed methods, especially not in an industry context. In our method, we exclusively use pure states from the Bloch sphere since this reduces the complexity of the application. 3 assures that our method with existing quantum techniques is applicable for nearest neighbour clustering. In contrast, the density matrices of mixed states and the normalised trace distance between the density matrices are used for binary classification in [34,35]. A crucial thing to consider here is to distinguish the contribution of the ISP from the quantum effects. We will see in Section 5 that the ISP seems to be the most important contributing factor. In [37], it is also proposed to encode classical information into quantum states using the ISP in the context of quantum generative adversarial networks. Their motivation for using the ISP is due to the fact that it is injective and can hence be used to uniquely represent every point in the 2D plane without any loss of information. On the other hand, angle embedding loses all amplitude information due to the normalisation of all points. A method to transform an unknown manifold into an n-sphere using ISP is proposed in [38] - here, however, the property of their concern was the conformality of the projection since subsequent learning is performed upon the surface. In [39], a parallelised version of [2] is developed using the FF-QRAM procedure [40] for amplitude encoding and the ISP to ensure a injective embedding.

In the method of Spherical Clustering [41], the nearest neighbour algorithm is explored based on the cosine similarity measure (). The cosine similarity is used in cases of information retrieval, text mining, and data mining to find the similarity between document vectors. It is used in those cases because the cosine similarity has low complexity for sparse vectors since only the non-zero co-ordinates need to be considered. For our case as well, it is in our interest to study with the cosine dissimilarity. This approach becomes particularly relevant once we employ stereographic embedding to encode the data points into quantum states.

1.2. Contribution

In this work, we first develop generalised stereographic embedding for hybrid quantum-classical kNN clustering as a better encoding that allows the quantum algorithm (Section 3) to outperform the accuracy and convergence of classical k-means algorithm in the absence of noise; in contrast, angle embedding introduces fundamental limitations to the accuracy not due to quantum noise. To validate this statement, we simulate this algorithm classically, which translates into an equivalent classical quantum-analogous stereographic kNN clustering algorithm (Section 4). One must note that we do not demonstrate that running the stereographic quantum kNN algorithm is more practical than the classical k-means algorithm in the NISQ context. We show that stereographic quantum kNN clustering converges faster and is more accurate than other hybrid quantum-classical kNN algorithms with angle or amplitude embedding. In parallel, the benchmarking of the classical stereographic kNN algorithm lets us claim that for the problem of decoding 64-QAM optical-fibre signals the generalised ISP and spherical centroid can allow for better accuracy and convergence.

The extensive testing upon the real-world, experimental QAM dataset (Section 2.1) revealed some significant results regarding the dependence of accuracy, runtime, and convergence performance upon the radius of projection, number of points, noise in the optical-fibre, and stopping criterion - described in Section 5. Noteworthy, we observe the existence of a finite optimal radius for the ISP (not equal to 1). To the best of our knowledge, no other work has considered a generalised projection radius for quantum embedding or studied its effect. Through our experimentation, we have verified that there exists an ideal radius greater than 1 for which accuracy performance is maximised. The advantageous implementation of the algorithm upon experimental data shows that our procedure is quite competitive. The fact that the developed quantum algorithm has an entirely classical analogue (with comparable time complexity to the classical k means algorithm) is a distinct advantage in terms of in-field deployment, especially compared to [2,9,18,34,35,36,39]. The developed quantum algorithm also has another advantage in the context of Noisy Intermediate-Scale Quantum (NISQ) realisations - it has the least circuit depth and circuit width among all candidates [2,9,36,39] - making it easier to implement with the current quantum technologies. Another significant contribution is our generalisation of the dissimilarity for clustering; instead of Euclidean dissimilarity (distance), we consider other dissimilarities which might be better estimated by quantum circuits (Appendix E). A somewhat similar approach was developed in parallel by [42] in the context of amplitude embedding. All previous approaches [2,9,36,39] only try to estimate the Euclidean distance. We also make the contribution of studying the relative effect of `quantumness’ and the ISP, something completely overlooked in previous works. We show that the quantum `advantage’ in accuracy performance touted by works such as [34,35,36,39] is in reality quite suspect and achievable through classical means. In the appendix, we describe a generalisation of the stereographic embedding - the Ellipsoidal embedding, which we expect to give even better results in future works.

Other secondary contributions of our work include:

• The development of a mathematical formalism for the generalisation of kNN to indicate the contribution of various parameters such as dissimilarities and dataspace (Section 2.4);
• Presenting the procedure and circuit for stereographic embedding using the Bloch embedding procedure, which consumes only $O\left(1\right)$ in time and resources (Section 3.1).

2. Preliminaries

In this section, for completeness, we touch upon some concepts and background information required to understand the paper. These concepts range from small general statements on quantum states (Bloch sphere, fidelity and Bell-state measurement in ) to the mathematical formalism of kNN (), and stereographic projection (Section 2.7). We begin by first describing the optic-fibre experimental setup used to collect the 64-QAM dataset, upon which the clustering algorithms were tested and benchmarked.

2.1. Optical-Fibre Setup

M-ary Quadrature Amplitude Modulation (M-QAM) is a simple and popular protocol for digital data transmission through analog communication channels. It is widely used in optical-fibre communication networks, and the decoding process of the received data often uses the k-nearest-neighbour algorithm to cluster nearby points. More details, including the description of the model used in the experiments, can be found in Appendix A. We now describe the experimental setup used to collect the dataset that is used for benchmarking the clustering algorithms.

The dataset contains a launch-power (laser-power feed into the fibre) sweep (four datasets collected at four different launch powers) of 80 km fibre transmission of coherent dual polarization (DP)-64QAM with a gross data rate of $960\times {10}^9$ bits/s. For the dataset, we assumed 15% overhead for forward error correction (FEC) and used 3.47% overhead for pilots and training sequences, so the net bit rate is $800\times {10}^9$ bits/s. Note that the pilots and training sequences are removed after the MIMO equalizer. An overview of the experimental setup [10,15] to capture this real-world database is shown in Figure 1. Four $120\times {10}^9$ Samples/s digital-to-analog converters (DACs) generate an electrical signal amplified by four 60 GHz 3dB-Bandwidth amplifiers. A tunable 100 kHz external cavity laser (ECL) source generates a continuous wave signal that is modulated by a 32 GHz DP-I/Q modulator. The receiver comprises an optical 90${}^{\circ }$-hybrid and four 100 GHz balanced photodiodes. The electrical signals are digitized using four 10-bit analog-to-digital converters (ADCs) with $256\times {10}^9$ Samples/s and 110 GHz. Subsequently, the raw signals are pre-processed by the receiver digital signal processing (DSP) blocks.

center/.style= align=center, block/.style = align=center, draw, fill=white, rounded corners, rectangle, minimum height=2em, minimum width=6em, EDFA/.style = align=center, draw, fill=white, regular polygon, regular polygon sides=3, minimum size=1.4cm, rotate=0, shape border rotate=-90, fiber/.style = align=center, draw, fill=none , circle, node distance= 2cm,color=red ,minimum size=20pt,

The datasets were collected in a very short time, corresponding to the memory size of the oscilloscope, which is limited. This is referred to as offline processing. At the receiver, the signals were normalised to fit the alphabet. The average launch power in watts can be calculated as follows:

$$\begin{array}{c}\hfill P_{\left(W\right)}=1W·{10}^{\left(P_{\left(dBm\right)}\right)/10)}/1000={10}^{(P_{\left(dBm\right)}-30)/10}\end{array}$$

There are four sets of published data with different launch powers, corresponding to different levels of non-linear distortions during transmission: $2.7$dBm, $6.6$dBm, $8.6$dBm, and $10.7$dBm. Each dataset consists of the `alphabet’ (initial analog transmission values), the error-corrected received analog values, and the true labels of the transmitted points. The data has been explained and visualised in detail in the Appendix A. To quantify the system performance of an amplified coherent optical communication system, one uses either a launch power sweep or an OSNR (Optical Signal to Noise Ratio) sweep. While the OSNR metric is used when the system is operating in the linear region, the launch power is the preferred metric to show the performance degradation in the nonlinear region since the induced nonlinear effects are directly proportional to the launch power. The signal-to-noise (SNR) ratio for each launch power can be computed using the following expression, where $\mathbf{z}$ are the received noisy signals and $\mathbf{x}$ are the noiseless target symbols (the launched signals):

$$\begin{array}{cc}\hfill SNR& =10{log}_{10}\left(\frac{\mathrm{mean}\left({∥\mathbf{z}∥}^2\right)}{\mathrm{mean}\left({∥\mathbf{z}-\mathbf{x}∥}^2\right)}\right).\hfill \end{array}$$

After obtaining this noisy real-world dataset, our task is to decode the received analog values into bit-strings. The kNN is the candidate of choice since it classifies datasets into clusters by associating an `average point’ (centroid) to each cluster. In our method, the objective of the clustering algorithm is first to identify, using the set of received signals, a given number M of centroids (one for each cluster) and then to assign each signal to the `nearest’ centroid. The second step is classification. This creates the clusters, which can then be decoded into bit signals through the process of demapping. Demapping consists of mapping the original transmission constellation (alphabet) to the current centroids and then assigning the bit-string label associated with that initial transmission point to all the points in the cluster of that centroid. This process completes the final step of the QAM protocol, translating the analog values to bit-strings read by the receiver. The size M of the constellation is known since we know beforehand which QAM protocol is being used. We also know the “alphabet”, i.e. the initial and ideal points at which the signals were transmitted.

2.2. Bloch Sphere

It is well known that all the qubit pure states can be obtained from the zero state using the unitary U [43]

$$\begin{array}{cc}\hfill |\psi (\theta ,\varphi )\rangle & =U(\theta ,\varphi )|0\rangle =cos(\theta /2)|0\rangle +e^{i\varphi }sin(\theta /2)|1\rangle \hfill \end{array}$$

(1)

where

$$\begin{array}{c}\hfill U(\theta ,\varphi )\left(\begin{array}{cc}cos\frac{\theta }{2}& -sin\frac{\theta }{2}\\ e^{i\varphi }sin\frac{\theta }{2}& e^{i\varphi }cos\frac{\theta }{2}\end{array}\right).\end{array}$$

(2)

These are unit vectors in the unit sphere of ${\mathbb{C}}^2$, but it is also well known that the corresponding density matrices are uniquely represented by the Bloch vectors

$$\mathbf{a}(\theta ,\varphi )(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta )$$

as points in the unit sphere $S^2\left(1\right)\subset {\mathbb{R}}^3$ [43] (the Bloch sphere) through the relation

$$\begin{array}{cc}\hfill \rho (\theta ,\varphi )& =\left|\psi (\theta ,\varphi )\rangle \langle \psi (\theta ,\varphi )\right|\hfill \end{array}$$

(3)

$$\begin{array}{cc}& =\left(\begin{array}{cc}{cos}^2\frac{\theta }{2}& -e^{-i\varphi }cos\frac{\theta }{2}sin\frac{\theta }{2}\\ e^{i\varphi }cos\frac{\theta }{2}sin\frac{\theta }{2}& {sin}^2\frac{\theta }{2}\end{array}\right)\hfill \end{array}$$

(4)

$$\begin{array}{cc}& =\frac{1}{2}\left(\begin{array}{cc}1+cos\theta & -e^{-i\varphi }sin\theta \\ e^{i\varphi }sin\theta & 1-cos\theta \end{array}\right)\hfill \end{array}$$

(5)

$$\begin{array}{cc}& =\frac{1}{2}\left(\mathbb{1}+\mathbf{a}(\theta ,\varphi )·\overrightarrow{\sigma }\right)\hfill \end{array}$$

(6)

where $\mathbb{1}$ is the identity matrix and $\overrightarrow{\sigma }=({\sigma }_x,{\sigma }_y,{\sigma }_z)$ is the vector of Pauli matrices

$$\begin{array}{cccccc}\hfill {\sigma }_x& =\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\hfill & \hfill {\sigma }_y& =i\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\hfill & \hfill {\sigma }_x& =\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).\hfill \end{array}$$

Regarding mixed states, notice that Eq. (6) is linear and thus, convex combinations of density matrices translate to convex combinations of Bloch vectors, meaning that the interior of the sphere represents the mixed states. Namely, the most general qubit quantum states can be represented by

$$\begin{array}{cccc}\hfill {\rho }_{\mathbf{a}}\equiv \rho \left(\mathbf{a}\right)& =\frac{1}{2}\left(\mathbb{1}+\mathbf{a}·\overrightarrow{\sigma }\right),\hfill & \hfill & {∥\mathbf{a}∥}_2\le 1.\hfill \end{array}$$

(7)

Finally, since the Pauli matrices are orthogonal operators under the Hilbert-Schmidt inner product, this inner product is easily computed as

$$\begin{array}{cc}\hfill Tr\left({\rho }_{{\mathbf{a}}_1}{\rho }_{{\mathbf{a}}_2}\right)& =\frac{1}{2}\left(1+{\mathbf{a}}_{\mathbf{1}}·{\mathbf{a}}_{\mathbf{2}}\right).\hfill \end{array}$$

(8)

which for pure states coincides with the fidelity.

Using the Bloch sphere representation of qubit quantum states also makes it easy to find orthogonal states and compute diagonalizations. Indeed, let $\mathbf{a}$ be a unit vector ($∥\mathbf{a}∥=\mathbf{a}·\mathbf{a}=1$), thus representing the pure state ${\rho }_{\mathbf{a}}=\frac{1}{2}\left(\mathbb{1}+\mathbf{a}·\overrightarrow{\sigma }\right)$, then the orthogonal state to ${\rho }_{\mathbf{a}}$ is simply the antipodal point

$$\begin{array}{c}\hfill {\rho }_{-\mathbf{a}}=\frac{1}{2}\left(\mathbb{1}-\mathbf{a}·\overrightarrow{\sigma }\right)\end{array}$$

(9)

which can be shown by computing the inner product as in Eq. (8)

$$\begin{array}{c}\hfill Tr\left({\rho }_{+\mathbf{a}}{\rho }_{-\mathbf{a}}\right)=\frac{1}{4}\left(1+\mathbf{a}·(-\mathbf{a})\right)=0.\end{array}$$

(10)

Hence, the Bloch eigenvectors for any Bloch vector $\mathbf{a}$ are $\pm \frac{\mathbf{a}}{∥\mathbf{a}∥}$, the two antipodal points where the line of $\mathbf{a}$ intersect the Bloch sphere. Namely, for any mixed quantum state corresponding to the Bloch vector $\mathbf{a}$, we can decompose the quantum state as

$$\begin{array}{cc}\hfill \frac{1}{2}\left(\mathbb{1}+\mathbf{a}·\overrightarrow{\sigma }\right)& =p\frac{1}{2}\left(\mathbb{1}+\frac{\mathbf{a}}{∥\mathbf{a}∥}·\overrightarrow{\sigma }\right)+(1-p)\frac{1}{2}\left(\mathbb{1}-\frac{\mathbf{a}}{∥\mathbf{a}∥}·\overrightarrow{\sigma }\right)\hfill \end{array}$$

(11)

$$\begin{array}{cc}& =\frac{1}{2}\left(\mathbb{1}+(2p-1)\frac{\mathbf{a}}{∥\mathbf{a}∥}·\overrightarrow{\sigma }\right)\hfill \end{array}$$

(12)

with

$$\begin{array}{ccc}\hfill 2p-1& =∥\mathbf{a}∥\hfill & \hfill \Rightarrow p=\frac{1}{2}(1+∥\mathbf{a}∥).\end{array}$$

(13)

In the next section, we discuss how we use the Bell-state measurement to estimate the fidelity between quantum states and exposit when this should be chosen over the SWAP test.

2.3. Bell-State Measurement and Fidelity

We use the Bell-state measurement to estimate the fidelity between two pure states. The Bell-state measurement is defined as the von-Neumann measurement of the maximally entangled basis

$$|{\varphi }_{ij}\rangle CNOT(H\otimes \mathbb{1})|ij\rangle ,$$

(14)

which by construction is equivalent to a standard basis measurement after $(H\otimes \mathbb{1})CNOT$ as displayed in Figure 2. This measurement can be used to estimate the fidelity as follows.

Lemma 1.

Let $|\psi \rangle$ and $|\chi \rangle$ be two qubit pure states and let $|{\varphi }_{11}\rangle CNOT(H\otimes \mathbb{1})|11\rangle$ (the singlet Bell state). Then

$$\begin{array}{cc}\hfill {\left|\langle {\varphi }_{11}|\left(|\psi \rangle \otimes |\chi \rangle \right)\right|}^2& =\frac{1}{2}{(1-|\langle \psi |\chi \rangle |}^2).\hfill \end{array}$$

(15)

Proof. 

Let us write the states as

$$\begin{array}{cccc}\hfill |\psi \rangle & ={\psi }_0|0\rangle +{\psi }_1|1\rangle \hfill & \hfill |\chi \rangle & ={\chi }_0|0\rangle +{\chi }_1|1\rangle ,\hfill \end{array}$$

(16)

Then the state before the standard-basis measurement is

$$\begin{array}{cc}\hfill |{\psi }_{\mathrm{out}}\rangle =\left(H\otimes \mathbb{1}\right)CNOT\left(|\psi \rangle \otimes |\chi \rangle \right)& =\frac{1}{\sqrt{2}}\left(\begin{array}{c}{\psi }_0{\chi }_0+{\psi }_1{\chi }_1\\ {\psi }_0{\chi }_1+{\psi }_1{\chi }_0\\ {\psi }_0{\chi }_0-{\psi }_1{\chi }_1\\ {\psi }_0{\chi }_1-{\psi }_1{\chi }_0\end{array}\right)\hfill \end{array}$$

(17)

and in particular, the probability of outcome $ij=11$ (i.e. simultaneous measurement of both qubits yields value `1’ on each qubit) can be written as

$$\begin{array}{cc}\hfill {\left|\langle {\varphi }_{11}|\left(|\psi \rangle \otimes |\chi \rangle \right)\right|}^2& ={\left|\langle 11|{\psi }_{out}\rangle \right|}^2=\frac{1}{2}{|{\psi }_0{\chi }_1-{\psi }_1{\chi }_0|}^2.\hfill \end{array}$$

(18)

The fidelity is obtained now by adding and subtracting ${\psi }_0^{*}{\psi }_0{\chi }_0^{*}{\chi }_0+{\psi }_1^{*}{\psi }_1{\chi }_1^{*}{\chi }_1$ and computing

$$\begin{array}{cc}\hfill {\left|\langle {\varphi }_{11}|\left(|\psi \rangle \otimes |\chi \rangle \right)\right|}^2& =\frac{1}{2}\left(1-{\psi }_1{\chi }_0{\psi }_0^{*}{\chi }_1^{*}-{\psi }_0{\chi }_1{\psi }_1^{*}{\chi }_0^{*}-{\psi }_0^{*}{\psi }_0{\chi }_0^{*}{\chi }_0-{\psi }_1^{*}{\psi }_1{\chi }_1^{*}{\chi }_1\right)\hfill \end{array}$$

(19)

$$\begin{array}{cc}& =\frac{1}{2}(1-|{\psi }_0^{*}{\chi }_0+{\psi }_1^{*}{\chi }_1{|}^2)\hfill \end{array}$$

(20)

$$\begin{array}{cc}& =\frac{1}{2}{(1-|\langle \psi |\chi \rangle |}^2),\hfill \end{array}$$

(21)

concluding the proof.    □

Lemma 1 is used to construct the quantum clustering algorithm in Section 3. We will use the quantum circuit of Figure 2 for fidelity estimation in the developed quantum algorithm.

Remark 1.

Since we are only interested in the $ij=11$ outcome and we are measuring qubits, the course-grained projective measurement defined by $\left|{\varphi }_{11}\rangle \langle {\varphi }_{11}\right|$ and $\mathbb{1}-\left|{\varphi }_{11}\rangle \langle {\varphi }_{11}\right|$ is sufficient for computing the inner product. The non-destructive version of this measurement is known as theSWAP test [25,26], first described in [24]. This test has been used extensively for overlap estimation in quantum algorithms [2]. The SWAP test requires one to only measure an ancilla qubit instead of the two input qubits, leaving them in the post-measurement state, which can be used later for other purposes. However, given the current limitations of NISQ technologies, storing quantum information for reuse is quite impractical; therefore, we prefer the destructive measurement version for overlap estimation. Namely, we use the Bell-state measurement instead of the SWAP test because the post-measurement state is unnecessary.

2.4. Nearest-Neighbour Clustering Algorithms

Clustering is a simple, powerful and well-known machine-learning algorithm that has been extensively used throughout the literature. In this section, we summarise some standard and basic notions introduced by clustering and define this class of heuristic algorithms precisely so that we can make clear the difference between regular clustering and the quantum and quantum-inspired clustering algorithms introduced in this paper. We first define the involved variables needed for the kNN.

Definition 1

(Clustering State). We define ak-Nearest-Neighbour Clustering State, orclustering statefor short, as a collection $(D,\overline{\mathbf{c}},\mathbb{D},d)$ where

• $\mathbb{D}$ is a space calleddataspacewith elements calledpoints.
• $D\subseteq \mathbb{D}$ is a subset calleddatasetconsisting of points calleddatapoints.
• $\overline{\mathbf{c}}=({\mathbf{c}}_1{\mathbf{c}}_2\dots {\mathbf{c}}_k)\subseteq {\mathbb{D}}^k$ is a list (of size k) of points calledcentroids
• $d:\mathbb{D}\times \mathbb{D}⟼\mathbb{R}$ is a lower bounded function calleddissimilarity function, ordissimilarityfor short.

Note that d does not have to be a distance metric. We now define the basic steps that are repeated in the clustering algorithm.

Definition 2

(Clusters and Centroid update). Let $(D,\overline{\mathbf{c}},\mathbb{D},d)$ be a clustering state. We define theclustersof the state as, for each $j=1,...,k$, the set

$$\begin{array}{c}\hfill C_j\left(\overline{\mathbf{c}}\right)=\left\{\mathbf{p}\in D|d(\mathbf{p},{\mathbf{c}}_j)\le d(\mathbf{p},{\mathbf{c}}_{\ell })\forall \ell =1,...,k,\mathbf{p}\notin \bigcup _{\ell <j}{C}_{\ell }\left(\overline{\mathbf{c}}\right)\right\}.\end{array}$$

(22)

We now definethe possible new centroidsof a subset $C\subseteq D$ as the set

$$\begin{array}{cc}\hfill \mathcal{P}\left(C\right)& \underset{\mathbf{x}\in \mathbb{D}}{argmin}\sum _{\mathbf{p}\in C}d(\mathbf{x},\mathbf{p})\hfill \end{array}$$

(23)

of all points minimising the total (and thus the average) dissimilarity. Then, we call acentroid updateany function ${\mathbf{c}}^{\mathrm{update}}:P\left(\mathbb{D}\right)\to \mathbb{D}$ (where P denotes the power set) of clusters such that ${\mathbf{c}}^{\mathrm{update}}\left(C\right)$ is a possible new centroid, namely such that ${\mathbf{c}}^{\mathrm{update}}\left(C\right)\in \mathcal{P}\left(C\right)$, for all $j=1,...,k$. We then define the following short-hand notation for the centroid update of $\overline{\mathbf{c}}$, namely the new list of centroids

$$\begin{array}{cc}\hfill {\overline{\mathbf{c}}}^{\mathrm{update}}\left(\overline{\mathbf{c}}\right)& =\left({\mathbf{c}}^{\mathrm{update}}\left(C_1\left(\overline{\mathbf{c}}\right)\right),...,{\mathbf{c}}^{\mathrm{update}}\left(C_k\left(\overline{\mathbf{c}}\right)\right)\right).\hfill \end{array}$$

(24)

We now define the general k-nearest-neighbour clustering algorithm.

Definition 3

(K-Nearest-Neighbour Clustering Algorithm (kNN)). Finally, we define aK-Nearest-Neighbour clustering algorithm (kNN)as a pair of clustering state and centroid update $(D,{\overline{\mathbf{c}}}_1,\mathbb{D},d,{\overline{\mathbf{c}}}^{\mathrm{update}})$. The kNN algorithm defines a sequence of clustering states $(D,{\overline{\mathbf{c}}}_i,\mathbb{D},d)$ via ${\overline{\mathbf{c}}}_{i+1}={\overline{\mathbf{c}}}^{\mathrm{update}}\left({\overline{\mathbf{c}}}_i\right)$ for all $i\in \mathbb{N}$ which we call the iterations (of the algorithm).

A point of note is that Eq. (23) implies that the new centroid is one of the points $\mathbf{x}$ in the dataspace that minimises the total (and hence the average) dissimilarity with all the points $\mathbf{p}$ in the cluster. Also, notice that this definition requires one to initialise or populate the list $\overline{\mathbf{c}}$ with initial values, i.e. the initial centroids${\overline{\mathbf{c}}}_1$ must be defined as a starting point for the clustering algorithm. The initial centroids can be assigned randomly or defined as a function of parameters such as the dataset.

Another comment about Eq. (23): in our case, we will see later in Section 3.2 that all choices of points from the set ${\mathcal{P}}_j$ will be equivalent. As in our algorithm, this freedom of choice can be exploited to reduce the amount of computation or for other optimisations.

Notice that Eqs. (23) and (24) implies that centroids are generally not part of the original data set; however, according to Eqs. (23) and (24) they must be restricted to the space in which the dataset is defined. Definitions involving centroids for which $\overline{\mathbf{c}}\notin {\mathbb{D}}^k$ are possible but are not used in this work.

One can see that any kNN can be broken down into two steps that keep alternating until a stopping condition (a condition which, when true, forces the algorithm to terminate) is met: a cluster update which updates the points associated with the newly calculated centroid, and then a centroid update which recalculates the centroid based upon the new points associated to it through its cluster. For the cluster update, the value of the centroid calculated in the previous iteration is taken, and its cluster set is constructed by collecting all the points in the dataset that are `closer’ to it than any other centroid. The `closeness’ is computed by using a pre-defined dissimilarity. In the next step, the centroids are updated by searching in the dataspace, for each updated cluster, a new point for which the sum of dissimilarities between that point and all points in the cluster is minimised.

This procedure will lead to different results if one changes the dissimilarity or the space of data points or both. In this paper, we explore the effects of changing this dissimilarity as well as the space of data points, and we shall explain it in the context of quantum states.

2.5. Euclidean Dissimilarity and Classical Clustering

It can be seen from the centroid update in Eqs. (22) and (23) that the dissimilarity plays a central role in the clustering algorithm. The nature of this function directly controls the first step of cluster update since the dissimilarity is used to compute the `closeness’ between any two points in the dataspace. It is also apparent that if the dissimilarity is changed in the centroid update, the points at which the minimum is achieved could also change.

The Euclidean dissimilarity $d_{\mathrm{e}}:{\mathbb{R}}^n\times {\mathbb{R}}^n\to \mathbb{R}$ is defined simply as the square of the Euclidean distance between the points:

$$\begin{array}{c}\hfill d_{\mathrm{e}}(\mathbf{a},\mathbf{b})={∥\mathbf{a}-\mathbf{b}∥}^2.\end{array}$$

(25)

For a finite subset $C\subset {\mathbb{R}}^n$, the minimisation of Eqs. (23) and (24) yields a unique point, reducing to the average point of the cluster:

$$\begin{array}{c}\hfill {\mathbf{c}}_{\mathrm{e}}^{\mathrm{update}}\left(C\right):=\underset{\mathbf{x}\in {\mathbb{R}}^n}{argmin}\sum _{\mathbf{p}\in C}{d}_{\mathrm{e}}(\mathbf{x},\mathbf{p})=\frac{1}{\left|C\right|}\sum _{\mathbf{p}\in C}\mathbf{p},\end{array}$$

(26)

which we call the Euclidean centroid update. This is the most typical case of centroid update, where the new centroid is updated as the mean point of all points in the cluster. This corresponds to the classic k-means clustering algorithm [44], which now can be defined as follows.

Definition 4

(n-Dimensional Euclidean Classical kNN (nDEC-kNN)). Ann-dimensional classical Euclidean kNN algorithmis any clustering algorithm with dataspace ${\mathbb{R}}^n$, Euclidean dissimilarity and cluster update as the average point, as in Eq. (26). Namely, any clustering algorithm of the form $(D,\overline{\mathbf{c}},{\mathbb{R}}^n,d_{\mathrm{e}},{\mathbf{c}}_{\mathrm{e}}^{\mathrm{update}})$.

The computation of the centroid through Eq. (26) instead of Eqs. (23) and (24) reduces the complexity of the centroid update step; such a reduced expression is used to compute the updated centroids rather than searching the entire dataspace for the minimising points during the centroid update.

2.6. Cosine Dissimilarity

In this work, we project the collected two-dimensional dataset (described in Section 2.1) into a sphere via the ISP. After this projection, the calculation of the centroids according to Eq. (26) would generally yield centroids which lie inside the sphere instead of on the $S^2\left(r\right)$ surface due to the convex nature of the sphere’s surface.

In our work, to use qubit pure states, we restrict the dataspace $\mathbb{D}$ to the sphere surface $S^2\left(r\right)$, forcing the centroids to lie on the surface of a sphere. This naturally leads to the question of what the proper reformulation of Eqs. (23) and (24) is, and whether a computationally inexpensive formula similar to Eq. (26) exists for this case as well. This question will be answered in Lemma 3. For this purpose, it is useful to first define the cosine dissimilarity [45] and see how it relates to the Euclidean dissimilarity.

Definition 5

(Cosine Dissimilarity). For two points, $\mathbf{a}$ and $\mathbf{b}$ in an inner-product space $\mathbb{D}$ thecosine dissimilarityis defined as:

$$\begin{array}{c}\hfill d_{\mathrm{s}}(\mathbf{a},\mathbf{b})=1-\frac{\mathbf{a}·\mathbf{b}}{∥\mathbf{a}∥∥\mathbf{b}∥},\end{array}$$

(27)

where $\mathbf{a}·\mathbf{b}$ is the inner product between the two points expressed as vectors from the origin, $∥\mathbf{a}∥$ is the norm of $\mathbf{a}$ induced by the inner product.

This is called cosine dissimilarity because when $\mathbf{a},\mathbf{b}\in {\mathbb{R}}^n$ the cosine dissimilarity $d_{\mathrm{s}}(\mathbf{a},\mathbf{b})$ reduces to $1-cos\left(\alpha \right)$, where $\alpha$ is the angle between $\mathbf{a}$ and $\mathbf{b}$. The cosine dissimilarity is also known sometimes as cosine distance (although it is not a distance), while $\frac{\mathbf{a}·\mathbf{b}}{∥\mathbf{a}∥∥\mathbf{b}∥}$ is well known as cosine similarity. This quantity, by construction, only depends on the direction of the vectors and not their magnitude. Said otherwise, we have

$$\begin{array}{c}\hfill d_{\mathrm{s}}(\mathbf{a},\mathbf{b})=d_{\mathrm{s}}(c\mathbf{a},\mathbf{b})=d_{\mathrm{s}}(\mathbf{a},c\mathbf{b})\end{array}$$

(28)

for any positive constant $c>0$. We also note that cosine dissimilarity of Eq. (27) can be related to the Euclidean dissimilarity of Eq. (25) if $\mathbf{a}$ and $\mathbf{b}$ lie on the n-sphere $S^n\left(r\right)\left\{\mathbf{s}\in {\mathbb{R}}^{n+1}|{∥\mathbf{s}∥}_2=r\right\}$ of radius r, as stated by the following lemma.

Lemma 2.

Let $d_{\mathrm{s}}$ and $d_{\mathrm{e}}$ be the cosine and Euclidean dissimilarities, respectively. Let ${\mathbf{s}}_1$, ${\mathbf{s}}_2$$\in S^n\left(r\right)$ be points on the n-sphere of radius r, then

$$\begin{array}{c}\hfill d_{\mathrm{e}}({\mathbf{s}}_1,{\mathbf{s}}_2)=2r^2{d}_{\mathrm{s}}({\mathbf{s}}_1,{\mathbf{s}}_2).\end{array}$$

(29)

Proof. 

Assuming ${\mathbf{s}}_1$, ${\mathbf{s}}_2$$\in S^n\left(r\right)$, Eq. (27) reduces to:

$$\begin{array}{cc}\hfill d_{\mathrm{s}}({\mathbf{s}}_1,{\mathbf{s}}_2)& =1-\frac{1}{r^2}{\mathbf{s}}_1·{\mathbf{s}}_2,\hfill \end{array}$$

(30)

then

$$\begin{array}{cc}\hfill 2r^2{d}_{\mathrm{s}}({\mathbf{s}}_1,{\mathbf{s}}_2)& =2r^2-2{\mathbf{s}}_1·{\mathbf{s}}_2={∥{\mathbf{s}}_1∥}^2+{∥{\mathbf{s}}_2∥}^2-2{\mathbf{s}}_1·{\mathbf{s}}_2={∥{\mathbf{s}}_1-{\mathbf{s}}_2∥}^2=d_{\mathrm{e}}({\mathbf{s}}_1,{\mathbf{s}}_2),\hfill \end{array}$$

(31)

concluding the proof.    □

From this, we can expect that the minimiser of the centroid update equation (23) computed using the cosine dissimilarity will closely relate to the Euclidean centroid update. However, the derivation is not straightforward since the Euclidean centroid update does not lie on the same sphere, but lies inside at a smaller radial distance. This is shown in the following lemma.

Lemma 3.

Let $C\subset S^n\left(r\right)$ be a finite set, then

$$\begin{array}{cc}\hfill {\mathbf{c}}_{\mathrm{s}}^{\mathrm{update}}\left(C\right):=& \underset{\mathbf{x}\in S^n\left(r\right)}{argmin}\sum _{\mathbf{p}\in C}{d}_{\mathrm{s}}(\mathbf{x},\mathbf{p})=r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}.\hfill \end{array}$$

(32)

We call this the cosineorspherical centroid update. In particular, thus

$$\begin{array}{cc}\hfill 1-1{\mathbf{c}}_{\mathrm{s}}^{\mathrm{update}}\left(C\right)=& =r\frac{{\mathbf{c}}_{\mathrm{e}}^{\mathrm{update}}\left(C\right)}{∥{\mathbf{c}}_{\mathrm{e}}^{\mathrm{update}}\left(C\right)∥}\hfill \end{array}$$

(33)

where ${\mathbf{c}}_{\mathrm{e}}^{\mathrm{update}}\left(C\right)=\frac{1}{\left|C\right|}{\sum }_{\mathbf{p}\in C}\mathbf{p}$ is the Euclidean centroid update of Eq. (26).

Proof. 

The second claim is trivial, we thus have to prove only the first claim. Given that $C\subseteq S^n\left(r\right)$, then according to Lemma 2, the cosine dissimilarity given in Eq. (27) reduces for all $\mathbf{a},\mathbf{b}\in C$ to:

$$\begin{array}{cc}\hfill d_{\mathrm{s}}(\mathbf{a},\mathbf{b})& =1-\frac{\mathbf{a}·\mathbf{b}}{∥\mathbf{a}∥∥\mathbf{b}∥}=1-\frac{1}{r^2}\mathbf{a}·\mathbf{b}.\hfill \end{array}$$

(34)

The minimisation in Eq. (23) can then be calculated for the cosine dissimilarity with a Lagrangian (see Eq. (37)) that satisfies Eqs. (23) and (24) at the minimising point. Namely, we have to find $\mathbf{x}\in {\mathbb{R}}^n$ that minimises

$$\begin{array}{c}\hfill f\left(\mathbf{x}\right)=\sum _{\mathbf{p}\in C}d(\mathbf{x},\mathbf{p}),\end{array}$$

(35)

subject to the restriction condition that assures that $\mathbf{x}\in S^n\left(r\right)$, that is

$$\begin{array}{c}\hfill g\left(\mathbf{x}\right)={∥\mathbf{x}∥}^2-r^2=0.\end{array}$$

(36)

Such a Lagrangian is expressed as

$$\begin{array}{cc}\hfill \mathcal{L}(\mathbf{x},\lambda )& =f\left(\mathbf{x}\right)-\lambda g\left(\mathbf{x}\right)=\sum _{\mathbf{p}\in C}{d}_{\mathrm{s}}(\mathbf{x},\mathbf{p})-\lambda ({∥\mathbf{x}∥}^2-r^2),\hfill \end{array}$$

(37)

where $\lambda$ is the Lagrangian multiplier. We then calculate the centroid update by employing the derivative criteria to Eq. (37).

$$\begin{array}{cc}\hfill 0& =\nabla \left({\sum }_{\mathbf{p}\in C}\left(1-\frac{1}{r^2}\mathbf{x}·\mathbf{p}\right)-\lambda {∥\mathbf{x}∥}^2+r^2\right)=-\frac{1}{r^2}\sum _{\mathbf{p}\in C}\mathbf{p}-2\lambda \mathbf{x}\hfill \end{array}$$

(38)

Therefore the following holds:

$$\begin{array}{cc}\hfill \mathbf{x}& =-\frac{1}{2\lambda r^2}\sum _{\mathbf{p}\in C}\mathbf{p}.\hfill \end{array}$$

(39)

Substituting Eq. (39) into the restriction in Eq. (36), we obtain the multiplier $\lambda$ as:

$$\begin{array}{cc}\hfill \left|\lambda \right|& =\frac{1}{2r^3}·∥\sum _{\mathbf{p}\in C}\mathbf{p}∥.\hfill \end{array}$$

(40)

Therefore, the critical point and minimising point ${\mathbf{c}}_{\mathrm{s}}$ is written as

$$\begin{array}{c}\hfill {\mathbf{c}}_{\mathrm{s}}=\frac{r}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}\sum _{\mathbf{p}\in C}\mathbf{p},\end{array}$$

(41)

as claimed.    □

We can observe that Lemma 3 implies that the minimiser obtained by restricting the point to lie on the surface of the sphere is the projection (from the origin) of the minimiser of the Euclidean dissimilarity into the sphere’s surface.

Corollary 1.

Let $C\subset S^n\left(r\right)$ be a finite set, the possible new centroids of C under cosine dissimilarity in ${\mathbb{R}}^3$ are

$$\begin{array}{cc}\hfill {\mathcal{P}}_{\mathrm{s}}\left(C\right):=& \underset{\mathbf{x}\in {\mathbb{R}}^n}{argmin}\sum _{\mathbf{p}\in C}{d}_{\mathrm{s}}(\mathbf{x},\mathbf{p})=\left\{r{\sum }_{\mathbf{p}\in C}\mathbf{p}:r>0\right\}.\hfill \end{array}$$

(42)

We call these thecosine possible new centroids.

Proof. 

We have

$$\begin{array}{cc}\hfill {\mathbf{c}}_{\mathrm{s}}^{\mathrm{update}}:=& \underset{\mathbf{x}\in {\mathbb{R}}^n\setminus \left\{0\right\}}{argmin}\sum _{\mathbf{p}\in C}{d}_{\mathrm{s}}(\mathbf{x},\mathbf{p})\hfill \end{array}$$

(43)

$$\begin{array}{cc}& =\underset{\mathbf{x}\in S^n\left(r\right),r>0}{argmin}\sum _{\mathbf{p}\in C}{d}_{\mathrm{s}}(\mathbf{x},\mathbf{p})\hfill \end{array}$$

(44)

$$\begin{array}{cc}& ={\left\{r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}\right\}}_{r>0}\hfill \end{array}$$

(45)

$$\begin{array}{cc}& =\left\{r{\sum }_{\mathbf{p}\in C}\mathbf{p}:r>0\right\},\hfill \end{array}$$

(46)

where the last equality follows from Eq. (28), namely

$$\begin{array}{c}\hfill d_{\mathrm{s}}\left(r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥},\mathbf{p}\right)=d_{\mathrm{s}}\left({\sum }_{\mathbf{p}\in C}\mathbf{p},\mathbf{p}\right)\end{array}$$

(47)

for all r and $\mathbf{p}$; thus making all the points $r{\sum }_{\mathbf{p}\in C}\mathbf{p}$, $r>0$ equivalent possibilities for the centroid update.    □

2.7. Stereographic Projection

The inverse stereographic projection (ISP), shown in Figure 3, is a bijective mapping

$$\begin{array}{c}\hfill s_r^{-1}:{\mathbb{R}}^n\mapsto S^n\left(r\right)\setminus \left\{N\right\}\end{array}$$

(48)

from the Euclidean space ${\mathbb{R}}^n$ into an n-sphere $S^n\left(r\right)\subset {\mathbb{R}}^{n+1}$ without the north pole N.

This mapping is interesting because of the natural equivalence between the 3D unit sphere $S^2\left(1\right)$ and the Bloch sphere of qubit quantum states. In this case, as displyed in Figure 3, the ISP maps a 2-dimensional point $\mathbf{p}=(p_x,p_y)\in {\mathbb{R}}^2$ into a three-dimensional point $s_r^{-1}\left(\mathbf{p}\right)=(s_x\left(\mathbf{p}\right),s_y\left(\mathbf{p}\right),s_z\left(\mathbf{p}\right))\in S^2\left(r\right)\setminus \left\{(0,0,r)\right\}$ through the following set of transformations:

$$\begin{array}{cccccc}\hfill s_x\left(\mathbf{p}\right)& =p_x·\frac{2r^2}{p_x^2+p_y^2+r^2}\hfill & \hfill s_y\left(\mathbf{p}\right)& =p_y·\frac{2r^2}{p_x^2+p_y^2+r^2}\hfill & \hfill s_z\left(\mathbf{p}\right)& =r·\frac{p_x^2+p_y^2-r^2}{p_x^2+p_y^2+r^2}\hfill \end{array}$$

(49)

$$\begin{array}{cccccc}\hfill & =p_x·\frac{2r^2}{{∥\mathbf{p}∥}^2+r^2}\hfill & \hfill & =p_y·\frac{2r^2}{{∥\mathbf{p}∥}^2+r^2}\hfill & \hfill & =r·\frac{{∥\mathbf{p}∥}^2-r^2}{{∥\mathbf{p}∥}^2+r^2}.\hfill \end{array}$$

(50)

The polar and azimuthal angles of the projected point are given by the expressions:

$$\begin{array}{cccc}\hfill \varphi \left(\mathbf{p}\right)& ={tan}^{-1}\left(\frac{p_y}{p_x}\right)\hfill & \hfill \theta \left(\mathbf{p}\right)& =2·{tan}^{-1}\left(\frac{r}{∥\mathbf{p}∥}\right)\hfill \end{array}$$

(51)

This information, particularly Eq. (51), will allow us to associate each point in ${\mathbb{R}}^2$ to a unique quantum state through the Bloch sphere. Still, the inverse stereographic projection does not need to be bound to the preparation of quantum states and can be used as a transformation between classical kNN algorithms. Indeed, we can stereographically project and then perform classical clustering on the 3D data, namely perform 3DEC-kNN as defined in Definition 4.

Definition 6 (3D Stereographic Classical kNN (3DSC-kNN)). Let $s_r^{-1}$ be an ISP, and let $(D,\overline{\mathbf{c}},{\mathbb{R}}^2,d_{\mathrm{e}})$ be a clustering state (recall, D is the dataset and $\overline{\mathbf{c}}$ are the initial centroids). We then define the 3D Stereographic Classical kNN (3DSC-kNN) as ($s_r^{-1}\left(D\right)$, $s_r^{-1}\left(\overline{\mathbf{c}}\right)$, ${\mathbb{R}}^3$, $d_{\mathrm{e}},{\mathbf{c}}_{\mathrm{e}}^{\mathrm{update}})$.

Here we apply $s_r^{-1}$ elementwise and thus $s_r^{-1}\left(\overline{\mathbf{c}}\right)=(s_r^{-1}\left({\mathbf{c}}_1\right),...,s_r^{-1}\left({\mathbf{c}}_k\right))$ for any list of centroids $\overline{\mathbf{c}}$ and $s_r^{-1}\left(D\right)=\left\{s_r^{-1}\left(\mathbf{p}\right):\mathbf{p}\in C\right\}$ for any set of points C, and where $C_j$ are the clusters as defined in Eq. (22).

Remark 2.

Derivations and further observations can be found in Appendix C. Of particular note, as explained in more detail in Appendix C.2, is that changing the plane’s distance from the centre of the sphere is equivalent to a change of radius. Therefore we can limit our analysis to projections where the centre of the plane is also the centre of the sphere without loss of generality.

3. Stereographic Quantum Nearest-Neighbour Clustering (SQ-kNN)

This section proposes and describes the quantum k (nearest-neighbour) clustering using stereographic embedding. In Section 4, we demonstrate an equivalent quantum-inspired (classical) version. Section 3.1 defines the method to convert the classical data into quantum states. In what follows, we describe how these states are manipulated so that we obtain an output that can be used to perform clustering, using the circuit of Section 2.3 for dissimilarity estimation. Section 3.2 defines the quantum algorithm in terms of , and Section 3.3 discusses the complexity and scalability of the algorithm. Section 3.4 discusses the SQ-kNN algorithm in the context of mixed states.

3.1. Stereographic Embedding, Bloch Embedding and Quantum Dissimilarity

For quantum algorithms to work on classical data, the classical data must be converted into quantum states. This process of encoding classical data into quantum states is also called embedding. The embedding of classical data into quantum states is not unique, and each technique’s pros and cons must be weighed in the context of a specific application. The process of data embedding is an active field of research. More details on existing embedding can be found in Appendix B.

Here, we propose the stereographic embedding as an improved embedding of classical vector $\mathbf{p}\in {\mathbb{R}}^2$ into quantum state using its stereographic projection. We can split stereographic embedding into two steps: inverse stereographic projection and Bloch embedding. We define Bloch embedding, a variation of angle embedding, as follows.

Definition 7

(Bloch embedding). Let $\mathbf{P}\in {\mathbb{R}}^3$. We define theBloch embedded quantum state, orBloch embeddingfor short, of $\mathbf{P}$ as the quantum state

$${\psi }_{\mathbf{P}}:=\frac{1}{2}\left(\mathbb{1}+\frac{\mathbf{P}}{∥\mathbf{P}∥}·\overrightarrow{\sigma }\right)$$

(52)

which is simply the pure state obtained using $\mathbf{P}/∥\mathbf{P}∥$ as Bloch vector.

At this point, we define this general embedding for general 3-dimensional points since this general form will yield the quantum dissimilarity defined next. We will also define this embedding in the context of the ISP in 9 below.

To obtain ${\psi }_{\mathbf{P}}$ the state can be encoded as explained in the preliminaries in Section 2.2, through . For Bloch embedding, the $\theta$ and $\varphi$ of would be the polar and azimuthal angles of $\mathbf{P}$ respectively. We now define the quantum dissimilarity as follows.

Definition 8

(Quantum Dissimilarity). For any two points ${\mathbf{P}}_1,{\mathbf{P}}_2\in {\mathbb{R}}^3$, we define thequantum dissimilarityas

$$d_{\mathrm{q}}({\mathbf{P}}_1,{\mathbf{P}}_2):=\frac{1}{2}(1-Tr\left({\psi }_{{\mathbf{P}}_1}{\psi }_{{\mathbf{P}}_2}\right)),$$

(53)

where ${\psi }_{\mathbf{P}}$ is the Bloch embedding of $\mathbf{P}$.

Notice that, as per this definition, the classical 2-dimensional points are embedded in pure states only. In Section 3.4, we consider Bloch embedding the centroids also into mixed states, showing that it does not give an advantage in our framework. This quantum dissimilarity can be obtained either with the SWAP test or with the Bell state measurement on ${\psi }_{{\mathbf{P}}_1}\otimes {\psi }_{{\mathbf{P}}_2}$ as described in Section 2.3. In our application, we use the Bell state measurement (depicted in Figure 2)), as we do not need the extra resources of the SWAP test that allow us to keep the post-measurement state. For more details, see Lemma 1.

By Eq. (8), the quantum dissimilarity is proportional to the cosine dissimilarity.1

$$d_{\mathrm{q}}({\mathbf{P}}_1,{\mathbf{P}}_2)=\frac{1}{2}(1-Tr\left({\psi }_{{\mathbf{P}}_1}{\psi }_{{\mathbf{P}}_2}\right))=\frac{1}{4}\left(1-\frac{{\mathbf{P}}_1}{∥{\mathbf{P}}_1∥}·\frac{{\mathbf{P}}_2}{∥{\mathbf{P}}_2∥}\right)=\frac{1}{4}{d}_{\mathrm{s}}({\mathbf{P}}_1,{\mathbf{P}}_2)$$

(54)

It is also proportional to the Euclidean distance for points on the same sphere (points with the same magnitude) as per Lemma 2. Namely, if ${\mathbf{s}}_1,{\mathbf{s}}_2\in S^2\left(r\right)$ then:

$$\begin{array}{cc}\hfill d_{\mathrm{q}}({\mathbf{s}}_1,{\mathbf{s}}_2)& =\frac{1}{8r^2}{d}_{\mathrm{e}}({\mathbf{s}}_1,{\mathbf{s}}_2).\hfill \end{array}$$

(55)

We can finally define stereographic embedding as follows.

Definition 9

(Stereographic Embedding). The stereographic embedding of a classical vector $\mathbf{p}\in {\mathbb{R}}^2$ consists of

• projecting the 2D point $\mathbf{p}$ into a point on the sphere of radius r in 3D space through the ISP:

  $$\mathbf{s}:=s_r^{-1}\left(\mathbf{p}\right)\in S^2\left(r\right)\subset {\mathbb{R}}^3;$$

  (56)
• Bloch embedding $\mathbf{s}$ into ${\psi }_{\mathbf{s}}={\psi }_{s_r^{-1}\left(\mathbf{p}\right)}$.

Comparing the distance estimate of the Stereographic embedding procedure (Eqs. (55) and (A34)) with that for the hybrid quantum-classical k-means with angle embedding (the `distance loss function’ described in [10]), we can see that the theoretical performance has been improved, since the estimate has been much improved with respect to the closeness to Euclidean distance. This leads us to expect a performance improvement of the SQ-kNN algorithm over the hybrid quantum-classical implementation of quantum k-means with angle embedding.

A very time-consuming computational step of kNN involves the repeated calculations of distances between the data set points meant to be classified and each centroid. In the case of the quantum kNN in [2], since angle embedding is not injective, many steps must be spent after estimating the fidelity to calculate the distance between the points using the norms. Even in [34,35,39], the norms of the points have to be stored classically, leading to much computational expense. Our method has the clear benefit of calculating the cosine dissimilarity directly through fidelity estimation. No further calculations are required due to all stereographically projected points having the same norm r in the sphere, and the existence of a bijection between the ISP and the original 2D datapoints, thus saving computational time and resources. In summary, Eqs. (54) and (55) portray a method to measure a dissimilarity that leads to consistent clustering involving pure states.

As one can see, in the case of stereographically projected points, $d_{\mathrm{q}}$ is directly proportional to the Euclidean dissimilarity between them. Since all the points after projection into the sphere have equal modulus r, and each projected point corresponds to a unique 2D data point, we can directly compare the probability of getting outcome $ij=11$ on the Bell-state measurement circuit for cluster assignment. This eliminates extra steps needed during computation to account for the different moduli of points on the 2-dimensional plane.

3.2. The SQ-kNN Algorithm

We now have all the building blocks to define the quantum clustering algorithm. The quantum part will be the dissimilarity estimation $d_{\mathrm{q}}$, obtained by embedding the data into quantum states as described in Section 3.1 and then feeding it into the quantum circuit described in to estimate an outcome probability. The finer details of distance estimation are further described in Appendix E. We can now formally define the developed algorithm building on the definition of clustering state (Definition 1), of clustering algorithm and cluster update provided by Definitions 2 and 3, of ISP as defined in Section 2.7 and of quantum dissimilarity $d_{\mathrm{q}}$ from Definition 8.

Definition 10 (Stereographic Quantum kNN(SQ-kNN)). Let $s_r^{-1}$ be the ISP, let $d_{\mathrm{q}}$ be the quantum dissimilarity and let $(D,\overline{\mathbf{c}},{\mathbb{R}}^2,d_{\mathrm{e}})$ be a clustering state (where D and $\overline{\mathbf{c}}$ are 2-dimensional dataset and initial centroids). We then define the Stereographic Quantum kNN (SQ-kNN) as the kNN clustering algorithm

$$\begin{array}{c}\hfill \left(s_r^{-1}\left(D\right),s_r^{-1}\left(\overline{\mathbf{c}}\right),{\mathbb{R}}^3,d_{\mathrm{q}},{\overline{\mathbf{c}}}_{\mathrm{q}}^{\mathrm{update}}\right)\end{array}$$

(57)

where ${\mathbf{c}}_{\mathrm{q}}^{\mathrm{update}}:={\sum }_{\mathbf{p}\in C}\mathbf{p}$.

The complete process of performing SQ-kNN in practice can be described in detail as follows.

• First, prepare to embed the classical data and initial centroids into quantum states using the ISP: project the 2-dimensional datapoints and initial centroids (in our case, the alphabet) into a sphere of radius r and calculate the polar and azimuthal angles of the points. This first step is executed entirely on a classical computer.
• Cluster Update: The calculated angles are used to create the states using Bloch embedding (Definition 7). The dissimilarity between the centroid and point is then estimated using the Bell-state measurement. Once the dissimilarities between a point and all the centroids have been obtained, the point is assigned to the cluster of the `closest’ centroid. This is repeated for all the points that have to be classified. The quantum circuit and classical controller handle this step entirely. The controller feeds in the classical values at the appropriate times, stores the results of the various shots and classifies the point to the appropriate cluster.
• Centroid Update: Since any non-zero point on the subspace of ${\mathbf{c}}_{\mathbf{s}}$ (see Corollary 1, Figure 4) is an equivalent choice, to minimise computational expense, the centroids are updated as the sum point of all points in the cluster - as opposed to the average, for example, which minimises the Euclidean dissimilarity (Eq. (26)).

Once the centroids are updated, Step 2 (Cluster Update) is repeated, followed once again by Step 3 (Centroid Update) until a decided stopping condition is fulfilled.

Compared to 2D quantum kNN clustering with angle or amplitude embedding, the differences with the SQ-kNN algorithm lie in the embedding and the post-processing after the inner-product estimation.

• The stereographic embedding of the 2D datapoints is done by inverse stereographic projecting the point into a sphere of a chosen radius and then producing the quantum state obtained by rescaling the sphere to radius one.
  In contrast, in angle embedding, the coefficients of the vectors are used as the angles of the Bloch vector (also known as dense angle embedding [47]), while in amplitude embedding, they are used as the amplitudes in the standard basis. For 2D vectors, amplitude embedding allows one to encode only one coefficient (instead of two) in one qubit, and sometimes angle embedding would also encode only one coefficient by using a single rotation (standard angle embedding [48]). Both angle and amplitude embeddings require the lengths of the vectors to be stored classically beside the quantum state, which is not needed in Bloch embedding.
• No post-processing is needed after the overlap estimation of stereographically embedded data as the obtained estimate is already a linear function of the inner product, as opposed to standard approaches using angle or amplitude encoding. Amplitude embedding also requires non-trivial computational time in the state preparation process. In contrast, in angle embedding, though the state preparation time is constant, recovering a useful dissimilarity (e.g. Euclidean) may involve many post-processing steps.

In short, stereographic embedding has the advantage of angle over amplitude embedding of being able to encode all values of a vector and low state preparation time, and the advantage of amplitude versus angle embedding in the recovery of the dissimilarity.

3.3. Complexity Analysis and Scaling

Let the ISP of a d-dimensional point $\mathbf{p}=[x_1{x}_2\dots x_d]$ into $S^d\left(r\right)$, using the point $(r,0,0,\dots ,0)$ (the `North pole’) as the projection point, be the point $\mathbf{s}=[s_0{s}_1\dots s_d]$. It is known that the Cartesian coordinates of $\mathbf{s}$ are given by:

$$\begin{array}{cccc}\hfill s_0& =r\frac{{∥p∥}^2-r^2}{{∥p∥}^2+r^2}\hfill & \hfill s_i& =\frac{2r^2{x}_i}{{∥p∥}^2+r^2}\forall i=1,\cdots ,d\hfill \end{array}$$

(58)

where ${∥\mathbf{p}∥}^2={\sum }_{j=1}^d{x}_j^2$. Hence one can see that the time complexity of projection for a single $d-$dimensional point is $O\left(d\right)$, provided we only need the Cartesian coordinates. However, for the Stereographic Embedding procedure, one would need to calculate the angles made by $\mathbf{s}$ with the axes of the space, making the time complexity of projection $O\left(poly\right(d\left)\right)$. Therefore, the total time complexity of the Stereographic Embedding for a $d-$dimensional dataset D of size $\left|D\right|=N$ and k centroids is given by $O\left(\right(k+N\left)poly\right(d\left)\right)$. We now specify two strategies for scaling our algorithm for higher dimensional datapoints.

3.3.1. Using Qubit-Based System

We consider the case where we have two d-dimensional vectors ${\mathbf{p}}_1,{\mathbf{p}}_2\in {\mathbb{R}}^d$, and we want to compute the quantum dissimilarity of vectors. If we have a qubit-based system and use dense angle encoding for encoding the stereographically projected point, we would encode the d calculated angles using $\frac{d}{2}$ qubits. Namely, for the $(d+1)$-dimensional projection of a d-dimensional point ${\mathbf{p}}_1$, one would get d angles $[{\theta }_1{\theta }_2\dots {\theta }_d]$ that specify the projected point ${\mathbf{s}}_1$ on $S^d\left(r\right)$. We then encode this vector using the same unitary as follows:

$$\begin{array}{c}\hfill |{\psi }_1\rangle =\underset{j\in (1,3,\dots ,d-1)}{⨂}|{\psi }_{1_j}\rangle =\underset{j\in odd\left(d\right)}{⨂}U({\theta }_j,{\theta }_{j+1},0)|0\rangle .\end{array}$$

(59)

If d is odd we can pad $[{\theta }_1{\theta }_2\dots {\theta }_d]$ with an extra 0 to make it even. The other point ${\mathbf{s}}_2=s_r^{-1}\left({\mathbf{p}}_2\right)$ will be encoded into the state

$$\begin{array}{c}\hfill |{\psi }_2\rangle =\underset{j^{\prime }\in (1,3,\dots ,d-1)}{⨂}|{\psi }_{2_{j^{\prime }}}\rangle =\underset{j^{\prime }\in odd\left(d\right)}{⨂}U({\theta }_{j^{\prime }}^{\prime },{\theta }_{j^{\prime }+1}^{\prime },0)|0\rangle .\end{array}$$

(60)

Now, to find the overlap between the states, one would have to perform the Bell-state measurement (Section 2.3) pairwise using $|{\psi }_{1_j}\rangle$ and $|{\psi }_{2_{j^{\prime }}}\rangle$ as inputs, i.e.

$$\begin{array}{c}\hfill {|\langle {\varphi }_{11}||{\psi }_{1_j}\rangle \otimes |{\psi }_{2_{j^{\prime }}}\rangle |}^2=\frac{1}{2}{(1-|\langle {\psi }_{1_j}|{\psi }_{2_{j^{\prime }}}\rangle |}^2)\end{array}$$

(61)

In the common practical case of the vectors being expressed in an orthogonal basis, one would only have to find the overlap for $j=j^{\prime }$. We would then have the quantum dissimilarity by adding up the individual probabilities

$$\begin{array}{c}\hfill d_{\mathrm{q}}({\mathbf{p}}_1,{\mathbf{p}}_2)=\sum _{j\in odd\left(d\right)}{|\langle {\varphi }_{11}||{\psi }_{1_j}\rangle \otimes |{\psi }_{2_j}\rangle |}^2\end{array}$$

(62)

This procedure has a time complexity of $O\left(d\right)$. It is important point to note that this quantum dissimilarity will no longer correspond directly to the inner product or Euclidean distance between either ${\mathbf{p}}_1$ and ${\mathbf{p}}_2$ or ${\mathbf{s}}_1$ and ${\mathbf{s}}_2$. With the strategy of pairwise overlap estimation, we see that if the number of shots to estimate $\left|\langle {\varphi }_{11}|\right(|{\psi }_{1_j}\rangle \otimes |{\psi }_{2_j}\rangle {|}^2$ is kept constant, the error in estimation will be $\propto d$. Hence, taking into account the increase in the number of shots to estimate the quantum dissimilarity with a given total error $ϵ$, the time complexity of this qubit implementation of overlap estimation between two points using SQ-kNN scales as $O\left({ϵ}^{-1}poly\left(d\right)\right)$. Hence for all points and clusters, the time complexity would be $O\left({ϵ}^{-1}kNpoly\left(d\right)\right)$.

It is shown in [49] that collective measurements are a better strategy than repeated individual measurements for overlap estimation. Although this is shown in [49] for estimating overlap between two states given the availability of multiple copies of the same states, similar collective measurement strategies could be applied in this case for better results. In conclusion, the time complexity of SQ-kNN for qubit-based implementation is

$$\begin{array}{cc}\hfill & O\left({ϵ}^{-1}kNpoly\left(d\right)\right).\hfill \end{array}$$

(63)

3.3.2. Using Qudit-Based System

Consider

$$|\mathbb{1}\rangle :=\sum _{i\in \left\{0,...,d-1\right\}}|ii\rangle .$$

Then, for any two real vectors $|\psi \rangle =\sum {\psi }_i|i\rangle$ and $|\varphi \rangle =\sum {\varphi }_i|i\rangle$, namely when ${\psi }_i,{\varphi }_i\in \mathbb{R}$ we have

$$\begin{array}{c}\hfill \langle \mathbb{1}|(|\psi \rangle \otimes |\varphi \rangle )=\sum {\psi }_i{\varphi }_i=\langle {\varphi }^{*}|\psi \rangle =\langle \varphi |\psi \rangle .\end{array}$$

Now, if we make a qudit Bell measurement2, one of the basis states of this von Neumann measurement will be

$$|\Phi \rangle \equiv {|\Phi \rangle }_d=\frac{1}{\sqrt{d}}|\mathbb{1}\rangle .$$

Thus, the inner product between two real vectors can still be measured with a Bell measurement, but the resulting probability of measuring outcome $|\Phi \rangle$ scales as

$${\left|\langle \Phi |(|\psi \rangle \otimes |\varphi \rangle )\right|}^2=\frac{1}{d}{\left|\langle \varphi |\psi \rangle \right|}^2;$$

meaning that, as the inner product remains constant going to higher dimensions, the number of shots needed to estimate the inner product with constant precision scales polynomially in the dimension. In contrast, such complexity for the SWAP test remains constant because the contribution of the fidelity to the outcome probability is not divided by the dimension d. This is why the SWAP test is usually considered for inner product estimation, even if in the case of qubits the Bell measurement is a simpler solution [50].

3.4. SQ-kNN and Mixed States

Instead of estimating the quantum dissimilarity, we can use the datapoints produced by the ISP to perform classical kNN clustering on the 3D projected data. We called this the 3DSC-kNN (3D Stereographic Classical kNN) in 6. This algorithm produces centroids that are inside the sphere. As previously pointed out, when computing the Euclidean 3D centroid on the data projected on the sphere, the result is a point inside the sphere rather than on the sphere itself.

In the Bloch sphere, internal points are mixed states, namely states with added noise. In contrast, the quantum algorithm (SQ-kNN) always produces pure centroids, namely points on the surface of the sphere. The only noiseless states are the pure states on the surface of the sphere, and thus the intuition is that arguably mixed states should not help. However, this is not immediately clear from the algorithm. Comparing 3DEC-kNN to SQ-kNN, it is thus natural to ask whether embedding the centroids into mixed states inside the Bloch sphere improves the accuracy.

Here, we show that the intuition is correct, namely that projecting into the pure state centroid is a better option. The reason is that while the quantum dissimilarity is proportional to the Euclidean dissimilarity for states in the same sphere, the same is not true for Bloch vectors with different lengths.

To allow for mixed state embedding, we can modify the definition of quantum dissimilarity (Eq. (55)) to produce mixed states whenever the 3D vector has a length of less than one. This results in the following new quantum dissimilarity.

Definition 11

(Noisy Quantum Dissimilarity). Let $B^2\left(1\right)=\left\{\mathbf{P}\in {\mathbb{R}}^3|∥\mathbf{P}∥\le 1\right\}$ be the ball of radius 1. We define thenoisy quantum dissimilarityas the function ${\tilde{d}}_{\mathrm{q}}:B^2\left(1\right)\times B^2\left(1\right)\to \mathbb{R}$

$$\begin{array}{cc}\hfill {\tilde{d}}_{\mathrm{q}}({\mathbf{P}}_1,{\mathbf{P}}_2):=& \frac{1}{2}\left(1-Tr\left({\rho }_{{\mathbf{P}}_1}{\rho }_{{\mathbf{P}}_2}\right)\right)\hfill \end{array}$$

(64)

where ${\rho }_{\mathbf{P}}$ is the quantum state of the Bloch vector $\mathbf{P}$ as in Eq. (7).

Now suppose we have a convex combination of pure states; namely, suppose we have

$$\begin{array}{cccc}\hfill \overline{\rho }& =\sum _i{p}_i{\rho }_{{\mathbf{P}}_i}\hfill & \hfill {∥\mathbf{P}∥}_i& =1\forall i.\hfill \end{array}$$

(65)

where $p_i$ is a probability distribution ($p_i>0$ and $\sum p_i=1$). By linearity, we have

$$\begin{array}{cccc}\hfill \overline{\rho }& =\rho \left(\sum p_i{\mathbf{P}}_i\right)\rho \left(\overline{\mathbf{P}}\right)\hfill & \hfill \overline{\mathbf{P}}& =\sum p_i{\mathbf{P}}_i,\hfill \end{array}$$

(66)

By convexity, $\overline{\mathbf{P}}$ will always lie in the sphere. Namely, we have $∥\overline{\mathbf{P}}∥\le 1$ and thus by linearity

$$\begin{array}{cc}\hfill {\tilde{d}}_{\mathrm{q}}\left(\overline{\mathbf{P}},\mathbf{P}\right)& =\sum p_i{\tilde{d}}_{\mathrm{q}}\left({\mathbf{P}}_i,\mathbf{P}\right)=\sum p_i{d}_{\mathrm{q}}\left({\mathbf{P}}_i,\mathbf{P}\right).\hfill \end{array}$$

(67)

The result in Eq. (67) can be interpreted as another two-step process: first, repeatedly performing the Bell-state measurement of each state ${\rho }_{{\mathbf{P}}_i}$ that makes up the cluster and ${\rho }_{\mathbf{P}}$ corresponding to the datapoint, to estimate each individual dissimilarity; and then, taking the weighted average of the dissimilarities according to the composition of the mixed state centroid. This procedure is clearly impractical experimentally and also no longer correlates to the cosine dissimilarity for mixed states.

Computing the diagonalization of $\overline{\rho }$ as per Eq. (11)

$$\begin{array}{cccc}\hfill \rho & =p\rho \left(\frac{\mathbf{P}}{∥\mathbf{P}∥}\right)+(1-p)\rho \left(\frac{-\mathbf{P}}{∥\mathbf{P}∥}\right)\hfill & \hfill p& =\frac{1}{2}(1+∥\overline{\mathbf{P}}∥)\hfill \end{array}$$

(68)

$$\begin{array}{cc}& =p{\psi }_{\overline{\mathbf{P}}}+(1-p){\psi }_{-\overline{\mathbf{P}}}\hfill \end{array}$$

(69)

(where $\psi$ is the Bloch embedding) makes the estimation more practical by reducing it to two estimations of $d_{\mathrm{q}}$, namely

$$\begin{array}{cc}\hfill {\tilde{d}}_{\mathrm{q}}\left(\overline{\mathbf{P}},\mathbf{P}\right)& =p{\tilde{d}}_{\mathrm{q}}\left(\frac{\overline{\mathbf{P}}}{∥\overline{\mathbf{P}}∥},\mathbf{P}\right)+(1-p){\tilde{d}}_{\mathrm{q}}\left(-\frac{\overline{\mathbf{P}}}{∥\overline{\mathbf{P}}∥},\mathbf{P}\right)\hfill \end{array}$$

(70)

(71)

The implementation portrayed at Eq. (71) simplifies the measurement procedure of the mixed state. Furthermore, instead of estimating $d_{\mathrm{q}}(\pm \overline{\mathbf{P}},\mathbf{P})$ separately, the estimation can be done directly by preparing $\psi \left(\mathbf{P}\right)$ with probability p and $\psi (-\mathbf{P})$ with probability $1-p$, and finally collecting all the outcomes in a single estimation, which requires a larger number of shots to achieve the same precision of estimation. Another issue is that the points $\overline{\mathbf{P}},-\overline{\mathbf{P}}$ have to be computed, which is quite time-consuming. This is true even for Eq. (67); however, a number of shots proportional to the number of Bloch vectors ${\mathbf{P}}_i$ in the cluster is needed for an accurate estimation. Regardless, linearity and convexity make it clear that using mixed states can only increase the quantum dissimilarity.

Namely, while in Euclidean dissimilarity, points inside the sphere can reduce the dissimilarity, the quantum dissimilarity is proportional to the Euclidean dissimilarity only for unit vectors and actually increases for points inside the Bloch sphere. Hence we conclude that the behaviour of 3DSC-kNN does not carry over to SQ-kNN.

4. Quantum-Inspired Stereographic $\mathbf{k}$-Nearest-Neighbour Clustering

We have detailed the developed quantum algorithm in the previous Section 3. This section develops the classical analogue to this quantum algorithm - the `quantum-inspired’ classical algorithm. A table summarising all the algorithms discussed in this paper, including the next one, can be found in Table 1. We begin by defining this analogous classical algorithm in terms of the clustering state (1), deriving a relationship between the Euclidean and spherical centroids given datapoints that lie on a sphere, and then proving our claim that the defined classical algorithm and previously described stereographic quantum kNN are indeed equivalent.

Recall from Lemma 3 that

$$\begin{array}{cc}\hfill {\mathbf{c}}_{\mathrm{s}}^{\mathrm{update}}\left(C\right):=& \underset{\mathbf{x}\in S^n\left(r\right)}{argmin}\sum _{\mathbf{p}\in C}{d}_{\mathrm{s}}(\mathbf{x},\mathbf{p})=r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}.\hfill \end{array}$$

(72)

Definition 12 (2D Stereographic Classical kNN(2DSC-kNN)). Let $s_r^{-1}$ be the ISP, and let $(D,\overline{\mathbf{c}},{\mathbb{R}}^2,d_{\mathrm{e}})$ be a 2D euclidean clustering state . We define the 2D Stereographic Classical kNN (2DSC-kNN) as

$$\begin{array}{c}\hfill \left(s_r^{-1}\left(D\right),s_r^{-1}\left(\overline{\mathbf{c}}\right),S^2\left(r\right),d_{\mathrm{s}},{\overline{\mathbf{c}}}_{\mathrm{s}}^{\mathrm{update}}\right).\end{array}$$

(73)

Remark 3.

Notice that due to the cluster update being cosine (${\overline{\mathbf{c}}}_{\mathrm{s}}^{\mathrm{update}}$) and 3, we can equivalently substitute $d_{\mathrm{s}}$ with $d_{\mathrm{e}}$, namely we can substitute it without changing the outcome of the cluster update. In our implementation, we use the Euclidean dissimilarity for simplicity of coding.

To expand upon Definition 12, for the quantum-inspired/classical analogue stereographic kNN, the steps of execution are as follows:

• Stereographically project all the 2-dimensional data and initial centroids into the sphere $S^2\left(r\right)$ of radius r. Notice that the initial centroids will lie on the sphere by construction.
• Cluster Update: Form the clusters using the method defined in Eq. (22), i.e. form all $C_j\left({\overline{\mathbf{c}}}_i\right)$. Here, $\mathbb{D}=S^2\left(r\right)$ and dissimilarity $d=d_{\mathrm{s}}(\mathbf{p},\mathbf{c})=\frac{1}{2r^2}{d}_{\mathrm{e}}(\mathbf{p},\mathbf{c})$ (Definition 5 and Lemma 2).
• Centroid Update: A closed-form expression for the centroid update was calculated in Eq. (41) $\left({\mathbf{c}}_{\mathrm{s}}^{\mathrm{updated}}=\frac{r}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}{\sum }_{\mathbf{p}\in C}\mathbf{p}\right)$. This expression recalculates the centroid once the new clusters have been formed. Once all the centroids are updated, Step 2 (cluster update) is then repeated, and so on, until a stopping condition is met.

Table 1. Summary of various kNN algorithms, where SC stands for “stereographic classical” (2D or 3D) and SQ for “stereographic quantum”. Here, D is the 2-dimensional dataset, $\overline{\mathbf{c}}$ are the 2-dimensional initial centroids (initial transmission points), $s_r^{-1}$ is the ISP into $S^2\left(r\right)$, and the dissimilarities $d_{\mathrm{e}}$, $d_{\mathrm{s}}$ and $d_{\mathrm{q}}$ are defined in Eq. (25) and Definitions 5 and 8, respectively. The option of using $d_{\mathrm{e}}$ instead of $d_{\mathrm{s}}$ in the 2DSC-kNN is due to Remark 3.

Algorithm	Reference	Dataset	Initial	Dataspace	Dissimilarity	Centroid
			Centroids			Update
		D	${\overline{\mathbf{c}}}_1$	$\mathbb{D}$	d	${\mathbf{c}}^{\mathrm{update}}$
2DEC-kNN	Definition 4	D	$\overline{\mathbf{c}}$	${\mathbb{R}}^2$	$d_{\mathrm{e}}$	$\frac{1}{\left|C\right|}{\sum }_{\mathbf{p}\in C}\mathbf{p}$
3DSC-kNN	Definition 6	$s_r^{-1}\left(D\right)$	$s_r^{-1}\left(\overline{\mathbf{c}}\right)$	${\mathbb{R}}^3$	$d_{\mathrm{e}}$	$\frac{1}{\left|C\right|}{\sum }_{\mathbf{p}\in C}\mathbf{p}$
2DSC-kNN	Definition 12	$s_r^{-1}\left(D\right)$	$s_r^{-1}\left(\overline{\mathbf{c}}\right)$	$S^2\left(r\right)$	$d_{\mathrm{s}}$ (or $d_{\mathrm{e}}$)	$r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}$
SQ-kNN	Definition 10	$s_r^{-1}\left(D\right)$	$s_r^{-1}\left(\overline{\mathbf{c}}\right)$	${\mathbb{R}}^3$	$d_{\mathrm{q}}$	${\sum }_{\mathbf{p}\in C}\mathbf{p}$

Table 1. Summary of various kNN algorithms, where SC stands for “stereographic classical” (2D or 3D) and SQ for “stereographic quantum”. Here, D is the 2-dimensional dataset, $\overline{\mathbf{c}}$ are the 2-dimensional initial centroids (initial transmission points), $s_r^{-1}$ is the ISP into $S^2\left(r\right)$, and the dissimilarities $d_{\mathrm{e}}$, $d_{\mathrm{s}}$ and $d_{\mathrm{q}}$ are defined in Eq. (25) and Definitions 5 and 8, respectively. The option of using $d_{\mathrm{e}}$ instead of $d_{\mathrm{s}}$ in the 2DSC-kNN is due to Remark 3.

Algorithm	Reference	Dataset	Initial	Dataspace	Dissimilarity	Centroid
			Centroids			Update
		D	${\overline{\mathbf{c}}}_1$	$\mathbb{D}$	d	${\mathbf{c}}^{\mathrm{update}}$
2DEC-kNN	Definition 4	D	$\overline{\mathbf{c}}$	${\mathbb{R}}^2$	$d_{\mathrm{e}}$	$\frac{1}{\left|C\right|}{\sum }_{\mathbf{p}\in C}\mathbf{p}$
3DSC-kNN	Definition 6	$s_r^{-1}\left(D\right)$	$s_r^{-1}\left(\overline{\mathbf{c}}\right)$	${\mathbb{R}}^3$	$d_{\mathrm{e}}$	$\frac{1}{\left|C\right|}{\sum }_{\mathbf{p}\in C}\mathbf{p}$
2DSC-kNN	Definition 12	$s_r^{-1}\left(D\right)$	$s_r^{-1}\left(\overline{\mathbf{c}}\right)$	$S^2\left(r\right)$	$d_{\mathrm{s}}$ (or $d_{\mathrm{e}}$)	$r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}$
SQ-kNN	Definition 10	$s_r^{-1}\left(D\right)$	$s_r^{-1}\left(\overline{\mathbf{c}}\right)$	${\mathbb{R}}^3$	$d_{\mathrm{q}}$	${\sum }_{\mathbf{p}\in C}\mathbf{p}$

4.1. Equivalence

We now want to show that the 2DSC-kNN algorithm of 12 is equivalent to the previously defined quantum algorithm using stereographic embedding (10). For that, we first define the equivalence of two clustering algorithms.

Definition 13

(Equivalence of Clustering Algorithms). Let $\mathcal{K}=(D,{\overline{\mathbf{c}}}_1,\mathbb{D},d,{\overline{\mathbf{c}}}_{\mathrm{update}})$ and ${\mathcal{K}}^{\prime }=(D^{\prime },{\overline{\mathbf{c}}}_1^{\prime },{\mathbb{D}}^{\prime },d^{\prime },{\overline{\mathbf{c}}}_{\mathrm{update}}^{\prime })$ be two clustering algorithms. They are said to be equivalent if there exists a transformation $t:\mathbb{D}\to {\mathbb{D}}^{\prime }$ such that it maps the data, initial centroids and centroid update, and clusters of $\mathcal{K}$ to the data, initial centroids and centroid update, and clusters of ${\mathcal{K}}^{\prime }$; namely if

• $D^{\prime }=t\left(D\right)$,
• ${\overline{\mathbf{c}}}_1^{\prime }=t\left({\overline{\mathbf{c}}}_1\right)$ and ${\overline{\mathbf{c}}}_{\mathrm{update}}^{\prime }=t\circ {\overline{\mathbf{c}}}_{\mathrm{update}}$,
• $C_j^{\prime }\left(t\left(\overline{\mathbf{c}}\right)\right)=t\left(C_j\left(\overline{\mathbf{c}}\right)\right)$ for all $j=1,...,k$ and any $\overline{\mathbf{c}}\in {\mathbb{D}}^k$.

where we apply t elementwise and thus $t\left(\overline{\mathbf{c}}\right)=(t\left({\mathbf{c}}_1\right),...,t\left({\mathbf{c}}_k\right))$ for any list of centroids $\overline{\mathbf{c}}$ and $t\left(C\right)=\left\{t\left(\mathbf{p}\right):\mathbf{p}\in C\right\}$ for any set of datapoints C, and where $C_j$ are the clusters as defined in Eq. (22).

Theorem 1.

SQ-kNN (Definition 10) and 2DSC-kNN (Definition 12) are equivalent.

Proof. 

By definition, let $(D,{\overline{\mathbf{c}}}_1,{\mathbb{R}}^2,d_{\mathrm{e}})$ be the 2D clustering state and thus giving us the SQ-kNN algorithm as

$$\begin{array}{c}\hfill \mathcal{K}=\left(S,{\overline{\mathbf{s}}}_1,{\mathbb{R}}^3,d_{\mathrm{q}},{\overline{\mathbf{c}}}_{\mathrm{q}}^{\mathrm{update}}\right)\end{array}$$

(74)

and the 2DSC-kNN clustering algorithm as

$$\begin{array}{c}\hfill {\mathcal{K}}^{\prime }=\left(S,{\overline{\mathbf{s}}}_1,S^2\left(r\right),d_{\mathrm{s}},{\overline{\mathbf{c}}}_{\mathrm{s}}^{\mathrm{update}}\right)\end{array}$$

(75)

where

$$\begin{array}{cccc}\hfill S:=& s_r^{-1}\left(D\right)\hfill & \hfill {\overline{\mathbf{s}}}_1:=& s_r^{-1}\left({\overline{\mathbf{c}}}_1\right)\hfill \end{array}$$

(76)

Let us use the notation $\widehat{\mathbf{p}}:=\frac{\mathbf{p}}{∥\mathbf{p}∥}$ and define the transform $t:{\mathbb{R}}^3\mapsto S^2\left(r\right)$ as $t\left(\mathbf{p}\right)=r\widehat{\mathbf{p}}$, which rescales any vector to have length r. Observe that trivially for all $\mathbf{p}\in S^2\left(r\right)$, $t\left(\mathbf{p}\right)=\mathbf{p}$ and thus $t\circ s_r^{-1}=s_r^{-1}$. Therefore

$$\begin{array}{cccc}\hfill t\left(S\right)& =S,\hfill & \hfill t\left({\overline{\mathbf{s}}}_1\right)& ={\overline{\mathbf{s}}}_1.\hfill \end{array}$$

(77)

Also, the equivalence of centroids is obtained since

$$\begin{array}{c}\hfill t\left({\mathbf{c}}_{\mathrm{q}}^{\mathrm{update}}\left(C\right)\right)=t\left(\sum _{\mathbf{p}\in C}\mathbf{p}\right)=r\frac{{\sum }_{\mathbf{p}\in C}\mathbf{p}}{∥{\sum }_{\mathbf{p}\in C}\mathbf{p}∥}={\mathbf{c}}_{\mathrm{s}}^{\mathrm{update}}\left(C\right).\end{array}$$

(78)

For the clusters, we prove the equivalence of the cluster updates as follows. We will use $d_{\mathrm{s}}(\mathbf{a},\mathbf{b})=4·d_{\mathrm{q}}(\mathbf{a},\mathbf{b})$ (54) and the fact that $d_{\mathrm{s}}$ and $d_{\mathrm{q}}$ are invariant under t, namely $d_{\mathrm{s}}\circ t=d_{\mathrm{s}}$ and $d_{\mathrm{q}}\circ t=d_{\mathrm{q}}$, or more explicitly

$$\begin{array}{cccc}\hfill d_{\mathrm{s}}(t\left(\mathbf{a}\right),\mathbf{b})& =d_{\mathrm{s}}(\mathbf{a},\mathbf{b})=d_{\mathrm{s}}(\mathbf{a},t\left(\mathbf{b}\right)),\hfill & \hfill d_{\mathrm{q}}(t\left(\mathbf{a}\right),\mathbf{b})& =d_{\mathrm{q}}(\mathbf{a},\mathbf{b})=d_{\mathrm{q}}(\mathbf{a},t\left(\mathbf{b}\right)).\hfill \end{array}$$

(79)

Let now $\overline{\mathbf{s}}\in {\left({\mathbb{R}}^3\right)}^k$. Then, using the above equations and that $t\left(\mathbf{s}\right)=\mathbf{s},t\left(S\right)=S$, we have

$$\begin{array}{cc}\hfill C_j^{\prime }\left(t\left(\overline{\mathbf{s}}\right)\right)& =\left\{\mathbf{p}\in S|d_{\mathrm{s}}(\mathbf{p},t\left({\mathbf{s}}_j\right))\le d_{\mathrm{s}}(\mathbf{p},t\left({\mathbf{s}}_{\ell }\right))\forall \ell =1,...,k,\mathbf{p}\notin \bigcup _{\ell <j}{C}_{\ell }^{\prime }\left(t\left(\overline{\mathbf{s}}\right)\right)\right\}\hfill \end{array}$$

(80)

$$\begin{array}{cc}& =\left\{\mathbf{p}\in S|d_{\mathrm{q}}(\mathbf{p},{\mathbf{s}}_j)\le d_{\mathrm{q}}(\mathbf{p},{\mathbf{s}}_{\ell })\forall \ell =1,...,k,\mathbf{p}\notin \bigcup _{\ell <j}{C}_{\ell }\left(\overline{\mathbf{s}}\right)\right\}\hfill \end{array}$$

(81)

where the change in the dissimilarity inequality has also transformed the calculation of $C_{\ell }^{\prime }\left(\overline{\mathbf{s}}\right)$ into the calculation of $C_{\ell }\left(\overline{\mathbf{s}}\right)$. We are now done, since $t\left(\mathbf{s}\right)=\mathbf{s}$ for $\mathbf{s}\in S$ and thus

$$\begin{array}{cc}\hfill 1-1C_j^{\prime }\left(t\left(\overline{\mathbf{s}}\right)\right)& =\left\{t\left(\mathbf{p}\right)\in t\left(S\right)|d_{\mathrm{q}}(\mathbf{p},{\mathbf{s}}_j)\le d_{\mathrm{q}}(\mathbf{p},{\mathbf{s}}_{\ell })\forall \ell =1,...,k,\mathbf{p}\notin \bigcup _{\ell <j}{C}_{\ell }\left(\overline{\mathbf{s}}\right)\right\}\hfill \end{array}$$

(82)

(83)

This concludes the proof    □

The following discussion provides a visual intuition of Theorem 1. In Figure 4, the sphere with centre origin (O) and radius r is the stereographic sphere into which the 2-dimensional points are projected, while the sphere with centre O and radius 1 is the Bloch sphere. The points ${\mathbf{p}}_1,{\mathbf{p}}_2,\dots ,{\mathbf{p}}_n$ are the stereographically projected points defining a cluster, corresponding to the previously used labels ${\mathbf{s}}_1,{\mathbf{s}}_2,\dots ,{\mathbf{s}}_n$. The centroid ${\mathbf{c}}_{\mathrm{e}}$ is obtained with the euclidean average in ${\mathbb{R}}^3$. In contrast, the centroid ${\mathbf{c}}_{\mathrm{s}}$ is restricted to be in $S^2\left(r\right)$ and equal ${\mathbf{c}}_{\mathrm{e}}$ rescaled to lie on this sphere. The quantum states $|{\psi }_{{\mathbf{p}}_1}\rangle ,|{\psi }_{{\mathbf{p}}_2}\rangle ,\dots ,|{\psi }_{{\mathbf{p}}_n}\rangle$ are obtained after Bloch embedding the stereographically projected points ${\mathbf{p}}_1,{\mathbf{p}}_2,\dots ,{\mathbf{p}}_n$, and $|{\psi }_{\mathbf{c}}\rangle$ is the quantum state obtained after Bloch embedding the centroid. The points marked on the Bloch sphere in Figure 4 are the Bloch vectors of the quantum states $|{\psi }_{{\mathbf{p}}_1}\rangle ,|{\psi }_{{\mathbf{p}}_2}\rangle ,\dots ,|{\psi }_{{\mathbf{p}}_n}\rangle$ and $|{\psi }_{\mathbf{c}}\rangle$.

One can see from Definition 7 that the origin, any point $\mathbf{p}$ on the sphere and $|{\psi }_{\mathbf{p}}\rangle$ are collinear. Hence, it can be seen that in the process of SQ-kNN clustering, the points on the stereographic sphere are projected radially into the sphere of radius 1. Once the labels were assigned in the previous iteration, the new centroid is computed, giving an integer multiple of the average point ${\mathbf{c}}_{\mathrm{e}}$, which lies within the stereographic sphere. Crucially, when we embed this new centroid into the quantum state for the quantum dissimilarity calculation of the next step, since we only use the polar and azimuthal angle of the point for embedding (see Definition 7), the prepared quantum state is also projected into the surface of the Bloch sphere - or, in other words, a pure state is prepared ($|{\psi }_{\mathbf{c}}\rangle$). Hence, we can see that all the dissimilarity calculations in the SQ-kNN will take place between points on the surface of the Bloch sphere, even though the calculated quantum centroid is contained outside the stereographic sphere. This argument also illustrates why any point on the ray $\overrightarrow{\mathbf{O}{\mathbf{c}}_{\mathrm{e}}{\mathbf{c}}_{\mathrm{s}}}$ can be used for the centroid update step of the stereographic quantum kNN; any chosen point on the ray, when embedded into a quantum state for dissimilarity calculations will reduce to $|{\psi }_c\rangle$.

In short, we know from Lemma 3 that $O,{\mathbf{c}}_{\mathrm{e}},$ and ${\mathbf{c}}_{\mathrm{s}}$ lie on a straight line. Therefore one can see that if the Bloch sphere is scaled by r, the point on the Bloch sphere corresponding to $|{\psi }_{\mathbf{c}}\rangle$ will transform to ${\mathbf{c}}_{\mathrm{s}}$, i.e. $0,|{\psi }_{\mathbf{c}}\rangle ,{\mathbf{c}}_{\mathrm{e}}$ and ${\mathbf{c}}_{\mathrm{s}}$ are all collinear. Eq. (55) shows that SQ-kNN clustering clusters points on a sphere as per Euclidean dissimilarity; that implies that simply scaling the sphere makes no difference to the clustering. Therefore, we conclude that clustering on the surface of the stereographic sphere $S^2$ (2DSC-kNN) is equivalent to the quantum algorithm with stereographic embedding (SQ-kNN).

4.2. Complexity Analysis and Scaling

As we showed in Section 3.2, the time complexity of ISP for calculating Cartesian coordinates of a $d-$dimensional vector is $O\left(d\right)$. Hence the total time complexity of projection for the 2DSC-kNN will be $O\left(\right(k+N\left)d\right)$ where $N=\left|D\right|$ is the total number of points, and k is the total number of centroids. Since the cluster update step uses Euclidean dissimilarity, it will take $O\left(kNd\right)$ time in total ($O\left(d\right)$ for each distance calculation, which is to be done for each pair of N points and k centroids). The centroid update expression (41) can be calculated in $O\left(Nd\right)$ making the total time for this step $O\left(kNd\right)$ since we have k centroids. Hence we have

$$\begin{array}{c}\hfill \mathrm{Time}\mathrm{complexity}\mathrm{of}2\mathrm{DSC}-\mathrm{kNN}\mathrm{algorithm}=O\left(kNd\right),\end{array}$$

(84)

on par with the classical k-means clustering algorithm, and at least polynomially faster than the stereographic quantum kNN (63) or Lloyd’s quantum clustering algorithm (taking into account input assumptions) [2,9,13].

5. Experiments and Results

We defined the procedure for SQ-kNN in Section 3. Section 3.1 introduces our idea for state preparation - projecting the 2-dimensional data points into a higher dimension. Section 3.1 details the hybrid quantum-classical method used for our process and then proves that the output of the quantum circuit is not only a valid but also an excellent metric that can be used for distance estimation between two points. Section 4 describes the quantum-inspired classical algorithm analogous to the quantum algorithm (2DSC-kNN). In this section, we test and compare the quantum-inspired rather the quantum algorithm for two main reasons:

• The hardware performance and availability of quantum computers (NISQ devices) is currently so much worse than that of classical computers that no advantage can likely be obtained with the quantum algorithm.
• The goal of this paper is not to show a “quantum advantage” in time complexity over the classical k-means in the NISQ context - it is to show that stereographic projection can lead to better learning for classical clustering and be a better embedding for quantum clustering. In particular, the equivalence between 2DSC-KNN and SQ-KNN proves that noise is the only limitation for the stereographic quantum algorithm to achieve the accuracy of the quantum-inspired algorithm.

All the experiments were carried out on a server with the following specifications: 2 Intel Xeon E5-2687W v4 chips clocked at 3.0 GHz (24 cores / 48 threads), 128GB RAM. All experiments are performed on the real-world 64-QAM data provided by Huawei (see Section 2.1 and Appendix A). Due to the extensive nature of testing and the large volume of analysis generated, we do not present all the figures in the following sections. Figures which sufficiently demonstrate general trends and observations have been included here. An exhaustive collection of all figures and other such analysis results, as well as the source code, real-world data, and raw data collected, can be obtained from the GitHub repository of the project.

The terminology used is as follows:

• Radius: the radius of the stereographic sphere into which the 2-dimensional points are projected.
• Number of points: the number of points upon which the clustering algorithm was performed. For every experiment, the selected points were a random subset of all the 64-QAM data (of a specific noise) with cardinality equal to the required number of points. The random subset is created using the numpy.random.sample() from the Python Numpy library.
• Number of runs: Since for each choice of parameters for each experiment we select a subset of points at random, we repeat each of the experiments many times to remove bias from the random choice and obtain stable averages and standard deviations for the collected performance parameters (described in another list below). This number of repetitions is the “number of runs”.
• Dataset Noise: As explained in Section 2.1, data was collected for four different input powers. The data is divided into four datasets labelled with powers 2.7, 6.6, 8.6 and 10.7 dBm.
• Natural endpoint: The natural endpoint of a clustering algorithm occurs when

  $$\begin{array}{cccc}\hfill C_j\left({\overline{\mathbf{c}}}_{i+1}\right)& =C_j\left({\overline{\mathbf{c}}}_i\right)\hfill & \hfill \forall & j=1,...,k\hfill \end{array}$$

  (85)

  i.e. when all the clusters remained unchanged (stay the same) even after the centroid update. It is the natural endpoint since if the clusters do not change, the centroids will not change either in the next iteration, leading to the same clusters (Eq. (85)) and centroids for all future iterations.

The algorithms that we test are:

• 2DSC-kNN: The quantum-analogue algorithm of Definition 12, the classical equivalent of the SQ-KNN and the most important candidate for our testing.
• 2DEC-kNN: The standard classical kNN of Definition 4 implemented upon the original 2-dimensional dataset ($n=2$), which serves as a baseline for performance comparison.
• 3DSC-kNN: The standard classical kNN, but implemented upon the stereographically projected 2-dimensional dataset, as defined in Definition 6. We again emphasise that in contrast to the 2DSC-kNN the centroid lies within the sphere, and in contrast to the 2DEC-kNN, the clustering takes place in ${\mathbb{R}}^3$. This algorithm serves as another control, to gauge the relative impacts of stereographically projecting the dataset versus restricting the centroid to the surface of the sphere. It is an intermediate step between the 2DSC-kNN and the 2DEC-kNN algorithms.

From these algorithms, we measure the following performance parameters (or KPIs, Key Performance Indicators):

• Accuracy: Since we have the true labels of the datapoints available, we can measure the accuracy of the algorithm as the percentage of points that have been given the correct label, i.e. symbol accuracy rate. All accuracies are recorded as a percentage.
• Symbol or Bit error rate: As mentioned in Appendix A, due to Gray encoding, the bit error rate is approximately $\frac{1}{6}$ of the symbol error rate, which in turn is simply one minus the accuracy. Although error rates are the standard performance parameter in channel coding, we decided to measure the accuracy instead, which is the standard performance parameter in machine learning.
• Accuracy gain: The gain is calculated as (accuracy of candidate algorithm - accuracy of 2-dimensional classical k-means clustering algorithm), i.e. it is the increase in accuracy of the algorithm over the baseline, defined as the accuracy of the classical k-means clustering algorithm acting on the 2D dataset for those number of points.
• Number of iterations: One iteration of the clustering algorithm occurs when the algorithm performs the cluster update followed by the centroid update (the algorithm must then perform the cluster update again). The number of times the algorithm repeats these two steps before stopping is the number of iterations. We use the number of iterations the algorithm requires to reach its `natural endpoint’ as a proxy for convergence performance. The lesser the number of iterations performed, the faster the algorithm’s convergence. The number of iterations does not directly correspond to time performance since the time taken for one iteration differs between all algorithms.
• Iteration gain: The gain in iterations is defined as (the number of iterations of 2D k-means clustering algorithm - the number of iterations of candidate algorithm), i.e. the gain is how many fewer iterations the candidate algorithm took than the 2DEC-kNN algorithm to reach its natural endpoint.
• Execution time: The amount of time taken for a clustering algorithm to give the final output (the final centroids and clusters) given the 2-dimensional data points as input, i.e. the time taken end to end for the clustering process. All times in this work are recorded in milliseconds (ms).
• Execution time gain: This gain is calculated as (the execution time of 2DEC-kNN k-means clustering algorithm - the execution time of candidate algorithm).
• Overfitting Parameter: (accuracy in testing−accuracy in training).

With these algorithms and variables, we perform two main experiments:

• The Overfitting Test: The dataset is divided into a `training’ and a `testing’ set, to characterise the clustering and classification performance of the algorithms.
• The Stopping Criterion Test: The iterations and other performance parameters are varied, to test whether and what kind of stopping criterion is required.

We see that the tested algorithms display some very promising and interesting results. We manage to get improvements in accuracy and convergence performance almost across the board, and we discover the very important optimisation parameters of the radius of projection and the stopping criterion.

5.1. Experiment 1: Overfitting

Here, the datasets were divided into training and testing data. First, a random subset of cardinality equal to the number of points was chosen from the dataset, and then 80% of the selected points were assigned as `Training Data’, while the other 20% was assigned as `Testing Data’.

In the training phase, the algorithms were first run on the training data with the maximum possible iterations set to 50 to keep an acceptable running time. The stopping criterion for all algorithms was chosen as the natural endpoint - the algorithm stopped either when the number of iterations hit 50, or when the natural endpoint was reached (whichever happened first). The final centroid co-ordinates (${\overline{\mathbf{c}}}_{\mathrm{last}\mathrm{iteration}}$) were recorded in the training phase, to be used for the testing phase, along with several performance parameters. The recorded performance parameters were the algorithm’s accuracy, the number of iterations taken and the execution time.

Once the training was over, the centroids calculated at the end of training were used as the initial centroids for the testing set datapoints, and the algorithm was run with the maximum number of iterations set to 1, i.e. the calculated centroids were then used to classify the remaining points as per the dissimilarity and dataspace of each algorithm. The recorded performance parameters were the algorithm’s accuracy and execution time. Once both the testing and training accuracy had been recorded, the overfitting parameter was also recorded.

For each set of input variables (just the number of points for 2DEC-kNN clustering, the radius and number of points for the 2DSC-kNN and 3DSC-kNN clustering), the entire experiment (training and testing) was repeated 10,000 times in batches of 100 to calculate reasonable standard deviations for every performance parameter.

There are several reasons for this choice of experiment:

• It exhaustively covers all the parameters that can be used to quantify the performance of the algorithms. We were able to observe very important trends in the performance parameters with respect to the choice of radius and the effect of the number of points (affecting the choice of when one should trigger the clustering process on the collected received points).
• It avoids the commonly known problem of overfitting. Though this approach is not usually used in testing the kNN due to its iterative nature, we felt that from a machine learning perspective, it is useful to know how well the algorithms perform in a classification setting as well.
• Another reason that justifies the training and testing approach (clustering and classification) is the nature of the real-world application setup. When transmitting QAM data through optical-fibre, the receiver receives only one point at a time and has to classify the received point to a given cluster in real-time using the current centroid values. Once a number of data points have accumulated, the kNN algorithm can be run to update the centroid values; after the update, the receiver will once again perform classification until some number of points has been accumulated. Hence, we can see that in this scenario, the clustering and the classification performance of the chosen method become important.

5.1.1. Results

We begin the presentation of the results of this experiment by first showing the characterisation of the 2DSC-kNN algorithm with respect to the input variables.

Figure 5 characterises the testing and training accuracy of the 2DSC-kNN algorithm acting upon the 2.7dBm dataset, i.e. classification and clustering performance respectively. Figure 6 portrays the same results in the form of a heat map, with a focus on the region of interest of the algorithm. These figures are representative of the trends of all four datasets.

Figure 7 characterises the convergence performance of the quantum algorithm - it shows how the number of iterations required to reach the natural endpoint of the 2DSC-kNN algorithm varies as the number of points and radius of projection changes. Once again, the figures for all the other datasets follow the same pattern as the included figures.

We then compare the performance of the 2DSC-kNN algorithm with that of the 3DSC-kNN and 2DEC-kNN algorithms.

Accuracy Performance

In all the figures in this section, the winner is chosen as the radius for which the maximum accuracy is achieved for the given number of points.

Figure 8 depicts the variation in testing accuracy with the number of points for all three algorithms along with error bars. As mentioned before, this characterises the performance of the algorithms in `classification’ mode, that is, when the received points must be decoded in real-time.

Figure 9 portrays the trend in training accuracy with the number of points for all three algorithms along with error bars. This characterises the performance of the algorithms in `clustering’ mode, that is, when the received points must be used to update the centroid for future re-classification or if the received datapoints are stored and decoded in batches.

Figure 8 and Figure 9 also plot the gain in testing and training accuracies respectively for the 3DSC-kNN and 2DSC-kNN algorithms. The label of the points in these figures is the radius of ISP for which that accuracy gain was achieved.

Iteration Performance

In all the figures in this section, the winner is chosen as the radius for which the minimum number of iterations is achieved for the given number of points.

Figure 10 shows how the required number of iterations for all three algorithms varies as the number of points increases. Figure 10 also displays the gain of the 2DSC-kNN and 3DSC-kNN algorithms in the number of iterations to reach their natural endpoints. The label of the points in these figures is the radius of ISP for achieving that iteration gain.

Time Performance

In all the figures in this section, the winner is chosen as the radius for which the minimum execution time is achieved for the given number of points.

Figure 11 puts forth the dependence of testing execution time upon the number of points for all three algorithms along with error bars. As mentioned before, these times are effectively the amount of time the algorithm takes for one iteration. This characterises the performance of the algorithms when performing direct classification decoding of the received points in real time.

This figure reveals the trend in training execution time with the number of points for all three algorithms, along with error bars. This characterises the time performance of the algorithms when performing clustering - that is, when the received points must be used to update the centroid for future re-classification or if the received datapoints are stored and decoded in batches.

Figure 12 graphs the gains in testing and training execution times for the 3DSC-kNN and 2DSC-kNN algorithms.

Overfitting Performance

Figure 13 exhibits how the overfitting parameter for the 2DEC-kNN k-means clustering, 3DSC-kNN and 2DSC-kNN algorithms vary as the number of points changes.

5.1.2. Discussion and Analysis

From Figure 5 we can see that there is an ideal radius > 1 for which maximum accuracy is achieved. This ideal radius is usually between 2 and 5 for our datasets. For a good choice of radius (>1), the accuracy increases monotonically (with an upper bound) with the number of points. In contrast, for a poor choice of radius (<1), the accuracy nosedives as the number of points increases. This is due to the clusters getting squished together near the North pole of the stereographic sphere (the point $(0,0,r)$). If one is dealing with a large number of points, the accuracy becomes even more sensitive to the choice of radius as the decline in accuracy for a bad radius is much steeper as the number of points increases. These observations hold for both training and testing accuracy (classification and clustering), regardless of the noise in the dataset. These observations are also well-reflected in the heatmaps, where one can see that the best training and testing performance is for $r=2$ to 3 and the maximum number of points. It would seem that choosing too large of a radius is not too harmful. This might hold true for the classical algorithms, but when the quantum algorithm is deployed, all the points will be clustered around the South pole of the Bloch sphere and even minimal noise in the quantum circuit will degrade performance. Hence, there is a sweet spot of radius to be chosen.

Figure 7 also shows that there is an ideal radius > 1 for which one needs the minimum number of iterations to reach the natural endpoint. This ideal radius is once again between 2 and 5 for our datasets. As the number of points increases, the number of iterations always increases. The increase is minimal for a good choice of radius, while for a bad choice, the convergence is very slow. For our experiments, we chose the maximum iterations as 50 hence the observed plateau at 50 iterations. If one is dealing with a large number of points, the convergence becomes more sensitive to the choice of radius. The increase in iterations for a poor choice of radius is much steeper. 2DSC-kNN algorithm and 3DSC-kNN algorithm display near-identical performance.

From Figure 8 we can see that both 2DSC-kNN algorithm and 3DSC-kNN algorithm perform better in accuracy than the 2DEC-kNN algorithm for all datasets. The advantage becomes more definitive as the number of points increases as the increase in accuracy moves beyond the error bar. We observe the highest increase in accuracy for the 2.7dBm dataset.

In Figure 9 one can see the noticeably better performance of the 2DSC-kNN algorithm and 3DSC-kNN algorithm over the 2DEC-kNN algorithm for all datasets than in the testing case (classification mode). Once again the 2.7dBm dataset shows the maximum increase. The advantage again becomes more definitive as the number of points increases as the increase in accuracy moves beyond the error bar. The 2DSC-kNN algorithm and 3DSC-kNN algorithm show an almost identical performance.

From Figure 8 and Figure 9 we can also see that almost universally for both algorithms, the gain is greater than 0, i.e. we beat the 2DEC-kNN algorithm in nearly every case! We can also see that the best radius is almost always between 2 and 5. Another observation is that the gain in training accuracy increases with the number of points. The figures further display how similarly the 3DSC-kNN algorithm and 2DSC-kNN algorithm perform in terms of accuracy, regardless of noise.

From Figure 10 it can be concluded that for low noise datasets, since the number of iterations is already quite low, there is not much gain or loss; all three algorithms perform almost identically. For high-noise datasets, however, both the 3DSC-kNN algorithm and 2DSC-kNN algorithm show significant performance improvement, especially for a higher number of points. For a high number of points, the improvement is beyond the error bars and hence very significant. It can be noticed that the ideal radius for minimum iterations is once again between 2 and 5. Here also the 3DSC-kNN algorithm and 2DSC-kNN algorithm perform similarly, with the 2DSC-kNN algorithm performing better in certain cases.

One learns from Figure 11 that most importantly, the 2DSC-kNN algorithm and 2DEC-kNN algorithm take nearly the same amount of time for execution in classification mode, and the 2DSC-kNN algorithm in most cases beats the 3DSC-kNN algorithm. Here too, the gain is significant since it is much beyond the error bar. The execution time increases linearly with the number of points, as expected. These conclusions are supported by Figure 12.

Since the 2DSC-kNN algorithm takes almost the same time, and provides greater accuracy, it is an ideal candidate to replace the 2DEC-kNN algorithm for classification applications.

Figure 11 also shows that all three algorithms take almost the same amount of time for training, i.e. in clustering mode. The 3DSC-kNN and 2DSC-kNN algorithms once again perform almost identically, almost always slightly worse than 2DEC-kNN clustering. Figure 12 supports these observations. Here execution time increases linearly with the number of points as well, as expected.

In Figure 13, all three algorithms have nearly identical performance. As expected. the overfitting decreases with an increase in the number of points.

5.2. Experiment 2: Stopping criterion

Based on the results obtained from the first experiment, we performed another experiment to see how the accuracy of the algorithms varies iteration by iteration. It was observed that the natural endpoint of the algorithm was rarely the ideal endpoint in terms of performance. Hence we wished to observe the performance of each algorithm as the number of iterations progressed.

In this experiment, the entire random subset of datapoints was used for the clustering algorithm. The algorithms were run on the dataset, and the accuracy of the algorithms at each iteration as well as the iteration number of the natural endpoint was recorded. The maximum number of iterations was once again 50. By repeating this 100 times for each number of points (and radius, if applicable), we obtained the general performance variation of each algorithm with the iteration number. The input variables were the number of points, the radius of the stereographic sphere and the iteration number; the recorded performance parameters were the accuracy and probability of stopping.

This experiment revealed that the natural endpoint was indeed a poor choice of stopping criterion and that the endpoint should be chosen per some “loss function”. It also revealed some important trends in the performance parameters which not only emphasised the importance of the choice of radius and number of points but also gave greater insight into the disadvantages and advantages of each algorithm.

5.2.1. Results

Characterisation of the 2DSC-kNN algorithm

Figure 14 depicts the dependence of the accuracy of the 2DSC-kNN algorithm upon the iteration number and projection radius for the 2.7dBm dataset. The figures for the rest of the datasets follow the same trends and are nearly identical in shape.

Figure 15 shows the dependence of the probability of the 2DSC-kNN algorithm reaching its natural endpoint versus the radius of projection and iteration number for the 10.7dBm dataset with 51200 points and for the 2.7dBm dataset with 640 points. Once again, the figures for the rest of the datasets follow the same trends and their shape can be extrapolated from the presented Figure 15.

Comparison with 2DEC-kNN and 3DS-kNN K-Means Clustering

Figure 16 portrays the gain of the 2DSC-kNN and 3DSC-kNN algorithms in the number of iterations to reach maximum accuracy for the 2.7 and 10.7dBm datasets. In these figures, a gain of `g’ means that the algorithm took `g’ fewer iterations than the classical k means acting upon the 2D dataset did to reach maximum accuracy.

Figure 17 plots the gain of the 2DSC-kNN and 3DSC-kNN algorithms in the maximum achieved accuracy for the 2.7 and 10.7dBm datasets. Here, a gain of `g’ means that the algorithm was $g\%$ more accurate than the maximum accuracy of the classical k means acting upon the 2D dataset.

Lastly, Figure 18 illustrates the maximum accuracies achieved by the 2DSC-kNN, 3DSC-kNN, and 2DEC-kNN algorithms for the 2.7 and 10.7dBm datasets.

5.2.2. Discussion and Analysis

Figure 14 shows that once again, there is an ideal radius for which maximum accuracy is achieved. The ideal projection radius is larger than one, in particular, it seems to be between two and five. Most importantly, there is an ideal number of iterations for maximum accuracy, beyond which the accuracy reduces. As the number of points increases, the sensitivity of the accuracy to radius increases significantly. For a bad choice of radius, accuracy only falls with an increase in the number of iterations and stabilises at a very low value. For a good radius, accuracy increases to a point as iterations proceed, and then stabilises at a slightly lower value. If the allowed number of iterations is restricted, the choice of radius to achieve the best results becomes extremely important. With a good radius one can achieve nearly the maximum possible accuracy with very few iterations. As mentioned before, this holds for all dataset noises. As the dataset noise increases, the iteration number at which maximum accuracy is achieved also expectedly increases. Since accuracy always falls after a point, choosing a stopping criterion is essential rather than waiting for the algorithm to reach its natural endpoint. An idea for the stopping criterion is to record the sum of the average dissimilarity for each centroid at each iteration and stop the algorithm if that quantity increases.

Figure 15 portrays that for a good choice of radius, the 2DSC-kNN algorithm approaches convergence much faster. For $r<1$ the algorithm converges much slower or never converges. As the number of points increases, the convergence rate for poor radius falls dramatically. For a radius greater than the ideal radius as well, the convergence rate is lower. As one would expect, the algorithm takes longer to converge as dataset noise increases. As mentioned before, if the number of iterations is severely limited, the choice of radius becomes very important. The algorithm can reach its ideal endpoint in very few iterations if the radius is chosen well.

Through Figure 16 we see that for lower values of noise, both algorithms do not produce much advantage in terms of iteration gain, regardless of the number of points in the data set. However, both algorithms significantly outperform the classical one at higher noise in the dataset and a high number of points. This effect is especially significant for the 2DSC-kNN algorithm. For the highest noise and all the points, it saves over 20 iterations compared to the 2DEC-kNN algorithm - an advantage of over 50%. One of the reasons for this is that at low noises, the algorithms already perform quite well, and it is at high noise with a high number of points that the algorithm is stressed enough to reveal the difference in performance. It should be noted that these gains are much higher than when the algorithms are allowed to reach their natural endpoint, suggesting another reason for choosing an ideal stopping criterion.

Figure 17 shows that for all datasets and numbers of points, the two algorithms perform better than 2DEC-kNN clustering. The 3DSC-kNN algorithm and 2DSC-kNN algorithms perform nearly the same, and the accuracy gain seems to stabilise with an increase in the number of points. Figure 18 supports these conclusions.

5.3. Overall Observations

Overall Observations from Experiment 1

• The ideal projection radius is greater than 1 and between 2 and 5. At this ideal radius, one achieves maximum testing and training accuracy, and minimum iterations.
• In general, the accuracy performance is the same for 3DSC-kNN and 2DSC-kNN algorithms - this shows a significant contribution of the ISP to the advantage as opposed to `quantumness’. This is a significant distinction, not made by any previous work.
• 2DSC-kNN and 3DSC-kNN algorithms lead to an increase in the accuracy performance in general, with the increase most pronounced for the 2.7dBm dataset.
• The 2DSC-kNN algorithm and 3DSC-kNN algorithm give more iteration performance gain (fewer iterations required than 2DEC-kNN) for high noise datasets and for a large number of points.
• Generally, increasing the number of points favours the 2DSC-kNN and 3DSC-kNN algorithms, with the caveat that a good radius must be carefully chosen.

Overall Observations from Experiment 2

• These results further stress the importance of choosing a good radius (2 to 5 in this application) and a better stopping criterion. The natural endpoint is not suitable.
• The results justify the fact that the developed 2DSC-kNN algorithm has significant advantages over 2DEC-kNN k-means clustering and 3DSC-kNN clustering.
• The 2DSC-kNN algorithm performs nearly the same as the 3DSC-kNN algorithm in terms of accuracy, but for iterations to achieve this max accuracy, the 2DSC-kNN algorithm is better (especially for high noise and a high number of points).
• The developed 2DSC-kNN algorithm and 3DSC-kNN algorithm are better than the 2DEC-kNN algorithm in general - in terms of accuracy and iterations to reach that maximum accuracy.
• The supremacy of the 2DSC-kNN algorithm over the 2DEC-kNN algorithm implies that a fully-quantum SQ-kNN algorithm would have an advantage over the fully-quantum k-means algorithm of [2].

6. Conclusion and Further work

This work considers the practical case of performing kNN on experimentally acquired 64-QAM data. This work has described the problem in detail and explained how the SQ-kNN and its classical analogue, the 2DSC-kNN clustering algorithm, can be used. The proposed processes and circuits, as well as the theoretical justification for the SQ-kNN quantum algorithm and the 2DSC-kNN classical algorithm, have been described in detail. Finally, the simulation results on the real-world datasets have been presented, along with relevant analysis. From the analysis, one can see that the classical analogue of the stereographic quantum kNN, the 2DSC-kNN algorithm, is something that should be considered for industrial implementation - the experiments provide proof of concept. It also shows the importance of choosing the projection radius and provides a very useful embedding for quantum machine learning algorithms - the generalised stereographic embedding. The theoretical advantage offered by the SQ-kNN algorithm over the hybrid quantum-classical quantum k-means algorithm has also been demonstrated. Another important inference from the obtained results is that the SQ-kNN algorithm offers a way to achieve the same advantage compared to the fully quantum k-means that 2DSC-kNN has over 2DEC-kNN - by using the stereographically projected quantum states. These results warrant the practical implementation and testing of both quantum and classical algorithms.

Quantum and quantum-inspired computing has the potential to change the way certain algorithms are performed, with potentially significant advantages. However, as the field is still in relative infancy, finding where quantum and quantum-inspired computing fits in practice is a challenging problem. Here, we have seen that quantum and quantum-inspired computing can indeed be applied to signal-processing scenarios and could potentially work well in the noisy quantum era as clustering algorithms that are relatively robust to noise and inaccuracy.

6.1. Future Work

One of the most important directions of future work is to experiment with more diverse datasets. More experimentation may also lead to more sophisticated methods of selecting the radius for ISP. A more detailed analysis of how to choose a radius of projection through analytical methods is another important direction for future work. A differential geometric analysis of the effects of the ISP on a square grid gives a rough intuition of why one needs an appropriate radius. A comparison with amplitude embedding is also warranted. The ellipsoidal projection (Appendix D) is another promising and novel idea that is to be explored further. In this project, two different stopping criteria for the algorithm were proposed and revealed a change in its performance; yet there is plenty of room to explore more possible stopping criteria.

Further directions of study include improved overlap estimation methods [49] and communication scenarios where the dimensionality of the data points is greatly increased. For example, when multiple carriers experience identical or at least systematically correlated phase rotations.

Another future work is to benchmark against sampling-based quantum-inspired algorithms. As part of a research analysis to evaluate the best possibilities for achieving a practical speed-up, we investigated the landscape of classical algorithms inspired by the sampling in quantum algorithms. Initially, we found that such algorithms have theoretical complexity competing with quantum algorithms, however only under arguably unrealistic assumptions on the structure of the classical data. As the performance of the quantum algorithms turns out to be extremely poor, this reopens the possibility that quantum-inspired algorithms can yield performance improvements while we wait for quantum computers with sufficiently low noise. Thus future work will also be a practical implementation of the quantum-inspired kNN [13], with the goal of testing the computational advantage over 2DEC-kNN, 3DSC-kNN, and 2DSC-kNN algorithms.

Appendix A.1. Description of 64-QAM Data

The various datasets we collected through the setup described in Section 2.1 are visualised in this section. As mentioned, there are four datasets with launch powers of $2.7,6.6,8.6,10.7$ dBm, corresponding to various noise levels. Each data set consists of three variables:

• `alphabet’: The initial analog values at which the data was transmitted, in the form of complex numbers, i.e. for an entry ($\mathrm{a}+i\mathrm{b}$), the transmitted signal was of the form $asin\left(\theta \right)+bcos\left(\theta \right)$. Since the transmission protocol is 64-QAM, there are 64 values in this variable. The transmission alphabet is the same irrespective of the non-linear distortions.
• `rxsignal’: The received analog values of the signal by the receiver. This data is in the form of a $52124\times 5$ matrix. Each datapoint was transmitted five times to the receiver, and so each row contains the values detected by the receiver during the different instances of the transmission of the same datapoint. The values in different rows represent unique datapoint values detected by the receiver.
• `bits’: This is the true label for the transmitted points. This data is in the form of a $52124\times 6$ matrix. Since the protocol is 64-QAM, each analog point represents six bits. These six bits are the entries in each column, and each value in a different row represents the correct label for a unique transmitted datapoint value. The first three bits encode the column, and the last three bits encode the row - see Figure A3.

In Appendix A.1, we can see the analog transmission values (alphabet) for all channels. In the subsequent four figures (Figure A2 to Figure A2), we can observe the transmission data for all the iterations for each channel. The first transmission data is represented as blue crosses, the second transmission as orange circles, the third transmission as yellow dots, the fourth transmission as purple stars, and the fifth transmission as green pluses. The de-mapping alphabet is depicted in Figure A3.

Figure A1. The analog alphabet (initial analog transmission values) for the data transmitted in all the channels. The real part represents the amplitude of the transmitted sine wave, and the imaginary part represents the amplitude of the cosine wave.

Figure A1. The analog alphabet (initial analog transmission values) for the data transmitted in all the channels. The real part represents the amplitude of the transmitted sine wave, and the imaginary part represents the amplitude of the cosine wave.

One can see from these figures that as the noise in the channel increases, the points are further scattered away from the initial alphabet. In addition, the non-linear noise effects also increase, causing distortion of the `shape’ of the data, most clearly visible in Figure A2 - especially near the `corners’. The birefringence phase noise also increases with an increase in the channel noise, causing all the points to be `rotated’ about the origin.

Once the centroids have been found and the data has been clustered, as mentioned before, we need to `de-map’ the analog centroid values and clusters to bit-strings. For this, we need a de-mapping alphabet which maps the analog values of the alphabet to the corresponding bit strings. The de-mapping alphabet is depicted in Figure A3. It can be seen from the figure that, as in most cases, the points are Gray coded, i.e. adjacent points differ in binary translation by only 1 bit. This helps minimise the number of bit errors per symbol error in case of misclassification or exceptionally high noise. In case a point is misclassified, with most probability, it will be assigned to a neighbouring cluster. Since the surrounding clusters differ by only 1 bit, it minimises the bit error rate. Due to Gray coding, the bit error rate is approximately $\frac{1}{6}$ of the symbol error rate.

Figure A2. The four datasets detected by the receiver with various launch powers corresponding to different noise levels: 2.7dBm (top left), 6.6dBm (top right), 8.6dBm (bottom left), 10.7 dBm (bottom right). All five iterations of transmission are depicted together.

Figure A2. The four datasets detected by the receiver with various launch powers corresponding to different noise levels: 2.7dBm (top left), 6.6dBm (top right), 8.6dBm (bottom left), 10.7 dBm (bottom right). All five iterations of transmission are depicted together.

Figure A3. The bit-string mapping and demapping alphabet.

Figure A3. The bit-string mapping and demapping alphabet.

Appendix B.1. Angle Embedding

Angle encoding [48,54,55] is one of the most fundamental forms of encoding classical data into a quantum state. Each data point is represented as a separate qubit. The $n^{th}$ classical real number is encoded into the rotation angle of the $n^{th}$ qubit. In its most basic form, this encoding requires N qubits to represent N dimensional data. It is quite cheap to prepare in terms of complexity – all that is needed is one rotation quantum gate for each qubit. This is one of the forms of encoding we have used to implement quantum kNN clustering. It is generally useful for quantum neural networks and other such QML algorithms. Angle encoding encodes N features into the rotation angles of n qubits where $N\le n$.

The rotations can be chosen as either $R_X\left(\theta \right)$, $R_Y\left(\theta \right)$ or $R_Z\left(\theta \right)$ gates. As a first step, each input data point is normalised to the interval $[0,\pi ]$. To encode the data points, a rotation around the y-axis is used. The angle of rotation depends on the value of the normalised data point. This creates the following separable state:

$$\begin{array}{c}\hfill |\psi \rangle =R_Y\left(x_0\right)|0\rangle \otimes R_Y\left(x_1\right)|0\rangle \otimes \dots \otimes R_Y\left(x_n\right)|0\rangle \end{array}$$

(A3)

$$\begin{array}{c}\hfill =\left(\begin{array}{c}cos{x}_0\\ sin{x}_0\end{array}\right)\otimes \left(\begin{array}{c}cos{x}_1\\ sin{x}_1\end{array}\right)\otimes \dots \otimes \left(\begin{array}{c}cos{x}_n\\ sin{x}_n\end{array}\right)\end{array}$$

(A4)

It can easily be seen that one qubit is needed per data point, which is not optimal. To load the data, the rotations on the qubits can be performed in parallel; thus, the depth of the circuit is optimal [48].

The main advantage of this encoding is that it is very efficient in terms of operations – only a constant number of parallel operations are needed regardless of how many data values need to be encoded. This is not optimal from a qubit point of view (the circuit is very wide), as every input vector component requires one qubit. Another related encoding, dense angle encoding, exploits an additional property of qubits (relative phase) to use only $n/2$ qubits to encode n data points. QRAM can be used to generate the more compact quantum state $\sum |i\rangle R_Y\left({\theta }_i\right)|0\rangle$

Appendix C.1. ISP for General Radius

In this appendix, the transformations for obtaining the Cartesian coordinates of the projected point on a sphere of general radius are derived, followed by the derivation of polar and azimuthal angles of the point on the sphere. First mentioned are three conditions that the point on the sphere must satisfy, and then follows the rest of the derivation. Refer to Figure A4 for a better understanding of the conditions and calculations.

Figure A4. ISP for a sphere of radius `r’. In this figure, the plane is the plane of angle $\varphi$ with respect to the X-axis perpendicular to the XY plane, i.e. the plane containing the 2D point and it projection.

Figure A4. ISP for a sphere of radius `r’. In this figure, the plane is the plane of angle $\varphi$ with respect to the X-axis perpendicular to the XY plane, i.e. the plane containing the 2D point and it projection.

• Azimuthal angle of the original point and the projected point must be the same, i.e. the original point, projected point, and the top of the sphere (the point from which all projections are drawn) lie on the same plane, which is perpendicular to the 2D plane.

  $$\begin{array}{cc}\hfill \Rightarrow \frac{s_y}{s_x}& =\frac{p_y}{p_x}\hfill \end{array}$$

  (A5)
• The projected point lies on the sphere.

  $$\begin{array}{cc}\hfill \Rightarrow s_x^2+s_y^2+s_z^2& =r^2\hfill \end{array}$$

  (A6)
• The triangle with vertices $(0,0,r)$, $(0,0,0)$ and $(p_x,p_y,0)$ is similar to the triangle with vertices $(0,0,r)$, $(0,0,s_z)$ and $(s_x,s_y,s_z)$:

  $$\begin{array}{cc}\hfill \Rightarrow \frac{\sqrt{p_x^2+p_y^2}}{r}& =\frac{\sqrt{s_x^2+s_y^2}}{r-s_z}.\hfill \end{array}$$

  (A7)

Using Eqs. (A7) and (A5), we get

$$\begin{array}{cc}\hfill s_x=p_x·\left(1-\frac{s_z}{r}\right)& \mathrm{and}{s}_y=p_y·\left(1-\frac{s_z}{r}\right)\hfill \end{array}$$

Substituting in Eq. (A6) we get,

$$\begin{array}{cc}\hfill s_z& =r\left(\frac{p_x^2+p_y^2-r^2}{p_x^2+p_y^2+r^2}\right)\hfill \end{array}$$

Hence one gets the set of transformations:

$$\begin{array}{cc}\hfill s_x& =p_x\left(\frac{2r^2}{p_x^2+p_y^2+r^2}\right)\hfill \end{array}$$

(A8)

$$\begin{array}{cc}\hfill s_y& =p_y\left(\frac{2r^2}{p_x^2+p_y^2+r^2}\right)\hfill \end{array}$$

(A9)

$$\begin{array}{cc}\hfill s_z& =r\left(\frac{p_x^2+p_y^2-r^2}{p_x^2+p_y^2+r^2}\right)\hfill \end{array}$$

(A10)

Calculating the polar angle ($\theta$) and azimuthal angle ($\varphi$):

$$\begin{array}{cc}\hfill \varphi & ={tan}^{-1}\left(\frac{s_y}{s_x}\right)\hfill \end{array}$$

(A11)

$$\begin{array}{cc}\hfill \Rightarrow \varphi & ={tan}^{-1}\left(\frac{p_y}{p_x}\right)\hfill \end{array}$$

(A12)

$$\begin{array}{cc}\hfill tan\left(\frac{\pi -\theta }{2}\right)& =\frac{\sqrt{p_x^2+p_y^2}}{r}\hfill \end{array}$$

(A13)

$$\begin{array}{cc}\hfill \Rightarrow \theta & =2·{tan}^{-1}\left(\frac{r}{\sqrt{p_x^2+p_y^2}}\right)\hfill \end{array}$$

(A14)

Appendix C.2. Equivalence of Displacement and Scaling

Refer to Figure A5. Here, $\mathbf{T}$ is the point from which all the projections originate; $\mathbf{O}$ and ${\mathbf{O}}^{\prime }$ are the centres of the projection spheres of radius $\mathrm{r}$ and $(1+\delta )\mathrm{r}$ respectively; $\mathbf{P}$ is the point on the 2D plane to be projected; and $\mathbf{S}$ and ${\mathbf{S}}^{\prime }$ are the ISP of P on the spheres of radius $\mathrm{r}$ and centre $\mathbf{O}$ and radius ${\mathrm{r}}^{\prime }$ and centre ${\mathbf{O}}^{\prime }$ respectively.

$$\begin{array}{c}\hfill |\overline{\mathrm{OT}}|=|\overline{\mathrm{OS}}|=\mathrm{r}\end{array}$$

(A15)

$$\begin{array}{c}\hfill \Rightarrow \angle OTS=\angle OST=\theta \end{array}$$

(A16)

$$\begin{array}{c}\hfill \Rightarrow \angle TOS=\pi -2\theta \end{array}$$

(A17)

Also,

$$\begin{array}{c}\hfill |\overline{{\mathrm{O}}^{\prime }\mathrm{T}}|=|\overline{{\mathrm{O}}^{\prime }{\mathrm{S}}^{\prime }}|=(1+\delta )\mathrm{r}\end{array}$$

(A18)

$$\begin{array}{c}\hfill \Rightarrow \angle O^{\prime }T{S}^{\prime }=\angle O^{\prime }{S}^{\prime }T=\theta \end{array}$$

(A19)

$$\begin{array}{c}\hfill \Rightarrow \angle T{O}^{\prime }{S}^{\prime }=\pi -2\theta \end{array}$$

(A20)

Hence

$$\begin{array}{c}\hfill \angle TOS=\angle T{O}^{\prime }{S}^{\prime }=\pi -2\theta \end{array}$$

Figure A5. ISP on a sphere displaced above the plane and a sphere centred at origin.

Figure A5. ISP on a sphere displaced above the plane and a sphere centred at origin.

Since both $\mathbf{S}$ and ${\mathbf{S}}^{\prime }$ lie on the same plane, which is a vertical cross-section of the sphere (plane perpendicular to the data plane and passing through the centre of both stereographic spheres), the azimuthal angle of both points is equal ($\varphi ={tan}^{-1}\left(\frac{p_y}{p_x}\right)$).

Hence one can see that the azimuthal and the polar angle generated by ISP on a sphere of radius $\mathrm{r}$ displaced above the 2D plane containing the points by $(1+\delta )\mathrm{r}$ is the same as the azimuthal and the polar angle generated by ISP on a sphere of radius $(1+\delta )\mathrm{r}$ centred at origin. This reduces the effective number of parameters that can be chosen for the embedding.