Efficient Combinatorial Optimization Using Quantum Annealing
Hristo N. Djidjev Guillaume Chapuis Georg Hahn Guillaume Rizk
quant-ph




Efficient Combinatorial Optimization Using Quantum Annealing
    Hristo Djidjev (PI)
Los Alamos National Laboratory
Guillaume Chapuis
Los Alamos National Laboratory
Georg Hahn
Imperial College, London, UK
Guillaume Rizk
INRIA, Rennes, France
    November 8, 2016
====================================================================================================================================================================================







The recent availability of the first commercial quantum computers has provided a promising tool to tackle NP hard problems which can only be solved heuristically with present techniques. However, it is unclear if the current state of quantum computing already provides a quantum advantage over the current state of the art in classical computing. This article assesses the performance of the D-Wave 2X quantum annealer on two NP hard graph problems, in particular clique finding and graph partitioning. For this, we provide formulations as Qubo and Ising Hamiltonians suitable for the quantum annealer and compare a variety of quantum solvers (Sapi, QBSolv, QSage provided by D-Wave Sys, Inc.) to current classical algorithms (METIS, Simulated Annealing, third-party clique finding and graph splitting heuristics) on certain test sets of graphs. We demonstrate that for small graph instances, classical methods still outperform the quantum annealer in terms of computing time, even though the quality of the best solutions obtained is comparable. Nevertheless, due to the limited problem size which can be embedded on the D-Wave 2X chip, the aforementioned finding applies to most of problems of general nature solvable on the quantum annealer. For instances specifically designed to fit the D-Wave 2X architecture, we observe substantial speed-ups in computing time over classical approaches.







§ INTRODUCTION


The availability of the first commercial quantum computers <cit.> has led to many new challenges in the areas of computing, physics and algorithmic design. In particular, the new devices provide a promising tool to solve NP hard problems which can only be solved heuristically with present techniques. However, it is still disputed in the literature if the available quantum computers indeed provide an advantage over classical methods <cit.>.

This article investigates the performance of the D-Wave 2X quantum computer on two NP hard graph problems, maximal clique finding and (edge-cut as well as ch-) partitioning. The first problem aims to find the largest complete (i.e., fully connected) subset of vertices for a given graph <cit.>. The latter typically aims to divide a graph into a pre-specified number of components in such a way that the number of edges in-between all components is minimized <cit.>. Another variant of graph partitioning investigated in this article is core-halo (CH) partitioning by <cit.>.

A variety of classical heuristics are known to solve both problems in reasonable time with high accuracy. The device developed by D-Wave Systems, Inc., is a quantum annealer which follows a different approach. Briefly speaking, a quantum annealer is a hardware realization of the classical simulated annealing algorithm, with the exception that it makes use of quantum effects to explore an energy or objective function landscape. Similarly to the random moves proposed in simulated annealing, quantum annealers use quantum fluctuations to escape local minima and to find low-energy configurations of a physical system. The quantum bits (qubits) are realized via a series of superconducting loops on the D-Wave chip. Each loop encodes both a 0 and 1 bit at the same time through two superimposed currents in both clockwise and counter-clockwise directions until annealing is complete and the system turns classical <cit.>.

Due to the particular architecture of the D-Wave chip, the device is designed to minimize an unconstrained objective function consisting of a sum of linear and quadratic qubit contributions, weighted by given constants. To be precise, it aims to minimize the Hamiltonian

    ∑_i ∈ V a_i q_i + ∑_(i,j) ∈ E a_ij q_i q_j

for given a_i, a_ij∈ℝ, where V = {1,…,N} and E=V × V <cit.>. The advantage of a quantum computer for minimizing an Hamiltonian of the aforementioned form lies in the fact that, although classically an NP hard combinatorial problem, its mapping of the problem onto a physical system allows to read off solutions by observing how the system arranged itself in an (energy) minimum. Notably, the computing time is independent of the problem size and unlike existing digital computers, which are basically deterministic, quantum annealing is fundamentally probabilistic: Due to the experimental nature of the minimisation procedure, solutions are not always identical.

The article is structured as follows. Section <ref> starts by introducing the Qubo (quadratic unconstrained binary optimization) and Ising Hamiltonians and contains a note on the D-Wave architecture. Moreover, it describes the quantum tools (Sapi, QBsolv and QSage, see <cit.>) as well as the classical algorithms (simulated annealing <cit.>,  <cit.> and  <cit.>) we make use of in this study.

Section <ref> investigates the maximum clique problem. It describes its Qubo formulation, a graph splitting algorithm as well as the parameters we used to solve the maximum clique problem on D-Wave. Moreover, Section <ref> presents a comparison study of D-Wave to for small graphs as well as graphs specifically chosen to fit the D-Wave qubit architecture. The comparison involves classical methods (quantum adiabatic evolution of <cit.>, specialized simulated annealing of <cit.>), the  clique heuristic and exact solver <cit.> and ) as well as the aforementioned solvers provided by D-Wave.

Section <ref> focuses on two flavors of graph partitioning, called edge-cut and CH-partitioning. It derives Hamiltonians in both Qubo and Ising formulation and proves conditions on the Hamiltonian weights which guarantee that the global minimum of the Hamiltonian is attained at the optimal partitioning solution. Using random graphs, Section <ref> evaluates all the aforementioned classical and quantum solvers in terms of computing time as well as the quality of the partitionings computed, both on the edge-cut and on the CH-partitioning problem. We also compare the embedding time vs. the anneal time on the D-Wave device.

We conclude with a discussion of our results in Section <ref>, arguing that a considerable quantum advantage is observable in instances only which are tailored to the D-Wave chip, but not in problem instances of general form.

In the entire paper, we denote a graph as G=(V,E), consisting of a set of n vertices V = { 1,…,n } and a set of undirected edges E. The cardinality of a set V is denoted as |V|.






§ METHODS


This section introduces quadratic integer programs, the type of problem the D-Wave device is capable of solving. We provide two different formulations of such quadratic programs, and show how to convert between the two. Next, this section briefly presents three tools provided by D-Wave Inc. to submit quadratic programs to the quantum computer. To this end, we treat their input format, the solver used and its interface, as well as available post-processing methods and output format. The section concludes with a summary of classical methods we will use in the forthcoming sections to compare any D-Wave results.



 §.§ Ising and Qubo formulations


The D-Wave device is designed to solve a certain type of problem only, called quadratic integer programs. It operates on n variables q_i, i ∈{1,…,n}, and allows linear and quadratic interactions, summarized as

    min_q_1,…,q_n(∑_i=1^n a_i q_i + ∑_1≤ i<j≤ na_ij q_i q_j)

where the linear weights a_i ∈ and the quadratic coupler weights a_ij∈ are pre-specified by the user. The aforementioned quadratic minimization program exists in two flavors: In a Qubo formulation, in which all variables take binary values, i.e. q_i ∈{0,1}, and as an Ising formulation operating on q_i ∈{-1,+1}. In both formulations, the form of the quadratic program (<ref>) is unchanged.

As a convention, Qubo problems are often stated in the literature as a lower- (or upper-) diagonal matrix of weights ( Q_ij )_i,j, where Q_ii=a_i encode the linear weights and off-diagonal elements encode the quadratic couplers. Ising problems are sometimes given as a vector h=(h_1,…,h_n)=(a_1,…,a_n) of linear weights and a lower- (or upper-) diagonal matrix J_ij = Q_ij for i<j (i>j) with couplers.



 §.§ Conversion of Ising to Qubo and vice versa


Qubo and Ising quadratic integer programs are equivalent. They only differ in the values the variables q_i are allowed to take, and hence allow to be converted from one formulation to the other. This is done as follows.

Suppose we are given a Qubo problem with linear weights Q_ii, i ∈{1,…,n}, and quadratic weights Q_ij, i<j. Then, an equivalent Ising problem is given by setting

    h_i    = Q_ii/2 + ∑_j=1^n Q_ij/4,
    
    J_ij   = Q_ij/4

for i ∈{1,…,n} and all i<j.

Likewise, any Ising problem given by linear weight vector h and coupler matrix J can be converted to a Qubo by defining

    Q_ii   = 2 ( h_i - ∑_j=1^n J_ij),
    
    Q_ij   = 4 J_ij,

for i ∈{1,…,n} and all i<j.



 §.§ Qubit architecture and chains of qubits


The D-Wave quantum computer operates on roughly 1000 qubits arranged on a particular graph architecture. The precise number of available qubits varies from chip to chip due to manufacturing errors, meaning that slightly less of the envisaged qubits are usually operational in each device. The qubits are arranged in a chimera architecture on the chip consisting of 4 × 4 bipartite graph cells, interconnected on a lattice of 12  12 cells.

The particular architecture of the qubits implies two important consequences: First, the chip design only allows for pairwise contributions of two qubits. Together with the linear contributions, this motivates the form of the objective function as a quadratic integer program in (<ref>) which can be solved with D-Wave.

Second, although in theory, (<ref>) includes pairwise interactions between any two qubits q_i and q_j, this is not the case in practice. Since any physical qubit on the D-Wave chip possesses edges only within its 4 × 4 bipartite graph cell, as well as possibly some further connections to qubits of adjacent cells, pairwise interactions are actually very limited. Therefore, in practice, a pairwise interaction specified in the Qubo or Ising input problem has to be re-routed through various qubits in case the two qubits involved are not already neighbors on the physical chip. These additional qubits used to connect two bits q_i and q_j which are actually not neighbors on the physical chip are called a chain.

The existence of chains has two vital consequences which will play an important role in Sections <ref> and <ref>. On the one hand, the need for chains uses up qubits which would otherwise be available in the quadratic program, meaning that for problems containing many pairwise interactions, less than the around 1000 qubits are available. In the worst-case of a complete graph, only 48 qubits can be used in (<ref>) when arbitrary pairwise couplers are desired. On the other hand, due to the probabilistic nature of the quantum annealing process, solutions returned by D-Wave are not always identical. This is a problem for all chains used in a particular instance, since all qubits in a chain encode the same single bit in (<ref>) by construction, and should hence have a consistent value in any solution. This phenomenon is called a broken chain. It is not clear which value a particular qubit shall be assigned if its chain is broken, making post-processing techniques necessary. Naturally, chains can be ensured to not break by assigning them higher coupler weights than in the actual quadratic program the user is interested in. In this case, however, the quantum annealer with be implicitly tuned towards ensuring that all chains are satisfied (i.e. possess one consistent value), at the expense of violating some of the couplers of the actual quadratic program, thus leading to infeasible solutions. The trade-off between ensuring satisfied chains and putting sufficient weight on the actual quadratic program the user wishes to solve requires extensive further research.



 §.§ D-Wave solvers


D-Wave Inc. provides several solvers to implement given instances of (<ref>), to submit them to the quantum computer, apply necessary post-processing steps to the solution it returns and to format the output. This section highlights the interfaces to the D-Wave computer used in this article.



  §.§.§ Sapi

Sapi interfaces for the programming languages C++ and Python are provided by D-Wave Inc. These interfaces allow to establish an internet connection to a D-Wave server, to submit quadratic programs in Ising and in Qubo format, and to post-process the solution.

Most importantly, the sapi interface allows the user to take own control over the D-Wave parameters: It allows to compute individual embeddings of the Ising or Qubo to the chimera physical graph structure, to choose the number of anneals the quantum device is supposed to perform, or to specify the type of post-processing. The quadratic program itself is provided as a vector (containing linear weights) and matrix (containing quadratic weights) in appropriate data structures, solutions are returned as vectors which can be evaluated or processed further in any way the user sees fit. Although there are sapi routines for solving both Ising and Qubo formulations, we noticed that the Ising routines seem to work a lot better. This is not an issue since sapi routines to convert between the equivalent Ising and Qubo formulations are provided.

Since the user has full control over the graph embedding, quadratic programs can be provided which are designed to perfectly fit the chimera architecture. Arbitrary graphs are limited to 48 vertices (qubits), the maximal graph size that can be embedded onto the chimera graph even when all pairwise interactions are present.

In summary, the sapi interface provides the lowest level control of the D-Wave system and is hence the best candidate for working directly on the quantum annealer.



  §.§.§ QBsolv

QBsolv is a high-level tool that aims to take care of the D-Wave parameters and embedding issues automatically. It is not a programming interface but rather given by stand-alone binaries for selected computer architectures.

QBsolv solves Qubo problems only. The user provides a Qubo in a specific file format which lists linear weights first and quadratic couplers afterwards, where for each quadratic coupler Q_ij, the convention i<j has to be satisfied.

Notably, instances submitted to QBsolv are not limited in size. They are automatically split up into subproblems submitted to D-Wave, and an extensive tabu search is applied to post-process all D-Wave solutions. QBsolv further allows to specify certain parameters such as the number of individual solution attempts, the subproblem size used to split up instances which do not fit onto D-Wave, a target value of the Hamiltonian that is sought or a timeout parameter.

The solution is either displayed in the terminal or written into an output file from which it has to be re-read for further processing.



  §.§.§ QSage

QSage is another programming interface which allows to establish a connection to the D-Wave server within textitC++ and Python code. In contrast to Sapi or QBsolv, it is a blackbox solver which does not require a Qubo or Ising formulation.

Instead, QSage is able to minimize any one dimensional function operating on a binary input string of arbitrary size. To this end, QSage internally samples random points from the blackbox objective and uses the resulting function values to interpolate a quadratic program which is then submitted to D-Wave. The minimum returned by D-Wave is then used as a proxy for the minimum of the objective function and thus as a starting point for a tabu search or for further expansions of the objective function. Although QSage does not require the objective function to be smooth, naturally, best results are obtained in practice for smooth landscapes.



 §.§ Classical solvers


Apart from the tools provided by D-Wave Inc., we employ classical solvers in our comparison. These are described in this section.



  §.§.§ Simulated annealing with random starting point


Simulated annealing (SA) is a widely used optimization algorithm <cit.>. It provides an effective way to minimize (maximize) virtually any objective function using pointwise evaluations only.

SA iteratively minimizes an objective function by proposing a sequence of random moves (modifications) to an existing solution. These moves are drawn randomly according to a user-specified distribution. In order to give the algorithm the ability to escape local minima, modifications which do not further minimize the currently found best solution are not always rejected, but accepted with a certain acceptance probability. This acceptance probability is proportional to the magnitude of the (unfavorable) increase in the objective function and antiproportional to the runtime. The latter gives SA the property to explore the objective function landscape at the beginning of a run while gradually turning into a steepest descent algoritm over time. In order to make the acceptance probability antiproportionally depend on the runtime, a so-called temperature function is employed. This is any decreasing function to zero. SA can be run with a variety of stopping criteria, for instance with a fixed number of iterations (as in our case), until no better solution has been found for a certain number of iterations, or when the temperature becomes vanishingly small.

Algorithm <ref> shows an implementation of SA adapted from <cit.> for solving graph partitioning instances given by an Ising Hamiltonian (see Section <ref>). It uses the common temperature function T(t) = log(1+t) <cit.>.




  §.§.§ METIS and Simulated Annealing


 <cit.> is a popular heuristic multilevel algorithm to perform graph partitioning based on a three-phase approach:

  
  * The input graph is coarsened by generating a sequence G_0, G_1, …, G_n of graphs starting from the original graph G=G_0 and ending with a suitably small graph G_n (typically less than 100 vertices). The sequence has the the property that G_i contains less nodes than G_j whenever i<j. To coarsen the graph, several vertices are combined during the coarsening phase into a multinode of the next graph in the sequence.
  
  * G_n is partitioned using some other algorithm of choice, for instance spectral bisection or k-way partitioning.
  
  * The partition is projected back from G_n to G_0 through G_n-1,…,G_1. As each of the finer graphs during the uncoarsening phase contains more degrees of freedom than the multinode graph, a refinement algorithm such as the one of Fiduccia-Mattheyses <cit.> is used to enhance the partitioning after each projection.

 allows the user to change several parameters, for instance the size of the subgraphs G_n, the coarsening or the partitioning algorithms.

We employ  in two ways. First, to calculate edge-cut partitionings of given graphs (Sections <ref> and <ref>) and to use its partitions as a starting point for core-halo (CH-) partitioning <cit.> investigated in Section <ref>.

For CH-partitioning, <cit.> provide an algorithm based on Simulated Annealing capable of refining a given partitioning in order to optimize it further with respect to a different metric. Whenever a CH-partitioning is sought in Section <ref>, we thus employ the Simulated Annealing approach of <cit.> with initial solution provided by .



  §.§.§ Gurobi


 <cit.> is a mathematical programming solver for linear programs, mixed-integer linear and quadratic programs, as well as certain quadratic programs.

We employ  to solve given Qubo problems (Ising problems can be solved as well, nevertheless  explicitly allows to restrict the range of variables to binary inputs, making it particularly suitably for Qubo instances). To this end, we programmed a simple converter which takes a given Qubo and outputs a quadratic program in  format.






§ THE MAXIMUM CLIQUE PROBLEM






 §.§ Introduction

The maximum clique problem is a famous NP-Hard problem, with applications in social networks, bioinformatics, data mining, and many other fields. Its formulation is quite simple, given an undirected graph  G = (V, E), a clique is a subset S of the vertices forming a complete subgraph, i.e. where every two vertices of S are connected by an edge in G. The clique size is the number of vertices in S, and the maximum clique problem is to find the largest possible clique in G. 

A previous study investigated the use of quantum adiabatic evolution to solve the clique problem <cit.>. However, the system was not available at that time and the study was therefore limited to small graphs (n ≤ 18). Its conclusion was that their quantum algorithm appeared to require only a quadratic run time. Our objective is to verify that assumption on the real system, on much larger graphs if possible.



 §.§ QUBO formulation

A QUBO is written as 

    ∑_i≥ j Q_ijx_ix_j

with  x_i∈{0,1} and i ∈ [0;N]. The Q_ij with i ≠ j are the quadratic terms, and Q_ii are the linear terms. Paper from A.Lucas provides this QUBO formulation of the clique problem <cit.>:

    H =  A ( K- ∑_v x_v )^2 + B [ K(K-1)/2 - ∑_(uv) ∈ E x_ux_v].


The ground state will be H=0 when a clique of size K exists. Extra variables are then added to find the largest clique in the graph. However, this formulation requires N=|V| variables on a complete graph, because of the ( K- ∑_v x_v )^2 term that will generate every posible quadratic term in the QUBO formula. We prefer a simpler formulation, derived from the maximum independent set (MIS) problem.

An independent set of a graph G is a set of vertices, where every two vertices are not connected by an edge in G. It is easy to see that an independent set of G = (V, E) is a clique in the complement graph H = (V, E), where E is the set of edges not present in G. Therefore looking for the maximum clique in G or the maximum independent set in H is the same problem. The QUBO formulation of the MIS problem on a graph H is straightforward, let x_i=1 if vertex i is in the independent set, and x_i=0 otherwise. We should have a penalty when two vertices in this set are connected by an edge in the graph, therefore we need some positive quadratic term M over all existing edges of the graph H. We want the largest independent set, therefore a negative weigth L on the linear terms. The penalty for an edge should be greater than the bonus for a vertex, so that it is detrimental to add a vertex with a single edge to the set, so for example L=-1 and M=2.

The QUBO for the MIS problem of a graph G = (V, E) is therefore:

    -∑_i ∈ V x_i + 2∑_(i,j) ∈ E x_ix_j.


This formulation requires |V| variables and |E| quadratic terms. Note that the ising graph (the set of edges and vertices with non zero weigth, also called the primal graph) required for this problem is the graph G itself.

When solving the maximum clique problem for a graph G_1 = (V,E) we need to complement the graph first to G_2= (V, E) and solve the MIS on G_2, thus requiring  |V| variables and |E| =  |V| (|V|-1)/2 - |E| quadratic terms. It will therefore be easier to solve the maximum clique problem on a very dense graph, that will generate a sparse complement graph.



 §.§ Methods

The system consists in a specific graph topology called the chimera graph. It is a lattice of 12x12 cells, where each cell is a 4x4 bipartite graph. can naturally solve ising problems where non-zero quadratic terms are represented by an edge in the chimera graph. For other problems, an embedding is needed, i.e a logical variable has to be represented by a chain of physical qubits, in order to increase the connectivity of the simulated variable. The largest complete graph that the can embed in theory has 1+4*12 = 49 vertices. It is slightly smaller (n=45) because of missing qubits coming from manufacturing defaults.

In our problem, the maximum number of connections needed may be up to a complete graph. The largest graph that we are sure to be able to solve on is therefore a 45-vertex graph. For larger graphs, it will depend on the edge density and topology.

To solve the MIS problem on any graph, we therefore require a divide-and-conquer strategy that first splits the input graph into smaller graphs of at most 45 vertices.



  §.§.§ Sapi parameters

Sapi stands for solver API, it is the lowest level of control we can have over the quantum annealer. We use here the sapi C client version 2.4.2. We first convert our QUBO problem to an ising problem. We use the pre-computed embedding of the full 45-vertex graph for the experiments in section <ref>. In section <ref> we do not need to compute an embedding. We use num_reads=500, postprocessing=SAPI_POSTPROCESS_OPTIMIZATION and provide the embedding to the “chains” parameters of the solver so that the post-processing is applied to the original problem without chains (with this parameter the solver first applies majority voting to resolve chains, before applying post-processing). Our experiments showed that this is the best way to deal with chains.



 §.§ Results




  §.§.§ Software solvers

Results from are compared to different solvers. A simulated annealing algorithm working on the ising problem (SA-ising), a simulated annealing algorithm specifically designed to solve the clique problem (SA-clique <cit.>), a software designed to find cliques in heuristic or exact mode (fmc <cit.>), and the gurobi solver. We also include in the comparison results coming from the Dwave post-processing heuristics only.



  
SA-ising
This is a simulated annealing algorithm working on an ising problem. The intial solution is a random solution, and a single move in the simulated algorihtm is the flip of one random bit.



  
SA-clique
We implememted the simulated algorithm designed to find cliques described by X.Gen et al. <cit.>. It is designed to find a clique of a specific size m, we therefore need to apply a binary search on top of it to find the maximum clique size. The main parameter is the α value controlling the temperature update at each step: T_n+1 = α T_n. They typically use α=0.9996. A value closer to 1 will yield a better solution but will increase computation time.



  
Fast Max-Clique Finder (fmc)
This algorithm is designed to find efficiently the maximum clique for a large sparse graph. It provides an exact and a heuristic search mode. We use the 1.1 version of the fmc software <cit.>.



  
Post-processing heuristics alone (PPHa)
The Dwave pipeline includes a post-processing step: first, if chains exist, a majority vote is applied to fix the broken chains. Then a local search is performed so that solutions become locally minimum (the raw solutions coming from dwave might not be locally minimum) <cit.>. For a given solution coming out of the pipeline, one might wonder what are the relative contributions of the system and of the post-processing step. For some small and simple problems, the post processing step alone might be able to find a good solution.

We try to answer this by applying post-processing step only, and see how it compares with the actual quantum annealing. However, this post-processing step runs on the dwave server and is not available separately. 

We therefore need some special procedure: we set a very high absolute chain strength (e.g. 1000 times greater than the largest weight in our ising problem), and turns on the “auto-scale” feature. Consequently, chains weights will be set to minimum value -1 and all other weights will be scaled down to 0. In this way, the quantum annealer will only satisfy the chains, and will not solve anything else. Each chain will not be connected to anything, and all linear terms will be zero, therefore each chain should get a random value -1 or +1. The post-processing step will therefore be fed with some random initial solution, and what we get at the end of the pipeline should be the result of post-processing step only. This will be referred to later as PPHa, post-processing heuristic alone.



  §.§.§ Results




  
Results on a small graph
We first analyze the behavior of the quantum annealer on one graph of size at most 45 vertices, i.e. the largest graph we are sure to be able to embed. 



We test our graph splitting method on random graphs with 500 vertices and increasing density. The density of the generated graph is controlled by the probability of the presence of an edge ranging from 0.1 to 0.4. Fig. <ref> shows the results of this experiment. We report here the number of solver calls required to solve the input graph -ie. the number of generated subgraphs- with respect to the edge presence probability. Each data point is the median value of ten runs with the standard deviation given as error bars. The number of solver calls follows an exponential curve with respect to the density of the input graph.



In a second experiment, we measure the impact of future generations of DWave systems on our clique finding approach for arbitrary large graphs. Assuming a similar chimera topology for future generations of DWave systems, doubling the number of available QBits will increase the size of the maximal complete subgraph that can be embedded by a factor √(2). Fig <ref> shows the evolution of the number of solver calls for increasing number of QBits. In this experiment the graph size is fixed at 500 vertices and the edge presence probability at 0.3. Each data point is the median of ten runs with standard deviation given as error bars. If we consider that the number of QBits doubles with each new generation, 7 generations of DWave machine are required to directly embed and solve an arbitrary 500 vertex graph.





  §.§.§ Artificial graph designed to fit chimera


The idea here is to try to see what is the speedup can achieve in the best possible scenario. This best possible scenario is a graph that fits well the chimera graph, i.e. a graph that requires a very few chains, and hopefully for which the MIS problem is not trivial.



  
Graph generation
We start with the actual chimera graph C = (V,E) of the system, in our case |V = 1100|  qubits (this may vary from system to system because of missing qubits).  This is the most dense 1100-vertex graph the system can solve. Note that this graph is not interesting for the MIS problem since it is bipartite, and one could argue that MIS on a bipartite graph is a trivial problem.

Consider now the graph C_1 obtained by contracting one random edge from C. An edge contraction is the removal of an edge v_1,v_2 and fusion of vertices v_1 amd v_2.  With 𝒩_1 and 𝒩_2 the set of neighbouring vertices from v_1 and v_2, the neighbours of the fused node are 𝒩_1∪𝒩_2 ∖{v_1,v_2}. C_1 has |V|-1 vertices and is called a minor of the chimera graph. By definition the ising formulation of the MIS problem[The primal graph of the ising MIS problem of a graph is the graph itself.] for this graph can be embedded: the embedding is given by the edge contraction. Moreover, if we add any edge to C_1, the resulting graph will not be embeddable into chimera, since C_1 already uses every possible qubits and edges of the chimera graph. We can therefore say that C_1 is one of the most dense graph of size |V| -1 than chimera can embed.

We can generalize after m random edge contractions: the resulting graph C_m will have |V| -m vertices, will be one of the most dense graph of size |V| -m chimera can embed, and will not be bipartite. This family of graphs C_m with 0<m<1100 is therefore a good candidate for the best case scenario experiment for the MIS problem: large graphs that can embed on chimera, and for which solution is not trivial. 
 


  
Experiment
We solve the MIS problem on the C_m family of graphs, using Sapi on the system, PPHa and the SA-ising software. We solve the equivalent maximum clique problem on the complement graphs using SA-clique and fmc. Figure <ref> shows the result.  We can see that for graph sizes up to 400, PPHa gets same result as Dwave. For these small graphs the problem is probably simple enough to be solved by the post-processing step alone. As expected, the simulated annealing designed specifically for clique is behaving better than SA-ising. fmc software is run in its heuristic mode. It is designed for large sparse graphs, and run here on very dense graphs, probably the reason of comparatively lower quality result. 

For large graphs (larger than 800 vertices), gives the best solution. Note that we do not know if that is the optimal solution.



  
Speedup
Since SA-clique seems to best the best candidate to compete with , we choose to compute the speedup relatively to SA-clique on this C_m family of graphs. The procedure is as follows: for each graph size, we run with 500 anneals and report the best solution. The runtime is the total qpu runtime for 500 anneals, i.e. approximately 0.15s. For SA-clique, we start with a low α parameter (i.e. a fast cooling schedule), and gradually increase alpha until SA-clique finds the same solution as . This alpha for which SA-clique finds the same solution gives us the SA-clique best execution time. The SA-clique algorithm is run on one CPU core of an Intel E8400 @ 3.00GHz.  Figure <ref> shows the speedup for different graph sizes (of the C_m family). For graph sizes lower than 200 vertices, is slower, but for larger grah sizes, it gets exponentially faster, up to a million speedup. 

Although this is the best possible case scenario for and not representative of what we would get with random graphs, this may be an evidence showing quantum behavior in the system.  is able to find very quickly a solution that is very difficult to obtain with all the classical solvers we tried.










§ GRAPH PARTITIONING


This section investigates a variety of methods for graph partitioning. Our focus is on two variations of the graph partitioning problem: Edge-cut partitioning and core-halo partitioning <cit.>. To this end, Section <ref> provides Hamiltonians in either Ising or Qubo formulation to carry out both types of partitioning. We then employ the methods outlined in Section <ref> to minimize those Hamiltonians on a variety of graphs in order to assess their performance in a simulation study (Section <ref>).



 §.§ Hamiltonian formulations




  §.§.§ Ising model for graph partitioning


<cit.> gives an Ising Hamiltonian to carry out standard edge-cut partitioning. However, this formulation is only suitable to partition a graph into two components. A generalisation to n components will be derived in the next section.

Edge-cut partitioning aims to divide the graph into partitions in such a way that the edge-cut, the number of edges connecting vertices belonging to different partitions, is minimized.

To this end, <cit.> considers an Ising spin s_v ∈{ +1,-1} per vertex v ∈ V indicating which of the two partitions v belongs to (one being the "+" partition and one being the "-" partition). He defines the Hamiltonian as H=A H_A + B H_B, where

    H_A    = ( ∑_i=1^n s_i )^2,
    
    H_B    = ∑_(u,v) ∈ E1-s_u s_v/2.

In the above, H_A is a penalty term enforcing equal balance of vertices in both partitions. Assuming an even number of vertices, H_A is zero if and only if the number of vertices in both partitions is equal. For an odd number of vertices, H_A ≥ 1. The term H_B is the objective function. It adds a penalty of 1 for any edge connecting two vertices which are not in the same partition and is thus equal to the edge-cut. Minimizing H will thus lead to a minimal edge-cut partitioning of the graph.

The weights A and B have to be chosen dependent on the input graph in such a way that it is never favorable to violate H_A in order to achieve a further decrease in the objective function H_B. To this end, it suffices to determine the minimal increase in the penalty for any given change while maximizing the decrease in the objective function. Consider the configuration of s_v, v ∈ V, at the global minimum (the equilibrium state) of H. For simplicity, assume without loss of generality that the number of vertices is equal and hence that H_A=0.

When changing the spin of one s_v, the change in H_A is Δ H_A ≥ 4 since s_v was previously compensated in equilibrium by another spin indicator. We now consider any possible decrease in H_B. In the best case, all the neighboring vertices of v did not belong to the partition v was in since in this case, a change in s_v will cause 1-s_u s_v/2 to change from one to zero for all w ∈ N(v), where N(v) is the set of neighboring vertices of v. Hence, the decrease in H_B can be bounded above by Δ H_B ≤Δ, where Δ denotes the maximal degree of any vertex in G. To prevent an increase in H_A to be ever favorable in order to achieve a decrease in H_B, we must set A Δ H_A ≥ 4A > Δ B ≥ B Δ H_B, thus A > Δ/4 B. Setting B=1 leads to the choice A = Δ/4+1.



  §.§.§ A generalisation to partitioning into multiple parts


The Hamiltonian considered in Section <ref> can easily be generalised to an arbitrary number of partitions. For this, let s_vi∈{ 0,1 } be a binary (qubo) indicator which encodes if vertex v belongs to partition i. Let K be the desired number of partitions.

Define H=A H_A + B H_B + C H_C, where

    H_A    = ∑_v ∈ V( ∑_k=1^K s_vk -1 )^2,
    
    H_B    = ∑_k=1^K ( ∑_v ∈ V s_vk -|V|/K)^2,
    
    H_C    = ∑_(u,v) ∈ E∑_k=1^K (1-s_uk s_vk).

The reasoning behind the above formulation is as follows. H_A is zero if and only if each vertex v only belongs to exactly one partition, that is if for each v ∈ V there exists only one k_v such that s_v,k_v=1. Since ∑_v ∈ V s_vk is the size of partition k ∈{1,…,K}, H_B enforces partitionings of equal size |V|/K. The term H_C is the objective function: It incentivises two vertices u and v to belong to the same partition, that is s_uk=s_vk=1 (for partition k, say), and thus minimizes the edge-cut between partitions.

As in Section <ref>, the weights A and B for the two penalty terms and the weight C for the objective function have to be chosen appropriately. Consider H in a global minimum and an arbitrary s_vk. Suppose first that s_vk=1 in the optimal solution and hence changes to 0. Since we assume optimality, H_A=H_B=0, meaning that a change of s_vk from 1 to 0 results in Δ H_A ≥ 1, Δ H_B ≥ 1 and Δ H_C ≥ 0 (in H_C, a change of any s_vk from 1 to 0 potentially only results in terms 1-s_uk s_vk switching from 0 to 1). This means that changing an indicator s_vk with value 1 in the optimal solution is never favorable to achieve a further decrease in the objective function.

Next, assume that s_vk changes from 0 to 1. In this case, Δ H_A ≥ 1 and Δ H_B ≥ 1 as before. However, a change of any s_vk from 1 to 0 potentially cancels up to Δ (the maximal degree of the graph) terms of the form 1-s_uk s_vk in H_C, hence resulting in a decrease of the objective function of up to Δ H_C ≤Δ≤ |V|. To prevent that violating the penalty terms H_A and H_B ever leads to a decrease in H_C, we must hence ensure that A Δ H_A+B Δ H_B ≥ A(1+1) > C |V| ≥ C Δ H_C, meaning that A+B > C |V|. Setting C=1 leads to the condition A+B>|V|. Since the constraints lead to an under-determined equation, more than one choice for A and B exists. For simplicity, splitting up the weight equally among both penalty terms (that is, A=B) leads to the choice A=B=1/2|V| +1.



  §.§.§ Qubo model for CH-partitioning


Apart from standard edge-cut partitioning, we consider a special flavor of the graph partitioning problem termed core-halo (CH) partitioning in <cit.>. In CH-partitioning, we aim to find a partitioning of the vertices V into K sets (partitions) in such a way that the (squared) sum of both the number of vertices c_i (the core) in each partition i ∈{ 1,…,K } as well as the number of neighboring vertices in other partitions h_i (the halo) is minimized:

    min_(c_i,h_i)_i=1,…,K∑_i=1^K (c_i+h_i)^2.


In order to design a Hamiltonian for this problem consider H=A H_A + B H_B + C H_C, where

    H_A    = ∑_v ∈ V( ∑_k=1^K c_vk -1 )^2,
    
    H_B    = ∑_v ∈ V∑_i=1^K ( ∑_w ∈ N(v) ∪{v} (h_vi - c_wi - z_(v,w),i)^2 ),
    
    H_C    = ∑_i=1^K ( ∑_v ∈ V h_vi)^2

and N(v) denotes the set of all neighbors of v ∈ V. We use two types of indicators to encode whether a vertex is in the core or in the halo of a given partition. Each variable c_vi is a binary (qubo) indicator encoding with a value of one that vertex v belongs to partition (core) i. Likewise for the halo indicators h_vi. The auxiliary variables z_(v,w),i for each edge (v,w) in the graph and each partition i have no interpretation. They are used to ensure that the penalty H_B takes the value zero in equilibrium.

The reasoning behind the above formulation is as follows. As each vertex can only belong to one core, as in Section <ref>, only one of the indicators c_v1,…,c_vK can take the value one for each v ∈ V. This is enforced in H_A which is zero if and only if there is exactly one i such that c_vi=1 for each v ∈ V.

Moreover, by definition of the CH-partitioning problem, a vertex v ∈ V belong to the halo of a partition j if and only if there is an edge in the graph connecting v to a vertex w in the core of partition j. We extend this definition and define v to be a halo vertex of partition j if either it is part of core j or if there is an edge connecting v to a vertex in core j. This is enforced in H_B: Given a v ∈ V and partition number i (see the two sums in H_B), h_vi is forced to be one whenever either c_vi=1 or c_wi=1 for any w ∈ N(v), the set of all neighbors of v. Naturally, one neighbor w (of v) in partition i suffices to cause h_vi=1. Therefore, h_vi can be one even if c_wi=0 for a given w, thus causing the penalty to not attain the value zero in equilibrium. We correct this by introducing variables z_(v,w),i which in the latter case will take the value 1 and ensure that H_B=0 in equilibrium. If h_vi - c_wi=0 already holds true for given v and w, there is no incentive for z_(v,w),i to be one since this increases the overall value of the Hamiltonian.

The objective function sought to be minimized in CH-partitioning is the squared sum of core and halo nodes. Since in the above definition, h_vi=1 if v is either a core or a halo vertex, it suffices to minimize ( ∑_v ∈ V h_vi)^2, the squared sum of halos for each partition i ∈{1,…,K}, expressed in H_C.

The weights A, B and C have to be chosen appropriately again to ensure the correct global minimum of H at a valid CH-partitioning. This is done as follows. Suppose H is in equilibrium at a given configuration of c_vi, h_vi and z_(v,w),i. First, assume we change a c_vi from 0 to 1. This will not happen in practice since c_vi=1 potentially only results in more h_wi=1 (for w ∈ N(v) ∪{ v}), thus increasing the objective function H_C. Assume now that c_vi changes from 1 to 0. Since H_A=0 in equilibrium, changing one c_vi results in Δ H_A ≥ 1. Although the objective function H_C does not depend on any c_vi, changing a c_vi can still lead to a decrease of the objective. This is due to the fact that switching c_vi from 1 to 0 potentially causes Δ (the maximal degree of the graph) halo indicators to not be forced to be one any more, hence leading to a maximal decrease of H_C of Δ H_C ≤( ∑_v ∈ V h_vi)^2 - ( ∑_v ∈ V h_vi - Δ)^2 ≤ |V|^2. Together, we see that A Δ H_A ≥ A > C |V|^2 ≥ C Δ H_C, hence A > C|V|^2. Setting C=1 leads to the choice A=|V|^2+1.

Likewise, consider changing a halo indicator h_vi out of equilibrium. Assume h_vi changes from 0 to 1. In this case, we would increase the penalty H_B as well as the objective function H_C, meaning this case is infeasible in practice. Assume now that we change h_vi from 1 to 0. In this case, we incur an increase in H resulting from the penalty term H_B of at least Δ H_B ≥ 1 (since we can assume that in equilibrium, c_vi=1 whenever h_vi=1). The term H_A is independent of h_vi and hence not affected. The objective function H_C will decrease by at most Δ H_C ≤ |V|^2 - (|V|-1)^2 in the worst case if one h_vi switches from 1 to 0, hence Δ H_C ≤ 2|V|-1 ≤ 2|V|. Together, we see that Δ H_B ≥ B > 2|V|C ≥Δ H_C, thus B > 2|V|C. Using C=1 set before results in the choice B = 2|V|+1.



 §.§ Simulation results


The following section shows simulation results for an evaluation of classical edge-cut graph partitioning on D-Wave (both using the Hamiltonian for two partitions given in Section <ref> as well as the one of Section <ref>) and for CH-partitioning (Section <ref>). We compare D-Wave to  <cit.>, to simulated annealing as well as to  (see Section <ref> for details and/or the precise algorithm). All D-Wave solvers described in Section <ref> are included in our study.

We start by describing our simulation setting and proceed to evaluate the D-Wave performance on edge-cut partitioning as well as CH-partitioning. Since graph which can be embedded on the D-Wave chip are rather small, we switch to graph partitioning into two components for maximal sized graphs. The section finishes with an evaluation of the embedding vs. computing time on D-Wave.



  §.§.§ Setting and choice of test systems


All graph partitioning methods were evaluated on simulated graphs consisting of 6 to 48 vertices. For all simulated graphs, we inserted random edges between all pairs of vertices with probability 0.9. These graphs remained fixed throughout the comparison.

The largest graph size we were able to use was limited by the size of the corresponding Hamiltonian that could still be embedded on the D-Wave chip. However, solely the Sapi method for D-Wave (see Section <ref>) has this constraint – the two other D-Wave solvers (QBSolv and QSage) are (classical) hybrid methods which can handle Hamiltonians of arbitrary size.

Unless otherwise stated, all evaluations presented hereafter compare the classical methods , simulated annealing (SA) and  to the D-Wave solvers Sapi, QBSolv and QSage. QBSolv and QSage are straightforward to use since they do not allow for significant parameter tuning. For Sapi, we re-compute the graph embedding onto the D-Wave chip in each run in order to maximize the Hamiltonian size solvable on D-Wave, always perform 1000 anneals, and use the answer_mode='raw' as well as the embeddings='chains' and broken_chains='minimize_energy' post processing modes.

Computing times were measured as averages of 1000 repetitions. For D-Wave, we reported the total_real_time provided by D-Wave's solve_ising method. This time does not include the time needed to compute embeddings of the Hamiltonian onto D-Wave's chimera chip.



  §.§.§ Edge-cut partitioning




We start by evaluating D-Wave's performance on classical edge-cut graph partitioning. For added difficulty we aim to split graphs into four balanced partitions. Since we aim to split graphs into more than two partitions we use the Qubo Hamiltonian derived in Section <ref> and apply classical and D-Wave solvers to the set of graphs fixed in Section <ref>.

Figure <ref> shows the best solution found among 1000 repetitions (for the classical methods) and 1000 anneals for Sapi and QSage. QBSolv generally produces very stable results, meaning that repeating computations with QBSolv generally does not lead to varying answers.

As shown in Figure <ref>, all methods perform comparably well on the test systems, especially for small graphs. For larger graphs (towards 20 vertices),  and QBSolv find the best partitionings (in terms of a minimal edge-cut), SA performs well for small graphs but increasingly worse for larger graphs. However, the implementation we use for SA is a generic simulated annealing algorithm not further optimized (for instance, by choosing specific proposals or a certain temperature function). QSage seems to perform slightly worse than the other methods throughout the range of test graphs.

Figure <ref> shows the corresponding average computing time to Figure <ref>. As can be seen from Figure <ref>, SA is the fastest method (on average), thus providing very good results at short times for small graphs, and less good results for larger graphs.

As QSage is run with a time-out criterion (set to 60 seconds in our experiments), and as it is not possible to determine for QSage when a global minimum has been reached, using up its full computing time in each run. The time-out of 60 seconds is an arbitrary choice (although for too short timethe s, results will start getting worse).

Notably,  is able to compute very good solutions in around 0.4 seconds throughout the comparison.

The time for Sapi is constant at 0.3 seconds. This is the overall time for computation and post-processing on the D-Wave server side for 1000 anneals from D-Wave. However, Sapi is unable to embed graph Hamiltonians of the form derived in Section <ref> beyond 14 vertices.

 allows to compute partitionings which are proven to be globally optimal. However, its running time is exponential, thus greatly limiting its range of applicability. We were not able to compute partitionings with  beyond 14 vertices.

Overall, we conclude from these results that Sapi/D-Wave does not yet offer a quantum advantage over classical methods since the results obtained with QBSolv make heavy use of classical post-processing algorithms.



  §.§.§ CH-partitioning




Next, we repeat the previous comparison for CH-partitioning (again into four partitions) using the same test graphs and the Hamiltonian derived in Section <ref>. The method abbreviated by SA in this section is a simulated annealing algorithm tailored to CH-partitioning proposed in <cit.> (see Section <ref>). We employ this SA algorithm both as a post-processing step to tune partitions computed by  (which minimize the edge-cut instead of (<ref>)) and as a stand-alone algorithm initialized with random partitions.

Figure <ref> shows simulation results. Precisely, we evaluate the partitions found for each method and graph size by plotting the sum of squares criterion (<ref>) (Section <ref>), where c_i is the partition (core) size for each partition i ∈{ 1,…,K}, and h_i is the number of neighbors in adjacent partitions (the halo). As visible from the plot, performance is more similar for all methods under investigation apart from QBSolv, which seems to perform poorly.

Notably,  with a SA post-processing step performs worse than using SA on its own with a random initial partitioning allocation. QSage finds very good CH-partitions, albeit having a very long 60 second running time (Figure <ref>). As before,  finds optimal partitions in the sense that they globally minimize the sum of squares criterion, even though problems become infeasible to solve exactly as the number of vertices increases. Unfortunately, Sapi is only able to embed the CH-partitioning Hamiltonian for 6 and 8 vertex graphs. This is due to the fact that CH-partitioning uses a large number of auxiliary variables for each partition and vertex (see Section <ref>), thus resulting in relatively large Hamiltonians even for small graph sizes.

Figure <ref> shows corresponding computing times. The three classical methods ( and SA, SA by itself and QBSolv) perform comparably with a decent running time of around one second. Due to the user-specified time-out criterion of 60 seconds, QSage again has a constant running time. Sapi again achieves a fast 0.3 second anneal time but is, however, only applicable to small graph instances.  achieves the fastest runtimes for small instances among all methods considered, nevertheless its runtime quickly falls behind the other methods (apart from Sapi).



  §.§.§ Edge-cut partitioning with Ising Hamiltonian




Section <ref> already investigated edge-cut partitioning with four components. It established that D-Wave's QBSolv performs as well as classical methods, though with the help of classical post-processing to enhance the initial solution returned by the D-Wave device. Unfortunately, the tool which provides closest quantum control, Sapi, was not able to embed graph instances of sufficient size.

This is due to two factors: First, splitting a graph into more partitions uses up more qubits to indicate which partition each vertex belongs to. Thus, evaluating graph partitioning with two partitions allows for a maximal number of vertices in the input graph. Second, the Hamiltonian for graph partitioning into an arbitrary number of partitions derived in Section <ref> uses as many qubits per vertex as there are partitions, meaning that even for two partitions, two qubits per vertex are used. This unnecessarily uses up qubits and hence limits the size of the input graph. For the special case of two partitions, the Ising formulation of Section <ref> allows to allocate one qubit only per vertex.

Since for small instances, classical (brute force) methods are at least as fast as sophisticated (quantum or classical) methods, we re-evaluate graph partitioning with Ising Hamiltonian to allow for a more pronounced quantum speedup.

Figure <ref> shows simulation results for graphs up to 48 vertices (qubits in Ising formulation), the largest graph size which can be embedded on the D-Wave chip. We depict the difference in the best edge-cut found by , simulated annealing (as classical solvers) in relation to D-Wave's Sapi, since Sapi offers closest control over D-Wave. The SA algorithm we employ in this section is the one given as Algorithm <ref> in Section <ref>. The figure shows that all methods are able to find the same minimal edge-cut for instances up to 18 vertices. For larger graphs, however, D-Wave's Sapi considerably outperforms classical simulated annealing and finds partitions which even improve upon the ones found by  by up to edge-cut two.

The corresponding computing times for all three methods are displayed in Figure <ref>. As before, it shows that D-Wave/ Sapi has a constant computing time of around 0.3 seconds, while  produces good solutions in a very short time. SA runs very fast for small graphs (up to 18 vertices) but performs considerably poorer in terms of computing time for larger instances.



  §.§.§ Timings for embedding and anneal steps



In order to solve a Hamiltonian on the D-Wave chip, it is necessary to embed it first on its chimera architecture. This has to be done once for each Hamiltonian and can be quite time consuming.

Using the set-up of Section <ref> (graph partitioning with Ising Hamiltonian into two partitions using our set of test graphs), we compare the timings to compute graph embeddings to the (constant) anneal times for 1000 samples from D-Wave. Figure <ref> reports mean timings (circles) for 10 repetitions together with error bars indicating the variance. It shows that the time to compute the embeddings roughly increases linearly with the graph (Hamiltonian) size. For small instances, the embedding time is negligible. However, for large instances the embedding time quickly exceeds the actual anneal (computing) time on D-Wave, thus limiting the capabilities of D-Wave.

Alternatively, a straightforward way to mitigate this problem would be to always use pre-computed embeddings for full graphs: D-Wave allows to use pre-compute embeddings for complete graphs of size up to 45 vertices, hence allowing to use these embeddings for all graphs up to that vertex size. Nevertheless, computing individual embeddings allows to use qubits which are physically close on the D-Wave chip, thus minimizing chain lengths which increases computation accuracy.



  §.§.§ Histogram of maximal clique sizes



Figure <ref> shows a histogram of maximal clique sizes found by D-Wave and a greedy algorithm. D-Wave is used with two options, once with chain strength 0 and once with chain strength -1.






§ CONCLUSIONS


This article evaluated the performance of the D-Wave 2X quantum annealer on two NP-hard graph problems, clique finding and graph partitioning.

We compared D-Wave's solvers QSage, QBSolv and Sapi to common classical solvers (quantum adiabatic evolution of <cit.>, specialized simulated annealing of <cit.>),  of <cit.>, and  for clique finding;  <cit.>, two simulated annealing approaches of <cit.> and  for graph partitioning) to determine if current technology already allows to observe a quantum advantage in these areas. We summarize our findings as follows.

  
  * Problems need to be of a certain (problem-dependant) size to see a performance gain with D-Wave. General problems can often not be embedded on the D-Wave chip; if they can be embedded, however, D-Wave returns solutions of comparable quality than classical ones, even though classical methods are usually faster for such small instances.
  
  * For random larger instances which were designed to fit the D-Wave architecture, a speed-up for D-Wave in comparison to classical solvers is observable. However, it is not straightforward to justify the correct time to actually measure on D-Wave.
  
  * The graph embedding needed to map a particular Hamiltonian to the qubits on the D-Wave chip is computationally expensive and hence time-consuming, thus considerably distorting the time measurements if added to the actual anneal (solver) time.
  
  * In comparison to our simulated annealing (SA) implementations, D-Wave perfroms considerably better. However, it is not clear if such a comparison is fair since our implementation of SA might not be sufficiently optimized and might be outperformed by any hardware solver, either classical or quantum.

Overall, we conclude that general problems which allow to be mapped onto the D-Wave architecture are typically still too small to show a quantum advantage. Such instances can also often be solved faster (and often even exact) with best available classical algorithms (returning solutions of at least the quality that D-Wave returns). In comparison to simple SA algorithms, D-Wave is considerably faster, even though we claim no optimization for our SA routines. Selected instances especially designed to fit D-Wave's particular chimera architecture can be solved orders of magnitude faster with a quantum annealer than with classical techniques.







[Balas and Yu, 1986]Balas1986
Balas, E. and Yu, C. (1986).
Finding a maximum clique in an arbitrary graph.
SIAM J. Comput., 15:1054–1068.

[Buluc et al., 2015]Buluc2015
Buluc, A., Meyerhenke, H., Safro, I., Sanders, P., and Schulz, C. (2015).
Recent advances in graph partitioning.
arXiv:1311.3144, pages 1–37.

[Bunyk et al., 2014]Bunyk2014
Bunyk, P., Hoskinson, E., Johnson, M., Tolkacheva, E., Altomare, F., Berkley,
  A., Harris, R., Hilton, J., Lanting, T., Przybysz, A., and Whittaker, J.
  (2014).
Architectural considerations in the design of a superconducting
  quantum annealing processor.
IEEE Transactions on Applied Superconductivity, 24(4):1–10.

[Childs et al., 2000]childs2000finding
Childs, A. M., Farhi, E., Goldstone, J., and Gutmann, S. (2000).
Finding cliques by quantum adiabatic evolution.
arXiv preprint quant-ph/0012104.

[D-Wave Sys, 2016a]DWaveDevelopperGuide2016
D-Wave Sys, I. (2016a).
D-wave developer guide.

[D-Wave Sys, 2016b]dwavepostprocess24
D-Wave Sys, I. (2016b).
D-wave post-processing guide.

[D-Wave Sys, 2016c]dwave2016
D-Wave Sys, I. (2016c).
Introduction to the d-wave quantum hardware.

[Del Moral et al., 2006]DelMoral2006
Del Moral, P., Doucet, A., and Jasra, A. (2006).
Sequential monte carlo samplers.
JRSS-B, 68(3):411–436.

[Denchev et al., 2016]Denchev2016
Denchev, V. S., Boixo, S., Isakov, S. V., Ding, N., Babbush, R., Smelyanskiy,
  V., Martinis, J., and Neven, H. (2016).
What is the computational value of finite-range tunneling?
Phys. Rev. X, 6:031015.

[Djidjev et al., 2016]Djidjev2016
Djidjev, H., Hahn, G., Mniszewski, S., Negre, C., Niklasson, A., and
  Sardeshmukh, V. (2016).
Graph partitioning methods for fast parallel quantum molecular
  dynamics.
CSC 2016, 1(1):1–17.

[Fiduccia and Mattheyses, 1982]FiducciaMattheyses1982
Fiduccia, C. and Mattheyses, R. (1982).
A linear time heuristic for improving network partitions.
In Proc. 19th IEEE Design Automation Conference, pages
  175–181.

[Geman and Geman, 1984]GemanGeman1984
Geman, S. and Geman, D. (1984).
Stochastic relaxation, gibbs distributions, and the bayesian
  restoration of images.
PAMI, 6(1):721–741.

[Geng et al., 2007]geng2007simple
Geng, X., Xu, J., Xiao, J., and Pan, L. (2007).
A simple simulated annealing algorithm for the maximum clique
  problem.
Information Sciences, 177(22):5064–5071.

[Gurobi Optimization, 2015]gurobi
Gurobi Optimization, I. (2015).
Gurobi optimizer reference manual.

[Johnson et al., 2011]Johnson2011
Johnson, M., Amin, M., Gildert, S., Lanting, T. Hamze, F., Dickson, N., Harris,
  R., Berkley, A., Johansson, J., Bunyk, P., Chapple, E., Enderud, C., Hilton,
  J., Karimi, K., Ladizinsky, E., Ladizinsky, N., Oh, T., Perminov, I., Rich,
  C., Thom, M., Tolkacheva, E., Truncik, C., Uchaikin, S., Wang, J., B., W.,
  and Rose, G. (2011).
Quantum annealing with manufactured spins.
Nature, 473:194–198.

[Kadowaki and Nishimori, 1998]KadowakiNishimori1998
Kadowaki, T. and Nishimori, H. (1998).
Quantum annealing in the transverse ising model.
Phys. Rev. E, 58(5355):1–15.

[Karypis and Kumar, 1999]KarypisKumar1999
Karypis, G. and Kumar, V. (1999).
A fast and high quality multilevel scheme for partitioning irregular
  graphs.
SIAM J. Sci. Comput., 20(1):359–392.

[King et al., 2015]King2015
King, J., Yarkoni, S., Nevisi, M., Hilton, J., and McGeoch, C. (2015).
Benchmarking a quantum annealing processor with the time-to-target
  metric.
arXiv:1508.05087, pages 1–29.

[Kirkpatrick et al., 1983]Kirkpatrick1983
Kirkpatrick, S., Gelatt Jr, C., and Vecchi, M. (1983).
Optimization by simulated annealing.
Science, 200(4598):671–680.

[Lucas, 2014]Lucas2014
Lucas, A. (2014).
Ising formulations of many np problems.
Frontiers in Physics, 2(5):1–27.

[Pattabiraman et al., 2013]pattabiraman2013fast
Pattabiraman, B., Patwary, M. M. A., Gebremedhin, A. H., Liao, W.-k., and
  Choudhary, A. (2013).
Fast algorithms for the maximum clique problem on massive sparse
  graphs.
In International Workshop on Algorithms and Models for the
  Web-Graph, pages 156–169. Springer.

[Rønnow et al., 2014]Ronnow2014
Rønnow, T.F. Wang, Z., Job, J., Boixo, S., Isakov, S., Wecker, D., Martinis,
  J., Lidar, D., and Troyer, M. (2014).
Defining and detecting quantum speedup.
Science, 345:420–424.


