Fast Algorithms for Stoquastic Spin Systems

Abstract: We establish a general framework for developing fast sampling and counting algorithms for stoquastic spin systems at high temperature. Our framework is based on a rapidly mixing Markov chain for polymer models and a subcritical percolation process for sampling individual polymers. We apply our framework to obtain fast algorithms for approximating the partition function and sampling from the thermal distribution of (1) general stoquastic spin systems, (2) ferromagnetic Heisenberg models, and (3) antiferromagnetic Heisenberg models on bipartite graphs. For the Heisenberg models, we obtain an improved bound on the inverse temperature by using their respective cycle and loop representations.

• The paper develops a polymer-dynamics and graphlet-sampling framework with unbiased weight estimators, achieving expected sampling time O(|G| log(|G|/ε)) and counting time O(|G|² ε⁻² log²(|G|/ε)).
• The paper improves high-temperature guarantees for general stoquastic systems to β ≤ 1/(e³d²Δ) with probabilistic estimators, while deterministic methods achieve β ≤ 1/(2e³d³Δ).
• The paper obtains dimension-independent thresholds for Heisenberg models—β ≤ 1/(e³Δ) ferromagnetically and β ≤ 1/(2e³Δ) antiferromagnetically on bipartite graphs—while leaving low temperatures, unbounded degree, and non-stoquastic systems open.

This paper develops a general algorithmic framework for approximating partition functions and sampling from thermal distributions of stoquastic spin systems at high temperature, achieving near-linear sampling runtimes. The framework combines two classical-algorithmic ingredients — the polymer dynamics Markov chain of Chen et al. and the graphlet sampling (percolation) technique of Blanca et al. — and adapts them to quantum polymer models arising from stoquastic Hamiltonians, where polymer weights are expressible as non-negative traces over exponentially large Hilbert spaces.

Framework: polymer dynamics with efficient weight sampling

The central abstraction is a graphlet polymer model on a bounded-degree graph $G$, with a non-negative weight function satisfying a decay condition $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$ for all polymers $\gamma$. Under this condition, together with a "vertex incompatibility condition" derived from the standard bound of at most $(e\Delta)^{k-1}$ graphlets of order $k$ containing a fixed vertex, the polymer distribution $\mu_\mathcal{P}$ is well-conditioned.

Two subroutines drive the framework:

• Polymer dynamics: a Glauber-type Markov chain on admissible polymer sets, shown to be ergodic with unique stationary distribution $\mu_\mathcal{P}$ and mixing time $O(|G|\log(|G|/\epsilon))$. Each update requires sampling from a local proposal distribution $\pi_v$.
• Graphlet sampler: a subcritical percolation process (a breadth-first search with edge-inclusion probability $p$ chosen so that $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$0) that produces perfect samples from an auxiliary local distribution. The percolation cluster size is stochastically dominated by a Galton–Watson total progeny with offspring $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$1, $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$2, which is what makes expected runtime constant when weight evaluation is cheap.

A key technical contribution is the treatment of $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$3-computable polymer weights, where computing or sampling a weight costs $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$4. The sampler rejects polymers larger than a threshold $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$5 before evaluating their weights, yielding expected runtime $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$6 when $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$7 and $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$8 otherwise. This cleanly separates the statistical decay parameter $w_\gamma \leq (1/(e^{r+2}\Delta))^{|\gamma|}$9 from the computational cost parameter $\gamma$0: fast algorithms require the weight decay to dominate the evaluation cost.

The paper also introduces probabilistic polymer models, in which each weight is replaced by a non-negative unbiased random estimator bounded by the same decay condition. This is essential for quantum applications, since exact trace evaluations are expensive but unbiased estimators from random operator products are not. The main results are then:

• An $\gamma$1-approximate sampling algorithm for $\gamma$2 in expected time $\gamma$3 when $\gamma$4.
• An $\gamma$5-approximate counting algorithm for $\gamma$6 via simulated annealing, in expected time $\gamma$7 under the same condition.

A truncation argument extends these guarantees to arbitrary (untruncated) polymer models at polynomial overhead when $\gamma$8, using that truncating at $\gamma$9 incurs only $(e\Delta)^{k-1}$0 error in both the partition function and total variation distance.

Application to general stoquastic systems

For a stoquastic spin system on a graph of maximum degree $(e\Delta)^{k-1}$1 with local dimension $(e\Delta)^{k-1}$2, the paper constructs a subgraph polymer model whose weights arise from inclusion–exclusion over connected components of edge subsets. Stoquasticity (all matrix elements of each interaction non-positive in the spin basis) ensures non-negativity of weights and of the conditional spin distributions $(e\Delta)^{k-1}$3, which is what permits probabilistic estimation.

Using the deterministic representation, the algorithms apply for $(e\Delta)^{k-1}$4; using the probabilistic representation based on Poisson-sampled edge sequences with Stirling-number correction factors, this improves to $(e\Delta)^{k-1}$5, since the estimator's $(e\Delta)^{k-1}$6 improves on the deterministic value $(e\Delta)^{k-1}$7. In both cases, sampling runs in $(e\Delta)^{k-1}$8 expected time and counting in $(e\Delta)^{k-1}$9. The improvement factor of roughly one power of $k$0 in the temperature threshold comes entirely from cheaper unbiased estimation rather than sharper combinatorics — a point the paper makes concrete through the explicit tradeoff in the truncated-model corollary, which allows any $k$1 with polynomially degraded runtime $k$2 when $k$3.

Heisenberg models via cycle and loop representations

The strongest quantitative results concern Heisenberg models, where specialized representations remove the dependence on $k$4 entirely, giving $k$5 for the ferromagnetic model and $k$6 for the antiferromagnetic model on bipartite graphs — both with $k$7.

For the ferromagnetic case, the paper uses Toth's cycle representation: transpositions $k$8 generate permutations of the vertex set, and the polymer weight decomposes into factors $k$9 over cycles $\mu_\mathcal{P}$0 of the induced permutation, each bounded via $\mu_\mathcal{P}$1. Sampling a permutation product and its cycle decomposition takes polynomial time in $\mu_\mathcal{P}$2, so weight estimators are computable in essentially linear time ($\mu_\mathcal{P}$3).

For the antiferromagnetic model on bipartite graphs, the Aizenman–Nachtergaele loop representation applies: projectors onto singlet states generate loop configurations on the space-time graph, with internal loops contributing factors of $\mu_\mathcal{P}$4 and external loops contributing $\mu_\mathcal{P}$5 weighted by winding numbers. Bipartiteness is essential here, as it controls the winding-number contribution and yields the same $\mu_\mathcal{P}$6 suppression used in the ferromagnetic bound. Loop decomposition is again polynomial-time computable.

In both Heisenberg cases, the inverse-temperature thresholds are independent of local dimension and match, up to constants, the best known high-temperature thresholds for related classical and quantum polymer methods, while the runtimes are substantially faster than prior cluster-expansion-based quantum algorithms.

Summary of results

Model	Inverse temperature condition	Sampling runtime	Counting runtime
General stoquastic (deterministic)	$\mu_\mathcal{P}$7	$\mu_\mathcal{P}$8	$\mu_\mathcal{P}$9
General stoquastic (estimator)	$\mu_\mathcal{P}$0	$\mu_\mathcal{P}$1	$\mu_\mathcal{P}$2
Ferromagnetic Heisenberg	$\mu_\mathcal{P}$3	$\mu_\mathcal{P}$4	$\mu_\mathcal{P}$5
Antiferromagnetic Heisenberg (bipartite)	$\mu_\mathcal{P}$6	$\mu_\mathcal{P}$7	$\mu_\mathcal{P}$8

All bounds hold for graphs of maximum degree at most $\mu_\mathcal{P}$9, and the sampling guarantee is in total variation distance while the counting guarantee holds with probability at least $O(|G|\log(|G|/\epsilon))$0.

Limitations and open questions

Several restrictions are intrinsic to the present approach. All results are confined to the high-temperature regime: the polymer weight decay condition fails below the stated $O(|G|\log(|G|/\epsilon))$1 thresholds, and nothing is claimed about critical or low-temperature behavior, where the relevant polymer expansions diverge. The framework also assumes bounded degree, inherited from the graphlet-counting lemma and the percolation analysis; extension to unbounded-degree graphs, achieved classically by Blanca et al., remains open here. For antiferromagnetic Heisenberg models, bipartiteness is required for the loop-representation argument, so the non-bipartite case is untreated. Finally, stoquasticity itself is load-bearing: it is what guarantees non-negativity of the inclusion–exclusion weights and hence the validity of the probabilistic estimators, and the paper offers no mechanism for handling signful Hamiltonians.

The conclusion identifies three specific open problems: determining the optimal inverse-temperature threshold for these methods, extending to low temperatures and unbounded degrees as has been done for classical polymer models, and obtaining comparably fast classical algorithms for general (non-stoquastic) quantum spin systems.

Conclusion

This work transfers the modern toolkit of fast classical polymer algorithms — rapidly mixing polymer dynamics combined with subcritical percolation sampling — to quantum stoquastic spin systems. The conceptual device that makes the transfer possible is the separation of statistical decay ($O(|G|\log(|G|/\epsilon))$2) from computational cost ($O(|G|\log(|G|/\epsilon))$3), together with unbiased, bounded weight estimators built from Poisson-sampled operator sequences. The resulting algorithms achieve $O(|G|\log(|G|/\epsilon))$4 sampling and quasi-quadratic counting runtimes, with dimension-free temperature thresholds for Heisenberg models obtained through the cycle and loop representations. The framework leaves open the natural extensions to lower temperatures, unbounded degree, and signful interactions.