Classical Simulation of High Temperature Quantum Ising Models

02/06/2020 ∙ by Elizabeth Crosson, et al. ∙ University of New Mexico 0

We consider generalized quantum Ising models, including those which could describe disordered materials or quantum annealers, and we prove that for all temperatures above a system-size independent threshold the path integral Monte Carlo method based on worldline heat-bath updates always mixes to stationarity in time O(n log n) for an n qubit system, and therefore provides a fully polynomial-time approximation scheme for the partition function. This result holds whenever the temperature is greater than four plus twice the maximum interaction degree (valence) over all qubits, measured in units of the local coupling strength. For example, this implies that the classical simulation of the thermal state of a superconducting device modeling a frustrated quantum Ising model with maximum valence of 6 and coupling strengths of 1 GHz is always possible at temperatures above 800 mK. Despite the quantum system being at high temperature, the classical spin system resulting from the quantum-to-classical mapping contains strong couplings which cause the single-site Glauber dynamics to mix slowly, therefore this result depends on the use of worldline updates (which are a form of cluster updates that can be implemented efficiently). This result places definite constraints on the temperatures required for a quantum advantage in analog quantum simulation with various NISQ devices based on equilibrium states of quantum Ising models.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Quantum transverse Ising models (TIM) have occupied a distinguished role in the study of many-body quantum systems sachdev2007quantum ; suzuki2012quantum ; tanaka2017quantum . The 1D TIM has been extensively studied as an exactly solvable model, which exemplifies statistical mechanical dualities by its relation to free spinless fermions and to the 2D classical Ising model RevModPhys.36.856 . In Hamiltonian complexity, ground states of TIM capture the hardness of the complexity class StoqMA that is on the border of quantum and classical complexity bravyi2017complexity ; cubitt2018universal and are universal for a broad class of stoquastic adiabatic computations albash2018adiabatic . Effective Ising interactions are also ubiquitous schauss2018quantum ; zhang2017observation in NISQ era preskill2018quantum devices. A general TIM on -qubits has the form


where the couplings are all real. Here denotes adjacency in the interaction graph which associates qubits with vertices and pairwise Hamiltonian terms with edges.

Every Hamiltonian of the form (1) is stoquastic bravyi2006merlin , which means that there is some choice of local basis in which all of the off-diagonal matrix elements of are real and non-positive. The computational basis matrix elements of satisfy the required property after conjugating by the 1-local unitary (this unitary is said to “cure the sign problem” marvian2019computational ; klassen2019two ), and so we take for each without loss of generality. Approximating the ground energy of a stoquastic local Hamiltonian problem can be done in the complexity class AM bravyi2006merlin while for general local Hamiltonians it is QMA-complete kempe2004complexity . The special case of frustration-free stoquastic adiabatic computation can be classically simulated in polynomial time bravyi2009complexity (though this result does not include any non-trivial transverse Ising models due to frustration caused by the anticommuting nature of Pauli and ), and the random walk used in that result has recently been applied to show that quantum probabilistically checkable proofs based on reductions that preserve the stoquastic property would imply  aharonov2019stoquastic . TIM of the form (1) with polynomially bounded coupling strengths are universal for bounded-degree stoquastic adiabatic computation bravyi2017complexity ; cubitt2018universal .

In 1977 Suzuki introduced a Markov chain Monte Carlo algorithm for approximating the partition function of quantum Ising models and other stoquastic Hamiltonians 

suzuki1977monte (this algorithm motivates special consideration models with restrictions on the signs of Hamiltonian matrix elements). Suzuki’s method, which is now called path integral Monte Carlo (PIMC), is based on relating the quantum partition function of interest to a partition function of a classical spin system suzuki1976relationship , and using a Markov chain Monte Carlo procedure levin2017markov to approximate properties of the latter. In the last decade rigorous polynomial-time upper bounds on the run time of Suzuki’s algorithm been obtained for 1D systems with power-law interactions at constant temperature crosson2018rapid , specific problems related to quantum annealing crosson2016simulated ; jiang2017scaling ; jarret2016adiabatic , and for ferromagnetic systems on arbitrary graphs for temperatures which are at least inverse polynomial small bravyi2015monte ; bravyi2017polynomial . These examples all adopt premises that preclude computational complexity obstructions to finding a fully polynomial-time approximation scheme (FPRAS) for the partition function. In contrast, for general low-temperature classical Ising systems with non-ferromagnetic interactions there can be no FPRAS for the partition function unless randomized polynomial-time is equal to NP jerrum1993polynomial . In the present work we treat arbitrary non-ferromagnetic interactions, but restrict the temperature to be sufficiently high that no complexity obstructions can occur. Algorithms for simulating high temperature quantum systems have been a subject of recent interest harrow2019classical ; kuwahara2019clustering , with those works finding complimentary domains of simulation to the algorithm presented here.

Main result. We establish the existence of a temperature threshold such that the PIMC method is guaranteed to yield an FPRAS for the partition function at any temperature above that threshold. For the same range of temperatures one can also use PIMC to approximately sample from distribution of computational basis measurements () with total variation distance error after heat-bath worldline updates. The runtime bound of

is optimal for a Markov chain method with updates that act locally on the (classical degrees of freedom associated with the) qubits 

hayes2005general .

In terms of the maximum coupling strength and maximum interaction degree these results hold whenever the inverse temperature satisfies


This result is derived from an analysis of the mixing time of the PIMC Markov chain with worldline updates with heat-bath transition probabilities. We show that the mixing time of this Markov chain is

when (2) is satisfied. Therefore the overall sampling algorithm runs in time where is the time it takes to perform a single heat-bath worldline update. In appendix A We analyze a standard implementation of these updates based on the cavity method krzakala2008path to show that . This algorithm appeals to the continuous imaginary-time limit of the quantum-to-classical mapping to obtain a run time that is independent of the Trotter number, and depends only on the expected number of jumps along the imaginary-time direction. In appendix B maximum number of these jumps is determined by a Poisson process with mean .

Ii Preliminaries

Path Integral Monte Carlo.

The PIMC method is based on Suzuki’s quantum-to-classical mapping suzuki1976relationship from a system of qubits to a system of classical spins, which are sometimes described as “replicas” of the original system that are coupled together ferromagnetically. For classical configurations we either write with or where . Individual spins are denoted by where and , and we may write or . In the form the are called “replicas” or “(imaginary) time slices”, while in the form the are called “worldlines.” The goal of the PIMC method is to sample from the following equilibrium distribution on the classical spins,


where is proportional to and

The distribution is the marginal distribution of on a single replica. See bravyi2015monte for a full derivation of the PIMC method.

The distribution 3 is sampled by generalized heat-bath updates dyer2004mixing (i.e. sampling a region of spins from the conditional distribution that fixes spins outside of that region) applied to one worldline at a time, which are called worldline heat-bath updates. If two configurations and differ only in the -th worldline, then the transition probability is

where the factor of is the probability of selecting worldline . Our proof makes use of the form


where the conditional energy function is

with .

Mixing times and path coupling.

Given a Markov chain with stationary distribution , transition matrix , and state space , let be the distribution that results from starting at the initial state and evolving for steps. We measure the distance from stationarity after steps as the total variation distance between and for the worst-case initial state,


and the mixing time of the chain is


The mixing time is an appropriate notion of convergence in this setting because the total variation distance also bounds the difference in expectation values of observables.

A powerful and versatile technique for bounding the mixing time of Markov chains is based on the notion of a coupling. A coupling of two probability distributions


is a pair of random variables

defined on the same probability space with distributed according to and distributed according to . In the analysis of mixing we seek to define a coupling where is distributed according to and is distributed according to , thereby analyzing two trajectories of the Markov chain starting from distinct initial states. If for each we have such a coupling with and , then the time it takes for the two copies of the chain to coincide can be used to upper bound the distance from stationarity,


In other words, the distance to stationarity at time is upper bounded by the probability that the two branches of the coupling have not coincided yet at time , for the worst-case possible pair of starting states. In applications of this method one uses the fact that the two branches of the coupling are defined on the same probability space (which, for the sake of intuition, can be thought of as shared access to random coin flips) to try to update them together as often as possible while still respecting the transition probabilities of each respective branch.

Instead of starting from an arbitrary pair of states , we apply a simplified version of this proof technique called path coupling which was originally proposed by Bubley and Dyer. Here one defines a path metric on pairs of states in and shows that an arbitrary pair of states beginning a distance 1 apart with respect to will come closer together (in expectation) after one step of the Markov chain.

To define a path metric on first consider a connected graph and define for each 111Note that in our usage will be the set of edges corresponding to transitions of , although this is not a requirement in the general method. In addition, one can assign distinct lengths to each edge but this is not required for our usage.. In general a path in from to is a sequence with , , and for each . From this structure (which is sometimes called a pre-metric) one defines a path metric on the entire set by


Suppose for each edge there is a coupling of and such that


for some then , where , which implies


Iii Proof of rapid mixing

Our path coupling applies to the state space graph with vertices and edges given by pairs of configurations that differ only on a single worldline. If we define .

Let the initial states differ at a single worldline . To simulate one step of the Markov chain a worldline is first chosen to be updated uniformly at random. In the case that is not an element of then the conditional distributions of worldline , and , are equal and we can update the worldlines to the same value in the coupling. Otherwise if then due to the influence of worldline and it is not always possible update the chains to the same value. Therefore the expected distance satisfies

By proposition 4.7 in levin2017markov 222The theorem states that for any two distributions , . the probability that they are not updated together in the optimal coupling is

Thus we turn to bounding

First we note that

As the lattice configuration only differs on worldline we have that

from which the bounds

follow. As a direct consequence of this we have

Putting these together implies

since . Therefore

since is normalized. Summarizing,


for any worldline , and so


as . Setting

in 10 we see that provided we have , as . The restriction is satisfied when


From one can obtain the weaker but simpler sufficient expression (2). From this analysis, a tighter but less transparent upper bound on inverse temperatures that suffice for rapid mixing can be obtained by requiring for

This form is particularly useful when the interaction degree is large, but most of the couplings are small, for example in TIM with qubits embedded in a spatial lattice and interaction strengths that decay as a power law with Euclidean distance.

Iv Concluding Remarks

These rapid mixing results depend crucially on the use of worldline updates that avoid the critical slowing down of the single-site Glauber dynamics. In particular, the high temperature regime is said to lead to a “classical freezing” (strong ferromagnetic coupling) along each worldline, and the generalized heat-bath updates avoid this by erasing and resampling an entire worldline at once. This is the reason our result does not follow from known general results on rapid mixing for high temperature classical systems (such systems are known to be rapidly mixing when ).

Since the PIMC method can also be applied to stoquastic Hamiltonians with more general off-diagonal terms than those in (1), it is natural to ask whether a similar rapid mixing result holds above some system-size independent temperature for such models (e.g. those containing terms of the form ). The reason our techniques fall short of addressing this case is that -local off-diagonal terms for create strong interactions between worldlines at all temperatures, implying that is near 1. Therefore an extension of these techniques to such cases would likely require considering more general cluster updates than those applied here (such as “worm algorithm” boninsegni2006worm updates).

V Acknowledgements

The research is based upon work partially supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.


  • [1] Subir Sachdev.

    Quantum phase transitions.

    Handbook of Magnetism and Advanced Magnetic Materials, 2007.
  • [2] Sei Suzuki, Jun-ichi Inoue, and Bikas K Chakrabarti. Quantum Ising phases and transitions in transverse Ising models, volume 862. Springer, 2012.
  • [3] Shu Tanaka, Ryo Tamura, and Bikas K Chakrabarti. Quantum spin glasses, annealing and computation. Cambridge University Press, 2017.
  • [4] T. D. Schultz, D. C. Mattis, and E. H. Lieb. Two-dimensional Ising model as a soluble problem of many fermions. Rev. Mod. Phys., 36:856–871, Jul 1964.
  • [5] Sergey Bravyi and Matthew Hastings. On complexity of the quantum Ising model. Communications in Mathematical Physics, 349(1):1–45, 2017.
  • [6] Toby S Cubitt, Ashley Montanaro, and Stephen Piddock. Universal quantum Hamiltonians. Proceedings of the National Academy of Sciences, 115(38):9497–9502, 2018.
  • [7] Tameem Albash and Daniel A Lidar. Adiabatic quantum computation. Reviews of Modern Physics, 90(1):015002, 2018.
  • [8] Peter Schauss. Quantum simulation of transverse Ising models with rydberg atoms. Quantum Sci. Technol, 3:023001, 2018.
  • [9] Jiehang Zhang, Guido Pagano, Paul W Hess, Antonis Kyprianidis, Patrick Becker, Harvey Kaplan, Alexey V Gorshkov, Z-X Gong, and Christopher Monroe. Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator. Nature, 551(7682):601, 2017.
  • [10] John Preskill. Quantum computing in the nisq era and beyond. Quantum, 2:79, 2018.
  • [11] Sergey Bravyi, Arvid J Bessen, and Barbara M Terhal. Merlin-Arthur games and stoquastic complexity. arXiv preprint quant-ph/0611021, 2006.
  • [12] Milad Marvian, Daniel A Lidar, and Itay Hen. On the computational complexity of curing non-stoquastic Hamiltonians. Nature communications, 10(1):1–9, 2019.
  • [13] Joel Klassen and Barbara M Terhal. Two-local qubit Hamiltonians: when are they stoquastic? Quantum, 3:139, 2019.
  • [14] Julia Kempe, Alexei Kitaev, and Oded Regev. The complexity of the local hamiltonian problem. In International Conference on Foundations of Software Technology and Theoretical Computer Science, pages 372–383. Springer, 2004.
  • [15] Sergey Bravyi and Barbara Terhal. Complexity of stoquastic frustration-free Hamiltonians. Siam journal on computing, 39(4):1462–1485, 2009.
  • [16] Dorit Aharonov and Alex Bredariol Grilo. Stoquastic PCP vs. randomness. In 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 1000–1023. IEEE, 2019.
  • [17] Masuo Suzuki, Seiji Miyashita, and Akira Kuroda. Monte Carlo simulation of quantum spin systems. i. Progress of Theoretical Physics, 58(5):1377–1387, 1977.
  • [18] Masuo Suzuki. Relationship between d-dimensional quantal spin systems and (d+ 1)-dimensional Ising systems: Equivalence, critical exponents and systematic approximants of the partition function and spin correlations. Progress of theoretical physics, 56(5):1454–1469, 1976.
  • [19] David A Levin and Yuval Peres. Markov chains and mixing times, volume 107. American Mathematical Soc., 2017.
  • [20] Elizabeth Crosson and Aram W Harrow. Rapid mixing of path integral Monte Carlo for 1D stoquastic Hamiltonians. arXiv preprint arXiv:1812.02144, 2018.
  • [21] Elizabeth Crosson and Aram W Harrow. Simulated quantum annealing can be exponentially faster than classical simulated annealing. In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 714–723. IEEE, 2016.
  • [22] Zhang Jiang, Vadim N Smelyanskiy, Sergei V Isakov, Sergio Boixo, Guglielmo Mazzola, Matthias Troyer, and Hartmut Neven. Scaling analysis and instantons for thermally assisted tunneling and quantum Monte Carlo simulations. Physical Review A, 95(1):012322, 2017.
  • [23] Michael Jarret, Stephen P Jordan, and Brad Lackey. Adiabatic optimization versus diffusion Monte Carlo methods. Physical Review A, 94(4):042318, 2016.
  • [24] SERGEY BRAVYI. Monte Carlo simulation of stoquastic Hamiltonians. Quantum Information and Computation, 15(13&14):1122–1140, 2015.
  • [25] Sergey Bravyi and David Gosset. Polynomial-time classical simulation of quantum ferromagnets. Physical review letters, 119(10):100503, 2017.
  • [26] Mark Jerrum and Alistair Sinclair. Polynomial-time approximation algorithms for the Ising model. SIAM Journal on computing, 22(5):1087–1116, 1993.
  • [27] Aram Harrow, Saeed Mehraban, and Mehdi Soleimanifar. Classical algorithms, correlation decay, and complex zeros of partition functions of quantum many-body systems. arXiv preprint arXiv:1910.09071, 2019.
  • [28] Tomotaka Kuwahara, Kohtaro Kato, and Fernando GSL Brandão. Clustering of conditional mutual information for quantum Gibbs states above a threshold temperature. arXiv preprint arXiv:1910.09425, 2019.
  • [29] Thomas P Hayes and Alistair Sinclair. A general lower bound for mixing of single-site dynamics on graphs. In 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05), pages 511–520. IEEE, 2005.
  • [30] Florent Krzakala, Alberto Rosso, Guilhem Semerjian, and Francesco Zamponi. Path-integral representation for quantum spin models: Application to the quantum cavity method and Monte Carlo simulations. Physical Review B, 78(13):134428, 2008.
  • [31] Martin Dyer, Alistair Sinclair, Eric Vigoda, and Dror Weitz. Mixing in time and space for lattice spin systems: A combinatorial view. Random Structures & Algorithms, 24(4):461–479, 2004.
  • [32] Massimo Boninsegni, Nikolay Prokof’ev, and Boris Svistunov. Worm algorithm for continuous-space path integral Monte Carlo simulations. Physical review letters, 96(7):070601, 2006.
  • [33] Edward Farhi, David Gosset, Itay Hen, AW Sandvik, Peter Shor, AP Young, and Francesco Zamponi. Performance of the quantum adiabatic algorithm on random instances of two optimization problems on regular hypergraphs. Physical Review A, 86(5):052334, 2012.

Appendix A Worldline Heat Bath Updates

To implement each worldline heat-bath we apply an algorithm introduced in [30] and improved upon in [33] for sampling the conditional distribution (the algorithm is suitable for any Hamiltonian where is diagonal in the basis and where ). For define and let . In order to avoid a dependence on the Trotter number it suffices to store and track a sequence of events at times in which neighboring spins in a worldline flip their value (“jumps”). Notice that the diagonal part of can be written

where and are operator valued functions involving . The influence of the neighboring worldlines on worldline at imaginary time can be deduced via . will be a piecewise constant function, switching values whenever a neighbor of worldline flips its value,

where the times correspond to times where a neighboring spin has flipped.

To generate a new path for worldline the algorithm first computes the value of as a function of imaginary time along the path. Next it generates boundary conditions at the imaginary time points where changes its value (i.e. a neighboring worldline undergoes a flip). This is done by sampling the distribution


is a matrix and .

Now that boundary conditions for regions of contstant spin have been chosen, the algorithm generates subpaths on each interval of length . These subpaths are described by a number of flips and and times

at which flips occur. The number of flips is restricted to being either even or odd depending on the boundary conditions chosen in the previous step. These subpaths are drawn from the probability density on configurations

and finally all subpaths are combined in order to get a full path.

To sample over paths in a region of constant field with fixed boundary conditions, with boundary condition at

, the algorithm draws waiting times from an exponential distribution. Starting with

, the waiting time until the next flip is drawn from

and set , repeating this process repeats until . At this point the path is output if it satisfies the boundary conditions, and otherwise it is discarded and the process is repeated until success.

Appendix B Bounding the Time per Update

The algorithm described in the previous section has a probabilistic run time. In order to upper bound this run time and guarantee an FPRAS for the partition function we introduce a failure condition as follows: if at any point during the PIMC method the number of jumps in any worldline exceeds some value (to be determined below) then terminate and output . We now turn our attention to choosing a suitable . In a given imaginary time region of constant local field the waiting time drawn for the flip is given by


where , switching its value after each flip. (i.e. ). In bounding the probability of a high number of flips we seek the maximum value of as when this is maximized the expected interarrival time between each flip is minimized. It is clear that and , so


is the exponential distribution corresponding to the highest rate of spin flips. Physically this corresponds to the neighboring spins being adversarially aligned at each moment in imaginary time. This interarrival process with the fixed maximum rate

corresponds to a Poisson process over the entire imaginary time interval with the same rate defined by

where the random variable is the number of spin flips in the imaginary time interval. Using a Chernoff bound we have:

Thus the probabililty of observing more than spinflips decreases exponentially witk when . For simplicity we take so that . Since each sample of (3) requires

steps of the Markov chain and assuming whatever quantity we are attempting to estimate requires

samples to compute the probability of failure during entire runtime will obey

for some constant . We wish to keep below a constant threshold . This is equivalent to obeying

and so . Each worldline heat-bath update as described in the previous section therefore runs in time .