Fast Solvers for Discrete Diffusion Models: Theory and Applications of High-Order Algorithms
Authors: Yinuo Ren
Haoxuan Chen
Yuchen Zhu
Wei Guo
Yongxin Chen
Grant M. Rotskoff
Molei Tao
Lexing Ying
Venue: NeurIPS 2025 poster ()
Keywords: Discrete diffusion model
fast solver
high-order scheme
Runge-Kutta method
Primary area: deep_learning
TL;DR: We made discrete diffusion model inference faster with high-order methods, both theoretically and empirically.

Abstract: Discrete diffusion models have emerged as a powerful generative modeling framework for discrete data with successful applications spanning from text generation to image synthesis. However, their deployment faces challenges due to the high dimensionality of the state space, necessitating the development of efficient inference algorithms. Current inference approaches mainly fall into two categories: exact simulation and approximate methods such as $\tau$-leaping. While exact methods suffer from unpredictable inference time and redundant function evaluations, $\tau$-leaping is limited by its first-order accuracy. In this work, we advance the latter category by tailoring the first extension of high-order numerical inference schemes to discrete diffusion models, enabling larger step sizes while reducing error. We rigorously analyze the proposed schemes and establish the second-order accuracy of the $\theta$-Trapezoidal method in KL divergence. Empirical evaluations on GSM8K-level math-reasoning, GPT-2-level text, and ImageNet-level image generation tasks demonstrate that our method achieves superior sample quality compared to existing approaches under equivalent computational constraints, with consistent performance gains across models ranging from 200M to 8B. Our code is available at https://github.com/yuchen-zhu-zyc/DiscreteFastSolver

Fast Solvers for Discrete Diffusion Models: Theory and Applications of High-Order Algorithms
Yinuo Ren
Haoxuan Chen l.*i
Yuchen Zhu
Wei Guo Lexing Ying!
Yongxin Chen
Grant M. Rotskoff
Molei Tao
Stanford University Georgia Institute of Technology fyinuoren, haoxuanc, rotskoff, lexing) ord. edd (yzhu738, wei.guo yongchen, mtao)@gatech. edu
Abstract
Discrete diffusion models have emerged powerful generative modeling frame- work for discrete data with successful applications spanning from text generation to image synthesis. However, their deployment faces challenges due to the high dimensionality the state space, necessitating the development of efficient infer- ence algorithms. Current inference approaches mainly fall into two categories: exact simulation and approximate methods such as T-leaping. While exact meth- ods suffer from unpredictable inference time and redundant function evaluations, r-leaping is limited by its first-order accuracy. In this work, we advance the latter category by tailoring the first extension of high-order numerical inference schemes to discrete diffusion models, enabling larger step sizes while reducing error. We rigorously analyze the proposed schemes and establish the second-order accuracy of the 0- O-Trrpezooiall method in KL. divergence. Empirical evaluations on GSMSK- level math-reasoning, GPT-2-level text, and ImageNeetlevel image generation tasks demonstrate that our method achieves superior sample quality compared to existing approaches under equivalent computational constraints, with consistent performance gains across models ranging from 200M to 8B. Our code is available
Introduction
Diffusion and flow-based models on discrete spaces [I] [00 have emerged a cornerstone of modern generative modeling for categorical data, offering unique advantages in domains where continuity assumptions fail. Unlike their continuous counterparts, discrete diffusion models inherently accom- modate data with discrete structures, e.g., language tokens, molecular sequences, tokenized images, and graphs, enabling principled generation and inference combinatorially complex spaces. These models have exerted a large impact on numerous applications, from the design of molecules [[]]].
proteins [12], and DNA sequences |IK [4] under biophysical constraints, to generation of high-fidelityy text [TS] and images [[G]] via autoregressive or masked transitions, etc.. Beyond stan-
dalone tasks, discrete diffusion models also synergize with methodologies, ranging from tensor net-
works [| to guidance mechanisms [TX-20]
Discrete diffusion models, despite their broad applicability, face a critical bottleneck: inference inefficiency. Current inference methods include: (1) exact simulation methods [2I]], which ensure unbiased sampling from the pre-trained model but suffer from unpredictable inference time and redundant score evaluations, resulting in poor scaling w.r.t. dimensionality; and (2) approximate
"Equal contribution Corresponding author
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

methods such r-leaping [22], which offer simple and parallelizable implementation but, due to their first-order accuracy. requires small step sizes to control discretization error, forcing stringent trade-off between speed and sample quality.
To address these limitations in possibly computationally constrained environments, we develop high- order numerical schemes tailored for discrete diffusion model inference. Drawing inspirations from acceleration techniques developed for ordinary differential equations (ODEs) [23], stochastic dif-
ferential equations (SDEs) [244 25]. chemical reaction simulations [26]]. and most recently continu-
ous diffusion [27 291, our work represents the first successful adaptation of high-order numerical
schemes to the discrete diffusion domain. Through careful design, these high-order schemes provide unprecedented efficient and versatile solutions for discrete diffusion model inference.
Our Contributions.
The main contributions of this paper are summarized follows:
We introduce the firs high-order numerical solvers for discrete diffusion model inference, namely the B-Runge-Kutta-2 (0-RK-2) method and the 8-Trapezoidal method;
We rigorously establish the theoretical properties of both methods, proving second-order con- vergence of the d-Trapezoidal method and conditional second-order convergence of the A-RK-2 method;
We empirically validate our theoretical results and demonstrate the superior performance the Trapezoidal method through comprehensive evaluations on large-scale text and image generation benchmarks.
1.1 Related Works
Here we briefly review related works and defer more detailed discussion to App.
Discrete Diffusion Models. Since their introduction, discrete diffusion models have undergone significant refinements, including the development of score-entropy loss [30] and flow-matching formulation [3T]. 5- These models generally fall into two categories based on their noise distribu-
tion: uniform [M]. MI and masked (absorbing state) [<K-t5] ES. [l]], each offering unique advantages in modeling discrete distributions. Recent theoretical advances have emerged through studies [36- B8).
High-Order Scheme for Continuous Diffusion Models. The development of high-order numer- ical schemes for solving ODEs and SDEs represents decades of research, as comprehensively re- viewed in [23, 39, [M]. These schemes have recently been adapted accelerate continuous dif- fusion model inference, encompassing approaches such as the exponential integrators [H] [3], Cllam--Bashfoe methods [29] 9. A5], Taylor methods [27 [6] and (stochastic) Runge-Kutta meth- ods |47 28 43-511 SI]..
High-Order Scheme for Chemical Reaction Systems. Regarding approximate methods for sim-
ulating compound Poisson processes and chemical reaction systems with state- dettedde intensities.
27. 55, 56).. For quick review of the problem setting and these methods, one may refer to [57] 58 The adaptation of
efforts have been made on the r-leaping method [52]], and its extensions [53, 54
these methods to discrete diffusion models presents unique challenges due to the presence of both time and state-inhomogeneuus intensities in the underlying Poisson processes.
2 Preliminaries
In this subsection, we review several basic concepts and previous error analysis results of discrete diffusion models.
2.1 Discrete Diffusion Models
In discrete diffusion models, one considers continuous-time Markov chain (CTMC) (æ:)osisT finite space X as the forward process. We represent the distribution of æ; by vector Pt E A [x] where A X| denotes the probability simplex in [ Given a target distribution Po, the CTMC

---

satisfies the following equation:
dp:
Q.pr. where Q: (Q)( x))
(2.1)
is the rate matrix at time t satisfying
① Q((x)))
Σ Qe(y, Vor X; (ii) Qe v) 2 0, Vx t yęe
Below we will use the notation Qi = Q. diag Qt- It can be shown that the corresponding backward process is of the same form but with a different rate matrix [59]:
dp. - Q.P,: ds
(.( P.(a) a. (z,
where (v,
(2.2)
- -y Q,),
Yr y-=
the rate matrix and &, denotes *r--s- The rate matrix Q: is often chosen to possess certain sparse structures such that the forward process converges to simple distribution that is easy to sample from. Popular choices include the uniform and absorbing state cases [30], where the forward process (2.1l converges to the uniform distribution on X and Dirac distribution, respectively.
Common training practice is to define the score function (or the score vector) as 81(x) (#e ()),,,))xxxx) == pe[m]
Pt_. for any æ E X, t E [0, and estimate by a neural network 8 (x), where
the parameters & are trained by minimizing the score entropy [60] 60] for some weights ei 2 0 as:
min /
Σ Qi(Ee. v) (sc(at, y) log sx(x.v) se v) + 5 (It. (f(())))
dt. (2.3)
Similar to the continuous case, the backward process is approximated by another CTMC
dg. de
Q% qs. with qo Pro and rate matrix Q where Q, (y, x) = 5% (z. V)Q, y) for any x + y E X. The inference is done by first sampling from pso and then evolving the CTMC accordingly. For simplicity, we drop the superscript & hereafter.
2.2 Stochastic Integral Formulation of Discrete Diffusion Models
Discrete diffusion models can also be formulated as stochastic integrals, which is especially useful for their theoretical analysis [[8]. In this section, we briefly recapitulate the relevant results and refer to App. B for the mathematical details. Below, we work on the probability space (St, B, IP) and denote the pairwise difference set of the state space X by DD = {x - y # y E X]. In this work, we focus on the case where [s] with d data dimensions and sites along each dimension.
We first introduce the Poisson random measure, a key concept in the formulation.
Definition 2.1 (Informal Definition of Poisson Random Measure). The random measure N [[] dv) on R++ x D called Poisson random measure with evolving intensity w.t.t measure 7 ID if, roughly speaking, the number of jumps of magnitude during the infinitesimal time interval (t, I + df) is Poisson distributed with mean >e(v)r(du)er
The forward process (2. can thus represented by the following stochastic integral:
21 = xo 2 + 1 uNN[)(dss, dr),
where the intensity X is defined (e) = (1(x (w) vT (w)) if x. (w) X and 0 otherwise. Here, the outcome w E Ω and -y- denotes the left limit of the càdlàg process de at time i with zi- 3op We will also omit the variable w, should be clear from context. The backward process in discrete diffusion models (2-2) can also be represented similarly as:
V. = 900 S/¹ J [N
(2.4)
+ E X where the intensity μ is defined as (v, (v) = F. --- +v) Ms- + vf if and 0 otherwise. During inference, T. = Do + So (d] du) is used instead of (2.4). where the estimated intensity f is defined by replacing the true score with the neural network estimated score & in (v, w). the following, we also denote the intensity ps (v, w) time by μs (M (,,)))) with slight abuse of terminology to emphasize its dependency through ys- (w).
3

---

Inference Process
7-Leaping
Perplexity
3000
NFE intersilty
() Pan
R-RK-2
Perplexity
2000
(ii) (1 26 Mon + 28 Pi
36
(i) Pan.
0-Trapezoidal
500
(ii a a cxafan.
100
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Us
0
Pp
(1 (1-0)44 0)
Ven+
Figure I: Leftt: Application of the uniformization algorithm to discrete diffusion models for text gen- eration. The æ-axis denotes the time of the backward process, and the y-axis denotes the frequency of jumps (NFE). Perplexity convergence occurs before the NFE grows unbounded. Right: Compar- ison between r-leaping and the proposed second-order schemes (0-RK-2 and A- T-aapezoidal).
3 Numerical Schemes for Discrete Diffusion Model Inference
Before introducing the proposed numerical schemes, we first review existing numerical schemes for discrete diffusion models, including exact simulation methods and the r-leaping method, and discuss their merits and limitations.
3.1 Exact Simulation Methods
Unlike in continuous diffusion models, where exact simulation is infeasible, discrete diffusion mod- els permit inference without discretization error. Notable examples of unbiased samplers include uniformization [[56] for the uniform state case and the First-Hitting Sampler (FHS) [2I] for the ab- sorbing state case. The main idea behind these methods is first sample the next jump time and then the jump itself. Theoretical analysis [GX] reveals that such schemes lack guarantees with finite finite computation budget, since the number of required jumps (and thus the inference time) follows random distribution with expectation ((d). This computational restriction may be less favorable for high-dimensional applications, such as generative modeling of DNA or protein sequences.
Furthermore, the absence of discretization error does not necessarily translate to superior sample quality, given the inherent estimation errors in neural network-based score functions. This limitation further amplified by the highly skewed distribution of jumps, with concentration occurring during the terminal phase of the backward process, when the neural network- score function exhibits the highest estimation error. This phenomenon stems from the potential singularity of the target distribution Pa, which induces singularities in the score function, making accurate neural network estimation particularly challenging during that phase (ef. Assump. 4.4 [[8]).
The left figure in Fig. M illustrates application of the uniformization algorithm to discrete dif- fusion inference for text generation, with detailed experimental parameters presented in Sec. b.d and App. D As the process approaches the target distribution (t - 7"), the number of jumps (in terms of the number of score function evaluations, NFE) grows unbounded, while perplexity improvements become negligible. This skew in computational effort leads to redundant function evaluations. Although early stopping is commonly adopted at T 5 for some small & 1 to allevi- ate this inefficiency, this approach introduces challenges in its selection, particularly under compu- tational constraints or when efficiency-accuue -accuracy trade-offs are desired. Moreover, the variable jump schedules across batch samples complicate parallelization efforts in exact methods, highlighting the need for more adaptable and efficient algorithmic solutions.
3.2
Approximate Method: r-Leaping Method
The r-leaping method [52 22 is widely adopted scheme that effectively addresses both dimen- sionality scaling and inference-time control challenges. This Euler-type scheme approximates the backward process with dependepn intensity ft, via the following updates:
Dr+ Te Σ vP (Po(v)A),
(3.1)

---

where A denotes the time step and P(-) denotes a Poisson random variable. In general, one may design different discretization schemes for r-leaping, and the summation in (B.I) parallelizable. underscoring the method's s flexibility and efficiency. We refer to Alg. B and App. B for detailed description of the r-leaping method for discrete diffusion model inference. Regarding convergence properties as the time discretization becomes increasingly refined, theoretical analyses by [22, [8] have established the error bounds of the r-leaping method, the results which are summarized in the following theorem. Further discussion can be found in App. B.
Theorem 3.1 (Thm. 4.7 in [[8]], Under a certain discretization scheme and technical assumptions, and given an e-accurate score function, the following error bound holds:
DKL (r-6) Σ exp(-T) + KT,
(3.2)
where ā < is the early stopping time, K controls the step size, and T is the time horizon. The notation < indicates the inequality holds up to a constant factor as -+
The error bound (3.2) decouples three error sources of the r-leaping scheme: the truncation error O(e="), the score estimation error and the discretization error (kd"). Similar to the case for the Euler method for ODEs and the Euler-Maruyama scheme for SDEs, the r-leaping method is a first-order scheme in terms of the discretization error (KT)).
4 Algorithms: High-Order Inference Schemes
A natural improvement of r-leaping is develop high-order schemes for discrete diffusion mod- els. As a foundational example, consider the second-order Kunge-Kutta (RK-2) method with two stages [23] for solving the ODE dæ, fe(x))dt. This method represents one of the simplest high- order numerical schemes:
i4+0- 2, + fa, 3 ū. + [(1 (.)) + + fio (a(+a\)) A.
(4.1)
This scheme reduces to the exact midpoint method when 2 and Heun's method when 0 1. The underlying intuition stems from the observation that for f E (B)) [(1 3 f(0) + 1
28 f(A)) A
offers second-order approximation of f f(r)dx contrast to f(0) which is only first-order. This approach has been successfully adapted for SDE simulation [24]] and continuous diffusion model inference [48, 28, 29 gy, [Il]. Notably, these methods enhance sample quality and computa- tional efficiency without requiring additional model training. making the development of high-order schemes for discrete diffusion inference both theoretically appealing and practically viable.
In this section, we propose two high-order solvers for inference in the discrete diffusion model. We will primarily focus on two-stage algorithms aiming for second-order accuracy. Specifically, we will introduce the 9-RK-2 and T-Trapezodaal methods. Throughout this section, we assume time discretization scheme (8:))e(o:N) with 0 So stopping time and use the shorthand notations * max +). For any
< SN T 8, where 8 is the early
(Sn: 8n+1] and [0 N 1]. we define [s] = Sn, P. (1 1--0 0)s,, fsn+1. A., = 8n+1 - 8n, and f-section points as Pn. (1 0)sn + 0sa41. We choose (dv) to be the counting measure ID.
4.1 O-RK-2 Method
We first present the R-RK-2 method, which is simple in de- sign and serves as a natural analog of the second-order RK method for ODEs [I in terms of time and state-dependent Poisson random measures, as a warm-up for the O-Trapezoidal method. We note that simi- lar methods have been proposed for simulating SDEs driven by Brownian motions Poisson processes, such the stochastic [24]] and the Poisson ods. A summary of this method is given in Alg. MI.
Algorithm 0-RK-2 Method
Input: To qo. E (0, 1), (NN Output: A sample Tox T.y
for n 0 to N do
2
Tp + Da, + Even uP (Ds.. (v)0An);
3
Dents + Von +
EvER uP (1 20 [(1 29 Pe.. + 1 Pr (V)A.):
20
end
[54] RK meth-
Intuitively, the 0-RK-2 method is two-stage algorithm that:
5

---

(i) Firstly, it rums 7-leaping with step size A obtains an intermediate state Vy. at the f-section point Pa. and evaluates the intensity 2e there; (ii) Then another step of r-leaping for full step A, is run using weighted sum of the intensities the current time point Sn and the 8-section point Px-
We emphasize that our method differs from the midpoint method proposed in [52] for simulating chemical reactions, in which the Poisson random variable in the first step is replaced by its expected value. Such modification is in light of the lack of continuity and orderliness of the state space.
4.2 0- TrTapezziial Method
As to be shown theoretically and empirically, the conceptually simple 9-RK-2 method may have limi- tations in terms of both accuracy and efficiency. To this end, we propose the following 0- O-Trapezoidal method, which is developed based on existing methods proposed for simulating SDEs [25] and chem- ical reactions [26]. Below, we introduce two parameters that will be used extensively later:
(1-
a1 20(1 oj
and 02 20(1-0)
with 0j
1.
The e-Trapezoidal method is summa- rized in Alg. M Intuitively, it sepa- rates each interval (5n, 9n++1 into two sub-intervals (8p, Pa] and (po, Sn+1], on which simulations are detached with different intensities designed in a balanced way.
Algorithm 2: 8- e-Trapezoidal Method
Input: Do go. 0 E (0, 1]. (Sns Pn)nP(unnNNN--- 7²" Output: sample TaN ~ INN
for n 0 to N do
Tp + Don + EVED vP (v)OAN))
Vonti t Fp
Compared to the 0-RK-2 method, the o-Trapezoidal method is also two- stage with an identical first step. The second step, however, differs in two major aspects:
EvED vP (la17 agHra) (v) 0))A.::
end
(I) The second step starts from the intermediate state Vp. instead of De.. and only runs for a frac- tional step (1 - 0) rather than a full step An;
(2) The weighted sum is comprised of an altered pair coefficients (a), -a2). performing an extrapolation instead of interpolation with coefficients (1 20) as in the 9-RK-2 method with 0 E [}: 1]. This feature will be shown to render the algorithm unconditionally second-order.
Following the common practice in the literature [22], we reject updates with multiple jumps along one dimension in both algorithms,, ensuring their well-posedness A simple analysis shows that rejection only happens with probability (kx), and we refer to further details in Rmk. C.4. We refer to Props. Kw and e:3 for the stochastic integral formulations of these algorithms. We provide visual comparison between the 8-RK-2 and the o-Trapezoidal method in the right figure of Fig. MI
Theoretical Analysis
In this section, we provide the theoretical results of the 0-Trapezoidal and --K-2 methods. The goal of this section to show that under certain conditions, both methods are second-order accurate, improving from the first-order accuracy of the r-leaping method (cf. Thm. (.I. Our theoretical analysis also reveals that the 4-Trapezoidal method is more robust to the choice of f than 6-RK-2, to be confirmed by empirical results in Sec. R.
5.1 Assumptions
For simplicity, we impose a periodic boundary condition on the state space X [S] i.e. embed the state space in the d-dimensional torus T to streamline the proofs (cf. Rmk. mc!)
Assumption 5.1 (Convergence of Forward Process). The forward process converges to the station- ary distribution exponentially fast, i.e., DKr(Prlp>> S cxp(-T).
This assumption ensures rapid convergence of the forward process, controlling error when termi- nated sufficiently large time horizon T, and is automatically satisfied in the masked state case

---

and the uniform state case, given sufficient connectivity of the graph (ef. [38]. The exponential rate aligns with continuous diffusion models (cf. [b]l].
Assumption 5.2 (Regularity of Intensity). For the true intensity (v ys- and the estimated in- tensity (v. ye-), ), holds almost everywhere w.r.t Ha(v, Us- (y(dv)p,- (dy.-) that: (1) Both in-
tensities belong to (20 T 6]); (2) Both intensities are upper and lower bounded on [0, T 8].
This assumes two key requirements of the scores: (I) the forward process maintains sufficient smoothness, which is achievable through appropriate time reparametrization; and (2) if and only if state y E X is achievable by the forward process and v is permissible jump therefrom, then both its true and estimated intensity are bounded, corresponding to Assumps. 4.3(1, 4.4, and 4.5 [38].
Assumption 5.3 (Estimation Error). For all grid points and 0-section points, the estimation error of the neural network-based score is small, i.e., for any E UnE(o:N (8n, Pa), we have
J. (Ho(z) log Pa(v) -1) + Re(ws) ((dv)) S <;;
(1) E
[v(v)
J I -
(2) E
(a (u))))) P c11
This assumption quantifies the proximity of the estimated intensity j to the true intensity μe after sufficient training. Compared with [[8]. the additional L part in (2) is required for technical reasons, which is similar to [62, 5ll]. In practice. such additional assumptions may be realized by adding extra penalty terms to the objective function during training.
5.2 Convergence Guarantees
The following theorem summarizes our theoretical guarantees for the 9- T-aapezzodal method:
Theorem 5.4 (Second Order Convergence of f- T-Tapezzadaal Method). Suppose O E (0. 1] and a17%, 0z3 s for all [O.T 8). then the following error bound holds for Alg. D un-
der Assumps. to b
DKL (psll11197 S exp(-I)) + + + En)) R2T,
where 5 is the early stopping time, max An, i.e., the largest stepsize, and Tr-8 is the distribution obtained by Alg. W as defined in Prop.
The complete proof is presented in App. KL divergence between the corresponding path measures, as established in Thm. C.S. and then decompose the integral in the log-likelling ligliiillhood and bound respectively, where the primary technique used is Dynkin's 's formula (Thm. C.m. With term-by-term comparison with Thm. B. we observe a significant improvement in the discretization error term from (KTT)) to (x²T). This confirms that the T-Tippezoidal method achieves second-order accuracy given sufficient time horizon T and accurate score estimation, with empirical validation presented in Sec. G
The outline is to first bound DKL (p& (T-88 by the
Theorem 5.5 (Conditional Second-Oree Convergence of 9-RK-2 Method). Suppose A E (0. b and 0 for all s E [0, - 6). then the following error bound holds for Alg.
(1 - 20 1+ +
7 Je
under Assumps. D to X
DKL(p6 RRK S exp( + (a + E1)T + KTT,
cRK where 5 is the early stopping time, (o::N An. i.e., the largest stepsize, and 9T-8 5 is the distribution obtained by Alg. B as defined Prop.
The proof of the theorem above is provided in App. The restricted range of f is caused by one specific error term (III.4) (C.9) that permits bounding with Jensen's 's inequality only when 0 E (0, 11
similar to its counterpart (II.4) (CI in the 6- Trapezoidal method. The limitation arises partially
because the weighted sum with coefficients (1 28: W becomes an extrupolation only if 28 < 0,
feature that naturally holds for all 0 E (0, 1) in the 0-Trapezoidal method. These theoretical
findings are consistent with the empirical observations in Fig. of App. D.33. where the performance of 0-RK-2 method clearly peaks when 0 e (0, J
7

---

Remark 5.6 (Comparison between Trapezoidal and RK-2 Methods). Trapezoidal methods were originally proposed by [2S]] as a minimal second- order scheme in the weak sense for simulating SDEs. In simulating chemical reaction contexts, [20] claimed that trapezoidal methods also achieve second-order convergence for covariance error apart from the weak error, a property not shared by midpoint (RK-2) methods. Our empirical results partly reflect these findings, while we defer theoretical investigation of covariance error convergence in discrete diffusion models to future work. Remark 5.7 (Remark on the Positivity of Extrapolated Intensity). Due to the nature of extrapolation, both our theorems require an additional assumption on the positivity of the extrapolated intensity, which is classically assumed in [22, [B], and resolving this issue is long-stton stand-inn open problem. The best result so far Prop. 5 [26], claiming clamping the intensity above 0 only causes an error of order (K²)). for any large integer p. We empirically demonstrate the validity of this assump- tion Tab. in practice through the text generation task (Sec. (.2) and find that positivity occurs for both methods with high probability over 95% approaching 100% with increasing NFE. We refer to further discussion in Rmk. C.O.
Experiments
Based on the theoretical analysis, expect the 0 O-Trapezoidal method to outperform the r-leaping method and the 0-RK-2 method in terms of sample quality, given the same number function evalu- ations. This section empirically validates the anticipated effectiveness of our proposed 0 Trapezzoidal method (Alg. ( through comprehensive evaluations across text and image generation tasks. Our comparative analysis includes established discrete diffusion samplers basclines, e.g., the Euler method [33]] r-leaping [22], Tweedie r-leaping [G0), First-Hitting Sampler (FHS) [2I].. Parallel De- coding [63]] and Semi-Autoregressive (Semi-AR) sampler [64]]. We benchmark both uniform and masked discrete diffusion models, with experiment details provided in App. D
6.1 15-State Toy Model
We first evaluate the performance of the 0 0-Trapezoidal method using a 15-state toy model (d 1, S 15). The target distribution is uniformly generated from 1 with rate matrix Q 1 I, where E is the all-one and the iden- tity matrix. This setup provides analyti- cally available score functions, allowing isolation and quantification of numerical errors introduced by inference algorithms. We apply both the d-Trapezoidal and the 0-RK-2 method to generate 106 samples and estimate the KI. divergence between the true ground truth and the generated distribution QT.
10-1
Dulovila)
10
Rope a 1181p
10
slpp -1937
Trapezo dal
6-8622 Filld
10
6 Trapazo dal tad
RC2 Fitted list 4 poins))
16
32
64
128
256
512
1024
Number of Steps
Figure 2: Empirical KL divergence between the true and generated distribution of the toy model vs. number of steps. Data are fitted with linear regression with 95% confidence interval by bootstrapping.
For fair comparison, we choose 0 = 1/2 for both methods, and the results are presented in Fig. While both methods exhibit super-linear convergence the total number of steps grows, the 0- Trapezoidal method outperforms the 8-RK-2 method terms of both absolute value and conver- gence rate, while the 6-RK-2 method takes longer to enter the asymptotic regime. Moreover, the fitted line indicates that the 8- T-Tapeezodall method approximately converges quadratically with re- spect to the step count, confirming our theoretical results.
6.2 Text Generation
For the text generation task, we employ the pre-trained score function from RADD [B as our base model for benchmarking inference algorithms. RADD is masked discrete diffusion model with
GPT-2-level text generation capabilities [65]] and is trained on the OpenWeeText dataset [66] with d = 1024 and S = 50258. Our comparative analysis maintains consistent computational resources across methods, quantified through the number of score function evaluations (NFE), and evaluates the sample quality produced by FHS, the Euler method, r-leaping, Tweedie r-leaping, Semi-AR.

---

9.30
10.2
11.50
Zoomiih
Tuler
PMS
9.52
8.53
9-Trapezoidal INFE
@TTrapczoddal INFE 64)
22.44
Paralid secestig
7:37
Frechet Incoption Distanee
$7.54
7:00
0.50
32
0.30
0.50
0.50
0.80
0.90
1.00
13.42
70.00
#-rrppezoidal INFL 120:
55.00
@rTrapozoidal INFE
32.20
47.50
40.00
7.85
52.30
0.70
0.00
0.90
3.00
Number of Function Evalustion (NFE)
Figure 4: Sampling quality vs. 0 E (0, 1] in
Figure 3: FID of images generated by differ- ent sampling algorithms vs. number of function evaluations (NFE). Lower values are better.
0- Trapezood method. Upper: Image genera- tion (FID). Lower: Text generation (perplexity). Lower is better.
and our proposed 0- 0-Trapezoidal method. We generate text sequences of 1024 tokens and measure their generative perplexity (computed with GPT-2 Large [[5i] following the evaluation protocol established in 655
Tab. 01 presents the results for both low (128)
Table 1: Generative perplexity (on GPT-2 large) of texts generated by different sampling algorithms. Lower values better, with the best in bold.
and high (1024) NFEs, with comprehensive re-
sults across additional NFE values in Tab. B The empirical results demonstrate that the A- Trapezoidal method consistently produces bet- ter samples within a fixed computation bud- get than existing popular inference algorithms. Notably, it outperforms Fuler and Tweedie т- leaping, two of the best-performing samplers adopted by RADD, by a large margin. It also consistently prevails over FHS, which performs exact simulation at high NFE (1024), support- ing again our observations that being free of discretization error does not necessarily imply better sampling quality. These results validate the practical efficiency and accuracy of Alg. Additional evaluation results, including unigram en- tropy and generative perplexities evaluated under LLaMA 3 [677], are detailed in App. D.21
Method
NFE = 128
NFE- 1024
FHS
122.732
109.406
Buler
86.276
44.686
Tweedie r-leap.
85.738
< 44.257
r-leaping Semi-AR 0-RK-2
52.366
< 28. 797
S 360.793
< 147.406
64.317
< 36.330
T-Trapezoidal
< 49.051
< 27.553
6.3 Image Generation
Our experiments on image generation utilize the pre-trained score function from MaskGIT [63. 68] the base model, which can be converted into masked discrete diffusion model by introducing noise schedule (see App. D.3)). MaskGIT employs a masked image transformer architecture trained on ImageNet [X-]] of 256 x 256 resolution, where each image amounts to sequence of 256 dis- crete image tokens following VQ-GAN tokenization [00] (d = 256, S 1025). We evaluate the O-Trapezoidal method against FHS, the Buler method, r-leaping, and parallel decoding under equiv- alent NFE budgets ranging from 4 to 64. Following the setting in [633], we generate 5 x 104 images and compute their Fréchet Inception Distance (FID) against the ImageNet validation split.
Fig. B reveals that 0- Trapezoidal method (Alg. D consistently achieves lower (and thus better) FID values compared to both the Euler method and r-lapping across all NFE values. While FHS and parallel decoding show advantages at extremely low NFE (< 8), their performance saturates with increased computational resources, making them less favorable compared to our rapidly converging method. Additional results, including generated image samples (Fig. M. are detailed in App. DI
Remark 6.1 (Algorithm Hyperparameters), We evaluate the performance of the B- T-Trapezoidal method across various and NFE values for both text and image generation lasks. As illustrated in Fig. A. we observe that the O-Trapezoidal method demonstrates robustness to 0, with flat land- scupe near the optimal choice. Our empirical analysis suggests that O E [0.3, 0.5) consistently yields competitive performance across different tasks.

---

6.4 Diffusion Large Language Model and Math Reasoning
To verify the effectiveness of the proposed method on a scale, we additionally benchmark its perfor-
Table Response accuracy on GSM8K with different NFEs. The best results are in bold.
mance on L.LDDA-Instrct [64]], an 8B masked diffusion L.LM with one
Accuracy (%)
NFE = 64
NFE = 128
NFE 256
Semi-AR (Conf.)
33.6
32.0
39.1
of the best language modeling per- formances among discrete diffusion- based LLMs. We examine its perfor- mance on GSMSK [Z[]], a math math-reasoning dataset consisting grade-school-level problems. We compare 0- o-Trapezoidal to the semi-autoregressive (Semi-AR) sampler therein, with both confidence- based (Conf.) and purely random (Rand.) remasking strategies. For each method, we generate a response of 256 tokens in a zero-shot prompting manner, with NFE ranging from 64 to 256, and re- port the accuracy in Tab. D T-Trapezoidal outperforms the Semi-AR sampler in the low NFE regime, where NFE is strictly smaller than the sequence length. At high NFE regime, 0-Trapezoidal exhibits a similarly competitive performance as other solvers. This observation accords with our claim that high-order samplers perform better with lower NFE budgets, and that these advantages persist when the model size is large. Further implementation details are available at App. D44.
Semi-AR (Rand.)
33.8
31.3
40.3
0-Trapezoidal
35.1
38.4
39.7
7 Conclusion and Future Works
In this work, we introduce the 9-RK-2 and 6- o-Trapezoidal methods as pioneering high-order numer- ical schemes tailored for discrete diffusion model inference. Through rigorous analysis based on their stochastic integral formulations, we establish second-order convergence of the 0- o-Trapezoidal method and that of the 9-RK-2 method under specified conditions. Our analysis indicates that the 6- Trapezoidal method generally provides superior robustness and computational efficiency compared to the 0-RK-2 method. Our empirical evaluations, spanning both 15-dimensional model with pre- cise score functions and large-scale text and image generation tasks, validate our theoretical findings and demonstrate the superiority performance of our proposed 8- 8-Traezzoddal method over existing samplers in terms of sample quality under equivalent computational constraints. Additionally, we provide a comprehensive analysis of the method's robustness by examining the optimal choice of the parameter in our schemes.
Future research directions include comparative analysis of these schemes and development of more sophisticated numerical approaches for discrete diffusion model inference, potentially developing inference methods of higher order [566 or incorporating adaptive step sizes and parallel sampling methodologies. From the perspective of applications, these methods may also show promise for tasks in computational chemistry and biology, particularly in the design of molecules, proteins, and DNA sequences. Moreover, it would also be interesting to explore the usage of higher-order numer- ical solvers in other problem settings involving the inference of diffusion models, such sampling via discrete diffusion models [722-Z5], inference-ting scaling of diffusion models [[ZG--1]] improv- ing the reasoning ability of diffusion large language models 1
93 991, etc.
Acknowledgments and Disclosure of Funding
YC is supported by the National Science Foundation under Grants No. ECCS-1942523, DMS- 2206576, and CMMI-2450378. GMR is supported by Google Research Scholar Award. YZ and MT are grateful for partial supports by NSF Grants DMS-1841 847802, DMS-2513699, DOE Grants NA0004261, SC0026274, and Richard Duke Fellowship. YR. HC and LY acknowledge support of the National Science Foundation under Award No. DMS-220816
10

---

References
[1] Ting Chen, Ruixiang Zhang, and Geoffrey Hinton. Analog bits: Generating discrete data using diffusion models with self-conditioning. arXiv preprint arXiv:2208. S.Ii:2200.042222, 2022.
[2] Jacob Austin, Daniel D Johnson, Jonathan Ho, Daniel Tarlow, and Rianne Van Den Berg. Structured denoising diffusion models in discrete state-spaces. Advances in Neural Informa- tion Processing Systems, 19411781-179999 2021.
[3] Sander Dieleman, Laurent Sartran, Arman Roshannai, Nikolay Savinov, Yaroslav Ganin, Pierre H Richemond, Arnaud Doucet, Robin Strudel, Chris Dyer, Conor Durkan, et al. Con- tinuous diffusion for categorical data, arXiv preprint arXiv:221 1,5089, 2022.
[4] Griffin Floto, Thorsteinn Jonsson, Mihai Nica, Scott Sanner, and Eric Zhengyu Zhu. Diffusion the probability simplex. arXiv preprint wrXiv:2309002530 2023.
[5] Emiel Hoogeboom, Alexey Gritsenko, Jasmijn Bastings, Ben Poole, Rianne van den Berg, and Tim Salimans. Autoregressive diffussion models. In International Conference on Learning Representationss 2022. URL https://openreviiot wtt.://wwww..ooouuuoo/aaaaa11110
[6] Emiel Hoogeboom, Didrik Nielsen, Priyank Jaini, Patrick Forré, and Max Welling, Argmax flows and multinomial diffusion: Learning categorical distributions. Advances in Neural Information Processing Systems, 34:12454-12465, 2021.
[7] Chenlin Meng. Kristy Choi, Jiaming Song, and Stefano Ermon. Concrete score matching: Generalized score matching for discrete data. Advances in Neural Information Processing Systems, 35:34532-34545, 2022.
[8] Pierre H Richemond, Sander Dieleman, and Arnaud Doucet. Categorical sdes with simplex diffusion. arXiv preprint arXiv:22 14784, 2022.
[9] Haoran Sun, Lijun Yu, Bo Dai, Dale Schuurmans, and Hanjun Dai. Score-based continuous- time discrete diffusion models. In The Eleventh International Conference Learning Rep- resentations, 2023. URL https://openrer///
[10] Javier E Santos, Zachary R Fox, Nicholas Lubbers, and Yen Ting Lin. Blackout diffusion: generative diffusion models in discrete-state spaces. In International Conference on Machine Learning, pages 9034-9059. PMLR, 2023.
[11] Thomas J Kerby and Kevin R Moon. Training-fre free guidance for discrete diffusion models for molecular generation. arXiv preprint arXiv:2409.07559 2024.
[12] Nathan Frey, Daniel Berenberg, Karina Zadorozhny, Joseph Kleinhenz, Julien Lafrance- Vanasse, Isidro Hotzel, Yan Wu, Stephen Ra, Richard Bonneau, Kyunghyun Cho, et al. Pro- tein discovery with discrete walk-jump sampling. arXiv preprint arXiv:2306. 12366 2023.
[13] Pavel Avdeyev, Chenlai Shi, Yuhao Tan, Kseniia Dudnyk, and Jian Zhou. Dirichlet diffusion score model for biological sequence generation. In International Conference on Machine Learning, pages 1276-1301. PMLR, 2023.
[14] Wei Guo, Yuchen Zhu, Molei Tao, and Yongxin Chen. Plug-and-play controllable generation for discrete masked models. arXiv preprint arXiv:2410.02143 2024.
[15] Do Huu Dat, Do Duc Anh, Anh Tuan Luu, and Wray Buntine. Discrete diffusion language model for long text summarization. arXiv preprint arXiv:2407.1099) 2024.
[16] Minghui Hu, Yujie Wang, Tat-Jen Cham, Jianfei Yang, and Ponnuthurai N Suganthan. Global context with discrete diffusion in vector quantised modelling for image generation. In Pro- ceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11502-11511, 2022.
[17] Luke Causer, Grant M Rotskoff, and Juan P Garrahan. Discrete generative diffusion mod- els without stochastic differential equations: a tensor network approach. arXiv preprint arXiv:2407.11333 2024.

---

[18] Hunter Nisonoff, Junhao Xiong, Stephan Allenspach, and Jennifer Listgarten. Unlocking guidance for discrete state-space diffusion and flow models. arXiv preprint arXiv:2406.01572. 2024.
[19] Xiner Li, Yulai Zhao, Chenyu Wang. Gabriele Scalia, Gokcen Eraslan, Surag Nair, Tommaso Biancalani, Aviv Regev, Sergey Levine, and Masatoshi Uehara. Derivative-free guidance in continuous and discrete diffusion models with soft value-based decoding, arXiv preprint AXXx:22408.082522 2024.
[20] Yair Schiff, Subham Sekhar Sahoo, Hao Phung, Guanghan Wang, Sam Boshar, Hugo Dalla- torre, Bernardo P de Almeida, Alexander Rush, Thomas Pierrot, and Volodymyr Kuleshov. Simple guidance mechanisms for discrete diffusion models. arXiv preprini arXiv:2412. 10193, 2024.
[21] Kaiwen Zheng. Yongxin Chen, Hanzi Mao, Ming-Yu Yu Liu, Jun Zhu, and Qinsheng Zhang. Masked diffusion models are secretly time-agnostic masked models and exploit inaccurate categorical sampling. arXiv preprint arXiv:2409.0223 0220009,, 2024.
[22] Andrew Campbell. Joe Benton, Valentin De Bortoli, Thomas Rainforth, George Deligiannidis, and Arnaud Doucet. continuous time framework for discrete denoising models. Advances in Neural Information Processing Systems,
2022.
[23] John Charles Butcher. The numerical analysis of ordinary differential equations: Kutta and general linear methods. Wiley-Intersciens 1987.
[24] Kevin Burrage and Pamela Marion Burrage. High strong order explicit runge-kuttr methods for stochastic ordinary differential equations. Applied Numerical Mathematics, 22(1-3):88 101, 1996.
[25] David F Anderson and Jonathan C Mattingly. A weak trapezoidal method for a class of stochastic differential equations. Communications in Mathematical Sciences, 91):301-33 2011.
[26] Yucheng Hu, Tiejun Li, and Bin Min. weak second order tau-leaping method for chemical kinetic systems. The Journal of chemical physics, 135(2), 2011.
[27] Hideyuki Tachibana, Mocho Go, Muneyoshi Inahara, Yotaro Katayama, and Yotaro Watan- abe. Quasi-taylor samplers for diffusion generative models based on ideal derivatives. arXiv preprint arXiv:2112.13333 2021.
[28] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps. Advances in Neural Information Processing Systems, 35:5775-5788 2022.
[29] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu. Dpm- solver++: Fast solver for guided sampling of diffusion probabilistic models. arXiv preprint arXiv:2211.01095, 2022.
[30] Aaron Lou, Chenlin Meng, and Stefano Ermon. Discrete diffusiin modeling by estimating the ratios of the data distribution. In Forty-first International Conference on Machine Learning. 2024.
[31] Andrew Campbell, Jason Yim, Regina Barzilay, Tom Rainforth, and Tommi Jaakkola. Gener- ative flows on discrete state-spaces: Enabling multimodal flows with applications to protein co-design. arXiv preprint arXiv:2402.044 049977 2024.
[32] Itai Gat, Tal Remez, Neta Shaul, Felix Kreuk, Ricky T. Q. Chen, Gabriel Synnaeve. Yossi Adi, and Yaron Lipman. Discrete flow matching. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024. URL https://openreviee wtwpeeeerreieeeeeee/rurnnnt net/forum? 1d=CIDKO3SV9p.
[33] Jingyang Ou, Shen Nie, Kaiwen Xue, Fengqi Zhu, Jiacheng Sun, Zhenguo Li, and Chongxuan Li. Your absorbing discrete diffusion secretly models the conditional distributions of clean data. arXiv preprint arXiv:2406.03336, 2024.
12

---

[34] Jiaxin Shi, Kehang Han, Zhe Wang, Arnad Doucet, and Michalis K Titsias. Simplified and generalized masked diffusion for discrete data. arXiv preprint wrXiv:2406004529 2024.
[35] Subham Sekhar Sahoo, Marianne Arriola, Yair Schiff, Aaron Gokaslan, Edgar Marroquin, Justin T Chiu, Alexander Rush, and Volodymyr Kuleshov. Simple and effective masked diffusion language models. arXiv preprint arXiv:2406.07552 2024.
[36] Hongrui Chen and Lexing Ying. Convergence analysis discrete diffusion model: Exact implementation through uniformization. arXiv preprint arXiv:2402. 08005, 2024.
[37] Zikun Zhang, Zixiang Chen, and Quanquan Gu. Convergence of score-based discrete diffu- sion models: A discrete-time analysis. arXiv preprint arXiv:2410.02321, 2024.
[38] Yinuo Ren, Haoxuan Chen, Grant M Rotskoff, and Lexing Ying. How discrete and continuous diffusion meet: Comprehensive analysis of discrete diffusion models via stochastic integral framework. arXiv preprint arXiv:2410.0366 2024.
[39] Peter Eris Kloeden and Eckhard Platen. Numerical solution of stochastic differential equa- tions. Stochastic Modelling and Applied Probability, Applications of Mathematics, Springer, 1992.
[40] Peter Eris Kloeden, Eckhard Platen, and Henri Schurz. Numerical solution of SDE through computer experiments. Springer Science & Business Media, 2012.
[41] Qinsheng Zhang and Yongxin Chen. Fast sampling of diffusion models with exponential integrator. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreviin http:://pppppaaeeeeeeeeeeeoro/ee
[42] Qinsheng Zhang, Molei Tao, and Yongxin Chen. gddim: Generalized denoising diffusion implicit models. In The Eleventh International Conference on Learning Representations, 2023. URL https:://pperreviex. www.://www.eeeeeeeeeeeeeeeeeeeee www.eeee///toumieeeoeeeeoeoooeee
1nk 1VKE9q
[43] Martin Gonzalez, Nelson Fernandez Pinto, Thuy Tran, Hatem Hajri, Nader Masmoudi, et al. Seeds: Exponential sde solvers for fast high-quality sampling from diffusion models. Ad- vances in Neural Information Processing Systems, 36, 2024.
[44] Shuchen Xue, Mingyang Yi, Weijian Luo, Shifeng Zhang, Jiacheng Sun, Zhenguo Li, and Zhi-Ming Ma. Sa-solver: Stochastic adams solver for fast sampling of diffusion models. Advances in Neural Information Processing Systems, 36, 2024.
[45] Qinsheng Zhang, Jiaming Song, and Yongxin Chen. Improved order analysis and design of exponential integrator for diffusion models sampling. arXiv preprint arXiv:2308.02157. 2023.
[46] Tim Dockhorn, Arash Vahdat, and Karsten Kreis. Genie: Higher-order denoising diffusion solvers. Advances in Neural Information Processing Systems, 35:30150-30166, 2022.
[47] Luping Liu, Yi Ren, Zhijie Lin, and Zhou Zhao. Pseudo numerical methods for diffusion models on manifolds. arXiv preprint arXiv:2202.09777. 2022.
[48] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-basee generative models. Advances in Neural Information Processing Systems, 35: 26565-26577, 2022.
[49] Kaiwen Zheng, Cheng Lu, Jianfei Chen, and Jun Zhu. Dpm-solver-v3: Improved diffusion ode solver with empirical model statistics. Advances in Neural Information Processing Sys- tems, 36:55502-555 555422, 2023.
[50] Gen Li, Yu Huang, Timofey Efimov, Yuting Wei. Yuejie Chi, and Yuxin Chen. Accelerating convergence of score-based diffusion models, provably. arXiv preprint 2rXiv:2403.03852. 2024.
[51] Yuchen Wu, Yuxin Chen, and Yuting Wei. Stochastic Runge- Kutta methods: Provable accel- eration of diffusion models. arXiv preprint arXiv:2410.04766, 2024.
13

---

[52] Daniel T Gillespie. Approximate accelerated stochastic simulation of chemically reacting systems, The Journal of chemical physics, 115(4):171--1733, 2001.
[53] Yang Cao, Linda R Petzold, Muruhan Rathinam, and Daniel T Gillespie. The numerical stability of leaping methods for stochastic simulation of chemically reacting systems. The Journal of chemical physics, 121( 2004.
[54] K Burrage and T Tian. Poisson runge-kutta methods for chemical reaction systems ad- vances in scientific computing and applications, 2004.
[55] Yucheng Hu and Tiejun Li. Highly accurate tau-leaping methods with random corrections. The Journal of chemical physies, 130(12), 2009.
[56] Markus Arns, Peter Buchholz, and Andriy Panchenko. On the numerical analysis of inhomo- geneous continuous-time markov chains. INFORMS Journal Computing, 22(3):416-433 432, 2010.
[57] Desmond J Higham. Modeling and simulating chemical reactions. SIAM review, 50(2):347- 368, 2008.
[58] Weinan E, Tiejun Li, and Eric Vanden-Eijnden. Applied stochastic analysis, volume 199. American Mathematical Soc., 2021.
[59] Frank P Kelly. Reversibility and stochastic networks. Cambridge University Press, 2011.
[60] Joe Benton, Yuyang Shi, Valentin De Bortoli, George Deligiannidis, and Arnaud Doucet. From denoising diffusions denoising markov models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 86(2):286-301 2024.
[61] Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis. Nearly d-linear convergence bounds for diffusion models via stochastic localization. In The Twelfth Interna- tional Conference on Learning Representations, 2024. URL https://openreview .ntt//
forum jrrum//d=rinnV3BBuu
[62] Sitan Chen, Sinho Chewi, Holden Lee, Yuanzhi Li, Jianfeng Lu, and Adil Salim. The prob- ability flow ode is provably fast. Advances in Neural Information Processing Systems, 36, 2024.
[63] Huiwen Chang, Han Zhang, Lu Jiang, Ce Liu, and William T Freeman. Maskgit: Masked generative image transformer. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11315- 11115-11325, 11325, 2022.
[64] Shen Nie, Fengqi Zhu, Zebin You, Xiaolu Zhang, Jingyang Ou, Jun Hu, Jun Zhou, Yankai Lin, Ji-Rong Wen, and Chongxuan Li. Large language diffusion models. arXiv preprint arXiv:2502.0999) 2025.
[65] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, al. Language models are unsupervised multitask learners. OpenAl blog, 188)99,
[66] Aaron Gokaslan and Vanya Cohen. Openwebtext corpus. http://Skyliin0077.gtt
to/DpenneeTexxtCrppss 2019.
[67] Aaron Grattafiori, Abhimanyu Dubey, Abhinav Jauhri, Abhinav Pandey, Abhishek Kadian, Ahmad Al-Dahle, Aiesha Letman, Akhil Mathur, Alan Schelten, Alex Vaughan, et al. The Llama 3 herd of models. arXiv preprint arXiv2240..21883 2024.
[68] Victor Besnier and Mickael Chen. A pytorch reproduction of masked generative image trans- former. arXiv preprint arXiv:2310. /4400, 2023.
[69] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large- scale hierarchical image database. 2009 IEEE conference on computer vision and pattern recognition, pages 248-255 Icee, 2009.
14

---

[70] Patrick Esser, Robin Rombach, and Bjorn Ommer. Taming transformers for high-resolution image synthesis. Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 12873-12883, 2021.
[71] Karl Cobbe, Vineet Kosaraju, Mohammad Bavarian, Mark Chen, Heewoo Jun, Lukasz Kaiser, Matthias Plappert, Jerry Tworek, Jacob Hilton, Reiichiro Nakano, et al. Training verifiers to solve math word problems. arXiv preprini urXiv:21 14168, 2021.
[72] Peter Holderrieth, Michael S Albergo, and Tommi Jaakkola. Leaps: A discrete neural sampler via locally equivariant networks. arXiv preprint arXiv:2502. 2025.
[73] Yuchen Zhu, Wei Guo, Jaemoo Choi, Guan-Horng Liu, Yongxin Chen, and Molei Tao. Mdns: Masked diffusion neural sampler via stochastic optimal control. arXiv preprint arXiv:2508. 10684, 2025.
[74] Zijing Ou, Ruixiang Zhang, and Yingzhen Li. Discrete neural flow samplers with locally equivariant transformer. arXiv preprint arXiv:2505. 17741, 2025.
[75] Wei Guo, Jaemoo Choi, Yuchen Zhu, Molei Tao, and Yongxin Chen. Proximal diffusion neural sampler. arXiv preprint arXiv:2510.3382/ 2025.
[76] Masatoshi Uehara, Yulai Zhao, Chenyu Wang, Xiner Li, Aviv Regev, Sergey Levine, and Tommaso Biancalani. Inference-time alignment in diffusion models with reward-guided gen- eration:
Tutorial and review. arXiv preprint arXiv:2501.0966 2025.
[77] Nanye Ma, Shangyuan Tong, Haolin Jia, Hexiang Hu, Yu-Chuan Su, Mingda Zhang, Xuan Yang. Yandong Li, Tommi Jaakkola, Xuhui Jia, et al. Scaling inference time compute for diffusion models. In Proceedings of the Computer Vision and Pattern Recognition Conference, pages 2523-2534 2534, 2025.
[78] Raghav Singhal, Zachary Horvitz, Ryan Teehan, Mengye Ren, Zhou Yu, Kathleen McKe- own, and Rajesh Ranganath. A general framework for inference-time scaling and steering of diffusion models. arXiv preprint ar Xiv:2501100645. 2025.
[79] Marta Skreta, Tara Akhound-Sadege Viktor Ohanesian, Roberto Bondesan, Alán Aspuru- Guzik, Arnaud Doucet, Rob Brekelmans, Alexander Tong, and Kirill Neklyudov. Feynman- kac correctors in diffusion: Annealing, guidance, and product of experts. arXiv preprint urXiv:2503.02875 2025.
[80] Haoxuan Chen, Yinuo Ren, Martin Renqiang Min, Lexing Ying, and Zachary Izzo. Solv- ing inverse problems via diffusion-based priors: An approximation-ffee ensemble sampling approach. arXiv preprint wxXv:25506.3399, 2025.
[81] Yinuo Ren, Wenhao Gao, Lexing Ying, Grant Rotskoff, and Jiequn Han. Driftlite: Lightweight drift control for inference-time scaling of diffusion models. arXiv preprint arXiv:2509.21653 2025.
[82] Sophia Tang, Yuchen Zhu, Molei Tao, and Pranam Chatterjee. Tr2-d2: Tree search guided trajectory-aware fine-tuning for discrete diffusion. arXiv preprint arXi::2509.25171, 2025.
[83] Meihua Dang, Jiaqi Han, Minkai Xu, Kai Xu, Akash Srivastava, and Stefano Ermon. Inference-time scaling of diffusion language models with particle gibbs sampling. arXiv preprint arXiv:2507.08300 2025.
[84] Vignav Ramesh and Morteza Mardani. Test-time scaling of diffusion models via noise trajec- tory search. arXiv preprint arXiv:2506.0314 2025.
[85] Vincet Jain, Kusha Sareen, Mohammad Pedramfar, and Siamak Ravanbakhsh. Diffusion tree sampling: Scalable inference-time alignment of diffusion models. arXiv preprint urXiv:2506.207011 2025.
[86] Tianlang Chen, Minkai Xu, Jure Leskovec, and Stefano Ermon. Rfg: Test-time scaling for diffusion large language model reasoning with reward-free guidance. arXiv preprint Sttiv:2509.26000 2025.
15

---

[87] Xiangcheng Zhang. Haowei Lin, Haotian Ye, James Zou, Jianzhu Ma, Yitao Liang, and Yilun Du. Inference- time scaling of diffusion models through classical search. arXiv preprint WwXiv:2305.336/1 2025.
[88] Mohsin Hasan, Marta Skreta, Alan Aspuru-Guzik, Yoshua Bengio, and Kirill Neklyudov. Discrete feynman-kac correctors. In 2nd Al for Math Workshop@ ICML 2025, 2025.
[89] Guanghan Wang, Yair Schiff, Subham Sekhar Sahoo, and Volodymyr Kuleshov. Remask- ing discrete diffusion models with inference-time scaling. arXiv preprint 2xXiv:5533.00077 2025.
[90] Cheuk Kit Lee, Paul Jeha, Jes Frellsen, Pietro Lio, Michael Samuel Albergo, and Francisco Vargas. Debiasing guidance for discrete diffusion with sequential monte carlo. arXiv preprint arXiv:2502.06079, 2025.
[91] Zijing Ou, Chinmay Pani, and Yingzhen Li. Inference-time scaling of discrete diffusion models via importance weighting and optimal proposal design. arXiv e-prints, pages arXiv- 2505, 2025.
[92] Jiacheng Ye, Shansan Gong, Liheng Chen, Lin Zheng. Jiahui Gao, Han Shi, Chuan Wu, Xin Jiang, Zhenguo Li, Wei Bi, et al. Diffusion of thought: Chain-of-thought reasoning in diffusion language models. Advances in Neural Information Processing Systems, 37: 37:105345- 105374, 2024.
[93] Jiacheng Ye, Zhihui Xie, Lin Zheng, Jiahui Gao, Zirui Wu, Xin Jiang, Zhenguo Li, and Ling- peng Kong. Dream 7b: Diffusion large language models. arXiv preprint arXiv:2508. 15487, 2025.
[94] Fengqi Zhu, Rongzhen Wang, Shen Nie, Xiaolu Zhang, Chunwei Wu, Jun Hu, Jun Zhou, Jianfei Chen, Yankai Lin, Ji-Rong Wen, al. Llada 1.5: Variance-reduced preference opti- mization for large language diffusion models. arXiv preprint arXiv:2505. 19223, 2025.
[95] Siyan Zhao, Devaansh Gupta, Qinqing Zheng, and Aditya Grover. d1: Scaling reasoning in diffusion large language models via reinforcement learning. arXiv preprint arXiv:2504.12211 12216, 2025.
[96] Chengyu Wang, Paria Rashidinejad, DiJia Su, Song Jiang, Sid Wang, Siyan Zhao, Cai Zhou, Shannon Zejiang Shen, Feiyu Chen, Tommi Jaakkola, et al. Spg: Sandwiched policy gradient for masked diffusion language models. arXiv preprint arXiv:2510.0955
2025.
[97] Siyan Zhao, Mengehen Liu, Jing Huang, Miao Liu, Chenyu Wang. Bo Liu, Yuandong Tian, Guan Pang, Sean Bell, Aditya Grover, al. sraintiing-guiee policy optimization for diffu- sion large language models. arXiv preprint arXiv:2509. /0396, 2025.
[98] Cai Zhou, Chenxiao Yang, Yi Hu, Chenyu Wang, Chubin Zhang, Muhan Zhang, Lester Mackey, Tommi Jaakkola, Stephen Bates, and Dinghuai Zhang. Coevolutionary continu- ous discrete diffusion: Make your diffusion language model a latent reasoner. arXiv preprint arXiv:2510.03206, 2025.
[99] Yuchen Zhu, Wei Guo, Jaemoo Choi, Petr Molodyk, Bo Yuan, Molei Tao, and Yongxin Chen. Enhancing reasoning for diffusion lims via distribution matching policy optimization. arXiv preprint 10.iv:2510.08233. 2025.
[100] Ilia Igashov, Arne Schneuing, Marwin Segler, Michael Bronstein, and Bruno Correia. Retro- bridge: Modeling retrosynthesis with markov bridges. arXiv preprint arXiv:2308. 162/2, 2023.
[101] Yang Li, Jinpei Guo, Runzhong Wang, and Junchi Yan. From distribution learning in training to gradient search in testing for combinatorial optimization. Advances in Neural Information Processing Systems, 36, 2024.
[102] Zhiqing Sun and Yiming Yang. Difusco: Graph-based diffusion solvers for combinatorial optimization. Advances in Neural Information Processing Systems, 36:3706-3731, 2023.
16

---

[103] Naoki Murata, Chieh-Hsin Lai, Yuhta Takida, Toshimitsu Uesaka, Bac Nguyen, Stefano Er- mon, and Yuki Mitsufuji. G2d2: Gradient guided discrete diffusion for image inverse prob- lem solving. arXiv preprint arXiv:2410. 14710, 2024.
[104] Wenda Chu, Yang Song, and Yisong Yue. Split gibbs discrete diffusion posterior sampling. arXiv preprint arXiv:2503.01 161, 2025.
[105] Sarah Alamdari, Nitya Thakkar, Rianne van den Berg, Alex X Lu, Nicolo Fusi, Ava P Amini, and Kevin K Yang. Protein generation with evolutionary diffusion: sequence all you need. BioRxiv, pages 2023-09, 2023..
[106] Patrick Emami, Aidan Perreault, Jeffrey Law, David Biagioni, and Peter St John. Plug & play directed evolution of proteins with gradient-based discrete meme. Machine Learning: Science and Technology, 4(2):025014, 2023.
[107] Dmitry Penzar, Daria Nogina, Elizaveta Noskova, Arsenii Zinkevich, Georgy Meshcheryakov, Andrey Lando, Abdul Muntakim Rafi, Carl De Boer, and Ivan V Kulakovskiy. Legnet: best- in-class deep learning model for short dna regulatory regions. Bioinformatics, 39(8)bbaa4457, 2023.
[108] Joseph L Watson, David Juergens, Nathaniel Bennett, Brian L Trippe, Jason Yim, Helen E Eisenach, Woody Ahern, Andrew J Borst, Robert J Ragotte, Lukas F Milles, et al. De novo de- sign of protein structure and function with rfdiffusion. Nature, 620(7976):
2023.
[109] John J Yang, Jason Yim, Regina Barzilay, and Tommi Jaakkola. Fast non-autoregressive inverse folding with discrete diffusion. arXiv preprint
2023.
[110] Hannes Stark, Bowen Jing, Chenyu Wang, Gabriele Corso, Bonnie Berger, Regina Barzilay, and Tommi Jaakkola. Dirichlet flow matching with applications to dna sequence design. arXiv preprint arXiv: 2:xxv:2402.0884 2024.
[111] Kai Yi, Bingxin Zhou, Yiqing Shen, Pictro Liò, and Yuguang Wang. Graph denoising dif- fusion for inverse protein folding. Advances in Neural Information Processing Systems, 36, 2024.
[112] Yiheng Zhu, Jialu Wu, Qiuyi Li, Jiahuan Yan, Mingze Yin, Wei Wu, Mingyang Li, Jieping Ye, Zheng Wang, and Jian Wu. Bridge-if: Learning inverse protein folding with markov bridges. arXiv preprint arXiv:2411.02/20 2024.
[113] Patrick Esser, Robin Rombach, Andreas Blattmann, and Bjorn Ommer. Imagebart: Bidirec- tional context with multinomial diffuusinn for autoregressive image synthesis. Advances in neural information processing systems, 34:3518-200 05322, 2021.
[114] Jose Lezama, Tim Salimans, Lu Jiang. Huiwen Chang, Jonathan Ho, and Irfan Essa. Dis- crete predictor-correct diffusion models for image synthesis. In The Eleventh International Conference on Learning Representations, 2022.
[115] Shuyang Gu, Dong Chen, Jianmin Bao, Fang Wen, Bo Zhang, Dongdong Chen, Lu Yuan, and Baining Guo. Vector quantized diffusion model for text-to-image synthesis. In Proceedings of the IEEE/CVF conference on Computer Vision and Pattern Recognition, pages 10696-10706, 2022.
[116] Ari Seff, Wenda Zhou, Farhan Damani, Abigail Doyle, and Ryan Adams. Discrete object generation with reversible inductive construction. Advances in neural information processing systems, 32, 2019.
[117) Chenhao Niu, Yang Song, Jiaming Song, Shengjia Zhao, Aditya Grover, and Stefano Ermon. Permutation invariant graph generation via score-hased generative modeling. In International Conference on Artificial Intelligence and Statistics, pages 4474 4484. PMLR, 2020.
[118] Chence Shi, Minkai Xu, Zhaocheng Zhu, Weinan Zhang, Ming Zhang, and Jian Tang. Graphaf: flow-based autoregressive model for molecular graph generation. arXiv preprint arXiv:2001. 09388, 2020.
17

---

[119] Yiming Qin, Clement Vignac, and Pascal Frossard. Sparse training of discrete diffusion models for graph generation. arXiv preprint arXiv:231..0212,2250 2023.
[120] Clement Vignac, Igor Krawezuk, Antoine Siraudin, Bohan Wang, Volkan Cevher, and Pas- cal Frossard. Digress: Discrete denoising diffusion for graph generation. arXiv preprint
urXiv:2209. 14734, 2022.
[121] Kilian Konstantin Haefeli, Karolis Martinkus, Nathanaël Perraudin, and Roger Watten- hofer. Diffusion models for graphs benefit from discrete state spaces. arXiv preprint arXiv:2210.01549, 2022.
[122] Yiming Qin, Manuel Madeira, Dorina Thanou, and Pascal Frossard. Defog: Discrete flow matching for graph generation. arXiv preprint arXiv:2410.0422 2024.
[123] Jun Hyeong Kim, Seonghwan Kim, Seokhyun Moon, Hyeongwoo Kim, Jeheon Woo, and Woo Youn Kim. Discrete diffusion schrödinger bridge matching for graph transformation. arXiv preprint 240.02400001500 2024.
[124] Naoto Inoue. Kotaro Kikuchi, Edgar Simo-Serra, Mayu Otani, and Kota Yamaguchi. Lay- outdm: Discrete diffusion model for controllable layout generation. Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 10167-10176. 2023.
[125] Junyi Zhang. Jiaqi Guo, Shizhao Sun, Jian-Guang Lou, and Dongmei Zhang. Layoutdiffu- sion: Improving graphic layout generation by discrete diffusion probabilistic models. In Pro- ceedings of the IEEE/CVF International Conference on Computer Vision, pages 7226-7236, 2023.
[126] Seunggeun Chi, Hyung-gun Chi, Hengbo Ma, Nakul Agarwal, Faizan Siddiqui, Karthik Ra- mani, and Kwonjoon Lee. M2d2m: Multi-motion generation from text with discrete diffusion models. arXiv preprint arXiv 2407v:2407.14502, 2024.
[127] Yunhong Lou, Linchao Zhu, Yaxiong Wang, Xiaohan Wang, and Yi Yang. Diversemo- tion: Towards diverse human motion generation via discrete diffusion. arXiv preprint arXiv:2309.01372. 2023.
[128) Donghhao Yang, Jianwei Yu, Helin Wang, Wen Wang, Chao Weng, Yuexian Zou, and Dong Yu. Diffsound: Discrete diffusion model for text-to-sound generation. IEEE/ACM Transac- tions on Audio. Speech, and Language Processing, 31:1720-1733. 2023.
[129] Sam Bond-Taylor, Peter Hessey, Hiroshi Sasaki, Toby P Breckon, and Chris G Willcocks. Unleashing transformers: Parallel token prediction with discrete absorbing diffusion for fast high-resolution image generation from rector-quantized codes. In European Conference on Computer Vision, pages 170-188. Springer, 2022.
[130] Zhicong Tang, Shuyang Gu, Jianmin Bao, Dong Chen, and Fang Wen. Improved vector quantized diffusion models. arXiv preprint arXiv:2205. 16007, 2022.
[131] Ye Zhu, Yu Wu, Kyle Olszewski, Jian Ren, Sergey Tulyakov, and Yan Yan. Dis- crete contrastive diffusion for cross-modal music and image generation. arXiv preprint urXiv:220.007711. 2022.
[132] Zhichao Wu, Qiulin Li. Sixing Liu, and Qun Yang. Detts: Discrete diffusion model with con- trastive learning for text-to-speech generation. In ICASSP 0024-2024 IEEE International Con- ference on Acoustics, Speech and Signal Processing (ICASSP), pages 11336-11340 IEEE, 2024.
[133] Jun Han, Zixiang Chen, Yongqian Li, Yiwen Kou, Eran Halperin, Robert E Tillman, and Quanquan Gu. Guided discrete diffusion for electronic health record generation. arXiv preprint arXiv:2404. 12314, 2024.
[134] Juntong Shi, Minkai Xu, Harper Hua, Hengrui Zhang, Stefano Ermon, and Jure Leskovec. Tabdiff:: a multi-modal diffusion model for tabular data generation. arXiv preprint arXiv:2410.20626 2024.
18

---

[135] Zhengfu He, Tianxiang Sun, Kuanning Wang, Xuanjing Huang, and Xipeng Qiu. Diffusion- bert: Improving generative masked language models with diffusion models. arXiv preprint arXiv:2211.15029 2022.
[136] Nikolay Savinov, Junyoung Chung, Mikolaj Binkowski, Erich Elsen, and Aaron van den Oord. Step-unrolled denoising autoencoders for text generation. arXiv preprint arXiv:2 12.06749, 2021.
[137] Tong Wu, Zhihao Fan, Xiao Liu, Hai-Tao Zheng, Yeyun Gong, Jian Jiao, Juntao Li, Jian Guo, Nan Duan, Weizhu Chen, et al. Ar-diffusion: Auto-regressive diffusion model for text generation. Advances in Neural Information Processing Systems, 01.39957-399777 2023.
[138] Shansan Gong, Mukai Li, Jiangtao Feng. Zhiyong Wu, and Lingpeng Kong. Diffuseq-v2: Bridging discrete and continuous text spaces for accelerated seq2seq diffusion models. arXiv preprint arXiv:2310.05793, 2023.
[139] Lin Zheng, Jianbo Yuan, Lei Yu, and Lingpeng Kong. A reparameterized discrete diffusion model for text generation. arXiv preprint vrXiv:2302.05753 2023.
[140] Kun Zhou, Yifan Li, Wayne Xin Zhao, and Ji-Rong Wen. Diffusion-nat: Self-prompting discrete diffusion for non-autoregressive text generation. arXiv preprint arXiv:2305.0m .44404.00.0004. 2023.
[141] Minkai Xu, Tomas Geffner, Karsten Kreis, Weili Nie, Yilun Xu, Jure Leskovec, Stefano Ermon, and Arash Vahdat. Energy-based diffusion language models for text generation. arXiv preprint arXiv:2410.2133:
2024.
[142] Yinuo Ren, Grant M Rotskoff, and Lexing Ying. A unified approach to analysis and design of denoising markov models. arXiv preprint arXiv:2504.01938, 2025.
[143] Yong-Hyun Park, Chieh-Hsin Lai, Satoshi Hayakawa, Yuhta Takida, and Yuki Mitsufuji. Jump your steps: Optimizing sampling schedule of discrete diffusion models. arXiv preprint arXiv:2410.07761, 2024.
[144] Yixiu Zhao, Jiaxin Shi, Lester Mackey, and Scott Linderman. Informed correctors for discrete diffusion models. arXiv preprint arXiv:2407.21243 2024.
[145] Zixiang Chen, Huizhuo Yuan, Yongqian Li, Yiwen Kou, Junkai Zhang, and Quanquan Gu. Fast sampling via discrete mon-markov diffusion models with predetermined transition time. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
[146] Lingxiao Zhao, Xueying Ding, Lijun Yu, and Leman Akoglu. Improving and unifying discrete&ccontimuus-iiing discrete denoising diffusion. arXiv preprint 20X.v.240..03701, 2024.
[147] Machel Reid, Vincent Josua Hellendoorn, and Graham Neubig. Diffuser: Diffusion via edit- based reconstruction. In The Eleventh International Conference Learning Representations, 2023.
[148] Satoshi Hayakawa, Yuhta Takida, Masaaki Imaizumi, Hiromi Wakaki, and Yuki Mitsu- fuji. Distillation of discrete diffusion through dimensional correlations. arXiv preprint arXiv:2410.08709, 2024.
[149] Ludwig Winkler, Lorenz Richter, and Manfred Opper. Bridging discrete and continuous state spaces: Exploring the ehrenfest process in time-continuous diffusion models. In Forty-first International Conference on Machine Learning, 2024. URL https://openrevieet net/
forum CEFBYYGCCUEE
[150] Harshit Varma, Dheeraj Nagaraj, and Karthikeyan Shanmugam. Glauber generative model: Discrete diffusion models via binary classification. arXiv preprint arXiv:2405. /70.. 2024.
[151] Nate Gruver, Samuel Stanton, Nathan Frey, Tim GJ Rudner, Isidro Hotzel, Julien Lafrance- Vanasse, Arvind Rajpal, Kyunghyun Cho, and Andrew G Wilson. Protein design with guided discrete diffusion. Advances in neural information processing systems, 36, 2024.
19

---

[152] Severi Rissanen, Markus Heinonen, and Amo Solin. Improving discrete diffusion models via structured preferential generation. arXiv preprint urXiv:2405.7788%, 2024.
[153] Sulin Liu, Juno Nam, Andrew Campbell, Hannes Stärk, Yilun Xu, Tommi Jaakkola, and Rafael Gómez-Bombarelli. Think while you generate: Discrete diffusion with planned de- noising, arXiv preprint arXiv:2410.06264. 2024.
[154] Oscar Davis, Samuel Kessler, Mircea Petrache, Avishek Joey Bose, et al. Fisher flow match- ing for generative modeling over discrete data. arXiv preprint 1xXiv:2405.14664 2024.
[155] GN Mil' shtejn. Approximate integration of stochastic differential equations. Theory of Prob- ability & Its Applications, 19(3):557-562, 1975.
[156] Assyr Abdulle and Stephane Cirilli. S-rock: Chebyshev methods for stiff stochastic differen- tial equations. SIAM Journal on Scientific Computing, 30(2):997-1014, 2008
[157] Evelyn Buckwar and Renate Winkler. Multistep methods for sdes and their application to problems with small noise. SIAM journal on numerical analysis, 44(2):779 803, 2006.
[158] Denis Talay and Luciano Tubaro. Expansion of the global error for numerical schemes solving stochastic differential equations. Stochastic analysis and applications, 8(4):483-509 1990.
[159] Kevin Burrage and Pamela M Burrage. Order conditions of stochastic uune-kuuta methods by b-series. SIAM Journal on Numerical Analysis, 38(5):122611666, 2000.
[160] Kevin Burrage and Tianhai Tian. Predictor-corrector methods of runge-kutta type for stochas- tic differential equations. SIAM Journal on Numerical Analysis, 40(4):1516-1537, 2002.
[161] Andreas Rössler. Runge-kutte methods for the numerical solution of stochastic differential equations. Shaker- Aachen, 2003.
[162] Andreas RöBler. Runge-kutta methods for the strong approximation of solutions of stochastic differential equations. SIAM Journal on Numerical Analysis, 48(3):922-99) 952, 2010.
[163] James M Foster, Goncalo Dos Reis, and Calum Strange. High order splitting methods for sdes satisfying a commutativity condition. SIAM Journal on Numerical Analysis, 62(1):500-532 2024.
[164] Lei Li, Jianfeng Lu, Jonathan Mattingly, and Lihan Wang. Numerical methods for stochastic differential equations based on gaussian mixtures. Communications in Mathematical Sci- ences, 19(6):1549-133 1577, 2021.
[165] Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. Advances in Neural Information Processing Systems, 32, 2019.
[166] Nima Anari, Sinho Chewi, and Thuy-Duong Vuong. Fast parallel sampling under isoperime- try. arXiv preprint arXiv:2401.090) 2024.
[167] Lu Yu and Arnak Dalalyana. Parallelized midpoint randomization for langevin monte carlo. arXiv preprint arXiv:2402. 14434, 2024.
[168] Changyou Chen, Nan Ding, and Lawrence Carin. On the convergence of stochastic gradient meme algorithms with high-order integrators. Advances in neural information processing systems, 28, 2015.
[169] Alain Durmus, Umut Simsekli, Eric Moulines, Roland Badeau, and Gaël Richard. Stochastic gradient richardson-romber markov chain monte carlo. Advances in neural information processing systems, 29, 2016.
[170] Xuechen Li, Yi Wu, Lester Mackey, and Murat Erdogdu. Stochastic runge-kutta accelerates langevin monte carlo and beyond. Advances in neural information processing systems, 32, 2019.
[171] Sotirios Sabanis and Ying Zhang. Higher order langevin monte carlo algorithm. Electron. Statist, ((22:3805-3850 2019.
20

---

[172] Wenlong Mou, Yi-An Ma, Martin Wainwright, Peter L Bartlett, and Michael I Jordan. High- order langevin diffusion yields an accelerated meme algorithm. Journal of Machine Learning Research, 22(42): 2021.
[173] Pierre Monmarché. dighhdimmnsional meme with standard splitting scheme for the under- damped langevin diffusion. Electronic Journal of Statistics, (5(2):4117-416) 2021.
[174] James Foster, Terry Lyons, and Harald Oberhauser. The shifted ode method for underdamped langevin meme. arXiv preprint orXiv:2101.03446 2021.
[175] Kevin Burrage, PM Burrage, and Tianhai Tian. Numerical methods for strong solutions of stochastic differential equations: an overview. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460(2041):373-402. 2004.
[176] Grigori N Milstein and Michael V Tretyakov. Stochastic numerics for mathematical physics, volume 39. Springer, 2004.
[177] Alfred Bortz. Malvin Kalos, and Joel Lebowitz. new algorithm for monte carlo simulation of ising spin systems. Journal of Computational physics, 17(1):10-18, 1975.
[178] Daniel T Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of computational physics, 22(():403 434, 1976.
[179] Daniel T Gillespie. Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry, 81(25):2340-366) 1977.
[180] Yang Cao, Dan Gillespie, and Linda Petzold. Multiscale stochastic simulation algorithm with stochastic partial equilibrium assumption for chemically reacting systems. Journal of Computational Physicx, 206(2):395-411, 2005.
[181] Yang Cao, Daniel Gillespie, and Linda R Petzold, The slow-scale stochastic simulation algorithm. The Journal of chemical physics, 122(1), 2005.
[182] Weinan E, Di Liu, Eric Vanden- Eijnde, et al. Nested stochastic simulation algorithm for chemical kinetic systems with disparate rates. The Journal of chemical physics, 123(19), 2005.
[183] Weinan E. Di Liu, and Eric Vanden- Eijnden.. Nested stochastic simulation algorithms for chemical kinetic systems with multiple time scales. Journal of computational physics, 221 (1):158-180, 2007.
[184) Michael A Gibson and Jehoshua Bruck. Efficient exact stochastic simulation of chemical systems with many species and many channels. The journal of physical chemistry A, 104(9): 1876-1889, 2000.
[185] David F Anderson. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. The Journal of chemical physics, 127(21), 2007.
[186] Casper HL Beentjes and Ruth Baker. Uniformization techniques for stochastic simulation of chemical reaction networks. The Journal of Chemical Physics, 150(15). 2019.
[187] Muruhan Rathinam, Linda R Petzold, Yang Cao, and Daniel T Gillespie. Stiffness in stochas- tic chemically reacting systems: The implicit tau-leaping method. The Journal of Chemical Physics, 119(24):12784 12794, 2003.
[188] Daniel T Gillespie and Linda R Petzold. Improved leap-size selection for accelerated stochas- tic simulation. The journal of chemical physics, 119(1))8229 823344222299928,,,-3 2003.
[189] Kevin Burrage, Tianhai Tian, and Pamela Burrage. A multi-scaled approach for simulating chemical reaction systems. Progress in biophysics and molecular biology, 85(2-3):217-234, 2004.
[190] Yang Cao, Daniel T Gillespie, and Linda R Petzold, Avoiding negative populations in explicit poisson tau-leaping. The Journal of chemical physics, 123(5), 2005.
21

---

[191] Anne Auger, Philippe Chatelain, and Petros Koumoutsakos. R-leaping: Accelerating the stochastic simulation algorithm by reaction leaps. The Journal of chemical physics, 125(8). 2006.
[192] Yang Cao, Daniel T Gillespie, and Linda R Petzold. Adaptive explicit-implice tau-leaping method with automatic tau selection. The Journal of chemical physics, 126(22), 2007.
[193] Basil Bayati, Philippe Chatelain, and Petros Koumoutsakos. D-leaping: Accelerating stochas- tie simulation algorithms for reactions with delays. Journal of Computational Physics, 228 (16):5908-590 5916, 2009.
[194] Yang Cao and Linda Petzold. Slow-scale tau-leaping method. Computer methods in applied mechanics and engineering, 197(43-44):342234775 2008.
[195] Zhouyi Xu and Xiaodong Cai. Unbiased T-leap methods for stochastic simulation of chemi- cally reacting systems. The Journal of chemical physics, 128(15), 2008.
[196] Mary Schl, Alexander V Alekseyenko, and Kenneth Lange. Accurate stochastic simulation via the step anticipation r-leaping (sal) algorithm. Journal of Computational Biology, 16(9): 1195-1208, 2009.
[197] Krishma A Iyengar, Lconard A Harris, and Paulette Clancy. Accurate implementation of leaping in space: The spatial partitioned-leapii algorithm. The Journal of chemical physics, 132(9), 2010.
[198] David F Anderson and Desmond J Higham. Multilevel monte carlo for continuous time markov chains, with applications in biochemical kinetics. Multiscale Modeling & Simulation, 10(1):146-179) 2012.
[199] Alvaro Moraes, Raúl Tempone, and Pedro Vilanova. Hybrid chernoff tau-leap. Multiscale Modeling & Simulation, 12(2):581-615, 2014.
[200] Jill Padgett and Silvana llie. An adaptive tau-leaping method for stochastic simulations of reaction-diffusion systems. AIP Advances, 6(3). 2016
[201] Jana Lipková, Georgios Arampatzis, Philippe Chatelain, Bjoern Menze, and Petros Koumout- sakos. S-leaping: an adaptive, accelerated stochastic simulation algorithm, bridging r-leaping and r-leaping. Bulletin of mathematical biology. (188):3044-308 2019.
[202] Muruhan Rathinam, Linda R Petzold, Yang Cao, and Daniel Gillespie. Consistency and stability of tau-leaping schemes for chemical reaction systems. Multiscale Modeling & Simu- lation, ((3)) 2005.
[203] Tiejun Li. Analysis of explicit tau-leaping schemes for simulating chemically reacting sys- tems. Multiscale Modeling & Simulation, 6(2):417-436. 2007.
[204] Yucheng Hu, Tiejun Li, and Bin Min. The weak convergence analysis of tau-leaping methods: revisited. Communications in Mathematical Sciences, 9()):65-9996 2011.
[205] David F Anderson, Desmond Higham, and Yu Sun. Complexity of multilevel monte carlo tau-leaping. SIAM Journal on Numerical Analysis, 52(6):3106-3117, 2014.
[206] Chuchu Chen and Di Liu. Error analysis for d-leaping scheme of chemical reaction system with delay. Multiscale Modeling & Simulation, 15(4):1797-182) 2017.
[207] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep un- supervised learning using nonequilibrium thermodynamics. In International Conference Machine Learning, pages 2256-2265. PMLR, 2015.
[208] Linfeng Zhang, Weinan E, and Lei Wang. Monge-ampire flow for generative modeling. arXiv preprint 1809 10188, 2018
[209] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
22

---

[210] Jonathan Ho. Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Ad- vances in neural information processing systems, 33:6840-683 2020.
[211] Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13450 2020.
[212] Yang Song, Conor Durkan, lain Murray, and Stefano Ermon. Maximum likelihood training of score-based diffusion models. Advances neural information processing systems, 34: 1415-1428, 2021.
[213] Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXivv2210.02747, 2022.
[214] Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. arXiv preprint wxXiv:209.030033 2022.
[215] Michael S Albergo and Eric Vanden-Eijnden Building normalizing flows with stochastic interpolants. arXiv preprint arXiv:2209. 155711, 15571,
[216] Michael S Albergo, Nicholas M Boffi, and Eric Vanden-Eijnded Stochastic interpolants: A unifying framework for flows and diffusions. arXiv preprint arXiv:2303.08797, 2023.
[217] Haoxuan Chen, Yinuo Ren, Lexing Ying, and Grant M Rotskoff. Accelerating diffusion models with parallel sampling: Inference sub-linear time complexity. arXiv preprint arXiv:2405. 2024.
[218] Zhenyu Zhou, Defang Chen, Can Wang, and Chun Chen. Fast ode-based sampling for dif- fusion models in around 5 steps. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7777-7786, 2024.
[219] Yilun Xu, Mingyang Deng, Xiang Cheng, Yonglong Tian, Ziming Liu, and Tommi Jaakkola. Restart sampling for improving generative processes. Advances in Neural Information Pro- cessing Systems, 36.76806-76838, 2023.
[220] Hanzhong Guo, Cheng Lu, Fan Bao, Tianyu Pang, Shuicheng Yan, Chao Du, and Chongx- uan Li. Gaussian mixture solvers for diffusion models. Advances in Neural Information Processing Systems, 36, 2024.
[221] Jonathan Heck, Emiel Hoogeboom, and Tim Salimans. Multistep consistency models. arXiv preprint arXiv:2403.06807 2024.
[222] Yang Song. Prafulla Dhariwal, Mark Chen, and Ilya Sutskever. Consistency models. arXiv preprint arXiv:2303.0144 2023.
[223] Yang Song and Prafulla Dhariwal. Improved techniques for training consistency models. arXiv preprint arXiv:2310. 14189, 2023.
[224] Cheng Lu and Yang Song. Simplifying, stabilizing and scaling continuous time consistency models. arXiv preprint arXiv:2410.1 1081, 2024.
[225] Eric Luhman and Troy Luhman. Knowledge distillation in iterative generative models for improved sampling speed. arXiv preprint arXiv:2 023388, 2021.
[226] Chenlin Meng, Robin Rombach, Ruiqi Gao, Diederik Kingma, Stefano Ermon, Jonathan Ho, and Tim Salimans. On distillation of guided diffusion models. Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 14297-14306. 2023.
[227] Tim Salimans and Jonathan Ho. Progressive distillation for fast sampling of diffuusion models. arXiv preprint arXiv:2202.0055: 2022.
[228] Vinh Tong, Trung-Dung Hoang, Anji Liu, Guy Van den Broeck, and Mathias Niepert. Learn ing to discretize denoising diffusion odes. arXiv preprint 2005 2024.
23

---

[229] Eric Frankel, Sitan Chen, Jerry Li, Pang Wei Koh, Lillian Ratliff, and Sewoong Oh. S4s: Solving for a diffusion model solver. arXiv preprint arXiv:2502.17423, 2025.
[230] Tim Dockhorn, Arash Vahdat, and Karsten Kreis. Score-based generative modeling with critically-damped langevin diffusion. arXiv preprini urXiv:2112.07068. 2021.
[231] Hongkai Zheng, Weili Nie, Arash Vahdat, Kamyar Azizzadenesheli, and Anima Anandkumar. Fast sampling of diffusion models via operator learning. In International Conference on Machine Learning, pages 42390 42402. PMLR, 2023.
[232] Xinyin Ma, Gongfan Fang, and Xinchao Wang. Deepeache: Accelerating diffusion mod- els for free. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 15762-15772, 2024.
[233] Valentin De Bortoli, Alexandre Galashov, Arthur Gretton, and Arnaud Doucet. Accelerated diffusion models via speculative sampling. arXiv preprint arXiv:2501.05370 2025.
[234] Alexia Jolicoeur-Martinenu Ke Li, Rémi Mitliagkas. Gotta go fast when generating data with score-based models. arXiv preprint arXiv:2105.14080, 2021.
Tal Kachman, and loannis
[235] Wenliang Zhao, Lujia Bai, Yongming Rao, Jie Zhou, and Jiwen Lu. Unipe: A unified predictor-correctorr framework for fast sampling of diffusion models. Advances in Neural Information Processing Systems, 36, 2024.
[236] Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari. Parallel sam- pling of diffusion models. Advances in Neural Information Processing Systems, 36, 2024.
[237] Hyungjin Chung, Jeongsol Kim, Sehui Kim, and Jong Chul Ye. Parallel diffusion models of operator and image for blind inverse problems. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 6059-6069, 2023.
[238] Zhiwei Tang, Jiasheng Tang, Hao Luo, Fan Wang, and Tsung- Hui Chang. Accelerating paral- lel sampling of diffusion models. In Forty-first International Conference on Machine Learn- ing. 2024.
[239] Jiezhang Cao, Yue Shi, Kai Zhang. Yulun Zhang, Radu Timofte, and Luc Van Gool. Deep equilibrium diffusion restoration with parallel sampling. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2824-2834, 2024.
[240] Nikil Roashan Selvam, Amil Merchant, and Stefano Ermon. Self-refining diffusson samplers: Enabling parallelization via parareal iterations. arXiv preprint arXiv:2412.08229, 2024.
[241] Saravanan Kandasamy and Dheeraj Nagaraj. The poisson midpoint method for langevin dynamics: Provably efficient discretization for diffusion models. arXiv preprint arXiv:2405.17068. 2024.
[242] Shivam Gupta, Linda Cai, and Sitan Chen. Faster diffusion- sampling with randomized midpoints: Sequential and parallel. arXiv preprint arXiv:2406.0080 00924, 2024.
[243] Philip Protter. Point process differentials with evolving intensities. In Nonlinear stochastic problems, pages 467-472. Springer, 1983.
[244] Bernt Øksendal and Agnes Sulem. Applied Stochastic Control of Jump Diffusions. Springer, 2019.
24

---

A Further Discussion on Related Works
In this section, we provide a more detailed literature review of both continuous and discrete diffusion models, well as several studies on the numerical methods for SDEs chemical reaction systems, which are highly related to our work.
Discrete Diffusion Models: Methodology, Theory, and Applications. Discrete diffusion and flow-based models [M--), B-IO. 22] have recently been proposed generalizations of continuous diffusion models to model discrete distributions.
Such models have been widely used in various areas of science and engineering, including but not limited to modeling retrosynthesis [Лооо], combinatorial optimization [[[II], [02],, solving in-
verse problems [Лετ) [e] and sampling high-dimensional discrete distributions 190] [7 designing
molecules, proteins, and DNA sequences
[1005, [3, [O6. 12, 407-IOX BT. IT I.I.
image
synthesis [TE--TS) text summarization [TS]] well as the generation of graph [ [23], lay-
out [1124, 125), motion 2 I22ZI, sound [22 [128], image [TG 22 [E]]] speech [[332]. electronic
health record [|εلر], tabular data [[K3]]] and text [135)-440 54) 35, [411, [E]] Inspired by the huge
success achieved by discrete diffusion models in practice, researchers have also conducted some
studies on the theoretical properties of these models, such as [66-688 [42]
An extensive amount of work has also explored the possibility of making discrete diffusion models more effective from many aspects, such as optimizing the sampling schedule [143], adding correc- tors | developing fast samplers [[45]] designing correctors based on information learnt by the model | simplifying the loss function for training [146]] adding editing- based refinements [[47]]
synergizing these models with other techniques and methodologics like distillation []48], Ehren- fest processes [ Glauber dynamics [[500], tensor networks [[ enhanced guidance mecha- nisms [[51]L I8- 20J, structured preferential generation [[152], the plan-and-denoiss framework [[56]] and alternative metries, e.g., the Fisher information metric [[L54]] However, to the best of our knowl- edge, existing work on accelerating the inference of discrete diffusion models is relatively sparse compared to the ones listed above, which makes it a direction worth exploring and serves one of the main motivations behind this work.
Numerical Methods for SDEs and Chemical Reaction Systems. Below, we review advanced numerical methods proposed for simulating SDEs and chemical reaction systems, which the main techniques adopted in our work. For the simulation of SDEs driven by Brownian motions, many studies have been performed widely applied to tackle problems in computational physics, optimization, and Monte Carlo sam- pling. Examples of such work include the Milstein method [[SS]] explicit methods [[LS6]) multistep methods [[57]], extrapolation-type methods [[58, 22]], stochastic Runge Kutta methods [241, [3%- 1621, splitting methods [1663], methods based on gaussian mixtures [[64]] randomized midpoint method [1655], parallel sampling methods [1666] [6*l] as well as high-order methods for stochastic gradient Markov Chain Monte Carlo [168. 169], underdamped and overdamped Langevin Monte
design more accurate numerical schemes, which have been
Carlo [70-174) For a more comprehensive list of related numerical methods, one may refer
to [39 176, 40. 58).
Regarding the simulation of chemical reaction systems, numerical methods can be categorized into two classes. The first class consists of exact simulation methods, which are similar to the Kinetic Monte Carlo (KMC) method [[ZZ]] developed for simulating spin dynamics and crys- tal growth in condensed matter physics. Examples of such methods include the Gillespie algo- rithm (or the Stochastic Simulation Algorithm. a.k.a. SSA) [[78], 179] and its variants for mul- tiscale modeling [[80 I83], the next reaction method and its variants [|4 185], uniformization- based methods [[]6]] etc. The second class of methods are approximate simulation methods, including but not limited to the r-leaping method [52]] and its variants [[8Z, 188, 533, 54, 18g- 195, 555, 199, 56, 1197, 26, 198 201]] For a subset of the methods listed above, numerical analysis has also been performed in many works (06) to justify their validity.
Continuous Diffusion Models: Methodology, Theory, and Acceleration. Continuous diffusion and probability flow- models [207)-216] have also been the most popular methods in generative modeling, with a wide range of applications in science and engineering. For a list of related work on the thcoretical studies and applications of these models, one may refer to the literature review
25

---

conducted in [2T72 B8]. Here we will only review studies on accelerating the inference of continuous diffusion models, which motivates our work.
An incomplete list of accelerating methods includes approximate mean direction solver [218]. restart
sampling 121%% gaussian mixture solvers [22[]], self-consistency
knowledge distilla-
tion [225 229J. combination with underdamped Langevin dynamics [ operator learning [
and more recently ideas from accelerating large language models (LLMs) like caching [2 and
speculative decoding [|33]]_ Among all the proposed accelerating methods, one major class of meth-
ods are developed based on techniques from numerical analysis like adaptive step sizes [234]] expo-
nential integrators [9T 2.]. predictor- cerrrctttrrrrrttorr solver [235]] Adams-Bashforth Bashforth methods [29 44 35]
Taylor methods [27 - Picard iteration and parallel sampling [236 2401 2177], (stochastic) Runge-
Kutta methods [97] 28, 48 5 and randomized midpoint method [2411. 22 In contrast, there have been fewer studies on the acceleration of discrete diffusion models via techniques from numerical analysis, which inspires the study undertaken in this paper.
B Mathematical Background
In this section, provide the mathematical background for the stochastic integral formulation of discrete diffusion models, the error analysis of the r-leaping method, and useful lemmas for the theoretical analysis of high-order schemes for discrete diffusion models.
B.1 Stochastic Integral Formulation of Discrete Diffusion Models
Throughout this section, we will assume that (Ω, FPP) is a probability space, is finite-state space, and denote the pairwise difference set of the state space by ID := {x x o E X}. We also assume that the pairwise difference set X is equipped with metric || -|. a finite measure 7 and o-algebra B.
As warm-up, we introduce the definition of the Poisson random measure for time-homogeneous counting process.
Definition B.1 (Poisson Random Measure C& Definition A.I]). The random measure (dl, du) on R+ + x DD is called a Poisson random measure w.r.. measure if it a random counting measure satisfying the following properties:
(i) For any B E B and <t. N((s,t) x B) ~ (tt(B)tt - s));
(ii) For any t 2 0 and pairwise disjoint sets [Bi c B, (B.) = N(O0 x B.)
are independent stochastic processes.
Then we define the Poisson random measure with evolving intensities. The term "evolving" refers to that the intensity is both time and sate-dependent. Definition B.2 (Poisson Random Measure with Evolving Intensity [[8, Definition A.3]). Suppose (y) is non-negative predictable process on R+ D x satisfying that for any < T T. I
Jo (e())dt < oc, a.s.
The random measure [[A](dtt, dv) on R+ is called Poisson random measure with evolving intensity ><(v) w.nt. measure if is random counting measure satisfying the following proper- ties:
(i) For any B E B and < < l,
[[)]((s,t x B) ~ P (.. >((o)((a))a):
(ii) For any t 2 0 and pairwise disjoint sets (Bi)ie() CB
[N.(A](B.) == [[A] x B.))
are independent stochastic processes.
26

---

Remark B.3 (Construction of Poisson Random Measure with Evolving Intensity). As discussed in Thm. A.4 in [58] and originally proposed by [243), the Poisson random measure with evolving intensity can be constructed in the following way.
product space (TT B B(RR), m), where m is the Lebesgue measure on E, and B(RR) is the Borel o-algebra on RR The Poisson random measure with evolving intensity >e(1) can be defined in the augmented measure space as
One first augments the (X. v) measure space
N[A t) x B) = J. JaJ (uS N dv, de),
(B.1)
where Nidr, du, dg) is the Poisson random measure on RR+ D R w.E. measure v(dy)dç.
The following theorem provides the change of measure theorem for Poisson random measure with evolving intensity, which crucial for the theoretical analysis of numerical schemes for discrete diffusion models.
Theorem B.4 (Change of Measure for Poisson Random Measure with Evolving Density [38
Thm. 3.3]). Let (dt, dv) be Poisson random measure with evolving intensity >c(P), and (v(r) positive predictable process on IR x D x Q. Suppose the following exponential process is local Fe-martingal:
[e == exp (J. log hy bog x dv) (1((()) [(nle)-))((h(aa)))
(B.2)
and another probability measure on (St, F) such that Q < with Radon- Cadon-Nikodym derivative dQ/dPl/. = [b).
Then the Poisson random measure [()](dt, du) under the measure Q is Poisson random measure with evolving intensity heg().
B.2 Error Analysis of r-leaping
The r-leaping method originally proposed by [52] and adopted for the inference of discrete diffusion models by [22]. A summary of the algorithm is given Alg. B. In this subsection, we provide a sketch for the error analysis of the r-leaping method when applied discrete diffusion models, which will be compared with that of high-order schemes later on.
Algorithm 3: 7 r-Leaping Method for Discrete Diffusion Model Inference
Input: Do 9n. 0 [0, 1], time discretization (sn, Pn)ne P, î as defined in Prop. C.2.
Output: A sample Vay dx
for n to N do
I Jon+i < Vo. Σ vP (v)
3 end
i.e. the state space is d- dimensional grid with S states along each dimension, we have log | d log S. Then we consider a simple time-homogeneous transition matrix Q₂ = Q that allows jumps between neighhoring states with equal probability. Specifically, we have
Proof of Thm. 1 As we are considering the case where X [S]d d
1, -2d,
| vll 1,
Q(v.)) x)
T y.
which can be verified to satisfy Assumption 4.3(i) in [BS] with C' = 1 and D D = 2d. Assump- tion 4.3(i)) is also satisfied, shown in Example 1. of [[8]]
Then we may apply Thm. 4.7 in [X] by using the required time discretization scheme according to the properties of the target distribution and plugging in the corresponding values of C, D). The result follows by scaling the transition matrix Q by 3. equivalent to scaling the time by d. D
27

---

C Proofs
In this section, provide the missing proofs in the main text. We will first provide the proofs of the stochastic integral formulations of high-order schemes for discrete diffusion models in App. Cal. Then we will provide the proofs of the main results for the o-Trapezoidal method in App. C and the 8-RK-2 method in App. OK We remark that the proof for the 6- o-Trapezodaa method requires more techniques and is more involved, to which the proof for the 6-RK-2 method is analogous. In App. C.4. we provide the detailed lemmas and computations omitted in the proofs of Thms. 5.4 and 5.5
C.1 Stochastic Integral Formulations of High-Order Schemes
In order to rigorously analyze the 0-RK-2 method, we need the following definition: Definition C.1 (Intermediate Process). We define the intermediate process V piecewisely on (Sn, Sat1l follows:
To = Don. + [ vN [Pe..] (dss (ds,dv),
(C.1)
where the intensity P... is given by Ps.. (v, -=_)= 3. (To.: Non + v v)Q2, (Ben: Do.. + v), i.e, Th is the process obtained by performing T-leaping from time Sn with intensity p.
The following proposition provides the stochastic integral formulation of this method.
Proposition C.2 (Stochastic Integral Formulation 6-RK-2 Method). The 0-RK-2 method (Alg. m equivalent to solving the following stochastic integral:
RK = Us
Vo
RRK
uN UNNNNT [Z: (ds, du),
(C.2)
in which the intensity Tul is defined as weighted sum
RK Pe (BK(v) = [ (v, + 20 (v, (p.)))
RRK
(C.3)
and the intermediate intensity Ø³ is defined piecewisely as
(v. = (0, Vo + +. Th +v),
(C.4)
with the intermediate process V. defined in (C.l)) for the corresponding interval. We will call g the interpolating process of the 6-RK-2 method and denote the distribution of V. by q
The following proposition establishes the stochastic integral formulation of the 8-Trapezoidal method, whose proof can be found in App. Dill
Proposition C.3 (Stochastic Integral Formulation of 0-Trapezoidal Method), The O-Trapezoidal method (Alg. ( is equivalent solving the following stochastic integral:
+ L'' J. N[ptrap]]ds dv)
Trrap = Yo
strap
(C.5)
where the intensity rr is defined piecewisely as Tr. (v) = Lesp. ((o) + 1.20. (aiPe, (v. (v, (2, (top))))
(C.6)
Above, 1 denotes the indicator function and the intermediate process V is defined in (C.I)) for the corresponding interval. We will call the process garap the interpolating process of the f Trappezziaal method and denote the distribution of T by grrup.
Proof of Props. C and [K Without loss of generality, we give the proof on the interval (Sns Sa41 for n E [0
1]. and the generalization to the whole interval [o, T' straightforward.
Notice that once we condition on the filtration F., and construct the intermediate process D, as spec- ified in (C.I)) along the interval (8n: 8n+1], 1) the intermediate intensity 7 and the piecewise intensity MLa do not evolve with time the interpolating processes 2RK (or Dz
gtrap, respectively) since
28

---

only depends on the state, the intensity at the beginning of the interval Sn, and other randommess that independent of the interpolating process.
Therefore, the stochastic integral on this interval can be rewritten as for the 0-RK-2 scheme that
JRK
RRK Ven
[N de)
Ven-1
+
Vsn
RRK
+
VN UNTARK Sut1l, de)
= -RK
+
uP (1) 8a))
and for the 0- o-Trapezoidal scheme that
grap = Vsn
+
D [N[p dv)
= Ven
+
v. uN[atrup] 8n+1),
v/ (/P(P..P())
sn)))(d)),
and the statement follows by taking (dv) as the counting measure.
Remark C.4 (Remark on Rejection Sampling and Periodicity Assumption). The rejection sampling procedure in both algoritiums (Algs. W and ( guarantees well-posednest in the rare scenarios where a large drawn value of Poisson random variables or multiple simultanouss jumps in one coordinate would result in an update out of the state space X = [S]] To enforce this, we simply allow at most one jump per update across the summation, for example, the update
Ton f Usn. + >>'V vP
VED
as the standard practice in the literature [22, B8). The indicator function p 20 in Alg. M is also used to ensure that only valid jumps from the current state Us., are considered, while in Alg. D this implicitly guaranteed by taking the positive part of ai a2p which implies the positivity of on and thus the validity of the jumps Ton. We point out that the single-jump rule is only convenient sufficient condition, one should notice that this condition is not necessary for the well- posedness of our algorithums, since our setting of the state space carries both orderliness and algebraic structure, and thus one could in principle admit multiple simultaneous jumps without ambiguity.
Over the full inference process, the total probability of rejection is at most (R).. Below, we give brief justification and we refer to Proposition A.14 4 in [ES] for a complete proof of this claim. During the update aforementioned, the probability of at least two jumps occurring is bounded by
P
Σ (v))An)))) WED
=1
Σ WED
Poon (v)04n) $1
=1 exp
--E Pon (v)8An
1+ Σ
Ton (v)0An
wETO
DED
2 Σ Thn (u)Osn
24
$ A2
VED
Summing O ( over N steps gives EN-1 A2 S RT, and identical argument applies to the second update in each iteration. Hence, the overall rejection rate is at most Okk).
When we impose periodic boundary conditions, X = [S]d equipped with convenient algebraic structure: addition and scalar multiplication are globally well-defined. In that case, Algs. W and D match exactly the stochastic integral formulations Props. C2 and K- This alignment removes the need for per-slep rejection, streamlines the application of the chunge-of-measure argument, and greatly simplifies the convergence proofs of Thms. B and D.3. Even without periodicity, those theorems hold with probability at least 1 OKk), shown above.
29

---

C.2 Convergence Analysis of the 6- 6-Trapezoidal Method
Strap Theorem C.5. Let Po:I-6 and Jo:T- be the path measures of the backward process with the stochastic inlegral formulation (2.3) and the interpolating process (C.5) of O-Trapezoidal method (Alg. 2) then it holds that
strapy DRL. DRI.(Fr-oll 6ol
strap S DKL. (Po:T -olT
<DKL(Po L
T- / (Ha(v) log Pe (V( Ho(v) (.oop()) r(dr)ds
(C.7)
uis pr (v)
where the intensity Tr is defined in (C.S), and the expectation taken w.r.. both paths generated by the backward process (2.3) and the randomness of the Poisson random measure used in the first step of each iteration of the algorithm, i.e., the construction of the intermediate process (I..)) which is assumed to be independent of that of the backward process.
Proof. First, we will handle the randomness introduced by the Poisson random measure in the first step of each iteration of the 8-Trapezoidal method. For the ease of presentation, encode the afore- mentioned randomness as a random variable Ç and suppose
still supported on the probability space F,, while being independent of the backward process. Then for each realization of the intermediate process g is constructed as in (C.l) and the corresponding intensity Ø is defined in (C.4).
Given the stochastic integral formulation of the backward process (2.4) and the interpolating process of the T-Tappezoidal method (CCS), have by Thm. Bæl that this particular realization of the path strup
measure 88 can be obtained by changing the path measure Po-T-s with the Radon-Nikodym derivative
trap
exp (-f- log Malv) N N[ul(ds, dv) + S/" (HL(V) - (M)p(()) ((v))ds) strap p(rup(()) (v)
Ze
1
i.e.
2rrap
trap all90:T-os5) = E log -T- DKL (Po:T
11
I
/ (Ha(v) log
Me (v)
(v( + p p)tap(()) (v)) r(dr)ds
=E
rrra (v)
Us
Then it is easy to see by the data processing inequality and the chain rule KL divergence that
11 trap S DKL (PoT Araaa ,) I E DKL (Pr 107 [)]
T- / (ualv) log (a(v)
(P.llb +
ma(()) + () r(dv)ds)
and the proof is complete.
In the following, we will provide the outline of the proof of Thm. 5.4 where we leave the proof of several lemmas and detailed calculations to App. C44 for the clarity of presentation.
Proof of Thm. 4 Throughout this proof, including the subsequent lemmas and propositions that will be detailed in App. Me) we will assume that (ya)aglo,r1 is process generated by the path measure Po.r of the backward process with the stochastic integral formulation (2.z)) and set it as the underlying paths of the expectation in (CT) as required by Thm. w< Especially, y, ~ P,, holds for any & E [0. T]. For simplicity, we will assume that the process ys is left-continuous at each grid point s; for E j0 NN, which happens with probability one.
We first consider the interval (sn; 8x+11 for n E [0 Px = Pn. Within this interval, we will denote its intermediate process as appeared in (Cll v and the corresponding intermediate intensity as appeared in (C.)) as 7.5 the following discussion, will assume implicitly that the processes are conditioned on the filtration F...
1], and thus we have [8] 8,, and
30

---

By the definition of the intensity (v) as specified in (C.6)
Tr 1 Ml. 1.2p. (a177, (=)
we can rewrite the corresponding part of the integral in (CM2)) as
1n+1
J. [Mo(V)log log
(s(v) - 1s(1)) + ( y(dv)ds M8 ptrap (v)
Pn ponti f I ( log HalP)
ps(v) 23 mpt)p())) r(dv)ds
(v)
Mis
on
Pn
Pm
(s(()) log Mo(v)
ua(v) + Msn (1)) r(dv)ds
Man (v)
(1)
p
Ma(v) + ar (v) aziisn (v)) rdvv)ds, J (palw)los an7%, (v) a2Pon (v)
(a(v)
Pn
(IT)
where the assumption that a or2 [a] term (IT) above.
for all E [0, T 8] is applied here for the second
Decomposition of the Integral. Next, we decompose the integral (I) and (II) into several terms, the magnitudes of which or combinations of which are to be bounded.
(i) The first term decomposed as
(I) = (I.1) + (T.2) + (I.3) + (T.4),
where each term defined as
(1.1) = Pn
S. (Psn (v)
Hea (v) Men (v) Me, (v) + Tis (v) %ddrdss
on
(I.2)
(Ho(v) log (v) - pa(v) - Msa (v) log Mog (v) + Mon (v)) (dv)dss
(T.3)
(Ho(V) (e),(() (log (a1 (v) - a2 (v)) - log Pon, (v)) 7ddv)ds,
(1.4)
Men log log (a1 (v) azZ (v)) 7(dv)ds
8m
PPu
(s [s) (anP% a2ā (v)) 7dv)dds.
(ii The second term is decomposed
(II) (II.1) + (II.2) + (II.3 + (II.4) + (II.5) + (II.6).
where each term is defined as
[ (Ho.(v)log
Hon (v) (v)
Mpn (v) + tr (v) 7(dv)ds
(1I.1) = a1
pn 1 [ Mon (Mn(()) log
Mop (v)
Ho (v) + Men (v) %ddr)ds,
Non (v)
Pm
8nt1
(11.2) =
(HolP) (Holv)logus(u)) (((v))
Pm
sutt
[ (arlMon (v) log log (v) Plp, (v)) ar( (V) log Ms, (v) Me,, (v))) (du)ds,
on
31

---

(11.3)
a1 (Fip (v) up (v)) 7(d))ds,
Pn
(TI.4)
(01l4pp log P.p. (v) 02//s,. (V) logTs, (v)) (dv)ds
Pa
(01l4pa - a2pan (v)) log (a1 (v) anpis, (v) (dv)ds,
Pr
11
(11.5) =
(a1lpp (v) a2H4,, (2) log (a1 (v) a2 (v) (d
(sut1
(a1llp (v) a2paan (1)) log (1 (v) azfis (v) (dv)ds
Pn
8111
(11.6)
(a1kpp (v) a2He,, (v)) log (1 (v) azū (v)) (dv)ds
Pa
[ [u(()) log (ai (v) a2i (v) (dv)ds..
Pn
Bounding the Error Terms. Then we briefly summarize the intuitions and related techniques used the bounds of the terms above, and the detailed calculations and proofs of the lemmas and propositions are deferred to App. C4
(i) Error due estimation error associated with the intensity: The terms (I.1) and (II.1) bounded by the assumption on the estimation error of the intensity Tta (Assump. (.3)), as
E [1.11) + (II.1)] < 0 + an 0)A.er 0Amer + T8 Aner S AnEI,
for any 0 E (0, 1].
The term (11.4) is bounded by Prop.
as
E (1L.4) S An
where Jensen's inequality is applied here based the convexity of the loss.
(ii) Error related the smoothness of intensity: By Cor. CIK. the terms (I.2) and (II.2) bounded by
E [L.22 + (I1.2) < A By Cor. the terms (I.4) and (II.6) are bounded by
(I.4) + (II.6) s A
Intuitively, the bounds on these terms closely relate to the properties of the jump process and quantify the smoothness assumption on the intensity H4 (Assump. (5.2), especially when the intensity does not vary significantly within the interval (8n, 8n+1]. The main technique used for bounding these terms is Dynkin's 's Formula (Thm. CaI)). The third-order accuracy here directly follows from the intuition provided Sec. based on numerical quadrature.
(iii) Error involving the intermediate process: The terms (II.3) and (II.5) are bounded by Prop. C.IS and Cor. CT respectively follows
[(I.3) & A + A3 and [((1)) X A + A
The term (1.3) is bounded by Prop. 6021 as below
[[33) < A3
The three terms above involve the intermediate process y and the corresponding interme- diate density p
32

---

In conclusion, by summing up all these terms, we have
pn++ (HLle) log (=(v)
Mo(v) + Po (M) 7(dv)ds
[1 (
( + e11 A + en S (eler + +e) +
Therefore, the overall error is bounded by first applying Thm. os and then the upper bound derived above to each interval (Sn; Snt1l, which yields
patrap
T-5
/ (Hale))os
Ma(v) Mo(v) p. (1) rddv)ds
<DKL(Pol@O0) +
Psg (u)
N-1
(PKL(ool
+ Σ (e1 +en) +
Sexp(-T) T(er en) + TTT
as desired.
Remark C.6 (Discussion the Positivity Assumption). In the following. will take the positivity assumption in Thm. 15 as an example, and the case ohe the O-RK-2 method is similar. In the statement of Thm. 5.44 have assumed that
a' 01P a22
(v)
in (C..)) for all 8 E T-, which allows us to replace (a177). a by the difference itself. (25) showed that this approximation is most of ( within the corresponding interval,
and [26] further proved that for any order Z 1, there exists sufficiently small step size A such
that this approximation is least p-th order, i.e., of order O(D)) for that step.
We give a brief justification of this assumption here. We consider the expectation of the difference itself, which is given by
=E E [a (v) =B 21 a1 (161 =1 11- 0(k),
ai (Ap, (v) (o. + a1 (v) [1))
where we used E (v) [v)
Kеп as established in and I K. as shown (C..). Therefore, as long as the step sizes An
[v) Ph sufficiently small, the positivity assumption is valid in the sense that the expectation of the difference is at least 1 - (K)).
C.3
Convergence Analysis of the 6-RK-2 Method
Here may again apply the data processing inequality and the chain rule of KL divergence to upper bound the error associated with the 6-RK-2 method. A statement of the upper bound is provided in Thm. ea below, whose proof is omitted here since it is similar to that of Thm. C above.
Theorem C.7. Let Po:T- and <RK To be the path measures of the backward process with the
stochastic integral formulation (2.4) and the interpolating process ((..2) of the 0-RK-2 method ( Alg. D) then it holds that
DRL(FT-slH 1loRk
S DKL (Po.T--Slll RRK 81190-T-8)
{. [ (po(i) log (v(
T-
(C.8)
(Dol(foll@o) +
RRK (v) Ma(p) (v) r(dz)ds
where the intensity BRK defined (C..A), and the expectation is taken w.r.l. both paths generated by the backward process (2.4) and the randomness of the Poisson random measure used in the first step of each iteration of the algorithmm, i.e. the construction of the intermediate process (K.,)), which is assumed be independent of that of the backward process.
33

---

Following the same flow as in the proof of Thm. 54 we will first provide an outline of the proof of Thm. E.S, and defer the proof of several key lemmas and detailed calculations App. C for the clarity of presentation. We will also comment on the differences that may lead to the less desirable numerical properties of the R-RK-2 method.
Proof of Thm. 1.51 In the following proof sketch, we will be using the same notation as in the proof of Thm. 5.4 and we will assume that the process ys is left-continuous at each grid point 8; for i E 0: N]. We also start by taking closer look at the integral within each interval (#n, 8n+1 1] for n E 0 -11, and denote the intermediate process as appeared in (aal) as y and the corresponding
intermediate intensity appeared in (C:4 as jG-.
As defined in (C.3), the intensity f (v) is given by
p. (v) - (1 T [ Tp. (v),
which helps us rewrite the corresponding part of the integral in (C..) as
1
S (Hale)
(v(v) (v) (a(() + pTkk(v)) (v)) r(d))ds
log log
(v)
An 18n+1
S.
Ma (x(1) log
(1 28 Man (V) + zal1.. (v) - Ms(()) + (1 20) Ps. (v) + 20 (0)) %ddr)ds.
(a(U)
(II)
Above we again use the positivity assumption that (1 20 above, just as what we have done in the proof and discussion of Thm. 5 above.
1 | +
20 Po 0 for the term (III)
Decomposition of the Integral. Then we perform similar decomposition of the integral as in the proof of Thm. 5.4 as follows:
(II) = (II.I) + (111.) + (III.3) + (III.4) + (II1.5) + (II1.6,,
where each term defined as
S (aal(v) log
Toon (nil112 - Man (v) + Piso (1) r(dv)ds
(111.1)
(1
28
Jnn
1 +
8n+i
[o (Ppa (H log
Pop. (v) Men(v) + Pen (v) 7(dv)ds,
(Mp (v)
26
8nt1
(III.2) =
(Ha(V) log - (u((v))(ddv))
p8tt1
[ (11 20 (Man (v) log Men (V) - Moa) + 20 (Mon (v) Mo. - Mon (v))) (dv)dss
[
(111.3)
1 (Pp,, (v) - Pon (v) 7(dv)ds,
28
s8+1 (. ((1 28 Man (v) log Poo.. (v) + Mp (V) log Pip. (v)) (dv)dds
(III.4)
nentl
[ ((1 21) Men (v) + 20 Mp (v)) log ((1 20 Menn (v) +
20 Pm. (v)) (dv))ds,
1
Jen
pants [ ((1 20) Ms, (v) + 20 up. (v)) log ((1 21)) P... (v) 20 pp. (v)) ((v)dss
(III.5)
1hnt1 S ((1 (1-1)) -1)) psn, (v) + 1/2 My (v)) log ((1 1/2 Mo., (v) + T. (v)) (dv)dss
p1n { ((1 1 Msa (v) + 2/ Mpn (v) log ((1 1 Pi.. (v) 17 (v)) (dv)
(111.6)
34

---

Sut1 [ [us())log log ((1 20) usn. (v) 20 (v)) (dd)dds.
Bounding the Error Terms. Then we briefly summarize the intuitions and related techniques used in the bound of the terms above.. Detailed calculations and proofs of the lemmas and propositions used here deferred to App. C4
(i) Error due to the intensity estimation: The terms in (II1.1 are bounded by the assumption on the estimation error of the intensity Pa (Assump. 5.33 as follows
IE [(IIL.)) (1 1 Ane + An
for any E (0, 1].
(ii) Error related to the smoothness of intensity: By Cors. e and C.17. the terms (TII.2) and (TII.6) are bounded by
E (111.2)) S A. and E [1L1.6)]
A
respectively.
(iii) Eyror involving the intermediate process: The term (III.3 and (III.5) are bounded in almost
the same way that of Prop. R and Cor. C. By simply altering the integral upper limits, we obtain that
(IL.3) & A + AM
[(IT.)) & A + A
The only term that cannot directly bounded based on results App. C is (III.4), which given by
S [ ((1 2/0 Men (v) log (v)
E [IIL.4)]
E
(v) + 20 Mp (v) log Ppn. (v)) ((v)dss
/ ((1 201 Msn (v) 2 up (v)) log ((1 1 io.m (v) + Pp. (v)) (dv)dss
(C.9)
Recall that the proof of its counterpart (Prop. and the extrapolation nature of the second step in the T-Traeezoddal method (C.G) to bound the error term. However, the same technique cannot be directly applied to the 0-RK-2 method for any 0 [0, 1). as the intensity TaRK is interpolation of the intensity u. when ( 1). Therefore, below we will first focus on the case when (0, 21-
we utilized the convexity of the loss function
To be specific, by the assumption on the estimation error (Assump. 1, we can reduce (CC) to
[ ((1 1) Peo (v) log Pa.. (V) + 28 Mo, (v) log Pon, (p)
26
1n-1 / ((1 (11- 2 Man, (v) + 28 Pp.. (v)) log ((1 28 Ms,, (v) + Pr, (v)) (dv)ds
(C.10)
which can then be upper bounded based Jensen's inequality and the convexity of the loss function for 0 E (0, 4-
Summing up the bounds of the terms above, we have
1 [ log (v) 1 (v)
(b(1)) To (w) r(dv)ds
<An(et +e1) + A3
2 S + su)+ + A3
35

---

Consequently, the overall error of the 8-RK-2 method is bounded by
Jo T-5
/ (na(w) [a(v) (T)(())
(s(()) (11)( rddu)ds
(PRi(Fol(00) +
N-1
2 (Rol9ollo)
+ Σ Anle +
7=
+ T.
which suggests that the 0-RK-2 also of second order when 0 E (0, 1l- For the other case when O E ( 1]. we will provide brief discussion in the remark below.
Remark C.8 (Discussions on the case when 0 E and thus not necessarily bounded. One may wonder if. despite being positive, this term is still of least second order. However, the answer seems negative. By applying the Dynkin's 's formula (Thm. Call and Cor. (.ll) to μs log Pa in the term (III.A,, we have that the first integral in (C.9) can be expanded as follows
(5,1)) (5. 1].. For 0 E (5 11, the term (CET) is positive
Sm+1
/ ((1 26 Men (v) log Po., (v) + 20lpp (v) log Po. (v)) 7(dv)ds
E
An
pen+1
(Jon (v) log Pen (v) + AD, (Ho, (v) log Pus, (v)) (dv)dds
26
pan+1
[ Mo (v) Mo, (v)ddddd + 0(42)
28
sn
/ Man log (v)((v) + 2 / C (1r log Pso (v)) (dv)) + 0(43).
An
Similarly, applying Dynkin's 's formula to the the following function
G.(,,Vs-) = (V, (1 10) Ms, (o) log 1 Lo vo-) + 11) Ma. (v, (U,
with Go(v, Us..) = Ms, yo.) log Dia, (v, )s.)) allows us to expand the second integral (0.9) as below
8m--
/ (tt (v) (v)+ + +1 (1- 11) Men. (v)) log P.. + ((1 1)) Po.. (v)) 7ddv)ds
E
[ Gon ((ow) + + / (Vs, + 0
A.
36

---

where
LGs, (V.
20 Dspan (v, Uon ) log Msn. yom) + 1
20 Msn (v, yon 20
Pso (1, yon)
Mon Mo, + +1') log i
Ds (v, Man + v0) (1 20) Msn (1,yon Man v\) 7(de/)
S Mon (v, Von. log Poo (v, Von. V..Ver))(dd (r(dvv)
T, (V, Mon.)
+ (1 20 1) Man (v, Us..)
20 Ds. (V, Von)
+ (1 1 [( Msn (V, Vs.. + v') log 1 (s Usn + +1/ + (1 20) Dsn (U.Ya. + +z)) 7(de)
(1 # [ Men (V, Von.) log Per (v. Ven M (dd'')
Dutten (v, Mon. Men (V, Won.
20 0, (v, Mo..) log Pen ye.) )
20 Msn (v, Men
1
Msn (u,yo, ys.. v') log (ua Us.. + (((T/))))
20
(1 28 S Hon HonlV (v, Von + +2') log (u Mon + +zddddd))
S Msn (V, (s.) To, (v, Us., Us. (Tdd)) - -(1 -1- 210) I Psn (v. Us. log Msn.. (v, Us, (dde))
20
This further implies that
1 _ (Men (v) log Pa. (v))
OCG., (y+.)
+ 22 / (Mamle Usn + v') log fully, Ysn v') Msa. (V, Us..) log Ds.. (V,ys.) 3) i(de").
Comparing the first and second order terms in the two expansions of the two integrals in above then implies that the term (III.4) is of at most second order.
C.4 Lemmas and Propositions
In this section, we provide the detailed proofs of the lemmas and propositions omitted in the proof of Thms. 5.4 and 5.5.
Error due to the Intensity Estimation. Apart from the terms (I.1) and (II.1) in the proof of Thm. 5.4 and the term (III.1) in the proof of Thm. 5.3 we also need to bound the error terms (II.4) in terms of the intensity estimation error, which is given by the following proposition. Notably, the following bound also utilizes the convexity of the loss function and the extrapolation nature of the second step in the Orappezzodal method (CG))
Proposition C.9. For the interval (sn, 8n+1 1] for n E [0 N - 11. we have the following error bound:
[ pntt
[ (V) (v) 022Me, (v) log Pu.. (v)) 7(d))ds
E (II44) =B
Pm
1811
[ (aaHoo. (v) - azttan (v)) log (ai (v) azton (v)) (dv))ds (C.11)
<AnCII
37

---

Proof. We first define and bound three error terms (II.4.1,, (IL4.2), and (II.43) with score estima- tion error (Assump. b.33 as follows:
V/ 1at1
J on (P (v) log P (v) Pip.. log
E [(T...1)) E
log Po.. (v)) (dd))ds
Say L Jon sut2
(v)) / | (v) - (p)( |log Pi (v)) rdl)dds
SS (Mal(v) Pr ( x An
pn
Similarly, we also have
[/ Pm
(v) [ 02 (1),, log Iis., (v) Ts, (w) Jo, (v)) (v)) (d))dds s SAncn
E [ITT.442))]
E
and
E [II43)))) =E [| Sn+1
(artHen (v) azglinn (v)) (a1 (v) azp (v)) ((v))ds
On
/ (arTppi()) a22 (v)) log (a17pp (v) artian (v)) 7 (dv ds
SAnEI1
The remaining term (1I.4.4) = (11.4) (11.4.1) (11.4.2) - (II.4.3) is then given by
(11.4.4)
(a1 (v) log Apn. (v) azñis. (v) log a.... (v)) (dd))ds
Pn (aut1
[ (as (v) - azpue (v)) log (a1 (v) arpin (v)) 7(dv)ds s 0,
where the last inequality follows from Jensen's 's incquality, i.e.,
xr log x 022 logy S (arx a22y) log(ar - a2y);
for 01, 02 2 0 and 01 - 022 1. Therefore, by summing up the terms above, we have
E (11.4)) < E (II.411) + (11.4.2) + (11.4.3) + (114.4)) 2 AnEI,
and the proof is complete.
Error Related to the Smoothness of Intensity. Below we first present the Dynkin's formula, which is the most essential tool for the proof of the error related to the smoothness of the intensity. Theorem C.10 (Dynkin's 's Formula). Let (Ue)ce(o.r) be the following process:
90 + S. WNTTu(dss, dv),
1/1 90
where (ds, dv) a Poisson random measure with intensity a of the form ys-). ). For any f C¹ ( 7] x X), we define the generator of the process (ye)re(o.r) as below
[e+r()e+)) Se(U)) yt
= =d.f+(e)+ / (fil+ - fell))
Cf(l))
lim 7-0+
T
(C.12)
Then we have that
E [f(()))] = fo())) + E [ Cfa(ya)ds)
38

---

Proof. The definition and the form of the generator L, as well as the Dynkin's 's formula are all well- known in the literature of jump processes. We refer readers to detailed discussions on these topics in [244]]
Here take (t) (t, va). ('osesm((.N)-) in the statement of Thm. 1.19 in [EE and replace the compensated Poisson random
(v, att, (())
(tt, 0, 7t,, (tt)),
measure [dt, dz) with the Poisson random measure N (dss, dv, dg) defined as Rmk. B Then we allowed to use the ordinary Poisson random measure instead of the compensated one since are working with a finite measure r(dv).
From Thm. 1.22 in [244], we have that
(f(()) = =.ff(()) + J. (fily (oosemtos))))) fe(w)) r(dv)dg
=.f((y)+ / (f(v+ -l.))
and the proof complete.
In many cases below, we will need the following first-order expansion of the expectation of the function fele)) by assuming the second-order smoothness of the function f.
Corollary C.I1. Suppose that the process (Wt)ee(o.r] and the generator L are defined as in Thm. C.I. If we further assume that f E (10, x X). then it holds that
IR. [f((x:)] = fo()n) + (V)) + 0(17)
Proof. We expand the function fs())) from / = as follows
E (f+()) =fo())) + TE J. c.s.cu)dss)
=so(yo) + IE L' (1ol60)) + L' CLo CIolua)do do
=oo(w)) + Cfo(vo)t + E [T (-" f2
where L the second-order generator of the process (We)etoor) defined as follows
( ( (v)
So(y)) / - (o = -sS(v) +2 (Dololy + (
(3 =
+ / (fo + v) fo())) (opo(v)
S (fo (foly -6(- fo (0
fo
which is bounded uniformly by constant based the assumption the smoothness of the function up to the second order and the boundedness of the measure 7(v)). Therefore, the second-order term above of magnitude Ot²2), and the proof is complete.
The following lemma provides general recipe for bounding combination of errors, which resem- bles standard analysis performed for numerical quadratures. In fact, the following lemma can be easily proved by Taylor expansion when the process (J0) is constant, i.e., Ye y. Cor. Call offers analogous approach to perform the expansion when the process (yt)et0or) is not constant.
39

---

Lemma C.12. For any function f E (20 r] x X) and the true backward process (Vi)eefo.r
defined in (Z.A). holds that
folvo)ds + [(a (a1 a22o((0)) ds L Is(ua)ds) S
E
Proof. Let L be the generator defined in Thm. cam By applying the Dynkin's 's formula (Thm. on and Cor. C.Il) to the function fe(V)) and plugging in the expression of the generator L, we have that
0 folyo)ds 02 S fo - a1 [[ for (yor) (oserrds_/A) - -f- fa(v)dds
dr
E8
0r foloo) a2(1 - -0) + a1l1 1111-1 0)T (fu(ye) + OrLfu(uo))
/-" (Jo(vo) + +LSolvo)) ds + 0
72 (0 - a1 0) + a ( 0) fo + a1 Lfo(wo) - (fo(1) + O(+9),
which is the order ((+²) by noticing that
0 - a2 0) + 0) 1- ( 02+(1-0)))
20(1-0) (1- - 0)=
0)0 2е(1 (1 0)0 1 0,
and the proof is complete.
We remark that in Thm. eam Cor. CAll, and Lem. Mi the smoothness of the function f implies
that its derivatives up to the relevant order are bounded by constants independent of the time step This condition is verified the subsequent proofs.
Then we are ready to bound some of the error terms in the proof of Thm. 5.4 with Lem. C.I2
Corollary C.13. For the interval (8n, (8n, 8011 E [0: N - 1], we have the following error bound:
IE ((1.2) (II.2)]
S( (ua(v) (og) ps (b)) r(d)dds
E
(Man (V) Usn Usm (v) Usa (v)) 7dvvdds
on
1844
f (ailHoa (v) log Me (v) up (v)) - az(Men (v) log Mea (v) Msn (v)))
SA3 3
Proof. The bound obtained by applying Lem. C.2 with being the function
(u,) = J (v( log (v))((d)))
Strictly speaking, (ys) is actually in the form of (V,-)), but the argument can be easily extended this case by assuming time continuity of the function f.
Corollary C.14. For the interval (8ns Sor) E [0 N - 1}, we have the following error bound:
[L [(1.4 (11.6)]
[ Msn. (v) (a10% (V) azū (V)) 7dd))ds
[( (a1/4t. - azpus (v)) log (11 (v) aziiann (v)) r(dv)ds
Pa
8
/ (o(() log (ai (V) azm (v)) r(de)ds SA.
40

---

Proof. Note that the intermediate process y defined in (C.II) is driven by a Poisson random mea- sure that independent of the Poisson random measure driving the process ys within the interval (8n, Sn+1]. Therefore, the error bound is obtained by
(1) Taking the expectation w.r.t. the intermediate process y and thus the intermediate intensity Fg and
(2) Then applying Lem. with f being the following function
fe(ya) (o( [log (a1 (v) azãan (v)) ((v)).
The result follows directly.
Now we turn to the emror term (II1.6) Thm. 5. for which we need the following variant of Lem. (
Lemma C.15. For any function E (² r] x X) and the true backward process (Ve)ke(o.r)
defined in (2.3), it holds that
((1 20) So(vo) + 20 sor (Jo+)) ds - L (.lu.Jds) z Z
E
Proof. The proof is similar that of Lem. Specifically, we let be the generator defined in Thm. will apply the Dynkin's formula (Thm. call and Cor. C.TT) to the function fe()) and plug in the expression of the generator C, which yields
((1 20) fo(vo) + (yo=)) ds - L fs(va)ds)
E
(1 20 (5o())) + $ / (Fo(so) + +T do f- (Folvo) + ds 0(73)
=0(73),
as desired.
Corollary C.16. For the interval (8ns 8141 [0: N 1]. we have the following error bound:
IE
[11I.))]]
[ (Halv)log ps(v)) ((do)ds
E
Su
[ ((1 (1- 10) (p),(v)) (v) Ms,, (v) (ss)) ) + 20 (Mp, (v) log pp (v) Mp, (v))) (dv)ds
8n
S 3 n
Proof. By applying Lem. C.IS with being the function
fa(v)) = [ (e(v) log uo (2)(dd))
we have that the result follows directly.
Corollary C.17. For any following error bound:
[0 N 1] and the corresponding interval (Sns Sn41), we have the
E [11L,))]]
[((1- --) Ms, (v) + 2Φ Pp.. (v)) log ((1 1/1) Na. (V) + to (V)) (ddv)ds
( (o( log ((1 - Pon (v) + 20 fy (v)) >ddz)ss
$ A4.
41

---

Proof. Following the arguments in the proof of Cor. C.4 the error bound is obtained by first taking the expectation w.r.t the intermediate process y and thus the intermediate intensity Pig and then applying Lem. I.I3 with being the function
fa(va)= / [s(U)EE [log ((1 20) Ps,, (V) + 29 pr (v))] 7(de),
as desired.
Error involving the Intermediate Process. Proposition C.18. For the interval (8ns Sn41] with n [0 N 11, we have the following error bound:
E (113)) = IB [
S (ipa (v) - (pp( (v)) 7dd))ds
SA A
Pm.
Proof. First, we rewrite the error term (11.3) as
18441
/
E [[I.3)]
=E
7(dr)ds
(C.13)
E
(v)]
E pm (1))) (dv))ds.
Pn
Then expand the integrand by applying the Dynkin's formula (Thm. eaul and Cor. M the function (v) w.rt. the intermediate process (y") and the process (ye (U))xc(snnen respec- tively follows
E [P7 (v) - -E [7p. (v)]
=E Pe.. (v) + Ae (v) + ((A)) (v) + Cf.. (v)An + O(A
=E C),,, (n) (AA)),
where the generators L and L are defined in w.r.t the processes (V (V")eE(sniPn)) and (Ne)se(sn.pa)) respectively, i.e., for any function f E (1 x X), we have
('f((y)=0,f5 = =.fx(v) [ (fa(v-v)- P.,, (v)
(C.14)
Cfs(y) = 0.Sx(v) (5.(u ( - fx()) n. (d)) Therefore, for the term E E[[(2- - L) (v)l] evaluated at Sn, we have
E [](C- L) (v) = E [1./ (Msn v) Poa (v)) (Pen (v) Msn (v)) (d))
(C.15)
$ E J/ 10=,, (1) Mo.. (v) (dv) S e11.
where used the assumption on the estimation error (Assump. 5.3 in the last inequality. Then can further reduce [<])) to
Sm++
S ( (v)
E Ppn (V))) Adv)ds $ Lrae (cnAn +0(A2) ds I cn + +
and the proof complete.
Corollary C.19. For the interval (Sn, snr [0: 1]. we have the following error bound:
[11.5) =B
(ailton (V) aglle, (v)) log (a Ap. (v) azzus,, (v)) 7ddv)ds
1a111
/ (aiHpp (v) a2 log (a1 (v) - axp (v)) 7(dv)dd
<A + A3,11
42

---

Proof. Since the two integrands in (II.5) only differ by replacing 8. (v) with Po.. (v), we have the following upper bound by using the assumption on the boundedness of the intensities (Assump. 12.2
(I))
[II.5)]
J (a1 (v) a2H4sn(1)
SE
an (v) azH (v) az 10 (v) - Pp (v) r(dr)ds
Jran I |monl()) My.. (v) r(de)do
SE
Pm
-AnLE 1E J. [Pn(())- -pp will r(dr)
(C.16)
Applying the same arguments as Prop. C. which uses the generators C and L defined in ( we can bound the RHS above as follows
E [1pp (v) - Pp.. (v)]l
-E (P.. (v) L Ps. (v) 0(43)) (v) (V),, + +04
(C.17)
SA,E L)),, (v)]] + O x F1 + O
where the last inequality follows from Substituting into (C.G then yields the desired upper bound.
Proposition C.20. For the interval 8n+1) with E (0: N 1}, we have the following error bound:
[1,mm
((HLlV) Pon (v)) (log (a1 (v) - a2ztt, (v)) - log Na., (v)) r(dv)ds
E [1.3)] =E
SA3 en
Proof. First, observe by Dynkin's formula (Thm. that
| Cua, ds + +0(4))
E [u)(V)) Ms, (v)]] E
x <
and also
(1)) P... (v)l) | C' ds 0
- An.
(C.18)
Secondly, applying the given assumption (Assump. yields
2 (II) on the boundedness of the intensities
E. [log (017%, (v) ar2 (v) log 7r. (v)
IR. [ (v) a2f (v) Pi (v)
nen
Jaipp (v) - 022 (v) -
(v)
(C.
aai pppp (v) - - [v)]
[1 (v) - pp. [v)l] + [1 (v) - Dia, (v)l SA, + A, + ( S An
where the last inequality follows from (C.17) proved above. Therefore, may further deduce that
E
[I.3)
Pm
<
E [1))(v) Mo (v)l
E log (a1Ai, (v) azRon (v)) log (a1 (v) - auzee, (v)) (dv)ds
<A3 where the first inequality due to the independency ys and y for Pn), and the proof complete.
43

---

D Details of Numerical Experiments
In Apps. DJ D we present additional numerical results for the 15-dimensional toy model, text generation, image generation, and diffusion LLM, respectively.
D.1 15- 15-Dimensional Toy Model
We first derive the closed-form formula of the marginal distributions p. in this model. Recall that the state space X = {1,2,-. 1.22....4 d} with d = 15, and the initial distribution is po E A". The rate matrix
any time Q E By solving () we see that
p: = e = ( po:
and therefore Pt converges to the uniform distribution Poo = 1 as 1 co. The formula of directly yields the scores (x(x) P
p=(a)'
During inference, we initialize at the uniform distribution = and run from time 0 to
12. The truncation eπrоr of this choice of time horizon is the magnitude of 10-12 reflected by DKL(P) ( and therefore negligible. The discrete time points form an arithmetic sequence.
We generate 10t samples for each algorithm and use np. bincount obtain the empirical distribu- tion Or as the output distribution. Finally, the KL divergence is computed by
(o DKI(Poll@n Σ po( log (i)
We also perform bootstrapping for 1000 times to obtain the 95% confidence interval the KL divergence, the results are shown by the shaded area in Fig. Q The fitted lines are obtained by standard linear regression on the log-log scale with the slopes marked beside each line Fig.
D.2 Text Generation
For text generation, we use the small version of RADD [X checkpoint trained with X-DCE loss. We choose an early stopping time 10-3 for stable numerical simulation. Since RADD is masked discrete diffusion model, we can freely choose the noise schedule (t) used in the inference process. We consider the following log-linear noise schedule used the model training.
(t)) = 1- (1 - c)t²
-
a(t) - (s)ds = log(1 1011--1-10 c)t)
(D.1)
where we choose e = 10-3
The score function su s2 t) used for computing the transition rate matrix can be computed from the RADD score model py using the following formula from [33]
e (1)
(x:)
1- e- e-(o) (x
(D.2)
where the model Pa is trained to approximate the conditional distribution of the masked positions given all unmasked positions. More specifically, let d be the length of the sequence and 21, be the vocabulary set (not including the mask token). Then given partially masked sequence x = (x',..,, x J. the model (x) outputs a x S matrix whose (f, element approximates XM
Px-paan (x x UM = x UM) when x' is mask, and 1 otherwisse. Here, x UM represents the unmasked portion of the sequence
We adopt uniform discretization of the time interval (8, 1]. For 8-RK-2 and A- @-rapeezoidal, we pick 0 = 1 We compare our proposed 6-RK-2 and 6- Traapezoidal with the Euler method, Tweedie r-leaping. r-lepping, and we present full results across all NFES ranging from 16 1024 in Tabs. and E. For each method, we generate 1024 samples and compute the average perplexities and uni- gram entropy on both GPT-2 large and I.LMMA 3,8 All the experiments are run on single NVIDIA A100 GPU.
44

---

Table Generative perplexity of texts generated by different sampling algorithms on GPT-2 large.
Lower values are better. with the best in bold.
Method
NFE 16
NFE 32
NFE 61
NFE 128
NFE
256
NFE 512
NFE 1024
FHS
S 307.426
< 186.594
< 141.625
122.732
S 113.310
< 113.026
109.406
Fuler
277.962
< 160.586
5. 111.597
$ 86.276
S 68.092
56.622
5 44.686
Tweedie r-leaping
277.133
160.248
110.848
85.738
70.102
160.944
44.257
r-leaping
126.835
< 96.3221
69.226
52.366
41.694
33.7889
28.797
Semi-AR
2857.469
1543.302
$ 741.184
360.793
222.303
1644162
1147.406
B-RK-2
127.363
109.351
86.102
< -61.17
49.816
40.375
-339971
d-Trapezoidal
$ 123.585
S 89.912
I 66.549
S 49.051
S 39.959
I 32.456
S 27.553
Table 4: Generative perplexity of texts generated by different sampling algorithms on LLaMA 3. Lower values better, with the best in bold.
Method
NFE
NFR
NFB
NFE 128
NFE
NFB
NFR 1024
FHS
S 342.498
210.742
$ 155.258
132.133
$ 127.526
$ 123.013
S 120.791
Euler
318
115.555
125.955
S 991.051
75.245
59.971
49.406
Tweedie r-leuping
2 316.744
- 1172.941
$ 1221.248
< $94.253
S 75.403
639.933
S 449.239
T-leuping
S 152.867
S 1117,930
5 $669900
< 68.000
53.664
$. 44.676
S. 238.293
Semi-AR
2696.883
16844.73
251.963
$
1162.93
0-RK-2
150-439
1322.900
107.066
80.742
63.277
52.5633
44.687
0-Trapezoidal
146.0277
< 113.260
< 83.456
66.071
54.307
S 44.2233
35.524
From the results in Tabs. to E. observe that 8-Trapezoidal almost outperforms all other ap- proaches and generates samplers with better perplexities across almost all NFEs. We also noticed that both the Buler method and Tweedie r-leaping perform similarly, and are beaten by a large margin by 0-RK-2 and r-leaping.
In Tab. B. present the percentage of positive extrapolated intensities for different algorithms across NFE values. This partially validates the assumption our theoretical analysis (Thms. 5.4 and 5.50 that the intensity remains positive throughout the sampling process.
D.3 Image Generation
For the image generation, we use the checkpoint of MaskGIT [63. [XS] reproduced in P Recall that the MaskGIT is a masked image model which, given partially masked sequence, outputs the conditional distributions of the masked positions given the unmasked portion, just like the model (-(-) in the aforementioned masked text model, RADD. Therefore, by similarly introducing time noise schedule (t) (for which we adopt the same log-linear schedule (D.II) in our experiment), we obtain masked discrete diffusion model akin to the RADD. The score function be computed accordingly using the model output as (12.2).
We choose an early stopping time = 10-3 and adopt uniform discretization of the time inter- val (8. 1] for 0-RK-2, O-Trapezoidal, r-leaping and the Euler method. For parallel decoding, we use linear randomization strategy in the re-masking step and arccos masking scheduler, as recommended in [63]. For each method, we generate 50k samples a class-conditioned way and compute its FID against the validation split of ImageNet. We use classifier-free guidance to enhance generation quality and set the guidance strength to w = 3.
We present the full results for NFE ranging from 4 to 64 in Fig. S All the experiments are run on NVIDIA A100. Notably, 8-Trapezoidal with 0 for extremely low NFE budgets. While 0-Trapezoidal with 0 1 in general demonstrates a less competitive performance, it converges to the same generation quality as = in the high NFE regime. We also noticed that when using extrapolation with f = 3. 6-RK-2 beats r-leaping for NFE larger than 8, which again accords with our theoretical prediction of its competitive performance in 0 E (0, regime.
is the best- best-performing method except
45

---

Table 5: Unigram entropy of texts generated by different sampling algorithms on LLaMA
Method
NFF NFF: 32
NFF.
NFF 128
NFF 256
NFR 512
NFF. 1024
FHS
7.843
7.793
7.748
7.712
7.714
7.716
7.717
Fuler
7.785
T.677
7.594
7.446
7.343
7.158
6.962
Tweedie r-leaping
7.786
7.675
7.564
7.453
7.345
7.151
6.970
7-lcaping
7.048
7.122
7.016
6.890
6.706
6.537
6.407
Semi-AR
8.019
8.056
7.994
7.908
7.836
7.771
7.810
0-RK-2
6.772
7.017
7.085
7.010
6.831
6.682
6.548
e-Trapezoidal
7.126
7.163
7.033
6.919
6.740
6.532
6.412
Table 6: Percentage of positive extrapolated intensities for different algorithms across NFE values.
Method
NFE 32
NFE 64
NFE 128
NFE 256
NFE 512
NFE 1024
0-RK-2
97.21 + 3.1
98.31 1: 2.0
98.01 + 1.3
99.27 = 0.9
99.44 1 0.7
99.52 = 0.6
0 o-rapezzoiaal
95.67 1 4.8
97.06 1 3.6
98.22 + 2.4
98.87 = 1.6
99.24 1 1.1
99.43 0.9
Zoom $n
11.50
FHiS
10.20
Edler
Para Decoding
Traprsioat 121
21.01
7.29
Reachet Distance
6,50
36.35
12.73
9.91
6.00
Number of Function Evalucaioon
Figure 5: FID of images generated by sampling algorithms vs. number function evaluations (NFE) with different parameter choices. Lower values are better.
To investigate the robustness of 0-RK-2 with respect to the choice of 0, we also benchmark its performance across multiple choices at NFE 32 and 64, and we present the results in Fig. G. We observe that the performance of 9-RK-2 has a flat landscape around the optimal 0 choices, which fall in the range [0.15, 0.4]. In general, shown by the curve, the method performs better when extrapolation is used to compute the transition rate matrix, confirming the correctness of our theory (Thm. 5.50 and our discussions. Similar to the behavior of 0-Trapezoidal method in Fig. EL the performance of A-RK-2 has a flat landscape around the optimal 9 choices, which typically falls in the range [0.3, 0.5). In general, shown by the curve, both methods exhibit better performances when extrapolations are deployed, which, once again, certifies the validity of our theoretical results.
Finally, we visualize some images generated with 0- o-Trapezoidal on 6 different classes in Fig. Ø 0- Trapezoidal consistently generates high-fidelity images that are visually similar to the ground truth ones and well aligned with the concept.
46

---

30
0-RK-2 (NFE 32)
88
6-RK-2 (NFE
64)
86-
84
82
80-
7.8
7.6
74
0.10
0.20
0.30
0.40
0.60
0.70
0.80
0.90
1.00
Figure 6: Sampling quality E (0, 1] in 0-RK-2 algorithm. Sampling quality is quantified through FID.
Figure 7: Visualization of samples generated by 6- Trappezoida.. Upper Left: Aircraft carrier
(ImageNet-Ik 1k class: 933); Upper Middle: Pirate (ImageNet-Ik lk class: 724); Upper Right: Vol- cano (ImageNet-Ik 1k class: 980); Lower Left: Ostrich (ImageNet- lk class: 009); Lower Middle:
Cheeseburger (ImageNet-14 1k class: 933); Lower Right: Beer bottle (ImageNet-Il 1k class: 440).
D.4
Diffusion Large Language Model and Math Reasoning
the checkpoint of LLaDA [64] after instruction tuning. LLaDA-Insstre 8B one of the otate-of-the-art diffusion-based LLMs trained natively based masked discrete diffusion models [33]]. For the math reasoning datasets, consider GSM8K,5 B standard benchmark for math reasoning consisting of clementary-scchol-eeee word problems. Similar RADD, since LLaDA also masked discrete diffusion model, can freely choose the noise schedule (()) used in the inference process. We adopt the same log-linear schedule described in (D.I), and choose 10-3 The computation of the score function using LLaDA follows the same procedure as is depicted in App.
For the evaluation of diffusion LLMs math reasoning datasets,
47

---

We choose an early stopping time = 10 10m3, and adopt a uniform discretization of the time interval (6, 1] for B-Trapezoidal, with 0 3. For Semi-AR with random remasking confidence-basee remasking strategies, we set the block length to 256, the same the generation length, to avoid exploiting LLaDA's inherent preference for auto-regressive generation order and ensure fair com- parison with our proposed method. However, we note that our proposed method can also be gener- alized to a blockwise sampling setting, and we include additional results here for proof of concept and to demonstrate that LLaDA indeed prefers a left-to-right generation order. We set the sampling temperature to 0.5 for evaluated inference methods. The evaluation of generated responses fol- lows a standard pipeline with in-evaluation harness, a standard LLM evaluation kit, following the instruction in the LLaDA repository. & Full results are presented in Tab. 7. All experiments are run on NVIDIA A100.
Table 7: Response accuracy on GSM8K with different NFEs.
Accuracy (%)
NFE 64
NFE 128
NFE 256
Semi-AR (Conf.)
33.6
32.0
39.1
Semi-AR (Rand.)
33.8
34.3
40.3
0-Trapezoidal
35.1
38.4
39.7
0-Trapezoidal (Blocksize 8)
50.4
57.4
62.5
48

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract and introduction accurately reflect the paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
The paper discusses the limitations of the work performed by the authors.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are discussed in the paper.
The authors are encouraged to create separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The au- thors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, c.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the ap- proach. For example, facial recognition algorithm may perform poorly when image resolution is low images taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and fairesss.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't it acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
49

---

Justification:
The paper provides the full set of assumptions and complete (and correct)
proof.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theo- rems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The paper fully discloses all the information needed to reproduce the main experimental results of the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or the contribution it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, other means that are appropriate the research performed.
specific model and empirical evaluation,
While NeurIPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results way to re- produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
50

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
The paper opens access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as deseribed in supplemental material.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (ittpa:)
//nips. cttccc/pp. htcp/pbllccgguiii hwwli//ggie/////
ssiionolicy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: The paper specifies all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: The paper reports error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
51

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the enror bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should prefer- ably report 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative еπτог rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
The paper provides sufficient information on the computer resources (type of compute workers, memory, time of execution) needed reproduce the experiments.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: The research conducted in the paper conforms, in every respect, with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., s if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There is no societal impact of the work performed.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
52

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), faimess considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: There is no such risk.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by re- quiring that users adhere usage guidelines restrictions access the model or implementing safery filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The creators original owners of assets, used in the paper, properly credited and the license and terms of use are explicitly mentioned and properly respected.
Guidelines:
The answer NA means that the paper does
existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
53

---

If assets are released, the license, copyright information, and terms of use in the pack- age should be provided. For popular datasets, paperswithcode com/datasets has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documenta- tion provided alongside the assets?
Answer: [Yes]
Justification: The assets introduced in the paper well documented and the docu- mentation is provided alongside the assets.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/model part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can cither create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the pa- per include the full text instructions given to participants and screenshots, if applicable, well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information the supplemental material is fine, but the main contri- bution of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, cura- tion, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsoureing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equiva- lent) may be required for any human subjects research. you obtained IRB approval, you should clearly state this in the paper.
54

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The paper does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy ltteas//Doariiscccc///ofeerrrerp for what should or should not be described.
55


Reviews:
[{"review_id": "nuhv2DbeDF", "reviewer": "Reviewer_61BM", "summary": "This work addresses  computational bottleneck in discrete diffusion models by developing higher-order solvers that significantly reduce the number of function evaluations (NFE) required for sampling. The authors introduce theoretical foundations for higher-order methods in the discrete setting and demonstrate substantial improvements in sampling efficiency for text generation tasks.", "questions": "Following Weakness\n\nIt would be informative for readers and broader community to see integration of higher order solvers with\n- Effect of Temperature/Confidence\n- Multi-token unmasking heuristics like in FastDLLMs, etc.\n- ReMasking\n- Predictor-Corrector Solvers\n- How well higher order solvers work with distilled discrete diffusion models based on e.g., SDTT Checkpoints\n\nAlso additional metrics like use state-of-the-art language models (LLaMA-3 or alternatives) to evaluate generated text quality independently of the diffusion model.\n\n**References** \n\n1. G. Wang et al. Remasking Discrete Diffusion Models with Inference-Time Scaling (https://arxiv.org/abs/2503.00307)\n2. J. Deschenaux et al. Beyond Autoregression: Fast LLMs via Self-Distillation Through Time (https://arxiv.org/html/2410.21035v1)\n3. C. Wu et al. Fast-dLLM: Training-free Acceleration of Diffusion LLM by Enabling KV Cache and Parallel Decoding (https://arxiv.org/html/2505.22618v1)", "limitations": "N/A", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n* The paper provides mathematical foundations for higher-order solvers in discrete diffusion models, extending beyond the first-order methods commonly used in practice.\n* The paper demonstrates meaningful reductions in NFE while maintaining or improving generation quality, addressing a key limitation of discrete diffusion models.\n* The paper is well-written, clearly organized, and accessible to readers.\n\n\n**Weakness**\nLimited Scope of Experimental Evaluation\n\nThe paper lacks comprehensive analysis of how different solver configurations perform under varying conditions, limiting understanding of when and why higher-order methods excel. This limits drawing insights and adoption guidelines for readers and broader community \n\n- Currently this work does not compare against several important recent methods and how effective higher order solvers are across these settings.\n- While the appendix mentions predictor-corrector solvers, the authors neither evaluate them experimentally nor discuss how higher-order methods could be integrated within predictor-corrector frameworks. \n- Remasking during sampling shows enhanced generation quality. Would be interesting to see higher order solvers with Remasking like in works like ReMDM\n- Absence of semantic quality metrics: No evaluation using ModernBERT-based MAUVE scores or similar metrics that assess semantic coherence and diversity (entropy)\n- Minimal Performance gain in case of Image Generation (but didn't consider this much for rating)", "was_revised": true, "final_justification": "I retain my original rating, while agree on the contribution of formulation of higher order solvers for discrete diffusion. \n\nHigher order integration/ODE solvers is well established literature with close relation to SGD variants too. \n\nWhile the work is valuable the contribution would be significantly enhanced and useful if we can understand implications of higher order solvers in discrete diffusion setting which authors do not address in anyway and largely adopt Numerical Integration convergence guarantees when applied to discrete diffusion setting", "rebuttal": "Thank you for your constructive review and recognition of our mathematical foundations, meaningful NFE reductions while maintaining generation quality, and clear presentation. We're particularly encouraged by your acknowledgment that our work addresses key limitations in discrete diffusion models and provides accessible theoretical insights.\n\nBeyond empirical results, our work establishes the first rigorous mathematical framework for high-order numerical methods in discrete diffusion inference, achieving provably second-order convergence guarantees that were previously unknown in this domain. Our $\\theta$-trapezoidal method maintains unconditional second-order accuracy across all $\\theta \\in (0,1]$, representing a fundamental theoretical advance over existing first-order $\\tau$-leaping methods. These theoretical contributions are validated through empirical evaluation across multiple domains, demonstrating that rigorous mathematical foundations translate directly into practical performance gains and open new research directions for discrete diffusion acceleration.\n\n---\n\n### Regarding the evaluation scope and integration with recent methods\n\nThank you for raising this important point about expanding experimental scope. We acknowledge that more comprehensive analysis across different solver configurations and integration with recent techniques would strengthen our work's practical impact. Our current work establishes fundamental theoretical frameworks for high-order methods in discrete diffusion models, with experiments strategically designed to validate theoretical predictions across multiple domains.\n\nRegarding the specific integration suggestions:\n- *Predictor-corrector frameworks*: Our high-order methods can naturally serve as enhanced predictor steps, with $\\theta$-trapezoidal's extrapolation nature particularly well-suited for this integration. The theoretical benefits are also interesting directions, similar to continuous counterparts [1].\n- *Remasking techniques*: ReMDM's remasking technique is intended as inference-time scaling for generation quality enhancement, where computation budget such as NFE is not a major concern. Our high-order methods serve a completely different purpose: reducing inference NFE while maintaining generation quality. Since ReMDM is concurrent with our submission, we'll investigate this integration in future work.\n- *Temperature and confidence*: Temperature is orthogonal to our main algorithm design and typically depends heavily on specific models and targeted tasks. For example, LLaDA [2] performs well under low sampling temperature, while DiffuCoder [3] works better under high temperature. We select best-performing temperatures for each model rather than ablating this choice. Confidence is also tangent to our approach as we don't adopt heuristics for determining token generation order.\n- *Multi-token unmasking*: Regarding advanced acceleration techniques such as those in fast-dLLM, incorporating KV-cache could further improve our algorithm's inference FLOPS. However, our focus is restricted to algorithmic level with efficiency quantified by number of function evaluations (NFEs), making fast-dLLM's techniques beyond our current scope. Since fast-dLLM is concurrent with our submission, we'll investigate this integration in future work. \n- *Distilliation*: Our high-order numerical schemes operate at inference level and should be compatible with various architectures and training approaches, including distilled models, as long as the distilled models satisfy our analytical assumptions. \n\nWe appreciate these diverse suggestions and agree they demonstrate intriguing directions for future exploration. We will include discussion of these relevant references and potential integration approaches in our revision to provide a more comprehensive view of the research landscape.\n\n\n### Regarding semantic quality metrics and evaluation depth\n\nThank you for pointing out this important issue. We adopted your advice and use Llama 3 for generative complexity computation, also reporting unigram entropy as a diversity metric. We re-benchmarked RADD performance under this evaluation framework and included LLaDA's semi-AR sampler as an additional benchmark.\n\nDue to time constraints, we report NFE up to 512 with each metric computed based on 512 samples. We will include full experimental results in the paper revision. For the semi-autoregressive (Semi-AR) method, we follow the original implementation proposed in LLaDA using the random remasking strategy. We found that when choosing confidence-remasking strategy for RADD, the generated samples consist of repeated sequences of most frequent tokens with extremely small entropy, so we omit those results.\n\nOur proposed algorithms perform best across nearly all NFEs under the new evaluation metrics. Other benchmarks such as Semi-AR or FHS either produce texts with higher perplexities or unnatural sample entropy, making them inferior in performance compared to $\\theta$-Trapezoidal. This again suggests the empirical advantage of high-order solvers in practice.\n\n\n| Generative ppl with Llama3                 | NFE = 16  | NFE = 32  | NFE = 64 | NFE = 128 | NFE = 256 | NFE = 512 |\n|:----------------------:|:---------:|:---------:|:--------:|:---------:|:---------:|:---------:|\n| FHS                    |342.498    |210.742    |155.258   |132.135    |127.526    |123.013    |\n| Euler                  |318.413    |175.555    |125.955   |91.051     |75.245     |59.971     |\n| Tweedie $\\tau$-leaping |316.744    |172.941    |121.248   |94.253     |75.403     |59.943     |\n| $\\tau$-leaping         |152.867    |117.930    |86.980    |68.090     |**53.664** |44.676     |\n| **$\\theta$-RK-2**          |150.439    |132.090    |107.066   |80.742     |63.277     |52.563     |\n| **$\\theta$-Trapezoidal**   |**146.027**|**113.260**|**83.456**|**66.071** |54.307     |**44.293** |\n| Semi-AR                |2696.883   |1684.973   |829.391   |410.177    |251.963    |183.599    |\n\n| Unigram entropy                 | NFE = 16 | NFE = 32 | NFE = 64 | NFE = 128 | NFE = 256 | NFE = 512 |\n|:----------------------:|:--------:|:--------:|:--------:|:---------:|:---------:|:---------:|\n| FHS                    |7.843     |7.793     |7.748     |7.712      |7.714      |7.716      |\n| Euler                  |7.785     |7.677     |7.594     |7.446      |7.343      |7.158      |\n| Tweedie $\\tau$-leaping |7.786     |7.675     |7.564     |7.453      |7.345      |7.151      |\n| $\\tau$-leaping         |7.048     |7.122     |7.016     |6.890      |6.706      |6.537      |\n| **$\\theta$-RK-2**          |6.772     |7.017     |7.085     |7.010      |6.831      |6.682      |\n| **$\\theta$-Trapezoidal**   |7.126     |7.163     |7.033     |6.919      |6.740      |6.532      |\n| Semi-AR                |8.019     |8.056     |7.994     |7.908      |7.836      |7.771      |\n\nAdditionally, we benchmarked our high-order solver on GSM8K with LLaDA to enhance empirical validation. Please refer to our response to Reviewer 5vgu for detailed results. We hope these additional experiments could alleviate your concerns about evaluation scope.\n\n### Regarding performance gains in image generation\n\nThe smaller gains in image generation compared to text reflect different domain characteristics. However, $\\theta$-trapezoidal consistently outperforms existing approaches across nearly all tested NFE budgets. Our theoretical analysis suggests high-order methods' benefits become more pronounced with accurate score estimation and longer time horizons, potentially explaining domain-specific improvement variations.\n\n---\n\nWe sincerely appreciate your detailed feedback and constructive suggestions. Our work provides a two-fold contribution to the field: establishing rigorous theoretical foundations for high-order methods in discrete diffusion models and demonstrating their practical effectiveness across multiple domains. The theoretical framework we establish, combined with empirical validation using state-of-the-art evaluation models, opens new avenues for efficient discrete diffusion inference that complement existing acceleration approaches.\n\nThe enhancements you suggest would indeed strengthen the practical applicability of our work, and we're committed to incorporating these improvements in future iterations. Given the novelty of our theoretical contributions, the consistent empirical improvements we demonstrate across different evaluation frameworks, and our willingness to address the evaluation scope concerns, we would be grateful if you would reconsider your evaluation of our work.\n\n---\n### References\n\n[1] Chen, Sitan, et al. \"The probability flow ode is provably fast.\" Advances in Neural Information Processing Systems 36 (2023): 68552-68575.\\\n[2] Nie, Shen, et al. Large Language Diffusion Models. arXiv:2502.09992.\\\n[3] Gong, Shansan, et al. DiffuCoder: Understanding and Improving Masked Diffusion Models for Code Generation. arXiv:2506.20639"}, {"review_id": "9AqsFrPqqX", "reviewer": "Reviewer_HxA8", "summary": "This paper introduces higher order methods for discrete diffusion models, enabling time stepping schemes with higher order convergence when compared to standard schemes. Explicit examples of higher order methods, such as RK-2 and Trapezoidal are given. Computational examples are provided.", "questions": "Can the techniques presented here be used to create even higher order methods, such as 3rd or fourth order methods? \n\nCan the authors comment on the influence that the noise schedule has on their method? Usually the constant on the order of a high order scheme has some dependence on the gradient of the function to be integrated. By choice of a clever schedule, would one be able to also help diminish the size of this constant.\n\nDo classical time stepping stability analysis apply? For standard Runge-Kutta methods, we can look at regions of stability, can we do anything similar here?", "limitations": "There is no limitations section of this paper, I would recommend that the authors add this to their submission.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\nThis paper is well motivated, it is true in the continuous setting that higher order methods have traditionally given a computational gain. As the discrete time diffusion is still continuous in its temporal component, this is a natural extension that should be explored. \n\nThe paper is well layed out and easy to follow. The paper starts with defining discrete diffusion models, shows current stepping methods and presents Theorem 3.1 revealing standard methods are first order. The authors then present their new schemes and Theorem 5.4 and 5.5, natural extensions of Theorem 3.1. \n\nThe paper well states its assumptions, with Assumptions 5.1,5.2 and 5.3.  In particular, it is nice that Assumption 5.2 is directly stated, showing that for their proposed higher order scheme to work, they need 2nd order regularity in the rate matrix. \n\nWeaknesses:\nI do not see any major weaknesses. It would be nice if some standard techniques from higher order schemes were also presented, such as Buchner tables for the construction of these methods if applicable.", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thorough and constructive review of our paper. We greatly appreciate your positive assessment of our work as well-motivated and well-laid out. We are particularly pleased that you found our theoretical contributions, including the second-order convergence results and the explicit statement of our assumptions, to be clearly presented. Your insightful questions and suggestions will help us improve the paper further.\n\n---\n\n### Regarding standard techniques from higher-order schemes\n\nThank you for this valuable suggestion. We agree that including discussion of standard techniques such as Butcher tableaux would enhance the paper's completeness. The reviewer's suggestion coincides with a classical work [1] that uses Butcher tableaux to derive various Poisson Runge-Kutta methods. In our revision, we will add a brief discussion in Section 4 explaining how these techniques relate to our methods and how certain modifications can be carried out for the discrete diffusion setting. \n\n### Regarding higher-order methods ($\\geq$ 3rd order)\n\nThank you for this insightful question. Indeed, constructing even higher-order methods is possible. As shown in the work [2] on numerical analysis of inhomogeneous continuous-time Markov chains, higher-order schemes can be developed. However, there are several considerations, including the smoothness requirement, score estimation accuracy, and computational trade-off, that may require case-by-case empirical studies. We will add a discussion of these points in Section 7 (Future Work) of our revision, noting this as an exciting direction for future research.\n\n### Regarding the influence of noise schedule\n\nThank you for this excellent observation. You are absolutely correct that the constant in our convergence bounds depends on the noise schedule and time parametrization. This dependency enters through the smoothness properties of the intensity functions (Assumption 5.2) and affects the constants hidden in our $\\mathcal{O}(\\kappa^2 T)$ bounds.\n\nAs noted in the continuous diffusion literature, particularly the work [3], the choice of noise schedule significantly impacts both training and inference. For discrete diffusion models, this influence is even more pronounced due to the discrete state space. The optimal noise schedule likely depends on the properties of the data distribution and the specific discrete state space structure. \n\nWe will add a remark after Theorem 5.4 discussing how the noise schedule affects the practical constants in our convergence bounds and reference relevant empirical studies. This is indeed an active area of research that deserves further investigation.\n\n### Regarding classical time-stepping stability analysis\n\nThank you for raising this important theoretical question. To the best of our knowledge, the classical notion of stability regions from deterministic numerical analysis does not directly translate to our stochastic setting. Instead of absolute stability regions, the relevant concepts for stochastic numerical schemes include weak convergence and moment bounds, and our analysis follows the framework established in the stochastic numerical analysis literature, such as the work [3] on tau-leaping methods. We will add a remark in Section 5.2 clarifying this distinction.\n\n### Regarding the limitations section\n\nThank you for pointing out this omission. We will add a dedicated Limitations section before Conclusion, discussing our theoretical assumptions, memory requirements for our proposed inference schemes, the dependence of practical performance on score accuracy, and current restriction to second-order methods, *etc..*\n\n---\n\nWe sincerely thank you again for your careful review and valuable feedback. Your questions have helped us identify several areas where we can clarify and strengthen our presentation. We will incorporate all the suggested improvements in our revision and believe these additions will make our paper more complete and accessible to readers.\n\n---\n### References\n\n[1] Burrage, Kevin, and Tianhai Tian. \"Poisson Runge-Kutta methods for chemical reaction systems.\" (2004): 82-96.\\\n[2] M. Arns, P. Buchholz, A. Panchenko, (2009) On the Numerical Analysis of Inhomogeneous Continuous-Time Markov Chains. INFORMS Journal on Computing 22(3):416-432.\\\n[3] Karras, Tero, et al. \"Elucidating the design space of diffusion-based generative models.\" Advances in neural information processing systems 35 (2022): 26565-26577.\\\n[4] Hu, Yucheng, Tiejun Li, and Bin Min. \"A weak second order tau-leaping method for chemical kinetic systems.\" The Journal of chemical physics 135.2 (2011)."}, {"review_id": "A3p2jpcwKO", "reviewer": "Reviewer_LqdK", "summary": "This paper addresses the challenge of slow sampling in discrete diffusion models. The authors propose and analyze several high-order sampling algorithms, which are novel extensions of the tau-leaping method. The primary goal is to accelerate the generation process, drawing a parallel to the successful application of high-order solvers for continuous diffusion models.", "questions": "1. Regarding Theorem 5.4, the error bound seems to suggest that performance could be further improved by setting the NFE higher than the sequence length (L). Is this a correct interpretation? If so, does it make practical sense to use more sampling steps than the number of tokens to be generated?\n2. A critical point of comparison arises at NFE = L. For a simple masking-based sampler like the one in LLaDA, setting NFE = L (i.e., unmasking one token per step) is analogous to autoregressive decoding. This process should, in principle, allow for exact sampling from the model's conditional distributions, resulting in zero sampling error. In this specific and practical regime: What are the theoretical and empirical advantages of using the proposed high-order methods over such a simple, exact-sampling baseline?", "limitations": "A key limitation is the lack of comparison to simpler, heuristic-based sampling algorithms (e.g., LLaDA's remasking). The practical benefits of the proposed high-order methods are not clearly established in regimes where these simpler methods might perform optimally (such as when NFE = Sequence Length), which could limit the adoption of this more complex approach.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper tackles a conceptually important and underexplored problem. While accelerating sampling for continuous diffusion models is a mature field, this work is one of the first to provide a rigorous exploration of high-order methods for the discrete domain.\n2. The theoretical analysis is a core strength. The proposed high-order variants of the tau-leaping method are well-motivated, and their derivation and convergence analysis are presented with mathematical rigor.\n3. The experiments are comprehensive and sufficiently validate the efficacy of the proposed algorithms. The results clearly demonstrate an acceleration in the sampling process when compared to baseline first-order methods across several tasks.\nWeakness: \nI think the paper misses a direct comparison to conceptually simpler methods like the confidence-ranked remasking technique proposed in LLaDA (Arxiv: 2502.09992). Such methods, while perhaps less theoretically grounded in diffusion theory, are highly practical and serve as a crucial baseline. A thorough analysis, particularly in the important inference scenario where the Number of Function Evaluations (NFE) equals the sequence length, is needed to demonstrate the practical advantages of this work's more complex approach over the strong, simple heuristics in LLaDA. See also my question below.", "was_revised": true, "final_justification": "I have read the author's response and I am satisfied with the clarifications. I will be keeping my score and recommending acceptance.", "rebuttal": "Thank you for your thorough and insightful review of our paper. We greatly appreciate your recognition of our theoretical contributions and the importance of addressing the underexplored problem of accelerating discrete diffusion model sampling. Your constructive feedback, particularly regarding comparisons with simpler methods and the practical regime where NFE equals sequence length, raises important points that will help strengthen our work.\n\n---\n\n### Regarding comparison with confidence-ranked remasking methods (*e.g.*, LLaDA)\n\nThank you for highlighting this important baseline. We want to clarify that we have indeed included such a comparison in our experiments through the \"parallel decoding\" baseline, which is conceptually equivalent to the confidence-ranked remasking technique proposed in LLaDA. Specifically, our parallel decoding implementation follows the same principle: at each step, we compute the logits for all masked tokens and then remask the tokens with the lowest confidence scores (see lines 1461-1465 in our paper).\n\nTo further resolve this issue, we have carried out additional experiments with LLaDA. We refer the reviewer to our response to Reviewer 5vgu for more results on confidence-ranked methods on RADD and LLaDA. The additional results confirm again the superior performance of our proposed algorithms. These experiments reflects that our proposed high-order solver consistently exhibits competitive performance across tasks of various type and models of various size, while other sampler heuristics such as semi-AR sampler might suffer from instability and have failure cases (such as on RADD).\n\nWe will incorporate these new experiments into our revision and add a clarification in the paper on these samplers. We will also expand our discussion to better contextualize how our high-order methods compare to these confidence-based heuristics in different sampling regimes.\n\n### Regarding NFE higher than sequence length (Theorem 5.4)\n\nThank you for this astute observation about our theoretical results. You are correct that Theorem 5.4 analyzes uniform discrete diffusion models where using more steps can indeed reduce discretization error. However, it's important to distinguish between two types of discrete diffusion models:\n- For *uniform* discrete diffusion models, the forward process converges to a uniform distribution over the entire state space, and multiple transitions per position are possible. In this setting, using $\\mathrm{NFE} > L$ can provide benefits by allowing finer time discretization.\n- For *masked* discrete diffusion models (like those used in our text and image experiments), each position typically transitions only once from masked to unmasked. Therefore, as you correctly note, using $\\mathrm{NFE} > L$ provides no additional benefit. Our experiments on text and image generation focus on the practical regime where $\\mathrm{NFE} \\leq L$, which is why we report results up to $\\mathrm{NFE} = 1024$ for sequence length 1024.\n\n\nWe will clarify this distinction more prominently in the revision to avoid confusion.\n\n### Regarding advantages over exact sampling when $\\text{NFE} = L$\n\nThis is an excellent question that goes to the heart of our contribution. While it's true that methods like LLaDA can perform exact sampling when $\\text{NFE} = L$ (unmasking one token per step), there are several important considerations that favor our high-order approach:\n\n1. The \"exact sampling\" property assumes perfect score function estimation. In practice, neural networks provide imperfect score estimates, and our high-order methods can better handle these estimation errors through their numerical properties. As discussed in Section 3.1, exact simulation methods suffer from highly skewed computational effort, with most jumps concentrated near the end of the generation process when score estimation errors are highest.\n2. Our methods provide flexibility in the NFE budget. While $\\text{NFE} = L$ might seem natural, practitioners often want to trade off quality for speed ($\\text{NFE} < L$) or achieve higher quality with more computation ($\\text{NFE} \\rightarrow L$). Our theoretical analysis shows that high-order methods achieve better quality-speed trade-offs across this entire spectrum, not just at the specific point $\\text{NFE} = L$.\n3. Our empirical results demonstrate practical advantages even when comparing against confidence-based methods. In Figures 4 and 7, the $\\theta$-trapezoidal method consistently outperforms parallel decoding (equivalent to LLaDA-style remasking) across various NFE budgets, including those approaching the sequence length.\n\nWe will expand our discussion in Section 3.1 to better articulate these advantages and add a more detailed comparison specifically at the $\\text{NFE} = L$ regime in our revision.\n\n---\n\nThank you again for your valuable feedback and thought-provoking questions. Your comments have helped us identify areas where we can clarify our contributions and better position our work relative to existing methods. We believe that addressing these points will significantly strengthen our paper and make it more useful to practitioners in the field. We hope our responses have addressed your concerns."}, {"review_id": "oK7PcK2jt2", "reviewer": "Reviewer_5vgu", "summary": "The paper introduces high-order numerical solvers (θ-RK-2 and θ-trapezoidal methods) for discrete diffusion models, aiming to accelerate sampling and inference compared to standard τ-leaping and exact simulation. Theoretical analysis demonstrates second-order convergence in KL divergence under certain conditions. Empirical results on text and image generation show better sample quality for the proposed methods.", "questions": "- Can the authors provide evidence that their methods remain effective on larger, more complex benchmarks and with a broader set of evaluation metrics?\n- How do the methods perform when the theoretical assumptions are not strictly satisfied? If the authors can address such problem, I am willing to improve my rating.", "limitations": "The main limitations regarding assumptions and empirical scope are discussed by the authors.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths and Weaknesses:  \n\nStrengths\n- The paper provides a rigorous theoretical analysis, with detailed stochastic integral formulations and error bounds.\n- The introduction of high-order schemes to discrete diffusion inference is novel and addresses an under-explored topic.\n- Empirical validation is thorough, with experiments on both text and image benchmarks.\n- The proposed methods are conceptually simple and easy to incorporate into existing diffusion frameworks.\n\nWeaknesses  \n- The experimental validation is limited: most results are on toy or small-scale models, with no evaluation on large-scale or widely recognized datasets (e.g., Geneval [1], MMMU [2]). The base models are relatively simple, and only basic baselines are considered.\n- Evaluation metrics are primarily FID and NFE, which may not be sufficient to fully support the claimed effectiveness of the methods.\n- Theoretical guarantees depend on strong assumptions (regularity, boundedness, positivity, accurate score estimation) that are difficult to ensure in practical scenarios. The robustness of the proposed methods when these conditions are not met is not addressed.\n\n[1] GenEval: An Object-Focused Framework for Evaluating Text-to-Image Alignment\n\n[2] MMMU: A Massive Multi-discipline Multimodal Understanding and Reasoning Benchmark for Expert AGI", "was_revised": true, "final_justification": "I have read the rebuttal. The authors gave a very detailed reply that addressed most of my concerns. I raise the rating to borderline accept after the discussion. However, I still suggest the authors conduct evaluation on the GenEval in the future.", "rebuttal": "Thank you for your constructive feedback. We appreciate your recognition of our rigorous theoretical analysis and novelty of introducing high-order schemes to discrete diffusion inference. We address your specific concerns below.\n\n---\n\n### Regarding experimental validation on larger benchmarks\n\nWe thank the reviewer for this comment. As our work presents a two-fold contribution combining rigorous theoretical analysis with empirical validation, our benchmarks were strategically selected to provide both rigorous theoretical verification and practical performance evaluation: The 15-dimensional toy model allows exact KL divergence to verify theoretical convergence rates, while the GPT-2 text generation and ImageNet image generation represent standard literature benchmarks that enable direct comparison with existing methods. \n\nTo enhance the empirical validation of the proposed method, we additionally benchmark its performance on GSM8K (a dataset of grade school math) with LLaDA [4], an 8B paramter masked diffusion language model. For each method, we generate a response of 256 tokens without using demonstration examples in the context, and the answer is evaluated with LM-evaluation-harness, a standard LLM evaluation kit. Due to time limit, we only compare $\\theta$-Trapezoidal to the semi-autoregressive sampler therein, considering both confidence-based remasking strategy and purely random remasking strategy. The results are reported below:\n\n| GSM8K Acc                 | NFE = 64  | NFE = 128  | NFE = 256 |\n|:----------------------:|:---------:|:---------:|:--------:|\n| Semi-AR (Confidence)    |33.6    |32.0   |39.1   |\n| Semi-AR  (Random)       |33.8    |34.3   |**40.3**|\n| **$\\theta$-Trapezoidal**   |**35.1**|**38.4**|39.7|\n\n$\\theta$-Trapezoidal outperforms the Semi-AR sampler in the low NFE regime, where NFE is strictly smaller than the sequence length. At high NFE regime, $\\theta$-Trapezoidal exhibits a similarly competitive performance as other solvers. This observation accords with our claim that high-order samplers perform better with lower NFE budgets. \n\nMoreover, the Semi-AR sampler, while performing decently on LLaDA, is a heurstic sampler not universally good-performing. As seen from the results in the reponse below, when apply to RADD (another diffusion LLM), semi-AR significantly underperforms, while our propsoed $\\theta$-Trapezoidal demonstrates stable, superior performance across models and NFEs, likely due to its theoretically-grounded effectiveness. \n\nWe hope this diminishes your concern on the scalability of our method. Exploring the practical performance on larger-scale text-to-image discrete diffusion models, *e.g.*, GenEval, involves substantial engineering efforts and computational resources beyond this rebuttal period. However, we view scaling to other larger models remains important future work building upon our theoretical foundations.\n\n### Regarding evaluation metrics\n\nThank you for pointing this out. To make the reported results more trustworthy, we additionally report the generative perplexity of RADD measured by Llama 3 (8B parameters in contrast to the 700M GPT-2 large), as well as the corresponding unigram entropy. We also include the Semi-AR sampler proposed in the LLaDA paper [4] as an additional benchmark. \n\nDue to time constraints, we only report NFE up to 512, and each metric is computed based on 512 samples. We will include full experiment results in the revision of the paper. For Semi-AR, we follow the orginal implementation proposed in LLaDA and use the random remasking strategy. We found that when choosing confidence-remasking strategy for RADD, the generate samples are consisted of repeated sequences of most frequent tokens, with extremely small entropy, and we thus omit the results in the table.\n\nFrom the table, we see that our algorithms still performs the best across nearly all NFEs under the new evaluation metric. Other benchmarks such as Semi-AR or FHS either produces texts with higher perplexities or unnatural sample entropy, making them inferior in performance. \n\n\n| Generative ppl with Llama3                 | NFE = 16  | NFE = 32  | NFE = 64 | NFE = 128 | NFE = 256 | NFE = 512 |\n|:----------------------:|:---------:|:---------:|:--------:|:---------:|:---------:|:---------:|\n| FHS                    |342.498    |210.742    |155.258   |132.135    |127.526    |123.013    |\n| Euler                  |318.413    |175.555    |125.955   |91.051     |75.245     |59.971     |\n| Tweedie $\\tau$-leaping |316.744    |172.941    |121.248   |94.253     |75.403     |59.943     |\n| $\\tau$-leaping         |152.867    |117.930    |86.980    |68.090     |**53.664** |44.676     |\n| **$\\theta$-RK-2**          |150.439    |132.090    |107.066   |80.742     |63.277     |52.563     |\n| **$\\theta$-Trapezoidal**   |**146.027**|**113.260**|**83.456**|**66.071** |54.307     |**44.293** |\n| Semi-AR                |2696.883   |1684.973   |829.391   |410.177    |251.963    |183.599    |\n\n| Unigram entropy                 | NFE = 16 | NFE = 32 | NFE = 64 | NFE = 128 | NFE = 256 | NFE = 512 |\n|:----------------------:|:--------:|:--------:|:--------:|:---------:|:---------:|:---------:|\n| FHS                    |7.843     |7.793     |7.748     |7.712      |7.714      |7.716      |\n| Euler                  |7.785     |7.677     |7.594     |7.446      |7.343      |7.158      |\n| Tweedie $\\tau$-leaping |7.786     |7.675     |7.564     |7.453      |7.345      |7.151      |\n| $\\tau$-leaping         |7.048     |7.122     |7.016     |6.890      |6.706      |6.537      |\n| **$\\theta$-RK-2**          |6.772     |7.017     |7.085     |7.010      |6.831      |6.682      |\n| **$\\theta$-Trapezoidal**  |7.126     |7.163     |7.033     |6.919      |6.740      |6.532      |\n| Semi-AR                |8.019     |8.056     |7.994     |7.908      |7.836      |7.771      |\n\n\n### Regarding theoretical assumptions and robustness\n\nThank you for raising these important questions regarding the assumptions underlying our theoretical analysis. We provide detailed responses below addressing the applicability and practical implications of each assumption:\n\n- *Regularity*: Our pursuit of second-order convergence naturally requires second-order continuity of the intensities for rigorous numerical analysis. This smoothness is typically guaranteed under standard forward process choices. However, controlling the Hölder constants involves complex interactions between data properties, time parameterization, and noise schedules—an active research area following empirical studies in the continuous case [1].\n- *Boundedness*: Given our early stopping scheme at time $T-\\delta$, boundedness of intensities is ensured under common forward process choices, particularly when the data distribution has a lower bound. This assumption holds in most practical scenarios and can be readily enforced through truncation and appropriate choice of $\\delta$.\n- *Positivity*: This assumption is standard in stochastic analysis literature [2], arising from the extrapolation nature of our scheme. The violation probability has been proven asymptotically small to arbitrary order [3]. We explicitly evaluated this probability for discrete diffusion models, with results summarized below:\n\n    | Method | NFE=64 | NFE=128 |  NFE=256 | NFE=512 |\n    |--------|-----------|-----------|-----------|------------|\n    | $\\theta$-RK-2 | 98.31 ± 2.0 | 98.01 ± 1.3 | 99.27 ± 0.9 | 99.44 ± 0.7  |\n    | $\\theta$-Trapezoidal | 97.06 ± 3.6 | 98.22 ± 2.4 | 98.87 ± 1.6 | 99.24 ± 1.1 |\n\n    These results demonstrate that positivity holds approximately 98-99% of the time. In practical scenarios, violations can be handled through simple cut-off without affecting performance, as confirmed by our empirical results.\n\n- *Score Accuracy*: Like all inference algorithms, our method depends on accurate score estimation. Generated samples contain three error sources: truncation (exponentially small), discretization (which our algorithm minimizes), and estimation (training-dependent). Our observations confirm that more accurate score estimation leads to more effective inference, approaching theoretical second-order convergence.\n\nOur empirical evaluation demonstrates robust performance across different settings, validating the practical effectiveness of our theoretical framework under realistic conditions. We will polish the paper to provide clearer practical guidance for implementation, including truncation strategies for positivity and early stopping for regularity. \n\n---\n\nWe appreciate your detailed feedback and hope our responses address your concerns. Our work establishes rigorous theoretical foundations for high-order methods in discrete diffusion models while demonstrating practical effectiveness through empirical validation. Our theoretical framework advances convergence understanding, while empirical results across multiple domains confirm practical utility. Given these clarifications and additional experiments, we would be grateful if you would reconsider your evaluation.\n\n---\n\n### References\n\n[1] Karras, Tero, et al. \"Elucidating the design space of diffusion-based generative models.\" Advances in neural information processing systems 35 (2022): 26565-26577.\\\n[2] David F Anderson and Jonathan C Mattingly. A weak trapezoidal method for a class of stochastic differential equations. Communications in Mathematical Sciences, 9(1):301–318, 2011. \\\n[3] Hu, Yucheng, Tiejun Li, and Bin Min. \"A weak second order tau-leaping method for chemical kinetic systems.\" The Journal of chemical physics 135.2 (2011).\n[4] Nie, Shen, et al. \"Large language diffusion models.\" arXiv preprint arXiv:2502.09992 (2025)."}]

Author rebuttal:

Meta-review: The work proposes a higher-order numerical inference scheme for discrete diffusion models. The paper provides theoretical justifications for their proposed scheme and provides empirical tests of their method. The main concern from two of the reviewers is the limited scope of the experimental evaluations. Specifically, Reviewer 61BM wrote that the work would be improved by testing how the higher-order solver interacts with the following:
> - Effect of Temperature/Confidence
> - Multi-token unmasking heuristics like in FastDLLMs, etc.
> - ReMasking
> - Predictor-Corrector Solvers
> - How well higher order solvers work with distilled discrete diffusion models based on e.g., SDTT Checkpoints

In practice, it is important to test the effects of these interactions. However, I believe the theoretical contributions are strong enough to warrant acceptance without additional empirical validation.
Decision: Accept (poster)