Quantum Annealing versus Digital Computing:An Experimental Comparison Quantum Annealing versus Digital Computing: An Experimental Comparison

The company D-Wave has introduced a special hardware for realizing the adiabatic quantum annealing approach. Their device basically finds approximate solutions of certain Ising spin glass problems on Chimera graphs. Possible interactions between qubits are defined by a Chimera graph whose nodes correspond to the qubits and whose edges give the possible interactions. The machine, “programmed” by setting the values $J_{ij}$ for the edges of the Chimera graph and $h_i$ for the qubits, finds Ising ($\pm 1$) configurations with small energy by quantum annealing.

In order to solve a problem with the help of this machine, one must write a program on a conventional front-end computer that transforms a given problem instance into a sequence of Ising problems on a Chimera graph. Since applications are rarely defined on Chimera graphs, a minor-embedding has to be employed [14], which is a challenge by itself. The Ocean Developers Toolkit (open source) has a big library of tools for problem transformation and for minor-embedding. However, it is not clear what the best strategy for a given problem is. Moreover, depending on the mapping, the experimental results of the machine may be very different. So, here, an experiment on a classical and a quantum machine would compare two different things. The D-Wave machine gives a heuristic solution of the Ising instances in a fraction of a second. It does not guarantee true optima, but conducting many runs on the same instance may increase the chance of finding them [48].

The outcome of many empirical comparisons was that the tested quantum annealers do not have enough qubits so far for solving large problems. In their very interesting report on benchmarking a smaller D-Wave machine, Parekh et al. [52] predict that this will change if the number of qubits reaches 2,048, “at which point exact classical algorithms will no longer be able to readily obtain optimal solutions for commonly used random quadratic unconstrained binary optimization (QUBO) benchmark instances on the Chimera graph.” This was a motivation for us to exploit the ability of exact solvers for instances of larger sizes in order to assess the solution quality of the current D-Wave 2000Q machine that has indeed 2,048 qubits.

Our contribution. We transform the Ising spin glass problem instances given to the D-Wave hardware into maximum cut problem instances and then use exact algorithms for solving them to optimality. We use exact approaches based on branch-and-cut [5, 12, 25, 26] and semidefinite programming [55]. All of these approaches were tailored for Chimera graphs and we compare the exact solution values to those obtained by the current D-Wave 2000Q.

We compare on the lowest possible level in the sense that we take our codes as simulators of the D-Wave hardware in order to circumvent all the messy details that come in when solving a problem by transforming it to a sequence of Ising instances for the D-Wave hardware on the classical front end. We consider this the only fair comparison, namely using universal approaches on both sides, the quantum and the classical machine side.

This allows for the assessment of the quality of solutions produced by the D-Wave hardware, rather than a compound of algorithms running on a conventional computer that may call an oracle for Ising instances on Chimera graphs. It has been a matter of discussion in the literature how well the D-Wave hardware performs at its native task, and our experiments shed some more light on this issue.

In order to determine the quality of the D-Wave solutions, we carefully select a set of 1,445 problem instances in the pursuit to cover the entire spectrum from “very easy” to “very hard,” and we invest up to 8,000 seconds of wall clock time for each instance. Our approach allows for treating both the D-Wave annealer and digital computers ranging from notebooks to massively parallel computers as “black boxes” that all address the same task, and the quality of the results as well as the wall clock times to achieve them are measured. Finally, we examine the competitiveness of a classical digital heuristic.

Since QUBO problems can be (linearly) transformed into Ising models, the D-Wave machine is capable of treating a wide class of optimization problems. So far, quantum annealing has been used for many combinatorial optimization problems such as, e.g., maximum independent set [52, 69], community detection [49, 52], clustering [8], graph isomorphism [15], and maximum 2-SAT [48]. Moreover, small real-world application problems have been solved with the help of the D-Wave quantum annealer [54, 56, 65, 66, 67].

It has been claimed that a D-Wave Two quantum annealer is 3,600 times faster than a conventional computer. This statement goes back to a computational study of McGeoch and Wang [48]. Although McGeoch and Wang [48] does not claim this in this generality, nevertheless, a race has been started in which researchers have compared the outcome of their classical approaches with those using quantum annealing by D-Wave. While McGeoch and Wang [48] have used the CPLEX quadratic programming software for solving the Ising models, Dash and Puget [24] have linearized the problem and solved the instances with the CPLEX mixed integer linear programming software. With this better approach they could observe a speed-up factor of only 17 for the McGeoch and Wang random Ising model instances. Meanwhile, there exist many comparisons in which the quantum annealing approach by D-Wave has been compared to classical algorithms, in particular on combinatorial optimization problems. Often, special purpose algorithms tailored to the considered optimization problem have been used in these comparisons with the result that the classical algorithms dominate the quantum annealing [1, 11, 38, 57, 59, 68]. The authors of these studies have pioneered the approach of avoiding minor embedding when testing the capabilities of the D-Wave machines. The article [1] is the most detailed comparison to date of D-Wave to various classical solvers, for a family of instances for which an optimal annealing time can be established on D-Wave.

Parekh et al. [52] state that “it is generally unknown which algorithm will perform best on a given instance of a combinatorial optimization problem”, and therefore they suggest to compare on the mapped instances instead. In this view, it does not seem fair to use a tailored algorithm and compare it with a universal algorithm. A conclusion drawn that the quantum annealing approach is worse would not be correct in general. Parekh et al. [52] and Coffrin et al. [18] have also suggested ideas to benchmark with quantum annealing.

Saket [58] has derived a PTAS for the Ising problem on the Chimera graph. His algorithm approximates the value of minimum energy in a Chimera graph with $n$ nodes within a factor of $(1-\epsilon)$ in time $O(n2^{32/\epsilon })$. Selby [61, 62] provides an exact algorithm based on dynamic programming that exploits the relatively low tree-width of Chimera graphs.

For our algorithms, we transform the Ising problem to a maximum cut problem as described below. The maximum cut problem has been shown to be NP-hard for general graphs by Karp [41] and APX-hard by Papadimitriou and Yannakakis [51]. Goemans and Williamson [32] have suggested a randomized constant factor approximation algorithm for instances with nonnegative weights that has been derandomized by Mahajan and Ramesh [46] and has performance guarantee 0.87856. There are a number of special cases for which the maximum cut problem can be solved in polynomial time. Among them are planar graphs [35, 44, 63], embedded one-planar graphs [21, 22], embedded graphs with a bounded number of edge crossings [16], graphs not contractible to $K_5$ [4], toroidal graphs with bounded weights [3, 31], graphs without long odd cycles [34], and instances where the edges with positive weights define a subgraph with bounded vertex cover number [47].

Early approaches based on integer linear programming models have been suggested by Barahona et al. [5, 6, 7]. A breakthrough for optimally solving Ising spin glass models was achieved in the 1990s by tailoring the algorithm to the special structure of grid graphs (see, e.g., [25, 26]). Based on this work there is a publicly available spin glass solver [64]. Whereas the cut polytope has been studied extensively for complete graphs [7], only recently, new studies on lifting and separation procedures have provided ideas for exploiting the structure of arbitrary graphs [12]. Our work is heavily based on this previous work.

Rendl et al. [55] have suggested an alternative approach based on semidefinite programming (SDP). An SDP-based branch-and-bound algorithm, the “Biq Mac Solver – Binary quadratic and Max cut Solver” [9, 55], is publicly available.

Many heuristics have been suggested for the maximum cut problem. An experimental comparison of popular approaches such as simulated annealing, tabu search, local search, and evolutionary computing-based approaches can be found in [28]. However, the mentioned study does not include approaches based on semidefinite programming and integer linear programming.

A fundamental problem in statistical solid state physics is the determination of ground states of Ising spin glasses. We can think of an Ising Model as a physical system in which the nodes, the spins, of a graph $G=(V,E)$ represent particles (atoms) and the edges correspond to the interactions between the particles. The given interactions $J_{ij}$ describe the forces acting between two spins, and the node weights $h_i$ describe the forces acting on the nodes (e.g., the effect of some external magnetic field). A classical physical state of the system describes the spin configurations for each spin, which can be either up ($+1$) or down ($-1$). Given a spin configuration $s=(s_1,s_2,\ldots ,s_n)\in \lbrace +1,-1\rbrace ^n$,

Beyond their importance in statistical physics, Ising models have many applications, e.g., in chemistry (motion of atoms), neuroscience (activity of neurons), and the layout of electronic circuits [5, 25]. The general IM has been shown to be NP-hard [2]. Choi has also shown NP-hardness for Ising models on Chimera graphs [17].

The D-Wave hardware is designed to find good approximate solutions to the Ising Spin Glass Ground State Problem. In this machine the spins correspond to the qubits. In the current version of the hardware, the possible interactions between qubits are represented by a so-called Chimera graph; i.e., a nonzero interaction $J_{ij}$ between qubits $i$ and $j$ is only possible if $ij$ is an edge in this graph. Furthermore, there is an external field that can act on each qubit $i$ (with an individual strength $h_i$).

Rather than giving a formal definition of a Chimera graph, we look at $C_4$ shown in Figure 1. The subscript $k$ in $C_k$ specifies a $k\times k$ grid of $K_{4,4}$ subgraphs. The Chimera graph $C_k$ has $8k^2$ nodes and $24k^2-8k$ edges. If there is a nonzero field, it can be represented by an extra node and extra $8k^2$ edges associated with the $h_i$-values, so that we obtain a graph with $n=8k^2+1$ nodes and $m=32k^2-8k$ edges. The D-Wave 2000Q machine uses a $C_{16}$ Chimera graph having $2,\!048$ nodes and $6,\!016$ edges. Including the external field, we have 2,049 nodes and 8,064 edges.

By quantum annealing the machine is capable of finding qubit (spin) configurations with low energy value.

Most real-world problems do not conform to the Chimera structure. Therefore, in nearly all cases it is required to embed the original problem into the hardware graph by mapping each original node to several qubits. This encompasses several problems. First, the problem, given two arbitrary graphs $G$ and $H$, to decide whether $G$ can be embedded into $H$ is NP-hard. It is assumed, but not yet proven, that this remains true when we fix $H$ to the Chimera graph. To obtain an embedding efficiently, heuristic approaches or precomputed generic embedding templates can be used. Nevertheless, for a fixed hardware the size of embeddable optimization problems is limited. In the worst case of a complete graph problem, the D-Wave 2000Q with its $C_{16}$ Chimera graph allows at most 64 fully connected original nodes. By partitioning on a classical computer and sending suitable parts to the D-Wave machine, as, e.g., suggested in [27], larger instances can be dealt with.

In turn, due to the embedding, the number of required qubits is increased. In order to implement the minor embedding, the qubits in a chain need to be coupled strongly in order to serve as a single logical variable. This can be done by assigning artificial negative weights on the edges in the chain. The absolute values of these weights must be large enough to favor chain synchronization and ensure the embedded problem is equivalent to the original one [17]. If they are set too small, the solution returned by the D-Wave machine might have broken chains. In this case postprocessing methods like majority voting could be applied to heuristically repair the solution.

If the absolute values of the weights are chosen too large, another problem of the D-Wave quantum annealer is exacerbated: its limited precision for specifying the coefficients of the Ising problem. In relation to a large weight, two coefficients might be too close to each other and thus are resolved to the same value by the machine. More precisely, the D-Wave hardware has a limited granularity of the coupling values $J_{ij}$ and field values $h_i$ in the following sense. Up to a scale factor, the available values are in the set

Any approach to a given combinatorial optimization problem using the D-Wave hardware must run a program on the conventional front end that may call the D-Wave hardware as an oracle for solving instances of (1).

Two prominent combinatorial optimization problems are especially close to the Ising problem (1): the maximum weight cut problem and the quadratic unconstrained binary optimization problem. The transformations are well known; see, e.g., [6]. We define both problems and summarize the most important transformations for the purpose of this work.

Maximum Weight Cut Problem. A cut in an undirected graph $G=(V,E)$ is defined by a subset $W\subseteq V$ and consists of the edges $\delta (W):= \lbrace ij\in E\mid i\in W,\,j\in V\setminus W\rbrace$. For edge weights the Maximum Weight Cut Problem or shortly MaxCut is to determine a cut with maximum total weight:

For the purpose of our exposition we prefer to represent edge sets $F\subseteq E$ by their characteristic vectors $\chi (F)\in \lbrace 0,1\rbrace ^E$, where $\chi _e=1$ if and only if $e\in F$. Then the problem reads

Solving (1) using a method for (2) is essential for our work and amounts to the following: Let $G=(V,E)$ be the Chimera graph plus an extra special node $v$ with connecting edges to all qubit nodes. Let $c_{ij}=J_{ij}$ for all Chimera edges and $c_{iv}=h_i$ for all extra edges. Then an optimum $x^*$ for (2) with value $c^*$ gives rise to an optimum solution $s^*$ for (1) with value $\sum _{ij\in E}c_{ij}-2c^*$ in which $s^*_i=1-2x^*_{iv}$.

Notice that the IM with zero field is trivial for nonnegative interactions, because the Chimera graph is bipartite. Hence, the value of the maximum cut is equal to the sum of all edge weights (and therefore the ground state energy is the negative sum of all edge weights).

Quadratic Binary Optimization Problem. The QUBO Problem is

Using (1) for solving an instance of (3) in which the nonzero entries $q_{ij}$ of the matrix $Q$ correspond to couplers of the Chimera graph amounts to the following: Given $Q$ and $q$, set $J_{ij}=q_{ij}$ and $h_{i}=\sum _{j=1}^{i-1}q_{ji}+\sum _{j=i+1}^nq_{ij}+2q_i$. If $s^*$ is an optimum solution of (1) with value $H(s^*)$, then $x^*_{i} = (s^*_i+1)/2$ is an optimum solution for (3) with value $(H(s^*)+\sum _{i=1}^{n-1}\sum _{j=i+1}^nq_{ij}+2\sum _{i=1}^nq_i)/4$.

Granularity Issues. With respect to granularity, an instance of MaxCut with data values in $\Gamma$ is not affected, but even for the closely related QUBO problem with all the $q_{ij}$- and $q_i$-values in $\Gamma$, the $h_i$-values resulting from the transformation above are not necessarily in $\Gamma$.

In our experiments we conservatively restrict the values using $\gamma =10$ (rather than $\gamma =30$), such that we have

We review our approaches to solve MaxCut to optimality.

Selby [61, 62] provides an exact and a heuristic algorithm for Chimera graphs based on dynamic programming. It is well known that MaxCut can be solved to optimality by dynamic programming in time $O(2^{t(G)}n)$, where $t(G)$ is the tree-width of the graph $G$ (see, e.g., [10]). Therefore, for graphs with small tree-width, dynamic programming can be an option as an exact optimization algorithm of practical use. Based on the fact that the tree-width of a Chimera graph $C_k$ is $4k$ (see, e.g., [13]), and so it is relatively small with respect to $|V|$, Selby implemented an exact algorithm, whose code is provided in [61], that exploits a tree decomposition of a Chimera graph. In [62] he reports on such an algorithm. The method is quite effective for small values of $k$; instances with $k\le 8$ can be solved regularly. However, for bigger values of $k$ the algorithm becomes impractical for both excessive time and memory requirements, which grow exponentially with $k$ in dynamic programming. For example, with randomly generated instances with $\pm 1$ weights and no node/edge failures, we estimated a time requirement of $2\cdot 10^{-6}\cdot 2^{3,47k}$ seconds and a memory requirement of $8,57\cdot 10^{-7}\cdot 2^{4,02k}$ Mbytes by a simple fitting based on the data collected from the experiments done for $k\le 8$.

We ran experiments with the algorithm on $C_8$ instances. Indeed, the times are comparable with those of the other exact algorithms. Sometimes they are the best among them. However, in the current implementation [61], only the ground state energy is computed, but not the ground state.

Selby [62] presents a heuristic for MaxCut on Chimera graphs that exploits the ability of solving MaxCut on low tree-width graphs to optimality but that can be used also for large Chimera graphs. The algorithm is of a randomized type and is based on the following technique, called subgraph sampling, inspired by earlier work [30, 37]. Here we briefly outline the algorithm.

We are given an instance of a weighted graph $G=(V,E)$. We define an ordered collection of subsets of $V$, i.e., ${\mathcal {S}}=\lbrace S_i\subset V\mid i=1,\ldots , t\rbrace ,$ with the following properties:

• Each $S_i$ is a large subset of $V$.
  
• The induced subgraph $G[S_i]$ is connected and has a tree-width bounded by a parameter $r$ that does not depend on $|V|$.
  
• $G$ is the union of the induced subgraphs $G[S_i]$, for $i=1,\ldots , t$.
  

The randomly ordered collection is built by setting the value of a parameter $w$ that controls the tree-width of its elements and by generating three random numbers: $x$ in the set $\lbrace 0,1\rbrace$, $y$ in the set $\lbrace 0,1,\ldots ,w+1\rbrace$, and $z$ in the set $\lbrace 0,1,\ldots ,k-1\rbrace$. Then the collection is given by

The inner algorithm. Given a weighted graph $G=(V,E)$, a node subset $S$ with tree-width at most $r$, and a spin configuration of the nodes in $V\!\setminus \! S$, the algorithm finds the spin configuration of the nodes in $S$ that yields the maximum cut among all spin configurations of the nodes in $S$.

The configuration of the nodes in $V\!\setminus \! S$ is never changed by the algorithm. These nodes act as an external magnetic field on the nodes in $S$. If the parameter $r$ is small, the algorithm is pretty fast in finding the optimum solution even if $S$ has a large size. Its basic task (executed several times) is to choose, among the 16 spin configurations related to some four-node sets (the h-nodes and the v-nodes), the one that yields the best objective function value. To avoid making non-deterministic decisions in case of a tie, at the beginning of the inner algorithm a suitably small breakeven-avoiding random perturbation is generated that is added to the value of each of the 16 configurations.

The outer algorithm. Given a weighted graph $G=(V,E)$, a node subset $S$ with tree-width at most $r$, and a spin configuration of the nodes in $V$, the algorithm generates a randomly ordered collection $\lbrace S_i\subset V\mid i=1,\ldots , t\rbrace$ of $t$ subsets of $V$. Then, for every $i\in 1,\ldots ,t$, it executes the inner algorithm for $S_i$ using the current spin configuration to set one of the nodes in $V\!\setminus \! S_i$. The solution generated by the inner algorithm becomes the current spin configuration used for the next run.

The heuristic. The heuristic starts by generating a random spin configuration of the nodes in $V$. Then the outer algorithm is executed repeatedly. Between two consecutive runs either a new spin configuration is randomly generated afresh or the best configuration obtained in the previous run is partially randomly perturbed. The latter is the default option. In the default setting the random perturbation is obtained by randomly generating a configuration for the nodes of 20% of the $K_{4,4}$ subgraphs, randomly chosen.

The algorithm is run until a certain condition is satisfied. For example, in our experiments, the algorithm was stopped when the elapsed wall clock time exceeded a prescribed duration.

In view of the fact that we have classical algorithms that are able to compute optimum solutions to the native problem that the D-Wave hardware solves heuristically, we have selected a benchmark set of problem instances and implemented increasingly sophisticated classical computer programs that give optimum solutions to the same instances that are given to the D-Wave hardware. Of course, we cannot expect to solve all instances to optimality in reasonable time, but when we fail, our classical approaches deliver a solution along with a quality guarantee.

The goal of our experimental study is a contribution to answering the following questions:

• How close are the D-Wave solutions to optimum or best known solutions?
  
• How well does D-Wave perform in terms of computation time in comparison to our classical exact approaches?
  
• Can classical heuristics compete with the D-Wave hardware in terms of quality and solution time?
  

The “classical” digital computers we use in our experiments range from notebooks to a computer cluster with 720 Intel E5-2690v2 3 GHz CPU cores and 4.5 TB of main memory; the quantum computer is a D-Wave 2000Q. We use all these devices as “black boxes” that all address the same task, and the quality of the results as well as the wall clock times to achieve them are measured.

In the following, we introduce our benchmark set of 1,445 instances. Then we explain the D-Wave parameter settings as well as the hardware and software for our exact approaches and the Selby heuristic.

The easiest instance sets for our exact approaches have been pga, mis, and k64ising, all of which could be solved within a few seconds. The reason for the pga instances is their small size. On the other hand, the mis instances belong to the largest tested instances. Before transformation to Ising instances on Chimera graphs, they are independent set instances on bipartite graphs that can be solved in polynomial time. Our approaches seem to be able to profit from this fact, while the D-Wave machine does not. For D-Wave, these instances belong to the hardest tested instances. This confirms observations made by Parekh et al. [52]. The k64ising instances have a nonzero field and belong to the sparsest instances in our study, which explains that general_maxcut solves them very quickly. Also for D-Wave, these instances are not hard, although the ground state cannot be found for 40% of the instances.

The instance sets of medium difficulty for our exact approaches have been rfr, mgw-c8-439, and mgw-c8-507. These were solved to optimality in a few minutes. Although the rfr instances belong to the largest tested instances, their weight distribution (weights of wider range) made branch-and-cut approaches easier. This is in line with observations that have been made in the context of solving Ising spin glass problems on grid graphs, where instances with a Gaussian distribution on the edge weights have turned out to be much easier than those with $\pm 1$ edge weights. The McGeoch-Wang instance sets mgw-c8-439 and mgw-c8-507 belong to the smallest with a nonzero field. They are also easy for D-Wave.

The most difficult instance sets for our exact approaches have been selby-c8, k64maxcut, lga, mgw-c16-2031, and selby-c16. The Selby instances have been constructed to be very hard (zero field and difficult weight distribution) and they have kept their promise. Nevertheless, we have been able to solve all selby-c8 instances to optimality due to their small size, but the selby-c16 have been the hardest instances of the entire study, both for our exact algorithms and for the D-Wave machine that had the largest average gap for this set. With an edge/vertex ratio of only 1.26, k64maxcut is the sparsest instance set, which might explain why general_maxcut has been able to solve some of them to optimality within 5 seconds, yet other than the sparse k64ising instances they have a zero field, making them harder as already observed in the study of Dash and Puget [24]. The lga instances have been constructed to be large versions of the pga instances, and they are much larger. Their edge/vertex ratio is roughly twice the one of the k64 instances, and their weight structure makes them so hard that we have not been able to solve them all to provable optimality. The mgw-c16-2031 instances are the largest and densest of the entire study, but nevertheless, we have been able to determine ground states for all of them.

For $1,\!423$ of the $1,\!445$ instances of this study, we have been able to compute proven ground states within the chosen wall clock time limits. For one more instance (selby-c8-06) we have certified optimality of our best solution with the help of chimera_maxcut without time limit. We have tried and failed so far to do the same for the remaining 21 cases. But for the latter, our best solutions have a proven high quality: for the 20 selby-c16 instances, the maximum relative deviations from the ground state energy range from 1.208% to 2.171%, and for the only remaining unsolved lga-09 instance, the maximum relative deviation is 0.354%.

Table 8 provides a summary of our findings. #DW-best gives the number of times D-Wave finds a true ground state or a spin configuration of the best known energy, and %DW-gap gives the maximum and average deviations from the best known energy in percent.

Before we draw conclusions, it is appropriate to recall and discuss our experimental setup. Our computational study is based on the “black box approach,” in which all classical methods and the D-Wave method are treated as black boxes whose results are measured as a function of the wall clock time provided to them. Our D-Wave parameter choice results in a fixed D-Wave black box wall clock time of 1.6 seconds.

Of course, we might have tested various variations within each black box. In particular, in the D-Wave black box we might have varied the parameters $T_{\text{anneal}}$ and $R$, keeping the 1.6-second total time or trying a shorter or longer total time. Even when keeping the 1.6-second total time, we might have adapted the internal parameter settings for particular instance sets in the quest of obtaining better performance. The limited resources did not allow for such variations. We would like to make clear that our discussion and conclusions are based on our fixed black box parameter choices that are in line with recommendations in the literature.

Based on the results presented in Section 9, we make an attempt at answering question (Q1):

• The energies of D-Wave solutions are always close to ground state energies; the maximum deviation we have observed is 0.647%. On the other hand, D-Wave misses true ground states very often.
  

Yet we have no short answer to question Q2:

• If we restrict to instance sets for which D-Wave reliably provides ground states, the classical methods do so faster (pga) or at roughly the same speed (mgw-c8-439). For the set mgw-c8-507, D-Wave finds all ground states at a faster speed. On the other hand, for the mis set, classical computation (on a notebook) outperforms D-Wave both in time and in quality. If we consider the 410 instances other than pga, then chimera_maxcut has determined proven ground states for 111 of these in less than 1.6 seconds. For these “easiest” instances beyond the pga instances D-Wave misses a true ground state for six instances, i.e., more than 5%. For all the other instance sets, digital computation takes longer or very much longer than 1.6 seconds yet produces better solutions in the majority of cases, as well as proofs of optimality or quality guarantees.
  

Concerning question (Q3), the only example heuristic we tried gives a clear picture: in less than 4.3 seconds, the Selby heuristic finds the best known solutions to all instances, while D-Wave fails to do so very often for the harder instance sets.

• There is at least one heuristic (Selby) that performs better than D-Wavefor all but 4 of the 1,445 instances both in terms of quality and time, and for the remaining 4 in terms of quality.
  

By quantum annealing, the D-Wave computer realizes a fast heuristic for certain Ising problems. It cannot be expected to deliver true optimum solutions, but the relative deviations will usually be very small: rhe maximum relative deviation from the ground state energy is 0.647% over all 1,424 instances in this study for which we have been able to compute proven ground states. For the remaining 21 instances the deviation could be worse.

We observed that for $C_8$/$C_{16}$ pairs like mgw-c8/mgw-c16 and selby-c8/selby-c16 as well as for the pair pga/lga, the quality of the D-Wave solutions deteriorates with increasing problem size.

The exact algorithms provide the necessary certificates with which we can evaluate D-Wave and the Selby heuristic. However, we should stress the fact that exact optimization requires a lot of time to prove optimality, and thus it is not fair to compare their times with the heuristic times, but even with this additional burden, the exact algorithms are faster than D-Wave on a large portion of the sample.

If one needs to compute proven optimum solutions, then general MaxCut solvers based on polyhedral combinatorics or semidefinite programming are the methods of choice.

When we started this study, we conjectured, “What is hard for our mathematical optimization solvers in terms of running time is hard for the quantum annealer in terms of solution quality.” But the outcome of our experiments shows that this is not always true.

The most surprising result is the remarkably excellent performance of the Selby heuristic. Like D-Wave, the Selby heuristic cannot give an approximation quality guarantee, but the experiments show that it can produce very good solutions in short computation time. We have demonstrated that although D-Wave obtains near-optimal solutions to most of the instances with 1.6 seconds, it misses out on the optimal solutions for many instances (224), and the Selby heuristic (running on 520 cores) is qualitatively better. In 1.6 seconds, it obtains the best-known solution (optimum or otherwise) for all but 4 instances (i.e., for 1,441 of the 1,445 instances), and for the 4 instances where it does not obtain the best-known solution, its solution quality is much better than the one obtained by D-Wave. In 4.6 seconds it obtains the best-known solution for all 1,445 instances.

In this sense, the Selby heuristic is always the winner in terms of quality and also in terms of time for all instances except four, which all belong to the selby-c16 set, for which D-Wave energy values are always inferior to those obtained by the Selby heuristic.

While our exact approaches as well as the Selby heuristic have been tailored to the special architecture of Chimera graphs, we expect that similar exact approaches will also be possible for future architectures of D-Wave machines, and also new “classical” heuristics adapted to future architectures might well emerge.

It seems not unfair to compare a classical network of 28 digital processors with 20 threads each with a quantum annealer that has 2,048 qubits. If the former outperforms the latter, it is questionable if a quantum annealer can outperform classical computers for solving combinatorial optimization problems.

Since our instances are publicly available, our study meets the important basic requirement that our results are falsifiable by using a future version of the D-Wave computer or any other quantum annealer. In particular, the emerging new D-Wave Pegasus architecture will easily accommodate our instances because the $C_{16}$ Chimera graph is a subgraph of the Pegasus graph.

It may well be (and we hope) that the exciting new quantum computer technology will make leaps in the future, but in our experiments we have certainly not observed superior performance of quantum annealing in comparison to “classical” methods. We hope that the scientific community, especially in the realm of quantum annealing, will take up the challenge to experiment with our instance set with known optimum solutions and find improved approaches to quantum annealing, including better parameter settings for current systems.

The authors would like to thank NASA Quantum Artificial Intelligence Laboratory for many valuable discussions and the opportunity to use the D-Wave 2000Q machine at NASA Ames. Work of the third author is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC-2047/1 – 390685813. Finally, this article has profited a lot from the advice of two anonymous referees.