I Introduction
Quantum transverse Ising models (TIM) have occupied a distinguished role in the study of manybody 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
(1) 
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 offdiagonal matrix elements of are real and nonpositive. The computational basis matrix elements of satisfy the required property after conjugating by the 1local 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 QMAcomplete kempe2004complexity . The special case of frustrationfree stoquastic adiabatic computation can be classically simulated in polynomial time bravyi2009complexity (though this result does not include any nontrivial 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 boundeddegree 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 polynomialtime upper bounds on the run time of Suzuki’s algorithm been obtained for 1D systems with powerlaw 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 polynomialtime approximation scheme (FPRAS) for the partition function. In contrast, for general lowtemperature classical Ising systems with nonferromagnetic interactions there can be no FPRAS for the partition function unless randomized polynomialtime is equal to NP jerrum1993polynomial . In the present work we treat arbitrary nonferromagnetic 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 heatbath 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
(2) 
This result is derived from an analysis of the mixing time of the PIMC Markov chain with worldline updates with heatbath 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 heatbath 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 imaginarytime limit of the quantumtoclassical mapping to obtain a run time that is independent of the Trotter number, and depends only on the expected number of jumps along the imaginarytime 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 quantumtoclassical 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,
(3) 
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 heatbath 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 heatbath 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
(4) 
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 worstcase initial state,
(5) 
and the mixing time of the chain is
(6) 
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
andis 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,(7) 
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 worstcase 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 ^{1}^{1}1Note 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 premetric) one defines a path metric on the entire set by
(8) 
Suppose for each edge there is a coupling of and such that
(9) 
for some then , where , which implies
(10) 
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 ^{2}^{2}2The 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,
(11) 
for any worldline , and so
(12) 
as . Setting
in 10 we see that provided we have , as . The restriction is satisfied when
(13) 
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 singlesite 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 heatbath 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 offdiagonal terms than those in (1), it is natural to ask whether a similar rapid mixing result holds above some systemsize 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 offdiagonal 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 W911NF17C0050. 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.
References

[1]
Subir Sachdev.
Quantum phase transitions.
Handbook of Magnetism and Advanced Magnetic Materials, 2007.  [2] Sei Suzuki, Junichi 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. Twodimensional 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, ZX Gong, and Christopher Monroe. Observation of a manybody dynamical phase transition with a 53qubit 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. MerlinArthur games and stoquastic complexity. arXiv preprint quantph/0611021, 2006.
 [12] Milad Marvian, Daniel A Lidar, and Itay Hen. On the computational complexity of curing nonstoquastic Hamiltonians. Nature communications, 10(1):1–9, 2019.
 [13] Joel Klassen and Barbara M Terhal. Twolocal 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 frustrationfree 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 ddimensional 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. Polynomialtime classical simulation of quantum ferromagnets. Physical review letters, 119(10):100503, 2017.
 [26] Mark Jerrum and Alistair Sinclair. Polynomialtime 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 manybody 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 singlesite 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. Pathintegral 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 continuousspace 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 heatbath 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
where
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 fromand 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
(14) 
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
where
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 bywhere 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 obeyfor some constant . We wish to keep below a constant threshold . This is equivalent to obeying
and so . Each worldline heatbath update as described in the previous section therefore runs in time .
Comments
There are no comments yet.