Discrete Diffusion Models via Evolving Variational Autoregressive Networks

# Discrete Diffusion Models via Evolving Variational Autoregressive Networks
Source: [https://arxiv.org/html/2609.27306](https://arxiv.org/html/2609.27306)
Kewen PanAffiliation:Institute of Fundamental and Frontier Sciences, University of Electronic Science and Technology of China, Chengdu 611731, ChinaYing TangCorresponding authors:[jamestang23@gmail\.com](mailto:[email protected])Affiliation:Institute of Fundamental and Frontier Sciences, University of Electronic Science and Technology of China, Chengdu 611731, ChinaAffiliation:School of Physics, University of Electronic Science and Technology of China, Chengdu 611731, ChinaAffiliation:Key Laboratory of Quantum Physics and Photonic Quantum Information, Ministry of Education, University of Electronic Science and Technology of China, Chengdu 611731, ChinaAffiliation:Non\-classical Information Science Basic Discipline Research Center of Sichuan Province, University of Electronic Science and Technology of China, Chengdu 611731, China

###### Abstract

Conventional score\-based diffusion models learn scores without representing normalized densities, whereas tractable normalized models support both sampling and direct likelihood evaluation\. A recent tensor\-network approach provides such a representation but is largely restricted to low\-dimensional lattices\. Here we introduce a discrete diffusion model that parameterizes normalized probability distributions using variational autoregressive networks\. Explicit Markov jump operators govern the forward noising and reverse denoising dynamics, extending discrete diffusion models with normalized distributions to spin systems on higher\-dimensional lattices\. We apply this framework to the two\- and three\-dimensional Ising models across ordered, critical, and disordered regimes, accurately computing thermodynamic quantities including free energy, energy, and magnetization\. We further integrate the framework with Monte Carlo sampling, using adaptive diffusion steps to maintain high acceptance rates even at low temperatures while enhancing sample diversity\. These results establish a neural\-network framework for the discrete diffusion model with normalized probability distributions\.

###### Keywords:

Discrete diffusion model, Markov chain Monte Carlo, Stochastic dynamics

## IIntroduction

Inspired by non\-equilibrium thermodynamics, diffusion models constitute a powerful class of generative models\[[43](https://arxiv.org/html/2609.27306#bib.bib1)\]\. They follow a two\-stage pipeline in which a forward process gradually corrupts samples from the target distribution toward a tractable prior \(Uniform, Gaussian, etc\.\), and a trained reverse process removes this corruption to generate new samples\[[3](https://arxiv.org/html/2609.27306#bib.bib47),[54](https://arxiv.org/html/2609.27306#bib.bib46),[23](https://arxiv.org/html/2609.27306#bib.bib6)\]\. Nevertheless, many classical diffusion models were originally designed for continuous\-valued data, such as images\[[23](https://arxiv.org/html/2609.27306#bib.bib6),[4](https://arxiv.org/html/2609.27306#bib.bib10),[24](https://arxiv.org/html/2609.27306#bib.bib11)\]and audio\[[39](https://arxiv.org/html/2609.27306#bib.bib12),[29](https://arxiv.org/html/2609.27306#bib.bib48)\], and cannot be applied directly to discrete\-state systems\.

To adapt to different data modalities, discrete diffusion models corrupt their inputs through explicit state transitions rather than additive Gaussian noise\[[2](https://arxiv.org/html/2609.27306#bib.bib7)\], and they have achieved substantial progress in numerous practical domains, such as natural language processing\[[32](https://arxiv.org/html/2609.27306#bib.bib8),[1](https://arxiv.org/html/2609.27306#bib.bib9),[56](https://arxiv.org/html/2609.27306#bib.bib25),[26](https://arxiv.org/html/2609.27306#bib.bib28),[49](https://arxiv.org/html/2609.27306#bib.bib40),[58](https://arxiv.org/html/2609.27306#bib.bib41)\], biological sequence generation\[[22](https://arxiv.org/html/2609.27306#bib.bib2),[38](https://arxiv.org/html/2609.27306#bib.bib33),[15](https://arxiv.org/html/2609.27306#bib.bib26),[48](https://arxiv.org/html/2609.27306#bib.bib32)\], and the study of spin systems\[[37](https://arxiv.org/html/2609.27306#bib.bib3),[42](https://arxiv.org/html/2609.27306#bib.bib4),[59](https://arxiv.org/html/2609.27306#bib.bib27)\]\. Nevertheless, the most widely used diffusion models are score based: they learn the gradient of the log density and therefore represent only unnormalized densities\. This creates three difficulties, which are especially challenging in statistical physics\. First, learning the score over the full support of a high\-dimensional, sparsely sampled distribution is hard; second, reaching the noise prior requires evolving the forward process for a long time, whereas finite\-time denoising in practice introduces reconstruction error\[[14](https://arxiv.org/html/2609.27306#bib.bib45)\]; and third, evaluating the likelihood—and hence the free energy—of generated samples from a learned score function is non\-trivial\.

Recent work on discrete diffusion models based on tensor networks directly confronted these limitations\[[10](https://arxiv.org/html/2609.27306#bib.bib5)\]\. By encoding the normalized distribution as a matrix product state \(MPS\) and the diffusion operators as matrix product operators, this approach obtains both a normalized distribution and an exact likelihood—precisely the quantities that score\-based models cannot provide\. This is an important advance that we build upon\. Its principal limitation, however, is representational: MPS is intrinsically adapted to quasi\-one\-dimensional geometries, and the associated computational cost grows prohibitively with the spatial dimension\. Indeed, on moderately sized two\-dimensional cylindrical Ising models \(8×308\\times 30and9×309\\times 30\) the tensor\-network magnetization per site already deviates from Monte Carlo benchmarks at the critical temperature, and a direct application to three\-dimensional lattices is not feasible\. These gaps motivate a normalized discrete diffusion model whose representation is not constrained by the lattice dimension\.

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig1.1.png)Figure 1:Discrete diffusion model by denoising normalized distribution of variational autoregressive networks\.Standard diffusion models adopt a forward noising process and recover samples via a learned score function, without yielding explicit and tractable normalized probability representations for physical simulations\. Prior work based on tensor networks addressed this by encoding normalized distributions as matrix product states, and representing the operators acting on them as matrix product operators\. Building on that idea, we replace tensor networks with VAN, which learns the normalized distribution at each diffusion step under a fixed transition operator\. Compared with tensor networks, our VAN method achieves improved performance in high\-dimensional systems and makes diffusion steps tractable\.In this work, we replace tensor networks with variational autoregressive networks \(VAN\)\[[52](https://arxiv.org/html/2609.27306#bib.bib17),[57](https://arxiv.org/html/2609.27306#bib.bib58)\], an approach that has proven effective for problems in chemical reaction networks\[[47](https://arxiv.org/html/2609.27306#bib.bib49),[50](https://arxiv.org/html/2609.27306#bib.bib50)\], computational biology\[[41](https://arxiv.org/html/2609.27306#bib.bib51)\], and quantum many\-body systems\[[9](https://arxiv.org/html/2609.27306#bib.bib52),[40](https://arxiv.org/html/2609.27306#bib.bib53),[33](https://arxiv.org/html/2609.27306#bib.bib54)\]\. Fig\.[1](https://arxiv.org/html/2609.27306#S1.F1)explains our research motivation and depicts the framework of our VAN‑based discrete diffusion model\. Unlike conventional approaches, which learn score functions and thus represent only unnormalized densities, and unlike the matrix\-product\-state construction, whose representation is tied to quasi\-one\-dimensional geometries, our framework retains an explicitly normalized distribution while using an autoregressive ansatz in place of a matrix product state\. This yields three concrete gains: \(i\) the likelihood is available at every diffusion step, so the variational free energy can be evaluated directly; \(ii\) the ansatz is agnostic to the lattice geometry, so the same network applies to two\- and three\-dimensional lattices alike, without introducing a bond dimension\. We validate the precise modeling capability of our method for spin systems through experiments on 2D and 3D Ising models\[[11](https://arxiv.org/html/2609.27306#bib.bib29)\], benchmarking thermodynamic observables against tensor network baselines and Wolff cluster Monte Carlo across ordered, critical, and disordered temperature regimes; and \(iii\) meanwhile, machine\-learning\-augmented Monte Carlo methods have attracted extensive research attention in recent statistical physics work\[[10](https://arxiv.org/html/2609.27306#bib.bib5),[16](https://arxiv.org/html/2609.27306#bib.bib13),[17](https://arxiv.org/html/2609.27306#bib.bib44)\]\. Adopting the same sampling paradigm utilized in the tensor network approach, we further integrate our VAN\-based discrete diffusion model with Markov chain Monte Carlo \(MCMC\), which produces non\-local updates that break the locality bottleneck of single\-spin\-flip sampling, sustains high acceptance rates down to low temperatures, and enhances sample diversity\.

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig1.2.png)Figure 2:Overview of our discrete diffusion model: training procedure, inference strategy, and application\.\(a\) The autoregressive mechanism and diffusion evolution enable tractable optimization of the loss function\. \(b\) We consider two denoising strategies for inference: an exact denoising scheme for the Euler‑discretized forward noising process, and an accelerated denoising scheme based on the tau\-leaping algorithm\. \(c\) We integrate the discrete diffusion model with the MCMC sampler to enable efficient sampling of spin configurations\.
## IIBackground

### II\.1Variational autoregressive networks

We consider the problem of sampling the target distributionPPdefined over a finite discrete state space𝒮=\{−1,\+1\}𝒟\\mathcal\{S\}=\\\{\-1,\+1\\\}^\{\\mathcal\{D\}\}\. Specifically, the target distribution adopted in our study follows the Boltzmann distribution

P⁡\(𝝈\)=1Z​e−β​E​\(𝝈\),P\(\\bm\{\\sigma\}\)=\\frac\{1\}\{Z\}e^\{\-\\beta E\(\\bm\{\\sigma\}\)\},\(1\)whereβ=1/T\\beta=1/Tdenotes the inverse temperature,E⁡\(𝝈\)E\(\\bm\{\\sigma\}\)is a known energy function, and𝝈=\(σ1,σ2,…,σ𝒟\)\\bm\{\\sigma\}=\(\\sigma\_\{1\},\\sigma\_\{2\},\\dots,\\sigma\_\{\\mathcal\{D\}\}\)denotes the spin configuration, with each discrete spin variable satisfyingσi∈\{−1,\+1\}\\sigma\_\{i\}\\in\\\{\-1,\+1\\\}\. The termZZdenotes the partition function, a critical normalization constant defined as the sum of exponential energy terms over all configurations of the system

Z=∑𝝈e−β​E​\(𝝈\)\.Z=\\sum\_\{\\bm\{\\sigma\}\}e^\{\-\\beta E\(\\bm\{\\sigma\}\)\}\.\(2\)The spin system has a state space containing2𝒟2^\{\\mathcal\{D\}\}distinct configurations, so the number of configurations grows exponentially with the system size𝒟\\mathcal\{D\}\. This exponential scaling of the configuration space means that the exact numerical computation ofZZbecomes computationally intractable for systems with large𝒟\\mathcal\{D\}\.

To address this limitation, the VAN achieves tractable modeling via a variational autoregressive factorization scheme in Fig\.[2](https://arxiv.org/html/2609.27306#S1.F2)\(a\), which decomposes the intractable joint probability distribution into a product of conditional probabilities in a prescribed sequential order

Pθ​\(𝝈\)=∏i=1𝒟Pθ​\(σi∣σ1,…,σi−1\),P^\{\\theta\}\(\\bm\{\\sigma\}\)=\\prod\_\{i=1\}^\{\\mathcal\{D\}\}P^\{\\theta\}\(\\sigma\_\{i\}\\mid\\sigma\_\{1\},\\dots,\\sigma\_\{i\-1\}\),\(3\)whereθ\\thetadenotes the learnable parameters\. Owing to its inherent autoregressive architecture, the VAN possesses two notable features for the Boltzmann distribution approximation\. First, it is capable of capturing long\-range spatial correlations among spin variables\. Second, the sequential factorization structure supports precise likelihood evaluation of generated configurationsPθ​\(𝝈\)P^\{\\theta\}\(\\bm\{\\sigma\}\)and autoregressive sampling, allowing independent configuration generation in a predefined order\.

The training objective of VAN is to optimize the parameterized variational ansatzPθP^\{\\theta\}to precisely fit the target Boltzmann distributionPP\. To quantitatively measure distribution mismatch and guide the optimization process, we utilize the Kullback\-Leibler \(KL\) divergence, a standard measure for evaluating the difference between two probability distributions\. Specifically, we employ reverse KL divergence, which only requires explicit knowledge of the energy function and allows model training using samples drawn directly from the model itself:

DKL\(Pθ\|\|P\)\\displaystyle D\_\{\\text\{KL\}\}\(P^\{\\theta\}\|\|P\)=∑𝝈Pθ​\(𝝈\)​log⁡\(Pθ​\(𝝈\)P⁡\(𝝈\)\)\\displaystyle=\\sum\_\{\\bm\{\\sigma\}\}P^\{\\theta\}\(\\bm\{\\sigma\}\)\\log\\left\(\\frac\{P^\{\\theta\}\(\\bm\{\\sigma\}\)\}\{P\(\\bm\{\\sigma\}\)\}\\right\)=β⁡\(F⁡\(θ\)−F\),\\displaystyle=\\beta\(F\(\\theta\)\-F\),\(4\)where

F=−1β​log⁡Z\\displaystyle F=\-\\frac\{1\}\{\\beta\}\\log Z\(5\)represents the free energy of the system\. Owing to the non\-negativity of the KL divergence, the variational free energy

F⁡\(θ\)=1β​∑𝝈Pθ​\(𝝈\)​\(log⁡Pθ​\(𝝈\)\+β​E​\(𝝈\)\)\\displaystyle F\(\\theta\)=\\frac\{1\}\{\\beta\}\\sum\_\{\\bm\{\\sigma\}\}P^\{\\theta\}\(\\bm\{\\sigma\}\)\\left\(\\log P^\{\\theta\}\(\\bm\{\\sigma\}\)\+\\beta E\(\\bm\{\\sigma\}\)\\right\)\(6\)serves as an upper bound of the exact free energyFF\. MinimizingF⁡\(θ\)F\(\\theta\)narrows the gap between this variational free energyF⁡\(θ\)F\(\\theta\)and the exact free energyFF, while driving the parameterized distributionPθ​\(𝝈\)P^\{\\theta\}\(\\bm\{\\sigma\}\)to gradually converge to the exact distribution of the system\.

Compared to tensor networks\[[30](https://arxiv.org/html/2609.27306#bib.bib34),[53](https://arxiv.org/html/2609.27306#bib.bib35)\], the key advantage of VAN is the ability to approximate complex target distributions with its superior generality\. Unlike MPS, which are inherently limited to one\-dimensional systems, the neural network\-based ansatz of VAN can be applied to systems with arbitrary topological structures\. For example, convolutional neural networks \(CNN\)\[[27](https://arxiv.org/html/2609.27306#bib.bib42)\]are highly suitable for modeling two\-dimensional systems, as their convolutional architecture efficiently captures local spatial structure in 2D data, similar to how it processes image information\.

VAN has proven to be a powerful tool for studying equilibrium statistical mechanics\[[35](https://arxiv.org/html/2609.27306#bib.bib39),[34](https://arxiv.org/html/2609.27306#bib.bib37),[5](https://arxiv.org/html/2609.27306#bib.bib38)\]\. Beyond equilibrium problems, it is also applicable to exploring the time evolution dynamics of nonequilibrium systems in statistical mechanics\[[46](https://arxiv.org/html/2609.27306#bib.bib23)\]\. By accurately approximating the normalized probability distribution at each time step, the VAN enables efficient calculation of the dynamical partition function — a core quantity for characterizing dynamical phase transitions and extracting essential dynamical observables\. This capability broadens the research scope of statistical mechanics, allowing researchers to reveal finite\-time dynamical phase diagrams, critical scaling laws, and emergent spatial structures that are challenging for conventional analytical and numerical approaches\. Discrete diffusion models are designed to describe unidirectional stochastic Markov processes and thus belong to the class of nonequilibrium dynamical systems\. Leveraging the strong capability of VAN in capturing nonequilibrium evolution, we integrate the above VAN\-based dynamic learning framework into discrete diffusion models\.

### II\.2Discrete noising process

We formulate the noising dynamics as a continuous time Markov chain \(CTMC\), which is a stochastic process on a discrete state space𝒮\\mathcal\{S\}with the Markov property\. The time evolution of the probability distributionPtP\_\{t\}satisfies the master equation

d​Ptd​t=W​Pt,\\displaystyle\\frac\{dP\_\{t\}\}\{dt\}=WP\_\{t\},\(7\)where the standard form of the Markov generatorWWis given by

W𝝈′​𝝈=\{w⁡\(𝝈′←𝝈\),if​𝝈′≠𝝈−R⁡\(𝝈\),if​𝝈′=𝝈\.\\displaystyle W\_\{\\bm\{\\sigma\}^\{\\prime\}\\bm\{\\sigma\}\}=\\begin\{cases\}w\(\\bm\{\\sigma\}^\{\\prime\}\\leftarrow\\bm\{\\sigma\}\),&\\text\{if \}\\bm\{\\sigma\}^\{\\prime\}\\neq\\bm\{\\sigma\}\\\\ \-R\(\\bm\{\\sigma\}\),&\\text\{if \}\\bm\{\\sigma\}^\{\\prime\}=\\bm\{\\sigma\}\.\\end\{cases\}\(8\)
We consider the bistochastic noising dynamics built from non\-interacting single\-spin flips\. Transitions are allowed only between configurations separated by a Hamming distanced=1d=1, i\.e\., configurations differing by a single spin flip\. The transition rates are defined as

w⁡\(𝝈′←𝝈\)=\{1,d⁡\(𝝈′,𝝈\)=10,otherwise\.\\displaystyle w\(\\bm\{\\sigma\}^\{\\prime\}\\leftarrow\\bm\{\\sigma\}\)=\\begin\{cases\}1,&d\(\\bm\{\\sigma\}^\{\\prime\},\\bm\{\\sigma\}\)=1\\\\ 0,&\\text\{otherwise\}\.\\end\{cases\}\(9\)The corresponding escape rates areR⁡\(𝝈\)=∑𝝈′≠𝝈w⁡\(𝝈′←𝝈\)=𝒟\.R\(\\bm\{\\sigma\}\)=\\sum\_\{\\bm\{\\sigma\}^\{\\prime\}\\neq\\bm\{\\sigma\}\}w\(\\bm\{\\sigma\}^\{\\prime\}\\leftarrow\\bm\{\\sigma\}\)=\\mathcal\{D\}\.This construction endowsWWwith two key properties: probability normalization is preserved under time evolution, and the probability distribution asymptotically converges to the uniform distribution, which serves as the stationary distribution satisfyingW​Puniform=0WP\_\{\\text\{uniform\}\}=0\.

Based on the master equation Eq\. \([7](https://arxiv.org/html/2609.27306#S2.E7)\), the evolution of the probability distribution fromt=t1t=t\_\{1\}tot=t2t=t\_\{2\}can be equivalently expressed as

Pt2θ=e\(t2−t1\)​W​Pt1θ=Qt2←t1​Pt1θ\.\\displaystyle\{P\_\{t\_\{2\}\}^\{\\theta\}\}=e^\{\(t\_\{2\}\-t\_\{1\}\)W\}\{P\_\{t\_\{1\}\}^\{\\theta\}\}=\{Q\}\_\{t\_\{2\}\\leftarrow t\_\{1\}\}\{P\_\{t\_\{1\}\}^\{\\theta\}\}\.\(10\)For a sufficiently small time intervalΔ​t=t2−t1\\Delta t=t\_\{2\}\-t\_\{1\}, the matrix can be approximated via the first‑order Euler expansion:eΔ​t​W≈𝕀\+Δ​t​We^\{\\Delta tW\}\\approx\\mathbb\{I\}\+\\Delta tW, where𝕀\\mathbb\{I\}denotes the identity matrix and𝒟​Δ​t≤1\\mathcal\{D\}\\Delta t\\leq 1\. A typical choice satisfying this bound isΔ​t=1/\(2​𝒟\)\\Delta t=1/\(2\\mathcal\{D\}\)\. Accordingly, the transition probability under this Euler approximation takes the form

pt2\|t1​\(y\|x\)=δ⁡\(y,x\)\+Δ​t​W​\(y,x\)\+o⁡\(Δ​t\),\\displaystyle p\_\{t\_\{2\}\|t\_\{1\}\}\(y\|x\)=\\delta\(y,x\)\+\\Delta tW\(y,x\)\+o\(\\Delta t\),\(11\)whereo⁡\(Δ​t\)o\(\\Delta t\)represents higher\-order infinitesimal terms vanishing faster thanΔ​t\\Delta t, andδ⁡\(y,x\)\\delta\(y,x\)is the Kronecker delta that takes value 1 whenx=yx=yand 0 otherwise\. We adopt this first‑order Euler approximation to construct the forward noising process of our discrete diffusion model, whose dynamics are governed by the single‑spin‑flip generatorWW\.

### II\.3Discrete denoising process

Efficient simulation of discrete stochastic processes is a fundamental challenge in statistical physics and computational modeling\. For our discrete diffusion model, which enables tractable probability flow between adjacent time steps, the Euler\-discretized reverse process can be analytically derived based on Bayes’ theorem\[[44](https://arxiv.org/html/2609.27306#bib.bib43)\]

πt1\|t2​\(x1∣x2\)=πt2\|t1​\(x2∣x1\)​Pt1​\(x1\)Pt2​\(x2\),t1<t2\.\\displaystyle\\pi\_\{t\_\{1\}\\mid t\_\{2\}\}\(x\_\{1\}\\mid x\_\{2\}\)=\\pi\_\{t\_\{2\}\\mid t\_\{1\}\}\(x\_\{2\}\\mid x\_\{1\}\)\\frac\{P\_\{t\_\{1\}\}\(x\_\{1\}\)\}\{P\_\{t\_\{2\}\}\(x\_\{2\}\)\},\\quad t\_\{1\}<t\_\{2\}\.\(12\)The Alg\.[1](https://arxiv.org/html/2609.27306#alg1)details the full implementation pipeline for the Euler\-discretized denoising update rules of our discrete diffusion model\.

Algorithm 1Stepwise Euler\-discretized denoising1:A transition rate

ww, time partition

0=t0<t1<⋯<tK0=t\_\{0\}<t\_\{1\}<\\cdots<t\_\{K\}\.

2:Draw

uK∼PtKθu\_\{K\}\\sim P\_\{t\_\{K\}\}^\{\\theta\}\.

3:for

k=Kk=Kdownto

11do

4:Set

uk−1=u\_\{k\-1\}=
5:

\{u,w\.p\.​Δ​tk​w​Ptk−1θ​\(u\)Ptkθ\(uk\)uk,w\.p\.​1−∑u≠ukΔ​tk​w​Ptk−1θ​\(u\)Ptkθ\(uk\)\\begin\{cases\}u,&\\text\{w\.p\. \}\\Delta t\_\{k\}w\\frac\{P\_\{t\_\{k\-1\}\}^\{\\theta\}\(u\)\}\{P\_\{t\_\{k\}\}\{\{\}^\{\\theta\}\(u\_\{k\}\)\}\}\\\\\[3\.0pt\] u\_\{k\},&\\text\{w\.p\. \}1\-\\sum\\limits\_\{u\\neq u\_\{k\}\}\\Delta t\_\{k\}w\\frac\{P\_\{t\_\{k\-1\}\}^\{\\theta\}\(u\)\}\{P\_\{t\_\{k\}\}\{\{\}^\{\\theta\}\(u\_\{k\}\)\}\}\\end\{cases\}
6:where

Δ​tk=tk−tk−1\\Delta t\_\{k\}=t\_\{k\}\-t\_\{k\-1\}and

u≠uku\\neq u\_\{k\}\.

7:endfor

8:return

u0∼Pt0θu\_\{0\}\\sim P\_\{t\_\{0\}\}^\{\\theta\}\.

The stepwise Euler\-discretized denoising method takes the form of a two\-stage algorithm outlined below:

1. \(i\)Firstly, it initializes the noisy statesuKu\_\{K\}\. The states are generated either by applying sequential flips to thet=0t=0samples following our noising rules or by direct sampling from the distributionPtKθP\_\{t\_\{K\}\}^\{\\theta\}\.
2. \(ii\)Then it iteratively reconstructs the latent statesuk−1u\_\{k\-1\}fromuku\_\{k\}\. The probability of updatinguku\_\{k\}is given byΔ​tk​w​Ptk−1θ​\(u\)/Ptkθ​\(uk\)\\Delta t\_\{k\}wP\_\{t\_\{k\-1\}\}^\{\\theta\}\(u\)/P\_\{t\_\{k\}\}^\{\\theta\}\(u\_\{k\}\)\. Instead of directly evaluatingPtkθ​\(uk\)P\_\{t\_\{k\}\}^\{\\theta\}\(u\_\{k\}\), we can compute it fromPtk−1θP\_\{t\_\{k\-1\}\}^\{\\theta\}using the noising rule to ensure that the probability of retaining the stateuku\_\{k\}remains positive\.

Although this stepwise denoising approach yields accurate results, it suffers from substantial computational overhead and low sampling efficiency, particularly for long denoising trajectories\. To balance numerical accuracy and computational cost, we incorporate the tau‑leaping into our discrete diffusion model\.

Tau\-leaping is an approximate stochastic simulation method designed to accelerate sampling of discrete\-state Markov processes, offering noticeable efficiency improvements for systems with frequent state transitions\. Originally proposed and widely employed in chemical physics, this technique accelerates simulations of large reaction networks with abundant reaction events\[[21](https://arxiv.org/html/2609.27306#bib.bib18),[7](https://arxiv.org/html/2609.27306#bib.bib20),[8](https://arxiv.org/html/2609.27306#bib.bib19)\]\. More recently, tau\-leaping has been extended to discrete diffusion models based on CTMC for an efficient approximation of the generative reverse process\[[6](https://arxiv.org/html/2609.27306#bib.bib16)\]\. Unlike the classic Gillespie algorithm\[[19](https://arxiv.org/html/2609.27306#bib.bib21),[20](https://arxiv.org/html/2609.27306#bib.bib22)\], which simulates only one transition event per iteration, tau\-leaping bundles numerous transition events in a single step to reduce computational burden\. This property makes tau\-leaping particularly suitable for high\-dimensional systems with rapid transition rates\.

Algorithm 2Tau\-leaping1:A transition rate

ww, time partition

0=t0<t1<⋯<tK0=t\_\{0\}<t\_\{1\}<\\cdots<t\_\{K\}, optimal leap duration

τ\\tau\.

2:Draw

uK∼PtKθu\_\{K\}\\sim P\_\{t\_\{K\}\}^\{\\theta\}, where

K=⌈tKτ⌉K=\\lceil\\dfrac\{t\_\{K\}\}\{\\tau\}\\rceil\.

3:Encoding map:

ϕ⁡\(uK\)=\(uK\+1\)/2\\phi\(u\_\{K\}\)=\{\(u\_\{K\}\+1\)\}/\{2\}
4:for

k=Kk=Kdownto

11do

5:for

d=1d=1to

𝒟\\mathcal\{D\}do

6:for

zd∈\{0,1\}∖\{ukd\}z^\{d\}\\in\\\{0,1\\\}\\setminus\\\{u\_\{k\}^\{d\}\\\}do

7:Draw

Md,zd∼Poisson​\(τ​w​Ptkθ​\(ukd,zd\)/Ptkθ​\(uk\)\)M\_\{d,z^\{d\}\}\\sim\\text\{Poisson\}\(\\tau wP\_\{t\_\{k\}\}^\{\\theta\}\(u\_\{k\}^\{d,z^\{d\}\}\)/P\_\{t\_\{k\}\}^\{\\theta\}\(u\_\{k\}\)\)\.

8:endfor

9:endfor

10:Set

x=uk\+∑d=1𝒟∑zd∈\{0,1\}∖\{ukd\}\(zd−ukd\)​Md,zdx=u\_\{k\}\+\\sum\_\{d=1\}^\{\\mathcal\{D\}\}\\sum\_\{z^\{d\}\\in\\\{0,1\\\}\\setminus\\\{u\_\{k\}^\{d\}\\\}\}\(z^\{d\}\-u\_\{k\}^\{d\}\)M\_\{d,z^\{d\}\}\.

11:Set

uk−1=ℳuk​\(x\)u\_\{k\-1\}=\\mathcal\{M\}\_\{u\_\{k\}\}\(x\), where

ℳuk​\(x\)=\{xd​δ𝒮​\(xd\)\+zd​\(1−δ𝒮​\(xd\)\)\}d∈\[𝒟\]\\mathcal\{M\}\_\{u\_\{k\}\}\(x\)=\\\{x^\{d\}\\delta\_\{\\mathcal\{S\}\}\(x^\{d\}\)\+z^\{d\}\(1\-\\delta\_\{\\mathcal\{S\}\}\(x^\{d\}\)\)\\\}\_\{d\\in\[\\mathcal\{D\}\]\}is a mapping from

ℤ𝒟\\mathbb\{Z\}^\{\\mathcal\{D\}\}to

𝒮𝒟\\mathcal\{S\}^\{\\mathcal\{D\}\}\.

12:endfor

13:returnDecoding map:

ϕ−1​\(u0\)=2​u0−1\\phi^\{\-1\}\(u\_\{0\}\)=\{2\\,u\_\{0\}\-1\}

Compared with the Alg\.[1](https://arxiv.org/html/2609.27306#alg1), the tau\-leaping method for our discrete diffusion model also adopts a two\-stage procedure and retains an identical first stage\. In its second stage, stepwise Euler\-discretized denoising is substituted with an approximate leap over a predetermined optimal leap durationτ\\tau\(Supplementary Appendix[C](https://arxiv.org/html/2609.27306#A3)\)\. Instead of simulating every transition separately, the tau\-leaping method assumes that the reverse transition rates and the stateuku\_\{k\}remain constant throughout the intervalτ\\tau\. We sample Poisson countsMd,zdM\_\{d,z^\{d\}\}for each coordinatedd, aggregate these counts and apply them to the stateuku\_\{k\}at the end of the intervalτ\\tauto form the statexx\.xxis then projected onto the valid discrete state space𝒮\\mathcal\{S\}via the mappingℳuk​\(⋅\)\\mathcal\{M\}\_\{u\_\{k\}\}\(\\cdot\)to produce the latent stateuk−1u\_\{k\-1\}\. Although the tau\-leaping method reduces computational cost, we revert to the stepwise Euler\-discretized denoising method when strict numerical accuracy is required\.

## IIINumerical results

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig2.5.png)Figure 3:Tensor networks and VAN\-based discrete diffusion models on 2D Ising model\.\(a\) Magnetization per site obtained via the tensor network approach\[[10](https://arxiv.org/html/2609.27306#bib.bib5)\], the Wolff cluster algorithm, and our method with averages taken over10510^\{5\}independent samples, on8×308\\times 30and9×309\\times 30lattices\. \(b\) Comparison of energy per site under the same simulation parameters\. Note that energy per site data from the tensor network approach are unavailable in the original reference\.The Ising model is defined on a discrete spin lattice with nearest\-neighbor pairwise interactions, whose Hamiltonian is given by

E\(𝝈\)=−J∑⟨i,j⟩σiσj,\\displaystyle E\(\\bm\{\\sigma\}\)=\-J\\sum\_\{\\langle i,j\\rangle\}\\sigma\_\{i\}\\sigma\_\{j\},\(13\)where⟨i,j⟩\\langle i,j\\rangledenotes nearest\-neighbor pairs of spins, and the ferromagnetic coupling strength is set toJ=1J=1\. In the thermodynamic limit, the two\-dimensional Ising model undergoes a continuous phase transition at the critical inverse temperatureβc=ln⁡\(1\+2\)/2\\beta\_\{c\}=\\ln\(1\+\\sqrt\{2\}\)/2\. The system remains disordered forβ<βc\\beta<\\beta\_\{c\}at high temperatures and develops spontaneous magnetic order forβ\>βc\\beta\>\\beta\_\{c\}at low temperatures\.

We first consider the two\-dimensional Ising model defined on a𝒟=L1×L2\\mathcal\{D\}=L\_\{1\}\\times L\_\{2\}lattice \(L2≥L1L\_\{2\}\\geq L\_\{1\}\) with cylindrical boundary conditions, where periodic boundary conditions \(PBC\) are applied along the first dimension and open boundary conditions \(OBC\) along the second dimension\. We adopt lattice sizes of8×308\\times 30and9×309\\times 30at the critical inverse temperatureβc\\beta\_\{c\}, which produce numerical biases in the tensor network approach compared to standard Monte Carlo benchmarks\[[10](https://arxiv.org/html/2609.27306#bib.bib5)\]\. To quantify the accuracy of the learned distributions, we define the absolute magnetization per site as

m=\|1𝒟​∑i=1𝒟σi\|,\\displaystyle m=\\left\|\\frac\{1\}\{\\mathcal\{D\}\}\\sum\_\{i=1\}^\{\\mathcal\{D\}\}\\sigma\_\{i\}\\right\|,\(14\)and compare both the magnetization and the energy per site against the reference data\. The comparison results in Fig\.[3](https://arxiv.org/html/2609.27306#S3.F3)demonstrate that our VAN method achieves higher accuracy in predicting magnetization per site than tensor networks, while producing reliable estimates for the energy per site\. We adopt results from the Wolff cluster algorithm \(Supplementary Alg\.[3](https://arxiv.org/html/2609.27306#alg3)\) as reference data in our work\.

We further perform systematic simulations on both 2D and 3D Ising models with full PBC: for 2D lattices, all horizontal and vertical boundaries are periodically connected; for 3D lattices, periodicity is imposed along all three spatial dimensions, which guarantees uniform nearest\-neighbor spin interactions across the entire system\. We compute the free energy per site, energy per site, and magnetization per site for lattices of size16×1616\\times 16\(2D\) and4×4×44\\times 4\\times 4\(3D\) over a wide temperature range spanning the low\-temperature ordered regime, the high\-temperature disordered regime, and the phase transition point\. Three\-dimensional lattices are substantially more demanding for tensor networks, whose computational cost and memory usage grow exponentially with the spatial dimension and severely hinder the accurate characterization of long\-range spin correlations\. This is precisely the regime in which a neural autoregressive representation remains tractable, and the resulting estimates are summarized in Fig\.[4](https://arxiv.org/html/2609.27306#S4.F4)\. Since the 3D Ising model has no exact solution for free energy, we computed the free energy using FlashVAN\[[57](https://arxiv.org/html/2609.27306#bib.bib58)\]—currently the state\-of\-the\-art method—and used the resulting values as a benchmark for comparison\.

## IVIntegrating discrete diffusion into MCMC proposal updates

From a theoretical perspective, perfect alignment between the learned VANPθP^\{\\theta\}and the target equilibrium distribution would render autoregressive sampling fully unbiased, eliminating the need for supplementary MCMC correction\. In practice, however, practical neural network approximations inevitably suffer from a distribution mismatch against the true distribution, due to limited training data, inherent architectural constraints of neural networks, and fundamental optimization barriers\. A typical symptom of this discrepancy is mode collapse, inducing biased results and poor sample diversity\. Furthermore, the learned distribution may correspond to a high temperature distribution, whereas our objective is to sample the target distribution at lower temperatures\. Instead of merely retraining the neural network to eliminate such distribution gaps, we introduce a proposal module based on a discrete diffusion model \(DDM\) embedded in the MCMC workflow, which enables efficient sampling of the desired target distribution\[[25](https://arxiv.org/html/2609.27306#bib.bib14),[12](https://arxiv.org/html/2609.27306#bib.bib57)\]\.

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig4.png)Figure 4:Learned discrete diffusion models for 2D and 3D Ising models at distinct temperatures\.\(a\) Free energy per site, \(b\) energy per site and \(c\) magnetization per site as functions of inverse temperatureβ\\beta\. The upper row corresponds to the𝒟=16×16\\mathcal\{D\}=16\\times 16two\-dimensional Ising model with temperaturesT=1/β=1/ln⁡\(1\+2\),1\.5000,2\.0000,2/ln⁡\(1\+2\),2\.5000,3\.0000,3\.5000,4\.0000,4/ln⁡\(1\+2\)T=1/\\beta=1/\\ln\(1\+\\sqrt\{2\}\),1\.5000,2\.0000,2/\\ln\(1\+\\sqrt\{2\}\),2\.5000,3\.0000,3\.5000,4\.0000,4/\\ln\(1\+\\sqrt\{2\}\)\. The lower row presents statistical results for the𝒟=4×4×4\\mathcal\{D\}=4\\times 4\\times 4three\-dimensional Ising model at temperaturesT=1/β=3\.5000,4\.0000,4\.5115,5\.0000,5\.5000,6\.0000T=1/\\beta=3\.5000,4\.0000,4\.5115,5\.0000,5\.5000,6\.0000\. All systems are simulated with periodic boundary conditions and all results are averaged over10510^\{5\}independent samples\.To ensure reliable MCMC sampling, we need to train a discrete diffusion model to recover corrupted samples and steer them toward reasonable data distributions\. Unlike standard discrete diffusion models, our method adopts a fundamentally different stepwise training scheme for modeling sequential distributions\. To adapt to this specialized training paradigm, we first formulate step\-by\-step loss functions to optimize the VAN parameters

loss⁡\(Ptθ\)=\{𝔼𝝈∼Ptθ​\[log⁡Ptθ​\(𝝈\)\+β​E​\(𝝈\)\],t=0𝔼𝝈∼Ptθ​\[log⁡Ptθ​\(𝝈\)−log⁡Pt​\(𝝈\)\],t\>0\.\\displaystyle\\mathrm\{loss\}\(P\_\{t\}^\{\\theta\}\)=\\begin\{cases\}\\mathbb\{E\}\_\{\\bm\{\\sigma\}\\sim P\_\{t\}^\{\\theta\}\}\\left\[\\log P\_\{t\}^\{\\theta\}\(\\bm\{\\sigma\}\)\+\\beta E\(\\bm\{\\sigma\}\)\\right\],&t=0\\\\ \\mathbb\{E\}\_\{\\bm\{\\sigma\}\\sim P\_\{t\}^\{\\theta\}\}\\left\[\\log P\_\{t\}^\{\\theta\}\(\\bm\{\\sigma\}\)\-\\log P\_\{t\}\(\\bm\{\\sigma\}\)\\right\],&t\>0\.\\end\{cases\}\(15\)Concretely, the training pipeline of our discrete diffusion model is implemented in two stages:

1. \(i\)We first train a VAN to fit the target distribution, which serves as the initial distribution of the diffusion process at stept=0t=0\.
2. \(ii\)Given the trained distributionPt−1θP\_\{t\-1\}^\{\\theta\}, we employ the transition ratewwand the time stepΔ​t\\Delta tto compute the supervised training targetPt​\(𝝈\)P\_\{t\}\(\\bm\{\\sigma\}\), which is required to train the VAN of the subsequent diffusion step\. Iterative execution of this procedure constructs the complete forward noising process in an explicit and deterministic manner\. For allt\>0t\>0, we initialize the VAN with the trained parameters from the previous step and fine\-tune it with a small number of optimization iterations to obtain the distributionPtθP\_\{t\}^\{\\theta\}\.

Equipped with the trained discrete diffusion model, we construct our adaptive MCMC sampler\. Fig\.[2](https://arxiv.org/html/2609.27306#S1.F2)\(c\) illustrates its complete framework and pipeline, and each MCMC iteration proceeds as follows\.

\(1\) All MCMC chains are initialized with spin lattice configurations𝝈\\bm\{\\sigma\}sampled autoregressively from the VAN\-learned distributionP0θP\_\{0\}^\{\\theta\}\. Under this autoregressive sampling scheme, the spin variableσi\\sigma\_\{i\}at each lattice siteiiis sampled from its conditional distribution given all previously sampled spins

σi∼Bernoulli⁡\(P0θ​\(σi∣σ1,…,σi−1\)\)\.\\displaystyle\\sigma\_\{i\}\\sim\\mathrm\{Bernoulli\}\\left\(P\_\{0\}^\{\\theta\}\\left\(\\sigma\_\{i\}\\mid\\sigma\_\{1\},\\dots,\\sigma\_\{i\-1\}\\right\)\\right\)\.\(16\)
\(2\) Then the initial spin configurations𝝈\\bm\{\\sigma\}are progressively corrupted through a sequence of forward diffusion stepsttto produce noisy states𝒖\\bm\{u\}, which disrupt the original structural features of the spin lattice\. Starting from these corrupted configurations𝒖\\bm\{u\}, we invert the noising process via either the Euler\-discretized denoising or the tau‑leaping method to reconstruct candidate spin states𝝈′\\bm\{\\sigma\}^\{\\prime\}\. The corresponding transition probability of this noising\-and\-denoising proposal process is given by

π⁡\(𝝈′∣𝝈\)\\displaystyle\\pi\(\{\\bm\{\\sigma\}\}^\{\\prime\}\\mid\{\\bm\{\\sigma\}\}\)=∑𝒖Pt\|0​\(𝒖∣𝝈\)​P0\|t​\(𝝈′∣𝒖\)\\displaystyle=\\sum\_\{\\bm\{u\}\}P\_\{t\\mid 0\}\(\{\\bm\{u\}\\mid\\bm\{\\sigma\}\}\)P\_\{0\\mid t\}\(\{\\bm\{\\sigma\}\}^\{\\prime\}\\mid\\bm\{u\}\)=∑𝒖Pt\|0​\(𝒖∣𝝈\)​Pt\|0​\(𝒖∣𝝈′\)​P0θ​\(𝝈′\)Ptθ​\(𝒖\)\.\\displaystyle=\\sum\_\{\\bm\{u\}\}P\_\{t\\mid 0\}\(\{\\bm\{u\}\\mid\\bm\{\\sigma\}\}\)P\_\{t\\mid 0\}\(\{\\bm\{u\}\\mid\{\\bm\{\\sigma\}\}^\{\\prime\}\}\)\\frac\{P\_\{0\}^\{\\theta\}\(\{\\bm\{\\sigma\}\}^\{\\prime\}\)\}\{P\_\{t\}^\{\\theta\}\(\\bm\{u\}\)\}\.\(17\)
\(3\) To draw samples that conform to the target equilibrium distribution, candidate spin configurations𝝈′\\bm\{\\sigma\}^\{\\prime\}are further accepted or rejected according to the standard Metropolis–Hastings acceptance probability, defined as

Acc\[𝝈→𝝈′\]\\displaystyle\\text\{Acc\}\\left\[\{\\bm\{\\sigma\}\}\\rightarrow\{\\bm\{\\sigma\}\}^\{\\prime\}\\right\]=min⁡\[1,Ptrue​\(𝝈′\)​π​\(𝝈∣𝝈′\)Ptrue​\(𝝈\)​π​\(𝝈′∣𝝈\)\]\\displaystyle=\\min\\left\[1,\\frac\{P\_\{\\text\{true\}\}\(\{\\bm\{\\sigma\}\}^\{\\prime\}\)\\pi\(\{\\bm\{\\sigma\}\}\\mid\{\\bm\{\\sigma\}\}^\{\\prime\}\)\}\{P\_\{\\text\{true\}\}\(\{\\bm\{\\sigma\}\}\)\\pi\(\{\\bm\{\\sigma\}\}^\{\\prime\}\\mid\{\\bm\{\\sigma\}\}\)\}\\right\]=min⁡\[1,e−β​E​\(𝝈′\)​P0θ​\(𝝈\)e−β​E​\(𝝈\)​P0θ​\(𝝈′\)\],\\displaystyle=\\min\\left\[1,\\frac\{e^\{\-\\beta E\(\{\\bm\{\\sigma\}\}^\{\\prime\}\)\}P\_\{0\}^\{\\theta\}\(\{\\bm\{\\sigma\}\}\)\}\{e^\{\-\\beta E\(\{\\bm\{\\sigma\}\}\)\}P\_\{0\}^\{\\theta\}\(\{\\bm\{\\sigma\}\}^\{\\prime\}\)\}\\right\],\(18\)where the simplification adopts the conditional probability formula in Eq\. \([17](https://arxiv.org/html/2609.27306#S4.E17)\)\. When employing tau‑leaping for denoising, a sufficient numerical accuracy must be maintained to ensure that Eq\. \([18](https://arxiv.org/html/2609.27306#S4.E18)\) is satisfied\. In Fig\.[5](https://arxiv.org/html/2609.27306#S4.F5), we investigate how acceptance and decorrelation vary with diffusion steps and target‑distribution temperature\. At fixed diffusion steps, a larger temperature gap lowers acceptance and leads to weaker decorrelation\. Interestingly, more diffusion step updates reduce acceptance, but the additional spin flips applied to candidate configurations enhance decorrelation upon acceptance\. We define the decorrelation coefficient as1−C1\-C, whereCCis the normalized correlation coefficient

C=∑d=1𝒟\(1N−1​∑i=1N\(σi,d\(k\)−σ¯d\(k\)\)​\(σi,d\(k\+1\)−σ¯d\(k\+1\)\)\)\(∑d=1𝒟Var​\(𝝈d\(k\)\)\)​\(∑d=1𝒟Var​\(𝝈d\(k\+1\)\)\)\\displaystyle C=\\frac\{\\displaystyle\\sum\_\{d=1\}^\{\\mathcal\{D\}\}\\left\(\\frac\{1\}\{N\-1\}\\sum\_\{i=1\}^\{N\}\\left\(\\sigma\_\{i,d\}^\{\(k\)\}\-\\bar\{\\sigma\}\_\{d\}^\{\(k\)\}\\right\)\\left\(\\sigma\_\{i,d\}^\{\(k\+1\)\}\-\\bar\{\\sigma\}\_\{d\}^\{\(k\+1\)\}\\right\)\\right\)\}\{\\displaystyle\\sqrt\{\\left\(\\sum\_\{d=1\}^\{\\mathcal\{D\}\}\\text\{Var\}\\\!\\left\(\\bm\{\\sigma\}^\{\(k\)\}\_\{d\}\\right\)\\right\)\\left\(\\sum\_\{d=1\}^\{\\mathcal\{D\}\}\\text\{Var\}\\\!\\left\(\\bm\{\\sigma\}^\{\(k\+1\)\}\_\{d\}\\right\)\\right\)\}\}\(19\)to evaluate statistical independence between successive sampled spin configurations\(𝝈\(k\),𝝈\(k\+1\)\)\\left\(\\bm\{\\sigma\}^\{\(k\)\},\\bm\{\\sigma\}^\{\(k\+1\)\}\\right\)\. Here,𝝈\(k\+1\)\\bm\{\\sigma\}^\{\(k\+1\)\}denotes the configuration obtained after the Metropolis–Hastings acceptance step\. Our adaptive strategy can effectively balance acceptance and decorrelation to enhance the quality of samples: we increase the diffusion timettif the acceptance probability exceeds the preset target threshold and decreasettotherwise\. The spin configurations obtained from the previous MCMC iteration are then adopted as initial states for the next iteration, and this sampling pipeline is run iteratively until all Markov chains achieve convergence\.

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig2.png)Figure 5:Sampling characteristics for different distributions and diffusion steps\.\(a\) Acceptance rate as a function of diffusion steps forN=128N=128independent configurations of the16×1616\\times 16two\-dimensional Ising model in the first MCMC iteration, each initialized at inverse temperatureβ=βc/2\\beta=\\beta\_\{c\}/2and quenched to target inverse temperaturesβ=0\.25,0\.30,0\.35\\beta=0\.25,0\.30,0\.35and0\.400\.40\(repeated over 40 independent runs\)\. \(b\) Decorrelation coefficient versus diffusion steps for the same set of configurations at each target temperature\. The results reveal two consistent trends across all tested temperatures: a larger temperature gap reduces both the acceptance and configuration decorrelation, while increasing the diffusion steps lowers the acceptance via a sharp drop at early steps but improves configuration decorrelation\.Our MCMC sampler exploits the two\-stage noising–denoising dynamics of discrete diffusion models to construct high quality proposal samples for spin lattice systems\. In our framework, the forward diffusion process introduces a sequence of unconstrained random spin flips to diversify the lattice structures, which effectively shifts samples to a “higher\-temperature” regime with a much more accessible state space for broad exploration\[[12](https://arxiv.org/html/2609.27306#bib.bib57)\]\. The subsequent reverse denoising process imposes structured constraints on the spin flips to guide the empirical distribution of the generated samples to align with the original distributionP0θP\_\{0\}^\{\\theta\}, reverting configurations to the “lower\-temperature” regime\. By restricting spin flips, this denoising procedure mitigates large deviations of candidate samples from the target distribution that arise from random spin flips, which moderately increases the acceptance rate in the subsequent Metropolis–Hastings acceptance step\. By evolving spin configurations through complete forward corruption and reverse denoising reconstruction, our method generates non\-local lattice perturbations, fundamentally breaking the spatial locality bottleneck of the conventional Local Metropolis Monte Carlo \(LMMC; Supplementary Alg\.[4](https://arxiv.org/html/2609.27306#alg4)\)\. Traditional LMMC relies on limited local updates, which frequently trap Markov chains in free\-energy local minima during long simulations and fail to fully explore the entire configuration space\.

Our MCMC method built on discrete diffusion models corresponds to the connected update proposed in prior work\[[10](https://arxiv.org/html/2609.27306#bib.bib5)\], where the current spin configuration𝝈\\bm\{\\sigma\}and the proposed candidate configuration𝝈′\\bm\{\\sigma\}^\{\\prime\}are statistically correlated\. Meanwhile, the prior work also presents a disconnected update scheme, which generates a fully independent candidate configuration𝝈′\\bm\{\\sigma\}^\{\\prime\}with no correlation to the original state𝝈\\bm\{\\sigma\}\. In our work, we adopt the neural MCMC \(NMCMC\) to implement such disconnected updates \(Supplementary Alg\.[5](https://arxiv.org/html/2609.27306#alg5)\)\. The generated candidate configuration is accepted or rejected following the acceptance probability

Acc\[𝝈→𝝈′\]\\displaystyle\\text\{Acc\}\\left\[\{\\bm\{\\sigma\}\}\\rightarrow\{\\bm\{\\sigma\}\}^\{\\prime\}\\right\]=min⁡\[1,Ptrue​\(𝝈′\)​π​\(𝝈∣𝝈′\)Ptrue​\(𝝈\)​π​\(𝝈′∣𝝈\)\]\\displaystyle=\\min\\left\[1,\\frac\{P\_\{\\text\{true\}\}\(\{\\bm\{\\sigma\}\}^\{\\prime\}\)\\pi\(\{\\bm\{\\sigma\}\}\\mid\{\\bm\{\\sigma\}\}^\{\\prime\}\)\}\{P\_\{\\text\{true\}\}\(\{\\bm\{\\sigma\}\}\)\\pi\(\{\\bm\{\\sigma\}\}^\{\\prime\}\\mid\{\\bm\{\\sigma\}\}\)\}\\right\]=min⁡\[1,e−β​E​\(𝝈′\)​P0θ​\(𝝈\)e−β​E​\(𝝈\)​P0θ​\(𝝈′\)\],\\displaystyle=\\min\\left\[1,\\,\\frac\{e^\{\-\\beta E\(\\bm\{\\sigma\}^\{\\prime\}\)\}P\_\{0\}^\{\\theta\}\(\\bm\{\\sigma\}\)\}\{e^\{\-\\beta E\(\\bm\{\\sigma\}\)\}P\_\{0\}^\{\\theta\}\(\\bm\{\\sigma\}^\{\\prime\}\)\}\\right\],\(20\)where the candidate state𝝈′\\bm\{\\sigma\}^\{\\prime\}is drawn from the marginal proposal distributionP0θP\_\{0\}^\{\\theta\}parameterized by the VAN\. Since𝝈′\\bm\{\\sigma\}^\{\\prime\}is statistically independent of the current configuration𝝈\\bm\{\\sigma\}, the transition kernel satisfiesπ⁡\(𝝈′∣𝝈\)=P0θ​\(𝝈′\)\\pi\(\\bm\{\\sigma\}^\{\\prime\}\\mid\\bm\{\\sigma\}\)=P\_\{0\}^\{\\theta\}\(\\bm\{\\sigma\}^\{\\prime\}\)\.

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig5.png)Figure 6:Monte Carlo quench\.\(a\) Schematic of the quench process, where spin lattice configurations transition from higher\-temperature initial states to lower\-temperature target states\. \(b\) Acceptance rates over 10 000 sampling iterations \(data points displayed at 10\-iteration intervals\), evaluated across 64 independent configurations of the𝒟=16×16\\mathcal\{D\}=16\\times 162D Ising model undergoing a quench fromT=2T=2toT=1\.5T=1\.5\. Three sampling schemes are compared: DDM \(discrete diffusion model; connected updates; adaptive acceptance target 0\.5\), NMCMC \(neural Markov chain Monte Carlo; disconnected updates\), and LMMC \(local Metropolis Monte Carlo\)\. \(c\) Comparison of magnetization‑based effective sample size \(ESS\) fraction and sample diversity, where the ESS fraction estimates the proportion of statistically independent samples in MCMC chains and sample diversity quantifies the number of distinct spin configurations explored during sampling\. Simulations adopt the same parameters as panel \(b\), and results are averaged over1010independent runs\.In Fig\.[6](https://arxiv.org/html/2609.27306#S4.F6), we present the performance of all sampling methods considered under Monte Carlo “quench” protocols\. Notably, equipped with non\-local update mechanisms realized via connected and disconnected update schemes, our two proposed methods sustain acceptance rates higher than LMMC throughout the low‑temperature regime, especially the connected update scheme\. A naive strategy to break the locality barrier of LMMC is to flip multiple random spins per iteration; unfortunately, this simple modification progressively reduces the acceptance rate\. To further quantify the quality of the samples about different methods, we evaluated two quantities in Fig\.[6](https://arxiv.org/html/2609.27306#S4.F6)\(c\): the magnetization‑based ESS fraction and sample diversity\. The ESS fraction is1/\(1\+2​∑k=1Kρ⁡\(k\)\)1\\big/\\left\(1\+2\\sum\_\{k=1\}^\{K\}\\rho\(k\)\\right\), whereρ⁡\(k\)\\rho\(k\)denotes the autocorrelation of the averaged absolute‑magnetization series at lagkkandKKis the lag at whichρ⁡\(k\)\\rho\(k\)first becomes non‑positive\. Sample diversity measures the number of distinct spin configurations among the samples\. All three methods are run for identical durations, and we compare the diversity of an equal number of samples collected within this fixed runtime\. Since NMCMC and LMMC generate more total samples in the same runtime, we randomly take a subset of their samples so that the sample count used for comparison equals that of DDM\. For a fair comparison, each method is run on the faster of our available CPU and GPU: DDM and NMCMC on a single NVIDIA GeForce RTX 4090 GPU \(24 GB, CUDA 12\.2, PyTorch 2\.8\), and LMMC on the CPU \(AMD Ryzen 9 7950X3D, 16 cores\)\.

In general, our methods explore a substantially broader configuration space than the single\-spin\-flip update of LMMC while attaining higher acceptance rates\. This generates more diverse samples and reduces the total number of iterations required for convergence, at the cost of increased computational time per iteration\. Even so, non\-local sampling strategies of this kind exhibit promising prospects for efficient equilibrium sampling\.

Lastly, we outline three key limitations and practical constraints of our framework\. First, the efficiency and reliability of the sampling depend heavily on how well the neural networks are trained, which may fail to capture low‑probability configurations that nevertheless bear physical significance\. Second, our framework brings higher computational costs than standard spin\-flip update algorithms, since generating candidate samples via the discrete diffusion model and computing their likelihood increases the cost of every iteration, although faster chain mixing can partially compensate for this extra overhead\. Third, it remains challenging to establish strict theoretical bounds for the entire MCMC process\.

## VDiscussion

In this work, we have constructed a new type of discrete diffusion model in which the normalized distributions at successive diffusion times are parameterized by an evolving VAN\. Compared with prior discrete diffusion models based on tensor networks\[[10](https://arxiv.org/html/2609.27306#bib.bib5)\], our framework shows favorable performance on high\-dimensional lattice systems and alleviates the dimensionality bottleneck commonly encountered by MPS\. We incorporate the tau\-leaping algorithm to accelerate the sequential denoising process\. Relative to the step\-by\-step denoising baseline, the tau‑leaping‑accelerated method reduces iterative computational overhead while maintaining acceptable numerical precision for most sampling scenarios, achieving a better trade\-off between computational cost and sampling accuracy\. Comprehensive numerical experiments on both two\- and three\-dimensional Ising models verify that our method can accurately capture thermodynamic observables including free energy per site, energy per site, and magnetization per site across diverse temperatures, delivering competitive or even better numerical accuracy than tensor\-network baselines, especially for challenging three\-dimensional spin systems\.

We have also shown how to integrate the discrete diffusion model with MCMC to establish an efficient sampling paradigm for discrete spin systems\. Our unified pipeline supports two complementary sampling strategies, namely connected update and disconnected update, which together alleviate the local trapping problem that plagues conventional LMMC methods\. Specifically, the connected update constructs correlated candidate configurations through sequential forward noising and reverse denoising trajectories\. By contrast, the disconnected update performed via NMCMC generates fully independent proposals through VAN sampling\. Interestingly, the connected update bears close connections to both LMMC and the disconnected update\. Within the connected\-update framework, the forward noising process can be interpreted as a sequence of single\-spin flips in LMMC, while the reverse denoising process acts as a weakened variant of the disconnected update and recovers the full disconnected update in the limitt→∞t\\rightarrow\\infty\.

Finally, we discuss potential extensions of our work for future research\. First, the current framework can be extended to more complex spin systems to verify the generalization ability of our diffusion model beyond classical ferromagnetic Ising models\. Second, in addition to the tau\-leaping adopted in this paper, advanced acceleration techniques such as diffusion model distillation and enhanced tau\-leaping schemes\[[36](https://arxiv.org/html/2609.27306#bib.bib15),[55](https://arxiv.org/html/2609.27306#bib.bib56)\]can be incorporated to accelerate the denoising process\. Third, efficient numerical simulation of such stochastic dynamical processes can be extended to advance quantum generative models\[[13](https://arxiv.org/html/2609.27306#bib.bib59)\]\. Overall, our framework explicitly captures the full distribution dynamics of diffusion, which suggests preliminary directions for subsequent work on model interpretability and architectural optimization\.

## Acknowledgments

We acknowledge Online Club Nanothermodynamica for helpful discussions\. This work is supported by Project 12322501, 12575035 of National Natural Science Foundation of China, and 2026NSFSCZY0124 of the Natural Science Foundation of Sichuan Province\. The HPC is supported by the Center for HPC at University of Electronic Science and Technology of China\.

## Data availability

The authors declare that the data supporting this study are available within the paper\. A PyTorch code implementation of the present algorithm is openly available on GitHub \[[GitHub Repository](https://github.com/cowenp/Discrete-Diffusion-Models-via-Evolving-Variational-Autoregressive-Networks.git)\]\.

## Appendix ANotation

The key mathematical notations employed throughout this work are defined in Table[1](https://arxiv.org/html/2609.27306#A1.T1)\.

Table 1:Explanations of mathematical symbols
## Appendix BMasked autoencoder for distribution estimation

![Refer to caption](https://arxiv.org/html/2609.27306v1/fig6.png)Figure 7:Feasibility validation of the tau\-leaping method for accelerating the denoising process\.\(a\) Results for a𝒟=16×16\\mathcal\{D\}=16\\times 162D Ising model: we investigate the effects of different tau leaping steps with a fixed total simulation length of 50 denoising steps\. Left y\-axis shows mean energy per site after denoising as a function of jump step, right y\-axis shows mean magnetization per site; solid curves are noised\-denoised samples using the tau‑leaping procedure \(30 repeats with 200 samples\), and dashed horizontal lines indicate baseline values computed from initial samples\. \(b\) Analogous experiment on a 4×4×4 3D Ising model\. Together the panels illustrate the behavior of tau\-leaping in our discrete diffusion model, and we select the jump step where the solid and dashed curves show the closest agreement as the optimal leap duration\.In this work, we adopt the masked autoencoder for distribution estimation \(MADE\)\[[18](https://arxiv.org/html/2609.27306#bib.bib36)\]to model two\- and three\-dimensional Ising models\. For a standard L\-layer multi\-layer perceptron \(MLP\), the forward propagation formulation is defined as:

𝒉l\+1=σ⁡\(𝑾l​𝒉l\+𝒃l\)\\bm\{h\}^\{l\+1\}=\\sigma\\left\(\\bm\{W\}^\{l\}\\bm\{h\}^\{l\}\+\\bm\{b\}^\{l\}\\right\)\(21\)where𝒉0=𝒙\\bm\{h\}^\{0\}=\\bm\{x\}refers to the flattened lattice spin vector,𝒉L\\bm\{h\}^\{L\}denotes the final output of the network, andσ⁡\(⋅\)\\sigma\(\\cdot\)is the nonlinear activation function\.

To impose autoregressive constraints on the spin lattice configurations, fixed block\-wise binary mask matrices are applied to the weight matrices of all linear layers via element\-wise Hadamard products, which blocks the redundant connection mappings between individual spins\. The forward propagation of the masked layers is then formulated as

𝒉l\+1=\(𝑾l⊙𝑴l\)​σ​\(𝒉l\)\+𝒃l,\\bm\{h\}^\{l\+1\}=\\left\(\\bm\{W\}^\{l\}\\odot\\bm\{M\}^\{l\}\\right\)\\sigma\\\!\\left\(\\bm\{h\}^\{l\}\\right\)\+\\bm\{b\}^\{l\},\(22\)where⊙\\odotdenotes the element\-wise Hadamard product and𝑴l\\bm\{M\}^\{l\}is a fixed binary mask of the same shape as𝑾l\\bm\{W\}^\{l\}\. Two categories of triangular masks are constructed for different network layers: the input projection layer adopts an exclusive strictly lower\-triangular mask,𝑻i​j=1\\bm\{T\}\_\{ij\}=1fori\>ji\>j, which excludes the diagonal and eliminates the self\-dependency of each spin; the remaining masked layers utilize an inclusive lower\-triangular mask,𝑻i​j=1\\bm\{T\}\_\{ij\}=1fori≥ji\\geq j, which keeps the diagonal and preserves a valid sequential autoregressive dependency order\. Since both patterns are lower\-triangular in the site index, no spin can receive information from a spin of larger index, which is consistent with the autoregressive factorization\.

Regarding network hyper‑parameters, we configure the depth and width conditioned on the lattice dimensions\. We setnet\_depth=3\\text\{net\\\_depth\}=3andnet\_width=32\\text\{net\\\_width\}=32for the two‑dimensional Ising model, whilenet\_depth=4\\text\{net\\\_depth\}=4andnet\_width=16\\text\{net\\\_width\}=16are adopted for the three‑dimensional Ising model\. Our VAN is compatible with Adam optimizer\[[28](https://arxiv.org/html/2609.27306#bib.bib55)\]and natural gradient descent\[[31](https://arxiv.org/html/2609.27306#bib.bib24)\], which can be incorporated into the parameter update phase to effectively accelerate convergence during model training\. In practical training configurations, Adam generally requires a small learning rate, typically set to10−410^\{\-4\}or10−510^\{\-5\}, while natural gradient descent allows a relatively higher learning rate of10−210^\{\-2\}\.

## Appendix COptimal leap duration for tau\-leaping

This appendix explains the principle for determining the optimal leap duration used in the tau\-leaping Alg\.[2](https://arxiv.org/html/2609.27306#alg2), which controls the trade\-off between numerical precision and computational efficiency in stochastic simulations of the denoising process\. In our discrete diffusion model, we discretize the continuous‑time Markovian dynamics of the forward noising process with fixed‑size jump intervals and adopt a VAN to learn the probability distribution at each time step, which means that the exact transition rates of the denoising process are kept constant within each individual jump interval\. In Fig\.[7](https://arxiv.org/html/2609.27306#A2.F7), we illustrate the procedure to identify the optimal jump\-step size of tau\-leaping applied to the denoising process of two\- and three\-dimensional Ising models\. Since the noising process relies on Euler approximation fitting, a greater deviation between the adopted valueτ\\tauand the optimal leap duration will lead to lower denoising accuracy\. For scenarios that require high precision, it can be difficult to select the appropriateτ\\tauto keep the final results within an acceptable accuracy range\. When we do not find a suitableτ\\tauto keep numerical errors within acceptable limits, Euler\-discretized denoising serves as a reliable alternative\.

## Appendix DSampling algorithm

This appendix presents the full pseudocodes of multiple baseline and reference sampling algorithms discussed in the main text, including the standard local Metropolis Monte Carlo \(LMMC\), the neural Markov chain Monte Carlo \(NMCMC\), and the Wolff cluster algorithm\. The Wolff cluster algorithm, proposed in 1989\[[51](https://arxiv.org/html/2609.27306#bib.bib30)\], is a highly efficient single\-cluster Monte Carlo method\. Its single\-cluster update distinguishes it from the closely related Swendsen–Wang \(SW\) algorithm\[[45](https://arxiv.org/html/2609.27306#bib.bib31)\]\. The primary advantage of the Wolff algorithm is its ability to suppress critical slowing down near second\-order phase transitions\. In traditional algorithms such as Metropolis sampling, autocorrelation times increase rapidly near phase transitions, drastically reducing sampling efficiency\. Samples generated by the Wolff cluster algorithm are used to compute mean energy and magnetization per site, which serve as baseline values for comparison in this work\.

Despite its efficiency for critical sampling, the Wolff cluster algorithm exhibits inherent limitations that restrict its general applicability across diverse spin systems and temperature regimes\. First, its underlying cluster\-building principle is theoretically tailored to ferromagnetic systems with uniform short\-range couplings\. Second, the cluster construction strategy intrinsically relies on specific symmetry properties and structured spin interactions, which further limits its adaptability to generic disordered or asymmetric spin systems\.

Algorithm 3Wolff Cluster AlgorithmStep 1\.Initialize the cluster and stack:

- •Uniformly select a random lattice siteii\.
- •Initialize the cluster set:𝒞=\{i\}\\mathcal\{C\}=\\\{i\\\}\.
- •Create an empty stack and push siteiionto the stack\.

Step 2\.Grow the cluster using the stack:

- •While the stack is not empty, pop a sitejjfrom the stack\.
- •Traverse each nearest neighbork∈Ujk\\in U\_\{j\}of sitejj\.
- •For each neighborkksatisfying𝝈k=𝝈j\\bm\{\\sigma\}\_\{k\}=\\bm\{\\sigma\}\_\{j\}andk∉𝒞k\\notin\\mathcal\{C\}: - –Addkkto the cluster𝒞\\mathcal\{C\}with probability pk=1−e−2​β,p\_\{k\}=1\-e^\{\-2\\beta\}, - –Push the successfully added sitekkinto the stack\.

Step 3\.Flip all spins in the cluster:

𝝈j→−𝝈j∀j∈𝒞\\bm\{\\sigma\}\_\{j\}\\rightarrow\-\\bm\{\\sigma\}\_\{j\}\\quad\\forall j\\in\\mathcal\{C\}

Algorithm 4Standard local Metropolis Monte Carlo \(LMMC\)Step 1\.Propose a new configuration by flipping the spin of a randomly chosen site:

𝝈i→−𝝈i;\\bm\{\\sigma\}\_\{i\}\\rightarrow\-\\bm\{\\sigma\}\_\{i\};
Step 2\.Calculate the energy difference between the new and current configurations:

Δ​E=E⁡\(𝝈′\)−E⁡\(𝝈\)\\Delta E=E\(\\bm\{\\sigma\}^\{\\prime\}\)\-E\(\\bm\{\\sigma\}\)
Step 3\.Accept the move with probability:

Acc\[𝝈→𝝈′\]=min\[1,e−β​Δ​E\],\\text\{Acc\}\[\\bm\{\\sigma\}\\rightarrow\\bm\{\\sigma\}^\{\\prime\}\]=\\min\\left\[1,e^\{\-\\beta\\Delta E\}\\right\],
Accept𝝈′\\bm\{\\sigma\}^\{\\prime\}as the new configuration𝝈⁡\(t\+1\)=𝝈′\\bm\{\\sigma\}\(t\+1\)=\\bm\{\\sigma\}^\{\\prime\}if the transition is accepted; otherwise, keep the original state𝝈⁡\(t\+1\)=𝝈\\bm\{\\sigma\}\(t\+1\)=\\bm\{\\sigma\}for the next iteration\.

Algorithm 5Neural Markov Chain Monte Carlo \(NMCMC\)Step 1\.Initialize the lattice spin configuration𝝈\\bm\{\\sigma\}and adopt a fixed pre\-trained neural generatorP0θ​\(⋅\)P\_\{0\}^\{\\theta\}\(\\cdot\)\.

Step 2\.Propose a fully independent candidate spin configuration𝝈′\\bm\{\\sigma\}^\{\\prime\}directly sampled from the neural generative distribution:

𝝈′∼P0θ​\(𝝈′\)\.\\bm\{\\sigma\}^\{\\prime\}\\sim P\_\{0\}^\{\\theta\}\(\\bm\{\\sigma\}^\{\\prime\}\)\.
Step 3\.Evaluate the Metropolis–Hastings acceptance probability based on the neural prior distribution and Boltzmann weight to satisfy detailed balance:

Acc\[𝝈→𝝈′\]=min\[1,P0θ​\(𝝈\)⋅e−β​E​\(𝝈′\)P0θ​\(𝝈′\)⋅e−β​E​\(𝝈\)\]\.\\text\{Acc\}\[\\bm\{\\sigma\}\\to\\bm\{\\sigma\}^\{\\prime\}\]=\\min\\left\[1,\\,\\frac\{P\_\{0\}^\{\\theta\}\(\\bm\{\\sigma\}\)\\cdot e^\{\-\\beta E\(\\bm\{\\sigma\}^\{\\prime\}\)\}\}\{P\_\{0\}^\{\\theta\}\(\\bm\{\\sigma\}^\{\\prime\}\)\\cdot e^\{\-\\beta E\(\\bm\{\\sigma\}\)\}\}\\right\]\.
Accept𝝈′\\bm\{\\sigma\}^\{\\prime\}as the new configuration𝝈⁡\(t\+1\)=𝝈′\\bm\{\\sigma\}\(t\+1\)=\\bm\{\\sigma\}^\{\\prime\}if the transition is accepted; otherwise, keep the original state𝝈⁡\(t\+1\)=𝝈\\bm\{\\sigma\}\(t\+1\)=\\bm\{\\sigma\}for the next iteration\.

## References

- \[1\]\(2025\)Block diffusion: interpolating between autoregressive and diffusion language models\.InInternational Conference on Learning Representations,Vol\.2025,pp\. 50726–50753\.External Links:[Link](https://proceedings.iclr.cc/paper_files/paper/2025/hash/7ede97c3e082c6df10a8d6103a2eebd2-Abstract-Conference.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[2\]J\. Austin, D\. D\. Johnson, J\. Ho, D\. Tarlow, and R\. Van Den Berg\(2021\)Structured denoising diffusion models in discrete state\-spaces\.Adv\. Neural Inf\. Process\. Syst\.34,pp\. 17981–17993\.External Links:[Link](https://proceedings.neurips.cc/paper/2021/hash/6f558439a07a94918a07773f523f87a-Abstract.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[3\]Y\. Bahri, J\. Kadmon, J\. Pennington, S\. S\. Schoenholz, J\. Sohl\-Dickstein, and S\. Ganguli\(2020\)Statistical mechanics of deep learning\.Annu\. Rev\. Condens\. Matt\. Phys\.11\(1\),pp\. 501–528\.External Links:[Link](https://www.annualreviews.org/content/journals/10.1146/annurev-conmatphys-031119-050745)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[4\]G\. Batzolis, J\. Stanczuk, C\. Schönlieb, and C\. Etmann\(2021\)Conditional image generation with score\-based diffusion models\.arXiv:2111\.13606\.External Links:[Link](https://arxiv.org/abs/2111.13606)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[5\]I\. Biazzo, D\. Wu, and G\. Carleo\(2024\)Sparse autoregressive neural networks for classical spin systems\.Mach\. Learn\.: Sci\. Technol\.5\(2\),pp\. 025074\.External Links:[Link](https://iopscience.iop.org/article/10.1088/2632-2153/ad5783)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p5.1)\.
- \[6\]A\. Campbell, J\. Benton, V\. De Bortoli, T\. Rainforth, G\. Deligiannidis, and A\. Doucet\(2022\)A continuous time framework for discrete denoising models\.Adv\. Neural Inf\. Process\. Syst\.35,pp\. 28266–28279\.External Links:[Link](https://papers.nips.cc/paper_files/paper/2022/hash/b5b528767aa35f5b1a60fe0aaeca0563-Abstract-Conference.html)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p4.1)\.
- \[7\]Y\. Cao, D\. T\. Gillespie, and L\. R\. Petzold\(2005\)Avoiding negative populations in explicit poisson tau\-leaping\.J\. Chem\. Phys\.123\(5\)\.External Links:[Link](https://pubs.aip.org/aip/jcp/article-abstract/123/5/054104/905859/Avoiding-negative-populations-in-explicit-Poisson?redirectedFrom=fulltext)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p4.1)\.
- \[8\]Y\. Cao, D\. T\. Gillespie, and L\. R\. Petzold\(2006\)Efficient step size selection for the tau\-leaping simulation method\.J\. Chem\. Phys\.124\(4\)\.External Links:[Link](https://pubs.aip.org/aip/jcp/article-abstract/124/4/044109/562210/Efficient-step-size-selection-for-the-tau-leaping?redirectedFrom=fulltext)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p4.1)\.
- \[9\]G\. Carleo and M\. Troyer\(2017\)Solving the quantum many\-body problem with artificial neural networks\.Science355\(6325\),pp\. 602–606\.External Links:[Link](https://www.science.org/doi/10.1126/science.aag2302)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[10\]L\. Causer, G\. M\. Rotskoff, and J\. P\. Garrahan\(2025\)Discrete generative diffusion models without stochastic differential equations: a tensor network approach\.Phys\. Rev\. E111\(2\),pp\. 025302\.External Links:[Link](https://doi.org/10.1103/PhysRevE.111.025302)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p3.1),[§I](https://arxiv.org/html/2609.27306#S1.p4.1),[Figure 3](https://arxiv.org/html/2609.27306#S3.F3),[§III](https://arxiv.org/html/2609.27306#S3.p2.1),[§IV](https://arxiv.org/html/2609.27306#S4.p8.1),[§V](https://arxiv.org/html/2609.27306#S5.p1.1)\.
- \[11\]D\. Chandler\(1987\)Introduction to modern statistical mechanics\.Oxford University Press,Oxford, UK\.Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[12\]H\. Chen, S\. Liu, and J\. Yang\(2026\)Markov chain monte carlo with diffusion paths\.arXiv:2607\.11631\.External Links:[Link](https://arxiv.org/abs/2607.11631)Cited by:[§IV](https://arxiv.org/html/2609.27306#S4.p1.1),[§IV](https://arxiv.org/html/2609.27306#S4.p7.1)\.
- \[13\]Z\. Cui, P\. Zhang, and Y\. Tang\(2026\)Quantum flow matching\.PRX Intelligence1,pp\. 013009\.External Links:[Document](https://dx.doi.org/10.1103/s6zj-vzdp),[Link](https://link.aps.org/doi/10.1103/s6zj-vzdp)Cited by:[§V](https://arxiv.org/html/2609.27306#S5.p3.1)\.
- \[14\]V\. De Bortoli, J\. Thornton, J\. Heng, and A\. Doucet\(2021\)Diffusion schrödinger bridge with applications to score\-based generative modeling\.Adv\. Neural Inf\. Process\. Syst\.34,pp\. 17695–17709\.External Links:[Link](https://proceedings.neurips.cc/paper_files/paper/2021/file/940392f5f32a7ade1cc201767cf83e31-Paper.pdf)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[15\]L\. de Groot, R\. J\. A\. Kuiper, and A\. BagheriA survey of discrete diffusion for text and genomic sequence generation\.InThe 37th Benelux Conference on Artificial Intelligence and the 34th Belgian Dutch Conference on Machine Learning,External Links:[Link](https://openreview.net/forum?id=ubOzVt7oMw)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[16\]L\. M\. Del Bono, F\. Ricci\-Tersenghi, and F\. Zamponi\(2025\)Performance of machine\-learning\-assisted monte carlo in sampling from simple statistical physics models\.Phys\. Rev\. E112\(4\),pp\. 045307\.External Links:[Link](https://doi.org/10.1103/PhysRevE.112.045307)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[17\]L\. M\. Del Bono, F\. Ricci\-Tersenghi, and F\. Zamponi\(2026\)Demonstrating real advantage of machine learning–enhanced monte carlo for combinatorial optimization\.Proc\. Natl\. Acad\. Sci\.123\(19\),pp\. e2534768123\.External Links:[Link](https://www.pnas.org/doi/10.1073/pnas.2534768123)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[18\]M\. Germain, K\. Gregor, I\. Murray, and H\. Larochelle\(2015\)Made: masked autoencoder for distribution estimation\.InInternational conference on machine learning,pp\. 881–889\.External Links:[Link](https://proceedings.mlr.press/v37/germain15.html)Cited by:[Appendix B](https://arxiv.org/html/2609.27306#A2.p1.1)\.
- \[19\]D\. T\. Gillespie\(1976\)A general method for numerically simulating the stochastic time evolution of coupled chemical reactions\.J\. Comput\. Phys\.22\(4\),pp\. 403–434\.External Links:[Link](http://web.mit.edu/endy/www/scraps/signal/JCompPhys%2822%29403.pdf)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p4.1)\.
- \[20\]D\. T\. Gillespie\(1977\)Exact stochastic simulation of coupled chemical reactions\.J\. Phys\. Chem\.81\(25\),pp\. 2340–2361\.External Links:[Link](http://web.mit.edu/endy/www/scraps/dg/JPC%2881%292340.pdf)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p4.1)\.
- \[21\]D\. T\. Gillespie\(2001\)Approximate accelerated stochastic simulation of chemically reacting systems\.J\. Chem\. Phys\.115\(4\),pp\. 1716–1733\.External Links:[Link](https://pubs.aip.org/aip/jcp/article-abstract/115/4/1716/451187/Approximate-accelerated-stochastic-simulation-of?redirectedFrom=fulltext)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p4.1)\.
- \[22\]N\. Gruver, S\. Stanton, N\. Frey, T\. G\. Rudner, I\. Hotzel, J\. Lafrance\-Vanasse, A\. Rajpal, K\. Cho, and A\. G\. Wilson\(2023\)Protein design with guided discrete diffusion\.Adv\. Neural Inf\. Process\. Syst\.36,pp\. 12489–12517\.External Links:[Link](https://proceedings.neurips.cc/paper_files/paper/2023/hash/29591f355702c3f4436991335784b503-Abstract-Conference.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[23\]J\. Ho, A\. Jain, and P\. Abbeel\(2020\)Denoising diffusion probabilistic models\.Adv\. Neural Inf\. Process\. Syst\.33,pp\. 6840–6851\.External Links:[Link](https://proceedings.neurips.cc/paper/2020/hash/4c5bcfec8584af0d967f1ab10179ca4b-Abstract.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[24\]J\. Ho, C\. Saharia, W\. Chan, D\. J\. Fleet, M\. Norouzi, and T\. Salimans\(2022\)Cascaded diffusion models for high fidelity image generation\.J\. Mach\. Learn\. Res\.23\(47\),pp\. 1–33\.External Links:[Link](https://jmlr.org/papers/volume23/21-0635/21-0635.pdf)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[25\]N\. T\. Hunt\-Smith, W\. Melnitchouk, F\. Ringer, N\. Sato, A\. W\. Thomas, and M\. J\. White\(2024\)Accelerating markov chain monte carlo sampling with diffusion models\.Comput\. Phys\. Commun\.296,pp\. 109059\.External Links:[Link](https://www.sciencedirect.com/science/article/pii/S0010465523004046?via%3Dihub)Cited by:[§IV](https://arxiv.org/html/2609.27306#S4.p1.1)\.
- \[26\]T\. Kathuria and S\. Kumar\(2026\)Leveraging pretrained language models as energy functions for glauber dynamics text diffusion\.arXiv:2605\.04291\.External Links:[Link](https://arxiv.org/abs/2605.04291)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[27\]N\. Ketkar and J\. Moolayil\(2021\)Convolutional neural networks\.InDeep learning with Python: learn best practices of deep learning models with PyTorch,pp\. 197–242\.External Links:[Link](https://link.springer.com/book/10.1007/978-1-4842-5364-9)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p4.1)\.
- \[28\]D\. P\. Kingma and J\. Ba\(2014\)Adam: a method for stochastic optimization\.arXiv:1412\.6980\.External Links:[Link](https://arxiv.org/abs/1412.6980)Cited by:[Appendix B](https://arxiv.org/html/2609.27306#A2.p3.1)\.
- \[29\]J\. Lemercier, J\. Richter, S\. Welker, E\. Moliner, V\. Välimäki, and T\. Gerkmann\(2025\)Diffusion models for audio restoration: a review\.IEEE Signal Processing Magazine41\(6\),pp\. 72–84\.External Links:[Link](https://ieeexplore.ieee.org/document/10819705)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[30\]M\. Levin and C\. P\. Nave\(2007\)Tensor renormalization group approach to two\-dimensional classical lattice models\.Phys\. Rev\. Lett\.99\(12\),pp\. 120601\.External Links:[Link](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.120601)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p4.1)\.
- \[31\]J\. Liu, Y\. Tang, and P\. Zhang\(2025\)Efficient optimization of variational autoregressive networks with natural gradient\.Phys\. Rev\. E111\(2\),pp\. 025304\.External Links:[Link](https://arxiv.org/abs/2409.20029)Cited by:[Appendix B](https://arxiv.org/html/2609.27306#A2.p3.1)\.
- \[32\]A\. Lou, C\. Meng, and S\. Ermon\(2023\)Discrete diffusion modeling by estimating the ratios of the data distribution\.arXiv:2310\.16834\.External Links:[Link](https://arxiv.org/abs/2310.16834)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[33\]D\. Luo, Z\. Chen, J\. Carrasquilla, and B\. K\. Clark\(2022\)Autoregressive neural network for simulating open quantum systems via a probabilistic formulation\.Phys\. Rev\. Lett\.128\(9\),pp\. 090501\.External Links:[Link](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.128.090501)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[34\]B\. McNaughton, M\. Milošević, A\. Perali, and S\. Pilati\(2020\)Boosting monte carlo simulations of spin glasses using autoregressive neural networks\.Phys\. Rev\. E101\(5\),pp\. 053312\.External Links:[Link](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.101.053312)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p5.1)\.
- \[35\]P\. Mehta, M\. Bukov, C\. Wang, A\. G\. Day, C\. Richardson, C\. K\. Fisher, and D\. J\. Schwab\(2019\)A high\-bias, low\-variance introduction to machine learning for physicists\.Phys\. Rep\.810,pp\. 1–124\.External Links:[Link](https://www.sciencedirect.com/science/article/pii/S0370157319300766)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p5.1)\.
- \[36\]Y\. Ren, H\. Chen, Y\. Zhu, W\. Guo, Y\. Chen, G\. M\. Rotskoff, M\. Tao, and L\. Ying\(2025\)Fast solvers for discrete diffusion models: theory and applications of high\-order algorithms\.arXiv:2502\.00234\.External Links:[Link](https://arxiv.org/abs/2502.00234)Cited by:[§V](https://arxiv.org/html/2609.27306#S5.p3.1)\.
- \[37\]S\. Sanokowski, W\. Berghammer, H\. Wang, M\. Ennemoser, S\. Hochreiter, and S\. Lehner\(2025\)Scalable discrete diffusion samplers: combinatorial optimization and statistical physics\.InInternational Conference on Learning Representations,Vol\.2025,pp\. 87053–87082\.External Links:[Link](https://proceedings.iclr.cc/paper_files/paper/2025/hash/d87eadbe6114b768e76c5d5e8eb39388-Abstract-Conference.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[38\]A\. Sarkar, Z\. Tang, C\. Zhao, and P\. K\. Koo\(2024\)Designing dna with tunable regulatory activity using discrete diffusion\.bioRxiv,pp\. 2024–05\.External Links:[Link](https://www.biorxiv.org/content/10.1101/2024.05.23.595630v1)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[39\]F\. Schneider\(2023\)Archisound: audio generation with diffusion\.arXiv:2301\.13267\.External Links:[Link](https://arxiv.org/abs/2301.13267)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[40\]O\. Sharir, Y\. Levine, N\. Wies, G\. Carleo, and A\. Shashua\(2020\)Deep autoregressive models for the efficient variational simulation of many\-body quantum systems\.Phys\. Rev\. Lett\.124\(2\),pp\. 020503\.External Links:[Link](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.124.020503)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[41\]J\. Shin, A\. J\. Riesselman, A\. W\. Kollasch, C\. McMahon, E\. Simon, C\. Sander, A\. Manglik, A\. C\. Kruse, and D\. S\. Marks\(2021\)Protein design and variant prediction using autoregressive generative models\.Nat\. Commun\.12\(1\),pp\. 2403\.External Links:[Link](https://www.nature.com/articles/s41467-021-22732-w)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[42\]N\. S\. Singh, M\. Mohseni\-Rajaee, S\. Niazi, and K\. Y\. Camsari\(2026\)From independent to correlated diffusion: generalized generative modeling with probabilistic computers\.arXiv:2603\.27996\.External Links:[Link](https://arxiv.org/abs/2603.27996)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[43\]J\. Sohl\-Dickstein, E\. Weiss, N\. Maheswaranathan, and S\. Ganguli\(2015\)Deep unsupervised learning using nonequilibrium thermodynamics\.InInternational conference on machine learning,pp\. 2256–2265\.External Links:[Link](https://proceedings.mlr.press/v37/sohl-dickstein15.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[44\]H\. Sun, L\. Yu, B\. Dai, D\. Schuurmans, and H\. Dai\(2022\)Score\-based continuous\-time discrete diffusion models\.arXiv:2211\.16750\.External Links:[Link](https://arxiv.org/abs/2211.16750)Cited by:[§II\.3](https://arxiv.org/html/2609.27306#S2.SS3.p1.1)\.
- \[45\]R\. H\. Swendsen and J\. Wang\(1987\)Nonuniversal critical dynamics in monte carlo simulations\.Phys\. Rev\. Lett\.58\(2\),pp\. 86\.External Links:[Link](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.58.86)Cited by:[Appendix D](https://arxiv.org/html/2609.27306#A4.p1.1)\.
- \[46\]Y\. Tang, J\. Liu, J\. Zhang, and P\. Zhang\(2024\)Learning nonequilibrium statistical mechanics and dynamical phase transitions\.Nat\. Commun\.15\(1\),pp\. 1117\.External Links:[Link](https://www.nature.com/articles/s41467-024-45172-8)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p5.1)\.
- \[47\]Y\. Tang, J\. Weng, and P\. Zhang\(2023\)Neural\-network solutions to stochastic reaction networks\.Nat\. Mach\. Intell\.5\(4\),pp\. 376–385\.External Links:[Link](https://www.nature.com/articles/s42256-023-00632-6)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[48\]C\. Wang, M\. Uehara, Y\. He, A\. Wang, A\. Lal, T\. Jaakkola, S\. Levine, A\. Regev, H\. Wang, and T\. Biancalani\(2025\)Fine\-tuning discrete diffusion models via reward optimization with applications to dna and protein design\.InInternational Conference on Learning Representations,Vol\.2025,pp\. 47871–47899\.External Links:[Link](https://proceedings.iclr.cc/paper_files/paper/2025/hash/771e09dd204ea339da0d8114c48afd21-Abstract-Conference.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[49\]D\. Wang, R\. Qiu, and Z\. Huang\(2026\)When to commit? towards variable\-size self\-contained blocks for discrete diffusion language models\.arXiv:2604\.23994\.External Links:[Link](https://arxiv.org/abs/2604.23994)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[50\]J\. Weng, X\. Zhu, J\. Liu, L\. Lü, P\. Zhang, and Y\. Tang\(2025\)Tracking large chemical reaction networks and rare events by neural networks\.arXiv:2512\.10309\.External Links:[Link](https://arxiv.org/abs/2512.10309)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[51\]U\. Wolff\(1989\)Collective monte carlo updating for spin systems\.Phys\. Rev\. Lett\.62\(4\),pp\. 361\.External Links:[Link](https://link.aps.org/doi/10.1103/PhysRevLett.62.361)Cited by:[Appendix D](https://arxiv.org/html/2609.27306#A4.p1.1)\.
- \[52\]D\. Wu, L\. Wang, and P\. Zhang\(2019\)Solving statistical mechanics using variational autoregressive networks\.Phys\. Rev\. Lett\.122\(8\),pp\. 080602\.External Links:[Link](https://doi.org/10.1103/PhysRevLett.122.080602)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1)\.
- \[53\]Z\. Xie, H\. Jiang, Q\. N\. Chen, Z\. Weng, and T\. Xiang\(2009\)Second renormalization of tensor\-network states\.Phys\. Rev\. Lett\.103\(16\),pp\. 160601\.External Links:[Link](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.103.160601)Cited by:[§II\.1](https://arxiv.org/html/2609.27306#S2.SS1.p4.1)\.
- \[54\]L\. Yang, Z\. Zhang, Y\. Song, S\. Hong, R\. Xu, Y\. Zhao, W\. Zhang, B\. Cui, and M\. Yang\(2023\)Diffusion models: a comprehensive survey of methods and applications\.ACM Comput\. Surv\.56\(4\),pp\. 1–39\.External Links:[Link](https://dl.acm.org/doi/pdf/10.1145/3626235)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p1.1)\.
- \[55\]Y\. Yao, H\. Zhou, A\. Han, W\. Huang, and M\. Sugiyama\(2026\)Accelerating discrete diffusion models with parallel\-in\-time sampling\.arXiv:2607\.00773\.External Links:[Link](https://arxiv.org/abs/2607.00773)Cited by:[§V](https://arxiv.org/html/2609.27306#S5.p3.1)\.
- \[56\]R\. Yu, Q\. Li, and X\. Wang\(2025\)Discrete diffusion in large language and multimodal models: a survey\.arXiv:2506\.13759\.External Links:[Link](https://arxiv.org/abs/2506.13759)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[57\]L\. Zhong, W\. Duan, J\. Liu, P\. Zhang, and Y\. Tang\(2026\)Scalable physics\-inspired transformers for spin glasses\.arXiv:2606\.22984\.External Links:[Link](https://arxiv.org/abs/2606.22984)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p4.1),[§III](https://arxiv.org/html/2609.27306#S3.p3.1)\.
- \[58\]H\. Zhou, S\. Roy, and R\. Gangadharaiah\(2026\)Steering without breaking: mechanistically informed interventions for discrete diffusion language models\.arXiv:2605\.10971\.External Links:[Link](https://arxiv.org/abs/2605.10971)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.
- \[59\]Y\. Zhu, W\. Guo, J\. Choi, G\. Liu, Y\. Chen, and M\. Tao\(2025\)Mdns: masked diffusion neural sampler via stochastic optimal control\.Adv\. Neural Inf\. Process\. Syst\.38,pp\. 35260–35308\.External Links:[Link](https://papers.neurips.cc/paper_files/paper/2025/hash/3289cc9cb9fdde172004497c61359544-Abstract-Conference.html)Cited by:[§I](https://arxiv.org/html/2609.27306#S1.p2.1)\.