Quantum-Enhanced MCMC

• Quantum-enhanced MCMC is a hybrid sampling paradigm that uses quantum proposals combined with classical Metropolis-Hastings corrections to sample from complex distributions like the Boltzmann distribution of Ising models.
• It leverages unitary evolution and tailored Hamiltonian dynamics to generate proposals that can improve spectral gaps and convergence rates compared to traditional local updates.
• Applications span optimization, many-body computations, and turbulence modeling, though challenges remain for unstructured sampling and overly delocalized operations.

Searching arXiv for papers on quantum-enhanced Markov chain Monte Carlo and related variants. Quantum-enhanced Markov chain Monte Carlo (QeMCMC, QEMC, or QE-MCMC) denotes a hybrid quantum-classical sampling paradigm in which a quantum device generates proposal moves for a Markov chain and a classical Metropolis-Hastings correction accepts or rejects those moves so that the chain targets a prescribed stationary distribution, most often the Boltzmann distribution of a classical Ising model (Layden et al., 2022). In the canonical construction, the current classical configuration is encoded as a computational-basis state, evolved under a Hamiltonian containing a problem term and a mixing term, measured to obtain a candidate configuration, and then classically corrected. Subsequent work has treated this template as a broader proposal-engineering framework: it has been extended to QAOA-type circuits, quantum annealing, continuous-space Hamiltonian Monte Carlo, coarse-grained subproblem updates, irreversible proposals, neural-surrogate proposals, constrained sampling, and application-specific optimization workflows, while also prompting sharp no-go results for worst-case unstructured instances and overly delocalized quenches (Orfi et al., 2024).

1. Foundational hybrid formulation

The original QeMCMC setting is sampling from the Boltzmann distribution of a classical Ising model with energy

$$E(\boldsymbol{s}) = -\sum_{j>k=1}^n J_{jk}s_js_k - \sum_{j=1}^n h_j s_j,$$

and target distribution

$$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$

As in ordinary Metropolis-Hastings, the transition kernel is decomposed into a proposal distribution $Q(\boldsymbol{s}'\mid \boldsymbol{s})$ and an acceptance probability $A(\boldsymbol{s}'\mid \boldsymbol{s})$, with off-diagonal transition probability

$$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$

QeMCMC changes only the proposal step: from the current basis state $\lvert \boldsymbol{s}\rangle$, a quantum processor applies a unitary $U$, measures in the computational basis, and thereby induces

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$

When the proposal is symmetric,

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$

the Metropolis-Hastings ratio simplifies to the usual Metropolis form

$$A(\boldsymbol{s}'\mid \boldsymbol{s})=\min\!\left(1,e^{-[E(\boldsymbol{s}')-E(\boldsymbol{s})]/T}\right).$$

This preserves detailed balance with respect to $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$0 while delegating proposal construction to the quantum device (Layden et al., 2022).

The standard proposal Hamiltonian is

$$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$1

with $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$2 diagonal in the computational basis and $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$3. The intended mechanism is that coherent evolution explores the energy landscape in superposition and proposes configurations that are simultaneously far in Hamming distance and small in energy difference, a combination that local classical updates rarely achieve. The initial superconducting-hardware implementation used this structure to sample hard low-temperature Ising instances and reported faster convergence than local and uniform classical baselines on the studied problems (Layden et al., 2022).

2. Proposal constructions, reversibility, and nonreversibility

A large part of the literature concerns how to construct quantum proposals that are either exactly reversible or correctable by a tractable classical acceptance rule. One line keeps proposal symmetry explicit. QAOA-MC uses a shallow QAOA-like ansatz

$$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$4

so that $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$5 and therefore $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$6 for arbitrary circuit parameters. This removes the need to evaluate proposal ratios explicitly and makes the acceptance step purely Metropolis. In numerical tests on spin-glass Boltzmann sampling, the fitted spectral-gap exponents were $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$7 for optimized QAOA-MC, $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$8 for uniform update, and $$\mu(\boldsymbol{s})=\frac{1}{\mathcal Z}e^{-E(\boldsymbol{s})/T}.$$9 for random circuits, which the authors interpret as an approximately quadratic improvement in convergence scaling (Nakano et al., 2023).

A second line transfers the same logic to continuous variables. Quantum Dynamical Hamiltonian Monte Carlo (QD-HMC) replaces classical leapfrog integration by coherent evolution under

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$0

and generates proposals with

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$1

The method preserves detailed balance by combining proposal symmetry with an explicit momentum-inversion operation

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$2

In this formulation, QeMCMC becomes a continuous-space analogue of Hamiltonian Monte Carlo in which the proposal-generation bottleneck is moved onto quantum hardware rather than onto a classical symplectic integrator (Lockwood et al., 2024).

A third line replaces gate-model dynamics by annealing dynamics. Quantum Annealing Enhanced Markov-Chain Monte Carlo (QAEMCMC) draws proposal states from the output of a quantum annealer, using the measured distribution as the proposal kernel inside a standard Metropolis-Hastings loop. On Sherrington-Kirkpatrick benchmarks, it reports larger spectral gaps, faster convergence of energy observables, and reduced total variation distance relative to classical local and uniform proposals (Arai et al., 12 Feb 2025).

Not all variants remain reversible. Irreversible quantum-enhanced Monte Carlo (IQEMC) introduces a state-dependent Hamiltonian

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$3

with energy-dependent sweep rate $Q(\boldsymbol{s}'\mid \boldsymbol{s})$4, so that generally

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$5

The construction is designed to preserve the Boltzmann stationary distribution through global balance rather than detailed balance. Its physical justification is Landau-Zener-like: fast driving at high energy promotes broad descents, while slow driving near low-energy states stabilizes local exploration (Cao et al., 22 Jun 2026).

3. Spectral gap, mixing-time theory, and no-go results

The central performance metric in this literature is the absolute spectral gap

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$6

which controls the mixing time through standard inequalities such as

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$7

Because the classical accept/reject layer preserves the target distribution, any quantum advantage in the proposal-based framework must appear through improved scaling of $Q(\boldsymbol{s}'\mid \boldsymbol{s})$8 or closely related conductance bounds. In the worst-case marked-item Hamiltonian

$Q(\boldsymbol{s}'\mid \boldsymbol{s})$9

Orfi and Sels prove that for any proposal generated by a unitary, and more generally by a unital channel, one has

$A(\boldsymbol{s}'\mid \boldsymbol{s})$0

which yields

$A(\boldsymbol{s}'\mid \boldsymbol{s})$1

in the low-temperature regime. For that unstructured sampling problem, unital quantum proposals therefore do not improve the asymptotic scaling over naive classical sampling (Orfi et al., 2024).

A complementary limitation arises from excessive delocalization. In the long-time limit of quench-based proposals,

$A(\boldsymbol{s}'\mid \boldsymbol{s})$2

and the gap can be upper-bounded in terms of the inverse participation ratio

$A(\boldsymbol{s}'\mid \boldsymbol{s})$3

When the quench Hamiltonian is ergodic, $A(\boldsymbol{s}'\mid \boldsymbol{s})$4 and $A(\boldsymbol{s}'\mid \boldsymbol{s})$5, so the proposal becomes equivalent to a uniform classical proposal. For the one-dimensional Ising chain, the optimal transverse field obeys

$A(\boldsymbol{s}'\mid \boldsymbol{s})$6

and the best achievable scaling is

$A(\boldsymbol{s}'\mid \boldsymbol{s})$7

matching local classical updates rather than outperforming them. In Sherrington-Kirkpatrick and 3-spin models, the best performance occurs only in an intermediate regime where the quench is neither perturbatively local nor fully ergodic (Orfi et al., 2024).

Later work shows that proposal quality can nevertheless be improved by controlling the localization-delocalization tradeoff. “Adiabatic dressing” replaces a sudden quench by a ramp-hold-ramp protocol

$A(\boldsymbol{s}'\mid \boldsymbol{s})$8

with time-symmetric ramps. In the periodic Ising chain at fixed $A(\boldsymbol{s}'\mid \boldsymbol{s})$9 and $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$0, the dressed protocol yields optimal $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$1 scaling whereas the undressed quench is exponentially bad. In SK and 3-spin models, the fitted exponents improve from $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$2 to $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$3 and from $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$4 to $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$5, respectively (Hsieh et al., 30 Mar 2026).

Irreversibility provides a different route to enlarging the gap. IQEMC reports larger average spectral gaps than reversible QEMC for all tested system sizes, a $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$6 spectral-gap improvement over QEMC at $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$7, and fitted decay slopes $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$8 for IQEMC, $$P(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}'\mid \boldsymbol{s})A(\boldsymbol{s}'\mid \boldsymbol{s}).$$9 for QEMC, and $\lvert \boldsymbol{s}\rangle$0 for annealing in the scaling law $\lvert \boldsymbol{s}\rangle$1 (Cao et al., 22 Jun 2026).

4. Scaling strategies, surrogates, and constrained proposals

One practical limitation of the original formulation is that the quantum register size scales with the problem size. Coarse Grained Quantum-enhanced Markov Chain Monte Carlo (CGQeMCMC) addresses this by updating only subsets of spins. In the “improved local groups” construction, the neglected environment is absorbed into effective fields

$\lvert \boldsymbol{s}\rangle$2

and the “multiple sample” version partitions $\lvert \boldsymbol{s}\rangle$3 spins into $\lvert \boldsymbol{s}\rangle$4 groups of size $\lvert \boldsymbol{s}\rangle$5, applies separate QeMCMC calls, and accepts or rejects only after all group updates are assembled. At $\lvert \boldsymbol{s}\rangle$6, the fitted gap exponents are $\lvert \boldsymbol{s}\rangle$7 for uniform proposals, $\lvert \boldsymbol{s}\rangle$8 for local group, $\lvert \boldsymbol{s}\rangle$9 for improved local group, and $U$0 for multiple sampling. The reported resource reduction is quadratic: the speedup persists while using only $U$1 simulated qubits, and in the 36-spin magnetization example only 6 simulated qubits are required (Ferguson et al., 2024).

Constrained problems have motivated blockwise and learned-proposal variants. For fixed-Hamming-weight Boltzmann sampling, one divide-and-conquer framework partitions the graph into two block decompositions, generates blockwise QAOA samples with an XY mixer

$U$2

and trains a conditional MADE surrogate

$U$3

to preserve the required block Hamming weight. The global feasible set is

$U$4

On random 3-regular graphs, the reported average speedup factors are about $U$5 over local Kawasaki dynamics and $U$6 over global Kawasaki dynamics; on an MNIST feature-mask problem with $U$7, the classification accuracy after 50 steps is $U$8 for the proposed method versus $U$9 for global Kawasaki (Kawamata et al., 22 Apr 2026).

A related surrogate program trains a generative neural sampler on QAOA output and then uses the learned density as an independent proposal,

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$0

Because the proposal density is explicit, the acceptance step uses

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$1

so symmetry constraints on the quantum circuit are no longer required. In low-temperature spin-glass sampling, this hybrid QAOA-plus-MADE construction reports about a $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$2 improvement in spectral gap over uniform proposals, and similar acceleration even without parameter optimization (Nakano et al., 2 Jun 2025).

5. Application domains

QeMCMC has increasingly been used as an optimization primitive rather than only as an exact sampler. For Maximum Independent Set, one implementation combines quantum proposals from a discretized annealing Hamiltonian with warm-starting, low-energy filtering, and parallel tempering. The workflow keeps only the best few quantum samples, randomly chooses among them, and then applies Metropolis-Hastings. On QOBLIB instances, the method recovered global optima for graphs up to 117 decision variables using 117 qubits on IBM hardware; on the 117-node instance, the reported hardware run on ibm_boston converged after 151 iterations using 5 replicas and 10,000 shots per iteration (Marshall et al., 5 Feb 2026).

In variational many-body computation, quantum-assisted variational Monte Carlo (QA-VMC) adapts the proposal mechanism to sampling the variational distribution

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$3

inside VMC. The proposal is generated by quantum time evolution under a chosen Hamiltonian $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$4, while the target remains the neural-network wave-function distribution. For the 10-site one-dimensional Fermi-Hubbard model at $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$5, the quantum proposal yields a spectral gap about an order of magnitude larger than ExcitationSD; for $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$6 at $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$7, it reduces the maximum absolute error and standard deviation by about a factor of 3 and corresponds to about a 9-fold increase in effective sample size (Chang et al., 28 Feb 2025).

Constrained discrete-gravity sampling is another domain. In causal set theory, QeMCMC has been adapted by adding a transitivity-closure penalty

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$8

so that valid causal sets form the ground space of the constraint Hamiltonian, and by deriving a qubit Hamiltonian for the Benincasa-Dowker action. For uniform sampling, the fitted spectral-gap decay rate is $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=\left|\langle \boldsymbol{s}'|U|\boldsymbol{s}\rangle\right|^2.$$9, versus classical rates around $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$0; for the BD-weighted sampler at $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$1, the reported quantum rate is $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$2 versus classical rates around 4 to 5 (Ferguson et al., 24 Jun 2025).

QE-MCMC has also been used as a generative module in turbulence modeling. In the turbulent boundary-layer study, a parametric quantum circuit constructs the proposal distribution for sampling three-component acceleration vectors from a height-dependent target. The method introduces an effective height-weighted spectral gap

$$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$3

and reports that this quantity significantly exceeds the classical MCMC value at the highest qubit numbers and resolutions for multivariate targets with cross-correlations. The module is reported to work reliably for $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$4 qubits per spatial dimension (Schindler et al., 14 Jun 2026).

Continuous-space inference remains a conceptual extension rather than an application-specific benchmark line. QD-HMC is explicitly motivated by machine learning and probabilistic inference, where the target distribution is known through its negative log density $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$5 and proposal generation may benefit from coherent simulation of continuous-space dynamics, although the reported experiments are noiseless exact statevector simulations rather than hardware demonstrations (Lockwood et al., 2024).

6. Limitations, controversies, and related quantum Monte Carlo lineages

The principal controversy concerns the status of quantum advantage. The original proposal reported numerically favorable scaling on small low-temperature spin-glass instances and experimental improvements on a superconducting processor, including an average quantum enhancement factor in the spectral-gap exponent of about $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$6 relative to the best classical baseline at $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$7 (Layden et al., 2022). Subsequent theory showed, however, that worst-case unstructured low-temperature sampling excludes asymptotic speedup for any unital proposal channel (Orfi et al., 2024), and that over-delocalized or ergodic quenches reduce the proposal to an effectively uniform update with no advantage over classical uniform proposals (Orfi et al., 2024).

Several practically strong variants also relax exact-sampler guarantees. The MIS algorithm biases the proposal distribution by filtering quantum shots toward good solutions, which improves optimization but explicitly breaks the exact reversibility assumptions of standard QeMCMC (Marshall et al., 5 Feb 2026). QD-HMC does not provide a full finite-resource error analysis and currently rests on noiseless statevector evidence (Lockwood et al., 2024). Coarse-grained and surrogate approaches introduce approximation error through reduced Hamiltonians, block decompositions, or learned proposal models, even though they preserve a Metropolis-Hastings correction at the outer level (Ferguson et al., 2024).

The term “quantum-enhanced MCMC” also borders adjacent but distinct quantum Monte Carlo programs. Quantum-walk and amplitude-estimation algorithms seek near-quadratic speedups for expectation estimation or partition-function approximation through coherent state preparation and quantum mean estimation, with bounded-variance complexity improving from $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$8 classically to $$Q(\boldsymbol{s}'\mid \boldsymbol{s})=Q(\boldsymbol{s}\mid \boldsymbol{s}'),$$9 quantumly (Montanaro, 2015). Hardware qMCMC implementations on trapped-ion devices encode the Markov chain itself into a quantum walk and then use phase estimation or quantum amplitude estimation, which is algorithmically different from the proposal-plus-classical-acceptance architecture of QeMCMC (Claudon et al., 9 Mar 2026). The inverse direction also exists: “Markov Chain Monte-Carlo Enhanced Variational Quantum Algorithms” uses classical MCMC to improve a variational quantum optimizer, rather than using quantum proposals to improve a classical Markov chain (Patti et al., 2021).

This suggests that quantum-enhanced MCMC is best understood not as a single algorithm with a uniform complexity story, but as a proposal-engineering framework whose performance depends on proposal symmetry, localization properties, irreversibility, constraint handling, and hardware-resource tradeoffs. Within the literature summarized here, the strongest negative results apply to worst-case unstructured sampling with unital proposals, whereas the strongest positive results arise in structured, constrained, or application-specific settings where proposal design, surrogate modeling, coarse graining, or nonreversible dynamics materially alter the spectral gap or the autocorrelation structure.