Quantum enhanced Markov chains require fine-tuned quenches

Sampling from a Boltzmann distribution of a classical system is a critical task across multiple fields. Consider a classical Hamiltonian, $H_{c}$, which characterizes a system of $N$ binary random variables $x_{i}=\pm 1$. The Boltzmann distribution,

describes the probability of the system occupying the configuration $x=\{x_{i}\}$ at a specified temperature $T=1/\beta$. As the partition function, $\mathcal{Z}$, is typically intractable, sampling from the Boltzmann distribution is often the only approach for estimating the system’s thermodynamic properties. This class of sampling problem is additionally a common subroutine in machine learning and combinatorial optimization, often emerging as a computational bottleneck in these methods.

Several quantum algorithms for this task, such as quantum walks and quantum simulated annealing, are known to give algorithmic speedups over classical methods for certain problems [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. However, these algorithms rely on fault-tolerant quantum computers with many qubits, making them infeasible on today’s devices. Current quantum devices are limited in size and performance, and few algorithms for sampling from a Boltzmann distribution are available within the constraints of current hardware [15, 16, 17, 18, 19].

One such algorithm, introduced in Ref. [17], the quantum-enhanced Markov chain Monte Carlo, demonstrates promising error resilience. This hybrid approach constructs a Markov chain over classical configurations, proposing new configurations based on measured quantum evolution. Next, a classical step decides whether to accept or reject the proposed configuration as the new state of the chain. Empirical tests reported in Ref. [17] indicated a polynomial speedup in sampling low-temperature spin glass problems compared to uniform and local classical MCMC methods. Additionally, this algorithm was as a subroutine in an adversarial quantum auto-encoder, showing a large temperature region with a larger spectral gap compared to its classical counterparts [20]. In Ref. [21], it was shown that the speedup is not generic across models as it is provably absent for some worst-case problems. In this work, we extend the analysis to a more practically relevant class of Hamiltonians and show generic performance behaviour related to the ergodicity of the quantum quench, allowing for the identification of favourable performance regimes.

Markov chain Monte Carlo (MCMC) is an extremely successful and widely used class of sampling algorithms. These methods rely on Markov chains designed to equilibrate to the desired probability distribution. Often, these chains are constructed with a proposal step followed by an acceptance step, facilitating easier engineering to ensure the conditions for convergence are met. Specifically, if the chain is in the state $y$, a new state $x$ is proposed with probability $Q(x|y)$. In the acceptance step, this new state is adopted as the new state of the chain, with probability $A(x|y)$, known as the acceptance probability. This process is outlined in Fig. 1(a).

The quantum-enhanced Markov chain Monte Carlo algorithm follows this two-step procedure, where only the proposal step requires quantum computation, as seen in Fig. 1(b). The current state of the Markov chain, a classical configuration, is encoded as a computational basis state $\ket{x}$. A unitary $U$ is applied to the state and then the system is measured in the computational basis, resulting in a new classical configuration. This method results in the following proposal probability,

A common acceptance rule, the Metropolis-Hasting probability,

depends only on the ratio of proposal probabilities. If $Q(y|x)=Q(x|y)$, the acceptance probability can be computed efficiently even if the proposal probability is intractable. With this construction, the chain satisfies the detailed balance condition. If $Q(x|y)>0$ $\forall x,y$, the chain is aperiodic and irreducible; ensuring convergence to the desired distribution [17].

Due to the classical acceptance step, this algorithm maintains convergence even with errors in the proposal step, making it well-suited for current noisy quantum devices. The only problematic errors are those that remove the needed proposal probability symmetry. However, mitigation methods such as state preparation and measurement twirling can be employed [17]. Although many choices of $U$ meet the symmetry requirement, they do not all lead to an algorithmic improvement over classical methods. In this work the focus is on evolution of the form $U=e^{-iHt}$ with the following time-independent Hamiltonian,

After initialization, Markov chains require a warm-in period, known as the mixing time, to equilibrate to its stationary distribution. Specifically, the mixing time is the number of steps required for the total variation distance to $\pi$ to be $\epsilon$-small. On difficult problems, chains suffer from slow convergence, making the mixing time the important metric in MCMC performance.

A Markov chain can be specified with a stochastic transition matrix $P$ whose elements $P(y,x)$ describe the probability of moving from the state $y$ to $x$. For the two-step construction,

As $P$ is stochastic, its eigenvalues have a magnitude at most 1, and the unique stationary distribution is associated with the eigenvalue 1. The repeated application of $P$ determines the distribution of the chain, therefore, the mixing time can be connected to the spectral gap $\delta=1-|\lambda_{2}|$ of $P$, formally:

where $\pi_{\min}=\min_{x\in S}\pi(x)$ [22]. The spectral gap is commonly used in performance analyses of MCMC algorithms. Exact spectral gap calculation is, however, limited to small systems as the transition matrix grows with the dimension of the state space and quickly becomes intractable to diagonalize.

The performance of the quantum-enhanced Markov chain Monte Carlo algorithm was probed through spectral gap calculations on the Sherrington-Kirkpatrick (SK) model in Ref. [17]. The two free parameters of the proposal strategy, $t$ and $h$, were randomly drawn at each MCMC step. This strategy was compared to local and uniform classical proposal strategies. Even without optimizing the free parameters, the algorithm exhibited favourable mixing times in the low-temperature regime over the classical counterparts. Although these numerical studies indicated a polynomial speedup, they are limited to small systems. It is unclear whether this advantage persists for larger systems, given the significant finite-sized effects known in the SK model.

The improvement in the mixing time is suggested to be due to the proposal strategy favouring states close in energy, as the quantum evolution depends on the classical Hamiltonian. These states can be difficult to propose classically as they can be decorrelated from the initial state due to the thermalization of the system after the quantum quench. Here, we focus on how the mixing time depends on the free parameters $t$ and $h$ as a strategy to probe the underlying mechanism responsible for the found advantage. Understanding this dependence is essential in knowing if a similar improvement can be obtained where this quantum evolution is classically simulatable. Through analysis of the bottlenecks of the quantum-enhanced method, we show that the speedup observed numerically can only be achieved when quenching to a non-ergodic quantum system, highlighting the difficulty in achieving a generic speedup.

The algorithmic performance of MCMC methods is often not formally derived but determined empirically. This poses a difficulty for the quantum-enhanced method as very limited system sizes are accessible. However, rigorous upper bounds for the mixing time of classical MCMC methods have been found in certain cases using analytical techniques [22]. One such technique is to analyze the bottlenecks of the chain, which restrict the chain’s movement from one section of states to another.

The equilibrium flow between two configurations $x$ and $y$ is defined as $E(x,y)=\pi(x)P(x,y)$. There is a natural graph corresponding to any MCMC method, depicted in Fig. 2, where the nodes represent the possible configurations and the edge weights indicate the equilibrium flow between these configurations. The mixing time of the chain is related to the equilibrium flow between a subset of configurations $S$, to its complement,

An upper bound can be placed on the spectral gap through minimization over possible subsets [22],

where $\pi(S)=\sum_{x\in S}\pi(x)$ is the Boltzmann weight of the subset $S$. Here, the minimization can be interpreted as finding the cut with the most restricted equilibrium flow, thus determining the chain’s mixing time.

At high temperatures, the Boltzmann distribution is close to uniform and, therefore, easy to sample. The difficulties arise at lower temperatures, where the distribution favours low-energy configurations. We will focus on this low-temperature regime, where the quantum-enhanced method numerically showed favourable mixing time. The equilibrium flow of any set $B$ to its complement provides a further upper bound on the spectral gap of this chain,

In the low-temperature regime, for some models, one can devise a reasonable choice of $B$ resulting in a tight bound on the spectral gap. This method allows for the behaviour of the MCMC to be probed for system sizes where exact spectral gap calculation is not accessible numerically [21].

The quantum quench in the quantum-enhanced method can propose configurations far from the previous state in Hamming distance. Ideally, this evolution causes the system to dephase within an energy window around the initial state’s energy, as then the proposed configurations are likely to be accepted. Difficult sampling problems, such as the low-temperature SK model, often have many nearly degenerate states. The energy windows can be exponentially small, suggesting the need for long time evolution. The proposal probability can be expanded in terms of the eigenstates of $H$,

In the long time limit, and in the absence of many degenerate gaps, the proposal probability relaxes to its equilibrium value, equal to the time-averaged probability [23]. Therefore, in the above form, only the terms where $n=m$ contribute,

It can be shown, as presented in Appendix A, that the spectral gap for this time-averaged proposal upper bounds the time-average of the gap. As a result, this limit has shorter mixing times than the average performance of the strategy of randomly selecting $t$ for each MCMC step. Numerical studies of the spectral gap dependence on time, presented in Appendix A, show the best mixing occurs once the system has reached equilibrium. The long-time limit, therefore will be the focus of the analysis. In addition, this regime alleviates the need for fine-tuning the evolution time.

In this limit, the bottleneck upper bound provides valuable insight into the properties of the quantum evolution that lead to favourable mixing. Consider a set $B$ where all configurations in $B$ have higher energy than those in $B^{c}$. An example of such a set partition within the Boltzmann distribution is illustrated in the left and rightmost panels of Fig. 3. As the Boltzmann measure weighs the equilibrium flow, the bound for this type of set $B$ is typically tight in the low-temperature limit. The bottleneck bound takes the following form,

with $\pi_{B}(x)$ the Boltzmann measure normalized over $B$. Two distributions are defined,

the probability of transitioning from the set $B$ to the $n$th quantum eigenstate, weighted by $\pi_{B}$ and

the probability that elements of $B^{c}$ move to the eigenstate $n$. Finally, $\overline{\pi}(B^{c})=|B^{c}|^{-1}\sum_{x\in B^{c}}\pi(x)$ denotes the average Boltzmann probability to be in set $B^{c}$. As Eq. (12) involves the product of these distributions, one achieves fast mixing if there is significant overlap between them. For example, Fig. 3 illustrates these distributions for two quantum proposal strategies that differ by field strength $h$. In the first strategy, shown in Fig. 3(a), the distributions $f(n)$ and $g(n)$ have minimal overlap, resulting in unfavourable mixing compared to in Fig. 3(b), where the overlap is significantly greater.

The purity of the two eigenstate distributions $f$ and $g$ limits the extent of their overlap and consequently, the spectral gap. This result follows from the Cauchy-Schwarz inequality applied to Eq. (12),

These distribution purities can be related to a physical quantity of the quantum evolution, namely the inverse participation ratio (IPR) of classical configurations in the eigenstate basis, defined as,

Specifically, applying Cauchy-Schwarz inequality again, the following bound on the gap is obtained,

The IPR is a measure of the localization of the classical configuration $x$ on the eigenstates of $H$. The minimum IPR corresponds to a completely delocalized state, with $\mathrm{IPR}(x)=1/2^{N}$, whereas a perfectly localized state gives the maximum IPR value of $1$. If the Hamiltonian under which the system is quenched is fully ergodic, the eigenstates are delocalized, and the spectral gap is upper bounded by $1/2^{N}$. Specifically, the proposal probability of Eq. (11) is $Q(x|y)=1/2^{N}$, identical to the strategy of picking configurations uniformly at random. This shows that the quantum quench must be non-ergodic to obtain a scaling advantage over the uniform classical strategy.

The IPR of many-body systems generally scales exponentially with the system size, a consequence of the tensor product structure. For instance, if the states $\ket{n}$ in Eq. (16) are product states where each qubit is rotated away from the classical z-axis by a small angle $\theta$, the IPR would still scale exponentially as ${\rm IPR}\sim\exp(-2\theta^{2}N)$. This immediately highlights the central difficulty in obtaining fast mixing with these quantum-enhanced Markov chains as they suffer from an orthogonality catastrophe, causing all terms within the square root in Eq. (17) to be exponentially small. A large IPR can be attained if the eigenstates of the quantum evolution are well localized on the classical configurations, specifically when the transverse field is perturbatively small. In this perturbative limit, there is no quantum advantage, as the quantum evolution can be replicated classically. Moreover, the likelihood of shifting spectral weight from $B$ to $B^{c}$ decreases with a smaller transverse field. This effect is not reflected in the bound of Eq. (17) but with Eq. (15).

Additionally, consider the following application of Jensen’s inequality,

where $E_{c}=|B_{c}|^{-1}\sum_{x\in B^{c}}H_{c}(x)$ and $F_{\beta}$ is the free energy of the system at inverse temperature $\beta$. Expression (15) can thus be further bound by,

where $F_{c}=E_{c}-T(S_{f}+S_{g})/2$ and $S_{f,g}$ is the Rényi 2-entropy of $f$ and $g$, respectively. This form is reminiscent of classical rate equations, where $F_{c}$ corresponds to the free energy of the transition state. At low temperatures, the Gibbs free energy $F_{\beta}$ approaches $E_{c}$, but the entropy of $f$ and $g$ remains finite, provided the transverse field is finite.

The upper bounds derived here link the mixing time of the quantum-enhanced method to the physical properties of the quantum evolution. Achieving fast mixing requires a balance of competing properties: the distributions $f$ and $g$ must have significant overlap but not too much entropy, as that results in an unfavourably small IPR. Specifically, the low-temperature gap is exponentially small and constrained by the IPR of the low-energy states, allowing regimes of slow mixing to be identified.

As with any MCMC method, the performance of the quantum-enhanced method is problem-dependent. For example, on an unstructured marked item problem, a quantum proposal method has no advantage over a uniform classical proposal [21]. Before investigating more complex spin glass systems, it is insightful to consider a simpler analytically tractable example to elucidate some points made in the previous section.

The performance of the quantum-enhanced MCMC method varies significantly across different systems and parameter regimes. In this work, this behaviour is explained through a bottleneck analysis of the chain, which links MCMC performance to the ergodicity of the quantum Hamiltonian. An upper bound on the spectral gap is established, dependent on the IPR of the classical configurations on the quantum eigenstates. An immediate result of this bound is that the quantum Hamiltonian under which the system is quenched needs to be non-ergodic to have improved scaling over the uniform classical method.

Since this analysis provides an upper bound on the spectral gap, it does not guarantee an improvement in scaling over classical strategies in the non-ergodic regime. Establishing lower bounds through bottleneck analysis is possible but more challenging for this class of problems. However, the bottleneck upper bound offers valuable physical insights into the properties that lead to favourable mixing allowing the identification of parameter regimes where a speedup is expected. We show that there is a generic behaviour across various models, with optimal performance occurring in the long time limit at an intermediate value of $h$. The exact value of this optimal transverse field depends on the specific localization properties of the model and is typically hard to determine. Due to the degeneracy of the quantum Hamiltonian’s transverse field, scaling improvement over the uniform strategy persists into the large $h$ limit. Effective Hamiltonians in this regime can be identified as having similar behaviour but with simpler experimental realization, potentially enabling large system-size demonstrations of this algorithm on today’s quantum devices.

The authors are grateful for ongoing support through the Flatiron Institute, a division of the Simons Foundation. D.S. is partially supported by AFOSR (Award no. FA9550-21-1-0236) and NSF (Award no. OAC-2118310). We acknowledge support from DARPA-STTR award (Award No. 140D0422C0035). The authors acknowledge useful discussions with Pooya Ronagh and Raymond Laflamme.