{'fraction_non_alphanumeric': 0.029482992290430804, 'fraction_numerical': 0.010620098202581706, 'mean_word_length': 4.66448527561037, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Annealing amid Local Ruggedness and Global Frustration

A recent Google study [Phys. Rev. X, 6:031015 (2016)] compared a D-Wave 2X quantum processing unit (QPU) to two classical Monte Carlo algorithms: simulated annealing (SA) and quantum Monte Carlo (QMC). The study showed the D-Wave 2X to be up to 100 million times faster than the classical algorithms. The Google inputs are designed to demonstrate the value of collective multiqubit tunneling, a resource available to D-Wave QPUs but not to simulated annealing. But the computational hardness in these inputs is highly localized in gadgets, with only a small amount of complexity coming from global interactions, meaning that the relevance to real-world problems is limited. In this study we provide a new synthetic problem class that addresses the limitations of the Google inputs while retaining their strengths. We use simple clusters instead of more complex gadgets and more emphasis is placed on creating computational hardness through frustrated global interactions like those seen in interesting real-world inputs. The logical problems used to generate these inputs can be solved in polynomial time [J. Phys. A, 15:10 (1982)]. However, for general heuristic algorithms that are unaware of the planted problem class, the frustration creates meaningful difficulty in a controlled environment ideal for study. We use these inputs to evaluate the new 2000-qubit D-Wave QPU. We include the HFS algorithm---the best performer in a broader analysis of Google inputs---and we include state-of-the-art GPU implementations of SA and QMC. The D-Wave QPU solidly outperforms the software solvers: when we consider pure annealing time (computation time), the D-Wave QPU reaches ground states up to 2600 times faster than the competition. In the task of zero-temperature Boltzmann sampling from challenging multimodal inputs, the D-Wave QPU holds a similar advantage as quantum sampling bias does not seem significant.


I. INTRODUCTION
Quantum annealers are designed to take advantage of quantum tunneling to find good solutions to hard optimization problems. When constructing a family of synthetic inputs to test the potential of a quantum annealing platform, one should therefore ensure that the inputs a) are such that solvers can benefit from quantum tunneling, and b) are hard optimization problems with global frustration.
For a solver to benefit from quantum tunneling, the energy landscape associated with the input must have tall, thin energy barriers. For an input to be computationally hard, the input must have constraints that interact with each other in nontrivial ways.
Quantum processing units (QPUs) developed by D-Wave Systems that use the quantum annealing algorithm have been commercially available since 2011. These QPUs solve Ising model inputs defined on the underlying working graph of the chip. There have been various efforts to evaluate the performance of the D-Wave systems using synthetic inputs generated randomly from different distributions, or input classes.
This study has two main contributions: to propose a new problem class ideal for evaluating D-Wave QPUs, and to use this problem class to evaluate the 2000-qubit D-Wave QPU.

A. Proposing a new problem class
Previous evaluations of D-Wave QPUs have used problem classes that benefit either too little or too much from quantum tunneling to be ideal for evaluating quantum annealers.
On one side of this spectrum we have problems such as random unstructured ±1 problems on the Chimera topology native to D-Wave QPUs. These were used by Rønnow et al. [1] in their evaluation of the D-Wave Two QPU in 2014, but they are now known [2] to lack a finite-temperature phase transition, meaning that quantum tunneling is unlikely to play a significant role when solving them.
On the other side of the spectrum, Denchev et al. [3] recently introduced an input class designed to benefit immensely from quantum tunneling. We refer to these inputs as Google problems. Their study showed a massive speed increase (up to 100 million times faster) of a D-Wave 2X system over simulated annealing (SA) and quantum Monte Carlo (QMC), also known as simulated quantum annealing. This provided strong evidence for the ability of quantum annealing to leverage quantum tunneling in a computationally relevant way. However, the spinglass backbones of the Google problems are easy to solve, meaning that a) the problems have limited relevance to real-world problems, and b) certain cluster-detecting algorithms can solve them with relative ease [4].
In this study we provide a problem class that aims to retain the advantages of Google problems while being more reflective of real-world problems. They are more reflective of real-world problems because, rather than relying too heavily on finely-tuned gadgets, they derive much of their computational hardness from larger spin-glass backbones with planted frustration.
Our problems are synthetic and are easy to solve using knowledge of the problem class. 1 More specifically, since the logical problems are Ising models on a planar lattice without fields, they are solvable in polynomial time [5]. However, frustration in the logical problems creates meaningful difficulty for heuristic methods that are unaware of the planted problem class. Further, the inputs have properties such as tunable ruggedness that make them useful for the evaluation of quantum annealing and classical approximations thereof. In this way, they are similar to Kauffman's NK model that has proved very useful in the analysis of evolutionary algorithms [6][7][8].

B. Evaluation of the 2000-qubit D-Wave QPU
We use this new problem class to evaluate the latestgeneration D-Wave QPU. We measure its performance in absolute terms and we analyze its response to the ruggedness parameters of the problem class.
The software competition we consider is much stronger than that considered by Denchev et al. [3], and includes GPU implementations of SA, QMC, and SVMC, and also includes Selby's implementation [9,10] of the Hamze-de Freitas-Selby (HFS) algorithm [9,11]. In the study of Mandrà et al. [4] that used a wide array of algorithms to solve Google problems, Selby's implementation of HFS was the fastest software solver in terms of both scaling and absolute speed.
We find that the D-Wave QPU is able to find ground states up to 2600 times faster than the software competition. We also consider the problem of sampling from ground states and find that the D-Wave QPU maintains a similar advantage and does not struggle to find a diverse set of optimal solutions. The remainder of the paper is organized as follows. In Section II we provide a description of the 2000-qubit D-Wave system and a history of D-Wave QPUs. In Section III we present the problem class analyzed in this paper and discuss the concept of ruggedness and its relevance to optimization problems. In Section IV we discuss the software solvers used in our evaluations, as well as notable solvers that were not suitable. In Section V we present our experimental results on optimization. In Section VI we present our experimental results on sampling from ground states. In Section VII we provide further discussion and conclude the paper.

II. D-WAVE QUANTUM PROCESSING UNITS
We start with an overview of D-Wave design features and introduce notation that will be used throughout. For details about underlying technologies see Bunyk et al. [12], Dickson et al. [13], Harris et al. [14], Johnson et al. [15] or Lanting et al. [16].
A. Ising minimization D-Wave annealing-based quantum processors are designed to find minimum-cost solutions to the Ising minimization (IM) problem, defined on a graph G = (V, E) as follows. Given a collection of fields h = {h i : i ∈ V } and couplings J = {J ij : (i, j) ∈ E}, assign values from {−1, +1} to n spin variables s = {s i } so as to minimize the energy function The spin variables s can be interpreted as magnetic poles in a physical particle system; in this context, negative J ij is ferromagnetic and positive J ij is antiferromagnetic, the optimal solution is called a ground state, and nonoptimal solutions are excited states. IM instances can be trivially transformed to Quadratic Unconstrained Boolean Optimization (QUBO) instances defined on integers s = {0, 1}, or to Maximum Weighted 2-Satisfiability (MAX W2SAT) instances defined on Booleans s = {true, false}, all of which are NP-hard.

B. Chimera topology
The native connectivity topology for the D-Wave QPU is based on a C 16 Chimera graph containing 2048 vertices (qubits) and 6016 edges (couplers).
A Chimera graph of size C s is an s × s grid of Chimera cells (also called unit tiles or unit cells), each containing a complete bipartite graph on 8 vertices (a K 4,4 ). Each vertex is connected to its four neighbors inside the cell as well as two neighbors (north/south or east/west) outside the cell: therefore every vertex has degree 6 excluding boundary vertices.
In this study, as in others, we vary the problem size using square subgraphs of the full graph, from size C 4 (128 vertices) up to C 16 (2048 vertices). Note that the number of problem variables n = 8s 2 grows quadratically with Chimera size. The reason we measure algorithm performance as a function of the Chimera size and not the number of qubits is that problem difficulty tends to scale exponentially with the Chimera size, i.e., with the square root of the number of qubits, since the treewidth of a Chimera graph C s is linear in s [17,18].
Because the chip fabrication and trapped magnetic flux leave some small number of qubits unusable, each processor has a specific hardware working graph H ⊂ C 16 . The qubit yield-the fraction of qubits that are operational-is typically around 98% for the 2000-qubit D-Wave system whereas 95% was typical for the D-Wave 2X. The working graph used in this study has 2035 working qubits out of 2048.
C. Quantum annealing D-Wave processors solve Ising problems by quantum annealing (QA) in the form proposed by Kadowaki and Nishimori [19]. The QA algorithm is implemented in hardware using a framework of analog control devices to manipulate a collection of qubit states according to a time-dependent Hamiltonian shown below.
QA carries out a gradual transition in time t : 0 → t a , from an initial ground state in H init , to a state described by the problem Hamiltonian H prob = The problem Hamiltonian matches the energy function (1), so that a ground state for H prob is a minimum-cost solution to E(s).
QA is closely related to adiabatic quantum computing (AQC). The AQC model of computation was proposed by Farhi et al. [20] who showed that if the transition is carried out slowly enough the algorithm will find a ground state (i.e., an optimal solution) with high probability.
Theoretical guarantees about solution times for quantum algorithms (found in [20]) assume that the computation takes place in an ideal closed system, perfectly isolated from energy interference from ambient surroundings. The 2000-qubit D-Wave chip is housed in a highly shielded chamber and cooled to near absolute zero; nevertheless, as is the case with any real-world quantum device, it must suffer some amount of interference, which has the general effect of reducing the probability of landing in a ground state. Thus, theoretical guarantees on performance may not apply to these systems. We consider any D-Wave QPU to be a heuristic solver, which requires empirical approaches to performance analysis.
The D-Wave QPU studied here contains 2035 active qubits (quantum bits) and 5912 active couplers made of microscopic loops of niobium connected to a large and complex analog control system via an arrangement of Josephson Junctions. Thermometry on the refrigerator of the D-Wave QPU and fits of single qubit measurements to a thermodynamic model indicate that T 15 mK. When cooled to temperatures below 9.3 K, niobium becomes a superconductor and is capable of displaying quantum properties including superposition, entanglement, and quantum tunneling. Because of these properties, the qubits on the chip behave as a quantum mechanical particle process that carries out a transition from initial state described by H init to a problem state described by H prob [13,16,21]. b. Anneal. Carry out the QA algorithm. Anneal time t a can be set by user to some value 5 µs ≤ t a ≤ 1000 µs.
c. Read. Record qubit states to obtain a solution; denote the elapsed readout time t r .
d. Repeat. Repeat steps b and c k times to obtain a sample of k solutions.
We define sample time t s and total time T as follows: For the D-Wave system studied in this paper, the median programming time t i is 9.5 ms and the median readout time t r is 123 µs.
In this study, both for software solvers and for the D-Wave QPU, we typically report annealing time rather than total time. Annealing time is the measure of the algorithm proper, and measuring total time often obscures trends in data. Scaling plots are particularly susceptible to this because the overhead of programming time makes scaling-typically presented on a semilog plot-look totally flat except for an uptick at the very largest problem sizes. Further, we are most interested in the future potential of D-Wave QPUs, and we expect that programming time and readout time will be reduced to small fractions of their current values; minimum annealing times will similarly be reduced, allowing us better control over the algorithm parameters. For reference, since many people will be interested in total wall clock time, rather than annealing time, a 1000 times speedup over software solvers in annealing time, typical for the D-Wave QPU in this study, translates roughly to a 30 times speedup in total wall clock time including programming and readout.
System characteristics of D-Wave QPUs such as yield can vary within a generation. If we compare this specific 2000-qubit D-Wave system to the specific D-Wave 2X QPU studied in 2015 [22], programming time has decreased by 20%, readout is three times faster, and yield has improved from 95% to 99%.

A. Ruggedness and clusters
Ruggedness is a feature of certain optimization problems-more specifically their energy landscapes-characterized by tall energy barriers and many local optima [23,24]. Typically, rugged problems are harder to solve, particularly with Markov chain Monte Carlo (MCMC) methods [25,26]. In the late 1980s, when ruggedness was first being explored in the context of evolutionary biology and bio-inspired computing, Kauffman's NK model was put forward as a model with tunable ruggedness inspired by genetic fitness functions under varying degrees of epistasis, or how many other genetic loci affect the fitness contribution of a given locus [6][7][8].
The tunable ruggedness of the NK model has proved very valuable in the study of optimization heuristics, particularly evolutionary algorithms [27].
Closely related to ruggedness is the analysis of spin overlap, in which landscape features are inferred from the distribution of overlap of two random states sampled from the Boltzmann distribution [28,29]. Tall, thin peaks in the spin overlap distribution tend to correspond to tall, thin energy barriers; the presence of these features correlates not only with ruggedness and classical hardness, but also with applicability of quantum annealing, since quantum tunneling is likely to be a useful computational resource in the presence of these tall, thin barriers. Zhu et al. [30] have used spin overlap features to predict whether a problem can be solved by QA more efficiently than by SA, showing promising preliminary results for optimization problems such as weighted partial MAX-2SAT, minimum vertex cover, satisfiability, graph partitioning, circuit fault diagnosis, and certain spin-glass instances [30,31]. This work points to the potential of quantum annealers to have a place in portfolio solvers [32] and hybrid algorithms running on heterogeneous computing systems alongside CPUs, GPUs, and other coprocessors [33,34].
To induce ruggedness using tall, thin energy barriers, Denchev et al. [3] used ferromagnetically coupled unit tiles as clusters. Flipping such a cluster in the absence of fields or external couplings involves jumping over or tunneling through an energy barrier that is 16 Ising units high and has a width of 8 in Hamming space. Denchev et al. [3] actually go beyond using single-tile clusters and use two-tile gadgets studied previously by Boixo et al. [21]. This gadget is made up of two clusters that form a deceptive trap to draw annealers into a local minimum using local fields; annealers must then go over or through an energy barrier to reach the gadget's ground state.
Instead of using two-cluster gadgets, we simply use single-cell ferromagnetic clusters to induce ruggedness, leaving us with a simpler problem class.

B. FCL problem generation
We create local ruggedness by treating unit cells of the Chimera graph as ferromagnetically-coupled clusters. We create global frustration by joining these clusters together using a problem generated on the logical graph of clusters. This creates an energy landscape that is macroscopically interesting and in which the clusters induce wells separated by tall, thin energy barriers.
The logical graph of clusters is a square lattice, with a logical 16×16 lattice of clusters spanning the working graph of a 2000-qubit D-Wave QPU. 2 The problems we generate on the logical graph are frustrated loop problems, constraint satisfaction problems first used in the evaluation of D-Wave QPUs by Hen et al. [35] and modified to allow precision limits by King et al. [36].
We refer to the final inputs as frustrated cluster loop (FCL) problems. For a given Chimera graph G C that may or may not have missing qubits or couplers, an FCL problem is generated from three parameters, α (the clauses-to-variables ratio), ρ (the range, or precision), and R ≥ ρ (the ruggedness) as follows: a. Define each unit cell as a logical spin if it has no missing qubits or couplers. Use c(v) to denote the logical spin index corresponding to qubit v.
b. Wherever all four couplers connecting two logical spins are present, define these couplers as a logical coupler.
c. Define the logical graph G L as the graph comprising the logical spins and logical couplers.
d. Generate a range-bounded frustrated loop problem Hamiltonian (h L , J L ) on G L using parameters α and ρ as per King et al. [36] (note that h L is the zero vector).
e. Define the native Chimera Hamiltonian (h C , J C ) with h C as the zero vector and J C as: It is worth repeating that these Hamiltonians have no fields (i.e., h L and h C are both zero vectors). Note also that in-tile couplings in J C are all −1 and inter-tile couplings take values in Since we ensure that ρ ≤ R, and ρ ≥ |j| for any logical coupling j, all couplings in J C are in the range The logical frustrated loop problems may be disconnected and have multiple components; we reject such disconnected inputs at generation time.
While these problems are large enough to span the entire working graph of the latest D-Wave QPUs, the repetition code inherent in logical couplers and spins makes them relatively robust to analog errors [37].

C. Problem class parameters
The FCL problem class has three parameters: the clauses-to-variables ratio α, the range ρ, and the ruggedness R. We would like to restrict our experiments to the most interesting region of the parameter space.
First we aim to determine the value of α that maximizes the difficulty of the logical problem. If α is too low, a problem is underconstrained and is easy to solve. If α is too high, the planted solution is expressed too strongly and the problem's features approach those of a ferromagnet, making it easy to solve. The difficulty of the logical problem depends only on α and ρ, not on R. For various values of ρ, we perform a sweep of α to determine the value that maximizes the hardness of the logical problem (see Figure 1).
The impacts of ρ are more nuanced. First, the value of ρ provides an upper bound on the limit of α because packing in more loops eventually raises the maximum coupler range. Second, for a fixed value of α, problems with a lower ρ value have their loops spread out more evenly over the logical spins. Finally, for the native problem, coupler values are scaled down by a factor of R ≥ ρ so that inter-tile couplings are in the range [−1, 1] (in-tile couplings are always −1). Thus higher values of ρ constrain R to be higher, and make problems more locally rugged relative to the global Hamiltonian. In Section V, we attempt to deconvolve the impacts of ρ and R.
The ability to tune the ruggedness of the inputs by varying R, either by specifying R = ρ and varying ρ, or by varying R independently, gives FCL problems an additional degree of utility, particularly when assessing the value of quantum tunneling and the potential of quantum annealing. Varying ρ and specifying R = ρ makes problem generation simpler by reducing the number of free parameters whereas fixing ρ and varying R allows us to isolate the impact of ruggedness without altering the complexity of the logical problem.

D. Confirming correlation between ruggedness and classical hardness
We expect to see a positive correlation between ruggedness and classical hardness. Here we characterize classical hardness using the decorrelation time of an MCMC procedure. To validate this assumption we measure the decorrelation time for a parallel tempering (PT) procedure that uses the Metropolis algorithm in combination with the standard replica exchange rule [38]; we use the autocorrelation of temperature as our measure of decorrelation [39]. For more details of this method, see Appendix A.

IV. SOFTWARE SOLVERS
The four software solvers we consider are GPU implementations of SA, QMC, and SVMC, and Selby's CPU implementation of HFS [9,10].
Recent studies of D-Wave QPUs have not included GPU-based software solvers despite the fact that SA is very amenable to GPU implementation [40]. The addition of GPU solvers is a significant raising of the bar in terms of software competition, and means that solvers that can be implemented on GPUs have taken a leap forward relative to solvers that cannot. Run on modern hardware, our GPU-based algorithm implementations are roughly 1000 times faster than the corresponding single-core CPU implementations.
Mandrà et al. [4] analyzed the performance of a diverse set of solvers on the inputs of Denchev et al. [3]. However the lack of GPU implementations of these solvers means that most are unlikely to be competitive in an absolute sense. Indeed, of the three classes of solvers that they study, only sequential algorithms, which they find to have the worst performance, have the massive parallelizability and low memory requirements that make efficient GPU implementations possible. It is also possible to implement SA in a field-programmable gate array (FPGA), but the additional speedup over GPU implementation is limited and generally not worth the increased cost of hardware.
In Appendix B we give further details of the software solvers and parameterizations we used. We also discuss algorithms we omitted because of prohibitive runtimes.

V. OPTIMIZATION
We measure the expected time to solution (TTS) of different solvers on the inputs, calculated as TTS = time per anneal ground state probability .
We consider only annealing times and exclude programming and readout times from our analysis as these are not part of the algorithms proper.
For a given value of ρ, we choose α to maximize the difficulty of the logical problem (see Figure 1). For each selection of ρ and R, we generate 100 FCL problems at each problem size and solve each problem with each solver.

A. Varying ruggedness via logical complexity
In our first experiment, we vary ρ and set R = ρ. In this case the only free parameter ρ controls both the ruggedness and the logical complexity of the inputs. Time-to-solution plots are shown in Figure 3. At the largest problem size, the D-Wave QPU is three orders of magnitude faster than the fastest software solver for each value of ρ. D-Wave's speedup over software peaks at 2600 times for ρ = 4.
As ρ increases, the impact of local ruggedness increases as the logical Hamiltonian is compressed relative to the local wells induced by the clusters. The performance of SA drops off sharply while the performance of DW and QMC declines gracefully. The performance of HFS decreases only very slightly. HFS is not affected by the local ruggedness because it is tailored to the Chimera topology and uses updates that contain entire clusters; the performance degradation is due to the slight increase in logical problem hardness.
All solvers except HFS have strictly convex scaling curves because the anneal lengths are optimized for the largest problem size and are too long for the smaller problems. HFS does not use fixed-length anneals and ends up using shorter anneals on smaller inputs.
Though true scaling is masked by the inability to optimize parameters for smaller inputs [1,41], we note that the performance of the D-Wave QPU scales at least as well as the software solvers between the two largest problem sizes.

B. Varying ruggedness by scaling
In our second experiment, we fix ρ = 3 and vary the ruggedness R. This keeps the logical complexity constant, allowing us to isolate the impact of ruggedness on the various solvers. Here we consider only the largest problem size having a 16 × 16 logical lattice.
Consistent with our findings when varying ρ, tuning the ruggedness directly by varying R increases difficulty dramatically for simulated annealing and less so for other solvers (see Figure 4). Excluding HFS, whose behavior is constant in this example, the work required by a solver essentially scales according to its quantumness. The D-Wave QPU is most capable of dealing with ruggedness. QMC-the most faithful classical simulation of quantum annealingcomes next, followed by SVMC, which is a meanfield approximation to QMC. Bringing up the rear is SA, a simulation of a fully classical process.
The improved scaling (versus ruggedness) of QMC HFS deals with these energy barriers using exponential brute force; therefore the parameter R does not affect its performance. Markers indicate medians (over 100 inputs) and error bars indicate 95% confidence intervals for the median.
over SVMC indicates that crucial information is being lost in the mean-field approximation. The improved scaling of QA (i.e., the D-Wave QPU) over QMC may indicate that QMC is failing to faithfully simulate the dynamics of the QA processor, or it may simply be an artifact of our inability to use faster D-Wave anneals. This bears further investigation using future D-Wave QPUs with faster annealing times, again utilizing the tunable ruggedness of FCL problems.

VI. SAMPLING
The ability of an Ising solver to sample diverse optima has both practical and theoretical importance. Ground state sampling in combinatorial problems is the basis for construction of space-efficient SATbased membership filters [42,43]. The associated complexity class, #P -the counting analog of NPhas been the subject of extensive research in theoretical computer science since the 1970s [44]. Sampling from the Boltzmann distribution, in which states with equal energy are sampled with equal probability, is of particular interest in machine learning. Boltzmann samples are used to train Boltzmann machines, a task known to be both hard and useful [45].
While machine learning applications typically depend on finite-temperature Boltzmann sampling, using near-optimal states as well as optimal states, we focus on zero temperature sampling to simplify our investigation. This saves us from having an additional input parameter β-the inverse temperature-that we would have to either set arbitrarily or determine empirically. Empirical estimation of β can be challenging [46] and basing the target β on the output of a solver would arguably give that solver an unfair advantage.

A. Sampling from all valleys
The expected time required for a solver to find all ground states of a problem is known, both in the equiprobable case and the biased case [47]. In the case of an Ising spin problem, ground states often lie in connected valleys in Hamming space, and given one ground state in the cluster it is easy to find the rest. We therefore adopt a more practical metric based on the time required to sample all valleys of ground states.
In ground states of FCL problems, all clusters have their spins in agreement; therefore the distance between any two ground states in the native Hamming space is a multiple of 8. However, ground states can be adjacent (i.e., differ by a single spin) in the logical space. We define a valley as a set of ground states that are connected in the logical Hamming space. While it is nontrivial to move from one state to another in the same valley because of the single tall, thin energy barrier, it can still be done with a modest amount of postprocessing. We also note that, since FCL problems do not have fields, ground states come in antipodal pairs, 3 and by extension so do valleys. We treat each pair of antipodal valleys as a single valley since it is trivial to move from one to the other.
With valleys defined in this way, we define the timeto-all-valleys (TTAV) metric as the expected amount of annealing time required to draw at least one sample from each valley. This metric captures the hardest part of sampling from these distributionsfinding ground states in every mode-and ensures a diverse set of solutions. With at least one sample from each valley, it is possible to find all ground states using only a modest amount of postprocessing. The TTAV metric is most meaningfully interpreted relative to TTS since hitting valleys directly depends on hitting ground states.

B. Mining for interesting valley structure
Sampling from all valleys is not always much harder than finding a single ground state-an input may have only a single valley or may have valleys that are all very close in Hamming space. We wish to generate inputs with multiple valleys that are wellseparated. Sampling from distributions with multiple, well-separated valleys is particularly hard [48] and has important applications such as classification using deep Boltzmann machines [49].
Because we define valleys as clusters of ground states in the logical space, analyzing the valley structure of an input is tractable. The largest logical graphs are 16 × 16 lattices having treewidth 16, so solving a logical problem using dynamic programming and returning some fixed number of ground states typically takes less than a second. This allows us to mine for inputs having interesting valley structure.
We quantify interesting valley structure using the distribution of spin overlap P (q) [28,29] at zero temperature, i.e., for two ground states sampled uniformly with replacement, what fraction of spins do they have in common? The random variable P (q) takes values in the range [−1, 1]. For inputs without fields the distribution is symmetric about zero; we can therefore consider the distribution of the absolute value P (|q|). We define the mean overlap as the expectation of P (|q|). Inputs with mean overlap near 1 tend to resemble ferromagnets-if there are multiple valleys they will be close together. Inputs with lower mean overlap tend to have valleys that are well-separated.
Inputs that are hard to sample from have multiple valleys that are well-separated. We mine for such inputs as follows. First we reject any input with more than 1000 ground states, as these slow down our analysis and may be too easy. Second, we reject any input that does not have at least 4 valleys since we want valley collection to be nontrivial. Finally, we reject any input with a mean overlap of 0.7 or higher since we want valleys to be well-separated.

C. Sampling results
We generated problems at the 16 × 16 lattice size with α = 0.85 and ρ = R = 6 and mined them for interesting valley structure as described above. We generated 50,000 inputs and rejected all but 74. This gave us an acceptance rate of roughly 0.15% of inputs. We sought to answer the question, after In certain situations quantum annealing in the transverse-field Ising model is subject to inherent sampling bias [50][51][52][53], although that does not prove to be a significant problem here. While the slope of the D-Wave curve is slightly less steep than the HFS curve, indicating that its samples might be less diverse, the D-Wave QPU still manages to outperform the competition by about three orders of magnitude.

KL-divergence of valley distributions
The TTAV results shown in Figure 5 fail to address a specific fear-that in a significant minority of inputs there are valleys that the D-Wave QPU would be simply unable to find due to quantum sampling bias.
To address this, we calculate for each (solver, problem) pair the KL-divergence between empirical valley distributions and exact valley distributions (i.e., relative valley sizes). KL-divergence is an asymmetric measure of the distance between two probability distributions; we calculate it such that it is infinite if the solver fails to see all valleys, i.e., where P (v) is the true Boltzmann probability of valley v and P (v) is the sample estimate of P (v), conditioned on samples being ground states. This KL-divergence measure includes two types of error. First, there is a distributional error, since each solver samples from a distribution that differs from the Boltzmann distribution. Second, there is a sample size error, since our sample estimate has finite size and therefore differs from the solver's true distribution. In this context it is appropriate to include both types of error. Figure 6 shows histograms of KL-divergence for the different solvers. For these FCL problems, fears of valleys suppressed by quantum sampling bias are unfounded. The D-Wave QPU has a superior KLdivergence distribution than any of the software solvers even when annealing for three orders of magnitude less time. On the single input for which the D-Wave KLD was infinite because at least one valley was never seen, it was also infinite for all other solvers.

Raw error on model marginals
The TTAV metric and valley distributions can be thought of as representing what sample quality would look like with postprocessing. We would also like a more raw metric that does not have this implicit postprocessing. For this we consider marginals of the zero-temperature Boltzmann distribution. Specifically, we consider the spin-spin expectations, i.e., for each coupler, what is the expected product of the two incident spins in the zerotemperature Boltzmann distribution? The L 1 error on these marginals (i.e., the empirical estimates of spin-spin expectations minus true expectations) is  a well-established metric of interest in the study of undirected graphical models (see, e.g., [54]).  Figure 7 shows the decay in error as more samples are taken. We measure errors in the logical space, so couplers within a cluster are ignored. Again, the D-Wave QPU is roughly three orders of magnitude faster than the fastest software solver. We expect that the error on marginals will decay to a certain point, then plateau, with the level of the plateau corresponding to the bias of the solver. As with the TTAV data, potential concerns about the bias of quantum annealers do not seem to play out here. The L 1 error of the D-Wave QPU is still decaying after 100,000 samples; at this point in time (0.5 s) the D-Wave QPU has achieved a median error of less than 1%, a value that software solvers fail to reach in 100 s to match.

VII. CONCLUSIONS
We have introduced a class of synthetic inputs on which to evaluate the performance of annealingbased QPUs, specifically the QPUs developed by D-Wave. This problem class is more representative of real-world problems and provides an alternative to the Google problems of Denchev et al. [3], which were more highly tuned to highlight the utility of quantum tunneling.
The D-Wave QPU is up to 2600 times faster than all software solvers considered and typically on the order of 1000 times faster at the largest problem size. These software solvers now include GPU implementations of SA, QMC, and SVMC as well as a CPU implementation of HFS, making the competition much stronger than that analyzed by Denchev et al. [3]. The set of software solvers we used was representative, but not exhaustive-in particular, we hope to include more of the solvers used by Mandrà et al. [4] in future studies.
These inputs have tunable ruggedness controlled either by the range parameter ρ or by the scaling parameter S. Of particular interest is the fact that QMC performed better relative to SA when the ruggedness is increased, and physical quantum annealing performed better still.
We also evaluated the 2000-qubit D-Wave QPU on the task of zero-temperature Boltzmann sampling, i.e., sampling uniformly from ground states. While concerns have been raised that quantum and analog sampling bias might make it difficult for quantum annealers to sample from Boltzmann distributions [50][51][52], there was little evidence for such a struggle in this study. In several metrics considered, the 2000-qubit D-Wave QPU maintains its speed advantage of roughly three orders of magnitude and provides sample diversity that is as good as, or better than, the software competition. initial sequence estimator [56], we average over the autocorrelations on the n available chains to reduce noise.
Appendix B: Details of software solvers 1. Classical hardware Table I contains the specifications for the classical processors used. CPU algorithms were run singlethreaded on one core each; multiple workers used cores concurrently for independent jobs. GPU algorithms are highly parallelized and each GPU job uses the entire GPU.

Included software solvers
Our experiments and analyses focus on four algorithms commonly used in comparisons with D-Wave QPUs-three that are highly amenable to GPU implementation and one that is highly tailored to the Chimera topology.

a. Simulated annealing
Simulated annealing [57] is a simulation of thermal annealing that is widely used as an optimization algorithm. It is the classical analog to quantum annealing. Since simulated annealing is a simple algorithm with very low memory requirements and a high degree of parallelizability, it is ideal for implementation on a GPU.

b. Quantum Monte Carlo
Quantum Monte Carlo, also known as simulated quantum annealing, is a classical approximation to quantum annealing. For the algorithm to work efficiently on a GPU, we implement the discrete time variant of QMC and fix the number of Trotter slices at 64 so that a worldline can be packed as bits in a word.
While the continuous time variant of QMC is a more faithful simulation of quantum annealing, in particular serving as a bias-free sampler that approaches the quantum Boltzmann distribution in the limit, it has been shown that discrete time QMC can have superior performance as an optimizer [58].

c. Spin vector Monte Carlo
Spin vector Monte Carlo (SVMC), also known as the O(2)-rotor model, is a mean-field approximation to QMC. SVMC can be thought of as falling between SA and QMC. Proposed for use as an approximation to D-Wave QPUs by Shin, Smith, Smolin, and Vazirani [59], it is also known as the SSSV algorithm. We use a GPU implementation of SVMC that is a minor modification of our implementation of SA.

d. Hamze-de Freitas-Selby
The Hamze-de Freitas-Selby (HFS) algorithm optimizes by repeatedly optimizing the spins in a lowtreewidth induced subgraph of the input, subject to the rest of the input being fixed. The subgraph over which the input is optimized changes at each step. The HFS algorithm is a greedy search algorithm in which reassignment of many variables is considered at once. We used Selby's implementation [9,10] that is heavily tailored to the Chimera topology; we modified this solver to return each stopping state for consistency with the other heuristic solvers. We note that the HFS algorithm cannot be efficiently implemented on GPU because the memory requirements are too high.

Excluded software solvers
In addition to these four algorithms, we considered several other software solvers that were prohibitively slow; due to limited time and resources it was not feasible to perform the long software runs needed to optimize parameters and determine ideal performance.
a. Nontailored HFS We tested an implementation of HFS that, rather than using subgraphs tailored to the Chimera topology, is topology-agnostic and generates subgraphs dynamically. This nontailored version of HFS performed far worse than Selby's tailored implementation, to the point where we failed to hit ground states in the largest problems. The failure of this algorithm highlights the extent to which Selby's HFS implementation, and specifically the hardcoded subgraphs to update, exploit the sparsity and modularity of the Chimera topology [4]. It is very likely that this type of exploitation will be impossible in future quantum annealer topologies [3].

b. Wolff cluster Monte Carlo
The Wolff algorithm [60] dynamically detects clusters of spins that should be flipped together. We used a modified implementation that considers the potential change in energy when deciding whether to flip a cluster, similar to Venturelli et al. [61]. This algorithm would, at first glance, be ideal for FCL problems due to the crucial role of clusters. However, finding clusters is slow and our CPU implementation was not competitive with other solvers. Note that the Wolff algorithm is not particularly amenable to GPU implementation; such implementations exist for topologies such as lattices [62] but they achieve only modest speedups over CPU implementations.

c. Parallel tempering
Parallel tempering runs multiple replicas of a Monte Carlo simulation at different temperatures in parallel, and can exchange information between replicas according to certain exchange rules. It is the algorithm of choice for approximately calculating features and statistics of an energy landscape when exact calculation is prohibitive. Our GPU implementation used 64 replicas, with replicas of a spin packed bitwise into a word, similar to our implementation of QMC. Parallel tempering is more powerful than simulated annealing, but in this case proved to be uncompetitive due to the increased cost of each step.

d. PT-ICM
In vanilla parallel tempering, replica exchange steps are combined with single-spin Monte Carlo updates. However, replica exchanges can be combined with other types of updates. Isoenergetic cluster move (ICM) updates can be combined with replica exchange and simple Monte Carlo updates. This variant was introduced by Swendsen and Wang [63] in their original parallel tempering paper and is sometimes called the PT-ICM algorithm. We implemented PT-ICM as described by Zhu et al. [64] and found that its scaling was similar to Selby's implementation of HFS, but absolute performance was an order of magnitude slower even when using a precomputed ideal temperature ladder (i.e., set of β values) specific to each input.
The advantage that PT-ICM has over Selby's implementation of HFS is that PT-ICM detects clusters dynamically and is not tailored to the underlying topology. Because of this, there may be a misconception that PT-ICM will be future proof against denser quantum processor topologies. However, problems on denser topologies will have smaller site-percolation thresholds [65], and the effectiveness of PT-ICM depends crucially on a large sitepercolation threshold [64]. Thus increased processor density will erode the computational value of cluster updates and consequently the advantage of PT-ICM over vanilla PT.

Parameter tuning
For the D-Wave QPU, optimal performance at all problem sizes was achieved at the minimum allowed annealing time of 5 µs. It therefore makes sense to optimize the parameters of software solvers only at the largest problem size; optimizing on a per-size basis would only make scaling look worse for the software solvers. Since we cannot properly optimize the performance of all solvers at all problem sizes, we focus on results at the largest problem size and take scaling results with a grain of salt.
For SA, we chose values of β that increase linearly from 0.01 to 3 and used 10 5 sweeps. For QMC, we used a fixed β of 30 and 10 4 sweeps with the transverse field A(t) decreasing linearly from 1 to 0 and the longitudinal field B(t) increasing linearly from 0 to 1. For SVMC we again used a β of 30 and the same annealing schedule, but used 10 5 sweeps. These values of β were chosen to optimize performance. For HFS, we used Selby's strategy GS-TW2 [9]. For the D-Wave QPU, we use the minimum annealing time of 5 µs.