Theory of Quantum Annealing of an Ising Spin Glass

See allHide authors and affiliations

Science  29 Mar 2002:
Vol. 295, Issue 5564, pp. 2427-2430
DOI: 10.1126/science.1068774


Probing the lowest energy configuration of a complex system by quantum annealing was recently found to be more effective than its classical, thermal counterpart. By comparing classical and quantum Monte Carlo annealing protocols on the two-dimensional random Ising model (a prototype spin glass), we confirm the superiority of quantum annealing relative to classical annealing. We also propose a theory of quantum annealing based on a cascade of Landau-Zener tunneling events. For both classical and quantum annealing, the residual energy after annealing is inversely proportional to a power of the logarithm of the annealing time, but the quantum case has a larger power that makes it faster.

The annealing of disordered and complex systems toward their optimal (or lowest energy) state is a central problem in statistical physics. The unknown ground state of a system can be approximated by slow-rate cooling of a real or fictitious temperature: The slower the cooling, the closer the approximation gets to the optimal solution (1, 2). Although this kind of standard classical annealing (CA) has been extensively investigated over the last two decades (1–3) and is routinely used in a variety of technological applications, such as chip circuitry design, any improved optimization algorithm would certainly be of enormous value.

Recent results on the spin ½ disordered Ising ferromagnet LiHo0.44Y0.56F4(4, 5) suggested that quantum annealing (QA) works better than CA. In QA, temperature is replaced by a quantum mechanical kinetic term; in the specific case, this term is a transverse magnetic field, Γ, mixing the up and down spin states at each site. Initially the quantum perturbation starts out so large in magnitude as to completely disorder the system even at zero temperature. When the transverse field is subsequently reduced to zero at some slow rate 1/τ, the system is “annealed” toward its ground state, much in the same way as when its temperature is reduced to zero in CA. The question is which of the two, CA or QA, works better, and how and why. Experimental comparison (4) of the properties displayed by the same system transported through two different routes in the [T, Γ] plane from the same initial state A (a classical high-Tstate) to the same nominal final state B (a low-T glassy state) indicates that QA, the “quantum route” from A to B, yields (with the same “cooling” rate) a state B apparently closer to the ground state than that yielded by CA, the classical route. The data, however, do not clarify why that should be so.

Theoretical suggestions and exemplifications of QA made by various groups over the past decade (6–10) have stimulated considerable interest in understanding the mechanism of QA better. A theoretical discussion of the relative merits of CA and QA is therefore desirable. For this, it is necessary to carry out a direct comparative test on a sufficiently representative benchmark system, such as a spin glass, and to lay the bases of a theory of the processes underlying QA. Even in the context of CA, open issues remain. For instance, the way the final annealed energy E final(τ) approaches the ground-state energy E GS as a function of the annealing rate 1/τ is still controversial. Whereas general theoretical arguments by Huse and Fisher (3) predict a slow logarithmic decrease of the residual energyɛ res(τ) =E final(τ) − E GS≈ log−ζ(τ), with ζ ≤ 2, early simulations (11) and more recent studies (12, 13) favor a different form such as power-law,ɛ res(τ) ≈ τ−α, or stretched exponential. The question remains whether the discrepancy between simulations and theory is real, or only apparent.

At the outset, inspired by Brooke et al.'s experimental system (4, 5), we selected the two-dimensional (2D) random Ising model as an appropriate realistic test case. This choice is dictated by the fact the 2D random Ising model is technically a polynomial problem (14)—whereE GS can be calculated up to sufficiently large lattice sizes (15) [thus avoiding an extra fitting parameter inɛ res(τ)]—which is nonetheless of prohibitively large complexity for any physical dynamics as a true glass (3, 16). On this system we carried out CA and QA, and found that QA is indeed faster. For CA we find thatɛ res(τ) deviates from a power law at large τ and is compatible with the Huse-Fisher lawɛ res(τ) ≈ log−ζ(τ). We then built a theory of QA of a spin glass based on the idea of a cascade of level crossings. Our theory suggests an asymptotic decayɛ res(τ) ≈ log−ζQA(τ) that is again logarithmic, as in CA, but is governed by somewhat different exponents that make it faster.

The Edwards-Anderson Hamiltonian of an Ising spin glass in transverse field isEmbedded Image(1)where nearest-neighbor spins 〈ij〉 of a cubic lattice in D dimensions interact with a random exchange couplingJij , Γ is the transverse field inducing transitions between the two states, ↑ and ↓, of each spin, and σi x, σi z are Pauli matrices at site i. The problem is to anneal this system as closely as possible to its classical (Γ = 0) ground state. In CA (1, 2), there is no quantum mechanics (Γ = 0): One starts with a sufficiently high temperatureT 0, which is then reduced linearly to zero in a time τ. In QA, T is instead fixed to zero or some small value, and one starts with a transverse field Γ0sufficiently large to throw the system into a disordered quantum paramagnetic state, decreasing Γ linearly to zero, again in a time τ.

Because real-time annealing is computationally out of the question for the large systems addressed here, we work with the fictitious “time” represented by the number of Monte Carlo (MC) steps. For our implementation of CA, we used a standard Metropolis MC; for QA we used a path integral MC (PIMC) (17, 18) scheme for a quantum system at a small finite temperature T. The 2D quantum Ising model is first mapped on a (2 + 1)D classical model consisting ofP copies (Trotter replicas) of the original lattice, with a nearest-neighbor uniform ferromagnetic coupling in the third (Trotter) direction J = −(PT/2) log tanh(Γ/PT), at temperature PT (17,18). When Γ is large, the replicas are only weakly coupled. As Γ decreases to zero, J increases, eventually forcing all replicas into the same configuration. For a given 2D lattice size L × L (Lup to 80) and for a fixed realization of theJij , drawn from a flat distribution in the interval (−2, 2), we got the exact classicalE GS (15) and then carried out repeated annealings (45 for L = 80), both CA and QA. The annealing parameters T (for CA) or Γ (for QA) were decreased stepwise from the initial value ofT 0 = 3 or Γ0 = 2.5 down to zero, with a total of τ MC steps per spin (18).

The residual energy ɛ res(τ) forL = 80 is plotted in Fig. 1 against τ, actually the MC computer time, for both CA and QA. QA appears definitely superior to CA, with a lower ɛ res(τ) for large τ. This theoretical finding is consistent with the experimental evidence of significantly faster relaxation times observed in QA (4). The τ dependence of our QA data does depend on the chosen values of P and T, particularly on the value of PT, whose optimal value appears to be aroundPT = 1. An increase of P for a fixed value of PT (Fig. 1, inset) ceases to be effective beyond a certain characteristic length (depending on PT) in the Trotter direction. The computational cost increases linearly withP, and the choice P = 20 (corresponding toT = 0.05) was found to be optimal up to the largest values of τ used.

Figure 1

Comparison of the residual energy per site for an 80 × 80 disordered 2D Ising model after classical and quantum annealing. The QA data shown correspond to the optimal value ofPT = 1, with T = 0.05 andP = 20 Trotter replicas. For fair comparison, the annealing time τ used in the QA has been rescaled (multiplied byP) so that points at the same τ require the same computer time (MCS, Monte Carlo steps ). The lower residual energy signifies that QA is superior to CA. Inset: τ-unrescaled QA data for the same system for increasing values of P. Note the satisfactory convergence for P = 20.

A feature evident in our ɛ res(τ) CA data is the small but consistent deviation from a pure power law, suggesting serious reconsideration of all the earlier power-law claims (11, 12). Because the slope (or apparent power) systematically declines for increasing τ, it is natural to ask whether it will asymptotically extrapolate to zero in accordance with the Huse-Fisher logarithmic law (3). If we write this law in the form ɛ res −1/ζ =A log(γτ), where A is a prefactor and γ is a rate constant, and replace the time with the number of MC steps, we can plot the CA data as in Fig. 2. The extrapolated behavior is indeed compatible with a straight line. However, as Fig. 2 shows, it is impossible to extract a value for the exponent ζ; in particular, we cannot establish whether ζ ≤ 2 (3) is any better, as might have been expected.

Figure 2

The same CA data as in Fig. 1, replotted (see text) so as to fall on a straight line if obeying the Huse-Fisher logarithmic law. Although this form is seen to be asymptotically compatible with the data, extraction of a value for the exponent ζ is impossible.

To shed some light on the actual asymptotic form of residual energy in QA, we begin with a cartoon of the instantaneous energy spectrum of Eq. 1 versus Γ (Fig. 3), suggested by small-systems exact diagonalizations. For sufficiently large initial Γ0 ≫ ∣Jij ∣, the ground state, generally nondegenerate (19), must have a finite excitation gap. Imagine following the Schrödinger evolution of an initial ground-state wave function ∣ΨΓ0(t = 0)〉 while reducing Γ gradually to zero as a function of time (10). The instantaneous gap of our disordered magnet will close as Γ decreases through the quantum phase transition at Γc (20–22). After that, ground-state level crossings begin. The arrows in Fig. 3 point to two crossings [really avoided crossings (19) because the problem possesses no symmetry]. Each instantaneous ground-state crossing is associated with tunneling of the whole system between two valleys, one broad but shallow and the other narrow but deep; such tunneling takes place when kinetic energy diminishes, and represents a major crisis in the otherwise quasi-adiabatic evolution caused by the time-dependent decrease of Γ(t).

Figure 3

Cartoon of the lowest instantaneous eigenvalues of a (finite-size) Ising glass as a function of the transverse field Γ, or of a generic complex system as a function of its zero-point kinetic energy Γ. Note the two avoided crossings of the ground state, marked by arrows and enlarged in the upper insets. Lower inset: Schematic of a Landau-Zener crossing. At each crossing the system will follow adiabatically the ground state only if Γ is reduced sufficiently slowly. The infinite system will exhibit an infinite cascade of crossings as Γ → 0.

For sufficiently slow annealing, each tunneling event can be treated as a Landau-Zener (LZ) problem (23) (Fig. 3, inset). The probability P(τ) that the system, starting in the lower state ∣b〉 at high Γ, will continue nonadiabatically onto the higher branch as Γ is reduced with time is given by P(τ) = exp(–τ/τc), where the characteristic tunneling time τc = (ℏ︀αΓ0)/(2πΔ2),ℏ︀ is Planck's constant divided by 2π, Δ is the tunneling amplitude between the two states ∣a〉 and ∣b〉 (whose splitting at crossing is 2Δ), and α is the relative slope of the two crossing branches as a function of Γ (23). One can estimate Δ ≈ exp[−dab /ξ(Γ)], wheredab is a suitable distance between states ∣a〉 and ∣b〉 [in the Ising case, the number of spins flipped in the tunneling process (21)], and ξ(Γ) is a typical wave function localization length, which must vanish as Γ → 0, ξ(Γ) ≈ Γϕ with some exponent ψ > 0. The tunneling time becomes exponentially large for small Γ, τΓ ≈ exp(2dab ϕ), and an exceedingly small width ≈ Δ of each tunneling event justifies treating the multiple crossings as a cascade of independent LZ events. Once the system fails [with a probability P Γ(τ) = exp(−τ/τΓ)] to follow the ground state at the LZ crossing occurring at Γ, it will eventually attain an average excitation energy E ex(Γ) > 0. If we letZ(Γ)dΓ be the number of LZ crossings taking place between Γ and (Γ + dΓ), the average residual energy can be estimated asEmbedded Image(2)where Γc marks the first level crossing. The large τ behavior of this function is dominated by the Γ → 0 behavior of Z(Γ)E ex(Γ), and of τΓ. If we assume that, for small Γ,Z(Γ)E ex(Γ) ≈ Γω and τΓ ≈ exp(Aϕ), we finally getɛ res(τ) ≈ log−ζQA(τ), with an exponent ζQA = (1 + ω)/ϕ. The exponents ω and ϕ are not obvious. A semiclassical expression (23) for the decay of a wave function inside a barrier suggests ϕ = ½. The average excitation energy attained by missing the ground-state “track” at Γ should scale as Γ2 for small Γ, because all eigenvalues start out as Γ2 for Γ → 0. The total number of LZ crossings occurring from 0 to Γ should not be larger than the number of classical states in the energy window (E GS, E GS + Γ), which is approximately equal to ρ(0)Γ [where ρ(0) is the density of classical states at the ground state energy (16)], so that the density of crossingsZ(Γ → 0) → ρ(0), at most. This yields ω = 2 as our most reasonable estimate.

We conclude that ζQA = (1 + ω)/ϕ can be as large as 6 for a spin glass, and that in any case it is above the classical bound ζ ≤ 2 (3). Hence, QA of the Ising spin glass is predicted to be again logarithmically accurate, not fundamentally different in this respect from CA. We therefore expect that a quantum computation based on QA will not transform a hard nonpolynomial (NP-complete) computational problem into a polynomial one. On the contrary, the above reasoning suggests that a logarithmically slow annealing also applies to the present 2D Ising case, which is not NP-complete (14). The slowing-down effect of the LZ cascade illustrated above is particularly severe in problems, like the Ising spin glass we have considered, where the classical spectrum has a gapless continuum of excitations above the ground state. Satisfiability problems, for which encouraging results were recently presented (10), differ from the Ising spin glass in that they possess a discrete classical spectrum and a finite excitation gap. We observe that in general a gap will cut off the LZ cascade precisely in the dangerous low-Γ region, which may eliminate the logarithmic slowing down of QA. Nonetheless, even in the gapless case, the advantage of QA over CA is far from negligible because of the generally larger exponent ζQA of the logarithm. To get an idea of the order of magnitudes involved, consider the relative increase of annealing time (τ′/τ) needed to improve the accuracy of a certain annealing, say with τ ≈ 106 (in appropriate units), by a factor of 10. In CA (ζ = 2), this would require (τ′/τ) ≈ τ[10(1/ζ) ]−1 ≈ 1013. In QA (ζ = 6), the same result would be accomplished with (τ′/τ) ≈ 102.8, an enormous saving of computer effort. Moreover, the PIMC version of QA is easy to implement on a parallel computer, which provides an extra advantage. Optimization by QA of a vast variety of problems, after a suitable fictitious kinetic term is identified case by case, is an open avenue and stands as a worthy challenge for the future.

  • * To whom correspondence should be addressed. E-mail: tosatti{at}


Stay Connected to Science

Navigate This Article