Solving large Maximum Clique problems on a quantum annealer

TL;DR

This paper explores methods to decompose large Maximum Clique problems into smaller subproblems suitable for quantum annealers, employing heuristics and graph analysis to improve solution efficiency.

Contribution

It introduces novel decomposition and pruning techniques tailored for quantum annealing to handle large-scale Maximum Clique problems.

Findings

Effective decomposition reduces problem size for quantum annealers

Pruning strategies eliminate non-contributing subproblems

Enhanced solution quality with reduced computational effort

Abstract

Commercial quantum annealers from D-Wave Systems can find high quality solutions of quadratic unconstrained binary optimization problems that can be embedded onto its hardware. However, even though such devices currently offer up to 2048 qubits, due to limitations on the connectivity of those qubits, the size of problems that can typically be solved is rather small (around 65 variables). This limitation poses a problem for using D-Wave machines to solve application-relevant problems, which can have thousands of variables. For the important Maximum Clique problem, this article investigates methods for decomposing larger problem instances into smaller ones, which can subsequently be solved on D-Wave. During the decomposition, we aim to prune as many generated subproblems that do not contribute to the solution as possible, in order to reduce the computational complexity. The reduction methods presented in this article include upper and lower bound heuristics in conjunction with graph decomposition, vertex and edge extraction, and persistency analysis.

Tables

Table 1. Table 1: TTS times as a function of the intended graph density (first row). The time t run subscript 𝑡 run t_{\text{run}} is the sum of averages for embedding ( T embedding subscript 𝑇 embedding T_{\text{embedding}} ) and QPU solve time ( T QPU subscript 𝑇 QPU T_{\text{QPU}} ).

density	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
p	0.007	0.003	0.008	0.015	0.020	0.003	0.003	0.007	0.012
$t_{\text{run}}$	1.938	2.057	2.317	3.382	6.736	6.981	4.975	4.333	3.429
TTS	0.127	0.296	0.141	0.105	0.152	0.948	0.746	0.273	0.137

Table 2. Table 2: Total solution time in seconds for DIMACS graphs based on a single run.

DIMACS graph	No. vertices	No. edges	$\omega \left(G\right)$	DBK algorithm	Chapuis et al. (2017)
johnson16-2-4	120	5460	8	3217	3217
keller4	171	9435	11	682	1310
p-hat300-1	300	10933	8	16	20
p-hat300-2	300	21928	25	794	2781
p-hat500-1	500	4459	9	61	77
p-hat700-1	700	60999	11	175	643
brock200-2	200	9876	12	72	295
brock200-3	200	12048	15	1185	1640
brock200-4	200	13089	17	2168	14735
hamming6-2	64	1824	32	4	17
hamming8-4	256	20864	16	5179	18253

Equations

$$H=\sum _{i\in V}{a}_i{x}_j+\sum _{(i,j)\in E}{a}_{ij}{x}_i{x}_j,$$

$$H=-A\sum _{i\in V}{x}_i+B\sum _{(i,j)\in \overline{E}}{x}_i{x}_j,$$

$$\text{TTS}=(T_{\text{embedding}}+T_{\text{QPU}})\frac{\log (0.01)}{\log (1-p)},$$

Full text

Keywords: Branch-and-bound; Decomposition; D-Wave; Graph algorithms; Maximum Clique; Optimization; Quantum annealing.

1 Introduction

Quantum annealers such as the commercially available ones manufactured by D-Wave (D-Wave Systems,, 2000) aim to solve quadratic unconstrained binary optimization (QUBO) problems by looking for a minimum-energy state of a quantum system implemented on the hardware. Notably, due to the stochastic nature of annealing, the approximations of the global minimum returned by D-Wave are not always identical, while the solving time is independent of the problem instance. The class of QUBO problems includes many well-known NP-hard problems (Lucas,, 2014), for which no classical polynomial algorithms are known. In this article, we focus on the NP-hard problem of finding a maximum clique in a graph.

The D-Wave quantum processing unit is designed to minimize the Hamiltonian

[TABLE]

where $V=\{1,\ldots,n\}$, $\{x_{i}~{}|~{}i\in V\}$ is a set of $n$ binary variables, and $E\subseteq V\times V$ describes quadratic interactions between them. The coefficients $a_{i}$ are the weights of the linear terms, and $a_{ij}$ are the weights of the quadratic terms (couplers). If in (1) we have $x_{i}\in\{0,1\}$ for all $i\in V$, the Hamiltonian is called a QUBO, whereas if $x_{i}\in\{-1,+1\}$ for all $i\in V$, it is called an Ising Hamiltonian. The Hamiltonian for the maximum clique problem we consider is a QUBO.

The maximum clique problem (MaxClique) is defined as the problem of finding a largest clique (a completely connected subgraph) in a graph $G=(V,E)$. The clique number is the size of the maximum clique. In Chapuis et al., (2017), the QUBO formulation of the maximum clique problem is given as

[TABLE]

provided $B/A$ is sufficiently large, where the choice $A=1$ and $B=2$ ensures equivalence of the maximum clique problem and the Hamiltonian in (2), and $\overline{E}$ is the edge set of the complement graph of $G$.

The D-Wave annealer has a fixed number of operational qubits which are arranged in a certain graph structure, consisting of a lattice of bipartite cells, on a physical chip Chapuis et al., (2017). D-Wave 2X has 1152 physical qubits, while the newest D-Wave 2000Q model has up to 2048 physical qubits. For many NP-hard problems, however, the connectivity structure of a QUBO (that is, the graph defined by its non-zero coefficients $a_{ij}$) does not match the specific connectivity of the qubits on the D-Wave chip. In particular, this is the case for general instances of MaxClique. This problem can be alleviated by computing a minor embedding of a complete graph onto the D-Wave architecture, thus guaranteeing that any QUBO connectivity can be mapped onto the D-Wave chip. In such an embedding, a connected set of physical qubits is identified to act as one logical qubit, which further reduces the number of available qubits (D-Wave,, 2016). Even if the connectivity structure required to solve a QUBO matches the one of D-Wave, the scaling of problem instances quickly exhausts the capacity of the D-Wave chip. In both cases, classical preprocessing can attempt to split up large problem instances into smaller ones, which can be solved on D-Wave. General purpose approaches for QUBO decomposition exist (Boros and Hammer,, 2002), even though they are not guaranteed to always work for arbitrary QUBO instances.

Special purpose techniques exist for decomposing certain NP-hard problems. This paper investigates a decomposition algorithm for finding a maximum clique in a graph (Chapuis et al.,, 2017), where the problem reduction does not result in a loss of information regarding the maximum clique of the original problem. The decomposition algorithm works by recursively splitting a graph at each level in such a way that the maximal clique is guaranteed to be contained in either of the two generated subgraphs. To reduce the computational complexity, we attempt in this article to identify and exclude those generated subgraphs that cannot contain a maximum clique. This is achieved by computing bounds on the maximal clique size at each iteration. The purpose of this article is to investigate several such techniques with the aim of demonstrating the viability of solving larger optimization problems with D-Wave than previously possible.

Decomposition has already been suggested to solve several NP-hard problems in Tarjan, (1985), and since used for NP-hard problems such as domination-type problems and graph coloring in Rao, (2008). For the Maximum Independent Set problem, an equivalent formulation of MaxClique, polynomial-time algorithms are known for many special graph classes (Dabrowski et al.,, 2011; Grötschel et al.,, 1988; Minty,, 1980), including several algorithms relying on graph decomposition (Giakoumakis and Vanherpe,, 1997; Courcelle et al.,, 2000).

The article is structured as follows. Section 2 considers the problem of decomposing a MaxClique instance into smaller subproblems which fit D-Wave. The efficiency of the decomposition algorithm relies on the computation of upper and lower bounds on the clique size of the subgraph it generates, and a variety of techniques for computing such bounds is presented in section 3, including subgraph extractions methods. Section 4 contains an experimental analysis of the sensitivity of our decomposition approach on the chosen subgraph extraction techniques and on the bounds of the clique size, as well as an analysis of the effectiveness of these techniques for very sparse and very dense graphs compared to earlier work. The article concludes with a discussion in section 5.

2 Maximum Clique Decomposition

We use the CH-partitioning introduced in Djidjev et al., (2015), see also Chapuis et al., (2017), in order to split a large input graph $G=(V,E)$ into smaller subgraphs on which a maximum clique is found. We next briefly review the method.

2.1 CH-partitioning

The CH-partitioning algorithm works as follows. We iteratively select a single vertex $v\in V$ and extract the subgraph $G_{v}$ containing all neighbors of $v$ and all edges between them. Afterwards, $v$ and all edges adjacent to $v$ are removed from $G$, thus creating another new graph denoted as $G^{\prime}$. This is visualized in figure 1. The clique number $\omega(G)$ of $G$ is equal to $\min(\omega(G_{v})+1,\omega(G^{\prime}))$.

Before decomposing $G^{\prime}$ and $G_{v}$ further in a similar fashion, our algorithm computes lower and upper bounds on the maximum clique sizes contained in $G^{\prime}$ and $G_{v}$. This is helpful for two reasons. First, if the current best lower bound is less (worse) than the lower bound of a newly generated subgraph, we can update the current best lower bound accordingly. Second, if the upper bound for $G^{\prime}$ or $G_{v}$ is less than the best current lower bound for the maximum clique, we can prune that subgraph and the branch of the decomposition tree associated with it. Our aim is to prune as many generated subproblems as possible, and section 3 highlights the approaches we explore.

The algorithm stops splitting any generated subgraph as soon as it can be embedded on the D-Wave hardware, meaning if the subgraph size is at most $46$ ($65$) vertices in case of D-Wave 2X (D-Wave 2000Q). The solution time of problems that fit the D-Wave hardware is essentially constant (anneals only without post-processing), so we are primarily concerned with the number of generated subproblems and the effectiveness of our pruning approaches. However, we do take into account that the success rate of annealing may depend on the graph structure and we compute different times-to-solutions for different graph types.

The proposed algorithm can be characterized as a branch-and-bound (B&B) algorithm, a class of algorithms often used to compute a solution of NP-hard optimization problems. As such, B&B has been previously used in exact classical algorithms for solving MaxClique, e.g., Pardalos and Rodgers, (1992); Carmo and Züge, (2012). However, these algorithms explore the search space in a quite different way than the one studied in this paper, as they work bottom up, starting with a single-vertex clique and expanding it as much as possible. In comparison, our algorithm is top-down, starting with the original graph and splitting it down to smaller graphs, which are examined for a maximum clique. Such structure allows for the generation of a set of smaller size MaxClique problems which fit in and are solvable on the quantum processing unit directly.

2.2 Vertex Choice

The algorithm of section 2.1 offers one more tuning parameter: the procedure for selecting the vertex $v$ that is used in each iteration to split the current graph instance $G$ into two new graphs. Possible choices include:

1. (i)

a vertex $v$ of lowest degree. 2. (ii)

a vertex $v$ of median degree. 3. (iii)

a vertex $v$ chosen at random. 4. (iv)

a vertex $v$ of highest degree. 5. (v)

a vertex $v$ in the $k$-core of lowest degree. The $k$-core of a graph $G=(V,E)$ is a maximal subgraph of $G$ in which every vertex has a degree of at least $k$. Specifically, we iteratively extract the $k$-core of sizes $k\in\{1,\ldots,|G|\}$ and stop at the first $k$ leading to a $k$-core which is not the original graph. We then choose $v$ as a vertex which was removed in that $k$-core reduction step. 6. (vi)

a lowest degree vertex $v$ whose extracted subgraph $G_{v}$, defined in section 2.1, is of lowest density.

In any of the above cases, if multiple vertices satisfy the selection criterion, the vertex $v$ which is extracted is chosen at random. The above vertex selection approaches are experimentally explored in section 4.1.

3 Pruning of subproblems

Next we describe two pruning approaches we use in order to reduce the number of subgraphs generated in the algorithm of section 2. Section 3.1 describes bounds which can be used to prune entire subgraphs. Section 3.2 describes reduction approaches which allow to reduce the size of generated subgraphs, and sometimes entirely remove subgraphs, before decomposing them further.

3.1 Bounds

We consider two approaches to compute upper bounds, and two approaches to compute lower bounds, on the maximum clique size $\omega(G)$ of each subgraph $G$ generated in the graph decomposition algorithm of section 2. As upper bounds we consider the following:

1. (i)

We use a greedy search heuristic to find an upper bound on the chromatic number of the graph. The chromatic number of a graph $G=(V,E)$, denoted as $\chi(G)$, is the minimum number of colors needed to color each vertex of $G$ such that no edge connects two vertices of the same color. Since each vertex in a clique must have a distinct color we know that $\chi(G)\geq\omega(G)$ (Knuth,, 1993), thus the chromatic number of $G$ is an upper bound on $\omega(G)$. Computing the chromatic number is NP-hard, so its exact computation would be intractable, but there are much better heuristics for its approximation compared with the ones for the clique number. Therefore, a greedy search heuristic for the chromatic number could provide an easily computable upper bound for the maximum clique number. In the simulations of section 4, we use two greedy algorithms: the greedy coloring algorithm greedy_color of the python module Networkx (Hagberg et al.,, 2008) as well as the upper bound algorithm of (Budinich,, 2003, eq. (4)) for graphs with densities larger than $0.8$. 2. (ii)

The Lovász number of the complement of a graph $G$, denoted as $\theta(\overline{G})$, is defined in (Lovász,, 1979) as an upper bound on the Shannon capacity of $\overline{G}$, and is sandwiched in-between the chromatic number $\chi(G)$ and the clique number $\omega(G)$, thus $\chi(G)\geq\theta(\overline{G})\geq\omega(G)$ (Knuth,, 1993). Therefore, $\omega(\overline{G})$ might improve upon the bound already given by $\chi(G)$, and cannot be worse. Although the Lovász number can be calculated in polynomial time (Knuth,, 1993), existing algorithms are not very scalable (due to large constants) since they involve solving a semidefinite programming problem. The python implementation to compute the Lovász number, adapted from (Stahlke,, 2013), becomes prohibitively slow for larger graphs: we only use it for subgraphs $|G|\leq 60$, and otherwise revert to the greedy search heuristic outlined above.

Likewise, we consider the following two lower bounds:

1. (i)

Any maximum clique heuristic can be used to provide a lower bound on the clique number, provided the result returned by the heuristic is a clique. We use the Fast Maximum Clique Finder of Pattabiraman et al., (2015, 2018) in heuristic mode. 2. (ii)

At any point in the decomposition tree, a lower bound on one of the current subbranches can be obtained simply by solving the other subbranch. This gives a lower bound for the maximum clique size in all nodes (and leaves) of the yet unsolved branch.

It should be noted that there are many other easily computable upper and lower bounds for maximum clique, maximum independent set, and chromatic number that generally apply to any type of graph (Soto et al.,, 2011; Elphick and Wocjan,, 2018). However, in experiments we conducted (not reported in this article) those generally applicable bounds proved to be quite conservative, with the exception of the bound of (Budinich,, 2003, eq. (4)) which we employ for high densities.

3.2 Reduction algorithms

We implement three types of reduction algorithms. The first two work directly on the subgraphs, the last one works with the QUBO formulation of MaxClique. Due to the inherent similarity between the first two techniques, we group them together as one reduction algorithm.

1. (i)

The (vertex) $k$-core algorithm can reduce the number of vertices of the input graph in some cases, and the edge $k$-core algorithm (Chapuis et al.,, 2017; Batagelj and Zaversnik,, 2011) can reduce the number of edges.

The (vertex) $k$-core of a graph $G=(V,E)$ was defined in section 2.2 as the maximal subgraph of $G$ in which every vertex has a degree of at least $k$. Therefore, if a graph has a clique $C$ of size $k+1$, then this clique $C$ must be contained in the $k$-core of $G$ and all vertices outside of the $k$-core can be removed.

The edge $k$-core of a graph $G$ is defined in Chapuis et al., (2017). It is easily shown that for two vertices $v$, $w$ in a clique of size $c$, the intersection $N(v)\cap N(w)$ of the two neighbor lists $N(v)$ and $N(w)$ of $v$ and $w$ has size at least $c-1$. Denoting the current best lower bound on the clique size as $L$ we can therefore choose a random vertex $v$ and remove all edges $(v,e)$ satisfying $|N(v)\cap N(e)|<L-2$, since such edges cannot be part of a clique with size larger than $L$. 2. (ii)

A second approach aims to extract information from the QUBO formulation of MaxClique in (2). For a general QUBO (1), persistency analysis of (Boros and Hammer,, 2002) allows us to identify the values that certain variables must take in any minimum (called strong persistencies) or in at least one minimum (called weak persistencies), as well as relations among them. We apply (2) to any newly generated subgraph, compute a persistency analysis for the resulting QUBO and remove a vertex $v$ from the subgraph if its corresponding variable $x_{v}\in\{0,1\}$ (indicating if $v$ belongs to the maximum clique or not) could be allocated a value. Furthermore, the count of solved variables which are guaranteed to be in the maximum clique for a particular subgraph can improve the lower bound on the clique number for the overall algorithm.

We use the vertex and edge $k$-cores algorithms with $k$ being set to the current best lower bound value for the clique number. This allows us to prune entire subgraphs that cannot contain a clique of size larger than the best current one. The persistency analysis allows to remove single vertices and their adjacent edges from subgraphs, as well as improve the lower bound on the clique number.

4 Experimental analysis

We investigate the effectiveness of the vertex choice of section 2.2 and the pruning techniques of section 3, and we compare an algorithm combining both to existing work. For each graph density in the interval $[0.1,0.9]$, incremented in steps of $0.1$, we generate random graphs of size varying from $65-110$ vertices and decompose them using only the partitioning outlined in 2.1. We then generate random graphs based on each subgraph’s density and size, and solve them on D-Wave. As a performance measure we report the Time-To-Solution measure (defined below) and/or the number of subgraphs generated, both as a function of the graph density. The Time-To-Solution (TTS) measure is defined in (Mandrà and Katzgraber,, 2018, eq. (3)) as

[TABLE]

where $T_{\text{embedding}}$ and $T_{\text{QPU}}$ are the average times in seconds (for $10^{4}$ anneals per graph) to compute an embedding onto the D-Wave architecture and to perform annealing, and $p$ is proportion of anneals which correctly found the maximum clique.

Table 1 shows the proportion of anneals which found the maximum clique on D-Wave, the average times $t_{\text{run}}$ and the resulting TTS times for each of the graph densities considered. We fix the TTS times given in table 1 and use them to report runtimes in the remainder of section 4.

As test data, we use instances of $G(n,p)$ random graphs (Erdös and Rényi,, 1960) generated with density $p$ and fixed graph size $n=80$, which the exception of figure 2 which uses selected DIMACS benchmark graphs (Bader et al.,, 2013). All results are averaged over $20$ repetitions unless otherwise stated.

To calculate TTS for each DIMACS graph tested, we decomposed the graph being considered until we had less than or equal to 500 subgraphs. We then randomly selected 60 of those subgraphs as representatives, determined their sizes and densities, and for each we solved a random subgraph of similar size and density on D-Wave. We then applied eq. 3 to the resulting data.

4.1 Vertex Choice

Figure 2 shows the total solution time to decompose our test graphs (to the point that each subproblem reaches a size of at most $46$ vertices, thus making it possible to compute the maximum clique on D-Wave) as a function of the graph density. We compare the six different vertex selection choices delineated in section 2.2. The lowest degree selection strategy seems to lead to the best solution time among the six approaches, where the advantage of the lowest degree selection becomes more pronounced as the graph density increases. Similar results are observed for the number of subgraphs.

4.2 Bounds

Figure 3 (left) shows the number of subgraphs generated by every combination of the two lower and upper bounds outlined in section 3.1. Figure 3 (right) shows the speedup in runtime when using the upper bounds of section 3.1 compared to not using any bounds. The speedup is defined as $T_{0}/T_{1}$, where $T_{0}$ ($T_{1}$) is the total solution time without (with) bounds. The right figure shows that the Lovász number bound is not an efficient bound in terms of solution time due to the preprocessing CPU time required to calculate it. However, the Lovász number leads to the lowest number of generated subgraphs (left figure). Over most densities tested, the combination of fast heuristic bounds given by upper bound (i) and lower bound (ii) (see section 3.1) seem to give a reasonable trade-off between total solution time and subgraph count. We will use this choice in section 4.4.

4.3 Reduction algorithms

Figure 4 (right) shows the total runtime to decompose our test graphs as a function of the density $p$. Two methods are compared: the $k$-core reduction, and the persistency analysis (section 3.2). Figure 4 (left) shows that the $k$-core reduction algorithm can significantly reduce the number of subgraphs, and therefore the total solution time, over more densities than the persistency analysis can, and that the latter yields improvements in runtime only for high densities with $p>0.8$ (figure 4, right).

4.4 Comparisons

We denote the decomposition algorithm of section 2 which uses the choice of fast heuristic bounds of section 3.1 for pruning subgraphs and the $k$-core of section 3.2 as the DBK algorithm (Decomposition, Bounds, $k$-core). Figure 5 shows the speedup of the DBK algorithm over the decomposition algorithm presented in Chapuis et al., (2017). To be precise, figure 5 displays $T_{\text{Chapuis}}/T_{\text{DBK}}$, where $T_{\text{DBK}}$ is the time the DBK algorithm takes for a particular graph instance, and analogously for $T_{\text{Chapuis}}$. The figure shows that both algorithms are quite similar for graph densities in the interval $[0.4,0.8]$, however DBK considerably improves upon Chapuis et al., (2017) the sparser or denser the graphs are.

For DIMACS graphs displayed in table 2, we observe that the DBK algorithm can yield substantial improvements in runtime especially for large graphs, and that it never performs worse than Chapuis et al., (2017).

5 Discussion

This article investigates subgraph extraction and pruning techniques for a decomposition algorithm designed to compute maximum cliques of a graph. These reduction techniques aim at reducing the number of recursively generated subgraphs for which the maximum clique problem needs to be solved. The overall runtime of our algorithm depends on the user’s choice of method to solve the decomposed maximum clique subproblems at leaf level of the decomposition. Our analysis shows the following:

The decomposition is sensitive to the vertex choice used to extract subgraphs. Extraction of the subgraph induced by the lowest degree vertex seems to lead to the lowest overall runtime. 2. 2.

A variety of lower and upper bounds exist to estimate the maximal clique size and help prune subgraphs which cannot contain the maximum clique. Of those investigated, the Lovász bound is most effective in pruning subgraphs but also most computationally intensive, whereas a heuristic upper bound leads to a good trade-off between subgraph reduction and runtime. 3. 3.

The $k$-core reduction prunes subgraphs more effectively than an alternative approach based on persistency analysis. 4. 4.

Combining the decomposition algorithm with subgraph extraction and fast heuristic bounds on the clique size for pruning leads to an improved algorithm which, despite the additional computations in each recursion, is faster at computing maximum cliques than a previous algorithm of Chapuis et al., (2017), with a considerable speed-up for very sparse and very dense graphs.

Future research could investigate the viability of trade-off algorithms, where differing decomposition and simplification strategies are used in response to different problem graph types. Another future improvement concerns the Lovász number: our experiments show that its bounds yield the best reduction in subproblem count, but also come at the highest computational cost (see figure 3). Additionally, the impracticality of applying the adapted python implementation (Stahlke,, 2013) on large graphs was not ideal. A faster way to compute the Lovász number would render those bounds usable in practice.