Efficient Parallel Ising Samplers via Localization Schemes

In a seminal work by [29], it was shown that for ferromagnetic Ising models, approximately sampling from the Gibbs distribution and approximately computing the partition function are both tractable in polynomial time. This result, along with celebrated breakthroughs in polynomial-time approximations of permanents [28, 30] and volumes of convex bodies [15], remarkably showcased the power of random sampling in approximately solving #P-hard inference problems on polynomial-time Turing machines.

Motivated by emerging applications in large-scale data and large models, parallel Gibbs sampling has recently drawn considerable attention [37, 3, 5, 1, 36, 4, 22, 20, 23, 21]. Despite these advances, the tractability of the ferromagnetic Ising model through parallel algorithms remains unresolved, in contrast to the successes of polynomial-time sequential algorithms.

In the theory of computing, the parallel complexity of the ferromagnetic Ising model is of foundational significance. In a seminal work, [39] asked whether an $\mathrm{𝖱𝖭𝖢}$ algorithm (with polylogarithmic depth on polynomial processors) exists for sampling bipartite perfect matchings, which would imply an $\mathrm{𝖱𝖭𝖢}$ approximation algorithm for the permanent of Boolean matrices. To this day, this remains a major open problem. [41] conjectured that no $\mathrm{𝖱𝖭𝖢}$ algorithm exists for this task, making it a rare example of a polynomial-time tractable problem suspected to be intrinsically sequential, yet not known to be P-complete. The challenge of parallelizing the ferromagnetic Ising sampler is closely tied to the problem of sampling matchings, as both problems were resolved sequentially using the same canonical path argument [28, 29].

In this work, we present an efficient parallel sampler in the $\mathrm{𝖱𝖭𝖢}$ class for the general ferromagnetic Ising model with nonzero external fields. To the best of our knowledge, this is the first $\mathrm{𝖱𝖭𝖢}$ sampler for Ising model beyond the critical phase-transition threshold.

One can view this result as a parallel counterpart to [29]. Through a standard reduction via non-adaptive annealing, this $\mathrm{𝖱𝖭𝖢}$ Ising sampler can be turned into an $\mathrm{𝖱𝖭𝖢}$ algorithm for approximating the Ising partition function $Z_{𝜷,𝝀}^{\mathrm{Ising}}$. Specifically, applying the work-efficient parallel annealing algorithm recently developed in [38], the parallel Ising sampler stated in Theorem 1 can be transformed to a randomized parallel algorithm which returns an $(1\pm ϵ)$-approximation of the Ising partition function in ${({ϵ}^{-\frac{1}{\log n}}\cdot \log n)}^{O_{\delta }(1)}$ parallel time on ${\tilde{O}}_{\delta }(m^{2}n/{ϵ}^2)$ machines.

The parallel Ising sampler in Theorem 1 is implied by a parallel sampler for the random cluster model (see Section 2.2 for definition). The correspondence between the Ising model and the random cluster model is due to the well-known Edwards-Sokal coupling [16] (see Lemma 6).

Theorem 1 follows from Theorem 2 via the Edwards-Sokal coupling, as formally described in Section 3. Our main technical contribution in this part is a parallel sampler for the random cluster model, as stated in Theorem 2, whose efficiency is proved in Section 4.

An Ising model can be formulated through the interaction matrix. Let $V$ be a set of $n$ vertices. The Gibbs distribution ${\mu }_{J,h}^{\mathrm{Ising}}$ of the Ising model, supported on ${\{\pm 1\}}^V$, is defined by the interaction matrix $J\in {ℝ}^{V\times V}$, which is a symmetric positive semi-definite matrix, and the external fields $h\in {ℝ}^V$:

Recent studies have established rapid mixing for Ising models in terms of their interaction norm [18, 6], showing that Glauber dynamics mixes rapidly when , whereas for ${\Vert J\Vert }_2>1$ sampling becomes intractable [34, 24, 31]. For parallel sampling, [36] proposed an $\mathrm{𝖱𝖭𝖢}$ algorithm under a stricter condition that requires $J$ to have a constant ${\mathrm{\ell }}_{\mathrm{\infty }}$-norm. However, under the broader condition , the best known parallel algorithms [32, 37] achieve only a depth of $\tilde{O}(\sqrt{n})$, highlighting a significant gap in achieving optimal parallelism.

In this work, we present an $\mathrm{𝖱𝖭𝖢}$ parallel sampler for the Ising model under the condition . To the best of our knowledge, this is the first $\mathrm{𝖱𝖭𝖢}$ Ising sampler that matches this critical threshold for rapid mixing in terms of the interaction norm.

The $\mathrm{𝖱𝖭𝖢}$ sampler stated in Theorem 3 is presented in Section 5.

Classic local Markov chains, such as Glauber dynamics, are inherently sequential, updating the spin of a single site at each step based on its neighbors. Recently, [36, 37] proposed a generic parallelization framework using correlated sampling (or universal coupling), enabling polylogarithmic-depth simulation of Glauber dynamics under a relaxed Dobrushin condition. When applied to the Ising model, this yields an $\mathrm{𝖱𝖭𝖢}$ sampler in the uniqueness regime $𝜷\in {(\frac{\mathrm{\Delta }-2}{\mathrm{\Delta }},\frac{\mathrm{\Delta }}{\mathrm{\Delta }-2})}^E$, where $\mathrm{\Delta }$ is the maximum degree. However, these conditions do not hold for the Ising models considered in this work.

Instead of directly parallelizing a local Markov chain, our approach focuses on simulating “global” Markov chains that update $O(n)$ spins in each step, satisfying the following:

• $\mathrm{■}$

  The chain mixes in polylogarithmic steps;

• $\mathrm{■}$

  Each step of the chain has an efficient parallel implementation in $\mathrm{𝖱𝖭𝖢}$.

These together ensure an $\mathrm{𝖱𝖭𝖢}$ sampler with polylog depth and polynomial total work.

For the ferromagnetic Ising model, this global chain corresponds to the field dynamics introduced by [10], while for the Ising model with a contracting interaction matrix, the global chain is the restricted Gaussian dynamics (also known as the proximal sampler) introduced by [33]. Both of these global Markov chains can be viewed as arising from different localization schemes for the Ising Gibbs sampling.

Proposed by [13], the localization schemes provide an abstract framework that generalizes high-dimensional expander walks and encompasses a wide range of stochastic processes. This framework can be interpreted in terms of noising and denoising processes, as discussed in [17, 9]. Let $X\sim \mu$ be drawn from the target distribution $\mu$. A noisy channel $D$, which is a Markov chain, is applied to $X$ to obtain a “noised” sample $Y\sim D(X,\cdot )$. The joint distribution of $(X,Y)$ is denoted as $\pi$. The denoising process $U$ represents the time reversal of $D$, where it draws $Z$ from the posterior distribution $U(Y,\cdot ):=\pi (\cdot \mid Y)$. The Markov chains associated with the localization scheme update the current state $X$ to a new state $Z$ according to the rules:

• $\mathrm{■}$

  Add noise to $X$ through the noisy channel: sample $Y\sim D(X,\cdot )$;

• $\mathrm{■}$

  Denoise $Y$ via the posterior distribution: sample $Z\sim U(Y,\cdot )$.

In particular, when the noisy channel $D$ corresponds to the continuous-time down walk channel, which drops each element in a set $X$ with a fixed probability, the localization process is called the negative-field localization, and the above associated global Markov chain is the field dynamics. Alternatively, when $D$ corresponds to a Gaussian noise channel, which adds Gaussian noise to the sample $X$, the localization process is called stochastic localization, and the resulting global Markov chain becomes the restricted Gaussian dynamics.

The localization scheme has previously proven remarkably effective in establishing rapid mixing of Glauber dynamics up to criticality [10, 11, 7, 13, 12, 8, 9]. Specifically, the localization process can effectively “tilt” the model parameters towards sub-critical directions. This allows the mixing properties established in sub-critical regimes to be conserved up to criticality, provided the associated global chains mix rapidly.

In the current work, we explore another perspective on localization schemes: leveraging the associated global Markov chains to achieve parallel efficiency in sampling. This task is highly non-trivial, as these chains were originally designed as analytical tools for studying mixing times, and their efficient parallel simulation poses significant challenges.

For the ferromagnetic Ising model with external fields, we leverage the field dynamics of the random-cluster model to design a parallel sampler. The field dynamics is the global Markov chain induced by negative-field localization. Each update consists of two steps:

1. 1.

   A noising step $Y\sim D(X,\cdot )$ which drops each element in $X$ with a fixed probability (as in Line 4 of Algorithm 1), and is straightforward to implement in parallel.

2. 2.

   A denoising step $Z\sim U(Y,\cdot )$, which requires sampling from a tilted posterior distribution $U(Y,\cdot )$ and is non-trivial. We simulate the denoising step via Glauber dynamics on $U(Y,\cdot )$, which corresponds to a sub-critical (low-temperature) random-cluster model, where rapid mixing is easier to establish.

A similar idea was employed in [12] to develop a near-linear time sequential sampler. Here, we further parallelize this Glauber dynamics on $U(Y,\cdot )$ using the parallel simulation algorithm introduced in [37], which yields Algorithm 2. A key challenge, however, is that despite being tilted to the low-temperature regime, the Dobrushin condition required in [37] to ensure efficient parallel simulation does not hold. To overcome this, we establish a new “coupling with the stationary” criterion, which significantly relaxes the Dobrushin condition and guarantees the efficiency of parallel simulation of Glauber dynamics. Altogether, this yields parallel random-cluster and ferromagnetic-Ising samplers with polylogarithmic depth and polynomial total work.

For the Ising model with a contracting interaction matrix $J$, where , we leverage the restricted Gaussian dynamics (also known as the proximal sampler) to design a parallel sampler. This dynamics defines a global Markov chain, which is the associated chain of the stochastic localization for the Ising model. Each update consists of two steps:

1. 1.

   A noising step $Y\sim D(X,\cdot )$, which introduces Gaussian noise to the current sample.

2. 2.

   A denoising step $Z\sim U(Y,\cdot )$, which resamples from a posterior distribution $U(Y,\cdot )$.

With appropriately chosen parameters, $U(Y,\cdot )$ becomes a product distribution, allowing an efficient parallel implementation of the denoising process. The non-trivial step is the noising step $Y\sim D(X,\cdot )$, which samples from a high-dimensional Gaussian distribution with mean $X$ and covariance $J^{-1}$ (as in Line 3 of Algorithm 3). A key observation is that the Gaussian distribution is log-concave, so the noising step can be efficiently approximated using recently developed parallelizations of Langevin Monte Carlo [2].

In a related prior work [5], an efficient parallel algorithm was developed for transport-stable distributions by implementing stochastic localization in discrete steps. However, their approach is fundamentally different from ours. Their method leverages the log-concavity of the tilted (or “noised”) distribution $D(X,\cdot )$, which enables efficient parallel sampling. In contrast, our approach relies on the rapid mixing of the proximal sampler, which also admits an efficient parallel implementation.

Given a finite nonempty ground set $E$, we will use boldface letters, such as $𝒑$, to denote a vector in ${ℝ}^E$. Given a vector $𝒂\in {ℝ}^E$ and function $f:ℝ\to ℝ$, we write $𝒃=f(𝒂)$ for the vector $𝒃\in {ℝ}^E$ satisfying that $b_e=f(a_e)$ for all $e\in E$.

Throughout the paper, we use $\log$ to denote the natural logarithm.

Let $G=(V,E)$ be a connected, undirected graph. The random cluster model on $G=(V,E)$ with edge probabilities $𝒑\in {[0,1]}^E$ and vertex weights $𝝀\in {[0,1]}^V$ is defined as follows. The random cluster distribution ${\mu }_{E,𝒑,𝝀}^{\mathrm{RC}}$ is supported on the power set $2^E$, where each $S\subseteq E$ is assigned a weight:

where $\kappa (V,S)$ is the set of connected components in the graph $(V,S)$. The distribution ${\mu }_{E,𝒑,𝝀}^{\mathrm{RC}}$ is then given by:

where the normalizing factor $Z_{E,𝒑,𝝀}^{\mathrm{RC}}:={\sum }_{S\subseteq E}{w}_{E,𝒑,𝝀}^{\mathrm{RC}}(S)$ is the partition function.

Let $\mathrm{\Omega }$ be a finite state space, and let ${(X_t)}_{t\ge 0}$ be a Markov chain over $\mathrm{\Omega }$ with transition matrix $P$. We use $P$ to represent the Markov chain for convenience. It is well known that a finite Markov chain $P$ converges to a unique stationary distribution $\mu$ if $P$ is irreducible and aperiodic (see [35] for the definitions of these concepts). The mixing time of a Markov chain $P$ is defined by

where $d_{\mathrm{TV}}(𝒑,𝒒):=\frac{1}{2}{\Vert 𝒑-𝒒\Vert }_1$ is the total variation distance.

We assume the concurrent-read exclusive-write ($\mathrm{𝖢𝖱𝖢𝖶}$) parallel random access machine ($\mathrm{𝖯𝖱𝖠𝖬}$) [27] as our model of computation, where concurrent writes to the same memory location are allowed, and an arbitrary value written concurrently is stored. The computational complexity is measured by the number of rounds (parallel time steps) and the number of processors (or machines) used in the worst case.

We use $\mathrm{𝖭𝖢}$ to refer to both the class of algorithms that run in polylogarithmic time using a polynomial number of processors and the class of problems solvable by such algorithms. $\mathrm{𝖱𝖭𝖢}$ denotes the randomized counterpart of $\mathrm{𝖭𝖢}$.