Temperature scaling law for quantum annealing optimizers

TL;DR

Quantum annealers operating at fixed finite temperatures face fundamental limitations in scalability, requiring temperature to decrease with problem size, following a specific scaling law for effective optimization.

Contribution

We derive a temperature scaling law for quantum annealers, showing the necessity of temperature reduction with problem size for scalable optimization.

Findings

Temperature must decrease at least logarithmically with problem size.

Experimental and simulation results support the derived scaling law.

Finite temperature limits the scalability of quantum annealing optimizers.

Abstract

Physical implementations of quantum annealing unavoidably operate at finite temperatures. We point to a fundamental limitation of fixed finite temperature quantum annealers that prevents them from functioning as competitive scalable optimizers and show that to serve as optimizers annealer temperatures must be appropriately scaled down with problem size. We derive a temperature scaling law dictating that temperature must drop at the very least in a logarithmic manner but also possibly as a power law with problem size. We corroborate our results by experiment and simulations and discuss the implications of these to practical annealers.

Figures

Equations

$$H=\sum _{⟨ij⟩}{J}_{ij}{s}_i{s}_j+\sum _i{h}_i{s}_i,$$

$$⟨e{⟩}_{\beta }=⟨H{⟩}_{\beta }/N,$$

$$c_{\beta }=\partial ⟨e{⟩}_{\beta }/\partial \beta =-N\left[⟨e^2{⟩}_{\beta }-⟨e{⟩}_{\beta }^2\right],$$

$$\beta ={\frac{\partial s}{\partial e}|}_{e=e^{\ast }}.$$

$$p_{\beta }(e)\approx \frac{{\text{e}}^{N{\Psi }_{\beta }(e^{\ast })}}{Z_{\beta }}\text{exp}\left[-\frac{N|{\Psi }_{\beta }^{\prime \prime }(e^{\ast })|}{2}(e-e^{\ast }{)}^2\right].$$

$$⟨e{⟩}_{\beta }=e^{\ast }\text{and}{c}_{\beta }=\frac{-1}{|{\Psi }_{\beta }^{\prime \prime }(e^{\ast })|}.$$

$${\sigma }_{\beta }(H)=N{\sigma }_{\beta }(e)=\sqrt{-N{c}_{\beta }},$$

$$⟨H{⟩}_{\beta }-E_0=-N{\int }_{\beta }^{\infty }\mathrm{d}\beta c_{\beta }.$$

$${\sigma }_{\beta }(H)\sim \sqrt{\frac{N}{{\beta }^{\alpha +2}}}\text{and}⟨H{⟩}_{\beta }-E_0=\frac{N}{{\beta }^{\alpha +1}}.$$

$$p_{\beta }(e)=\sqrt{\frac{N}{2\pi |c_{\beta }|}}{\mathrm{e}}^{-N\frac{[e-⟨e{⟩}_{\beta }{]}^2}{2|c_{\beta }|}}.$$

$$p_{\beta +\delta \beta }(e)=\frac{1}{\mathcal{Z}}{p}_{\beta }(e){\mathrm{e}}^{-N\delta \beta e},\mathcal{Z}={\int }_{-\infty }^{\infty }\mathrm{d}e{p}_{\beta }(e){\mathrm{e}}^{-N\delta \beta e}.$$

$$p_{\beta+\delta\beta}(e)=\frac{1}{\cal Z}\mathrm{exp}\big{[}-\frac{N}{2|c_{\beta}|}[e-\langle e\rangle_{\beta}-c_{\beta}\delta\beta]^{2}\big{]}\ ,$$

$$\mathrm{Gaussian}\mathrm{hypothesis}:⟨e{⟩}_{\beta }=⟨e{⟩}_{\beta =0}+\beta c_{\beta },\frac{\mathrm{d}{c}_{\beta }}{\mathrm{d}\beta }=0.$$

$$\frac{\mathrm{d}c_{\beta}}{\mathrm{d}\beta}=N^{2}\big{[}\langle e^{3}\rangle-3\big{\langle}e^{2}\rangle\langle e\rangle+2\langle e\rangle^{3}]=N^{2}\big{\langle}\,[e-\langle e\rangle]^{3}\,\big{\rangle}\,.$$

$$e=⟨e{⟩}_{\beta }+\sqrt{\frac{-c_{\beta }}{N}}\eta .$$

$$⟨\eta ⟩=0,⟨{\eta }^2⟩=1.$$

$$⟨{\eta }^3⟩=\frac{1}{\sqrt{N}}\frac{1}{[-c_{\beta }{]}^{3/2}}\frac{\mathrm{d}{c}_{\beta }}{\mathrm{d}\beta }.$$

$$H_{\mathrm{QA}}(s)=A(s)H_{\mathrm{D}}+B(s)H$$

$${\beta }_i={\left(\frac{{\beta }_{63}}{{\beta }_0}\right)}^{i/63}{\beta }_0,i=0,1,\dots ,63$$

$$Z=\sum _{k=0}^N\left(\genfrac{}{}{0ex}{}{N}{k}\right)e^{-\beta \left(k-N/2\right)}={\left[2\cosh (\beta /2)\right]}^N$$

$$\mu /N=-\frac{1}{2}\tanh (\beta /2)$$

$$\sigma /\sqrt{N}=\frac{1}{2}\mathrm{sech}(\beta /2)$$

$$\beta =-\ln \left(1-p_0^{-1/n}\right)$$

$$\beta =\ln (N)-\ln (-\ln p_0)+\frac{1}{2N}\ln p_0+\dots$$

$$Z=e^{\beta N}\left[1+(2^N-1)e^{-\beta }\right]$$

$$\mu =-N+1-\frac{e^{\beta }}{2^N-1+e^{\beta }}$$

$$\sigma =e^{\beta /2}\frac{\sqrt{2^N-1}}{2^N-1+e^{\beta }}$$

$$\beta =\ln (2^N-1)-\ln \left(p_0^{-1}-1\right)$$

$$\beta =N\ln 2-\ln \left(p_0^{-1}-1\right)-2^{-N}+\dots$$

$$H(s)=-\frac{1}{2}(1-s)\sum _{i=1}^N{\sigma }_i^x-\frac{1}{2}s\sum _{i=1}^N{\sigma }_i^z$$

Full text

Introduction.— Quantum computing devices are becoming sufficiently large to undertake computational tasks that are infeasible using classical computing Bremner et al. (2017); Gao et al. (2017); Bremner et al. (2016); Morimae et al. (2014); Broome et al. (2013); Spring et al. (2013); Boixo et al. (2016). The theoretical underpinning for whether such tasks exist with physically realizable quantum annealers remains lacking, despite the excitement brought on by recent technological breakthroughs that have made programmable quantum annealing (QA) Finnila et al. (1994); Brooke et al. (1999); Kadowaki and Nishimori (1998); Farhi et al. (2000); Santoro et al. (2002) optimizers consisting of thousands of quantum bits commercially available. Thus far, no examples of practical relevance have been found to indicate a superiority of QA optimization, i.e., to find bit assignments that minimize the energy, or cost, of discrete combinatorial optimization problems, faster than possible classically Young et al. (2008, 2010); Hen and Young (2011); Farhi et al. (2012); Rønnow et al. (2014); Hen et al. (2015); Boixo et al. (2013); Albash et al. (2015). Major ongoing efforts continue to build larger, more densely connected QA devices, in the hope that the capability to embed larger optimization problems would eventually reveal the coveted quantum speedup QEO ; Tolpygo et al. (2015a, b); Jin et al. (2015); King et al. (2017).

Understanding the robustness of QA optimization to errors that reduce the final ground state probability is critical. In this work, we consider perhaps the most optimistic setting where the only source of error is due to nonzero temperature. We analyze the theoretical scaling performance of ideal fixed-temperature quantum annealers for optimization. We show that even in the case where annealers are assumed to thermalize instantly (rather than only in the infinite runtime limit), the energies, or costs, of their output configurations would be computationally trivial to achieve (in a sense that we explain). We further derive a scaling law for QA optimizers and provide corroboration of our analytical findings by experimental results obtained from the commercial D-Wave 2X QA processor Johnson et al. (2010); Berkley et al. (2010); Harris et al. (2010); Bunyk et al. (Aug. 2014); King et al. (2015) as well as numerical simulations (our results equally apply to ideal thermal annealing devices). We discuss the implications of our results for both past benchmarking studies and for the engineering requirements of future QA devices.

Fixed-temperature quantum annealers.— In the adiabatic limit, closed-system quantum annealers are guaranteed to find a ground state of the target cost function, or final Hamiltonian $H$, they are to solve. The adiabatic theorem of quantum mechanics ensures that the overlap of the final state of the system with the ground state manifold of $H$, approaches unity as the duration of the process increases Kato (1950); Jansen et al. (2007). For physical quantum annealers that operate at positive temperatures ($T>0$), there is no equivalent guarantee of reaching the ground state with high probability. For long runtimes, an ideal finite-temperature quantum annealer is expected to sample the Boltzmann distribution of the final Hamiltonian at the annealer temperature Venuti et al. (2015).

In what follows, we argue that even instantly-thermalizing quantum annealers 111The existence of small minimum gaps prior to the end of the anneal suggests that it is extremely unlikely that the Gibbs state before crossing these gaps would have a larger overlap with the ground state manifold of the final Hamiltonian than after. Therefore, measuring the system midway through a quantum annealing process will generically yield lower success probabilities than measurements taking place at the end of it Altshuler et al. (2010); Farhi et al. (2012); Knysh (2016). We give two analytical examples in the Supplemental Information. are severely limited as optimizers due to their finite temperature. For concreteness, we restrict to annealers for which i) the number of couplers scales linearly with the number of qubits $N$ 222This is equivalent to having a bounded degree connectivity graph., ii) the coupling strengths are discretized and are bounded independently of problem size, and iii) the scaling of the free energy with problem size is not pathological, i.e., that our system is not tuned to a critical point. Other than the above standard assumptions, our treatment is general (we discuss the performance of quantum annealers when some of these conditions are lifted later on). For clarity, we consider optimization problems written in terms of a Hamiltonian of the Ising-type

[TABLE]

where $\{s_{i}=\pm 1\}$ are binary Ising spin variables that are to be optimized over, $\{J_{ij},h_{i}\}$ are the coupling strengths between connected spins and external biases, respectively, and $\langle ij\rangle$ denotes the underlying connectivity graph of the model. The discussion that follows however is not restricted to any particular model.

Under the above assumptions, the ground state energies, denoted $E_{0}$, of any given problem class, scale linearly with increasing problem size (i.e., the energy is an extensive property as is generically expected from physical systems) while the classical minimal gap $\Delta=E_{1}-E_{0}$ remains fixed. It follows then 333The analysis is based on the equivalence between the Canonical and the Microcanonical Ensembles of Statistical Mechanics. This equivalence is reviewed in many places, see e.g., Ref. Martín-Mayor (2007). that the thermal expectation values of the intensive energy

[TABLE]

and specific heat

[TABLE]

remain finite as $N\to\infty$ for any fixed inverse-temperature $\beta=1/T$. The intensive energy is discretized in steps of $\Delta/N$, yet its statistical dispersion $\sigma_{\beta}(e)=\sqrt{-c_{\beta}/N}$ is much larger. Treating $e$ as a stochastic variable, for large enough values of $N$ it can be treated as a continuous variable as the ratio of discretization versus dispersion is $\sqrt{-\Delta^{2}/(c_{\beta}N)}$ decaying to zero for large $N$. From the Boltzmann distribution it follows that the probability density of $e$ goes as $p_{\beta}(e)={Z_{\beta}}^{-1}\text{e}^{N(s(e)-\beta e)}\,,$ where $Z_{\beta}=\sum_{n}g_{n}\mathrm{e}^{-\beta E_{n}}$ is the partition function, $g_{n}$ is the degeneracy of the $n$-th level, i.e., the number of microstates with $H(\{s_{i}\})=E_{n}$, satisfying $2^{N}=\sum_{n\geq 0}g_{n}$, and $s(e)$ is the entropy density 444Equivalently, $Ns(e)$ is the logarithm of the number of microstates with intensive energy $e$.. The linear combination $\Psi_{\beta}(e)=s(e)-\beta e$ plays the role of a large-deviations functional for $e$. The most probable value of $e$, which we denote by $e^{\ast}$, is given by the maximum of $\Psi_{\beta}$. Solving $\Psi_{\beta}^{\prime}(e^{*})=0$, we find 555A relation best known as the second law of thermodynamics $\mathrm{d}e=T\mathrm{d}s$.

[TABLE]

Close to $e^{*}$, $\Psi_{\beta}$ can be Taylor-expanded as $\Psi_{\beta}(e)\approx\Psi_{\beta}(e^{*})-\frac{|\Psi_{\beta}^{\prime\prime}(e^{*})|}{2}(e-e^{*})^{2}$, from which it follows that

[TABLE]

The probability density is thus approximately Gaussian in the vicinity of $e^{*}$, although deviations from the Gaussian behavior are crucial 666An energy probability density that is precisely Gaussian implies that the energy density is a linear function of the inverse temperature $\beta$ and hence the specific heat is a constant. We elaborate on this point in the Supplemental Information.. Moreover, in the limit of large $N$, we find

[TABLE]

Therefore, the probability of finding by Boltzmann-sampling any energy $e<e^{*}$ (equivalently, $E<e^{*}N$) is exponentially suppressed in $N$, scaling in fact as $\text{exp}[-N(\,\Psi_{\beta}(e^{*})-\Psi_{\beta}(e)\,)]$. We thus arrive at the conclusion that even ideal fixed temperature quantum annealers that thermalize instantaneously to the Gibbs state of the classical Hamiltonian are exponentially unlikely to find the ground state since $e^{\ast}>e_{0}\equiv E_{0}/N$.

We now corroborate the above derivation by runs on the commercial DW2X quantum annealer Johnson et al. (2010); Berkley et al. (2010); Harris et al. (2010); Bunyk et al. (Aug. 2014). To do so, we first generate random instances of differently sized sub-graphs of the DW2X Chimera connectivity graph Choi (2008, 2011) and run them multiple times on the annealer, recording the obtained energies 777The reader is referred to the Supplemental Information for further details.. Figure 1 depicts typical resultant residual energy ($E-E_{0}$) distributions. As is evident, increasing the problem size $N$ ‘pushes’ the energy distribution farther away from $E_{0}$, as well as broadening the distribution and making it more gaussian-like. In the inset, we measure the departure of $\langle H\rangle_{\beta}$ from $E_{0}$ and the spread of the energies $\sigma_{\beta}(H)$ over 100 ‘planted-solution’ Hen et al. (2015) instances per sub-graph size as a function of problem size $N$ 888Details of these instances as well as similar results obtained for other problem classes are given in the Supplemental Information.. For sufficiently large problem sizes, we find that the scaling of $\langle H-E_{0}\rangle_{\beta}$ is close to linear while $\sigma_{\beta}(H)$ scales slightly faster than $\sqrt{N}$. While the slight deviations from our analytical predictions suggest that the DW2X configurations have not fully reached asymptotic behavior999The D-Wave processors are known to suffer from additional sources of error such as problem specification errors Martin-Mayor and Hen (2015); Zhu et al. (2016) and freeze-out before the end of the anneal Amin (2015b) that prevent thermalization to the programmed problem Hamiltonian., they exhibit a trend that closely matches our assumptions with the agreement getting better with growing problem sizes.

Given the scaling of the mean and standard deviation, we conclude that fixed-temperature quantum annealers will generate energies $e$ with a fixed distance from $e_{0}$, or in terms of extensive energies, configurations obtained from fixed-temperature annealers will have energies concentrated around $E=(1-\epsilon)E_{0}$ for some $\epsilon>0$ and $E_{0}<0$.

One could now ask what the difficulty is for classical algorithms to generate energy values in the above range. This question has been recently answered by the discovery of a polynomial time approximation scheme (PTAS) for spin-glasses defined on a Chimera graph  Saket (2013) (and which can be easily generalized to any locally connected model), where reaching such energies can be done efficiently 101010By no means however, is it meant that PTAS is able to thermally sample from a Boltzmann distribution of the input problem. In fact, it should be clear that PTAS does not.. While the scaling of the PTAS with $\epsilon$ is not favorable, scaling as $c^{1/\epsilon}$ for some constant $c$, in practice there exist algorithms (e.g., parallel tempering that we discuss later on) that are known to scale more favorably than PTAS.

Scaling law for quantum annealing temperatures.— In light of the above, it may seem that quantum annealers are doomed to fail as optimizers as problem sizes increase. We now argue that success may be regained if the temperature of the QA device is appropriately scaled with problem size. Specifically, we address the question of how the inverse-temperature $\beta$ should scale with $N$ such that there is a probability of at least $q$ of finding the ground state.

An estimate for the required scaling can be given as follows. From the above analysis, it should be clear that the probability of finding a ground state at inverse temperature $\beta$ will not decay exponentially with system size only if the ground state falls within the variation of the mean energy, specifically if

[TABLE]

is comparable to

[TABLE]

The third law of thermodynamics dictates that the specific heat $c_{T}\equiv\mathrm{d}\langle e\rangle/\mathrm{d}T$ goes to zero when $T\to 0$. Assuming a scaling of the form $c_{T}\sim T^{\alpha}$, or equivalently, $-c_{\beta}\sim\beta^{-\alpha-2}$, gives

[TABLE]

For a power-law specific heat, it thus follows that the sought scaling is $\beta\sim N^{1/\alpha}$. If on the other hand $c_{\beta}$ vanishes exponentially in $\beta$, the inverse-temperature scaling will be milder, of the form $\beta\sim\log{N}$.

To illustrate the above, we next present an analysis of simulations of randomly generated instances on Chimera lattices (we study several problem classes and architectures, see the Supplemental Information). To study the energy distribution generated by a thermal sampler on these instances, we use parallel tempering (PT) Geyer (1991); Hukushima and Nemoto (1996), a Monte Carlo method whereby multiple copies of the system at different temperatures are simulated 111111Details of our PT implementation can be found in the Supplemental Information.. In Fig. 2, we show an example of how the energy distribution of a planted-solution instance changes with $\beta$. The qualitative behavior is similar to what we observe with increasing problem size, whereby decreasing $\beta$ (increasing the temperature) pushes the energy distribution to larger energies and makes it more gaussian-like.

The behavior of the specific heat $c_{\beta}$ as the inverse-temperature $\beta$ becomes large is shown in Fig. 3. At large sizes, the scaling becomes $c_{\beta}\propto\exp(-\Delta\beta)$ as expected (here, $\Delta=4$ is the gap). Based on our predictions above, this should mean that if for a fixed $q$, the minimum $\beta^{\ast}$ such that $p_{\beta^{\ast}}(E_{0})\geq q$ falls in this exponential regime, then we should observe a scaling $\beta^{\ast}\propto\log N$. Indeed, the inset of Fig. 3, which shows simulation results of $\beta^{\ast}$ versus $N$, exhibits the expected $\log N$ behavior 121212Similar scaling behavior for other classes of Hamiltonians, specifically 3-regular 3-XORSAT instances and random $\pm 1$ instances, is also observed, and we give the results for these instances in the Supplemental InformationSupplemental Information..

While for problem classes with a fixed minimum gap $\Delta$, one may naively expect $c_{\beta}$ to vanish exponentially in general, implying that a logarithmic scaling of $\beta$ will generally be sufficient as our simulations indeed indicate, it is important to note that two-dimensional spin glasses are known to exhibit a crossover between an exponential behavior to a power law Jörg et al. (2006); Thomas et al. (2011); Parisen Toldin et al. (2011); Fernandez et al. (2016). This crossover is characterized by a constant $\theta\approx 1/2$, whereby the discreteness of the gap $\Delta$ is evident only for sizes $N^{\theta/2}\ll\beta$. Beyond $N^{\theta/2}\sim\beta$, the $2d$ system behaves as if the coupling distribution is continuous Thomas et al. (2011); Parisen Toldin et al. (2011) at which point the system can be treated as if with continuous couplings, for which the specific heat $c_{T}$ scales as $T^{\alpha}$ with $\alpha_{c}=2\nu$ Jörg et al. (2006), where $\nu=3.53(7)$ Fernandez et al. (2016). Therefore, for an ideal quantum annealer operating beyond the crossover, a scaling of $\beta\sim N^{1/(2\nu)\approx 0.14}$ is required. We may thus expect the same crossover to appear for instances defined on the Chimera lattice, which is $2d$-like. Interestingly, for the temperature scaling shown in the inset of Fig. 3, a power-law fit $\beta\sim N^{\alpha}$ with $\alpha=0.19\pm 0.05$ is almost indistinguishable from the logarithmic one, with a power that is consistent with the $2d$ prediction.

Suboptimal metrics for optimization problems.— For many classically intractable optimization problems, when formulated as Ising models, it is crucial that solvers find a true minimizing bit assignment rather than low lying excited states. This is especially true for NP-complete/hard problems Lucas (2014) where sub-optimal costs generally correspond to violated constraints that must be satisfied (otherwise the resultant configuration is nonsensical despite its low energy). Nonetheless, it is plausible to assume the existence of problems for which slightly sub-optimal configurations would still be of value Vinci and Lidar (2016). We thus also study the necessary temperature scaling for cases where the target energies obey $E_{\mathrm{T}}\leq E_{0}+\delta(N)$ with $\delta(N)$ scaling sub-linearly with problem size. In the inset of Fig. 3, we plot the required scaling of $\beta$ for $\delta(N)=\mathrm{const}$ and $\delta(N)\propto\sqrt{N}$. In both cases we find that a logarithmic scaling is still essential, albeit with smaller prefactors.

Conclusions and discussion.— We have shown that fixed temperature quantum annealers can only sample ‘easily reachable’ energies in the large problem size limit, thereby posing fundamental limitation on their performance. We derived a temperature scaling law to ensure that quantum annealing optimizers find nontrivial energy values with sub-exponential probabilities. The scaling of the specific heat with temperature controls this scaling: if $\beta$ lies in the regime where the specific heat scales exponentially with $\beta$, then the inverse-temperature of the annealer must scale as $\log N$. However, further considerations are needed because of a possible crossover behavior in the specific heat with temperature and problem size. For Chimera graphs, because of their essentially two-dimensional structure, this may lead to a crossover to power law scaling. Little is known about this crossover in three dimensions or for different architectures, so this concern may not be mitigated by a more complex connectivity graph.

Our results shed important light on benchmarking studies that have found no quantum speedups Boixo et al. (2014); Rønnow et al. (2014); Hen et al. (2015); Martin-Mayor and Hen (2015); Zhu et al. (2016), identifying temperature as a relevant culprit for their unfavorable performance. Our analysis is particularly relevant for both the utility as well as the design of future QA devices that have been argued to sample from thermal or close-to-thermal distributions Amin (2015a), calling their role as optimization devices into question.

One approach to scaling down the temperature with problem size is the (theoretically) equivalent scaling up of the overall energy scale of the Hamiltonian. However, the rescaling of the total Hamiltonian is also known to be challenging and may not represent a convenient approach for a scalable architecture. An alternative approach is to develop quantum error correction techniques to effectively increase the energy scale of the Hamiltonian by coupling multiple qubits to form a single logical qubit Pudenz et al. (2014, 2015); Vinci et al. (2015); Matsuura et al. (2015); Vinci et al. (2016); Matsuura et al. (2017) in conjunction with classical post-processing Chancellor (2017, 2016); Karimi and Rosenberg (2017); Karimi et al. (2017) or to effectively decouple the system from the environment Jordan et al. (2006); Bookatz et al. (2015); Jiang and Rieffel (2017); Marvian and Lidar (2017).

Our results reiterate the need for fault-tolerant error correction for scalable quantum annealing, however they do not preclude the utility of quantum annealing optimizers for large finite size problems, where engineering challenges may be overcome to allow the device to operate effectively at a sufficiently low temperature such that problems of interest of a finite size may be solved even in the absence of fault-tolerance. Our results only indicate that this ‘window of opportunity’ cannot be expected to continue as devices are scaled without further improvements in the device temperature or energy scale.

While our arguments above indicate that fixed-temperature quantum annealers may not be scalable as optimizers, the current study does not pertain to the usage of quantum annealers as samplers  Amin (2015a); Adachi and Henderson (2015); Benedetti et al. (2016), where the objective is to sample from the Boltzmann distribution. The latter objective is known to be very difficult task (it is #P-hard Papadimitriou (1995); Valiant (1979); Dahllöf et al. (2005)) and little is known about when or if quantum annealers can provide an advantage in this regard Campos Venuti et al. (2016).

I Deviations to the Gaussian probability density

In the main text, we indicated that deviations from a Gaussian distribution for the marginal of the classical Boltzmann probability (at inverse temperature $\beta$) for the energy density $p_{\beta}$ is crucial. To see why this is the case, let us consider what happens when the probability density is exactly Gaussian:

[TABLE]

The probability density at any other inverse temperature $\beta+\delta\beta$ can be obtained as Falcioni et al. (1982); Ferrenberg and Swendsen (1989)

[TABLE]

Eq. (11) is fully general, and we are not assuming $\delta\beta$ to be small. Now, let us plug the Gaussian probability (Eq. (10)) into Eq. (11). We find

[TABLE]

with ${\cal Z}=\sqrt{\frac{N}{2\pi|c_{\beta}|}}$. Comparing Eqs. (10) and (12), we are led to the conclusion that the Gaussian probability implies that the energy density is a linear function of $\beta$ and that the specific heat is constant:

[TABLE]

Of course, Eq. (13) is grossly in error, because in the limit of large $\beta$ (i.e. zero temperature) $\langle e\rangle_{\beta}$ should reach the Ground State energy-density [rather than diverge as wrongy implied by Eq. (13)]. In fact, the specific heat is not constant. A straightforward application of the fluctuation-dissipation theorem tells us that

[TABLE]

We can introduce $\eta$, the fluctuating part of the energy (regarded as a stochastic variable):

[TABLE]

In combination with the fluctuation-dissipation theorem, $c_{\beta}=-N\langle[e-\langle e\rangle]^{2}\rangle$, we have

[TABLE]

Furthermore, Eq. (14) implies

[TABLE]

However, if $\eta$ be a normal variable $N(0,1)$ as demanded by Eq. (10), we would have $\langle\eta^{3}\rangle=0$ and not what we have in Eq. (16). Hence, convergence to the main traits of the Gauss distribution law (symmetry under $\eta\leftrightarrow-\eta$, for instance) happens at a rate proportional to $1/\sqrt{N}$.

II The DW2X experimental quantum annealing optimizer

II.1 Description of the processor

The experimental results shown in the main text were taken on a 3rd generation D-Wave processor, the DW2X ‘Washington’ processor, installed at the Information Sciences Institute - University of Southern California (ISI). The processor connectivity is given by a $12\times 12$ grid of unit cells, where each unit cell is composed of 8 qubits with a $K_{4,4}$ bipartite connectivity, forming the ‘Chimera’ graph Choi (2008, 2011) with a total of 1152 qubits. Due to miscalibration, there are only 1098 operational qubits on the ISI machine. This is illustrated in Fig. 4.

The device implements the quantum annealing protocol given by the time-dependent Hamiltonian:

[TABLE]

where $H_{\mathrm{D}}=-\sum_{i}\sigma_{i}^{x}$ is the standard transverse field driver Hamiltonian, $H$ is the Ising Hamiltonian [Eq. (1) of the main text], and $A(s),B(s)$ are the annealing schedules satisfying $A(0)\gg B(0)$, $A(1)\ll B(1)$, and $s\equiv t/t_{f}\in[0,1]$ is the dimensional time annealing parameter. The predicted functional form for these schedules is shown in Fig. 5.

II.2 Details of the experiment and additional results

The randomly generated instances tested on the D-Wave processor were run with 20 random gauges Boixo et al. (2013) with 5000 reads per gauge/cycle for a total of $100,000$ anneals per instance. The annealing time chosen for the runs was the default $20\mu$-sec. We further corroborated the analytical derivations discussed in the main text using experiments on the commercial DW2X processor on randomly generated bi-modal $J_{ij}=\pm 1$ instances. As with the planted-solution instances, we first generate random instances of differently sized sub-graphs of the DW2X Chimera connectivity graph Choi (2008, 2011) and run them multiple times on the annealer, recording the obtained energies. Figure 6 depicts the resultant residual energy ($E-E_{0}$) distributions of a typical instance. As is evident, increasing the problem size $N$ ‘pushes’ the energy distribution farther and farther away from the ground state value, as well as broadening the distribution and making it more gaussian-like. In the inset we measure the departure of $\langle H\rangle_{\beta}$ from $E_{0}$ and the spread of the energies $\sigma_{\beta}(H)$ over 100 random bi-modal instances per sub-graph size as a function of problem size $N$. For sufficiently large problem sizes, we find that the scaling of $\langle H-E_{0}\rangle_{\beta}$ is almost linear while $\sigma_{\beta}(H)$ scales slightly faster than $\sqrt{N}$. The results are slightly worse than the analytical prediction but conform to the general trend.

III Simulation Methods

III.1 Instance generation

For the generation of instances in this work we have chosen one problem class to be that of the ‘planted solution’ type—an idea borrowed from constraint satisfaction (SAT) problems. In this problem class, the planted solution represents a ground-state configuration of the Hamiltonian that minimizes the energy and is known in advance. The Hamiltonian of a planted-solution spin glass is a sum of terms, each of which consists of a small number of connected spins, namely, $H=\sum_{j}H_{j}$Hen et al. (2015). Each term $H_{j}$ is chosen such that one of its ground-states is the planted solution. It follows then that the planted solution is also a ground-state of the total Hamiltonian, and its energy is the ground-state energy of the Hamiltonian. Knowing the ground-state energy in advance circumvents the need to verify the ground-state energy using exact (provable) solvers, which rapidly become too expensive computationally as the number of variables grows. The interested reader will find a more detailed discussion of planted Ising problems in Refs.Hen et al. (2015); King et al. (2015).

For the random $\pm 1$ instances on Chimera, we randomly (with equal probability) assign a value $\pm 1$ to all the edges of the Chimera graph. While the ground state energy for these instances is not known with 100% certainty, we ran the Hamze-Freitas-Selby algorithm (HFS) Hamze and de Freitas (2004); Selby (2014) for a sufficiently long time such that we were confident of having found the ground state for these instances.

For the 3-regular 3-XORSAT instances, for each spin, we randomly pick three other spins to which to couple. All couplings are picked to be antiferromagnetic with strength $1$. Because all terms in the Hamiltonian are of the form $+\sigma_{i}^{z}\sigma_{j}^{z}\sigma_{k}^{z}$, the ground state is simply that all-spins-down state.

III.2 Parallel tempering

For the planted-solution instances, we first ‘warmed-up’ our parallel tempering simulation with $5\times 10^{5}$ (for the smaller sizes) to $2\times 10^{6}$ (for the larger sizes) swaps with 10 Monte Carlo sweeps per swap. The temperature distribution is picked as follows:

[TABLE]

with $\beta_{0}=20$ and $\beta_{63}=0.1$. After the warm-up, we sample the energy after every 50 swaps in order to minimize correlation between the energies. We use a total of $10^{4}$ sample points, from which we extract the energies at different quantiles. In order to ensure that we have reached a thermal or near-thermal distribution, we performed the following check. The $10^{4}$ sample points are divided into three blocks: (a) $5\times 10^{3}$ samples from the last half of the samples; (b) $2.5\times 10^{3}$ samples from the second quarter of the samples; (c) $1.25\times 10^{3}$ samples from the second eighth of the samples. We then calculated the specific heat using the samples from each block separately; if the system has sufficiently thermalized and the samples are sufficiently uncorrelated, we expect to observe no change in the specific heat for the three sets of samples within the error bars. We show the results of this test in Fig. 7, where we indeed observe no significant difference.

IV Results for Planted-solution Instances with a Target Energy

In Fig. 8, we supplement the results presented in the main text with the scaling of $\beta$ when the target energy need not be the ground state, specifically $E_{\mathrm{T}}=E_{0}+\delta$. We consider three cases: (i) a constant about the ground state, $E_{\mathrm{T}}=E_{0}+8$, (ii) a square-root scaling above the ground state, $E_{\mathrm{T}}=E_{0}+\sqrt{N/2}$, and a linear scaling above the ground state $E_{\mathrm{T}}=E_{0}+(4+N/32)$. The specific values were picked so that the three cases would have the same target energy at the smallest size of $N=128$. If we fit all curves with a logarithmic dependence on $N$, we observe a similar scaling for the cases of $\delta=\mathrm{constant}$, and the case of $\delta\propto\sqrt{N}$ still exhibits a logarithmic scaling but with a milder coefficient. For the case of $\delta\propto N$, the required $\beta$ approaches a constant for sufficiently large problem sizes.

V Results for the 3-reg 3XORSAT and random $\pm 1$ Chimera Instances

Here we provide the equivalent plots to Fig. (3) of the main text but for the 3-regular 3-XORSAT (Fig. 10) and random $\pm 1$ instances (Fig.  11). The random $\pm 1$ instances were warmed-up with up to $24\times 10^{6}$ PT swaps depending on their size, while the XORSAT instances were warmed-up for with up to $200\times 10^{6}$ swaps depending on their size. For both, as in the planted-solution case, $10^{4}$ samples were taken with one sample after every 50 PT swaps. We perform the same thermalization test as for the planted-solution instances, and we observe no significant difference for the different blocks of samples (see Fig. 9).

We note that for both of these classes of instances, the $\beta$ values required fall in the regime where the scaling of the specific heat with $\beta$ is not yet exponential. The scaling behavior of $\beta^{\ast}$ is consistent with both a $\log N$ and a $N^{1/\alpha}$ behavior.

VI Scaling laws for temperatures: analytical examples

Let us consider the simple case of non-interacting spins in a global magnetic field. This case is particularly relevant if the initial state of the quantum annealer is prepared as the thermal state of the standard driver Hamiltonian $-\frac{1}{2}\sum_{i=1}^{N}\sigma_{i}^{x}$ with no overall energy scaling. The partition function is given by:

[TABLE]

Note that each energy spectrum has a degeneracy that grows polynomially with $N$. The mean energy is given by:

[TABLE]

and the standard deviation is:

[TABLE]

The ground state probability on a thermal state is then given by $p_{0}=\frac{e^{\beta N/2}}{Z}$, which we can then invert to write the inverse-temperature as:

[TABLE]

If we pick $p_{0}$ to be some small but fixed (independent of system size) number and take the large $N$ limit, we find that

[TABLE]

Therefore, we find that for this simple problem, in order to maintain a constant ground state probability while the system size grows, we must scale the inverse-temperature logarithmically with system size.

A Grover search problem Grover (1997); Roland and Cerf (2002) on the other hand yields the worst case scaling. In this case, we take a single state to have energy $-N$, while the remaining states have energy $-N+1$. The partition function is given by:

[TABLE]

with mean energy:

[TABLE]

and standard deviation

[TABLE]

Unlike our other local example the $\sigma$ does not scale as $\sqrt{N}$. The ground state probability is given by $p_{0}=1/Z$. Inverting this for $\beta$, we find:

[TABLE]

Again, for a fixed and small $p_{0}$, expanding for large $N$, we get:

[TABLE]

Therefore, in this case, $\beta$ must grow linearly with $N$ in order to maintain a constant $p_{0}$. Note of course that the Grover Hamiltonian is highly non-local as it contains $N$-body terms.

VII Final ground state weight in Gibbs distributions

We consider here the weight of the final ground state on the Gibbs distributions along a quantum annealing protocol. Let us consider a system of $N$ decoupled qubits evolving under:

[TABLE]

The probability of the final ground state, the state $|0\rangle^{\otimes N}$, in the instantaneous Gibbs distribution is given by:

[TABLE]

where $|\lambda_{0}(s)\rangle$ and $|\lambda_{1}(s)\rangle$ are the instantaneous ground state and first excited state for the single qubit system with eigenvalues $-\lambda(s)$ and $\lambda(s)$ respectively. Let us define $\Lambda(s)\equiv|\langle 0|\lambda_{0}(s)\rangle|^{2}=1-|\langle 0|\lambda_{1}(s)\rangle|^{2}$, so we can rewrite our expression as:

[TABLE]

We therefore have:

[TABLE]

Note that $\Lambda^{\prime}(s)>0,2\Lambda-1>0\forall s$, and $\lambda^{\prime}(s)<0$ for $s<0.5$. We can therefore ask, is it possible for $\frac{d}{ds}p(s)=0$ for $s<0.5$? Because of the exponential factors, for large $\beta$ the first term will dominate the second term, and this will not occur. Let us therefore consider the small $\beta$ case. If we expand the exponentials, we find

[TABLE]

However, an explicit evaluation of the expression in parenthesis gives 1, i.e.:

[TABLE]

Therefore, even in the high temperature limit, $\frac{d}{ds}p(s)$ remains positive. Numerically, we can confirm that $\frac{d}{ds}p(s)$ remains positive. Therefore, we can conclude that $p(s)$ is monotonically increasing and achieves its maximum value at $s=1$.

For the Grover problem, we take Roland and Cerf (2002)

[TABLE]

where $|\phi\rangle$ is the uniform superposition state and $|m\rangle$ denotes the ‘marked’ state which is the ground state at $s=1$. The spectrum is such that only the instantaneous ground state and first excited state have non-zero weight on the marked state for $s<1$. These two states can be written as:

[TABLE]

with eigenvalues $\frac{1}{2}(1-\Delta(s))$ and $\frac{1}{2}(1+\Delta(s))$ respectively and

[TABLE]

The probability of the final ground state in the instantaneous Gibbs distribution is given by:

[TABLE]

For small $\beta$, one has:

[TABLE]

and it is clear from Eq. (37b) that this expression is positive for all $s$. Numerically, we can confirm that $\frac{d}{ds}p(s)$ remains positive for $s\in[0,1]$. Therefore, we can conclude that $p(s)$ is monotonically increasing and achieves its maximum value at $s=1$.