A Quantum Annealing-Based Approach to Extreme Clustering

TL;DR

This paper introduces a distributed quantum annealing method for extreme clustering, enabling faster and near-optimal clustering of large datasets with many clusters, outperforming traditional algorithms in speed.

Contribution

The paper presents a novel distributed quantum annealing approach for extreme clustering, demonstrating its efficiency and optimality under certain conditions.

Findings

Achieves near-optimal clustering assignments

Operates significantly faster than classical algorithms

Effective for large-scale datasets with many clusters

Abstract

Clustering, or grouping, dataset elements based on similarity can be used not only to classify a dataset into a few categories, but also to approximate it by a relatively large number of representative elements. In the latter scenario, referred to as extreme clustering, datasets are enormous and the number of representative clusters is large. We have devised a distributed method that can efficiently solve extreme clustering problems using quantum annealing. We prove that this method yields optimal clustering assignments under a separability assumption, and show that the generated clustering assignments are of comparable quality to those of assignments generated by common clustering algorithms, yet can be obtained a full order of magnitude faster.

Tables

Table 1. Table 1: Dataset statistics (left) and runtime comparison of extreme clustering algorithms in seconds (right). PERCH-C (“collapsed-mode”) was run, as it outperforms standard PERCH. The parameter L 𝐿 L sets the maximum number of leaves (see  [ 18 ] for an explanation). BiDViT selected the values ε 0 = 30 subscript 𝜀 0 30 \varepsilon_{0}=30 and ε 0 = 0.5 subscript 𝜀 0 0.5 \varepsilon_{0}=0.5 , such that a percentage of the nodes collapsed in the initial iteration, for the Covertype and the MNIST datasets, respectively. The mean and standard deviation were computed over five runs.

Name	Description	Cardinality	Dimension
MNIST	handwritten images	$60$ K	$784$
MNIST-2D	$t$-SNE of the above	$60$ K	$2$
Covertype	forest data	$581$ K	$54$
grid-2D	synthetically generated	$100$ K	$2$
grid-3D	synthetically generated	$1000$ K	$3$

Equations

$$\underset{S\subseteq X}{\mathrm{maximize}}\omega (S)\text{subject to}S\text{ is ε-separated}.$$

$$C_x=\{y\in X:d(x,y)\le d(x^{\prime },y)\text{ for all }{x}^{\prime }\in U\}.$$

$$\underset{\begin{array}{c}x,y\in C_i\\ 1\le i\le k\end{array}}{\max }d(x,y)<\underset{\begin{array}{c}x\in C_i,y\in C_j\\ 1\le i\ne j\le k\end{array}}{\min }d(x,y),$$

$$\epsilon \in \left(\underset{\begin{array}{c}x,y\in C_i\\ 1\le i\le k\end{array}}{\max }d(x,y),\underset{\begin{array}{c}x\in C_i,y\in C_j\\ 1\le i\ne j\le k\end{array}}{\min }d(x,y)\right],$$

$$\underset{s\in \{0,1{\}}^n}{\mathrm{maximize}}\sum _{i=1}^n{s}_i{w}_i\text{subject to}\sum _{i=1}^n\sum _{j>i}{s}_i{N}_{ij}^{(\epsilon )}{s}_j=0.$$

$$\underset{s\in \{0,1{\}}^n}{\mathrm{maximize}}\sum _{i=1}^n{s}_i{w}_i\text{subject to}{s}_i+s_j\le 1,\text{for }i,j\text{ such that }{\overline{N}}_{ij}^{(\epsilon )}=1.$$

$$\omega (S)\le {\left(\overline{{\mathrm{deg}}_w}(G)+1\right)}^{-1}\omega (S^{\ast }),$$

$$\underset{s\in \{0,1{\}}^n}{\mathrm{maximize}}\sum _{i=1}^n{s}_i{w}_i-\lambda \sum _{i=1}^n\sum _{j>i}{s}_i{N}_{ij}^{(\epsilon )}{s}_j.$$

$$\underset{s\in \{0,1{\}}^n}{\mathrm{minimize}}{s}^{T}Qs,$$

$$v^{T}Qv$$

$$=s^{T}Qs-\sum _{i\ne k}{s}_i{\lambda }_{\sigma (i,k)}{N}_{ik}^{(\epsilon )}+w_k,$$

$$v^{T}Qv\le s^{T}Qs-{\lambda }_{k\ell }+w_k.$$

$$\underset{s\in \{0,1{\}}^n}{\mathrm{minimize}},-s^{T}w\text{subject to}{s}^T{\overline{N}}^{(\epsilon )}s=0.$$

$$s^T{\overline{N}}^{(\epsilon )}s=0\iff s^T\left(\Lambda \circ {\overline{N}}^{(\epsilon )}\right)s=0,$$

$$p_1(s_1)\le p_1(s_2)=p_2(s_2)\le p_2(s_1)=p_1(s_1).$$

$$w_S(x)=\omega (C_x)=\sum _{y\in C_x}w(y).$$

$$\underset{S\subseteq V}{\mathrm{maximize}}\omega (S)\text{subject to}S\text{ is independent},$$

$${\mathcal{S}}_{\text{CH}}(C_1,\dots ,C_{n_{\text{c}}})=\left(\frac{n-1}{n_{\text{c}}-1}\right)\frac{{\sum }_{k=1}^{n_{\text{c}}}|C_k|\Vert c_k-c{\Vert }_2^2}{{\sum }_{k=1}^{n_{\text{c}}}{\sum }_{x\in C_k}\Vert x-c_k{\Vert }_2^2},$$

$${\mathcal{S}}_{\text{DB}}(C_1,\dots ,C_{n_{\text{c}}})=\frac{1}{n_{\text{c}}}\sum _{k=1}^{n_{\text{c}}}\underset{j\ne k}{\max }\frac{S_k+S_j}{\Vert c_k-c_j{\Vert }_2},$$

$$\underset{s\in \{0,1{\}}^n}{\mathrm{minimize}}\sum _{i=1}^n{s}_i{c}_i\text{subject to}\sum _{j=1}^n{N}_{ij}^{(\epsilon )}{s}_j\ge 1,i=1,\dots n.$$

$$0\le {\xi }_i\le \left(\sum _{j=1}^n{N}_{ij}^{(\epsilon )}\right)-1,\text{ for }i=1,\dots ,n.$$

$$\underset{\begin{array}{c}s\in \{0,1{\}}^n\\ 0\le \xi \le (N^{(\epsilon )}1)-1\end{array}}{\mathrm{minimize}}\sum _{i=1}^n{s}_i{c}_i+\lambda \sum _{i=1}^n{\left[\left(\sum _{j=1}^n{N}_{ij}^{(\epsilon )}{s}_j\right)-1-{\xi }_i\right]}^2.$$

$${\xi }_i=\sum _{k=0}^{\lfloor m_i\rfloor }{b}_k^{(i)}{\gamma }_k^{(i)},\text{for }i=1,\dots ,n,$$

$$0\le \sum _{k=0}^{\lfloor m_i\rfloor }{b}_k^{(i)}{\gamma }_k^{(i)}\le \left(\sum _{j=1}^n{N}_{ij}^{(\epsilon )}\right)-1,$$

$$\underset{\begin{array}{c}{b}^{(i)}\in \{0,1{\}}^{\lfloor m_i\rfloor +1}\\ s\in \{0,1{\}}^n\end{array}}{\mathrm{minimize}}\sum _{i=1}^n{s}_i{c}_i+\lambda {\left[\left(\sum _{j=1}^n{N}_{ij}^{(\epsilon )}{s}_j\right)-1-\sum _{k=0}^{\lfloor m_i\rfloor }{b}_k^{(i)}{\gamma }_k^{(i)}\right]}^2.$$

Full text

††footnotetext: ∗ T. Jaschek and J. S. Oberoi have contributed equally to this research. Address correspondence to: tim.jaschek@1qbit.com

1 Introduction

Traditionally, clustering approaches have been developed and customized for tasks where the resultant number of clusters $k$ is not particularly high. In such cases, algorithms such as $k$-means++ [2], BIRCH [37], DBSCAN [9], and spectral clustering produce high-quality solutions in a reasonably short amount of time. This is because these traditional algorithms scale well with respect to the dataset cardinality $n$. However, in most cases, the computational complexity of these algorithms, in terms of the number of clusters, is either exponential or higher-order polynomial. Another common issue is that some of the algorithms require vast amounts of memory.

The demand for clustering algorithms capable of solving problems with larger values of $k$ is continually increasing. Present-day examples involve deciphering the content of billions of web pages by grouping them into millions of labelled categories [31, 34], identifying similarities among billions of images using nearest-neighbour detection [35, 23, 36]. This domain of clustering, where $n$ and $k$ are both substantially large, is referred to as extreme clustering [18]. Although there is great value in perfecting this type of clustering, very little effort towards this end has been made by the machine learning community. Our algorithm is, in fact, such an effort. Its output is a clustering tree, which can be used to generated multiple clustering assignments (or “levels”) with varying degrees of accuracy (i.e., coarseness or fineness) of the approximation. Generating such a tree is not uncommon for clustering algorithms. Consider, for example, hierarchical clustering algorithms which generate binary clustering trees. Clustering trees are useful tools for dealing with real-world data visualization problems. Our algorithm, the Big Data Visualization Tool, or BiDViT, provides this functionality.

BiDViT employs a novel approach to clustering problems, which is based on the maximum weighted independent set (MWIS) problem in a graph induced by the original dataset and a parameter we call the radius of interest or neighbourhood parameter, which determines a relation of proximity. The use of such a parameter has been successfully employed in density-based spatial clustering of applications with noise (DBSCAN) [9]. The MWIS problem can be transformed into a quadratic unconstrained binary optimization (QUBO) problem, the formulation accepted by a quantum annealer. An alternative way to address the underlying problem is to use a heuristic algorithm to approximate solutions to the MWIS problem. Quantum annealing and simulated annealing have been applied in centroid-based clustering [20, 30] and in density-based clustering [21]. However, the approaches studied are not capable of addressing problems in the extreme clustering domain.

We prove that, under a separability assumption on the ground truth clustering assignment of the original dataset, our method identifies the ground truth labels when parameters are selected that are within the bounds determined by that assumption. We provide runtime and solution quality values for both versions of our algorithm, with respect to internal evaluation schemes such as the Calinski–Harabasz and the Davies–Bouldin scores. Our results suggest that BiDViT yields clustering assignments of a quality comparable to that of assignments generated by common clustering algorithms, yet does so a full order of magnitude faster.

2 The Coarsening Method

Our algorithm is based on a combinatorial clustering method we call coarsening. The key idea behind coarsening is to approximate a set $X\subset\mathbb{R}^{d}$ by a subset $S\subseteq X$ such that, for any point $x\in X$, there exists a $y\in S$ such that $\|x-y\|_{2}<\varepsilon$, for some parameter $\varepsilon>0$. In this case, we say that $S$ is $\varepsilon$-dense in $X$ and call $\varepsilon$ the radius of interest. This concept is not restricted to subsets of Euclidean spaces and can be generalized to an arbitrary metric space $(M,d)$. For example, our coarsening method can be used for clustering assignments on finite subsets of Riemannian manifolds with respect to their geodesic distance, for instance, in clustering GPS data on the surface of the Earth when analyzing population density. In what follows, we assume that $X=\{x^{(1)},\ldots,x^{(n)}\}$ is a dataset consisting of $n$ $d$-dimensional data points, equipped with a metric $d:X\times X\rightarrow[0,\infty)$. Finding an arbitrary $\varepsilon$-dense subset of $X$ does not necessarily yield a helpful approximation. For example, $X$ itself is always $\varepsilon$-dense in $X$. However, enforcing the additional constraint that any two points in the subset $S$ must be separated by a distance of at least $\varepsilon$ yields more-interesting approximations, often leading to a reduction in the number of data points (one of our primary objectives). We call such a set $\varepsilon$-separated. Fig. 1 shows a point cloud and an $\varepsilon$-dense, $\varepsilon$-separated subset. The theorem that follows shows that a maximal $\varepsilon$-separated set $S$ of $X$ is necessarily $\varepsilon$-dense in $X$. Let $B(x,r)$ denote the open metric ball with respect to $d$, with centre $x$ and radius $r$.

Theorem 2.1**.**

Let $S$ be a maximal $\varepsilon$-separated subset of $X$ in the sense of set inclusion. Then the following properties must be satisfied.

• i)

We have the inclusion $X\subseteq\bigcup_{x\in S}B(x,\varepsilon)$.

• ii)

For every $y\in S$, it holds that $X\nsubseteq\bigcup_{x\in S\setminus\{y\}}B(x,\varepsilon)$.

• iii)

The sets $B(x,\varepsilon/2)$ for $x\in S$ are pairwise disjoint.

In particular, $S$ is a minimal $\varepsilon$-dense subset of $X$.

Proof 2.2**.**

Note that $i)$ is equivalent to $S$ being $\varepsilon$-dense in $X$ and that, in combination with $ii)$, is equivalent to $S$ being a minimal with respect to this property. To prove $i)$, let $S$ be a maximal $\varepsilon$-separated subset of $X$ and assume, in contradiction, that $S$ is not $\varepsilon$-dense in $X$. Then we could find $x\in X$ such that $d(x,y)\geq\varepsilon$, for every $y\in S$. Hence, $S\cup\{x\}$ would be $\varepsilon$-separated, which is in contradiction to the maximality of $S$. To prove $ii$), we fix a point $x\in S$. Since $S$ is $\varepsilon$-separated, $d(x,y)\geq\varepsilon$ for any $y\in S$ and, thus, $S\setminus\{x\}$ is not $\varepsilon$-dense in $X$. Property $iii$) follows from the triangle inequality.∎

Note that a maximal $\varepsilon$-separated subset does not refer to an $\varepsilon$-separated subset with fewer than or equally as many elements as all other $\varepsilon$-separated subsets but, rather, to an $\varepsilon$-separated subset that is no longer $\varepsilon$-separated when a single data point is added. Contrary to  Thm. 2.1, a minimal $\varepsilon$-dense subset does not need to be $\varepsilon$-separated. Consider the set $X=\{1,2,3,4\}\subset\mathbb{R}$, and let $d$ be the Euclidean distance on $\mathbb{R}$. Then, $S=\{2,3\}$ is $\nicefrac{{3}}{{2}}$-dense in $X$ but not $\nicefrac{{3}}{{2}}$-separated. Also note that an $\varepsilon$-separated subset is not necessarily an $\varepsilon$-coreset, which is a weighted subset whose weighted $k$-means cost approximates the $k$-means cost of the original set with up to an accuracy of $\varepsilon$ [12, 4].

In the following, we assume that $X$ is equipped with a weight function $w:X\rightarrow\mathbb{R}_{+}$. We call $w_{i}=w(x^{(i)})$ the weight of $x^{(i)}$ and gather all weights in a weight vector $w\in\mathbb{R}^{n}_{+}$. It will be clear from the context whether we refer to a weight function or a weight vector. The weight of a set $S\subseteq X$ is given by $\omega(S)=\sum_{x\in S}w(x)$. We have already argued that maximal $\varepsilon$-separated subsets yield reasonable approximations. However, such subsets are not unique. We are thus interested in finding an optimal one, that is, one that captures most of the weight of the original dataset. In other words, we are interested in solving the optimization problem

[TABLE]

If we impose unit weights, the solution set to this optimization problem will consist of the maximal $\varepsilon$-separated subsets of $X$ with a maximum number of elements among all such subsets. The term “maximal” refers to set inclusion and the “maximum” refers to set cardinality. Since $w(x)>0$ for all $x\in X$, a solution $S^{*}$ to P0 will always be a maximal $\varepsilon$-separated subset and, therefore, by Thm. 2.1, $\varepsilon$-dense. In Sec. 3.6, we show that this problem is equivalent to solving an MWIS problem for a weighted graph $G^{\varepsilon}(X,E^{\varepsilon},w)$, depending solely on the dataset $X$, the Euclidean metric $d$, and the radius of interest $\varepsilon$. Thus, the computational task of finding a maximal $\varepsilon$-separated subset of maximum weight is NP-hard [26, 17].

Every subset $U\subset X$ gives rise to a clustering assignment $\mathcal{C}=\{C_{x}\}_{x\in U}$. This assignment is given by

[TABLE]

Data points that are equidistant to multiple representative points are assigned to only one of them, uniformly at random. Typically, larger values of $\varepsilon$ result in smaller cardinalities of $\mathcal{C}$. The following corollary summarizes properties of $\mathcal{C}$ when $U$ is $\varepsilon$-separated, and can be readily verified.

Corollary 2.3**.**

Let $\mathcal{C}$ be the clustering assignment generated from a maximal $\varepsilon$-separated set $S\subset X$. Then, the following properties are satisfied:

• i)

The clusters in $\mathcal{C}$ are non-empty and pairwise disjoint.

• ii)

The cluster diameter is uniformly bounded by $2\varepsilon$, i.e., $\sup_{x\in S}\textup{diam}(C_{x})\leq 2\varepsilon$.

• iii)

For all $x\in S$, it holds that $\max_{y\in C_{x}}d(x,y)<\varepsilon$.

Notice that these properties are not satisfied by every clustering assignment, for example, the ones generated by $k$-means clustering. They are desirable in specific applications, such as image quantization, where a tight bound on the absolute approximation error is desired. However, they are undesirable if the ground truth clusters have diameters larger than $2\varepsilon$. More details on the clustering assignment are provided in Sec. 3.

One could argue that prior to identifying a maximum weighted independent set and using it to generate a clustering assignment, a dataset should be normalized. However, normalization is a transformation that would result in chunks not being defined by metric balls, but rather by ellipsoids. In particular, such a transformation would change the metric $d$. We assume that the metric $d$ already is the best indicator of proximity. In general, one can apply any homeomorphism $f$ to a dataset $X$, apply our clustering algorithm to the set $f(X)$, and obtain a clustering assignment by applying $f^{-1}$ to the individual clusters.

A common assumption in the clustering literature is separability—not to be mistaken with $\varepsilon$-separability—of the dataset with respect to a clustering $\mathcal{C}$. The dataset $X$ is called separable with respect to a clustering $\mathcal{C}=\{C_{1},\ldots C_{k}\}$ if

[TABLE]

that is, if the maximum intra-cluster distances are strictly smaller than the minimum inter-cluster distances. The following theorem shows that, if $\varepsilon$ is chosen correctly, our coarsening method yields the clustering assignment $\mathcal{C}$.

Theorem 2.4**.**

Let $X$ be separable with respect to a clustering $\mathcal{C}=\{C_{1},\ldots C_{k}\}$. Then, for any

[TABLE]

our coarsening methods yields the correct clustering assignment.

Proof 2.5**.**

To simplify our notation, we denote the lower and upper bounds of the interval in (3) by $l$ and $r$, respectively. By the separability assumption, this interval is non-empty. One can see that, for any admissible choice of $\varepsilon$, any two points from different clusters are $\varepsilon$-separated. Indeed, for $x\in C$ and $y\in C^{\prime}$, it holds that $d(x,y)\geq r\geq\varepsilon$. Furthermore, if a point $x$ in a cluster $C$ is selected, then no other point $y$ in the same cluster can be selected, as $d(x,y)\leq l<\varepsilon$. Therefore, every solution $S\subseteq X$ to P0 is a union of exactly one point from each cluster. Using the separability of $X$ with respect to $\mathcal{C}$, we can see that the clustering assignment induced by 1 is coincident with $\mathcal{C}$. ∎

In practice, the separability assumption is rarely satisfied, and it is challenging to select $\varepsilon$ as above (as this assumes some knowledge about the clustering assignment). However, Thm. 2.4 shows that our coarsening method is of research value, and can potentially yield optimal clustering assignments.

We have developed two methods, which we refer to as the heuristic method and the quantum method, to address the NP-hard task of solving P0. The heuristic method loosens the condition of having a maximum weight; it can be seen as a greedy approach to P0. In contrast, the quantum method explores all different maximal $\varepsilon$-separated subsets simultaneously, yielding one that has maximum weight. The quantum method is based on the formulation of a QUBO problem, which can be solved efficiently using a quantum annealer like the D-Wave 2000Q [7] or a digital annealer such as the one developed by Fujitsu [10].

3 The Algorithm

Let $X=\{x^{(1)},\ldots,x^{(n)}\}\subset\mathbb{R}^{d}$ denote a dataset of $n$ $d$-dimensional data points. Note that, mathematically speaking, a dataset is not a set but rather a multiset, that is, repetitions are allowed. Our algorithm consists of two parts: data partitioning and data coarsening, the latter of which can be further subdivided into chunk coarsening and chunk collapsing.

3.1 Data Partitioning

In general, the computational complexity of distance-based clustering methods is proportional to the square of the dataset cardinality, as all pairwise distances must be computed. This bottleneck can be overcome by dividing the dataset and employing distributed approaches [4, 28, 12], yielding a result different from the one we would obtain when applying clustering methods on the entire dataset. However, its slight imprecision results in a significant computational speed-up.

A partition $\mathcal{P}$ of $X$ is a collection of non-empty disjoint sets $P_{1},\ldots,P_{k}\subset X$ such that $X=\bigcup_{P\in\mathcal{P}}P$. Elements of partitions are typically referred to as blocks, parts, or cells; however, we refer to them as chunks. The partitioning is intended to be homogeneous: every extracted chunk has an equal number of data points (there might be minor differences when the cardinality of the chunk to be divided is odd). The upper bound on the number of points desired in a chunk is referred to as the maximum chunk cardinality $\kappa$. To determine $\kappa$, one should take into account the number of available processors, their data handling capacity, or, in the case of a quantum annealer, the number of fully connected qubits.

To break the data into chunks, we employ a modified version of the well-known “median cut” algorithm, which is frequently used in colour quantization. First, we select an axis of maximum variance. We then bisect the dataset along the selected axis, say $\ell$, at the median $m$ of $\{x^{(1)}_{\ell},\ldots,x^{(n)}_{\ell}\}$ in such a way as to obtain two data chunks $P_{1}$ and $P_{2}$ whose cardinalities differ by at most one (in the case where $n$ is odd) and which satisfy $P_{1}\subseteq\{x\in X:x_{\ell}\leq m\}$ and $P_{2}\subseteq\{x\in X:x_{\ell}\geq m\}$. We cannot simply assign $P_{1}=\{x\in X:x_{\ell}\leq m\}$ and $P_{2}=X\setminus P_{1}$, as these sets might differ drastically in cardinality. For example, when $x_{\ell}^{(1)}=\ldots=x_{\ell}^{(n)}$, this assignment would imply that $P_{1}=X$ and $P_{2}=\emptyset$.

By using $P_{1}$ and $P_{2}$ in the role of $X$, this process can be repeated iteratively, until the number of data points in the chunk to be divided is less than or equal to the maximum chunk cardinality $\kappa$, yielding a binary tree of data chunks. After $s$ iterations, this leaves us with $2^{s}$ chunks $P^{(s)}_{k}$ such that $X=\bigcup_{1\leq k\leq 2^{s}}P^{(s)}_{k},$ where the union is disjoint. Fig. 1 provides a visualization.

3.2 Chunk Coarsening

The goal of a data coarsening step is, for each chunk, to find representative data points such that their union can replace the original point cloud, while maintaining the original data distribution as accurately as possible.

Let $P=\{x^{(1)},\ldots,x^{(n)}\}$ be a chunk and $\varepsilon>0$ be the radius of interest. In what follows, we assume that all the data points are pairwise different. Practically, this can be achieved by removing duplicates and cumulatively incrementing the weight of the representative point we wish to keep by the weight of the discarded duplicates. The radius of interest $\varepsilon$ induces a weighted graph $G^{\varepsilon}=(P,E^{\varepsilon},w_{P})$, where $P$ is the vertex set, the edge set $E^{\varepsilon}$ is given by the relation $\sim_{\varepsilon}$ defined by $x\sim_{\varepsilon}y$ if and only if $d(x,y)<\varepsilon$ for all $x,y\in P$, and the weight function $w_{P}:P\rightarrow\mathbb{R}_{+}$ is the restriction of $w$ to $P$. For each data point $x^{(i)}$, we denote its weight $w_{P}(x^{(i)})$ by $w_{i}$.

For each data point $x^{(i)}$, we introduce a binary decision variable $s_{i}$ that encodes whether $x^{(i)}$ is used in a possible set $S^{*}$. Furthermore, we define the neighbourhood matrix $N^{(\varepsilon)}$ (or similarity matrix) of the graph $G^{\varepsilon}=(P,E^{\varepsilon},w_{P})$ by $N_{ij}^{(\varepsilon)}=1$ if $x^{(i)}\sim_{\varepsilon}x^{(j)}$, and $N_{ij}^{(\varepsilon)}=0$ otherwise. Problem P0 can then be posed as a quadratically constrained quadratic program (QCQP) given by

[TABLE]

Here, the inner summation of the constraint does not need to run over all indices, due to the symmetry of $N^{(\varepsilon)}$. The matrix form of P1 is given by maximizing $s^{T}w$ subject to the constraint $s^{T}\overline{N}^{(\varepsilon)}s=0$, where $\overline{N}^{(\varepsilon)}$ is the upper triangular matrix of $N^{(\varepsilon)}$ having all zeroes along the diagonal. As explained in Sec. 3.6, P1 is equivalent to the NP-hard MWIS problem for $G^{\varepsilon}=(P,E^{\varepsilon},w_{P})$, and thus is computationally intractable for large problem sizes. Note that P1 can be written as the 0–1 integer linear program (ILP)

[TABLE]

We present two methods we have devised to address P1.

3.2.1 The Heuristic Method

We wish to emphasize that the heuristic method does not provide us with a solution to P1. Rather, the aim of this method is to obtain an $\varepsilon$-separated subset $S$ with a high—but not necessarily the maximum—weight $\omega(S)$. The seeking of approximate solutions to the MWIS problem is a well-studied subject [3, 14, 16]. Typically, researchers employ greedy algorithms, LP-based algorithms (using the relaxation of P2), or semi-definite programming (SDP) algorithms; see [16] for an analysis.

We employ a classic greedy algorithm due to its simplicity and low computational complexity. In each step, we add the data point that locally is the best choice in the sense that the ratio of the weight of its neighbourhood to its own weight is as small as possible. Prior to the execution of the step, we remove the point and its neighbours from the set of candidates. Pseudocode of the greedy algorithm is provided in Algorithm 1. Before we state a theoretical result on the approximation ratio of this algorithm, we define the weighted degree $\deg_{w}(v)$ of a vertex $v$ in a weighted graph $G=(V,E,w)$ and the weighted average degree of $G$ as $\deg_{w}(v)=\omega(N_{v})/w(v)$ and $\overline{\deg_{w}}(G)=\sum_{v\in V}w(v)\deg_{w}(v)/\omega(V),$ respectively, where $N_{v}=\{u\in V:u\sim v\}$ is the neighbourhood of vertex $v$ [16].

Theorem 3.1**.**

Algorithm 1* has an approximation ratio of $\overline{\deg_{w}}(G)+1$, i.e.,*

[TABLE]

for any solution $S^{*}$ to P0 and any output $S$ of the algorithm. Moreover, the bound in (4) is tight.

Proof 3.2**.**

A proof is given in [16, Thm. 6]. ∎

3.2.2 The Quantum Method

In contrast to the heuristic method, the QUBO approach provides an actual (i.e., non-approximate) solution to P1. We reformulate the problem by transforming the QCQP into a QUBO problem.

Using the Lagrangian penalty method, we incorporate the constraint into the objective function by adding a penalty term. For a sufficiently large penalty multiplier $\lambda>0$, the solution set of P1 is equivalent to that of

[TABLE]

One can show that, for $\lambda>\max_{i=1,\ldots n},w_{i}$ every solution to P3 satisfies the separation constraint [1, Thm. 1]. Instead, we use individual penalty terms $\lambda_{ij}$, as this may lead to a QUBO problem with much smaller coefficients, which results in improved performance when solving the problem using a quantum annealer. Expressing P3 as a minimization, instead of a maximization, problem and using matrix notation yields the problem

[TABLE]

where $Q_{ij}=-w_{i}$ if $i=j$, $Q_{ij}=\lambda_{ij}$ if $N_{ij}^{(\varepsilon)}=1$ and $i<j$, and $Q_{ij}=0$ otherwise. Solutions to P4 can be approximated using heuristics such as simulated annealing [32], path relinking [25], tabu search [25], and parallel tempering [38]. Before solving P4, it is advisable to reduce its size and difficulty by making use of logical implications among the coefficients [11]. This involves fixing every variable that corresponds to a node that has no neighbours to one, as it necessarily is included in an $\varepsilon$-dense subset.

The following theorem show that P1 is equivalent to P4 for a suitable choice of $\lambda_{ij}$, for $1\leq i<j\leq n$.

Theorem 3.3**.**

Let $\lambda_{ij}>\max\{w_{i},w_{j}\}$ for all $1\leq i<j\leq n$. Then, for any solution $s\in\{0,1\}^{n}$ to P4, the corresponding set $S\subseteq X$ is $\varepsilon$-separated. In particular, the solution sets of P1 and P4 coincide.

Proof 3.4**.**

We generalize the proof of [1, Thm. 1] and show that every solution $s$ to P4 satisfies the separation constraint $\sum_{i=1}^{n}\sum_{j>i}s_{i}N_{ij}^{(\varepsilon)}s_{j}=0$. Assuming, in contradiction, that the opposite were to be the case, we could find a solution $s$ and indices $k$ and $\ell$ such that $1\leq k<\ell\leq n$ and $s_{k}=s_{\ell}=N_{k\ell}^{(\varepsilon)}=1$. Let $e_{k}$ denote the $k$-th standard unit vector, and let $v=s-e_{k}$. Then,

[TABLE]

where $\sigma:\mathbb{N}^{2}\rightarrow\mathbb{N}^{2}$, defined by $\sigma(i,k)=(\min(i,k),\max(i,k))$, orders the index accordingly. This technicality is necessary, as we defined $\lambda_{ij}$ only for $1\leq i<j\leq n$. As $N_{k\ell}^{(\varepsilon)}=s_{\ell}=1$, we have $\sum_{i\neq k}s_{i}\lambda_{\sigma(i,k)}N_{ik}^{(\varepsilon)}\geq\lambda_{\sigma(\ell,k)}$, and thus

[TABLE]

Therefore, as $\lambda_{k\ell}>\max\{{w_{k},w_{\ell}\}}\geq w_{k}$, it holds that $v^{T}Qv<s^{T}Qs$, which is absurd, as, by assumption, $s$ is a solution to P4.

We now show that the solution sets of P1 and P4 coincide. Note that P1 is equivalent to the optimization problem

[TABLE]

Let $s_{1}$ and $s_{2}$ be solutions to P5 and P4, respectively. We denote the objective functions by $p_{1}(s)=-s^{T}w$ and $p_{2}(s)=-s^{T}w+s^{T}\left(\Lambda\circ\overline{N}^{(\varepsilon)}\right)s$, where $\Lambda$ is the matrix defined by $\Lambda_{ij}=\lambda_{ij}$ for $1\leq i<j\leq n$, and zero otherwise, and the term $\Lambda\circ\overline{N}^{(\varepsilon)}\in\mathbb{R}^{n\times n}$ denotes the Hadamard product of the matrices $\Lambda$ and $\overline{N}^{(\varepsilon)}$, given by element-wise multiplication. Then, as $\lambda_{ij}>\max\{w_{i},w_{j}\}$ for $1\leq i<j\leq n$, by the observation above, both $s_{1}$ and $s_{2}$ satisfy the separation constraint. Since $s$ and $\overline{N}^{\varepsilon}$ are coordinate-wise non-negative and $\lambda_{ij}>\min_{k=1,\ldots,n}w_{k}>0$ for $1\leq i<j\leq n$, it holds that

[TABLE]

thus, if $s$ satisfies the separation constraint, then $p_{2}(s)=p_{1}(s)$. Using this observation, and that $s_{1}$ and $s_{2}$ minimize $p_{1}$ and $p_{2}$, respectively, we have

[TABLE]

Hence, the inequalities in 9 must actually be equalities; thus, the solution sets of the optimization problems coincide. ∎

Problem P4 can be transformed to an Ising spin model by mapping $s$ to $2s-1$. This form is desirable because the ground state of the Hamiltonian of an Ising spin model can be determined efficiently with a quantum annealer.

3.3 Chunk Collapsing

Having identified a maximal $\varepsilon$-separated subset $S\subseteq P$, we collapse the vertices $P\setminus S$ into $S$, meaning we update the weight of each $x\in S$ according to the weights of all $y\in P\setminus S$ that satisfy $x\sim_{\varepsilon}y$. We aim to assign each $y\in P\setminus S$ to a unique $x\in S$ by generating a so-called Voronoi decomposition (depicted in Fig. 1) of each chunk $P$, which is a partition, where each point $x\in P$ is assigned to the closest point within a subset $S$. More precisely, we define the sets $C_{x}$ as in 1, for each $x\in S$. By construction, $C_{x}$ contains all vertices that will be collapsed into $x$, in particular, $x$ itself. We then assign the coarsened chunk $S$ a new weight function $w_{S}$ defined by

[TABLE]

In practice, to prevent having very large values for the individual weights, one might wish to add a linear or logarithmic scaling to this weight assignment. In our experiments, we did not add such a scaling.

3.4 Iterative Implementation of BiDViT

BiDViT repeats the procedure of data partitioning, chunk coarsening, and chunk collapsing with an increasing radius of interest, until the entire dataset collapses to a single data point. We call these iterations BiDViT levels. The increase of $\varepsilon$ between BiDViT levels is realized by multiplying $\varepsilon$ by a constant factor, denoted by $\alpha$ and specified by the user. In our implementation we have introduced a node class that has three attributes: coordinates, weight, and parents. We initialize BiDViT by creating a node_list containing the nodes corresponding to the weighted dataset (if no weights are provided then the weights are assumed to be the multiplicity of the data points). After each iteration, we remove the nodes that collapsed into representative nodes from the node_list and keep only the remaining representative nodes. However, we append the removed nodes to the parents of the representative node. The final node_list is a data tree, that is, it consists of only one node, and we can move upwards in the hierarchy by accessing its parents (and their parents and so on); see Fig. 2. Two leaves of the data tree share a label with respect to a specific BiDViT level, say $m$, if they have collapsed into centroids which, possibly after multiple iterations, have collapsed into the same centroid at the $m$-th level of the tree. For the sake of reproducibility, we provide pseudocode (see Algorithm 2).

It is worth noting that, at each level, instead of proceeding with the identified representative data points, one can use the cluster centroids, allowing more-accurate data coarsening and label assignment.

3.5 Complexity Analysis

Our analysis shows that every interation of the heuristic version of BiDViT has a computational complexity of $\text{O}(dn\log(n/\kappa)+dn\kappa)$. Note that $\kappa\ll n$.

The order of complexity of the partitioning procedure is $\text{O}(dn\log(n/\kappa))$. To see this, note that there are at most $\lceil\log_{2}(n/\kappa)\rceil$ partitioning stages and in the $s$-th stage we split $2^{s-1}$ chunks $P_{i}$, where $i=1,\ldots,2^{s-1}$. Let $n_{i}$ denote the number of data points in chunk $P_{i}$. Finding the dimension of maximum variance has a complexity of $\text{O}(dn_{i})$ and determining the median of this dimension can be achieved in $\text{O}(n_{i})$ via the “median of medians” algorithm. Having computed the median, one can construct two chunks of equal size in linear time. Since $\sum_{1\leq i\leq 2^{s-1}}n_{i}=n$, a partitioning step is $\text{O}(dn\log(n/\kappa))$. Any division of a chunk is independent of the other chunks at a given stage; thus, this procedure can benefit from distributed computing.

The order of complexity for the collapsing process is $\text{O}(dn\kappa)$, as computing the neighbourhood matrix of a chunk is $\text{O}(d\kappa^{2})$ and the heuristic selection procedure is $\text{O}(\kappa^{2})$. The number of chunks is bounded from above by $\lceil n/\kappa\rceil$. This yields a complexity of $\text{O}((n/\kappa)(d\kappa^{2}+\kappa^{2}))=\text{O}(dn\kappa)$. As data coarsening in each chunk is independent, with $\lceil n/\kappa\rceil$ parallel processors available the complexity reduces to $\text{O}(d\kappa^{2})$.

3.6 Relation to the MWIS Problem

The process of identifying a maximal $\varepsilon$-separated set of maximum weight is equivalent to solving the MWIS problem for the weighted graph $G^{\varepsilon}=(P,E^{\varepsilon},w_{P})$. Let $G=(V,E,w)$ be a weighted graph. A set of vertices $S\subseteq V$ is called independent in $G$ if no two of its vertices are adjacent or, equivalently, if $S$ is a clique in the complement graph. This corresponds to the separation constraint mentioned earlier, where two vertices are adjacent whenever they are less than a distance of $\varepsilon$ apart. The MWIS problem can be expressed as

[TABLE]

and is NP-complete for a general weighted graph [17], yet, for specific graphs, there exist polynomial-time algorithms [29, 19]. Note that the QUBO formulation of the MWIS problem in [1, 13] is related to the one in P4.

If all weights are positive, a maximum weighted independent set is necessarily a maximal independent set. A maximal independent set is a dominating set, that is, a subset $S$ of $V$ such that every $v\in V\setminus S$ is adjacent to some $w\in S$. This corresponds to our observation that every maximal $\varepsilon$-separated subset is $\varepsilon$-dense.

4 Results

The datasets used to demonstrate the efficiency and robustness of our approach are the MNIST dataset of handwritten digits [22], a two-dimensional version of MNIST obtained by using $t$-SNE [27], two synthetic grid datasets, and a dataset called Covertype containing data on forests in Colorado [5]. The synthetic grid datasets are the unions of 100 samples (in the 2D case) and 1000 samples (in the 3D case) drawn from $\mathcal{N}(\mu_{ij},\sigma^{2})$ with means $\mu_{ij}=(10i+5,10j+5)$ and a variance of $\sigma_{2}=4$ for $0\leq i,j\leq 9$ in the 2D case and the natural extension in the 3D case. Dataset statistics are provided in Table 1. In addition to our technical experiments, explained in the following sections, a practical application of BiDViT for image qunatization is illustrated in Fig. 7. All experiments were performed using a 2.5 GHz Intel Core i7 processor and 16 GB of RAM.

4.1 Low-Range Clustering Domain

Although BiDViT has been specifically designed for extreme clustering, it yields accurate assignments for low values of $k$. Fig. 3 shows the clustering assignment of BiDViT on the 2D grid dataset and on MNIST. The results are obtained by manually selecting a BiDViT level. In the grid dataset, every cluster is identified correctly. In the MNIST dataset, all clusters are recognized, except one. However, as our algorithm is based on metric balls, and some datasets might not conform to such categorization, there are datasets for which it cannot accurately assign clusters. This is true for most clustering algorithms, as they are able to recognize only specific shapes.

4.2 Extreme Clustering Capability

To evaluate the performance of BiDViT on high-dimensional datasets in the extreme clustering range, we used the Calinski–Harabasz score [6] and the Davies–Bouldin score [8]. These clustering metrics are internal evaluation schemes, that is, their values depend solely on the clustered data, not requiring the ground truth label assignment for the dataset. Such schemes must be viewed as heuristic methods: their optimal values do not guarantee optimal clusters but provide a reasonable measure of clustering quality. Detailed analyses have been conducted on the advantages and shortcomings of internal clustering measures [24, 15]. In the extreme clustering scenario, where the objective is to obtain an accurate approximation of the entire dataset instead of categorizing its elements, no true labels are given and thus external evaluation schemes (ones based on the distance to a ground truth clustering assignment) do not qualify as success measures.

Let $C_{1},\ldots,C_{n_{\text{c}}}$ denote a total of $n_{\text{c}}$ detected clusters within a dataset $X$ with $n$ data points. The Calinski–Harabasz score $\mathcal{S}_{\text{CH}}$ of a clustering is defined as a weighted ratio of the inter-cluster squared deviations to the sum of the intra-cluster squared deviations. More precisely, $\mathcal{S}_{\text{CH}}$ is given by

[TABLE]

where $c_{k}$, for $k=1,\ldots,n_{\text{c}}$ are the cluster centroids, and $c$ is their mean. High values of $\mathcal{S}_{\text{CH}}$ are indicative of a high clustering quality. The Davies–Bouldin score $\mathcal{S}_{\text{DB}}$ is the average maximum value of the ratios of the pairwise sums of the intra-cluster deviation to the inter-cluster deviation. The score is defined as

[TABLE]

where $S_{i}=\sum_{x\in C_{i}}\|x-c_{i}\|/|C_{i}|.$ Low values of $\mathcal{S}_{\text{DB}}$ indicate accurate clustering.

Fig. 4 shows $\mathcal{S}_{\text{CH}}$ and $\mathcal{S}_{\text{DB}}$ of clustering assignments obtained with BiDViT and Mini Batch $k$-means clustering [33] for different values of $k$ on the Covertype dataset. Due to their high computational complexity with respect to $k$, many common clustering algorithms could not be applied. Remarkably, $\mathcal{S}_{\text{CH}}$ values are quite similar, indicating that the cluster assignments generated by BiDViT are of comparable quality even though the runtime of our algorithm is significantly shorter. For $\mathcal{S}_{\text{DB}}$, our algorithm outperforms the others for lower values of $k$, and is comparable for large values. One explanation for the slightly weaker performance of BiDViT with respect to $\mathcal{S}_{\text{CH}}$ is that BiDViT aims to minimize the non-squared distances, whereas $\mathcal{S}_{\text{CH}}$ rewards clustering methods that minimize squared distances. Similarly, this explains BiDViT’s advantage for $\mathcal{S}_{\text{DB}}$.

4.3 Runtime Comparison

In our experiments, we observed that, with respect to the total runtime, even the heuristic version of BiDViT restricted to a single core outperforms common clustering methods in the extreme clustering domain. Fig. 5 shows the runtime required by different clustering algorithms for the Covertype dataset. For the implementation of methods other than BiDViT, we used the publicly available sklearn.clustering module for Python. To generate the plots, we ran the entire BiDViT procedure, then applied classical algorithms for the same values of $k$. The results suggest that, in the extreme clustering domain, the runtime of BiDViT is an order of magnitude faster than that of the agglomerative methods against which it was compared, and multiple orders of magnitude faster than that of $k$-means and Mini Batch $k$-means clustering. The dataset cardinality was restricted to 20,000 points to obtain results for other methods, whereas BiDViT is capable of handling the entire dataset comprising 581,000 points.

We then compared the runtime of BiDViT to PERCH (“Purity Enhancing Rotations for Cluster Hierarchies”), a hierarchical algorithm for extreme clustering [18], to our knowledge the only other algorithm designed to solve extreme clustering problems. We restricted both algorithms to using a single core. Table 1 shows that BiDViT performs an order of magnitude faster than PERCH. However, they solve somewhat different problems: whereas BiDViT aims to gradually coarsen a dataset by finding $\varepsilon$-separated, $\varepsilon$-dense subsets, PERCH maximizes the dendrogram purity, a measure of the clustering tree’s consistency [18]. The clustering tree generated by PERCH is binary and thus enormous, allowing for much finer incremental distinctions between clustering assignments. In contrast, the tree generated by BiDViT is more compact, as multiple data points can collapse into the same representative point. When comparing dendrogram purities, we expect PERCH to outperform BiDViT; when comparing Davies–Bouldin scores at a given level, we expect the opposite. We did not test these hypotheses, as dendrogram purity is an external evaluation scheme, that is, it requires a clustering assignment to use for comparison, which is not available in unsupervised machine learning.

4.4 Results for the Quantum Version of BiDViT

We tested a prototype of BiDViT on a D-Wave 2000Q quantum annealer, a machine that has 2048 qubits and 5600 couplers. According to D-Wave Systems, the computer uses 128,000 Josephson junctions and was the most complex superconducting integrated circuit built to date when introduced in January of 2017 [7].

Before solving the QUBO problems, we applied preprocessing techniques, reducing their size and difficulty [11]. This proved effective and eliminated a great many variables. In most cases, we observed a size reduction of over 60%.

For the quantum version of BiDViT, we observed higher-quality solutions and a significant speed-up for BiDViT, when compared to common clustering methods. Both observations are based on results shown in Fig. 6.

However, the heuristic version of BiDViT and the common clustering algorithms were executed on a classical device that has a limited computational capacity, whereas the D-Wave 2000Q is a highly specialized device. Running these algorithms on a high-performance computer might lead to an equivalent degree of speed-up.

5 Conclusion

We have developed an efficient algorithm capable of performing extreme clustering. Our complexity analysis and numerical experiments show that if the dataset cardinality and the desired number of clusters are both large, the runtime of BiDViT is at least an order of magnitude faster than that of classical algorithms, while yielding a solution of comparable quality. With advances in quantum annealing hardware, one can expect further speed-ups in our algorithm and size of dataset that can be processed.

Independent of BiDViT, our coarsening method, based on identifying an $\varepsilon$-dense, $\varepsilon$-separated subset, is valuable in its own right—it is a novel approach to clustering which is not limited solely to the extreme clustering domain.

Further investigation of our coarsening approach is justified, as we have identified a domain for the radius of interest (in Thm. 2.4) such that, under a separability assumption, every solution to P0 (i.e., every maximum weighted $\varepsilon$-separated subset) yields the optimal clustering assignment.

Acknowledgements

We thank Saeid Allahdadian, Nick Condé, Daniel Crawford, and Austin Wallace for contributing to an earlier version of the algorithm. We thank Maliheh Aramon, Pooja Pandey, and Brad Woods for helpful discussions on optimization theory. The implementation of the QUBO preprocessing techniques was performed jointly with Brad Woods and Nick Condé. Inderpreet Singh contributed to the figure on image quantization. Victoria Wong assisted with graphical editing of two figures. Partial funding for this work was provided by the Mitacs Accelarate internship initiative.

Supplementary Information to the Paper, “A Quantum Annealing-Based Approach to Extreme Clustering”

Appendix A: An Alternative Coarsening Method

In certain situations, a user might not want the approximating set to be $\varepsilon$-separated but instead might be interested in finding an $\varepsilon$-dense subset with a minimum number of elements, or, more generally, the minimum cost for some cost function $c:X\rightarrow\mathbb{N}$. Finding such a set can be realized in a very similar way to the quantum method of BiDViT. In fact, the only modifications needed would be to Sec. 3.2 in the paper, where we introduce the concept of chunk coarsening.

Let $P=\{x^{(1)},\ldots,x^{(n)}\}$, and let $N^{(\varepsilon)}$ and $s_{i}$, for $i=1,\ldots,n$, be defined as in Sec. 3.2. Analogously to the weight vector $w$, we define a cost vector $c$ by $c_{i}=c(x^{(i)})$ for each $x^{(i)}\in P$. The problem of finding an $\varepsilon$-dense subset $S\subseteq P$ of minimum cost can then be expressed as follows:

[TABLE]

The constraints in P7 enforce the condition that for each solution (corresponding to a subset), every point in $P$ is represented by at least one of the points from the selected subset. The subset will not necessarily be $\varepsilon$-separated, but it will be $\varepsilon$-dense.

In the same way that finding an $\varepsilon$-separated subset of maximum weight corresponds to the MWIS problem, finding an $\varepsilon$-dense subset of minimum cost corresponds to the minimum weighted dominating set (MWDS) problem, which is equivalent to a weighted version of the minimal set covering (MSC) problem. Consider a set $U$ and subsets $S_{j}\subseteq U$ and $j\in J$, where $J$ is some set of indices, such that $U=\bigcup_{j\in J}S_{j}$. The MSC problem then consists of finding a subset $J_{0}\subseteq J$ such that the property $U\subseteq\bigcup_{j\in J_{0}}S_{j}$ is satisfied, and $J_{0}$ is of minimum cardinality with respect to this property. For example, if $U=\{a,b,c,d,e\}$, $S_{1}=\{a,c\}$, $S_{2}=\{a,d\}$, and $S_{3}=\{b,d,e\}$, then the solution to the MSC problem is given by $J_{0}=\{1,3\}$, as none of the subsets cover $U$, but the union $S_{1}\cup S_{3}$ does. The general MSC problem is known to be NP-hard [17]. By defining $S_{j}=B(x^{(j)},\varepsilon)\cap P$ for $j=1,\ldots,n$ in the above setting, one can see that we have solved a weighted version of the MSC problem.

To transform P7 into a QUBO problem, we convert the inequality constraints to equality constraints by adding integer slack variables. Note that the $i$-th constraint is satisfied if and only if there exists some $\xi_{i}\in\mathbb{N}_{0}$ such that $\sum_{j=1}^{n}N_{ij}^{(\varepsilon)}s_{j}-1=\xi_{i}$. In fact, given that $s\in\{0,1\}^{n}$, we can see that the $\xi_{i}$ must satisfy the bounds

[TABLE]

Thus, by dualizing the equality constraints, P7 can be expressed as a QUBO problem

[TABLE]

We will now describe how substituting a binary encoding for each of the $\xi_{i}$, for $i=1,\ldots,n$, in P8 yields the desired QUBO formulation. For each $i=1,\ldots,n$, the $(N^{(\varepsilon)}\mathbf{1})_{i}$ possible states of $\xi_{i}$ can be encoded by $\lfloor m_{i}\rfloor+1$ binary variables $b^{(i)}_{0},\ldots,b^{(i)}_{\lfloor m_{i}\rfloor}$, where $m_{i}=\log_{2}(N^{(\varepsilon)}\mathbf{1})_{i}$. The encoding has the form

[TABLE]

where $\gamma_{k}^{(i)}\in\mathbb{N}$ are fixed coefficients that depend solely on the bounds of 13. If we were to select $\gamma_{k}^{(i)}=2^{k}$ for $k=0,\ldots,\lfloor m_{i}\rfloor$, then, if $m_{i}\notin\mathbb{N}$, $\xi_{i}$ could assume states that do not satisfy these bounds. We can avoid this situation by manipulating the coefficient $\gamma^{(i)}_{\lfloor m_{i}\rfloor}$ of the final bit $b_{\lfloor m_{i}\rfloor}^{(i)}$ such that $\sum_{k=0}^{\lfloor m_{i}\rfloor-1}2^{k}+\gamma^{(i)}_{\lfloor m_{i}\rfloor}=(N^{(\varepsilon)}\mathbf{1})_{i}-1$. This may lead to a situation where there are multiple valid encodings for the same integer, but it will always hold that

[TABLE]

where $\gamma^{(i)}_{k}=2^{k}$ for $k<\lfloor m_{i}\rfloor$. Substituting the binary encoding into P8 yields the following QUBO formulation:

[TABLE]

One can show that the solution set of this QUBO problem is equivalent to the one for P7 for $\lambda>n\|c\|_{\infty}$. We have not investigated whether this bound is sharp. Note that our QUBO formulation is similar to the one described in [26], but uses a different encoding.

The number of binary variables in the QUBO formulation of this problem depends on the binary encoding of $\xi$. If the vertex degree in $G^{\varepsilon}$ is uniformly bounded from above by a constant $\eta>0$, then each $\xi_{i}$ can be encoded with fewer than $\lfloor\log_{2}(\eta)\rfloor$ binary variables. Therefore, the number of variables in the QUBO polynomial will be at most $n(1+\lfloor\log_{2}(\eta)\rfloor)$. In the worst case, that is, when there is vertex that is a neighbour of every other vertex, the polynomial would still comprise fewer than $n(1+\lfloor\log_{2}(n)\rfloor)$ variables.