Quantum annealing with more than one hundred qubits

Methods

Quantum annealing was performed on the D-Wave One Rainer chip installed at the Information Sciences Institute of the University of Southern California. The device has been described in extensive detail elsewhere harris_experimental_2010_1 ; Harris2010 ; berkley2010scalable . After programming the couplings, the device was cooled for $2.5\,s$, and then $1000$ annealing runs were performed using an annealing time of $t_{f}=5\,\mu s$. Annealing is performed at a temperature of $0.4\,\textrm{GHz}$, with an initial transverse field starting at $A(0)\approx 5\,\textrm{GHz}$, going to zero during the annealing, while the couplings and local fields are ramped up from near zero to about $B(t_{f})\approx 5\,\textrm{GHz}$ at the end of the schedule. Details of the schedule and results for longer annealing times are provided in the supplementary material.

Simulated annealing was performed using the Metropolis algorithm with local spin flips with codes optimised for the $\pm 1$ couplings used as test problems. A total of $N_{\rm updates}$ flips per spin were attempted, increasing the inverse temperature $\beta=1/k_{B}T$ linearly in time from $0.1$ to $3$. The simulated quantum annealing simulations were performed in both discrete and continuous time path integral quantum Monte Carlo simulations with cluster updates along the imaginary time direction to account for the transverse field, combined with Metropolis rejection sampling for the Ising interactions (see the supplementary material for details).

The classical spin dynamics model replaces the quantum spins by $O(3)$ classical unit vectors $\vec{M}_{i}$, where the sign of the $z$-component of each spin is the value of the Ising variable. The spins are propagated via the equations of motion $\frac{\partial\vec{M}_{i}}{\partial t}=\vec{H}_{i}(t)\times\vec{M}_{i}$, where the time-dependent field $\vec{H}_{i}(t)$ acting on spin $i$ is a sum of a decaying transverse field (along the $x$ direction) and a growing coupling term (along $z$): $\vec{H}_{i}(t)\equiv(1-t/t_{f})h\hat{e}_{x}-(t/t_{f})\sum_{j}J_{ij}M^{z}_{i}\hat{e}_{z}$. The initial condition is to have all spins perturbed slightly from alignment along the $x$ direction: $(-\sqrt{1-\delta_{i}^{2}-\eta_{i}^{2}},\delta_{i},\eta_{i})$, where $|\delta_{i}|$, $|\eta_{i}|<0.1$.

Gauge averaging was performed on the device by using gauge symmetries to obtain a new model with the same spectrum. This was achieved by picking a gauge factor $a_{i}=\pm 1$ for each qubit, and transforming the couplings as $J_{ij}\rightarrow a_{i}a_{j}J_{ij}$ and $h_{i}\rightarrow a_{i}h_{i}$. Success probabilities $s_{g}$ obtained from $G$ gauge choices were arithmetically averaged for the correlation plots and as $\overline{s}=\prod_{g=1}^{G}(1-s_{g})^{1/G}$ for the scaling of total effort (see supplementary material for a derivation).

The ground state energies were obtained using exact optimisation algorithms, an exact version of belief propagation using Bucket Sort dechter1999bucket and a related optimised divide-and-conquer algorithm described in the supplementary material.

Acknowledgements We acknowledge useful discussions with M.H. Amin, M.H. Freedman, H.G. Katzgraber, C. Marcus, B. Smith and K. Svore. We thank Lei Wang for providing data of spin dynamics simulations, G. Wagenbreth for help optimizing the belief propagation code and P. Messmer for help with optimizing the GPU codes. Simulations were performed on the Brutus cluster at ETH Zurich and on computing resources of Microsoft Research with the help of J. Jernigan. This work was supported by the Swiss National Science Foundation through the NCCR QSIT, by ARO grant number W911NF-12-1-0523, by ARO MURI grant number W911NF-11-1-0268, by the Lockheed Martin Corporation, by DARPA grant number FA8750-13-2-0035, and by NSF grant number CHE-1037992. MT acknowledges hospitality of the Aspen Center for Physics, supported by NSF grant PHY-1066293. The initial planning of the experiments by MT was funded by Microsoft Research.

Author Contributions MT, JM and DL designed the experiments and wrote the manuscript, with input from all other authors. SB and ZW performed the experiments on D-Wave One. SI, TR and MT wrote the simulated classical and quantum annealing codes and TR, SI, MT and DW performed the simulations. SB and TR wrote the bucket sort code and divide and conquer codes. TR, SI, MT, SB, ZW and DL evaluated the data. All authors contributed to the discussion and the presentation of the results.

Supplementary material for “Quantum annealing with more than one hundred qubits”

I Overview

Here we provide additional details in support of the main text. Section II shows details of the chimera graph used in our study and the choice of graphs for our simulations. Section III expands upon the algorithms employed in our study. Section IV presents additional success probability histograms for different numbers of qubits and for instances with magnetic fields, explains the origin of easy and hard instances, and explains how the final state can be improved via a simple error reduction scheme. Section V presents further correlation plots and provide more details on gauge averaging. Section VI gives details on how we determined the scaling plots and how quantum speedup can be detected on future devices. Finally, section VII explains how the spectral gaps were calculated by quantum Monte Carlo (QMC) simulations.

II The chimera graph of the D-Wave device.

The qubits and couplers in the D-Wave device can be thought of as the vertices and edges, respectively, of a bipartite graph, called the “chimera graph”, as shown in figure 1. This graph is built from unit cells containing eight qubits each. Within each unit cell the qubits and couplers realise a complete bipartite graph $K_{4,4}$ where each of the four qubits on the left is coupled to all of the four on the right and vice versa. Each qubit on the left is furthermore coupled to the corresponding qubit in the unit cell above and below, while each of the ones on the right is horizontally coupled to the corresponding qubits in the unit cells to the left and right (with appropriate modifications for the boundary qubits). Of the 128 qubits in the device, the 108 working qubits used in the experiments are shown in green, and the couplers between them are marked as black lines.

For our scaling analysis we follow the standard procedure for scaling of finite dimensional models by considering the chimera graph as an $L\times L$ square lattice with an eight-site unit cell and open boundary conditions. The sizes we typically used in our numerical simulations are $L=1,\ldots,8$ corresponding to $N=8L^{2}=8,32,72,128,200,288,392$ or $512$ spins. For the simulated annealers and exact solvers on sizes of $128$ and above we used a perfect chimera graph. For sizes below $128$ where we compare to the device we use the working qubits within selections of $L\times L$ eight-site unit cells from the graph shown in figure 1.

In references Choi1 ; Choi2 it was shown that an optimisation problem on a complete graph with $\sqrt{N}$ vertices can be mapped to an equivalent problem on a chimera graph with $N$ vertices through minor-embedding. The tree width of $\sqrt{N}$ mentioned in the main text arises from this mapping. See Section VI.1 for additional details about the tree width and tree decomposition of a graph.

III Classical algorithms

Simulated annealing (SA) is performed by using the Metropolis algorithm to sequentially update one spin after the other. One pass through all spins is called one sweep, and the number of sweeps is our measure of the annealing time for SA. Our highly optimised simulated annealing code, based on a variant of the algorithm in Ref. J.Stat.Phys.44.985 ; Comput.Phys.Commun.59.387 , uses multi-spin coding to simultaneously perform 64 simulations in parallel on a single CPU core: each bit of a 64-bit integer represents the state of a spin in one of the 64 simulations and all 64 spins are updated at once. A similar code for GPUs uses 32-bit integers and additionally performs many independent annealing runs and updates many spins in parallel in multiple threads.

The performance of our codes on the classical reference hardware is shown in Table 1. We use high-end chips at the time of writing, an 8-core Intel Xeon E5-2670 “Sandy Bridge” CPU and an Nvidia Tesla K20X “Kepler” GPU. To find a ground state of our hardest 108-spin instances with a probability of 99%, this translates to a median annealing $32\mu$s on a single core of the CPU, $4\mu$s on eight cores, and 0.8$\mu$s on the GPU, which should be compared to $15\mu$s pure annealing time on the D-Wave device for the same problems.

III.2 Simulated quantum annealing

For simulated quantum annealing we use both a continuous time algorithm and a discrete time algorithm.

The continuous time algorithm Eur.Phys.J.B.9.233 constructs segments of a world line in the (imaginary) time direction and flips them using the Metropolis algorithm. Specifically, we pick a random site and introduce new cuts in the time direction via a Poisson process. Then we calculate the overlaps with the neighbouring sites and use these overlaps to calculate the Metropolis acceptance ratios $P_{\text{Met}}$ and flip a segment with probability $P_{\text{Met}}/2$. We cannot grow the cluster along the space directions as in Ref. Eur.Phys.J.B.9.233, , which connects segments into larger clusters, since such cluster updates are inefficient in frustrated models like our spin glass.

In order to implement the fastest possible simulated quantum annealing code we also implemented a discrete time algorithm similar to that outlined in Ref. PhysRevB.66.094203, . However, unlike Ref. PhysRevB.66.094203, . we again used cluster updates along the imaginary time direction, typically with 64 time slices.

To verify that the discrete time code produces results similar to the continuous time algorithm – i.e., that the error in the imaginary time direction is small – we show a correlation plot between the discrete- and continuous time in figure 2A). The scatter observed in this plot is a measure for the dependence of success probabilities on details of the simulated quantum annealing implementation. We also show correlations for a continuous time Monte Carlo using two different sets of random seeds for initial configurations and updates in figure 2B). The scatter of points is within the 3% ($1/\sqrt{1000}$) error expected for the success probabilities when performing 1000 annealing runs per instance.

In Table 2, we summarise the performance of these two codes for typical cases using the linear schedule of figure 3A).

We have performed simulated quantum annealing with two different schedules shown in figure 3: A) a linear schedule (referred to as schedule I) where the transverse field is ramped down linearly in time and the Ising couplings are ramped up linearly, and B) the schedule used in the D-Wave device (referred to as schedule II). The performance was similar in both cases. For the scaling plots we used the linear schedule I at an optimised temperature ranging from $T=0.33$ to $T=1$ depending on system size, which gives slightly better performance than schedule II.

For the correlation plots we use the continuous time code (CTQ) with schedule II – the average over slightly different schedules for the individual qubits on the device – but at up to ten times lower temperature than the device temperature of 20mK (0.4 GHz). The simulated quantum annealer requires a temperature about three times lower than the nominal temperature of the device to exhibit a clear bimodal distribution. This can be explained and motivated as follows. When the transverse field is strong the quantum Monte Carlo updates mimic the quantum tunneling taking place in the device, however when the transverse field becomes smaller these Monte Carlo updates turn into local spin flips of a classical (thermal) annealer. The device in this regime, on the other hand, has high tunneling barriers between the two states of the qubits of the device that suppress thermal tunneling over the barrier. To achieve a similar suppression of thermal excitations we need to lower the temperature in the quantum Monte Carlo simulations by at least a factor two. Lowering by more than a factor of ten does not significantly change histograms or correlations, indicating that at the chosen temperature the simulated quantum annealer is dominated by quantum tunneling and not thermal effects.

III.3 Exact solvers

We investigated four exact solvers, akmaxsat akmaxsat , the biqmac algorithm biqmac used in the spin glass server sgserver , exact belief propagation using bucket sort dechter1999bucket and a related divide-and-conquer algorithm. The latter is specifically designed for the chimera graph, generalising divide-and-conquer for the square lattice. We consider a chimera graph of $M\times M$ 8-spin unit cells ($8M^{2}$ spins). For each possible configuration of the $4M$ spins on the left side of the unit cells in the first row (which couple vertically) we find and store the optimal configuration of the remaining $4M$ spins with effort $4M$, giving a total effort of $4M\,2^{4M}$. Finding the $2^{4M}$ optimal configurations of the next row of unit cells, one builds up the solution row by row, scaling as ${\cal O}\left(M^{2}\,2^{4M}\right)$. Since $M\propto\sqrt{N}$, where $N$ is the number of spins, this scaling demonstrates explicitly the claim made in the main text that exact solutions scale no worse than ${\cal O}(\exp(c\sqrt{N}))$, as for the bucket sort algorithm. We present the scaling of these algorithms in section VI.1.

IV Success probability histograms

In this section we provide additional experimental and simulation data, complementing the data shown in the main text.

In addition to the histogram for $108$ qubits shown in the main text, we also show histograms for $58,84,87$ and $108$ qubits without local fields in figure 4 and with local fields in figure 5. By comparing these two figures, one notes that the cases with local fields are in general easier for the D-Wave device. We also verified that this is the case for the simulated annealers - see figure 6.

To check that the bimodality is not due to faulty couplers or qubits we performed the following analysis for the hardest instances of $N=108$ spin problems where the D-Wave device did not find the ground state with one gauge choice. For each of these instances we looked at the lowest energy configurations reached and determined which spins differ compared to the closest ground state configuration. Closeness is measured in terms of the Hamming distance, which is the number of spins that have to be flipped to reach a ground state configuration. We did not observe a strongly peaked distribution which would have indicated a singly faulty qubit or coupler.

IV.2 Simulated annealing and simulated quantum annealing

Figures 7, 8 and 9 show the success probability histograms for $N=108$ spin problems for simulated annealing, simulated quantum annealing with different number of annealing sweeps (updates per spin), and simulated quantum annealing at different temperatures. While in the simulated classical annealer the distribution is always unimodal and shifts towards larger success probabilities upon increasing the annealing time, the simulated quantum annealer becomes more strongly bimodal when increasing annealing times.

Figures 10 and 11 show the success probability histograms for simulated quantum annealing for instances without and with local fields when increasing the problem size. As the problem size is increased the weight in the peak at low success probability increases. While for the instances without local fields the bimodality of easy and hard instances vanishes at about $N=288$ spins, it remains up to about $512$ spins for instances with local fields. This indicates that there are more easy instances here than in the case without fields, and hence that instances with fields can be used to test for quantumness in a device scaled up to more qubits. However, once the “easy” problems disappear beyond $512$ spins, and with them the bimodality, alternative analysis methods will be necessary to test for quantumness.

Supplementing the last two figures, we show the evolution of the success probability histograms for different numbers of spins for a simulated classical annealer in figure 12. It is seen that, in contrast to simulated quantum annealing where the weight in the two peaks of the bimodal distribution changes, here at fixed annealing time the unimodal peak gradually shifts to smaller success probabilities.

IV.3 Hardness, gaps and free qubits

Figures 3 and 4 of the main text show that the “hard” instances typically exhibit small gaps during the evolution and often get trapped in excited states with a large Hamming distance to the closest ground state. The correlation between small gaps and large Hamming distance can be understood with a simple perturbative argument in the regime of small transverse fields $-\Gamma\sum\sigma^{x}$, where most of the small gap avoided level crossings appear. At an avoided level crossing between two states with Hamming distance $d$, $d$ spins need to be flipped to adiabatically follow the ground state. In perturbation theory, the tunneling matrix element between the two states is of order $\Gamma^{d}$, exponentially small in the Hamming distance. The small matrix element not only poses problems for adiabatic evolution but also suppresses quantum Monte Carlo updates that connect the two states. This common origin of the hardness in both quantum annealing and simulated quantum annealing explains the observed correlations despite the different underlying dynamics (deterministic vs stochastic).

An intuitive understanding of how such small gap avoided level crossings can arise can be obtained by considering degenerate states that are connected by single spin flips. The free qubits of a state are the qubits that can be flipped without changing the energy. From perturbation theory, a small transverse field $\Gamma$ breaks the degeneracy of each free qubit PhysRevB.63.224401 ; amin2009 ; Boixo2012 . In the simplest case, degenerate states form a hypercube and have the same free qubits. If the unperturbed energy was $E_{0}$, the lowest energy state of a hypercube with $F$ free qubits will then have energy $E_{0}-\Gamma F$ to first order in perturbation theory. If the low energy excited states have more free qubits than the ground states, their energy will be lowered more, resulting in avoided level crossings and small gaps. As seen in figure 13, hard instances tend to have more free qubits in low energy excited states. This phenomenon has been previously observed in the random subcubes model, a specially constructed toy model Bapst2013 . It can also be seen that on the device the experimental average of free qubits for first excited states is higher than the unweighted average over all excited states. This is not the case in simulated annealing.

In the spin glasses with $\pm 1$ couplings chosen here, we can expect to find a significant number of free qubits per state. This is because many qubits have six couplings, that can cancel each other. We also expect more low energy excited states than ground states, and consequently some low energy excited states will have more free qubits than most ground states. As argued above, both physical and simulated quantum annealing are sensitive to this problem: many spins updates are needed to move between the states at both sides of the avoided crossing. On the other hand, classical simulated annealers, not having a transverse field, do not suffer from these avoided crossings. This might explain why simulated quantum annealing scales slightly worse than classical annealing for the $\pm 1$ spin glasses considered in this work.

IV.4 Improving the final state by fixing single spin errors

Single spin flip errors can be fixed with a linearly scaling effort by checking once for each spin whether the total energy can be lowered by flipping it, and then flipping the spin that lowers the energy the most. To illustrate the effect of this procedure we show in figure 14A), the success probability histogram for $N=108$ spin instances with local fields, with and without error reduction. Without error reduction only ground states (Hamming distance $d=0$) count as a successful annealing run, while with error reduction also the runs giving states a Hamming distance $d=1$ away from the ground state will, after error reduction, end up in the ground state. It is clear from figure 14A) that while “hard” instances do not benefit from this error reduction scheme, such single spin flip errors are a dominant error source for the “easy” instances and their failure rate is substantially reduced using error reduction. The most likely source for such errors are thermal excitations or readout errors.

To understand the dependence on problem size $N$ and hardness, we consider the percentiles of the success probability distribution as a function of $N$ in figure 14B) to E). We see that the success probability decreases upon increasing $N$, i.e., the probability of such errors increases for larger problems. Independently of $N$, single-bit error reduction is ineffective in increasing the success probability for the “hard” instances (high percentiles). It becomes more and more effective as problems become easier (low percentiles), improving the success probability to very nearly 1, independently of $N$, at the $1\%$ percentile.

Variants of this simple error reduction scheme can be derived that, also with linear scaling, fix not only one single spin errors, but also “disconnected” single spin errors (i.e., flips of spins that are not connected by one of the couplings). One such example is zero temperature Glauber dynamics, flipping randomly selected spins as long as they lower the energy. However, we saw only minimal improvements using such alternative schemes and sometimes worse results when the “wrong” spins are flipped and one reaches a different local minimum. The reason may be that disconnected single spin errors are not common in the problem sizes we studied, but that may change in larger problems.

V Correlations

To better understand the correlations we here show copulae in addition to the correlations shown in the main text. Copulae

factor out the marginal distributions $f_{1}(x_{1})$ and $f_{2}(x_{2})$ (e.g., the strong bimodal distribution) from the joint probability function $f(x_{1},x_{2})$ describing the correlation density. For independent sets the copula density is $c(x_{1},x_{2})=1$ while for perfect correlations it is a delta function $c(x_{1},x_{2})=\delta(x_{1}-x_{2})/f_{1}(x_{1})$.

To plot copula densities of success probabilities we replace each success probability by its rank after sorting the success probabilities divided by the number of values, and then plot the correlation densities of these normalised ranks. In figure 15 we plot copulae corresponding to the correlations shown in the main text, enhancing the visibility of the correlations especially in the corners.

In figure 16 we show the copulae between the hardness for akmaxsat and the D-Wave device and the simulated classical and quantum annealers. We do not observe any correlations.

V.2 Reproducibility

To verify reproducibility of the QA data we performed experiments on $N=108$ spin instances three times, with a month between the first and the second two repetitions. The correlations (see figure 17) show that the device is stable over the time of a month.

Strong deviations are seen for a small fraction of the instances. These are most likely due to 1/f noise or “programming errors” when flux quanta are loaded into the programmable magnetic memories to program a specific set of couplings into the device. Since these programming errors will limit the correlations between any model and the device, no better correlations than shown in this figure can be expected between our simulations and the device.

V.3 Gauge averaging

The spectrum of an Ising spin glass is invariant under a gauge transformation that changes spins $\sigma^{z}_{i}\rightarrow a_{i}\sigma^{z}_{i}$, with $a_{i}=\pm 1$, when at the same time changing the couplings $J_{ij}\rightarrow a_{i}a_{j}J_{ij}$ and the local fields as $h_{i}\rightarrow a_{i}h_{i}$. While the simulated annealers are invariant under such a gauge transformation, the calibration of the D-Wave device is not perfect and breaks the gauge symmetry. Different gauges hence realise slightly different physical systems with different success probabilities. We average over gauges to reduce calibration uncertainties.

Panel A of figure 18 shows a scatter plot of the hardness of instances in the simulated quantum annealer (SQA) and the D-Wave device (DW). The high density in the lower left corner (hard for both methods) and the upper right corner (easy for both methods) confirms the similarities between the quantum device and a simulated quantum annealer. However, a small percentage of instances appears in the other corners (hard for one method but easy for the other). We attribute these to calibration errors of the couplers in the device which vary by about 10% Harris2010 . To test for calibration errors we performed a gauge transformation on each instance to realise a different encoding of the same spin glass instance on the device. The correlations between two different encodings (gauge choices), shown in panel B, turn out to be comparable to those between the SQA and the device, demonstrating that the observed deviations can be attributed to calibration errors.

To minimise calibration errors, we then performed annealing with multiple encodings of the same instance related by gauge transformations and averaged the success probabilities. The resulting correlations between the simulated quantum annealer and the gauge-averaged results from the D-Wave device, shown in panel C (identical to the plot shown in the main text), are significantly improved compared to a single gauge choice, with a substantial reduction of the number of extreme outliers. These correlations are even better than those between a single embedding and the gauge averaged results on the device (see panel D).

In figure 19 we show correlations between the arithmetic average over eight gauges compared against the average over eight other gauges. Gauge averaging significantly increases correlations compared to single gauge choices. For the marginal distributions (see figure 20), the number of hard instances (weight of the peak close to zero) is reduced by gauge averaging, but the distribution remains bimodal even after averaging over many gauge choices. While calibration errors enhance the bimodality, the convergence to a bimodal distribution after averaging many gauges shows that the bimodality is intrinsic and not solely due to calibration errors.

For scaling plots of the total effort we use geometric means of failure rates instead of arithmetic means of success probabilities. If the probability for finding a ground state for a specific percentile and gauge choice $g$ is denoted by $s$ then the probability of achieving a ground state at least once in $R$ repetitions of the annealing is

Splitting the $R$ repetitions into $R/G$ repetitions for each of $G$ gauge choices and denoting the success probabilities for a specific gauge choice by $s_{g}$ the total success probability becomes

which can be written in a form similar to equation (2)

by using the geometric mean of the failure rates to define $\overline{s}$ as

The higher percentiles grow much more rapidly in the absence of gauge averaging, as can be seen by comparing figure 21 to figure 4A) in the main text. Without gauge averaging the device actually fails to find the ground states for the 5% hardest instances in 1000 repetitions and those percentiles are thus not shown.

V.4 Correlations for the simulated annealer

To complement figures 2 and 3 in the main text we show correlations between the D-Wave device and a simulated classical annealer in figure 22 and the energy-success distribution $p(s,\Delta)$ of the simulated classical annealer in figure 23. As can already be expected from the different success distributions (figure 1 of the main text), the classical annealer does not correlate well with the D-Wave device and $p(s,\Delta)$ is significantly different.

VI Scaling

In this section we present the scaling of the time to solution for the various exact solvers used. Timings are from our reference CPU, except for the biqmac algorithm where the timings data from the spin glass server was used.

Figure 24 shows the time scaling of the exact solvers for instances without fields. Two of the algorithms, an exact belief propagation algorithm using bucket sort dechter1999bucket and the divide and conquer algorithm presented in section III.3 do not depend on the specific instance. The time to find a solution is the same for all instances with and without field. These algorithms make use of the structure of the chimera graph with a tree width of $\sqrt{N}$ Choi1 ; Choi2 and thus show the expected scaling ${\cal O}(\exp(c\sqrt{N}))$.

The tree width of a tree decomposition of a graph is the size of the largest vertex set of the tree (minus one according to some definitions Choi1 ; Choi2 ). A tree decomposition of a graph $G$ is a tree whose nodes are sets of vertices of $G$, and that satisfies the following properties:

Every vertex of $G$ is in at least one node of the tree.

For every edge $(v,w)$ in $G$ there is at least one node of the tree which includes both $v$ and $w$.

The nodes of the tree that contain any given vertex $v$ of $G$ form a connected subtree.

The cost of doing exact belief propagation (or dynamical programming) on a tree decomposition of a graph $G$ is exponential in the tree width. Belief propagation proceeds from the leaves down. A tree node $t$ with a set of vertices $\{v_{1},\ldots,v_{l}\}$ will calculate the minimum of the subgraph corresponding to the tree above $t$, conditional on all possible assignments to $\{v_{1},\ldots,v_{l}\}$. The cost of this calculation is exponential in $l$, and the tree width is the largest $l$. A tree decomposition for chimera graphs with tree width ${\cal O}(\sqrt{N})$ can be seen in Ref. klymko_adiabatic_2012, . Because a complete graph with ${\cal O}(\sqrt{N})$ vertices can be minor-embedded in the same chimera graph Choi2 , this scaling is optimal for exact belief propagation.

The time to solution of the other two solvers, akmaxsat akmaxsat , and the biqmac algorithm biqmac used in the spin glass server sgserver depends on the specific instance. The scaling for the various percentiles from the easiest (1% percentile) to the hardest (99% percentile) is shown in figures 25 and 26 for instances with and without local fields.

While the biqmac algorithm scales as ${\cal O}(\exp(c\sqrt{N}))$, akmaxsat scales worse and – unlike the other solvers – does not benefit from the limited tree width of the chimera graph. Unfortunately, the spin glass server constrained us to instances of up to $200$ spins. It would be interesting to see if a crossover is present, i.e., to check whether the best exact code scales better than ${\cal O}(\exp(c\sqrt{N}))$.

VI.2 Optimising the total annealing time

For the purpose of scaling comparisons we consider the total annealing time needed to find a ground state with a probability of 99%. Using Eq. (2) we find that the number of repetitions to find the ground state at least once with probability $p=0.99$ is

where $s$ is the probability for a given percentile. The total annealing time is $Rt_{f}$. Note that this expression assumes uncorrelated repetitions. On the D-Wave device we have seen statistically significant correlations between repetitions in 15% of the instances. In those instances the observed positive autocorrelations increase the length of “runs” of consecutive failures or successes on average by about a factor of two. This will result in an increase of the number of repetitions $R$ required on the device, as given by equation (6) above. The simulated annealers, on the other hand, show no detectable correlations.

In the main text we argued that extrapolating the total annealing time on the D-Wave device at a fixed annealing $t_{f}$ is tempting but misleading. It can only be used to estimate the performance for slightly larger problem sizes, but will not give the true asymptotic scaling. To see this, consider figure 27 where we show the scaling of the median time for a simulated annealer for various fixed annealing times (similar behaviour is observed for a simulated quantum annealer). The scaling at fixed annealing time $t_{f}$ increases only modestly for small system sizes, and then suddenly shoots upwards once $t_{f}$ becomes too short for the problem size.

What needs to be done instead is to optimise the total annealing time for each problem size $N$. To do this we vary $t_{f}$, measure the success probability $s$, calculate the required number of repetitions $R$, and plot the total annealing time $Rt_{f}$ as a function of $t_{f}$. The required number of repetitions $R$ diverges when the annealing time $t_{f}$ is too short, causing diabatic transitions. When the annealing time is too long it becomes disadvantageous to perform repetitions since a single run already optimises the success probability; however, Eq. (6) always yields $R\geq 1$. We thus pick the optimal $t_{f}$ as the time that minimises $Rt_{f}$ and use it in the scaling plots.

In figure 28 we show typical data for simulated classical and quantum annealing (where time is measured in number of spin flips for SA and spin updates for SQA) and for the D-Wave device (measured in $\mu$s). For the simulated classical and quantum annealer we observe that both the optimal time and the required number of repetitions $R$ increase exponentially. Note that a minimum is not found for the D-Wave device but instead a monotonic increase, indicating that even the fastest annealing time of the hardware of $t_{f}=5\mu$s is slower than the optimal time. All timings reported for the device are thus only upper bounds on the optimal time.

While the total annealing time for the D-Wave device can simply be given in microseconds, the annealing times of the simulated classical and quantum annealers depend on compiler options and the specific machine used. We thus give the total effort in terms of the number of attempted spin flips for the simulated classical annealer. For the simulated quantum annealer we specify the total effort as the number of spins updated multiplied by the inverse temperature $\beta$, to account for the complexity of updating the imaginary time world lines of the spins of length $\beta$.

Scaling plots of the total effort for instances with fields, complementing the results for instances without fields shown in the main text, are shown in figure 29. In panel A, we show the scaling of the D-Wave device. Again simulated classical annealing scales slightly better than simulated quantum annealing. Comparing the scaling for instances without fields (figure 4 in the main text) and instances with local random fields (figure 29), we see that problems with fields are not only easier for small problem sizes but the total annealing time also scales better when going to larger problem sizes.

VI.3 Detecting quantum speed up

Quantum speedup of a hardware quantum annealer can be detected by comparing the scaling of the total annealing times to that of the simulated classical or quantum annealer. To draw valid conclusions about a speedup one must ensure that the experimental annealing times are optimal. As we pointed out in the previous subsection, this is not the case for the current device: the annealing time of $5\mu$s is suboptimal, as demonstrated in figure 28. It follows that the inferred total annealing times are only upper bounds, and those bounds are worse for smaller problem sizes, which leads to scaling plots with an underestimated scaling. To see the pitfall it suffices to consider extremely long fixed annealing times where a single repetition $R=1$ might be enough. We would then see a constant time needed to find the ground state.

In other words, since the fastest possible annealing time $t_{f}=5\mu$s on the D-Wave device is longer than the optimal time for the considered problem sizes, our experimental data is in the initial transient regime of modest increase, and thus cannot be used for reliable extrapolation or determination of quantum speedup.

The optimal annealing time defined in the previous section (see figure 27) increases exponentially with problem size for SA. It is expected to increase exponentially with problem size also for QA. Therefore, assuming that the minimum programmable annealing time on future D-Wave devices does not increase too rapidly, we expect that perhaps already for a device with $N=512$ spins, or possibly $N=2048$ spins, the optimal annealing time $t_{f}$ for a single run will be within reach. This will allow us to determine the optimal total annealing times $Rt_{f}$ for QA. To detect quantum speedup one should then compare to the total effort of the simulated annealers divided by the number of spins $N$. This division is necessary to compensate for the trivial parallelism of the hardware annealer, which updates $N$ spins in parallel, since an analog classical annealing device would have the same parallelism. For completeness we provide the scaling plots of total effort in units of sweeps (spin flips divided by $N$) in figure 30.

VII Gap calculation

We finally describe how the excitation gaps are obtained using a method similar to that of Refs. Kashurnikov1999 ; Young2008 . For simplicity we consider the transverse field Ising spin glass Hamiltonian without fields,

where $\Gamma$ is the strength of a transverse magnetic field and for our instances $J_{ij}=\pm 1$.

The excitation gap can be obtained from the connected correlation function in imaginary time

where $\hat{O}$ is an observable with non-vanishing matrix element between the ground state and first excited state. The correlation function is given by a sum:

where $\Delta_{i}$ is a gap to the $i$th eigenvalue above the ground state; $i=0$ corresponds to the lowest gap. Only one term survives when $\tau$ is large enough:

Thus $\Delta\equiv\Delta_{0}$ can be obtained by fitting $C(\tau)$ at large values of $\tau$. We use periodic boundary conditions in the imaginary time direction. In this case, $C(\tau)$ can be efficiently calculated at discrete points using the Fast Fourier Transformation.

The observable $\hat{O}$ must be chosen carefully because for a poorly chosen observable, $c_{0}$ can be much smaller than $c_{1},c_{2},\ldots$ leading to very small values of $C(\tau)$ (comparable to statistical noise) at large $\tau$ and making the gap extraction very difficult. This issue is discussed in Ref. PhysRevE.85.036705, . We use a simple observable, the local magnetisation $\sigma_{j}^{z}(\tau)$ and its correlation function, given by

where the sum runs over all the spins $N$.

We use the continuous time algorithm described in section III.2 and “anneal” the system from $\Gamma=3$ to small values of $\Gamma$ in steps of $0.1$. The simulation at a transverse field $\Gamma_{i}$ is started from the final configuration obtained in the previous simulation $\Gamma_{i-1}=\Gamma_{i}-0.1$, except for the initial easy simulation at a strong transverse field $\Gamma=3$, where the simulation is started from a random configuration. $200000/\Gamma$ Monte Carlo sweeps (one Monte Carlo sweep consists of $2N$ site updates) are performed for equilibration before measurements. Simulations are done at inverse temperatures $\beta=100$ and $\beta=200$.

Extraction of the gap for many instances and many values of $\Gamma$ can be a cumbersome task. To automate the process, we use the following approach. We obtain the gap by fitting the correlation function to the exponential function given by Eq. (10) in the range from $\tau_{0}$ to $\tau_{0}+5/J$, where $\tau_{0}$ is chosen empirically in such a way that $C(\tau)-C(\beta/2)>s_{0}$ for all $\tau<\tau_{0}$, $s_{0}=0.008f(\Gamma)$ and $f(\Gamma)=0.2+0.3\Gamma$. Let us denote this gap as $\Delta_{0}$. To obtain the error bars, we calculate another two gaps $\Delta_{1}$ and $\Delta_{2}$ by performing two extra fits with $s_{1}=0.006f(\Gamma)$ and with $s_{2}=0.011f(\Gamma)$. Then the error bar is just $\max(\Delta_{2}-\Delta_{0},\Delta_{0}-\Delta_{1})$. This procedure can be fully automated. We show an example of such a fit in figure 32.

In figure 4 of the main text we showed the excitation gap in units of temperature. In figure 33 we instead show the gap in units of $B(t)$, so that the gap shown is a property of the model only and not of the annealing schedule. While the physical temperature is constant ($17\,\textrm{mK}$), it increases as a function of $\Gamma$ when shown in units of $B(t)$. The gap and temperature are both shown in figure 33 for the two instances presented in figure 4 of the main text, as well as for four “easy” (left column) and four “hard” (right column) additional instances. Note that the gap closes trivially around $\Gamma=2.2$, related to a global $Z_{2}$ symmetry breaking. Once the gap becomes very small it can no longer be detected by our procedure since the decay of $C(\tau)$ becomes too slow and is indistinguishable from a constant. Our procedure then picks up the gap to the next excited state, which results in an apparent jump of the gap to a bigger value. Generally, all gaps shown here are upper bounds for the gap to the lowest excited state.

The gap always closes also in the limit of a weak transverse field $\Gamma\longrightarrow 0$, where multiple ground states become degenerate. Some of the instances, however, have an additional small gap (relative to the temperature) that can be associated with an avoided level crossing. These instances are “hard” for both the D-Wave device and SQA. Usually, the instances that do not have such a small gap are “easy” for both the D-Wave device and SQA. However, note that the last “easy” instance shown in figure 33 does have a small gap. This shows that a small gap is not a sufficient condition for making a problem “hard”. The fact that this last instance is “easy” can have multiple reasons, for example a second avoided level crossing with a small gap that takes the annealing back to a ground state or an avoided level crossing with a small Hamming distance from the ground state, from where thermal relaxation is not difficult.

References