Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
Authors: Huanjian Zhou
Masashi Sugiyama
Venue: ICML 2025 spotlightposter (2025)
Keywords: parallel sampling
log-concave sampling
diffusion model
Primary area: probabilistic_methods->monte_carlo_and_sampling_methods
TL;DR: We propose new parallel algorithms for log-concave sampling and diffusion models, and rigorously prove that our algorithms enjoy  $\log d$ iteration complexity.

Abstract: Sampling from high-dimensional probability distributions is fundamental in machine learning and statistics. As datasets grow larger, computational efficiency becomes increasingly important, particularly in reducing *adaptive complexity*, namely the number of sequential rounds required for sampling algorithms. While recent works have introduced several parallelizable techniques, they often exhibit suboptimal convergence rates and remain significantly weaker than the latest lower bounds for log-concave sampling.
To address this, we propose a novel parallel sampling method that improves adaptive complexity dependence on dimension $d$ reducing it from $\widetilde{\mathcal{O}}(\log^2 d)$ to $\widetilde{\mathcal{O}}(\log d)$.  Our approach builds on parallel simulation techniques from scientific computing.

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
Huanjian Zhou 12 Masashi Sugiyama" 2
Abstract
For Problem (a), specifically log-concave sampling, start- ing from the seminal papers of Dalalyan & Tsybakov (2012), Dalalyan (2017), and Durmus & Moulines (2017), there has been a flurry of recent works on proving non- asymptotic guarantees based on simulating a process which converges to т over time (Wibisono, 2018; Vempala & Wibisono, 2019; Mou al., 2021; Altschuler & Talwar, 2023). Moreover, these processes, such as Langevin dy- namics, converge exponentially quickly to under mild conditions (Dalalyan, 2017; Mou et al., 2021; Bernard et al., 2022). Such dynamics- based algorithms for Problem (a) share a common feature with the inference process of SGMs that they are actually a numerical simulation of an initial- value problem of differential equations (Hodgkinson et al., 2021). Thanks to the exponentially fast convergence of the process, significant efforts have been conducted on discretiz- ing these processes using numerical methods such as the forward Euler, backward Euler (proximal method), exponen- tial integrator, mid-point, and high-order Runge- Kutta meth- ods (Vempala & Wibisono, 2019; Wibisono, 2019; Shen Lee, 2019; Li et al., 2019; Oliva & Akyildiz, 2024).
Sampling from digensimensiona probability dis- tributions is fundamental in machine learning and statisties. As datasets grow larger, computational efficiency becomes increasingly important, par- ticularly in reducing adaptive complexity, namely the number of sequential rounds required for sam- pling algorithms. While recent works have intro- duced several parallelizable techniques, they often exhibit suboptimal convergence rates and remain significantly weaker than the latest lower bounds for log-concave sampling. To address this, we propose novel parallel sampling method that im- proves adaptive complexity dependence on dimen- sion reducing it from Öllog d) to Ö(log d). Our approach builds on parallel simulation tech- niques from scientific computing.
1. Introduction
We study the problem of sampling from probability distri- bution with density x(x) x exp(-f(x)) where f Rd + R is a smooth potential. We consider two types of set- ting. Problem (a): the distribution known only up to a normalizing constant (Chewi, 2023), and this kind prob- lem is fundamental in many fields such as Bayesian infer- ence, randomized algorithms, and machine learning (Robert et al., 1999; Marin et al., 2007; Nakajima et al., 2019). Problem (b): known as the score-based generative models (SGMs) (Song& & Ermon, 2019), we are given an approx- imation of V log Tr, where n, is the density of a specific process at time I. The law of this process converges to π over time. SGMs are state-of-the-an in applications like image generation (Ho et al., 2022a; Dhariwal & Nichol, 2021), audio and video generation (Ho et al., 2022b; Yang et al., 2023a), and inverse problems (Song et al., 2022).
Furthermore, in recent years, there have been increasing in- terest and significant advances in understanding the conver- gence of inherently dynamics-based SGMs (Bortoli, 2022;
Chen et al., 2023c; Lee et al., 2023; Pedrotti et al., 2024;
Chen et al., 2023b; Tang & Zhao, 2024; Li & Yan, 2024). Notably, polynomial-tne been established (Chen et al., 2023c;b; Benton et al., 2024; Liang al., 2025), and various discretization schemes for SGMs have been analyzed (Lu et al., 2022a;b; Huang et al., 2025).
convergence guarantees have
The algorithms underlying the above results are highly se- quential. However, with the increasing size of data sets for sampling, we need to develop a theory for algorithms with limited iterations. For example, the widely-used denoising diffusion probabilistic models (Ho et al., 2020) may take 1000 denoising steps to generate one sample, while the eval- uations of a neural network-based score function can be computationally expensive (Song et al., 2021).
'The University of Tokyo, Tokyo, Japan "Center for Advanced Intelligence Project, RIKEN, Tokyo, Japan. Correspondence to: Huanjian Zhou <zhou ms www.k..uttokk.... tokyo..c.jp....
As a comparison, recently, the (naturally parallelizable) Pi- card methods for diffusion models reduced the number of steps to around 50 (Shih al., 2024). Furthermore, in terms
Proceedings of the 4nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
Table 1. Comparison with existing parallel methods for strongly
of the dependency on the dimension d and accuracy e, Picard methods for both Problems (a) and (b) were proven to be able to return an e-accurate solution within (log" (d/=²)) iterations, improved from previous (d /=") with some a, b > 0. However, the Ologg" (d/=²)) adaptive complex- ity' may not be yet optimal for both Problem (a) and (b).
log-coneave sampling. The symbol represents that the results hold under a weaker condition, the log-Sobolev incquality.
Work
Adapdive
Space
dynamics
Measure
Complexity
Complexity
(Sben Lee, 2019. Theorem 4)
W,
5(lg" (:9)
o(*")
underdamped Langevin diffusion
(Yu & Dalla Jalaayyn, 2024, Covolary 2) underdamped Langovin diffusion
W2
Ó(lor" (#))
((
This motivates our investigation into the question:
(Anari et al., 2024. Theorem 15)
TV
O(lo""
o(*"))
underdamped Langevin diffusion
Can we achieve logarithmic adaptive complexity for both log-concave sampling and sampling for SGMs?
(Anari et il., 2024, Theorem 13)
overdamped Langevin diffiusion
KL
6(kgg (A))
Ō($)
Theorem 4.2 overdamped Langevin diffusion
KL
O log (4))
o(E)
Our Contributions
Table 2. Comparison with existing parallel methods for sampling for diffusion models.
In this work, we propose a novel sampling method that employs a highly parallel discretization approach for con- tinuous processes, with applications to the overdamped Langevin diffusion (Chewi, 2023) and the stochastic dif- ferential equation (SDE) implementation of processes in SGMs (Chen et al., 2024) for Problems (a) and (b), respec- tively.
Works Implementation
Measure
Adaptive Complesity
Complexity
Space
(Chen al., 2024, Theorem 3.5)
TV
"(lgg" ( ))
d 2/2
ODE
mechod
(Grppa et al., 2025, Theorem 8.13)
ODE/ Parallel midpoint methed (Chen al., 2024, Theorem 3.3)
TV
Ō ( ))
ō 2
KL
Ö Ing" (+4)
8
SDE Picard method Theorem 5.4
SDE Paralel Picand method
KL
6 log (,5 ))
Faster parallel log-concave sampling. We first present an improved result for parallel sampling from strongly log-concave and log-smooth distribution. Specifically, we improve the upper bound from Ô (log" (=)) (Anari et al., 2024) to Ô (log (1)), with slightly scaling the number of processors and gradient evaluations from O d to
Gupta et al., 2025), with slightly increasing the number of the processors and gradient evaluations and the space com- plexity for SDEs. We summarize the comparison in Table 2. Similarly, the better space complexity of the ordinary differential equation (ODE) implementations attributed to the smoother trajectories of ODEs, which are more readily discretized.
O 2 log (=)).
Compared with methods based on underdamped Langevin diffusion (Shen & Lee, 2019; Yu & Dalalyan, 2024; Anari et al., 2024), our method exhibits higher space complexity². This is primarily because underdamped Langevin diffusion typically follows smoother trajectory than overdamped Langevin diffusion, allowing for larger grid spacing and consequently, a reduced number of grids. We summarize the comparison in Table 1. on the adaptive complexity and discretization schemes for overdamped Langevin diffusion.
2. Problem Set-up
In this section, we introduce some preliminaries and key ingredients of log-concave sampling and diffusion models in Sections 2.1 and 2.2, respectively. Following this, Section 2.3 provides an introduction to the fundamentals of Picard iterations.
In this paper, we will focus
2.1. Log-concave Sampling
Problem (a) (Sampling task). Given the potential function
f D + R. the goal of the sampling task is to draw a sample from the density T = Z exp(-f), where Zy == Jo exp(-f(x))ddr is the normalizing constant.
Faster parallel sampling for diffusion models. We then present an improved result for diffusion models. Specifi- cally, we propose an efficient algorithm with θ (log (#) adaptive complexity for SDE implementations of diffusion models (Song& & Ermon, 2019). Our method surpasses all the existing parallel methods for diffusion models hav- ing O (log" (=) adaptive complexity (Chen et al., 2024;
Distribution and function class. If f is (strongly) convex, the density Ty is said to be (strongly) log-concau If f is twice-differentiable and V < 3I (where x denotes the Loewner order and is the identity matrix), we say the potential f is B-smooth and the density n is --ooohh
'Adaptive complexity refers to the minimal number of sequen- tial rounds required for an algorithm to achieve desired accuracy, assuming polynomially many queries can be executed in parallel each round (Balkanski & Singer, 2018).
We define relative Fisher information of probability density w.r.t. as FI(p][T) = I log(p/>))* and the Kull- back-Leibler (KL) divergence of p from ι as (p = Ep log(p/)). If h is a-strongly log- lonccveevve, then the fol- lowing relation between KI. divergence and Fisher informa-
"We note, in this paper, that the space complexity refers to the number of words (Cohen-Addad et al., 2023; Chen et al., 2024) instead of the numher of bits (Goldreich, 2008) to denote the approximate required storage.
2

---

Paraliel Simulation for Log-coneave Sampling and Score-based Diffusion Models
tion holds:
2.2. Score-based Diffusion Models
1 =I( for all probability measures p. (L(p]()) - 2a
Sampling for diffusion models. In score-based diffusion models, one considers forward process (2)))e in R governed by the canonical Ceastein-hhlebeeck (OU) pro- cess (Ledoux, 2000):
Langevin Dynamics. One of the most commonly-used dynamics for sampling is Langevin dynamics (Chewi, 2023),
} al + dB,, dæ, - 2
which is the solution to the following SDE,
Co ~ 4o>
I E [0, T), (I)
dæ dx))dt + V2dB,,
where Oo is the initial distribution over Rd The correspond-
ing backward process (Ze)re(o.r) in I& follows an SDE defined as
where (B)) is standard Brownian motion in IRd. If x exp(-f) satisfies strongly log-concavity, then the law of the Langevin diffusion converges exponentially fast to T (Bakry et al., 2014).
dØ 1 + V Vlogū.(ā)) dt + dB, E [0, To ~ Pa = (o
Score function for sampling task. We assume the score function $ Rd + R is a pointwise accurate estimate of VS. i.e., [|>())) - (=)] < 8 for all æ E Rd and some sufficiently small constant 8 R.+.
(2)
where ((:.)) represents the normal distribution over R In practice, the score function Voo P, (wx)) estimated by neural network (NN) sf R Rd where 0 is the parameters of NN. The backward process is approximated by
Measures of the output. For two densities and n, we define the total variation (TV) as
dy, By, + s((1) di d ddB, LE [0, yo ~ V()..I.).
TV(e, π) = sup(p(E) T(E) E is an event].
We have the following relation between the KL divergence and TV distance, known as the Pinsker inequality,
Problem (b) (Sampling task for SGMs). Given the learned NN-based score function 8f the goal is to simulate the approximated backward process such that the law of the output is close qo-
TV(p, S (L(pllT).
Distribution class. For SGMs, we assume the data density Po has finite second moments and is normalized such that covp, (Xn) = Ep [(x.- [ [=oo) - (z0)) = Ia. Such a finite moment assumption is standard across previous theoretical works on SGMs (Chen et al., 2023a;b;e) and we adopt the normalization to simplify true score function-related computations as Benton et al. (2024) and Chen et al. (2024) did.
We denote by W2 the Wasserstein distance between p and n. which is defined as
Wh x) inf (( [IX- -YI]]
II
where the infimum is over coupling distributions II of (X, Y) such that X ~ p, Y ~ n. If ग is a-strongly log- concave, the following transport-entropy inequality, known as Talagrand's To inequality, holds (Otto & Villani, 2000) for all p Pa )IR
OU process and inverse process The OU process and its inverse process also converge to the target distribution ex- ponentially fast in various divergences and metrics such as the 2-Wasserstein metric W2; see Ledoux (2000). Further- more, the discrepancy between the terminal distributions of the backward process (2) and its approximation version (3) scales polynomially with respect to the length of the time horizon and the score matching error. (Huang et al. (2025, Theorem 3.5) or setting the step size h 0 for the results in Chen et al. (2023a;b;c)).
,, i.e., with finite second moment,
W2 R) S (L(P)TT).
2
Complexity. For any sampling algorithm, we consider the adaptive complexity defined as unparallelizable evaluations of the score function (Chen et al., 2024), and use the notion of the space complexity denote the approximate required storage during the inference. We note, in this paper, that the space complexity refers to the number of words (Cohen- Addad et al., 2023; Chen et al., 2024) instead of the number of bits (Goldreich, 2008) to denote the approximate required storage.
Score function for SGMs. For the NN-based score, we assume the score function is accuraaee bounded and Lip- schitz; we defer the details in Section 5.2.
3

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
2.3. Picard Iterations
et al., 2025). With Olog d) Picard iterations for Öllog d) time slices, the total adaptive complexity of their algorithms is (log" d). However, while sequential updating of each time slice is not necessary for simulating the process, it remains unclear how to parallelize across time slices for sampling to obtain (log d) time complexity.
Consider the integral form of the initial value problem, at = To + J So (ag)dd + v2B, The main idea of Picard iterations (Clenshaw, 1957) is to approximate the difference over time slice [kn bn+1] as
Tla+1 - -bhp
Algorithmic novelty: parallel methods across time slices. Naively, if we directly update all the grids simultaneously, the Picard iterations will not converge when the total length is T = Ō (log d). Instead of updating all time slices together or updating the time slice sequentially, we update the time slices in a diagonal style as illustrated in Figure 1. For any j-th update at the n-th time slice (corresponding the rectangle in the n-th column from the left and the j-th row from the top in Figure 1), there will be two inputs: (a) the right boundary point of the previous time slice, which has been updated j times, and (b) the points on the girds of the same time slice that have been updated j j-1 times, Then we perform times Picard iterations with these inputs, where the hyperparameter P depends on the smoothness of the score function. This repetition is simple but essential for preventing the accumulation of matee-na errors, which could otherwise grow exponentially w.r.t. the length of the time horizon. The main difference compared to the existing Picard methods is that for a fixed time slice, the starting points in our method are updated gradually, whereas in existing methods, the starting points remain fixed once processed.
Mu+1 fa fa( + (2 - B,n))
M
Wist tTa, (1) )ds +
V2(Bt... Bt.):
with a discrete grids of M collocation points Lt, toat I tT., S +T., S ..... +T. -t,-1- We update the points in a wave-like fashion, which inherently allows for parallelization: for m' = =I...., M,
m
p+1 to +Tomm
= Z0
+
Σ
Wmm Sent Tn (2f mm
m
222BBBB ( +Tamm B...)
Various collocation points have been proposed, including uniform points and Chebyshev points (Bai & Junkins, 2011). In this paper, however, we focus exclusively the simplest case of uniform points, and extension to other cases is future work. Picard iterations are known to converge exponentially fast and, under certain conditions, even factorially fast for ODEs and backward SDEs (Hutzenthaler et al., 2021).
Time Dequthr N Öikgaa time slios
initial dist, mo
Lungst diel.
3. Technical Overview
N
We adopt the time splitting for the time horizon used in the existing parallel methods (Shen & Lee, 2019; Gupta et al., 2025; Chen et al., 2024; Anari et al., 2024; Yu & Dalalyan, 2024). With same time grids, our algorithm, however, depart crucially from prior work in the design of parallelism across the time slices, and the modification for controlling the score estimation error. Below we summarize these algorithmic contributions and technical novelties.
2
N
Pleaed Depth: J - (loag) d) depei
N+2
J+1
N+J
Figure 1. Illustration of the parallel Picard method: each rectan- gle represents an update, and the number within each rectangle indicates the index of the Picard iteration. The approximate time complexity is N - J = Ö(lo d).
Recap of existing parallel sampling methods. Existing works for parallel sampling apply the following generic dis- cretization schemes (Shen & Lee, 2019; Gupta et al., 2025;
Challenges for convergence. Similar to the arguments for sequential methods or parallel methods with sequentially up- dating the time slices, we use the standard techniques such as the interpolation method or Girsanov's 's theorem (Vem- pala & Wibisono, 2019; Oksendal, 2013; Chewi, 2023) to decompose the total error w.r.t. KL into four components: (i) convergence error of the continuous process, (ii) dis- cretization error, (iii) parallelization error, and (iv) score estimation error (See Eq. (4), Lemma B. and Lemma C.4). For (i) the convergence error of the continuous processes,
Chen et al., 2024; Anari et al., 2024; Yu & Dalalyan, 2024). At a high level, these methods divide the time horizon into many large time slices and each slice is further subdivided into grids with small enough step size. Instead of se- quentially updating the grid points, they update all grids at the same time slice simultancously using exponentially fast converging Picard iterations (Alexander, 1990; Chen et al., 2024; Anari et al., 2024), or randomized midpoint methods (Shen& & Lee, 2019; Yu & Dalalyan, 2024; Gupta
4

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
their exponential convergence rates allow this error to be effectively controlled by setting the total time length to (log $ ). regardless of the specific discretization scheme. For (ii), the discretization етrоr scales approximately as where h denotes the time slice length, M the number M.
To address the initial challenge related the truncation error, we choose the Picard depth as J = O(N log d). We first bound the error of the output for each update with respect to its inputs as L S |Lا 1 + b where a and are constants. By carefully choosing the length of the time slices, we can ensure that < along the Picard iteration direction. Consequently, the truncation error will converge if the iteration depth J is sufficiently large, such that a b' is sufficiently small. This requirement implies that J = O(N + log d).
d. A
of discretization points per slice, and h NI the grid resolu- tion (See Eq. (4), and Lemma C.5). Setting .7 x /d2 ensures the discretization error remains within O(E²). The technical challenges rise from controlling the remaining two errors, which we summarize below.
(iii) Parallelization error:
the parallelization error primar- ily arises from updating with 1(x)-1) instead of s(x)) in the Picard iteration, in contrast to the sequential method, where j indexes the steps along the Picard direction. In existing parallel methods, the sequential update across time slices benefits the convergence of truncation errors, E I|| -ل-1 |]. Assuming the truncation errors in the previous time slice have converged, its right boundary serves as the starting point for all grids in the current O(L)-enggh time slice which results in an initial bias of (d). Subse- quently, by performing (log d) exponentially fast Picard iterations, the truncation error will converge. However, in our diagonal-style updating scheme across time, the trunca- tion error interacts with inputs from both the previous time slice and prior updates in the same time slice. Consequently, the convergnee loop that holds in sequential updating no longer holds.
To mitigate the additional score estimation error for Problem (a), we perform P Picard iterations within each update. The interaction between the truncation error and additional score estimation error can be expressed as T S aT + P + c6², where a, b, c are constants. To ensure the total score estimation error remains bounded, it is necessary to have a, b < 1, which guarantees convergence along both the time and Picard directions. By the convergence of the Picard iteration, we can achieve b < 1. For a, the right boundary point of the previous time slice, and prior updates within the same time slice introduce discrepancies in the truncation error. For the impact from the previous time slice, we make use of the contraction of gradient decent to ensure convergence. However, since the grid gap scale as 1/d, the contraction factor is close to 1. Consequently, we have to minimize the impact from prior updates within the same time slice, which scales as ((1) by repeating P = log O11) Picard iterations for each update.
(iv) Score estimation error: If the score function itself is Lipschitz continuous (Assumption 5.3 for Problem (b)), no additional score matching error will arise during the Picard iterations. This allows the total score estimation error to remain bounded under mild conditions (Assumption 5.1). However, for Problem (a), since it is the velocity field ff instead of the score function s that is Lipschitz, additional score estimation errors will occur during each update. For the sequential algorithm, these additional score estimation errors are contained within the bias-convergence loop, ensur- ing the total score estimation error remains to be bounded. Conversely, for our diagonal-style updating algorithm, the absence of convergence along the time direction causes these additional score estimation errors to accumulate expo- nentially over the time direction.
Balance between time and Picard directions. We note that the Picard method, despite being the simplest approach for time parallelism, has achieved the state of art perfor- mance in certain specific settings. On the one hand, the continuous processes need to run for least Olog d) time. To ensure convergence within every time slice, the time slice length have to be set (1), resulting in necessity for at least O(log d) iterations. On the other hand, with a proper initialization (d), Picard iterations converge within (log d) iterations. Our parallelization balances the conver- gence of the continuous diffusion and the Picard iterations to achieve the improved results.
Related works in scientific computation. Similar paral- lelism across time slices has also been proposed in scientific computation (Gear, 1991; Gander, 2015; Ong & Schroder, 2020). especially for parallel Picard iterations (Wang, 2023). Compared with prior work in scientific computation, our approach exhibits several significant differences. Firstly, our primary objective differs from that in simulation. In sampling, we aim to ensure that the output distribution closely approximates the target distribution, whereas simu- lation seeks to make some points on the discrete grid closely match the true dynamics. Second, our algorithm differs sig-
Technical novelty. Our technical contributions address these challenges by the appropriate selection of the number of Picard iterations within each update P and the depth of the Picard iterations J. We outline the details of the choices below.
In the following, we assume that the truncation error at the n-th time slice and the j-th iteration scales with L3 and that the additional score estimation error for each update scales with 82.
5

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
nificantly from that of Wang (2023). In our algorithm, each update takes the inputs without the corrector operation. Fur- thermore, we perform P Picard iterations in each update to prevent error accumulation over time T = Ólog d). How- ever, these two fields are connected through the sampling strategies that ensure each discrete point closely approxi- mates the true process at every sampling step.
Algorithm 1 Parallel Picard Method for sampling
1: Input: 3y ~ μo. approximate score function 8 R Vf. the number of the iterations in outer loop J, the number of the iteration in inner loop P, the number of time slices N, the length of time slices h, the number of points on each time slices M.
2: for n = .....,N - 1 do M (in parallel) do
3: for m =
4. Parallel Picard Method for Strongly Log-concave Sampling
4:
Buh+m/Msk = Buh + +O, (mh/ M))a),
5:
x² -,- 2o, for .,...,J,
6:
2b.m = x 1
hm
M (x -1,M) +
In this section, we present parallel Picard methods for strongly log-concave sampling (Algorithm 1) and show it holds improved convergence rate w.r.t. the KL divergence and total variance (Theorem 4.2 and Corollary 4.3). We illustrate the algorithm in Section 4.1, and give a proof sketch in Section 4.3. All the missing proof's can be found in Appendix B.
v (2 Butt,
7: end for
8: end for
9: for k = 1,..,, N do
10: for j = 1...., min fk: - 1, and m = 1........ M
(in parallel) do
11:
let n j. x
x²/
and x
xx
n,m
12:
for p 1..., Pdo
4.1. Algorithm
13:
x
h
Σ
s
(( +
Our parallel Picard method for strongly log-concave sam- pling is summarized in Algorithm 1. In Lines 2-7, we gener- ate the noises and initialize the value at the grid via Langevin Monte Carlo (Chewi, 2023) with a stepsize h (1). In Lines 8-26, the time slices updated in diagonal manner within the outer loop, illustrated in Figure 1. In Lines 12-14 and Lines 21-23, we repeat P Picard iterations for each update to ensure convergence.
MI
m/=0
v V2(Bnh+mh/Mr - Bur):
14:
end for
15: xy x
16: end for
17: end for
18: for k = N 1,.... N + J do
19: for n = max (0, k - .J....,, N and m
Remark 4. Parallelization should be understood as evalu- ating the score function concurrently, with each time slice potentially being computed in an asynchronous parallel man- ner, resulting in the overall P(N + J) + N adaptive com- plexity.
1,..., M (in parallel) do
20:
let k-n, n, 20 x
1. and
z0.0 n,m
21:
for p p= 1 1...., Po
m-1
22:
x 0
h
Σ
s(
s(acig)))))) +
M
m/=0
/2(BnA+mM/M Bar).
4.2. Theoretical Guarantees
23: end for
24: xy.n x/I.
The following theorem summarizes our theoretical analysis for Algorithm
25: end for
26:
end for
Theorem 4.2. Suppose π is a-strongly log-concave and 3-log-smooth, and the score function 8 is 6-accurate. Let k= B/a. Suppose
27:
Return: N-1.
To make the guarantee more explicit, we can combine it with the following well-known initialization bound, see, e.g., Dwivedi et al. (2019, Section 3.2).
Rd gZ
(LL e2
3h- 0.1,
M Z
N 2 10x log
P2 2log n
Corollary 4.3. Suppose that π =exp(-> (pp(-1 -f) is a-strongly log-concave and B-log-smooth, and let 3/a. Let x² be the minimizer of f. Then, for po Na'" , it holds that KL [ < 2 log K. Consequently, setting
5 0.2V
+ 4,
and J-N2 log (N3/ h88h + (Hol(H)) + 122d
2²
then Algorithm I runs within N + (N J)P iterations with at most MN queries per iteration and outputs a sample with marginal distribution p such that
1
M kd
d log
h
= log
1033
Va (W m), TV(o.))
KL( 2
2 2og K
max
< 25.
8 8 0.2
P>
+ 4,
3
6

---

Parallel Simulation for Log-coneave Sampling and Score-based Diffusion Models
and J- - O (log lod K
ensures sufficiently long time horizon the sufficiently long time horizon T = Nh. Olog _ densely spaced grids with gap h/M = ( /d) and small score matching error, respectively, we can conclude that
E²
then Algorithm I runs within N + (N
On d
iterations with most MN O e log 1 queries per iteration and outputs sample with marginal distribution such that
e KL((oll) + dh
M + 82 Se
va W2(e, R), TV( sV KL(P,T) x 2
Thus, we will focus on proving the convergence of the truncation error & in the Picard iterations, and avoiding the additional accumulation of the score estimation error during Picard iterations as discussed before.
max
< 2c.
Remark 4.4. The main drawback of our method the sub- optimal space complexity due its application over- damped Langevin diffusion which has a less smooth trajec- tory compared to underdamped Langevin diffusion. How- ever, we anticipate that our method could achieve com- parable space complexity when adapted to underdamped Langevin diffusion.
Considering that the truncation error expands most ex- ponentially along the time direction, but diminishes expo- nentially with an increased depth of the Picard iterations, convergence can be achieved by ensuring that the depth of the Picard iterations surpasses the number of time slices as J 2 + Ollog = with initialization error bounded by (d) (the second part of Corollary B.7 and second part of Corollary B.9).
Remark 4.5. Regarding the condition number K, our method achieves the same adaptive complexity of ((k) as both the state-of-the-art sequential method and existing paral- lel approaches (Anari et al., 2024; Yu & Dalalyan, 2024; Altschuler Chewi, 2024). Whether parallelization can improve the dependence on the condition number remains an open question, which leave for future work.
Due to the non- non-Lipschitzzess, of the score function, we can only bound & by quantity - (Lemma B.5 and Lemma B.8), where A represents the trunca- tion error the right boundary of the previous time slice. Here, the coefficients are given by: = 1 0.1 Sh + (xx)3352)) b = (((32))))) and c = (kk5h?2) Intuitively, a comes from the contraction of the gradient mapping with an additional term from the Picard direction, b reflects convergence along the Picard direction, and c accounts for the accumulation of score estimation error å over time length h, with an additional scaling by due to Young's inequality. To control the growth of the score error, it is essential that the coefficients a and b remain strictly less than one. Setting P = Olog x) is sufficient to ensure this condition.
Remark 4.6. When the number of computation cores, de- noted by f, limited, the adaptive complexity of our al- gorithm is Ö( ( log² _ This matches the state-of-the- art adaptive complexity of Anari et al. (2024) and Yu & Dalalyan (2024) when f $ xd However, our algorithm 2 with potential applications as demonstrated in Nishihara et al. (2014), De Souza et al. (2022), Hafych et al. (2022) and Glatt-Holtz et al. (2024).
achieves improved adaptive complexity when é 2 nd
4.3. Proof Sketch of Theorem 4.2: Performance Analysis of Algorithm 1
5. Parallel Picard Method for Sampling of Diffusion Models
The detailed proof of Theorem 4.2 is deferred to Appendix B. As discussed in Section 3, by interpolation methods (Anari et al., 2024), we decompose the error w.r.t. the KL diver- gence into four error components (corollary B.4):
In this section, we present parallel Picard methods for diffu- sion models Section 5.1 and assumptions Section 5.2. Then we show it holds improved convergence rate w.r.t. the KL divergence (Theorem 5.4). All the missing details can be found in Appendix C.
(L((
dh M
-ONN
Convergence of Langevin dynamics
discretization error
5.1. Algorithm
N c ((x)
Ex-n +
82 Score estimation error
(4)
Due space limitations, the detailed methodology for the parallelization of Picard methods for diffusion models is provided in Appendix C.I and Algorithm 2. It keeps same parallel structure as that illustrated in Figure 1. Notably, it exhibits the following distinctions in comparison to the parallel Picard methods for strongly log concave sampling presented in Algorithm 1:
paralizaaton
where ε represents the truncation error of the grids at n-th time slice after j update. For the right terms, with the choice of N- = O(lo 2), M = Odhh/e2) and 5 S E, which
7

---

Assumption 5.1 ((L2((), tx]) 52-accurate learned score). The learned NN-based score of is da-accurate in the sense of Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
Since the score function itself is Lipschitz, there will not be additional score matching error during Picard iterations. As result, we perform single Picard iteration in one update, i.e., P=11
with total 2N + Ō 2) adaptive complexity and dM = Ō space complexity for parallelizable 62- accurate score function computations.
Instead of uniform discrete grids, we employ a shrinking step size discretization scheme towards the data end, and the early stopping technique which is unvoidable show the convergence for diffusion models (Chen et al., 2024).
Remark 5.5. Compared to existing parallel methods, our method improves the adaptive complexity from (log" ()) to O _ ). Its main drawback is suboptimal space com- plexity due to the less smooth trajectory in SDE implemen- tations, but we believe it can achieve comparable space complexity when adapted to ODE implementations.
We show the details in Appendix C.I;
We use an exponential integrator instead of the Euler- Maruyama Integrator in Picard iterations, where an additional high-order discretization error term would emerge (Chen al., 2023a), which we believe would not affect the overall (log d) adaptive complexity with parallel sampling.
Remark 5.6. When the number of computation cores & is lim- ited, the adaptive complexity is Ō dd log" _ matching the state-of-the-an result of Chen et al. (2024) for e E 2 However, employing a large batch size 2 nd in diffusion models may not always yield substantial benefits (Shih et al., 2024; Li et al., 2024b;a).
e²
5.2. Assumptions
Our theoretical analysis of the algorithm assumes mild con- ditions regarding the data distribution's 's regularity and the approximation properties of NNs discussed in Chen al. (2024). These assumptions align with those established in previous theoretical works, such as those described by Chen et al. (2023c;a;;; 2024).
Remark 5.7. We note that the uniformly Lipschitz assump- tion (Assumption 5.3) may be too strong. particular, the required Lipschitz constant can become quite large even becoming unbounded near the zero point (Salmona et al., 2022; Yang et al., 2023b). In this case, to ensure conver- gence of the Picard iterations under these conditions, the quantity Lie h2 must be sufficiently small. This require- ment implies that the length of each time slice, hn, should scale as (1 (()) Consequently, the number of time slices becomes N O (I2 log d), leading to an overall iteration complexity of N (L, log d). We also believe our algo- rithm, which is based on diagonal-style update, is robust to this assumption by adaptively adjusting the length of the time slices.
N- Ma
Eg
0 fn
0=
m==0
V log in + Fa.m (ā +n m|1?
< 82.
Assumption 5.2 (Regular and normalized data distribu- tion). The data density Po has finite second moments and normalized such that po (20) I...
6. Discussion and Conclusion
In this work, we proposed novel parallel Picard methods for various sampling tasks. Notably, we obtain aa-accurate sample w.r.t. the KL divergence within Ö (log 2 ). which represents a significant improvement from Ö (log² 4 for diffusion models. Furthermore compared with the existing methods applied to the overdamped Langevin dynamics or the SDE implementations for diffusion models, our space complexity only scales by a logarithmic factor.
Assumption 5.3 (Bounded and Lipschitz learned NN-based score). The learned NN-based score function sy has a bounded C² norm, i.e. with Lipschitz constant L..
. |||| ()Ill S M
5.3. Theoretical Guarantees
Theorem 5.4. Under Assumptions 5. 5.2, and 5.3, given the following choices of the order of the parameters
Our study opens several promising theoretical directions. First, as an analogue to simulation methods in scientific com- puting, it highlights the potential of leveraging alternative discretization techniques for faster and more efficient sam- pling. Another direction is exploring smoother dynamics to reduce space complexity in these methods.
h h 01), N N =
0 (og d )
M M=0(4 log
E2
g2
Nd 8 SSE, and J O +log -²
T T= T= log 5z
the parallel Picard algorithm for diffusion models (Algo- rithm 2) generates samples from satisfies the following error bound,
Lastly, although our highly parallel methods may introduce engineering challenges, such as the memory bandwidth, we believe our theoretical works will motivates the empirical development of parallel algorithms for both sampling and diffusion models.
dT
(L, S de T
dT + M
+2 + 82 S 2
(5)
8

---

Paralicl Simulation for Log-coneave Sampling and Score-based Diffusion Models
Acknowledgment
Chen, H., Lee, H., and Lu, J. Improved Analysis of Score- based Generative Modeling: User-Friendly Bounds un- der Minimal Smoothness Assumptions. In International Conference on Machine Learning. ICML 2023, 23-29 July 2023, Honolulu, Hawaii, USA, volume 202 of Pro- ceedings of Machine Learning Research, pp. 4735-4763. PMLR, 2023a.
The authors thank Sinho Chewi for very helpful conversa- tions. HZ was supported by International Graduate Program of Innovation for Intelligent World and Next Generation Artificial Intelligence Research Center. MS was supported by JST ASPIRE Grant Number JPMJAP2405.
Impact Statement
Chen, H., Ren, Y., Ying, L.., and Rotskoff, G. M. Accelerat- ing Diffusion Models with Parallel Sampling: Inference at Sub-Linear Time Complexity. In The Thirty- Annual Conference on Neural Information Processing Systems, 2024.
This work focuses on the theory of accelerating the sampling via parallelism. As far we can see, there is no foresecable negative impact on the society.
References
Chen, S., Chewi, S., Lee, H., Li, Y., Lu, J., and Salim, A. The probability flow ODE is provably fast. In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S. (eds.). Advances Neural Information Processing Systems 36: Annual Conference on Neural Information Processing Systems 2023, NewrlPS 2023, New Orleans, LA, USA, December 10 16, 2023, 2023b.
Alexander, R. Solving Ordinary Differential Equations Nonstiff Problems. Siam Review, 1990.
Altschuler, J. and Talwar, K. Resolving the Mixing Time of the Langevin Algorithm its Stationary Distribution for Log-Coo LLog-Concave Sampling. In Proceedings of Thirty Sixth Conference on Leurning Theory, 2023.
Chen, S., Chewi, S., Li, J., Li, Y., Salim, A., and Zhang, A. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In The Eleventh International Conference Learning Repre- sentations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023, 2023c.
Altschuler, J. M. and Chewi, S. Shifted Composition III: Lo- cal Error Framework for KL Divergence. arXiv preprint arXiv:24/2.17997 2024.
Anari, N., Chewi, S., and Vuong, T.-D. Fast parallel sam- pling under isoperimetry. Proceedings Thirty Seventh Conference on Learning Theory, pp. 161-185, 2024.
Chewi, S. Log-concave sampling. Book draft available at https://chewisinho. github. io, 2023.
Bai, X. and Junkins, J. L. Modified Chebyshev-Picee PPicard iteration methods for orbit propagation. The Journal of the Astronautical Sciences, 2011.
Chewi, S., Erdogdu, M. A., Li, M., Shen, R., and Zhang, M. S. Analysis of Langevin Monte Carlo from Poincaré Log-Sobolev. Foundations of Computational Mathe- matics, 2024.
Bakry, D., Gentil, I., Ledoux, M., et al. Analysis and geome- try of Markov diffusion operators, volume 103. Springer, 2014.
Clenshaw, C. The numerical solution of linear diffcrential equations in Chebyshev series. In Mathematical Proceed- ings of the Cambridge Philosophical Society. Cambridge University Press, 1957.
Balkanski, E. and Singer, Y. The adaptive complexity of maximizing a submodular function. In Proceedings of the 50th annual ACM SIGACT Symposium on Theory of Computing, 2018
Cohen-Addad, V., Woodruff, D. P., and Zhou, S. Streaming Fuclidean k-median and k-means with o(log n) Space. In 2023 IEEE 64th Annual Symposium on Foundations of Computer Science (FOCS). IEEE, 2023.
Benton, J., Bortoli, V. D., Doucet, A., and Deligiannidis, G. Nearly d-Linear Convergence Bounds for Diffusion Models via Stochastic Localization. The Twelfth Inter- national Conference on Learning Representations, ICLR 2024,, 2024.
Dalalyan, A. S. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Jour- nal of the Royal Statistical Society: Series B (Statistical Methodology), 2017.
Bernard, É., Fathi, M., Levitt, A., and Stoltz, G. Hypocoer- civity with schur complements. Annales Henri Lebesgue, 5:523-557, 2022.
Bortoli, V. D. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022. ISSN 2835-8856.
Dalalyan, A. S. and Tsybakov, A. B. Sparse regression learn- ing by aggregation and Langevin Monte-Carie Journal of Computer and System Sciences, 2012.

---

Paralicl Simulation for Log-concave Sampling and Seore-based Diffusion Models
De Souza, D. A., Mesquita, D., Kaski, S., and Acerbi, L. Parallel MCMC without embarrassing failures. In Interna- tional Conference Artificial Intelligence and Statistics, pp. 1786-1804. PMLR, 2022.
Hodgkinson, L., Salomone, R., and Roosta, Implicit Langevin algorithms for sampling from log-concave den- sities. Journal of Machine Learning Research, 22(136): 1-30, 2021.
Dhariwal, and Nichol, A. Diffusion models beat gans on image synthesis. Advances in neural information processing systems, 34:8780-8794, 2021.
Huang, D. Z., Huang, J., and Lin, Z. Convergence Analysis of Probability Flow ODE for Score- BereeBBsed Generative Models. IEEE Transactions on Information T'eeory 71 (6):4581-4601, 2025.
Dragomir, S. S. Some Gronwall type inequalities and appli- cations. Science Direct Working Paper, 2003.
Hutzenthaler, M., Kruse, T., and Nguyen, T. A. On the speed convergence of Picard iterations of backward stochastic differential equations. arXiv preprint arXiv:2107.018.m 2021.
Durmus, A. and Moulines, E. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. 2017.
Dwivedi, R., Chen, Y., Wainwright, M. J., and Yu, B. Log- concave sampling: Metropolis- Maassigo algorithms are fast. Journal of Machine Learning Research, 2019.
Ledoux, M. The geometry of markov diffusion genera- tors. In Annales de la Faculté des sciences de Toulouse: Mathématiques, 2000.
Gander, M. J. 50 years of time parallel time integration. In Multiple Shooting and Time Domain Decomposition Methods: MuS-TDD, Heidelberg, May 6-8, 2013, pp. 69-113. Springer, 2015.
Lee, H., Lu, J., and Tan, Y. Convergence of score-based generative modeling for general data distributions. In In- ternational Conference on Algorithmic Learning Theory. PMLR, 2023.
Li, G. and Yan, Y. Adapting to unknown low-dimensional structures in score-basee diffusion models. In The Thirty- eighth Annual Conference on Neural Information Pro- cessing Systems, 2024.
Gear, C. Waveform methods for space and time paral- lelism. Journal of Computational and Applied Math- ematics, 1991.
Glatt-Holtz, N. E., Holbrook, A. J., Krometis, J. A.. and Mondaini, C. F. Parallel MCMC algorithms: theoretical foundations, algorithm design, case studies. Transactions of Mathematics and Applications, t(2)00a, 2024.
Li, M., Cai, T., Cao. J., Zhang, Q., Cai, H., Bai, J., Jia, Y., Li, K., and Han, S. Distrifusion: Distributed parallel infer- ence for high-resolution diffusion models. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp- 7183-7193 2024a.
Goldreich, O. Computational complexity: a conceptual perspective. ACM Sigact News, 2008
Li, S.,, Yang, L., Jiang, X., Lu, H., Di, Z., Lu, W., Chen, J., Liu, K., Yu, Y., Lan, T., et al. SwiftDiffusion: Efficient Diffusion Model Serving with Add-on Modules. arXiv preprint arXiv:2407.0203/ 2024b.
Gupta, S., Cai, L., and Chen, S. Faster Diffusion-bari based Sampling with Randomized Midpoints: Sequential and Parallel. International Conference on Learning Represen- tations, 2025.
Li, X., Wu, Y., Mackey, L., and Erdogdu, M. A. Stochastic Runge- Kutta Accelerates Langevin Monte Carlo and Be- yond. Advances in neural information processing systems, 32, 2019.
Hafych, V., Eller, P., Schulz, O., and Caldwel, A. Paralleliz- ing MCMC sampling via space partitioning. Statistics and Computing, 32(4):56, 2022.
Liang, Y., Ju, P., Liang, Y., and Shroff, N. Broadening Target Distributions for Accelerated Diffusion Models via Novel Analysis Approach. The Thirteenth International Conference on Learning Representations, ICLR, 2025.
Ho, J., Jain, A., and Abbeel, P. Denoising diffusion proba- bilistic models. Advances in Neural Information Process- ing Systems, 2020.
Ho, J., Saharia, C., Chan, W., Fleet, D. J., Norouzi, M., and Salimans, T. Cascaded diffusion models for high fidelity image generation. Journal of Machine Learning Research, 23(47)11-33, 2022a.
Lu, C., Zhou, Y., Bao, F., Chen, J., and Li, C. A Fast ODE Solver for Diffusion Probabilistic Model Sampling in Around 10 Steps. Proc. Adv. Neural Inf. Process. Syst., New Orleans, United States, pp. 1-31, 2022a.
Ho, J., Salimans, T., Gritsenko, A., Chan, W., Norouzi, M., and Fleet, J. Video diffusion models. Advances in Neural Information Processing Systems, 35:8633-8644 3646, 2022b.
Lu, C., Zhou, Y., Bao, F., Chen, J., Li, C., and Zhu, J. Dpm- solver++: Fast solver for guided sampling of diffusion probabilistic models. arXiv preprint arXiv:2211.0109 2022b.
10

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
Marin, J.,M., Robert, C. P., al. Bayesian Core: Practi- cal Approach to Computational Bayesian statistics, vol- ume 268. Springer, 2007.
Song, Y. and Ermon, S. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
Mou, W., Ma, Y.A.., Wainwright, M. J., Bartlett, P. L.., and Jordan, M. I. High-order Langevin diffusion yields an accelerated MCMC algorithm. Journal of Machine Learning Research, 2021.
Song, Y., Shen, L.,, Xing, L., and Ermon, Solving inverse problems in medical imaging with score-based genera- tive models. In International Conference on Learning Representations, 2022.
Nakajima, S., Watanabe, K., and Sugiyama, M. Variational Bayesian learning theory. Cambridge University Press, 2019.
Tang, W. and Zhao, H. Contractive diffusion probabilistic models. arXiv preprint urXiv:2401133111 2024.
Vempala, S. S. and Wibisono, A. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. In Advances Neural Information Processing Systems 32, pp. 8092-8104, 2019.
Nishihara, R., Murray, I., and Adams, R. P. Parallel MCMC with generalized elliptical slice sampling. The Journal of Machine Learning Research, 15(1):0087-2112, 2014.
Wang. Y. Parallel Numerical Picard Iteration Meth- ods. J. Sci. Comput., 95(1):27, 2023. doi: 10.1007/ $10915-023-02156-Y URL https: //ti:./////.i////ooo org/10.
Oksendal, B. Stochastic differential equations: an intro- duction with applications. Springer Science Business Media, 2013.
Oliva, F. V. and Akyildiz, O. D. Kinetic Interact- ing Particle Langevin Monte Carlo. arXiv preprint arXiv:2407.05790, 2024.
Wibisono, A. Sampling as optimization in the space of measures: The Langevin dynamics as a composite opti- mization problem. In Conference on Learning Theory, 2018.
Ong, B. W. and Schroder, J. B. Applications of time paral- lelization. Computing and Visualization in Science, 23: 1-15, 2020.
Wibisono, A. Proximal Langevin algorithm: Rapid convergence under isoperimetry. arXiv preprint arXhv:1911.01469 2019.
Otto, F. and Villani, C. Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequal- ity. Journal of Functional Analysis, 2000.
Yang, R., Srivastava, P., and Mandr, S. Diffusion probabilis- tic modeling for video generation. Entropy, 25(10):1469, 2023a.
Pedrotti, F. Maas, J., and Mondelli, M. Improved Conver- gence of Score-Basee Diffusion Models via Prediction- Correction. Transactions Machine Learning Research, 2024. ISSN 2835-8856.
Yang. Z., Feng, R., Zhang, H., Shen, Y., Zhu, K., Huang, L., Zhang, Y., Liu, Y., Zhao, D., Zhou, J., et al. Lips- chitz singularities in diffusion models. In The Twelfth International Conference on Learning Representations, 2023b.
Robert, C. P., Casella, G., and Casella, G. Monte Carlo statistical methods, volume 2. Springer, 1999.
Yu, L. and Dalalyan, A. Parallelized Midpoint Ran- domization for Langevin Monte Carlo. arXiv preprint arXiv:2402. 2024.
Salmona, A., De Bortoli, V., Delon, J., and Desolneux, A. Can push-forward generative models fit multimodal dis- tributions? Advances in Neural Information Processing Systems, 35:10766-10777 2022.
Zhang, Q. and Chen, Y. Fast sampling of diffusion models with exponential integrator. In The Eleventh International Conference on Learning Representations, 2023.
Shen, R. and Lee. Y. The Randomized Midpoint Method for Log-Concave Sampling. Advances in Neural Infor- mation Processing Systems, 32, 2019.
Shih, A., Belkhale, S., Ermon, S., Sadigh, D., and Anari, N. Parallel sampling of diffusion models. Advances in Neural Information Processing Systems, 36, 2024.
Song, J., Meng, C., and Ermon, S. Denoising diffusion implicit models. In International Conference on Learning Representations, 2021.

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
A. Useful tools
A.1. Girsanov's Theorem
Following the notation introduced in Chen et al. (2024, Appendix A.2), we consider probability space (Q, p) on
which (E)) is Wiener process in R with the filtration (Fe 120. For an Itô process zr(w) satisfies dz)(w) all,w)dt + (w),,, we denote the marginal distribution of 2 by Pt, and the path measure of the process 2e by Ptytt
Definition A.1. Assume is the Borel o-algebra on IRd For any 0 S 4 < 12, we define V as the class of functions f(t,w): : +00) x f2 + which is F, measurable, Fe-adapted for any 2 and satisfies
exp (T, ((t.)att) < +00,
AA > 0.
For vectors and matrices, we say belongs 1 or vmxn if each component of the vector or each entry of the matrix belongs V.
Theorem A.2 ((Chen et al., 2024, Corollary A.4)). Let axt,w) E Vmh E(t, E vmx" and (wi(w) be Wiener process on the probability space (2, F.g). -)- For E [0, 7), suppose zi(w) is an ltô process with the following SDE: (z)(w) = (tt, E(t,
(6)
and there exist processes w(t,w) E v and ((t,w) E v such that:
1. (t,,w) = axt,, w) - ((),w);
2. The process Me(w) as defined below is martingale with respect the filtration {FiJe2o and probability measure q:
Mc(w) exp (-f. (''(ss.o) 8(s,w) w)
((s,w)Tdw,(e dw, 1// 100 [0((s:))ldds)
then there exists another probability measure p F) such that:
1. p < q with the Radon-Nikodym derivative
d2 dq (w) Mr(w),
2. The process we(w)) as defined below a Wiener process on (9t F.p):
wr(w) = wr(w) L' 8(s,w)ds,
3. Any continuous path in C(((,12) 12], R") generated by the process zt satisfies the following SDE under the probability measure
dZ+(w) = B(t, + Ett, ()).),(w).
Corollary A.3 ((Chen et al., 2024, Corollary A.5)). Suppose the conditions in Theorem A.2 hold, then any anyta, [0, with 11 < l2, the path measure of the SDE (7) under the probability measure in the sense of Ph1tz = - (-)) is absolutely continuous with respect the path measure of the SDE (6) in the sense of Itrit -2( Moreover, the KL divergence between the two path measures is given by
[ n (1((t,))l
(L (1g))sss) = (L(,,, ) + Eunoll
2
A.2. Comparison Inequalities
Theorem A.4 (Gronwall inequality (Dragomir, 2003, Theorem 1)). Let x, Vi and x be real continuous functions defined in Ja, b). x(t) 2 0f E [a,b]. b). We suppose that on [a, we have the inequality
x(t) S (tt)+ S' x(s)x(s)ds.
Then
x(t) s (t) + L (s)) exp [ x(w)aa) ds.
12

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
A.3. Help Lemmas for diffusion models
Lemma A.5 ((Chen al., 2023a, Lemma 9)), For do La) and pr is the distribution of the solution the forward process ((2), we have
KL S de-
B. Missing Proof for Log-concave Sampling
denotes KL = KL(H,, where Ha represents the law of a., define the truncation error from prior update
as, E max E
||| | - 2h. !|?)
and truncation error at right boundary point as, A := ===[]l [[[[== - -
M
B.1. One Step Analysis of K From KL's Convergence to Picard Convergence
In this section, we use the interpolation method analyse the change of KL along time direction, which will be bounded by discretization error and score enror.
Lemma B.I. Assume 3h S 0.1. For any j = 1,
= I,..., 1, we have
0.53dh
KL S (xxp(-1.2c/))Kr
+ h.4 h.15th.
M
Furthermore, for initialization part, ie., j = 0, n = 0, O,...,N N 1, we have
KL exp (-a(n + 1)) (L(PolI 832
Remark B.2. the first equation, the term exp(-1.2a.)) KLL characterizes the convergence of the continuous diffusion. reverts to the classical
Additionally, the second and third terms quantify the discretization error. Adopting P = 0 and scenario, where the discretization error approximates ((dd), as discussed Section 4.1 of Chewi (2023). Moreover, the second term is influenced by the density of the grids, while the third term dependent on the convergence of the Picard iterations. The fourth term accounts for the score error.
Proof. We will use the interpolation method and follow the proof of Theorem 13 in Anari al. (2024). For j E [J]. = .....,, NN
and m 0. M - 1, it is easy to see that
h
xxd a.m
s( + V2(Bnh++(++1)/ Buh+thh/M)-
M
Let æ, denote the linear interpolation between w
2,,+1 and wh.m i.c., for E [nh+ + nh (m++1ihh 1ihk
let
mh s(xg)ft))
(t-nh nh
V2(B, Bah-mh/M).
M
Note that s<)); is constant vector field. Let ur be the law of xt. The same argument as in Vempala & Wibisono (2019,
Lemma 3/Equation 32) yields the differential inequality
FI( [|)) + TE ())()))))))))))))))))))) s(ar).P-1), V log (2(
3,KL(u)|lm)
π (xx)
3 FI(ue (T) - |fx;))
s(x) ((2)
(8)
where used (a, < ||aa|12 -|b|| 2|b|| and E of have (L(
V log
Me (2:) ((x;)
FI(pll(T). For the first term, by strongly-log-coneanl
< 2m F1 (T) For the second term, have
V f(x)) s( ( 1|
E
2E
[ v ((( >|
2E [ (1t) s(ar..P-1) 11
< 2
||20
- -.PP- |1(] 28
(9)
13

---

Parallcl Simulation for Log-concave Sampling and Score-based Diffusion Models
Moreover,
E
||
n.m
P-P s 2E
11]
|| x
=mmmmll 11
+ 21E
n,m
11 !1?)
(10)
For the first term, which will be influenced by density of grids, we have
E
x
mh
mh M
- nh-
EE
s (122))1)1)1
+ d
t
nh
M
( }1? 11
mh
E
+
t - nh
M
M
2h
E
||
1)1
dh
M
M
M
h/2
E
P.P-1
[1]]
+
4h2
E
|IV )
28 dh
M +
(1I)
M
p,m
M
M
Taking Bh < 10 and combining (10) and (11), have
E
||
P.i.P-11
||]
4.4
M2 E [V
2.25
M
+
1.1dh
M
2.21E
|| -| kh
2n.
P...P-1 P-1
(12)
For the first term, we recall the following lemma.
Lemma B.3 (CChwwi et al., 2024, Lemma 16)).
E [[V Fl( [| 2,dd
Combining (8), (9), (12) and Sh <
10 we have for E [JI = O,..., 1, m 0.... M - 1, and
1
EE nh + má M nh (m++ 1ihh
0,KL (14 (TT
s FI(uelTT)
3
[|V]((x)) (
111"
s 3
FI( + 2
|
--P-111 11]
28
S FI(
3
h.
E
)|
4.43282//2
2.23t +
+ 4.4.3²
2.P-111 -11] + 282
M
M2
M
n,m
0.1
0.182
dh
< 4
+
VV()))
+4.4,2
+02
M
M2
M
3
+ 0.1
0.15
2.2
< FI(Mell=)=
M (FI(ee
[|)) 23d)
4.4 +82
M2
M
- - 1.20KL(44]I)) + 0.53d
+ &. + 2.18
M
Since this inequality holds independently of m, we integral from -= nh to (n 1)h,
KL < 1xL + 0.5,dh
+ 4.4,2 + 2.18
M
As for j = 0, actually, Line 4-7 performs Langevin Monte Carlo with step size h, by Theorem 4.2.6 in Chewi (2023), we have
8dh,32
KL s exp (-anh) KL.
with < h 1b
14

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
Corollary B.4. Assume 3h S 0.1. We have
N-1 Σ
0.5,d
2.18²
KL S e-1.2a(N-1)6 (KL(H))1))444337))) +
e 1.20(n
4.4
+
aM
a
==1
Furthermore, E has uniform bound En-n S + 5506 we have
0.53d
2.50
KL S se" -122 ( + +,3KE +
o.
Proof. By Lemma B.1, we decompose KL KL as
N-1
0.5,3dh
KL -1 Se-120(N-1))KL
Σ
E 1-220(n-1) 1)h
+ +4.48 +2.10% M
n=1
Se Se-1.2aNN11) 11th (KL(HO]ITT
+ h44
(E + 50082 1/2)
0.58dh
+ -2.182h M
+
1- expxp(- (exx(-1.22x))
Se 1.2a(NN1)4 (KL(HolIT) 4.4,8 + 1.1
1.1 0.53dh ah M
4.4 +
ah
1.1
+ ah
2
1120 (KL( + +4443%10)) 5m,38 + 0.63d
2862
aM a
where the third inequality holds since < x < 0.4, we have 1.1 1.1 exxeepp112..0 - x > 0. It is clear that ah < 8h < 0.1.
B.2. One Step Analysis of A
In this section, we analyze the one step change of A first.
Lemma B.5, Assume 3h to and P > 2log + 4. For any 2......J, 1....,,, N 1. we have
A 0.005
1
10x 1282 + 4.4 (
A A 4.4
+ 10x
M
M
Furthermore, for 1, n I,....
1, we have
A 1 ( M + 10x)
KL
(582)2 d. 0.4,3221
a
Proof. Decomposition when j 2 2. In fact, for j E [J], n O...,,N- - 1,
O,..., M 1, and 1,..., P. it
easy to see that
p m
zdP nmm
+ v2(Bnu++(+11/// Butt+mhiMa).
M
n.m
For any j 2...., 1, N - 1, the contraction of x(x) = (2023)), for any m 1,...,M M. we have,
fx) (Lemma 2.2 in Altschuler & Talwar M
h
E
P.
1-1
E
J.P. i.P
nnm-1
M
h
swx (P))))
j-I.P
M
h
[
n
1.P-11
15

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
$
(1 + B)
J.P
[f(a
n,m-1
2
11-1P n.
mm
n,
M
M
h
h
+ 2 +
E
M
S.P
m,m-1
M
V
f(
+
M
V
2 n.
M
f 2
2,m-1
+
2 +
E
M
h
f
p,m
M
8
n,m
) +
M
h
f
n.m
M
h
s(
P111P-11
n.m
(1+1)) (1
ah
M
2
E
p1.m-1
2 mn.m-1
I,P
2
+ 4
"
4
M2
hr2
6²
4
1/
E
p
PPP1 |/]
4
13212
E
--1-P-11
M
M
n,m-1
By setting
17
ah M
10 we have
27
E
J.P
x
11
M.M
R.M
ah
M
x-11. 2
h2
M
E
n.0
n,
+
44 44
M 82
+
M
4 + 4
h²h2
E
P.P
-1
2
M
n,
n,m-1
m=1
M
4
h2h22
1-P
--1.P-1 2
4 +
17
M
E
11.m11
n,
m=1
S exp(-aah)
+
(4+
4
h
82
1+ 4
h2
E9 +
4
4
hhh
11-3
7)
M
9)
M
1)
M
(++
4
112
B2h
4
/2 E1-1
S (1 - 0.1ah) A
+
82
4
E
4
17
M
M
M
(1 0.01
1
A 2 +
10k)
11882 4
10x
K
M
M
+
10k)
(13)
M
In the following, we further decompose E/e For any n = O 1. E 2........... P. and m 1,..., M, we can decompose E
|| -h
[|3]]]
as follows. By definition (Line 12 or 18 in Algorithm 1), we have
n,m
||
'1"]
k,m
n.m
h²2
m- 1
m-
E
Σ s(arhp-))
M2
n,
ml=0
m²=0
h2m
M
m-1
E
s
a p - s( (c00hip0 2-2
nimm
1
mi!
12m
3
E
f( (1) (. - V (x
27
+E
V
>1
27
M2
mn 0
E
a ( 2
ppppp- 2
< 23
max
E
+ 682
(14)
m' Il......M M
16

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
Furthermore,
E
2.1
23.0
n.m
23,m3
h
m-1
h
m-
E
1-1
Σ
n 1.M
n-1,M
M
h2mm
m-1
- ||
S 2E
n- 1--M 1.M
x n-1.M
M
E
s
22 + 6, + 1282
(15)
Combining (14) and (15), we have
P-P 27
Σ =E
n,m-
J.P
nm
S
2.0.03 A n-1 + 0.03 - 6..
(16)
Substitute it into (13), we have for any
=2, .,.....Inn=
1.
0.01
1
0.03*) A 1/-1 + +44 (
A S
1
8
+10k
+ 10xx 1282
M
M
+ 4.4
1
+
10x
18112
(17)
M
0.005
+ 4.4 ( + 10k 12822 444 (1 10n 82
A
(18)
k
R-
M
where the second inequality holds since P> 2 + implies ( + 10k) 0.03p s 0.005
3
Decomposition when j 1. When j similarly, we have for p I...,,,.
x 1p n,m+1
(
s (1-1-- + V2(Bm) Bat+mh/MM):
M
and
x n,
xb
s s(x ( + (22(B)h Bnh+th/MM)
n.m
M
Thus by the contraction of ((x) = x
fxx) (Lemma 2.2 in Altschuler & Talwar (2023)), we have
M
E
D n.m+1
1.P
h
h (
E
D n.m
L.P
s(x
s
M
M
h
h
27
<
(1 + +1]
L.P
(x
M
M
2+ 2
h
h
h M
h
E
V
(
[ n 1.
M
M
M
2 2 E
h
h M
h M
h
12
s
1.P-
V
s (.._1
1
M
M
< <11 (1+m) m))
(1
ah
E
x
I.P
2
+ 4+
82
M
2
+ +(4+ 4+ 4
7
182A2
M
E
1.P
1.
1112]
+
+44+
44+
n
4
5²h²
M
E
||0
n.m
x
1,M 112
17

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
For third term IE ||
111-11
we have
I.P
1-122
n,m
2
E
M
h
m
8 8(x
(10 1.111
n,m!
8(28 n,m!
1-P2
m
=0
mh?
M
s( 1.
s s n2 P-2))
|1/
E
M
$ 3,32 h2
max
E
a P
n.m"
2 n,m
1.P-2 2
+ 162
m -0..,,,
Thus
E ||
S 0-03
E
n.m
1,1
1,0 p.m'
+ 6.20%
(19)
O....
For E
|1
x p.
1.0
[12] by definition, we have
1.1
1.0
||
2 n.m
a
m,
m
E n
1,M
M
T 1n.-1,M 1.M
M
Σ
8 x
1.M)
m
mi 0
$ 2E
n 1.M
2 1. 112]
+ +2E
h M
m
m-1
s(x
n.m")
M
h
m
s(x -
m
12mm
m
$ 2E
n2 1.M
2 12-11.011 113]
+
M
E
s (f
s(xx 1.11 ) |1(
m' =0
< 2E
|| 20 -.
+ 1/
max
E
10 n,m"
x0 20-11,M11 11 + 1282
(20)
E[M]
For E
|1
w
-0-1,,r by definition of ahh (Line in Algorithm 1), we have
2 p,m
1MM 19
2
12m2
|||5(10)) 1,M, ( [1²] M
dhm
E
M²
282h
2
f [12) dh
482 23d
282h2
2h
(1 + dh
8 22 +
41.2dd +
KL - + dhi,
(21)
where the last inequality is implied from the following lemma, (Vempala & Wibisono, 2019, Lemma 10)
4
23d + (L ITT).
Combining (19), (20), and (21), and P > have
E ||| n,m
22.P-111 [1]]
18

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
< 0.033--1
maxx m' ml=......M M
E
2 1,1 n.m"
2 1.0
n,m!
+ 6.
+ 8/3112
0.003
24 n-1 + 6,321
[4h23d 22
KL dh + 1200
6.21282
22. 1..00-111/...000.00 + 6.3 +.00ddh 0.01
KL
0
(22)
n-1
By setting 77 = 33
10kM we have
E
n.M
ah
M
M
E
1
n,0
I.P.
2n.0
+ 4+ |
11
+++
8212
M
4
82h²
+
1
4
2. 2.0.03 P.
A + 6. + 0.01dh + 0.01
KL n-
9
M
4
1
8,32/2
+
4
2 3d
+
222822
KL + dh
7
M
0.01
+ 4
10ks
0.03P)
A 1 + (1
+ 10ks
552h2 + 6,3 dh + +0.433%2
KL
M
M
(1
1
+ + 0.4 KL
< - +
+ 10k)
M
where the last inequality holds since P> 2
implics ( + 10x) 0.03" <
0.005
When n = 0, the update is identical to the Picard iteration shown in Anari et al. (2024), thus we have the following lemma. Lemma B.6 (Lemma 18 in Anari et al. (2024)). For j l....,,.. J, we have
A/6 $ 0.03 A h.2
with A max E || xo ||2] < 18 KL (L))ll) + 1.4dh + 25
m=
Corollary B.7. For n = I,...,N
N
1, we have
A S n ( M
5.118P22 + 0.5 g2h2 goh (L(Hul(T) + 10n²3"an²)
+ 10k)
a
Furthermore, for j 1,...,J and n = 0, we have
A < 0.03 4.2
a KL(Hal|T)) + 1.4 0.03 dh 6.78
Proof. By Lemma B.6, we have
A 5.003 -/-1 1-1 6.28
$ 0.03
A + 2
S.0
KL(uol|T) 1.4dh + 282/22 +20"/3 h2
< 0.03IP 4.3%/²
(L(Pa]IT)
+14 P.03'P dh 6..7
a
19

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
Combining Lemma B.1 and Lemma B.5, we have
KL.
A
< 0
+
Σ
M + 10k
55 + 6, + 0.4,3%
a
t=1
$ 0
+
+ 10k)
(58² + 63 dh3))
M
11
1
0.4 8
exp (-amh) (L(uol|=) + 8,3
+ Σ
+ 10x
M
a
1
56 + dh 3 + 0.4 3212
$ 4 +
+ 10k
KL(Holl=T + +322K² B²
M
5.16%h2 + 0.5 hh
Sn
n
+ 10k²
(L(Hol|T) + 10x² 3222
M
B.3. One Step Analysis of E/
In this section, we analyze the one step change of E Lemma B.8. For any j = ................................ J. n = 1..... N 1, we have
E S 0.03 " + 2 2.00 g 182
Furthermore, for n 1, ,,..., N 1, we have
h2h²
E, 2.0003 P-1
6.3 + 0.01dh + 0.01 + -
KL n-
Proof. By (16), the first inequality holds. By (22), the second inequality holds.
Corollary B.9. For n = I,...,N - 1, we have
N
sh E (..58 + + KL( [()) 0.1 dh)
Proof. Combining Lemma B. Lemma B.8 and Corollary B.7, we have
E 2- 0.03 P-1 A 6.3h + 0.01dh 0.01 8h2
KL
a
2- P.03. - n- + 6. + 0.01dh
8422
eep (xp(-a(++1))) 1)4) (L(oll 833th
+
22.0.03 AAL n-1 + 20 + 0.02kdh + 0.01 h2
KL(HolIT))
2 22.0.03P-1
(1 1
+ 10k
5.162 + 0.5
32
KL(Hol(m) + 10x² 30n23²an))) )n²))
M
a
8212
+ 26 0.02kdh + 0.01
(L(Ho](T)
a
n 0.06 h.112 + 0.5 82
(L(polIT) + 0.1x2
1211
6. + 0.02kdh + 0.01
KL(Hol|TT)
8² Sn 5.55 0.1 (L(oll) 0.1
a
where the fifth inequality holds since >> 2 + implies ( + 10k) P. 0.03.
20

---

Parallel Simulation for Log-coneave Sampling and Score-based Diffusion Models
B.4. Proof of Theorem 4.2
We define energy function as
14 - + KEj-1
We note that 2- 0.003 12 + 78 2 Eg By Lemma B.5 and Lemma B.8, can decompose L2 as
LA = A n-1 + KE
1 0.005
1
1 M
A n-2 + 4.4
+ 10k
1282 + 4.4
+ 10k
M
+ <0.033
N +2 0.03 E +8
1 0.005
0.005
2 + 1-
8J
K 0.03 0.03--1A/-1 + -=00333 -x-003--1..
-KEj- -
R
k
+ 4
0.005 =
- + (x 0.13 - 56
(23)
K
Combining P
k + 4 implies K. 00.33 < 0.04, we recursively bound L as
+ 2) 0.04
(R. 0.03 (1 0.005 k
n-1
n n--++
1+ b) j-b
0) I
T2
1
j 2
0.001 n-b
0.04 165m62112
0 2 b=2
0.04
11
2)
1
(x. 0..0003
0.005 0103--115-6
n-1
n 1 + j 6) j b
I
j --2
16
b=2
280)
(24)
For the first term
Σ 0.04 0.043-2(m-92))
j-2
L. first bound 1. To do so, we first bound A follows. Combining Lemma
B.5 and Corollary B.9, we have
0.005
(1
[11
1
A3 < 1
A2
4.4
+ 10x
12802 + 4.4
+ 10x)
K
M
M
82h²
< A - + 48. 20
+
488.4x3311
5.58% + +0.1 11
(L(HOllT) + 0.1n2dh)
a
8212
< A + 48. 48.4
55.58
0.1
(L(Ho|[=) + 0.1k*dh
a
8212
< + 18
55.50
+ 0.1 KL(PolIT) + 0.1x
2P 4
KL(POl|T (T) + 1.4 O. dh + 6
0.1 821²2 KL(μoo| [T) + 0.1
55.55% 18.2
a
48.4n33222n 67 0.2 32 KL(HolI=) + 2
a
21

---

Parallcl Simulation for Log-concave Sampling and Score-based Diffusion Models
Thus
L = - + KEE
3%1
S 0.499k(a - 1)2 67.26" + 0.2
KL [1) + 0.2x22hh
KL(uolITT) + +0.1x2dh) )
+ a 5.500 h2 + 0.1
+0.1
h2
S na² 12 0.2
(L(uo]|=T 0.2x²ah)
Thus by ((n) 1 cm
for m 2 > 0, we have
a Σ 0.04-2 2
aa++ -
L
j-2
11 a +j - 22 J-2
0.04 22--
L3
j j-2
2=2
<
Σ
0.04 2 22-411² L[2
a=2
a
321²
0.3/ 2na²
33982 + 0.2 a
(L(pol| + 0.2x2
82h²
< 0.33 kn
h9 + 0.2
KL(HoolT)+0.24 + 0.2xx
(25)
a
For the second term
(R- 0.03 (1
0.005 x
n-1 (n-1+j-b) LG, we first bound L Firstly, for E combining j-b
Corollary B.7 and Corollary B.9, we have
gb-1 2.00.03 P2
2- 0.03 E gh-2 + 7821/2 2-.0.03 (0.03 432 4/22
(L(Hol[T)++ + 1-4 0..00. 1)p dh dh 6.78m2
+..0.03 8 1
0.01 48 (L(Holl=)+000100)
<2-0.03 + 0.03
+ 22
a
(0.08-
00.01 h8
<22-0.03")--21- +
(2 23-0.03")
(L(Hol|T) + 0.01dh + +7.16%1²
a
0.01"0.03' (0.03~-1 0.01 4.
<(2-0.03P)--22 2-2 +
(L(pol[=) + + + 7.182/2
821²
< (2-0.039))-
5.58 2 +011
a KL(Hol| + d.
0.00 4
0.03
KL(HO 0.02dhh 2.282h2
a 82 0..1
$0.03'
KL(Hol( + 0.1dh) + 2.
a
As for 4/ we have
A < 0.03hP 4/2222 4
KL |T]] + 1.4 0.00 dh - h.
22

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
Thus, bound the first term
1 4 + KET-!
- 0.03 482h2
(L((oll=) + 1.4 0.03 dh + 6.7
0..1 (0.1 82h
+x0.03'
(L((oo|=)) + + h.
(0. 8h² KL(wo + +..2dh
<m0.03*
148%
a
Thus by (*)) I ( (em) x2
for
0, and
(n+t) 1
x²
1-mm+1] m++++1
(1-2)=+1 we have
(1-2)
(N. -.00003
P.03P-192
11000 (1 0.005 n-1
+ - b)
) L
j b
0.04
0.005
n-1
n
1 + -) (NO.03* (0.2 (0.2 8h2 j-b b
(L((ollT) + 0.2dh
a
K 0.03 0.03--190~6 (1
0.005
-n- 1 b)
0) 1402/2
j-b
b=2
(n (n-1+3 --)
-b)x
b2
Σ 0.04
0.2
KL || + +022dh
j -
a
b=2
(R: 0.03P-195-1
(1
0.005
1 + - n
1452h2
j b
b=2
1-2
1+ n )
<
Σ
e
1+
0.2
KL [|]] + 0.2dh
a
j-2
R 0.03P-10 (1
0.005 K
n 1 1+ a-1
+
Σ
1462
i
1=
0.2 (0.2892
821/2
< ~0.111e-11 Ri
(L(Hol|T) + 0.2dh
0.005 n-1
+
(6.6+ + +
(1 - 110m
82h²
<0.11'e-116 K
0.2
(L(Hol(T) + 0.2dh) + x1 (1-3-0.03P-1) (6.6 7.9k)02²
11221
(L(HollT)+++ + 1.6dh + +200h²) )) + 20k02
50.111
2.2k
where the second-to-be inequality implied by 8 ( + 10k) 0.03P < 0.005
23

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
Combing (24) and (25), we bound Lg
1 -
n1
0.04
a
2
L
Σ
(x-
0.005 0.03 (1
n-1
nn n 1 --b
-
j
b=2
+ 68000x282
321.
x 0.33~2kn* kn
190022 + 0.2
(L()oll= + -.22k2th
1/82
+ +
2.2K
KL(uol(T) + 1.6dh + 202 + 20x=2 68000x 2
00.3* 20"=1 kn
(4162)²
1.8 dh 0.5h/KL + 68020x2521z 26.
Since 8 (M 10x) 0.03P=
0.005 implies x=00.0 S 0.003, we have
< 0.03 P-1 1/ h
22-0.03 (0.33-2 2 e kn (418 1.8 0.5xkLL(HOlIT)) + 68020x252 + 78
S 2 18 + HKL((HOTTT) + kdh) + 2
Thus when J - N 2 log (N3 [nPhtunk(nolltttttt)
we have for any n 0, ......N
ex4 - 5x3
+ 11
Recall
KL S 1N (KL(HolIT) + h44 + 3x38+ + 0.63d
288²
a.M +
thus when 8² aga 29
.M >
xd ş2,
and N > 10k log KL
we have
0.6,dd
2882
KL x- I -1.2 - (KL(Mo)IT) + +4.482haa) + 5k.8
a.M +
a
(0 4,3212
$ Se 1.2 (LLHll)) + 4.4,3%
KL (L(pol| (T) + 1.4 0.03 "h +0.7877))
+ 55.EE+ + 0.63d
2802
aM
$ 11.20 11220(N-1 [T) + 5kBE + 0.63d
2952
aM
a
< 3
C. Missing Details for Sampling for Diffusion Models
In this section, we begin by presenting the algorithm details in Appendix C.I. In Appendix C.2, following the approach of Chen et al. (2024), we apply Girsanov's Theorem and the interpolation method to decompose the KI. divergence and bound the discretization error, accounting for the influence of the step size scheme and the estimation error of the score function. Finally, in Appendix C.3, we analyze the additional parallelization error and derive the overall error bound.
C.1. Algorithm
In the parallel Picard method for diffusion model, we use the similar parallelization across time slices as illustrated in Figure 1. In Lines 2-6, we generate the noises and fix them. In Lines 7-10, we initialize the value at the grid via sequential method with a stepsize h2 = ((). In Lines 12-21, we update the grids diagonally, using the exponential integrator in Lines 14 and 19 instead of the Fuler- Maruyama scheme. The step size scheme also differs from that used for log-concave
24

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
sampling. Here, follow the discretization scheme with early stopping and exponential decay described in Chen et al. (2024, Section 3.1.2).
Algorithm 2 Parallel Picard Iteration Method for diffusion models
1: Input: Do To (0, Ia), the learned NN-based score function (+), the depth of Picard iterations J, the depth of inner Picard iteration P. and discretization scheme (T, (Nm))=1 and (n.m)ne(0:NN11)mm000
2: for = , N - do
3: for m = 0, ...,, M (in parallel) do
4: Sn.m ~ NN In)
5: end for
6: end for
7: for 0, ...,, N do
8: for = 0, .... Mn (in parallel) do
9:
V 1,M Do. for j = 0. J.
Tm.m
V Vn.Ta.
In-1.7m.h
m
'nmm
2(c5mmmm 1) - ( + Vetam 16m
(26)
102
10: end for
11: end for
12: for k = I..... N do
13: for j = 1,.... min (k- 1, } and m O,... Mn (in parallel) do
14:
let k and 0n V n
F fn.m
Vn,0
m-
(27)
2( 1) 1) (Frroomr + Vetam - -T5m²
+
m'=0
15:
end for 16: end for
17: for &- = N + 1....,N N + J do
18:
for n = max (0, t-- J}. N and m = O...., M, (in parallel) do
19: let = k - and Uno y M.M.
V Tn.mam
Fn,0
m-
(28)
2( 2(e 1))lnttnm (Distamm )+ + Ve - IImm
+
m 0
20:
end for end for
21:
22: Return: of
UN-1-N--1
Stepsize scheme. We first present the stepsize schedule for diffusion models, which is the same as the discretization scheme in Chen et al. (2024). Specifically, we split the the time horizon into N time slices with length h, S h
N = 2((1, and
T
a large gap grid (t,)) with f, = h;- For any n E [O - 1), we further split the n-th time slice into grid (mn,, M=o
with Tn,0 0 and M, hn- We denote the step size of the m-th step in the n-th time slice Ca,m Tn,mt1 Tn,ms and let the total number of grids in the n-th time slice M,. The grids Ma scheduled follows,
1. for the first N - 1 time slice, we use the uniform discretization: for 0, O,.... N 2 and m 0,..., M 1.
h = h, Cn.m e, and M,, = M = h
25

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
2. for the last time slice, apply early stopping and exponential decay:
hN-1 = h 6, 6N-1.m S SEAE TN- TN-1.m+1..
In(r)
We also define the indexing functions as follows: for T E [tns tntil, we define In(T) E N such that
Σ
Cn.j.
S T <
j=1
1 2 eng. We further define a piecewise function g such that (T(T) = Σ
1.fr)
En.j and thus we have LT(T) = LT/E] and
j
(T) LT/ele.
Exponential integrator for Picard iterations. forward Guler-Maruyami scheme for Picard iterations, we use the following exponential integrator scheme (Zhang & Chen, 2023; Chen et al., 2024). Specifically, In n-th time slice [t,, in + Tn,M,, J. for each grid in + Tn.m. we simulate the approximated backward process (3) with Picard iterations as
Compared with Line 12 and Line 18 in Algorithm 1, where use
Fm.m +1+1 Yn 1,
Vn.ua.
2(e5n - -1)s t
tT (Fn-I,t + Ve'nmm - 16m
192 0
We note such update also inherently allows for parallelization for m 1, I,..., Mn.
C.2. Interpolation Processes and Decomposition of KL Divergence
Following the proof flow in Chen et al. (2024), define the following processes:
1. the original backward process,
[ Ft + log (f. + dw., with To pri
dē
(29)
2. the approximated backward process,
dy: こリ + 8. 8(())) di dw., with Yo ~ NO Ia);
1
3. the interpolation processes (i (T,,,)) rE(.M) time & Fe, for any fixed O..... N 1, I,...,J,
over T [0, h] conditioned on the filtration of the backward SDE (29) up to
2 ( + sl.ton(t) (Di-1 (Binn)))t))())) (w)) dr dwe,+r(w)
1
d (w)
(30)
with Vi. (ww) Vh 1. 1 (w).
4. the initialization process,
1.0 tot(()))
(h (M (wi) dr + dww,.r(w))
do. (w)
5 Uin. (w) s
(31)
with Vto = Do and Vt y
Remark C.I. The main difference compared to the auxiliary process defined in Chen et al. (2024) is the change of the start point across each update.
We can demonstrate that the interpolation processes remain well-defined after parallelization across time slices.
26

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
Lemma C.2. The auxiliary process (F
(W)) (e))mn addpeed for any 1, ....... .1 and ,....,
Proof. Since the initialization Vi,,. (w) satisfies
(w)) = (w.) + toa()) (Vl (( dr + dww, +T(w),
1 +0
we can claim Pi.,.r (w) is aadpted Assume ye..r -ddpted by In(T) S T, the Itō integral So +,,4,,77) (n ddr is well-defined and -rdapied Therefore SDE
dy' (w) 1 / (w) + tl. (Umn) dr dwe,+r(w)
has a unique strong solution (Vi,. (W)) that is also Fr,, addptee.. The lemma is established through induction.
Finally, the following lemma shows the equivalence of our update rule and the auxiliary process, i.e., the auxiliary process an interpotation of the discrete points.
Lemma C.3. For any .,...,N. N - 1, the update rule ((26) in Algorithm and the update rule ((27) (28)) equivalent to the exact solution of the auxiliary process (31), and (30) respectively, for any = I,....J,, and [0, hn).
Proof. Due to the similarity, we only prove the equivalence of the update rule ((26). The dependency on w will be omitted in the proof below.
For SDE (30), by multiplying & on both sides then integrating on both side from 0 to 5, we have
M.
e --
Pl
2
AAnn m++1
2
battamm
(Dlnttamm + /- e dwo, +rl
m=
Thus
2
Von.o
2
TAra.m Thrm.mt1
-1)
over mmmmtl) 2
8 -mm (winttmm)
m=0
M
TATmm+1
dw,,, +r'
tr'.
By Itō isometry and let = Tn,m we get the desired result.
C.2.1. DECOMPOSITION OF KL DIVERGENCE
Similar as the analysis Section B.2 of Chen al. (2024), we conclude the following lemma by Corollary A.3.
Lemma C.4. Assume 6 (T,w)= o tn.ts decomposition,
@--1 ( (w)) log Pr, + (w)). Then we have the following one-step
1 hnn ||8.,, (T,w))12 dt
KL l1 S KL (1T...)
Buong
2 J0
27

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
Now, the problem remaining is reduced bound the following discrepancy,
| (T,
Am
V llgg P (i) dr
Pn
S
log Poo
V log +r (T (w)) dr Aen (w)
8 ( (T) V Vlog Pt,, (n)(t) (T) (w)) dr
(
:B
hn (
se := (w)
dr ( (wW))
(32)
0
where A,,, (w) measures the discretization error, Bi, (w) measures the estimation error of score function, and Ct. (w) measures the error by Picard iteration.
C.2.2. DISCRETIZATION ERROR AND ESTIMATION ERROR OF SCORE FUNCTION IN EVERY TIME SLICE
The following lemma from Benton et al. (2024); Chen et al. (2024) bounds the expectation the discretization error A...
Lemma C.5 (Discretization error (Benton et al., 2024, Section 3.1) (Chen et al., 2024, Lemma B.7)). We have for n 00: N
Ewwilr [Ae, (w)] edhn,
and
EuwBlren [Aex_, (w)) S edlogn
where i) is the parameter for early stopping.
The following lemma from Chen al. (2024) bounds the expectation of the estimation error of score function, Bt,, and we restate the proof for the convenience.
Lemma C.6 (Estimation error of score function (Chen et al., 2024, Section B.3)).
Σ E
Bl [Be,]] < 52.
n
Proof. By Assumption 5.1 and the the fact that the process (w) follows the backward SDE with the true score function under the measure i. we have
Euw
[Bn, (w)l
<. Euw-pl Fmm
bt (V) - V log Dr + (o) (o (r)
(w)) dr
Mn
s +x + (Vr.. (w)) V log Dr +on(r) (Di (()))) | dr
= EEw Fen
En.m
n=1 m==
N-1 Mr
V (t) + (w))) dr
-Eu-fils
Cn.m
n=0 m=
< 82
28

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
C.2.3. ANALYSIS FOR INITIALIZATION
By setting the depth of iteration as J 1 in Chen et al. (2024), our initialization parts (Lines 4-7 in Algorithm 2) and the initialization process ((31)) are identical the Algorithm I and the the auxiliary process (Definition B.1) in Chen et al. (2024). We provide a brief overview of their analysis and reformulate it to align with our initialization. Let
A², (w) =
V log Pt., ta. (T) (o (r. (T)( - V log ( (w)) dr
and
B1, (w) //0 Nnn || s ts. (T) (+) V log tt. ( on (sWW)) dr
Lemma C.7 (Lemma B.5 or Lemma B.6 with K 1 in Chen et al. (2024)). For any 0, .N - 1, suppose the 0 initialization o follows the distribution of Tt. ~ Pen. if th In < 0.5, then the following estimate
20 olo [1 (w)
sup Eu rE(D,A.,
e ha (Ma + 2d)
+ Ge hn Eum [A (w) + BI, (w)
Furthermore, the 49, (w) and Bi, (w) can be bounded as Lemma C.8 ((Chen et al., 2024, Lemma B.7)). We have for E 10: N 2]
Euwilyrn [A (w) S edhns
and
AAl_ (w)) S edlog
Bu Fen
where 1) is the parameter for early stopping. Lemma C.9 (Chen et al., 2024, Section B.3)). EnE Buw
BB (w)] S 63.
Thus we have the following uniform bound for our initialization. Corollary C.10. With the same assumption in Lemma C.7, we have
| (w) ow) Sd.
sup Bun O rE(0.
sup
C.3. Convergence of Picard iteration Similarly, define
E/ = sup Eum (w) - P.. (w]]
rE(O,ha.]
and
A Euwniee
11 (w) Pirtaa (w)
Pi | 00 Un-1,ta. 17) We note that by Corollary C.10, El 2 d.
Furthermore, we let E, sup
sup TE 0.
O N
Lemma C.11 (One-step decompasition of 8.). Assume 2h hn I 0.01 and 2hn S 2. For any j ,,
N..,NN- n n= 0, O....., 1, we have
0.
Furthermore, for j = 1, n= I,..., N - 1, we have
Euripe 00 Pla (w) Pla - (1o))((
E < 2 + 0.01
sup rEj0,ka]
29

---

Paralicl Simulation for Log-concave Log-concave-ampllnggg and Score-based Diffusion Models
Proof. For each w E Ω conditioned on the filtration Ft.,, consider the auxiliary process defined in the previous section,
1 2 tenr (w) + sloto (+)))) dr + (w:,+r((w,
d (w)
and
1 V (w) + s to ( (2)))) dr + +lw,++((w)
dig (w)
2
We have
d (w) Binn
2 Vinnt (W) - - (w)) + o sl (Blhoo si sout) (Dims dr.
Then we can calculate the derivative # [w) Vi. as
dr Fin.(w)) (w) Vi. (w))
12
1 ( bl))(w)) + 0 (ottt (Bimbon 2 (o)) (Dh (Bihio.c)))s))]
2 (8 (w) - : (w)
2
By integrating from T, we have
V (w) - (w))
(l (w) (w) dr'
(6 tontt (D (w) s s ton(t)) (3i (w)) dr' (w) 8 (w)) l| dr' + L'' o6g /
2 V (w) Tz (w))
9 slatom (Di (r))))))))()))) || dr'
2
< 2 (w)) : ())) dt' + L3 / | - (
(ri) (w)(()) dr'
By Theorem A.4, and T..P (ww V -i (w), we have
(w) Vi (w]] S Lz S | - "ro dr' + e² A
By taking expectation, for all r e hn]
BuwpFrnn Di (w) ( 27 A
T 1" Bundli Viin (r)(ww) - (Tr) (w))[( dr'
< 23 T sup Eun Fert (w) 8/ ( -'E(
Thus
sup Euniiiin 1 (w) (w) 1
rE(o,hal
< 2h A 13g a 8gg-11
30

---

Paralicl Simulation for Log-concave Sampling and Score-based Diffusion Models
For consider the following processes,
2 yt tn., (w) + +gn+))(7) ( ((ln.) dr + dw.,,t(ww)
1
d (w)
and
yr (w) + 8 (n (D 1st 1.30 (w)) dr dd +r (w
d (w)
Similarly, we have
Bumplien (w) Vinrr(w]] rE(0,An]
sup
L (w)) p 15th 15 1,30 (w))(()
S 2n A + 12 2m hm
sup re(0.h..)
BuwiLinn
Lemma C.12 (One-step decomposition of AW). Assume L2 hn S 0.01 and and 22 S 2. For any j = ,,
N- n , 1, we have
A s 34 n- + 1.48
Furthermore, for j = 1, n 1, I,...,N N 1, we have
Un.r - V 1,7,, - 19.
2
A E 3 + 0.4 sup rE(0.An)
Euwili
For n = 0, we have A & 0.32A 1- and A $ sup
EunijFi rE(0.ho)
1l (w) - (50.0((w))
Proof. By definition of y TM (w) we have
hn
he
e
y
bn In.a2
Tn.ld
m-
Vn.0
- Fna
+
mmmm +1
2 -1) 80 tia Cor.th o sln TT, ( mn
m
m
m
<
Ph.0
+
Σ mm==0
e
mmmmtta +1
2(e -
M-1 Σ
2
vn,
vn.0
1--
+ 32e2
n.m
M
bn + Tn,m
(
S tn (Fn. ]]1"
m ml=0
2
P.. P1.0
+ 32h rE(O,bm]
sup
L
where the second inequality is implied by that T 1 2r when 1. By taking expectation, and the assumption that 2he hn < 0.1 and 2hn 2, we have
AJ 2
c- Ap
|10 Planta.ss
10oo 3i.D 11"
|
0 2
2Eu-iilFi,
+ 32h 3 L²
sup Eu Fin re(0,Am)
24 + 0.328/-1
31

---

Parallcl Simulation for Log-concave Sampling and Score-based Diffusion Models
Thus
4 s $32, + -.48
In the remaining part, we will bound A By definition, we have
ha 11
(w)
yon. T (w)
m
Immn' 11 2(E -1) 100 tr.m (I (In-litnn lit
On.0
Vn 1.T..
+
e
tt VIn- (Fn-1.
m =0
M-1 Σ
<
|n.o
1.7mmr
2
+
Tim'tt 2(e5nmm -1)
2 tn 19n- (nn1i 8 +ram (I (fn-l
m'==0
Vn.0
In
1.
2
+
32h / r
sup
10 Fn-1,
-n-1, |
E
where the second inequality is implied by that ex - - 2x when I < 1. Thus with L2g hp S 0.01 and 2h S 2, we have
Rn A = BuNDLE
--
Denitnee
Den.Tn.
l']
| 20 vn.o
< 21
In-1,7m
+ +-3212, L3 sup rE(0,Am]
Dn 1n11,7m..4 Fn-1,
< <21 n-
0.32 sup
E
100 Fn.r
yn 1,
=E(U.Am.]
Let LAІ 2n + 0.011 We note that L3 2 E/ Thus for and 2 2,
E = 2 n + 0.018/-11
0.48% < + 0.011
< + +
(33)
We recursively bound L3 as
2
n aaa j
10000-0----1 (n-1+5- b)
L <
(0.01)-21600'
L + Σ (0.01)
160
j-b 0) L.
j 2
Bound for Σ (1000))-110000000-----
j-2 L2g.
Firstly, we bound 1,2 To so, by Lemma 12, we bound A as
A < 331
+ 4E, 33
+ Σ 4.3'8 S Σ 3'', < 3 2
2=0
2=0
and by Lemma C. I, bound E
& 24 0.18, < 33
Furthermore, by Lemma C.12, bound A3
カ-1
43 < 3A2 --1 0.48) < 3" 3
Σ
z'el 3igl_ 5 0.32 3r + 3
28 2m283"n8. n,,
32

---

Parallel Simulation for Log-concave Sampling and Score-based Diffusion Models
Thus
T 24 + 0.011 I 28 33.
Furthermore, by (m))
em
for m 2 n2 > 0, we have
(" a (001)-2100000009-2
1161))-1600--
1
j j-2
a=2
j-2
a + 2
$ (0.01))-2(28 160"n2)e
E1
j j-2
s
160"
Bound for (0.01)* 160 1600 ("- LG. By Lemma C.11, we have
j-b
& 20.01 -1 + 2
1-1
s (0.01) & Σ
(0.01)* (0.01)$24))) 2
Combining the fact that A < 0.332--E,, Ey, we have
Thus
1 22/ + +0011
2 0..20-2--16 E3 0.01 (b - 1) 0.32"-16 E3
77 b- 0.32
1- (m+ 1) (=+++) B. 1,++1)
Furthermore, by Σ
(*+t)) x² 1=0
S [1
1 [1-x)oII+ we have
(1
n (n-1+j- b)
n1 n-1 /)
(0.01) 160
(n-1+3-1 b)
n-1 b)7, 7. b. 0.32
(0.01)* 160--1
b
22 0.87 440" -1
jE,
Combining the above two results, have
3 < (e2 2.01))-((28 1600, 22- -40
If 45N > log NE; for any 0, ....NN N
E x²
(34)
C.3.1. OVERALI ERROR BOUND
By the previous computation, we have
(L (l@m)
< KL ((t..) ) + Eanglann
2 2 16 (T,
< KL (T...) + -3Ew
[A (w) + (w]) + h 8/
33

---

Paraliel Simulation for Log-concave Sampling and Score-based Diffusion Models
Combining Lemma A.5, Corollary C.10, and Eq. (34), have
(L (tatt)
N-1
x KL +3
Σ
(Bw-Bbrn (w) + Bt, (w)] + I. T
n=0
x de ed(TT + -loo ) + 62 +
with parameters J - 45N > Olog N h = 0(1), N Ollog $) T = =Ollog (t) O(d- log 10) M Olds log 3), logn) & T.
34


Reviews:
[{"review_id": "mX3Q6wSadV", "reviewer": "Reviewer_qiak", "summary": "Parallel sampling methods propose to speed up sampling by more efficiently simulating diffusions. Prior work splits up the simulation interval into $\\log(d)$ chunks and performs $\\log(d)$  iterations on each chunk sequentially. This yields an overall complexity of $\\log(d)^2$. This paper proposes to remove the need for the $\\log(d)$ sequential steps by communicating between the chunks diagonally as illustrated in their Figure 1. Two close settings of simulating a diffusion with an approximate score are analyzed in the paper: the overdamped langevin diffusion for sampling from un-normalized densities and the reverse diffusion of an OU process. In both cases, the authors show that the $\\log(d)$ sequential steps are unnecessary if the chunks communicate diagonally.\n\n## update after rebuttal\n\nThe authors have fixed their misstated theorem and have clarified their contribution. I think their diagonal update scheme is interesting and I agree that it does introduce some differences with prior analyses. I raise my score as a consequence.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful questions and appreciation of our work, particularly recognizing the clarity of the writing and **the effectiveness of the diagonal communication scheme for parallel sampling**. Below, we address your questions and comments:\n\n**\"...gradient mapping contraction...missing the strong-log-concavity assumption?\"**\n\nYou're absolutely right—the contraction property used in Lemma B.5 indeed requires the strong log-concavity of the target distribution. We will revise Theorem 4.2 to explicitly state this assumption and clarify its role in ensuring convergence under the diagonal update scheme. Additionally, we will update all related claims and discussions throughout the paper to consistently reflect this assumption.\n\n**\"...by closely extending the work of Anari et al 2024... build on the proofs of Chen et al 2024... (understandable) closeness to prior work...\"**\n\nWhile we adopt several standard tools in parts of our analysis—such as interpolation arguments (e.g., VW19 or Section 4.2 of Che23), and Girsanov-based comparisons in Appendix C.2 for the diffusion model—we emphasize that our work addresses analytical challenges not present in prior literature.\n\nIn particular, our **diagonal update scheme** introduces **changing start points**,  which presents analytical challenges not encountered in prior work. As discussed in Section 3, this required new analysis to control error propagation along the Picard iteration and time direction. As Reviewer KnEh pointed, the recurrence\n$$\nL_n^j = a L_{n-1}^j + b L_n^{j-1} + c \\delta^2\n$$\ncaptures the core of our theoretical analysis. We appreciate this insight and will provide additional intuition for this recurrence, highlighting its importance and distinguishing our analytic contributions from previous work.\n\n[VW19] Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices.\n\n[Che23] Sinho Chewi. Log-concave sampling.\n\n[CRYR24] Accelerating Diffusion Models with Parallel Sampling: Inference at Sub-Linear Time Complexity\n\n**\"The L∞ score approximation... unnecessary notational burden... does not capture any realistic setting...\"**\n\nWe would like to clarify that the $L^\\infty$ approximation follows the same formulation as in [ACV24]. As discussed in their paper, this form of approximation is natural in the context of discrete sampling, where approximate scores can be efficiently computed in $O(1)$ iterations using access to a weighted counting oracle or Laplace transform. Furthermore, in spin glass models, the score function can be approximated in the $L^\\infty$ norm using approximate message passing as shown in [AM22].\n\nMoreover, since the final accumulated error in our method depends on the score approximation error at each queried point with varying weights, a natural refinement is to define a weighted version of Assumption 5.1. Also, another alternative relaxation is to adopt a bias-variance decomposition, as in Assumption 4 of [BCESZ22], though it is still formulated in an $L^\\infty$-style metric.\n\n[ACV24] Fast parallel sampling under isoperimetry. \n\n[BCESZ22] Towards a Theory of Non-Log-Concave Sampling: First-Order Stationarity Guarantees for Langevin Monte Carlo\n\n[AM22] Sampling from the Sherrington-Kirkpatrick Gibbs measure via algorithmic stochastic localization\n\n**\"The Lipschitz constant of the score... affects the stepsize $h$...Could the authors expand on the possible dimension dependence on the uniform Lipschitz assumption and its effects on their bound?\"**\n\nWe acknowledge that the uniform Lipschitz assumption in Problem (b) is indeed strong; however, it is a common simplification adopted in prior works [CRYR24], [CCLL22]. As discussed in [SBDD22], the Lipschitz constant must be large to capture multimodal distributions, reflecting a trade-off between model expressivity and training stability. In particular, in settings with well-separated modes, the score's Lipschitz constant can grow significantly, even becoming unbounded near the zero point [YFZS24].\n\nTo ensure convergence of the Picard iterations under these conditions, the product $L_s^2 e^{h_n} h_n$ must be sufficiently small. This requirement implies that the length of each time slice, $h_n$, should scale as $O(1/L_s^2)$. Consequently, the number of time slices becomes $N = O(L_s^2 \\log d)$, leading to an **overall iteration complexity of $O(L_s^2 \\log d)$**. \n\nFurthermore, our regrouping of time grids provides a flexible framework that could help mitigate Lipschitz singularities in other diffusion models, such as E-TSDM, as studied in [YFZS24]. We will expand on this point in the revised version.\n\n[CRYR24] Accelerating Diffusion Models with Parallel Sampling: Inference at Sub-Linear Time Complexity\n\n[CCLL22] Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions\n\n[YFZS24] Lipschitz Singularities in Diffusion Models\n\n[SBDD22] Can Push-forward Generative Models Fit Multimodal Distributions?", "claims_and_evidence": "The paper claims are the following:\n1. Theorem 4.2: We can sample from smooth distributions verifying the log-sobolev inequality in $\\log(d/\\epsilon^2)$ steps by simulating the overdamped langevin diffusion with the diagonal communication scheme in Figure 1. Evidence: The authors provide proofs for this claim by closely extending the work of Anari et al 2024.\n2. Theorem 5.4: We can simulate the reverse diffusion with a learnt, uniformly Lipschitz score, in $\\log(d/\\epsilon^2)$ steps. Evidence: The authors closely build on the proofs of Chen et al 2024 to show that their diagonal communication scheme obviates the need for $\\log(d)$ sequential steps.", "methods_and_evaluation_criteria": ".", "theoretical_claims": "Issues with Theorem 4.2: It appears that theorem 4.2 is misstated. A crucial ingredient seems to be the gradient mapping contraction which is used in lemma B.5. This only holds for strongly log-concave distributions. If my understanding is correct, Theorem 4.2 does not apply to the broader class of smooth LSI distributions. If my understanding is wrong, I kindly ask the authors to clarify.", "experimental_designs_or_analyses": ".", "relation_to_broader_scientific_literature": "The related work is discussed well. The work consists of very close modifications of prior analyses to incorporate the diagonal communication scheme.", "essential_references_not_discussed": ".", "other_strengths_and_weaknesses": "The paper is very well written. The statements of problems a and b could be made more precise by explicitly stating the simulated diffusion. \nThe only minor weakness I see is the (understandable) closeness to prior work and the lack of discussion on the assumptions:\n- The $L_\\infty$ score approximation for problem a is in my view an unnecessary notational burden and does not capture any realistic setting. Approximate scores are either approximate in weaker metrics or are more common in the problem b setting.\n- The uniform Lipschitz assumption in problem b is very strong. The Lipschitz constant of the score can (or, in cases with well separated mode, always) depends on the dimension and it would be beneficial to discuss it as it affects the stepsize $h$ if my understanding is correct.", "other_comments_or_suggestions": ".", "questions_for_authors": "1. Could the authors clarify if their theorem 4.2 is misstated? Is it missing the strong-log-concavity assumption?\n2. Could the authors expand on the possible dimension dependence on the uniform Lipschitz assumption and its effects on their bound?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "R9bYQzL7Ir", "reviewer": "Reviewer_Q5PH", "summary": "The paper obtains novel rates for the parallel complexity of sampling, both in the gradient oracle (MCMC) setting and the score-based denoising setting. The rates are O(log d/eps), which is sharp (in epsilon). The proof is based on a refined Picard iteration scheme, with careful analysis in order to obtain the rate.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful comments and positive evaluation of our work, particularly your recognition of **the sharp iteration complexity in terms of $\\varepsilon$ and the rigor of our analysis**. We also appreciate your conclusion that the work merits acceptance. Below, we address your questions and comments in detail:\n\n**\"A more thorough discussion of condition number (kappa) dependence in both time and space would be appreciated, with reference to prior work.\"**\n\nThank you for pointing this out. We will revise the paper to include a more explicit discussion of how the condition number $\\kappa$ affects both the step size and the overall iteration complexity. In particular:\n\n- **For iteration complexity:**\nOur method achieves an iteration complexity of $O(\\kappa)$, which matches the state-of-the-art sequential query complexity established by [AC24]. In the parallel setting, we believe it is inherently difficult to break the $O(\\kappa)$ barrier using the Picard method, as the length of each time slice must scale inversely with the smoothness in order to preserve non-expansiveness along the Picard direction.\n\n- **For space complexity:**\nOur method builds on the regrouping of scheduled grids along the time direction, and such schedual was proposed by [VW19]. Since we update all grids simultaneously in the worst case, the space complexity is given by $d N_{\\text{seq}} = d \\kappa^2$, where $N_{\\text{seq}} = \\kappa^2$ is the number of time steps required in the sequential setting.\n\n[AC24] Shifted Composition III: Local Error Framework for KL Divergence, Jason M. Altschuler, Sinho Chewi, 2024\n\n[VW19] Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted Langevin algorithm: isoperimetry suffices, 2019\n\n**\"How difficult is it to handle the underdamped Langevin dynamics? What are the barriers and why can the current analysis not immediately give results?\"\"**\n\nApplying our method to underdamped Langevin dynamics, as in [ACV24], yields $\\log(d/\\varepsilon^2)$ iteration complexity and improved space complexity $d^{3/2}/\\varepsilon^2$. However, due to the failure of the triangle inequality for KL divergence, our current analysis cannot establish optimal iteration complexity for KL convergence—our main goal—so we leave this as future work.\n\n[ACV24] Fast parallel sampling under isoperimetry. Nima Anari, Sinho Chewi, Thuy-Duong Vuong, 2024\n\n**\"Can the Lipschitz score estimation assumption be removed, similar to the work from Benton et al.?\"**\n\nWe believe not, at least for Picard-type methods, as the Lipschitz condition is crucial to ensure convergence along the Picard direction.", "claims_and_evidence": "The authors claim a parallel complexity of O(log d/varepsilon) for MCMC under standard assumptions, and for score-based models (SDE variant) under Lipschitz score (which is one of the standard settings). The authors provide detailed proofs for their claims.", "methods_and_evaluation_criteria": "This is not applicable to the paper.", "theoretical_claims": "As discussed, the theoretical claims relate to the parallel complexity, which ultimately arises from a detailed analysis of a Picard iteration scheme. The results seem rigorous on a brief skim.", "experimental_designs_or_analyses": "The paper is not empirical in nature, and so does not contain experiments.", "relation_to_broader_scientific_literature": "The primary references are discussed; this result is related to the parallel complexity of MCMC samplers studied in Amari et al. and subsequent works. Whereas those works obtained $\\log^2(d/\\varepsilon)$ rates, the present rate obtains the same with a single logarithm.", "essential_references_not_discussed": "As far as I am aware, the most relevant references have been covered.", "other_strengths_and_weaknesses": "The primary weakness of this result is that the analysis feels somewhat like an incremental result, utilizing the same Picard scheme (although improved) compared to prior work.\n\nNonetheless, the improvement in logarithmic factor (which is highly important in the parallel setting) guarantees a tight rate, and the work will be useful for future follow-ups concerning the underdamped Langevin or ODE variant of the score-based generative model.\n\nThus, I think this work has merit and should be accepted as it stands. Nonetheless, if the authors are capable it would be helpful to supplement this work with any additional results.", "other_comments_or_suggestions": "Line 129 is unfinished.\nLine 425: parallalizable -> parallelizable", "questions_for_authors": "A more thorough discussion of condition number (kappa) dependence in both time and space would be appreciated, with reference to prior work.\n\nHow difficult is it to handle the underdamped Langevin dynamics? What are the barriers and why can the current analysis not immediately give results?\n\nCan the Lipschitz score estimation assumption be removed, similar to the work from Benton et al.?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lri3Wnma3r", "reviewer": "Reviewer_KnEh", "summary": "The proposed paper introduces a novel parallel sampling technique that significantly enhances the time complexity of both sampling under isoperimetry and score-based diffusion model inference problems from $O(\\log^2 d)$ to $O(\\log d)$. The primary algorithmic innovation lies in the parallel sampling across time slices, rather than sequentially updating each time slice. Furthermore, the authors provide a corresponding lower bound, demonstrating the optimality of their proposed algorithm.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed and positive feedback. We greatly appreciate your recognition of **our novel parallel sampling technique, the clarity of the paper’s structure, the soundness of the theoretical analysis, and the potential of our approach to substantially improve time complexity and contribute meaningfully to the field**. Below, we address your specific suggestions and questions:\n\n\n\n**\"...provide additional intuitive explanations and justifications for the validity of the claim  $L_n^j\\leq aL_{n−1}^j+bL_{n}^{j−1}$...selection of $J$, $N$, and $P$...\"**\n\nWe appreciate your thoughtful observation that this recurrence forms the core of our analysis. Recall we denote the truncation error at the $n$-th time slice and the $j$-th iteration as $L_n^j$ in the squared sense.\n\n- **For Problem (b):** by Young’s inequality and the definition of the exponential scheme, we can bound the change of $L_n^j$ in one update as\n$$\nL_n^j = a L_{n-1}^j + b L_n^{j-1}\n$$\nwith $a = O(e^{h/2})$ and $b = O(L^2 e^{h/2} h^2)$. Thus, by choosing the time step $h$ sufficiently small relative to the Lipschitz constant $L$, we ensure convergence along the Picard direction (indexed by $j$). Specifically, this choice guarantees that $b < 1$ and $a > 1$, with $a$ remaining bounded by a constant.\n\n- **For Problem (a):** since the score itself is not Lipschitz smooth, an additional score estimation error term $c \\delta^2$ appears in the one-step update of $L_n^j$. To ensure the total score estimation error remains bounded, it is necessary to have $a, b < 1$. Similarly, by Young’s inequality and the definition of the Euler–Maruyama scheme, we can bound the change of $L_n^j$ as\n$$\nL_n^j = a L_{n-1}^j + b L_n^{j-1} + c \\delta^2\n$$with$$\na = 1 - 0.1 \\frac{\\beta h}{\\kappa} + O(\\kappa)(3 \\beta^2 h^2)^P,$$$$\nb = O((\\beta^2 h^2)^P), \\quad\nc = O(\\kappa \\delta^2 h^2)\n$$Here, intuitively, $a$ comes from the contraction of the gradient mapping with an additional term from the Picard direction, $b$ reflects convergence along the Picard direction, and $c$ accounts for the accumulation of score estimation error $\\delta$ over time length $h$, with an additional scaling by $\\kappa$ due to Young’s inequality.\n\n- **On the choice of $J$, $N$, and $P$**:\nTo ensure convergence along the Picard direction, we choose the step size as $h \\approx 1/\\beta$. The total time horizon is $\\alpha \\log(d/\\varepsilon^2)$ for Problem (a), and $\\log(d/\\varepsilon^2)$ for Problem (b), leading to a total of $N = O(\\log(d/\\varepsilon^2))$ time slices in both cases.\nSince Picard iterations converge exponentially fast, it is sufficient for each time slice to undergo $O(\\log(d/\\varepsilon^2))$ updates to ensure the desired overall accuracy. This implies that the total number of diagonal iterations satisfies $J = N + O(\\log(d/\\varepsilon^2))$.\nFor the parameter $P$, which controls the number of steps within each blockwise Picard update, we require $a < 1$ in the recurrence relation to ensure contraction in Problem (a). This imposes the condition $P = \\Theta(\\log \\kappa).$ \n\n**\"Can the space complexity be improved from $O(d^2)$ to $O(d^{3/2})$?\"**\n\nYes—if only total variation (TV) convergence is required, applying our method to underdamped Langevin dynamics or ODE-based implementations of diffusion models is sufficient to achieve a space complexity of $O(d^{3/2})$. However, extending this to obtain similar guarantees under KL divergence remains an open question, which we leave as an important direction for future work.\n\n\n**\"Is it possible to extend this parallel sampling technique to discrete diffusion models?\"**\n\nThis is a great question. At a high level, our method operates by regrouping discretized grids along the time horizon and updating them in a diagonal fashion. We believe this diagonal scheduling approach offers a general framework for enabling parallelism along the time direction. While adapting the method to discrete diffusion models may require specific modifications, we believe the core strategy remains applicable and represents a promising direction for future work.", "claims_and_evidence": "The assertions regarding the algorithmic and technical novelty, as well as the enhancement in time complexity, are substantiated, well-defined, and compelling.", "methods_and_evaluation_criteria": "The proposed parallel sampling algorithm appears to be logically sound and presents a promising avenue for enhancing the time complexity, thereby making a significant contribution to the field.", "theoretical_claims": "I have reviewed the appendices for the proofs, and they appear to be sound.", "experimental_designs_or_analyses": "This is theoretical work and no experiment is involved.", "relation_to_broader_scientific_literature": "The work is well-positioned in the literature. The $O(\\log d)$ time complexity is clearly an improvement from $O(\\log^2 d)$ in [1] in the context of diffusion models and in [2] in the context of sampling.\n\n[1] Chen, Haoxuan, et al. \"Accelerating diffusion models with parallel sampling: Inference at sub-linear time complexity.\" Advances in Neural Information Processing Systems 37 (2024): 133661-133709.\n[2] Anari, Nima, Sinho Chewi, and Thuy-Duong Vuong. \"Fast parallel sampling under isoperimetry.\" The Thirty Seventh Annual Conference on Learning Theory. PMLR, 2024.", "essential_references_not_discussed": "The literature review is comprehensive.", "other_strengths_and_weaknesses": "The paper is well-structured, effectively summarizing the related works, proposing the novel algorithm, and presenting results in a mathematically rigorous manner. The authors have also included numerous explanatory remarks throughout the paper to facilitate comprehension.", "other_comments_or_suggestions": "It would be more beneficial to provide additional intuitive explanations and justifications for the validity of the claim $L_n^j \\leq a L_{n-1}^j + b L_n^{j-1}$. It appears that only some intuition for Problem (a) is presented in the main content, which still lacks certain details (“make use of the contraction of gradient descent”). Consequently, readers may require a thorough exploration of the proofs to gain a comprehensive understanding of the algorithm. In my opinion, this formula is the essence of the algorithm and even the entire paper. It would be preferable to explicitly express $a$, $b$, and (in the presence of $\\delta^2$) $c$ in terms of the length of the time slice, Lipschitz constant, and other relevant parameters. Furthermore, it would be beneficial to provide additional discussions on the selection of $J$, $N$, and $P$ (which is not adequately explained in the main content).", "questions_for_authors": "- Can the space complexity be improved from $O(d^2)$ to $O(d^{3/2})$?\n- Is it possible to extend this parallel sampling technique to discrete diffusion models?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QapunANkjN", "reviewer": "Reviewer_aCjg", "summary": "The authors of this paper study two separate but related problems: sampling from an isoperimetric distribution and also sampling from a score-based diffusion model.  Under\nsome assumptions on the target distribution $\\pi$---e.g. it satisfies a\nlog-Sobolev inequality (isoperimetric) and is $\\beta$-log smooth---the\nauthors demonstrate an approach to generate samples arbitrarily close to being\nsampled from the isoperimetric target or a score-based diffusion model that only takes\n$O(\\log(d/\\epsilon^2))$ iterations, while the current state of the art takes\n$O(\\log^2(d/\\epsilon^2))$ iterations. This result assumes that the number of\ncores available to parallelize computation is $\\Theta(d/\\epsilon^2)$, but\nthey provide analogous results in the case they are bottlenecked on the\nnumber of cores. Their results also typically require more space than some\nalternatives, which they attribute to the more volatile paths of the\noverdamped Langevin diffusion vis-a-vis the underdamped Langevin\ndiffusion. E.g., alternatives usually require\n$\\tilde{O}(d^{1.5}/\\epsilon^2)$ space while this approach requires\n$\\tilde{O}(d^2/\\epsilon^2)$.\n\nThe paper is purely theoretical in nature, and to prove their above result,\nthey build upon previous work that uses Picard iterations to subdivide the\nsampling problem across time. However, rather treating each subdivided\nsegment as independent and parallelizing across the different stages, the\nauthors propose a method that iterates over segments based on both the previous time slice and previous stage of the current time slice to propagate their approximations.\n\n## update after rebuttal\nBased on the authors feedback and other reviews, I'm inclined to keep my score of a weak accept mostly out of my lack of confidence.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their detailed and constructive feedback. We are encouraged by your recognition of the core contribution—**the diagonal reordering of Picard iterations to achieve improved parallel iteration complexity, the soundness of the theoretical analysis, \nand the potential applicability of our approach in high-compute environments (e.g., with sufficient RAM and cores)**. We address the specific concerns and questions raised below:\n\n**\"Are there any ways to empirically validate whether the approaches outlined here hold? If so, could they be added to the paper and compared against alternatives?\"**\n\nWe agree that empirical validation is important. While we are not aware of existing implementations of our diagonal update scheme, we believe it could be adapted to diffusion-based samplers by modifying their update scheduling. We will add a discussion on this potential direction in the revised version.\n\n**Regarding typos:**\nWe will carefully proofread the manuscript and correct all grammatical and typographical errors to improve the clarity and overall presentation in the revised version.", "claims_and_evidence": "There are no empirical studies to back up their claim. That being said, the\nauthors provide a sketch of the proof in Section 3 and also a\ndetailed proof of their main claims in the appendix. E.g., Theorem 4.2, which establishes\nthe superior rate (in terms of iterations) for approximately sampling from an isoperimetric\ndistribution, is proved in Appendix B, while Theorem 5.4, which establishes a similar result but for approximately sampling from a score-based diffusion model, is tackled in Appendix C.", "methods_and_evaluation_criteria": "There are no empirical analyses conducted in this paper.", "theoretical_claims": "No, I only skimmed the proofs in the appendix, these were not checked\nline-by-line.", "experimental_designs_or_analyses": "Again, there are no empirical results in this paper.", "relation_to_broader_scientific_literature": "The Algorithm 1 presented in the paper could provide a way to speed up the\nalgorithms that sample from diffusion-based methods when the number of\ncompute resources (e.g. RAM and cores) is high. If demonstrated empirically,\nthis could have a substantial effect on the genAI industry.", "essential_references_not_discussed": "To the best of my knowledge, no.", "other_strengths_and_weaknesses": "The paper's core contribution seems to be reordering the manner in which the\nPicard iterations are taken to approximate sampling from the path of the\nSDE. This appears to be an interesting idea, although the significance is a\nbit tempered by the extra space complexity needed to execute this approach.\n\nThe paper is only theoretical in nature, so it is hard to state its\npractical significance. It's case would be greatly bolstered if there were empirical\nwork demonstrating its efficacy. That being said, it appears to produce results on\nthe Pareto frontier of time and space complexity for sampling from an SDE.\n\nThe paper is fairly well written, although there are a few mistakes here and\nthere.", "other_comments_or_suggestions": "L128: capitalized \"The\"\n\nL144: Might have missed it, but what's V?\n\nL183: Is the integral correct? Seems like the $f_t$ should be $f_s$?\n\nL194: \"girds\"", "questions_for_authors": "[Q1] Are there any ways to empirically validate whether the approaches\noutlined here hold? If so, could they be added to the paper and compared\nagainst alternatives?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper presents a novel parallel sampling algorithm that reduces the iteration complexity from O(\log^2 d) to O(\log d) for sampling from isoperimetric and diffusion-based models. The key innovation is the diagonal reordering of Picard iterations, which removes the need for O(\log d) sequential steps.

All reviewers agree that the paper makes a strong theoretical contribution to the field of parallel sampling. The reduction in complexity is both significant and impactful. While an empirical evaluation would strengthen the paper, the novelty and rigor of the work are clear.
Decision: Accept (spotlight poster)