Quantum versus classical annealing of Ising spin glasses

See allHide authors and affiliations

Science  10 Apr 2015:
Vol. 348, Issue 6231, pp. 215-217
DOI: 10.1126/science.aaa4170

Benchmarking quantum simulation

Finding a solution to a problem often amounts to an optimization problem and thus can be recast in terms of the lowest-energy state of a system. To find such ground states, mathematical methods based on annealing were developed. To reach the ground state more quickly than with the earlier classical methods, a quantum-mechanical approach was proposed; however, the evidence for quantum speed-up is contradictory. Heim et al. show that the results depend on how the problem is described and how the optimization routine is implemented. This development should be valuable for benchmarking quantum machines.

Science, this issue p. 215


Quantum annealers use quantum fluctuations to escape local minima and find low-energy configurations of a physical system. Strong evidence for superiority of quantum annealing (QA) has come from comparing QA implemented through quantum Monte Carlo (QMC) simulations to classical annealing. Motivated by recent experiments, we revisit the question of when quantum speedup may be expected. Although a better scaling is seen for QA in two-dimensional Ising spin glasses, this advantage is due to time discretization artifacts and measurements that are not possible on a physical quantum annealer. Simulations in the physically relevant continuous time limit, on the other hand, do not show superiority. Our results imply that care must be taken when using QMC simulations to assess the potential for quantum speedup.

With the first archaeological records dating back more than 6000 years (1), thermal annealing is likely to be the oldest optimization method in human history. First heating a material and then letting it cool down slowly can relieve internal stresses and allow the material to achieve a lower-energy state. Inspired by thermal annealing, the simulated annealing (SA) algorithm (2) was proposed to find the ground states of combinatorial optimization problems—in particular, Ising spin glasses described by the HamiltonianEmbedded Image (1)Here, the N spins si can take the values ±1. Spins si and sj on lattice sites i and j are coupled by the Ising term Jij, and the hi are local fields. Nonconvex optimization problems—such as finding the ground state of this Ising spin glass (3), job scheduling (4), circuit minimization (5), and chain optimization (6)—find applications in many areas. These problems are all nondeterministic polynomially complete (7), which implies polynomial time mapping from one problem to the other. Thus, any method to efficiently find solutions to the Ising spin glass problem would provide an efficient way of solving other important problems.

Applying the Metropolis algorithm (8), Kirkpatrick et al. demonstrated that using SA—i.e., simulating the process of cooling Ising spin glasses in a Monte Carlo simulation—is an excellent method to minimize Hc (2). Starting in a random state at high temperature, the system is slowly cooled. Thermal excitations allow the escape from local minima and relaxation into a low-energy state with energy equal or near that of the ground state E0 (9). We will refer to the difference between the final energy E and E0 as the residual energy Eres = E − E0. Quantum annealing (QA) (1014) uses quantum tunneling instead of thermal excitations to escape from local minima, which can be advantageous in systems with tall but narrow barriers, which are easier to tunnel through than to thermally climb over. To perform QA of Ising spin glasses, an additional kinetic term is added, usually by applying a transverse magnetic field. The time-dependent Hamiltonian of QA is given byEmbedded Image (2)where Embedded Image and Embedded Image are Pauli z and x operators, respectively. The transverse field Γ(t) is initially much larger then the couplings, Embedded Image, and the spins start out aligned in the x direction. Γ(t) is slowly reduced to zero such that at the end of the annealing process, we recover the Hamiltonian of the initial Ising spin glass problem. On a perfectly coherent quantum device, this algorithm (13) would find the ground state of the spin glass in question with high probability, provided that the annealing time ta is sufficiently long to stay adiabatically in the ground state (15, 16). QA can also be performed at nonzero temperature—for example, on spin glass materials (17) or in programmable devices (18).

Although the simulation of unitary time evolution scales exponentially with N, a variant of QA can be efficiently performed on a classical computer using stochastic dynamics in a quantum Monte Carlo (QMC) simulation (19, 20). In this simulated quantum annealing (SQA) algorithm, the partition function of the quantum Ising model in a transverse field is mapped to that of a classical Ising model in one higher dimension corresponding to the imaginary time direction (21), as shown in Fig. 1. Details of the algorithms are discussed in the supplementary materials (22).

Fig. 1 Illustration of computational resources used in SA and DT-SQA.

(A) A lattice of classical spins used in SA. (B) DT-SQA maps the transverse field Ising model to a classical system representing imaginary time paths of the quantum spins in an additional imaginary time direction. We show an example of M = 4 time slices with discrete time steps Δτ = β/M. (C) With similar computational effort, we can perform SA on M independent replicas.

Strong evidence for QA being superior to classical annealing comes from a comparison of the performance of SQA and SA (14, 19, 20). Upon increasing the annealing time, the residual energy in a two-dimensional (2D) Ising spin glass was seen to drop faster in SQA than in SA, indicating that quantum tunneling may be advantageous in finding low-energy states. However, recent studies of the performance of D-Wave devices (D-Wave Systems, Burnaby, Canada), which are designed to use superconducting circuits to realize QA, failed to demonstrate indications of quantum speedup (23). Furthermore, in contrast to the work in (19, 20), no advantage of SQA over SA was observed.

To understand these contradictory findings, we first confirm the results of Santoro et al. (19) by observing better scaling of SQA compared with SA in Fig. 2 [and in more detail in fig. S1 in (22)]. However, these simulations were performed with a small number of time slices M and a corresponding large time step Δτ = β/M = 1, which we refer to as a discrete time SQA (DT-SQA) simulation. Here, β = 1/kBT (kB, Boltzmann’s constant; T, temperature) is the inverse temperature. Discrete time steps incur time discretization errors of order Embedded Image. To obtain accurate thermal averages for the quantum system, one needs to extrapolate results to Δτ → 0 or perform a continuous time simulation (CT-SQA) that works directly in the limit Δτ → 0 (24).

Fig. 2 Residual energy Eres dependence as a function of annealing time ta.

Graph shows the residual energy for SA, DT-SQA, and CT-SQA for the square-lattice Ising spin glass instance of (19) with 6400 spins. The annealing time is in units of MCS, corresponding to one attempted update per spin. Eres and error bars (indicating 1-σ statistical error) are obtained by averaging over 40 annealing runs.

Repeating the simulations using CT-SQA, we see an entirely different behavior. Although the performance is improved for short ta, the residual energy saturates for longer ta at a level higher than that reached by SA. Whereas the time discretization error in DT-SQA does not affect its use as a classical optimization algorithm, it does not reflect a physical quantum system, for which the continuous time limit is relevant. Hence, the circumstances under which SQA outperforms SA depend on whether we use SQA as a quantum-inspired classical algorithm or as simulation of a physical system.

Understanding the role of time discretization in SQA is important both to estimate the performance of experimental quantum annealers as well as to tune SQA as a classical optimization algorithm. To achieve this, we went beyond the single instance of (19) and studied 1000 random spin glass instances on an 80-by-80 square lattice with periodic boundary conditions. We use the same distribution of uniform couplings Embedded Image and hi = 0. The exact ground-state energies E0 are obtained using the spin glass server (25).

In Fig. 3A, we show Eres as a function of ta for various M. As expected, for M → ∞, DT-SQA converges toward the continuous time limit. We confirm the behavior already indicated in Fig. 2: Although the initial scaling is better in CT-SQA, lower residual energies are reached using a large time step. Comparing M = 16 and 64, a lower residual energy of 2 · 10−3 is found for M = 16 compared with 5 · 10−3 for M = 64, despite the computational effort being four times smaller.

Fig. 3 Convergence and temperature dependence of SQA.

(A) Convergence of DT-SQA toward the continuous time limit M = ∞ obtained by CT-SQA. Shown is the average of 1000 different disorder realizations annealed at an inverse temperature of β = 20. The lowest-energy configuration along the imaginary time axis was taken to calculate the residual energy. (B) Temperature dependence of DT-SQA with a constant number of M = 64 time slices. Lowering the temperature increases the time step Δτ = β/M. This reduces the initial drop in energy but allows us to ultimately find a final configuration with lower energy.

Analyzing the residual energies as a function of temperature with a constant number of time slices in Fig. 3B leads to a similar observation. For β < 20, the DT-SQA results match well with the CT-SQA results [shown in fig. S4 in (22)], indicating that M = 64 is sufficient to converge to the continuous time limit. At lower temperatures, deviations from CT-SQA are seen, and the larger Δτ = β/M allows us to find states with lower energy than in CT-SQA.

When using SQA as a classical optimization algorithm, we can search the final configuration for the time slice with the lowest energy instead of averaging over all time slices. However, we also have to take into consideration the increased computational effort of QMC simulations compared with SA and multiply the number of Monte Carlo steps (MCS) by M for DT-SQA and by β for CT-SQA. Plotting the residual energy as a function of total computational effort in Fig. 4A, we find that with suitable chosen β and M, DT-SQA outperforms SA, agreeing with (19).

Fig. 4 Performance of SQA as a classical optimizer versus a physical system.

(A) SQA as a classical optimizer: When choosing a suitable annealing temperature and time step, DT-SQA scales better than SA, consistent with the results of (19). (B) SQA as a simulation of a physical system, using large enough M to be converged to the continuous time limit. Here we average the final energy over imaginary time instead of picking the lowest-energy configuration.

To use SQA as a classical optimization algorithm, it is thus advantageous to use a small Δτ for short ta, because the continuous time limit has a more rapid initial decrease of Eres. When annealing for longer ta, a lower temperature and larger Δτ are preferred to reach lower asymptotic residual energies. To reach the lowest energies, rather large Δτ values of order unity are needed, where the system consists of few moderately coupled individual replicas instead of a more tightly coupled continuous path of configurations.

Simulated quantum annealing simulations can escape a local minimum through path-integral configurations corresponding to tunneling events. These configurations spend only short (imaginary) times in high-energy configurations on the barrier, thus avoiding the full cost of the barrier. Although the Monte Carlo dynamics in SQA is different from unitary dynamics of QA, they both use a form of quantum tunneling to escape local minima and are thus similarly sensitive to the nature of barriers. Hence, SQA may indicate whether a physical QA device is expected to show an advantage over SA. A good correlation of SQA [and a mean-field version thereof (26)] has been seen for spin glass problems on a D-Wave One device (27).

To use SQA to predict where QA might outperform SA, we have to take the physical continuous time limit and use either CT-SQA or DT-SQA with large enough M. Figure 4B shows the slightly higher residual energy obtained this way. For short ta up to ta = 105 MCS, SQA still outperforms SA when choosing an appropriate temperature, but the asymptotic scaling is better for SA.

In all of our simulations, the residual energy saturates at some point in SQA but not in SA, indicating that the simulations consistently get stuck in local minima for SQA. This may be understood by the more deterministic dynamics of QA, which prefers a subset of low-energy states (28). Repeating SQA, one may thus get consistently stuck in similar local minima. SA, on the other hand, starts in a random state at high temperatures and therefore explores the configuration space more evenly. The more deterministic nature of SQA also explains the counterintuitive result that for some choices of parameters (see Fig. 4), the Eres may increase when annealing more slowly. Similar to the work in (29), perturbing a quantum annealer (for example, by annealing faster) can excite a system out of a local minimum and thus help to ultimately uncover a better solution.

Reinvestigating evidence for QA outperforming classical annealing for spin glass instances, we thus find that the performance advantage previously observed for path-integral QMC annealing compared with SA for 2D spin glasses (19, 20) is due to large imaginary time steps in the path integral and choosing the lowest energy over all time slices. When taking the limit of continuous time and measuring the average energy, the advantage vanishes. These results are also consistent with recent arguments that 2D spin glasses should not see any quantum speedup in QA (30). It will now be important to explore whether 3D spin glasses or spin glasses with long-range couplings, where barriers are taller than in 2D, show superiority for QA. There, and for problem instances derived from hard application problems, it will once more be essential to investigate both discrete and continuous time simulations. CT-QMC simulations can estimate the potential of QA to outperform SA, but ultimately the unitary dynamics of QA will need to help it outperform DT-SQA, which is by itself an efficient classical optimization algorithm. Although the answer to this question will have to wait for experiments on improved QA devices, SQA simulations performed in the correct way can guide the design of these future experiments.

Supplementary Materials

Materials and Methods

Fig. S1 to S4

References (3134)

References and Notes

  1. See the supplementary materials on Science Online.
  2. Acknowledgments: We thank H. G. Katzgraber, G. Santoro, and I. Zintchenko for useful discussions and G. Santoro for providing the spin glass instance used in (19). This work was supported by the Swiss National Science Foundation through the National Competence Center in Research QSIT and by the European Research Council (ERC) through ERC Advanced Grant SIMCOFE. M.T. acknowledges the hospitality of the Aspen Center for Physics, supported by NSF grant 1066293. The spin glass server (25) was used to obtain the ground states for our problem instances. The original data used to create the figures can be obtained from the corresponding author.
View Abstract

Navigate This Article