Polynomial Convergence of Riemannian Diffusion Models
Authors: Xingyu Xu
Ziyi Zhang
Yorie Nakahira
Guannan Qu
Yuejie Chi
Venue: ICLR 2026 Poster (2026)
Keywords: diffusion model
Riemannian manifold
gradient estimate
Primary area: learning theory
TL;DR: We prove that discretized Riemannian diffusion model converges in total variation distance within polynomially many steps.

Abstract: Diffusion generative models have demonstrated remarkable empirical success in the recent years and are now considered the state-of-the-art generative models in modern AI. These models consist of a forward process, which gradually diffuses the data distribution to a noise distribution spanning the whole space, and a backward process, which inverts this transformation to recover the data distribution from noise. Most of the existing literature assumes that the underlying space is Euclidean. However, in many practical applications, the data are constrained to lie on a submanifold of Euclidean space. Addressing this setting, de Bortoli et al. (2022) introduced Riemannian diffusion models and proved that using an exponentially small step size yields small sampling error in Wasserstein distance, provided the data distribution is smooth and strictly positive.
In this paper, we prove that a *polynomially small stepsize* suffices to guarantee small sampling error in total variation distance, without any assumption on the smoothness or positivity of the data distribution. Our analysis only requires mild and standard curvature assumptions on the underlying manifold. The main ingredients in our proof are Li-Yau estimate for log-gradient of heat kernel, and Minakshisundaram-Pleijel parametrix expansion for perturbed heat equation. Our approach opens the door to a sharper analysis of diffusion models on non-Euclidean spaces.

Published as a conference paper at ICLR 2026
POLYNOMIAL CONVERGENCE OF RIEMANNIAN DIF- FUSION MODELS
Xingyu Xu, Ziyi Zhang, Yorie Nakahira, Guannan Qu, Yuejie Chi Department of Electrical and Computer Engineering
Carnegie Mellon University Pinsburgh, PA 15213, USA (xingyuxu, ziyizhan, ynakahir, gqu, yuejieceeeddeddre. cmu edu
ABSTRACT
Diffusion models have demonstrated remarkable empirical success in recent years and are considered one of the state-of- state-of-the-art generative models in modern AI. These models consist of forward process, which gradually diffuses the data dis- tribution to noise distribution spanning the whole space, and a backward pro- cess, which inverts this transformation recover the data distribution from noise. Most of the existing literature assumes that the underlying space is Euclidean. However, many practical applications, the data are constrained to lie on sub- manifold of Euclidean space. Addressing this setting, De Bortoli et al. (2022) introduced Riemannian diffusion models and proved that using an exponentially small step size yields small sampling error in the Wasserstein distance, pro- vided the data distribution is smooth and strictly positive, and the score estimate is To accurraee. In this paper, greatly strengthen this theory by establishing that, under L2-accurate score estimate, a polynomially small stepsize suffices to guarantee small sampling error in the total variation distance, without requiring smoothness or positivity of the data distribution. Our analysis only requires mild and standard curvature assumptions on the underlying manifold. The main in- gredients in our analysis are Li-Yau estimate for the log-gradient of heat kernel, and Minakshisundaram -aejee parametrix expansion of the perturbed heat equa- tion. Our approach opens the door to sharper analysis of diffusion models on -un-Euclddenn spaces.
INTRODUCTION
Initially introduced by Sohl-Dickstein et al. (2015) and later advanced by Song & Ermon (2019); Ho et al. (2020); Dhariwal & Nichol (2021), diffusion model has become one of the bedrocks in generative modeling across variety of application domains such vision, video, speech, and many others. On a high level, diffusion models generate samples from a target distribution by operating on two stochastic processes:
1. A forward process
add noise
add noise
add noise Xx.
Xp
X1
where X, is sampled from the target distribution po in Tid and XT resembles pure noise.
2. A reverse process
Yy
denoise
Yy-1
denoise
demoise
Yo:
where Yr starts from pure noise, and gradually removes the noise, so that Yo, we recover new sample from a distribution close to Po-
The reverse process built to recover the target data distribution by step-wise reversing the forward process, with goal of matching the probabilities Y, ~ X, in distribution for t {T..... 1}. Lever- aging the theory of backward stochastic differential equations (SDE) (Anderson, 1982; Haussmann & Pardoux, 1986), this be formally achieved soon the long the score function, i.e.,

---

Published as a conference paper at ICLR 2026
the log- gre-gradient of the the marginal density of the forward process, becomes available, which can be estimated via score matching (Hyvärinen & Dayan, 2005).
Tremendous recent progresses have been made in understanding the convergence of diffusion mod- els in the Euclidean space, e.g. Lee et al. (2023); Chen et al. (2023); Benton et al. (2024); Li et al. (2024), which establish near-tight polynomial iteration complexities of discrete-time samplers un- der L2-accurate score estimates and mild assumptions of the data distribution. These convergence results provide strong justifications the empirical success of diffusion models for generating from complex multi-modal distributions.
Many scientific domains, however, are intrinsically non-Fuclidean; examples include orientations on (3)3), directions on spheres, toroidal angles, articulated poses, and symmetric positive definite (SPD) matrices are naturally modeled on Riemannian manifolds (Piggott & Solo, 2016: Muniz al., 2022). Recently, there has been an increasing interest in effectively sampling from distributions sup- ported on manifolds and providing theoretical guarantees (Girolami Calderhead, 2011: I: Gatmiry & Vempala, 2022; Li & Erdogdu, 2023; Guan et al., 2025). Although sampling on manifolds has been studied extensively (Cheng et al., 2023), extending diffusion models to manifolds requires careful treatments to incorporate the manifold constraints into both the time-inhomogeneous forward and reverse processes, with selected attempts in De Bortoli et al. (2022); Huang et al. (2022); Lou et al. (2023); Liu et al. (2023); Fishman et al. (2023).
One notable development is De Bortoli et al. (2022), who introduced Riemannian Score Score-Based Gen- erative Models (RSGMs) with convergence guarantees in the Wasserstein distance. Specifically, they established a time-reversal diffusion process for geometric Brownian motion on manifolds, which can be similarly learned via score matching (Hyvärinen & Dayan, 2005). While groundbreaking. their convergence bound suffers from a few caveats: 1) it requires an exponentially small stepsize, leading to possibly exponential iteration complexity in some of the manifold parameters; 2) requires L -accuraee score estimates, which are impractical in deep learning; and 3) the data distri- bution required to be smooth and strictly positive on compact manifolds. This naturally raises the following questions:
Can we achieve polynomial iteration complexity for manifold diffusion models using Læ-accurate score estimates under mild data assumptions?
1.1 OUR CONTRIBUTION
We provide discrete-time analysis of the RSGM sampler in De Bortoli et al. (2022), assuming Læ-accurate score estimates. Under mild geometric conditions of the manifold without assuming smooth or strictly positive data densities, we establish that polynomial stepsizes suffice for accurate sampling on manifolds total variation (TV). More precisely, for some & > 0, the TV error between the distribution of the output Yo and ps. that is, an approximation to po with early-stopping time > 0, obeys
TV Law(Yo)) S + Eccore + VET poly(d,8-1)
as long as the horizon satisfies T R (d + log(d/E)). Here, d is the dimension of M, A is the spectral gap of -AM, h > is the stepsize, and Escore is the L2 score estimation error.' This bound suggests that a polynomial stepsize is sufficient: take T x (d d + log(d/e)) and h 2² -yT then the TV error is bounded by E Escore after iteration complexity of
N T/h = poly 5
(d, (AAee))
This conveys much more benign message about the efficiency of Riemannian diffusion models, compared with the iteration complexity in De Bortoli et al. (2022) that scales exponentially with the dimension d, under relaxed assumptions both the data distribution and the score estimates.
Techniques. Our proof highlights three ingredients: (i) high-probabilltty Yau gradient bounds for the manifold heat kernel together with early stopping to control || V log Prll without assuming positivity/smoothness of po: (ii) a localization scheme that "freezes" drifts across nearby tangent
For simplicity, we omitted polynomial dependence manifold geometry parameters and poly-log depen- dence The complete version be found in Theorem
2

---

Published as a conference paper at ICLR 2026
Work
Structure
Metric
Iteration complexity
Data distribution
Benton et al. (2024)
Euclidean
TV
Ō(d/E²)
bounded moment bounded support bounded moment smooth, strictly positive None (early stopping)
Li et al. (2024)
Euclidean
TV
((poly(d)/))
Li & Yan (2025)
Buclidean
TV
O(d/e) (exp(O(d)//e-/a)
De Bortoli et al. (2022)
Manifold
W P
Ó( (Polyid))
This work
Manifold TV
XPE"
"The original Wp error in De Bortoli et al. (2022) is stated in the form of Ô(Ce A + vh). where C' is defined in Proposition C.6 therein, which in turn is specified by Urakawa (2006, Proposition 2.6) as the supremum of ( 26, y). where H is the heat kernel on the manifold. In general, the best estimate for this is due to Yau (Li & Yau, 1986), which gives C' I col)) achieve e-error, we must set Ao (St(d) + loge and h = e -2 g², then the iteration complexity T/h has the claimed form.
Table I: Comparison of the current theoretical guarantees diffusion probabilistic models Eu- clidean spaces and manifolds. Here, A > 0 is the spectral gap of the ssplace-Beltrant operator.
spaces but preserves continuous Brownian motion (BM), separate the effects of discretizing scores and BM; and (iii) quantitative estimates for Minakshisunddram-PPaae one-step deviations between the manifold heat flow and its discretized proxy. These components allow us to handle the discretization errors sharply to avoid exponential dependence.
parametrix that controls
1.2 RELATED WORKS
Non-asymptotic convergence for Euclidean diffusion models. Early convergence analyses of diffusion models require Los -cccurate score estimates (De Bortoli et al., 2021). For stochastic sam- plers such as DDPM (Ho al., 2020), early bounds under Lipschiz/smmotthnes assumptions of the data distribution admit an O iteration complexity in the total variation distance assuming L2-accurate score estimates (Chen et al., 2023), with subsequent analyses relaxing the Lipschitz assumption yet retaining the same complexity (Lee et al., 2023; Benton et al., 2024; Li al., 2024). More recently, Li & Yan (2025) has improved the iteration complexity to Ö For determinis- tic samplers, Chen et al. (2023) established polynomial convergence with exact scores, and Li et al. (2024) established a convergence rate of OT under L2-accurate scores. See Beyler & Bach (2025); Liang et al. (2024); Li & Jiao (2024) for additional analyses that established convergence in the Wasserstein distance and improved discrete-time rates. Several works (Li Yan, 2024; Liang al., 2025; Huang al., 2024; Potaptchik et al., 2024) also developed non-asymptotic convergence rates of diffusion models under the manifold hypothesis, suggesting diffusion models are adaptive to low-dimensional structures. This line of work should not be confused with ours, where the diffusion process is designed specifically to be constrained on the manifold.
Sampling on Riemannian manifold. Cheng al. (2022; 2023) analyzed the geometric Eu- ler Maruyama (EM) discretization for homogeoomge SDEs, and proved a polynomial com- plexity guarantee under dissipative-distann geometric assumptions the manifold. See also Bharath et al. (2025) for follow-ups. Guan et al. (2025) proposed a Riemannian proximal sampler with con- vergence guarantees under the log-Sobolev inequality. Various sampling algorithms also studied for a related problem known as sampling from constrained spaces (Srinivasan et al., 2024; Ahn & Chewi, 2021). Nonetheless, convergence analyses of Riemannian diffusion models under general data distributions remain highly limited, with De Bortoli et al. (2022) being the only prior work with non-asymptotic convergence rates.
2 BACKGROUNDS
2.1 DIFFUSION MODELS ON EUCLIDEAN SPACE
We briefly recall diffusion processes Rd Let (W.) be standard Brownian motion in Rdd.
3

---

Published as a conference paper at ICLR 2026
Forward SDE and Fokker-Planck. Given drift term b, R R2, the forward process (Xt)kejo,T) solves the Itō SDE
dX, (x,)) dt + dW. Xo Po
Let Pt denote the law of X,, then Pt satisfies the Fokker-Planck equation d,p, (bbpp)) Ap.. In the driftless setting where b, = 0, the marginal Po * (p: is Gaussian smoothing of pu with kernel (z) (2xl)= dd/2
exp ( /22))
Score and reverse process. The score of the forward process at time t is defined as sx(2) = V log Pr (z)). The time-reversal identity (Anderson, 1982) yields reverse-time process (Y)) whose marginals match those of (Xx)1e(0,7] which solves the reverse-time SDE (note that t flows from to O):
dy, = (-) s, + dW,.
T t, Yr ~ pxr
Discretization with approximate score. Let 0= to < t < s.. .....t =tx=11 T, where & - 4--1 =: be a time grid. In practice, exact score function is often unavailable. Instead, we use an approxima- tion & trained via score matching (Hyvärinen & Dayan, 2005). The Euler-Maruyama discretization of the reverse-time SDE above is given by
yk yk l [- be, (yk) Szn (uk)]
Viga.
gk NO, Ta).
In our driftless setting, this reduces
Yk-1 yv h (yk) + Vhin. 9n ~ ~(0,
2.2 GEOMETRY AND NOTATION
We assume some familiarity with Riemannian geometry. and make use of standard notation. Please refer to Jost (2017); Petersen (2006) for more in-depth treatment. In particular, use o; E, etc., to index coordinate representation of tensors, and assume Einstein's 's summation convention. Let (M,g) be a connected, compact d-dimensional Riemannian manifold, with geodesic distance p(:s)) and volume measure μ.. We assume μ(M) = 1. The Levi-Civita connection is denoted by V, and the Laplace-Belle Beltrami operator by
AMI VVV
We use T.M for the tangent space at x and use exp, T.M- - M for the exponential map and log, for its local inverse the normal neighborhood of Furthermore, define geodesic diameter of (M,9)
Diam(M) = sup p(z,y); a.JE.
where -(: :) is the geodesic distance induced by 9. We further denote Rm as the Riemannian curvature tensor. Geodesic ball centered at x with radius is denoted (r)..
We use the total variation (TV) and the Kullback-Leibe (KL) distance to measure the discrepancy between two distributions
/ JM dp dq), KL(P 9) () / (log dp)
TV(p.9)
dp.
dq
2.3
HEAT FLOW, BROWNIAN MOTION, AND DIFFUSION ON M
We also recall the setup for SDE and diffusion processes on Riemannian manifolds introduced in De Bortoli et al. (2022); Cheng et (2023). Let (W,)ezo be standard Brownian motion in Rd and U, Rdd T T. any orthonormal frame x. The Geometric Brownian motion solves dX, Ux, dW,, where o denotes Stratonovich integral, and its transition density (T 3) with respect to μ solves the heat equation
(-.) = Ampe(s.).
4

---

Published as a conference paper at ICLR 2026
Algorithm 1 Riemannian Score-Based Generative Models (RSGM)
1: Manifold (M.g), g): score (z();; early stopping time 5 > O: reverse time grid ó to t
< ty = T; step size = tg t--1; initial ØN ~ μ (uniform distribution);
2: for k E
...., (N, 0, do
3: Choose an orthonormal frame Us at RA, which is linear map from Rd to M. Gk + UnEa E TyMM.
4:
Ek ~ NO,, LT) in Rd,
5:
bk (Yk) E MY
6:
AK f hbp + VhGR Tr
7:
if | h1/4 then
8:
Y + expy, (Ax)
9:
else
10:
Yk-1~
11:
return Yo
Equivalently, Brownian motion can be defined abstractly the solution to the martingale problem for the operator L.M. Concretely, for any E C ([0, oc) x M), the process
Mf = f(1, S(0, Xo) - - (a.
0, 0.++
(A [(, X,)ds
2
is martingale with respect to the natural filtration of X. More generally, forward diffusion process with drift is given by
dX, bt( dt + UX, oW
with Lakker-Plancee equation d.P: = - V(bipt) + AAMPt- Note that in this setting, the following process is martingale for smooth f:
Mf = f(t, f(0, Xo)
(be, Vf) ((... (br,
(Amf) (s, X,)ds.
Let Pt denote the density of X, w.r.t μr, and define the score st identity on manifolds yields a reverse SDE:
V log The time-reversal
xx, (X,, + V lo (X-))dtt + Ux, ddW, == t,
Xx pxr.
In practice, the score log Pt approximated by a trained neural network (z()..
Last but not least, note that compact manifolds, -AM admits spectral gap A > 0. Any initial distribution mixes to the uniform distribution f along the heat flow with rate -
MAIN RESULT
In this section, for completeness, we first introduce the RSGM algorithm in De Bortoli al. (2022).
Then, we offer our polynomial convergence guarantee in Theorem 1. For simplicity, we use driftless forward process:
dx, = Ux o dW, The reversallrsal identity yields the reverse-time SDE
Xo Po-
dY, Vog (Y;) dt + Uy o dW,, T t, YT Pr-
(2)
In Algorithm 1, we provide an outline discretized reverse-time SDE on Riemannian manifold, modified from De Bortoli al. (2022). In each reverse step &o {N,..., 0,0.. thonormal frame UR at then sample Gaussian noise &L and lift it to the tangent space Ty, using the orthonormal frame, obtaining G3 E MM Afterwards. we propose tangent update A = h&,,, (yx) + vhGx and the project to the manifold using the exponential map. To prevent the update from exitting the injective radius, we perform rejection sampling step that rejects exceed- ingly large update. The algorithm terminates at & = and returns the final iterate you In this way, ensure every update is well-defined in normal coordinates during the algorithm.
we select an or-
Before presenting the main theorem, we formalize the assumptions needed for the convergence guarantee.
5

---

Published as a conference paper at ICLR 2026
Assumption (Regularity). Let (M.g) be connected, compact d-dimensional Riemannian man- ifold. We assume the following conditions on M:
(A1) Positive injectivity radius: there exists some constant K 2 such that the injective radius 1/K.
(A2) Uniform curvature bounds: for the same constant K (which can be enlarged if necessary), we huve
max Diami(M), ||Rm]] [RRaal/.,,////Iliaaaaa 1V
K.
(A3) Regularity of score estimates: there exists polynomial poly K), such that
(x))l S (d, K) (IIV Pog (2)]l + ())), Vo M.
In Assumption I, we made the standard "bounded geometry" assumption; similar assumptions also occur in Cheng al. (2022); De Bortoli al. (2022). A positive injective radius ensures that have sufficient room to operate on the tangent spaces as a proxy of operating on manifolds. since for every x M, the exponential map exp,, diffeomorphism on the geodesic ball within injective radius. Bounds on Riemannian tensors rule out pathological cases, which helps to control the error propagation along the reverse diffusion. Lastly, compactness ensures positive spectral gap AM with À1 O, which is necessary to guarantec that the forward process mixes. The mild assumption (A3) on the score estimates avoids excessively large drifts in diffusion, and can be implemented easily in practice by clipping. In addition the above, also need standard assumption the score estimation error (Chen et al., 2023).
Assumption (Score estimation error). There exists Escore > 0 such that
(tx to- 13,, (Yr.) V llogtt, (112 I scoree
With the above assumptions, we are now ready to present our main convergence guarantee for RSGM, outlined in the following TV-accuracy bound.
Theorem 1. Assume Assumptions I and hold, There exists some universal constant C, C' > O (dlog( (dd)+ +K+ogl NN )), then the output Yo of Algorithm
such that the following holds. IfT 2 A1 obeys
TV(ps, Law (Yo)) < E S'core VhT poly(d, (, K,8-)),
where h is the discretization step size, A > 0 the mixing rate of the geometric Brownian Motion on M, i.e., the smallest eigenvalue of - AA in (2))
A few remarks are in order.
The error bound decomposes cleanly into three terms: E results from mixing of the heat semigroup at the spectral gap λ1 Escore captures error from imperfect score estimation, and VhT poly (d, K, 5 5-1) is the discretization error controlled by the step size and curvature. Con- sequently, choosing T x >T (dlog d + log(d/=)) and h poly(d,K,5 Ty7, then the TV error is bounded by Escore after polynomially many iterations
Iteration complexity.
N T/h poly(d, poly(d,K,ō-1)
(A15)2
Compared to prior convergence rates the Wasserstein metric (De Bortoli et al., 2022), which re- quire exponential complexity, we achieve polynomial convergence of Riemannian diffusion models for the first time. Nonetheless, we emphasize that TV and Wasserstein distances are incomparable with each other in general, and our guarantee complements prior Wasserstein results (De Bortoli et al., 2022) by ensuring distributional closeness different notion with much smaller number of iterations.

---

Published as a conference paper at ICLR 2026
Possible improvements. We note that the bound established in Theorem holds under very mild geometric assumptions, requiring only constraints on the injective radius and Riemannian curvature. The purpose of this study is demonstrate that, in the manifold setting, the exponential blow-up in 7 can be avoided and polynomial complexity can be achieved. To keep the exposition as simple possible and to clearly highlight the key ideas, have not attempted optimize the current bound on the degree of the polynomial. Potential approaches for sharper bounds include: (i) better design of discretization schedule, possibly adaptive the manifold geoemetry, and more careful computation of discretization error, such as those in Li & Jiao (2024), Benton et al. (2024) (notably, the dependence on 5 might be improved to poly-logarithmic in this way); (ii) a tailored analysis for TV error that does not rely on Pinsker's inequality, like those in Li & Yan (2025), may also be extended to manifolds; (iii) tighter version of our Minakshisundarame Pleijel parametrix bound. We leave these improvements as future work.
4 PROOF OUTLINE
Throughout the proof, we assume that
h 1
poly(d, K, K
since otherwise the bound in Theorem 1 would be trivial (recall that TV distance is always bounded by 2). We start by recalling the sequence considered in RSGM. Let (Y...E( be given by YN ~ a and for any k E NV 13:
expy, hin (Vi) + vhGe drawn from 14:
| (V)) + VhGal $//4
Yk-1-
otherwise.
This defines sequence of probability transition kernels Re For simplicity, we denote this by Kg. Let qk the law of Yk. We have
@o = qN RN R2.
Similarly, the probability transition kernel from time la to lx-1 in (2) is denoted by Keyde or KR in short. We have
Po = PNNKNNN- -.. .I.
Our goal would be to bound TV(po, 2a) in Theorem 1, by decomposing the total emror into four components:
(initialization error) (score error) (drift discretization error) (BM simulation error).
More concretely:
Initialization error arises from initializing YN with μ instead of the true marginal p.;
Score error arises from imperfect score estimation;
Drift discretization error arises from approximating the continuous-time drift 6,(Y)) by its "time-frozen" counterpart St, (Y,,);
Brownian motion (BM) simulation error is a distinctive feature of the manifold set- ting. Unlike Euclidean space where the transition kernel of Brownian motion over [tx, fx-1] is exactly Gaussian with variance (f4-t+-)) - the transition kernel of manifold- valued Brownian motion cannot be simulated exactly by any discrete-time process, even after time discretization. This inherent inexactness gives rise to this final error term.
The first two components are relatively easier to bound using well-established tools: bounds of heat flow (Urakawa, 2006) and Girsanov transform (Chen et al., 2023). For the drift dis- cretization error, recent techniques developed in the Euclidean setting (Benton et al., 2024) can also be adapted with modifications that account for the manifold curvature. However, the last component the Brownian motion simulation error represents the core challenge in the manifold setting, which fundamentally denies direct extension of Fuclidean analysis.
mixing rate
7

---

Published as a conference paper at ICLR 2026
Step I. Constructing auxiliary kernels via localization. In view of this, we first introduce an intermediate random process that separates the drift discretization error from the BM simulation error. Constructing such process, however, involves additional technicality. In particular, the frozen drift St, (Yt,,) is a vector in the tangent space Ty, and is therefore only well-defined at the fixed point Yg- This poses compatibility issue: Brownian motion evolves continuously the manifold, it immediately departs from Ys, rendering the frozen drift ill-defined. Careful geometric considerations are thus required to reconcile the piecewise-constan drift approximation with the intrinsic curvature of the manifold.
In our analysis, this is handled using localization by the construction of an auxiliary sequence of transition kernels Kat These kernels do not appear in the algorithm itself; they serve solely as an analytical tool to facilitate the proof. These kernels expose the behavior of the time-reverse SDE (2) when the estimated score 8, is frozen to be constant vector field in between discretization steps, meanwhile keeping the continuous Brownian motion.
Let [0, x)) - [0, 1] be smooth cutoff function, i.e., is decreasing, 7/ u.11 = and 7/14, = 0. Such function can be chosen such that bny + Iny" + in" < 100. Recall that the injective radius of M is lower bounded by 1/K, and the curvature is upper bounded by K. Define
Cw
4r
(r) = 17
r20 0,
(4)
K Kd
w2
where ew > O is a small universal constant. We have nwl = 1 and nwlw = 0. For >> 0, x, y E M, define the following vector field on M:
St,a(1) (d expx))oR,, y (n v)) 3e (x)) T.M.
is the "constant" velocity field (c(z) (z) in normal coordinates, which rep- resents our idea of freezing the drift term for time period. The d exp,, in the formula is responsible for identifying M. with T,M. On the other hand, the cut-off function Th necessary to keep all our discussions restricted to the injective radius, as to avoid pathologies of cut locus.
Intuitively speaking, (-)(-) (-)
With this mind, we are ready to define Kath as the transition kernel from time tx t-- of the reverse-time SDE
dY, P Nr,Y (Yr) + Uy, o dW, == -t, T [tN-1,
(5)
and in addition,
aux = PN Kaux K ... Kaux Pk
k= N, 1,
Step II. Decomposing different sources of error. We now decompose
TV(po, go) - TV(po, pith) TV(p go) S v 2KL ||
pāxx) + TV(p aux .90),
where the last inequality used Pinsker's inequality. To control KL || rox)), we further introduce the counterpart of Ftrr using the exact score function V log.p:
(v) (deex.))og. (p,x, log)((z)) E T.M.
We apply Girsanov's 's theorem (Hsu, 2002) to compare (5) with (2), in a way that is standard recent literature (Chen et al., 2023; De Bortoli et al., 2022). Denote the path law of the solution of (2) by Law(Y), , and the path law of the solution of (5) by Law(Yaau).. Girsanov's theorem asserts that the KL divergence (L(Law(Y)) Law(Y (Yrau)) upper bounded by the expectation of the squared norm of the difference between the drift terms in the two SDEs. More concretely.
N
logp:(Y)) - Pox. (x3)||= dt.
KL(Law(Y)) || Law aux)) <
Generally speaking. it is more natural to parallel transport identify different tangent spaces. However, this would later lead more complicated treatment the perturhed heat equation with variable drifts. We choose to use parallelism normal coordinates instead for simplicity.
In its classical form, Girsanov's 's theorem requires integrability such Novikov's condition to hold. In our setting, this he hypassed with localization argument as in Chen al. (2023).
8

---

Published as a conference paper at ICLR 2026
Since Po and p are marginals of Law(YY and Law(Yaau) respectively time t to. by post- processing inequality, we have
(P p
E
V Vlogpt(Y))
(v))
N
V log Po (Y) Finy, dt
+2 Σ
B Lu, (Y) Fi (x)
1k-1
fx-1
drift discretization
matching
It remains to decompose TV qo). To isolate the initialization error, we introduce
On PN KN K1.
By triangle inequality and post-processing inequality, we have
TV(p qo) S TV(pb ) + TV qo) S TV(po aux 9o +
TV(px.
BM simulation
initialization
errors. By design, QN ,4, and TV(pN. gN) = TV(px.pL p)) This is known as the mixing rate of heat flow in total variation norm, and has well-established bounds, e.g., Urakawa (2006). The score-matching error, on the other hand, can be controlled with an analysis on the distortion on the Riemannian metric in normal coordinates. We compile the bounds into the following lemma.
Step III. Controlling initialization and
Lemma 1. There exists universal constant C 0, such that whenever T 2 1, we have
TV(pN, QN) cC(K- - (T)
and
E Be. (Yi) Fiyy (Yi) S 22
Step IV. Controlling drift discretization error with Itt/Stratonovich calculus and Li-Yau esti- mates. The drift discretization error defined in (6) has a similar form the discretization error for the Euclidean setting (Benton et al., 2024), though additional complications arise due to non- constant St.y.. The idea study the time derivative of E V log (Y.) - -tiri, (Y:)lI: which in view of θ, log
A APt (negative sign due to reverse time) involves space derivatives of logP; up to third order. Fortunately, after applying Itò/Stratonovich calculus to simplify the ex- pression, key property in the proof of the Euclidean setting carries over: third-order derivatives of log Pt cancel out. The remaining first and second-order derivatives can be controlled by Li- estimates on the log gradient of the heat kernel. We obtain
Lemma 2. Under the assumptions in Theorem I, there is a universal constant C > O such that
Cc K
Σ
B log - Pixy, (X:) <
2 N.
53
Step V. Controlling BM simulation error using parametrix estimates. Our approach is inspired by the following consequence of post-pri post-processing inequality and Pinsker's inequality:
TV 46) S (KL
(6)
KL (1 BuaakuK aux Ka).
This leads to compare the kernel K. and Rg. In normal coordinates, the Fokker-Plance equation shows that these two are the solutions of the heat equations with the Euclidean Laplacian and with the manifold saplace-Beltram operator. We utilize the Minakshisundaram Pleijel parametrix theory (Berline et al., 2003) in geometric analysis for this comparison, and establish a quantitative bound in polynomially small radius and polynomially short time (cf. Lemma 20).

---

Published as a conference paper at ICLR 2026
5 CONCLUSION
We developed discrete-time theory for Riemannian diffusion models showing that a polynomial stepsize suffices for TV-accurate sampling under mild geometric conditions. In particular, our re- sults show that choosing stepsize polynomially small in manifold parameters achieves any pre- scribed TV target without exponential blow-ups in dimension or curvature. This complements prior Wasserstein-type guarantees which require exponentially many steps. Several important future di- rections remain open.
Sharper bounds. For simplicity, we did not attempt to establish sharp bounds for the error terms in our analysis, and is likely that the degree of the polynomial in the bound could be improved significantly by refining our analysis, and some of the polynomial dependencies can be improved logarithmic ones (Benton et al., 2024; & Jiao, 2024).
Analysis of deterministic samplers. We focused on DDPM-style stochastic samplers in our analysis. For practical purpose, it is also tempting to develop an analogous theory for DDIM-style deterministic samplers (Song et al., 2021; et al., 2024).
Conditional sampling. Our theory was for unconditional diffusion models. Applications like solving inverse problems require conditional sampling, which calls for both new algo- rithm design and new theoretical analysis (Xu & Chi, 2024).
ACKNOWLEDGMENT
The work of X. Xu and Y. Chi is supported part by Air Force Office of Scientific Research under FA9550-25-1-0086 1-0060, and by National Science Foundation under ECCS-21266442577778 Z. Zhang, Y. Nakahira, and G. Qu are supported in part by NSF CAREER Award 2339112, NSF Award 2512805, NSF Award 2442948, and the Pennsylvania Infrastructure Technology Alliance. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the United States Air Force.
REFERENCES
Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-Langevvin algorithm. Advances in Neural Information Processing Systems, 34:28405-28411 -24118, 2021.
Brian D.O. Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Ap- plications, 12(3):313-326, 1982. ISSN 0304-4149. doi: https://bb..or//10.0066044444444 90051-5.
Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis. Nearly $dØ-linear convergence bounds for diffusion models via stochastic localization. The Twelfth International Conference on Learning Representations, 2024.
Nicole Berline, Ezra Getzler, and Michele Vergne. Heat kernels and Dirac operators. Springer Science & Business Media, 2003.
Eliot Beyler and Francis Bach. Convergence of deterministic and stochastic diffusion-mont sam- plers: A simple analysis in Wasserstein distance. arXiv preprint arXiv:2508.032 2025.
Karthik Bharath, Alexander Lewis, Akash Sharma, and Michael V Tretyakov. Sampling and esti- mation on manifolds using the langevin diffusion. In Advances Neural Information Processing Systems, 2025.
Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R. Zhang. Sampling is easy as learning the score: theory for diffusion models with minimal data assumptions. arXiv preprint arXiv:2209.11215. 2023.
Xiang Cheng, Jingzhao Zhang, and Suvrit Sra. Efficient sampling on Riemannian manifolds via Langevin MCMC. Advances in Neural Information Processing Systems, 35:5995-6006, 2022.
10

---

Published as a conference paper at ICLR 2026
Xiang Cheng, Jingzhao Zhang, and Suvrit Sra. Theory and algorithms for diffusion processes riemannian manifolds, 2023.
Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to score-based generative modeling. Advances in neural information processing systems, 34: 14417695-17709 2021.
Valentin De Bortoli, Emile Mathieu, Michael John Hutchinson, James Thornton, Yee Whye Teh, and Arnaud Doucet. Riemannian score-based based generative modelling. In Alice H. Oh, Alekh Agar- wal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022.
Prafulla Dhariwal and Alexander Quinn Nichol. Diffusion models beat GANs image synthesis. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, 2021.
Nic Fishman, Leo Klamer, Valentin De Bortoli, Emile Mathieu, and Michael Hutchinson. Diffusion models for constrained domains. arXiv preprint arXiv:2304.05364 2023.
Khashayar Gatmiry and Santosh S. Vempala. Convergence of the riemannian langevin algorithm. ArXiv, abs/2204.10818, 2022.
Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73:123
214, 03 2011. doi: 10.1 (1111111677788.111000000000
Alfred Gray. Tubes, volume 221. Springer Science & Business Media, 2003.
Yunrui Guan, Krishnakumar Balasubramanian, and Shiqian Ma. Riemannian proximal sampler for high-accuracy sampling on manifolds. arXiv preprint arXiv:2502.07265, 2025.
Richard S Hamilton. Matrix harnack estimate for the heat equation. Communications in analysis and geometry, 1(1):113-126, 1993.
Qing Han and Qi S Zhang. An upper bound for hessian matrices of positive solutions of heat equations. The Journal of Geometric Analysis, 26(2):715-749, 2016.
Ulrich G. Haussmann and Étienne Pardoux. Time reversal of diffusions. Annals of Probability, 14: 1188-1205, 1986.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and Lin (eds.), Advances Neural Infor- mation Processing Systems, volume 33, pp. 6840-6851. Curran Associates, Inc., 2020.
E.P. Hsu. Stochastic Analysis on Manifolds. Graduate studies in mathematics. American Mathemat- ical Society, 2002. ISBN 0180821808023,
Chin-Wei Huang, Milad Aghajohari, Joey Bose, Prakash Panangaden, and Aaron C Courville. Rie- mannian diffusion models. Advances in Neural Information Processing Systems, 35:2750-2761, 2022.
Zhihan Huang, Yuting Wei, and Yuxin Chen. Denoising diffusion probabilistic models are optimally adaptive unknown low dimensionality, arXiv preprint arXiv:2410.18784, 2024.
Aapo Hyvärinen and Peter Dayan, Estimation of non-noomalized statistical models by score match- ing. Journal of Machine Learning Research, 6(4), 2005.
Jürgen Jost. Riemannian Geomelry and Geometric Analysis. Universitext. Springer, Cham, 7 edition, 2017. ISBN 978-3-319-61857 doi:
Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence of score-based generative modeling for general data distributions. Proceedings of the 34th International Conference on Algorithmic Learning Theory, volume 201 of PMLR, pp. 946-985, 2023.

---

Published as a conference paper at ICLR 2026
Gen Li and Yuchen Jiao. Improved convergence rate for diffusion probabilistic models. In The Thirteenth International Conference on Learning Representations, 2024.
Gen Li and Yuling Yan. Adapting to unknown low-dimensional structures in based diffusion models. Advances in Neural Information Processing Systems, 37:1 16297-16311,,3,, 2024.
Gen Li and Yuling Yan. O(d/t) convergence theory for diffusion probabilistic models under minimal assumptions. In The Thirteenth International Conference on Learning Representations, 2025.
Gen Li, Yuting Wei, Yuxin Chen, and Yuejie Chi. Towards non-asymptotic convergence for diffusion-based generative models. In International Conference on Learning Representations (ICLR), 2024.
Junfang Li and Xiangjin Xu. Differential harnack inequalities riemannian manifolds i: Linear heat equation. Advances in Mathematics, 226(5):4456 4491, 2011. ISSN 0001-8708. doi: https: 1ttr://1.1016/..aiii...
Mufan Bill Li and Murat A. Erdogdu. Riemannian langevin algorithm for solving semidefinite programs. Bernoulli, 29:3093 3113, 2023.
Peter Li and Shing Yau. On the parabolic kernel of the schödinger operator. Acta Mathematica, 156: 153- 07 1986. doi: 10. 10.1007/FF02399003
Jiadong Liang, Zhihan Huang, and Yuxin Chen. Low-dimensional adaptation of diffusion models: Convergence in total variation. arXiv preprint arXiv:2501.12982, 2025.
Yang Liang, Peng Ju. Yinghin Liang, and Ness Shroff. Non-asymptotic convergence of discrete-time diffusion models: New approach and improved rate. arXiv preprint arXiv:2402.13901, 2024.
Guan-Horng Liu, Tianrong Chen, Evangelos Theodorou, and Molei Tao. Mirror diffusion models for constrained and watermarked generation. Advances in Neural Information Processing Systems, 01-22898-42111 2023.
Aaron Lou, Minkai Xu, Adam Farris, and Stefano Ermon. Scaling Riemannian diffusion models. Advances in Neural Information Processing Systems, 36:80291-80308 2023.
Michelle Muniz, Matthias Ehrhardt, Michael Günther, and Renate Winkler. Higher strong order methods for linear itô sdes on matrix lie groups. BIT Numerical Mathematics, 62, 01 2022. doi: 10010077/115533021-0009999
Peter Petersen. Riemannian Geometry, volume 171 of Graduate Texts in Mathematics. Springer, New York, 2 edition, 2006. ISBN 978-0-387-29246--. doi: 00.1007/978--0-88-22993---
M. Piggott and Solo. Geometric euler-marus maruyama schemes for stochastic differential equations in so(n) and se(n). SIAM Journal on Numerical Analysis, 54:2490 2516, 08 2016. doi: 10. 137/ 15M1019726.
Yury Polyanskiy and Yihong Wu. Information theory: From coding learning. Cambridge univer- sity press, 2025.
Peter Potaptchik, Iskander Azangulov, and George Deligiannidis. Linear convergence of diffusion models under the manifold hypothesis. arXiv preprint 1rXiv:2410.09046 2024.
Daniel Revuz and Mare Yor. Continuous martingales and Brownian motion, volume 293, Springer Science & Business Media, 2013.
Richard M. Schoen, and Shing-Tung Yau. Lectures on differential geometry. Conference proceed- ings and lecture notes in geomelry and topology, 1994.
Jascha Sohl- Dicksteiin Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In Francis Bach and David Blei (eds.), Proceed- ings of the 32nd International Conference on Machine Learning, volume 37 Proceedings of Machine Learning Research, pp. 2256-2265, Lille, France, 07-09 Jul 2015. PMLR.
12

---

Published as a conference paper at ICLR 2026
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In Interna- tional Conference on Learning Representations, 2021.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In Neural Information Processing Systems, 2019.
Vishwak Srinivasan, Andre Wibisono, and Ashia Wilson Fast sampling from constrained spaces ing the Metropolis-adjuutee mirror Langevin algorithm. In The Thirty Seventh Annual Conference on Learning Theory, pp. 4593-4635, PMLR, 2024.
Hajime Urakawa. Convergence rates equilibrium of the heat kemells on compact riemannian manifolds. Indiana University mathematics journal, pp. 259-288, 2006.
Xingyu Xu and Yuejie Chi. Provably robust score-based diffusion posterior sampling for plug- and-play image reconstruction. Advances Neurul Information Processing Systems, 37:36148- 361 84, 2024.
A NUMERICAL EXPERIMENTS
In this section, we verify the results in Theorem 1 on compact manifolds by measuring the exit probability in the reverse steps of Algorithm I, which is extensively used in the proof of Lemma 19 to ensure the convergence of Algorithm I, and the TV distance between the target distribution (a Gaussian mixture) and the recovered distribution.
Reset probability on S and T We start by examining the total reset probability on the unit 2-sphere S and the 2-torus T We the backward process in Algorithm with different stepsizes each setting with po being Gaussian mixture (see below for definition), and record the fraction of trials whose tangent update As = St, (Yx) + VT.Gr violates || < 1/44 at least once among all steps. On both manifolds, Figure la shows a clear linear trend of the logarithm of reset probability of against h 1-1/2 which can be obtained by Gaussian tails. This confirms that the rejection sampling has no practical impact on the performance of Algorithm 1.
High-ddienniionl torus. We extend this experiment the d-dimensional flat torus T for differ- ent stepsizes. Figure 1b reports the logarithm of the reset probability versus h- for d {2, 4, 8). Increasing raises the baseline reset rate, yet the slope of the decay remains essentially un- changed-resed resets remain exponentially rare h O, which aligns with our analysis.
TV accuracy on To with warped Gaussian mixture. Finally, we assess the distributional accu- racy in TV for warped Gaussian mixture target T d 21, Here, the warped Gaussian distribution is defined to be the push-forward of the Gaussian distribution by the universal covering e22ma.), )) and warped Gaussian mixture is similarly the push-forward of Gaussian mixture in Rd The result is depicted in Figure 2, which confirms that the total variation decays fast with the increase of the number of steps.
F Tdd given by (....... za) + (m
B PRELIMINARIES
We first introduce some tools we use in the rest of the proof.
Lemma 3 (Pinsker's inequality, Polyanskiy & Wu (2025)). For any probability distributions q on M, we have
TV q) S /2KL(p
Lemma 4. Let be a vector field on M. In a local coordinate on M, we have
0.43 = (Vau)² T
Toyu4. Uy-
Here D ony is the Christoffel symbol, defined as
8 07 = 9 9 (Dagrs + a,gas 0sgaz).
13

---

Published as a conference paper at ICLR 2026
Reset et probabiitty on different manifolds
Reset probability for torus of different dimensions
Fiuzion
Muuiin
[am4)
Ing ewser pww.abiit
-12
7.5
10.0
12.5
15.0
17.5
20.0
22.5
25.0
12.5
17.5
20.0
22.5
25.0
27.5
b-in
h-10
(b)
Figure 1: Reset probabilities on spheres and tori. In Figure la, we examine the relationship between h 1/2 and the log of the reset probability of Algorithm 1 on both sphere S and torus T under the reset rules of Algorithm 1. In both cases,
that the reset probability decays exponentially, confirming the conclusion in (43). In Figure 1b, we examine the same statistics on high-dimensional tori, and find increasing d only shifts the curves to the right but leaves the exponential decay rate in h 1/2 unchanged.
Goussian-mixture target
1.0
d-3
0.6
0.2-
Number
Figure 2: TV distance on T with a warped Gaussian-mixture target. The total variation is estimated with kernel density estimator.
Lemma 5 (Metric distortion in normal coordinates). There exist coefficients c, C >> 0 polynomial in d and constant in other parameters, such that the following holds. Let E M. In the normal coordinates (0a) at x, for any E M such that p(æ, y) I c/ K, we have
llo(v) - 1|| [1 CKd (z.1), ||0 < CKp(C v), 10009 - CK.
Proof. This a quantitative version of the well-known Taylor expansion of g in normal coordinates (cf. Berline al. (2003, Proposition 1.28)):
1
(ns(exp,. = 5,, ( + [(]][m]] [VRm)))] 3
||u||| c/K.
Let (e, )a= be an orthonormal basis M,M, normal coordinates at x are defined by identifying a point exp, (a) with its coordinate vector (z²). Denote by O, the coordinate vector fields.
14

---

Published as a conference paper at ICLR 2026
Step I: Representation by Jacobi fields.
Define the geodesic segment with unit parameter
[0, 1]: 7(8) expx
d 7(0 2; 7(1) v, {(s)
ds (
Then v and
[+ss| | rr.
Fix an orthonormal basis (ea) of Tf, M and parallel transport it along y to obtain an orthonormal frame (Ea(s) along ?
V. E.E (s) 0, (a( Ca.
Let (s) be the Jacobi field along corresponding to varying the initial point in direction eg in the normal coordinate chart, i.e.
JØ(s) == (se E ET(s)
Equivalently, Jg is the unique Jacobi field solving
ViJa + Rm (Ja, 9)4 = 0,
J=(0) -0 VaJa(0) = eg.
(7)
Note that in normal coordinates, 0z1= = 6B and the geodesic variation expa (s(r EeG)) yields (7)- Write (s(s) in the parallel frame:
JØ(s) Σ (sa(s) (s)(;);
and let ((s) Radd be the matrix with entries (ss(s). Since E,, parallel, (7) becomes the matrix
Jacobi equation
('(s) + J(s) = O, where the curvature matrix R(s) is defined by
((O) = 0, J'(0) =
(R(s) u) = (Rm Σ up Ex)) En(s)
From |Rm|| S K and |+| = r we have
R[(s]] < K1)()) = Kr², 8 s [0,1]
(9)
In normal coordinates, the coordinate vector fields at y = (11 are
0g y = d(exp.) (es) Jz(1).
Since the frame at = 1 is orthonormal, the metric coefficients are
g (v) = (8ply, - (J,(),())))) Σ
(an(1) (1) ())) By
Hence, matrices,
gly) = J(1)
(10)
Step II: Control of J(1) I via Grönwall inequality. From (8), integrating twice and using J(0) = 0. ('( = I, we get the exact Volterra equation
J(s) = sI [ (s- (s-
8 8 [0, 1].
(11)
Taking operator norms and using (9) gives for & E [0,
[[s)]] s Kr² / (s- - [[(G]] dr.
A standard Grönwall argument yields
[[(s]] Cs and ||((s)]| c
for all E [0, 1], provided <c/VKK
(12)
15

---

Published as a conference paper at ICLR 2026
Now subtract s/ in (11):
((s) - -sl sl== =- -==], (s-))
Using (9) and (12),
[J(1) - S. ([->)|R[]]]]]G]]] dr S Sc (1-)) / SKK2
Combine with (10):
slv) 1= (J(1 = 1)" + (J() 1) + (J(1) 1) "T (J(1) 1),
so
(l( /|| S. CIJ(1) /|| CK'2.
(13)
Step III: Control of first derivatives. We first control (1) function of the coordinate t. Let +> (s) denote the Jacobi matrix for the geodesic 7x(8) expx (su). Differentiate the ODE
(8) w.r.t n²: J,J" + RØ,J = J,,R),, (aJ(0) 00 da) = 0.
(14)
k)) y in the parallel frame. Varying changes both y and ;; the corresponding variation field Va(s) := da (o(s) along y is itself Jacobi field with V. (0) 0, (O) Ca. hence by the same estimate as (12)
Bound on O²R. Recall R(s) represents the operator
[Va Cs.
(.). < C.
(15)
Using the product rule on Rm(- y and our Assumption 1, together with |il = r and (15), one obtains the uniform operator bound
[()]I < C [Rm|| (s)) 1+ + |mm|| H [[0+)())]]) s C(K + K K.r.1 s CKr,
(16)
for all [0, (since s 1). Here used 0, V.Va.
Now solve (14) by the same Duhamel principle: integrating twice with zero initial data gives dr.
0.J(s) L- 7)
(17)
Using (9), (16), and (12), we obtain
[s),](s]| s Kr² ( (a-)) (+) (s-))((ss(())))) S Kr² J( J (s-T) [1..
Apply the same Grönwall comparison before (now with forcing term CKr s²)) conclude, for r S c/vK,
||a,((1)1I S CKr.
(18)
Finally differentiate g = JJJ:
o.g = (0.J) J(O.J)
so by (12) (at 1) and (18),
[(ar(r)) s 2 ((OI & C(Kr) 1 s CKr.
(19)
Step IV: Control of second derívatives. Differentiate (14) once more: e + R R. = (00..R) J (0,, 0zJ - (0gR) a.
(20)
with zero initial data at s = O.
Bound R. Under Assumption 1, R. can be bounded uniformly by CK on [0, 1] as follows: expanding the second parameter derivative Rm(-, () produces terms of the schematic form
(VRm)(V) .4 (04), Rm(;84) (04),
(VRm)(aV) -4
16

---

Published as a conference paper at ICLR 2026
and also terms involving V (OV), all of which are controlled using || S 1. V < 1 and the fact that each appearance of y contributes a factor r. Concretely, one shows (using (15) for both Va Vi and the same Jacobi estimates for their derivatives) that
[[),)((s)Il & CK
for all E [0, 1].
(21)
Now apply Duhamel's principle to (20) with zero initial data:
(ss)(s) = - [ (s ---r) r) (Roast+ (DaaR) J+ (0aR)0gJ (0gR) (r)) dt. Take norms and use (9), (21), (12), (16), (18):
[[aa]( $ Kr² /-- (s (s- (()))(()))) C ' [
(U(T] (T ' (s-T) dr.
Since |J(T)]| Ct, the second integral is bounded by CKs², and the third is bounded by CK2,23² 22 provided c/v e/. Thus for 1,
[+)s)( - Kr² (((-->) (o--))l1ag (1aa](T)] d CKK.
Applying Grönwall argument yields
[[0a)(11]|| CK. CK..
(22)
Finally differentiate g JJJ twice:
Dagg (0asJ) JJ J (DasJ) + (0.J) (0gJ) + (0:J) (0.J).
Hence by (12), (18), (22) (and K2,² < CK for r S k/vK),
|[a s C'Il,,s - [|((1)] C(Ioa(()) | S CK CK+ CCKr) SCK.
(23)
This completes the proof.
Lemma 6. Fix x M. Define
Jar u) = del dexp, = d (exp,, u).
There exist universal constants c, C > O, such that for u T&M with ||u|| S kā. have the following bound on J [xx u):
[(x,2) -1 -1/5 [CdK||u22
(24)
In particular, we have
1 < [z u) 2, ||u|| S Kd'
Proof. Work in normal coordinates at so that exp, Baue (0, 1/K) C T.M - Broa (x. 1/K) is diffeomorphism and 9t, (0) 6,y. m (0)
From Lemma 5, we know that lgg u) - /|l - CK uu the region ||u|| S c
we have
Kd.
[la(exp, u) /|| S
Therefore, by Taylor expansion of determinants, know
det g(exp, u) 1 = |det(1 + g(exp, u) 1) - 1
< Ctr(o(eex, u)
< Cd |la(exp, u) I|l
< Cd CKI|
This concludes the proof by adjusting C' if necessary.
17

---

Published as a conference paper at ICLR 2026
The metric distortion bound implies that geodesic is almost straight line, in sufficiently small normal neighborhood. The following quantitative bound shall be useful.
Lemma (Geodesics are almost straight in small balls), There exist coefficients e, C 0 polyno- mial in d and constant in other parameters, such that the following holds. Fix any M and let < r & K. Let y, 2 E (r) and y be the unit- unit-speed geodesic connecting y to 2. Write
y(s) = exp, (r(s)) T&M = R
for its representation in normal coordinates at x. Then:
(i) (Almost constant velocity)
sup ((s) 5(0) I CKr². sE(0,8)
(ii) (Almost linear trajectory)
sup =E(o,
((() v(0) - s(0) & CKr
In words, in normal coordinates at x. any geodesic segment contained in Br(r) deviates from the Euclidean line segment connecting its endpoints by at most O(Kr²) in position and O(Kr2) in direction.
Proof. Work normal coordinates x. By bounded geometry and the choice of r, the metric coefficients satisfy, in view of Lemma 5, that
( - I|| < CK
|19(())|| - CKM.
y rr
which implies the Christoffel symbols obey
[(w)] S Ckkbl CKr.
The coordinate representation y(s) of the geodesic satisfies the geodesic equation
y (s) + (y(s) y (s)g² (s) =
Since is unit-speed and gvy) is uniformly equivalent to the Euclidean metric on ly] S r, we have ()(s) x Consequently,
[ū(s)| < CKr
for all s E [0,
Integrating once gives
(s)) - ((0) < / [ū(w) du s CKrs SKrssCKrr
proving (i). Integrating again yields
(u(s) y(0) so(0)) ' [ [[(w)) dw S CKr 3 < SCK'
which proves (ii).
We spell out the constants in a few classical inequalities in geometric analysis. Lemma 8 (Schoen et al. (1994, Thm. 4.6)). Let (M,g) g) be complete Riemannian manifold with Ric(M) -K for some K 2 O. Let H H y, 1) be the heat kernel, i.e., the fundamental solution of (A & )u(x,1) 1) 0. Then, for every 8sch >> O and
1,
, (2
H (1, < C'osch d, a) Va (v() (/3)) VV()) - 1/2
exp (4 dsch.)
+ C1 8seh Kl
where V.(R)= ({BB(R)) C(Ssch, d, + (gch)) da exp (te). and C1 ad
18

---

Published as a conference paper at ICLR 2026
To unleash the power of Lemma 8, we need the following lower bound on volume of geodesic balls assuming bounded geometry.
Lemma (Günther's comparison theorem). Under Assumption we have
(2TTd//2 (a(r) T(d/2)
[-"
K)
dt, OSrS
VK
Here D is the Gamma function. In particular, when S c/ K for some small universal constant C > O, we have
d/22
1
(2(r) (d/2)
da/2
Proof. We observe that ||Rm|| < K implies that the sectional curvature is upper bounded by K, since by definition sec(u, v) Rm(u, v, u, v). The first inequality follows from the classical form of Günther's 's comparison theorem, see for example Gray (2003, Theorem 3.17). The second inequality follows from the elementary bound that sin(x) 2 2 for E [0, c). where is a small universal constant, and the crude bound ((x) S =- for 2 1.
complementary lower bound to Lemma 8 is as follows.
Lemma (Li & Xu (2011, Thm. 1.5)). Let (M,g) ) be complete, possibly with Ric(M) 2 K. For the (Neumann) heat kernel H 1, and all v, 1 M, 1> 0,
(2Kt)a/2
p(x,y)?
Kt Ktcoth(Ktt
((t, x. v) 2 (4Tt) d/2
(e2kt
2Kt
1)d/4
exp
4f
(1+
Kt
H(t, 2 (4Tt) d/2 exp
p(z, (
(1+ ktt
d
Kt
(25)
4t
4
The above bounds for heat kernel translates seamlessly to Pt, since pe is convolution of po with (t, .,)). We formalize this in the following lemma.
Lemma 11. We have
inf (t,,,,)) S inf pr(x) S sup (t( S sup (t,, v).
a.
=EM
zEM Z.E
This follows from taking infimum and supremum respectively in the formula (Duhamel prin- ciple)
Proof.
(v()) 1 po (1,
The following lemma compiles a few follow-ups of Li-Yau estimates (Hamilton, 1993; Han & Zhang, 2016) with constants made explicit.
Lemma 12. Under Assumption we have Han-Zhang 's inequality
P2pt 3 CHz (++ + (r sup
Pt
Pt
On the other hand, we also have Hamilton's 's Harnack inequality
V2P:
(Vp.)
1 2 CHam 1 log
sup supPt/22
2 PPlogp:
g
Pt
Pt.
and
V
sup 1t/2
V log
S C
+ KK
log
p
Pt
Here C > 0 is universal constant, Chz Cd,,, Ctam CdK?
19

---

Published as a conference paper at ICLR 2026
Proof. The last inequality follows from Hamilton (1993, Theorem 1.1). The proof of the rest two inequalities require tracing the proofs of Hamilton (1993); Han & Zhang (2016). The details would be too tedious reproduce here, we leave pointers relevant proofs for interested readers.
To prove the second inequality, we trace the proof Hamilton (1993, Theorem 4.3) to see that A > 0 is such that
AMP Pt
A
Pt/2
d log
Pt
and
Vpll
A
sup Pt/2
S
d log
Pt
then
VPPc
(Vpi) (Vp))(P)))
21 9 CA[]] + [Rml]))
Pr
Here C > 0 is an absolute constant. By Assumption 1, this implies
CHam S CHCHN
Tracing the proof the main theorem in Han & Zhang (2016, Page 9), we can see
CHZ S CLy(1 K) CdK,
where CLY is the maximum of the coefficients before t-1 and K in Li & Yau (1986, Theorem 1.2). This was explicitly defined as Cd therein, by setting
there. This completes the proof. G
Lemma 13. For any three points x, y, E M, we have
2)
p(z,))²
p(x,y)²
AI (0,
1
The equality is attainable at some point 2, on the minimum-length geodesic from x to y. Moreover, x. y, z are within < 1/ poly K) distance to each other, then the function
((2) = p(z,2)2 2(2.8)2 p(x,8)2
-
is 1- =(1-1) -strogly convex in the normal coordinates at x (or y). where C > 0 is universal constans.
Proof. The first inequality follows from Cauchy-Schward inequality and triangle inequality:
p1x, [)(x,2)² 2)² p(z,))² 1-t
(1 1 1)
2 (p(x, 2) + p(z,u))² 2 pla, y)
Let 7 : [0, + M be the constant-speed, minimum-length geodesic from x to y. then straightforward to check that
2(z, p(z,,u)² p( v)², 2. ===)()))
1-t
where X2 E (0, 1) solves the quadratic equation
(1=
1-
Proving the strong convexity requires Lemma and Lemma 6, which implics that p(z, 2)2 and p(z, 8 are both (1 - Cdl K (2)))-strongly convex in the normal coordinates. The desired conclusion then follows from
1 1 1
t t
1(1-t)
The proof is completed.
20

---

Published as a conference paper at ICLR 2026
C INITIALIZATION ERROR
We require the following result from Urakawa (2006, Proposition 2.6). Lemma 14. Denote H(t, x, y) the heat kernel on M. Assume
A = sup sup (d12 (1,
1g zEM
Then for any probability distribution po, its evolution along heat flow d.p: = AAm satisfies
TV( Vee (1-(1-)) 121.
We combine this with the Li-Yau upper bound (Lemma 8) to obtain
Lemma (First part of Lemma 1). Under Assumption 1, there exists universal constant C such that
TV(pN, μ) & C(K+dlogKd)) (T--)
Proof. Plug the bound in Lemma 9 into Lemma 8 and use the fact that supx.y HI x, y) is decreas- ing (by convolution inequality), we obtain
A s (Cd/K)&ecK
S exp(C' CKd)o
for some universal constants C, C" > 0. Note that we absorbed dlogg K into K + dllog d. The desired claim follows.
D SCORE MATCHING ERROR
We now prove the second inequality in Lemma 1. Lemma 16. Under the same assumptions as in Theorem we have
N
Eli (Y)) - Hion (VI) S 22
Proof. This is relatively straightforward. Notice that in normal coordinates, by Lemma have
1 (Yi) Fi.y, (Y2) = (a.((Y) log (S,,, (Yr,) log Pt, (Yi)))
- |g)() || (Yr,) V Po, (Yr, )|
< & (Yr,) V Voggt, (12
for p(Y,, Yu.) S e/K, and is 0 otherwise due to our cutoff Tw. Therefore
El (Yi) - Fi (Y:)1² < 2 (Yr,) V log (Y,,))1
and consequently,
N
ElI.Sryyy, (Yi) Piv (Y))
N
tt
1 E (Vr,) VlogPry (12 202
t
k=1
This proves the claim, desired.
21

---

Published as a conference paper at ICLR 2026
E DISCRETIZATION ERROR
Lemma 17. Under the same assumptions as in Theorem and assuming (3) without loss of gener- ality, there is universal constant C such that for ty. S S tke we have
Cd K
E log (V:) P Nn (V)11
(La 1).
f3
Proof. For convenience, set the reverse time
T 1k
The main challenge is that Li-Yo Lau estimates provide sharp uniform control up to second-order derivatives of log Pt. but a naïve calculation of the difference V logg (Y)) - - -YyY, (Yi) involves third-order derivatives. More precisely, straightforward Taylor expansion will introduce factor of 0, V log Pt, which, by reverse-time heat equation d-p: = 2 A MPt, contains third-order derivatives of Pt. We bypass this difficulty by making use of Itô's calculus show that third-order derivatives cancel out; this is inspired by Benton et al. (2024), where a similar strategy was employed to the Euclidean setting.
Step I. Applying Itó/Stratonovish formula. We first compute O, V log Pt. Since the forward heat equation d.pt- APt, have
Oepr
AMP:
a v logp:
O logp, =-V
V
Pt
Pt
Use the manifold quotient rule and the identity AM log
AMPL
V
logp:|1 to rewrite
pt
AMP:
AM + || 1ogg.112
Pt
Therefore, have
1
A AMPi
V
VAM log V 1
2
Pt
VAM log
22 Pogg V log
2
On the other hand, it is straightforward calculate
V
(VP))
V
logge
Pt
Pi
Now, from Itô's formula, we know
dV log (Yi) (0, v (o + V log Pa (Yi) (oggo Uy, dw.) A log))
2 VAM
v2 logg: V log + V2 log V p.dr
D Uv, dW, + Am logp,d7
(AnV- 7 + log Uy, dWW
2
1
Ric log )dr + 2 logpe Uy, odW,
Ric (Vf.-), and Ric" de- notes the (1, 1)-tensor obtained by raising one index in Ricci curvature. Notice here the cancellation of order-order derivatives.
where the last line follows from Bochner's identity (AMV-VAM) Ric
On the other hand,
d
N, (Yr) = V. .VV.(Voog + Uy, o dW,) + AmM (Yi)dr.
22

---

Published as a conference paper at ICLR 2026
In normal coordinates, Nu is constant vector field inside B(O,w/3), therefore we have (ef.
Lemma 4):
9 ty Ye (8)) I aay a Ve (x)) = V Pog Pas (Yu), P()), S w/3.
and similarly, when p(Y,, Y(,) S w w/3, have
A AmTrr tue (Y)) - V log Pru (Yr) -aal)-g--- get (ar T
r +
B7) V log Pte (Vr)).
Step Bounding the coefficients. Combine the above formulas with the estimates given in Lemma 5, we obtain for some universal constant C >> 0:
Ric V logpt
CKIVV Po Pa (Y:)II,
V. Yo
C Pog (Y,,)l,,
p(Y,, Yox) S w/2, p(Y,, Y..) S w/2.
AM Pu
CK*
V log (Yr,)l,,
Outside the geodesic ball By, (w/2), the field Fr,Y, is non-zero only inside By, (w). Between these two balls, we have to take into account the radial derivative of the cutoff function T whose first order derivative is bounded by Cw and second order derivative by C by construction of Ty (recall that lnyl + |7/1 S 100). Apply Lemma and Lemma 5 again, this time we bound
V Mi
- CKWw-1I log (Vr,)l,,
A
WCK2
| Pog ptu (Yr,)) )l
We now apply Itô's formula manifold (i.e., take expectation and invoke the martingale property (1)) and collect the above bounds to obtain (cf. Benton et al. (2024))
d dr
V (Y.) Fixx, (1:)
$ CK² (BI|V pr (Y;) E Po, Pt.. bog
+ +CK 2E log (Y log (Y
(Y)))))))) y,,,
CK2 (BIVVogn
Pt, Pt,, (Y, (o
WCK logp: + B lo + Logg /(p((Y. )n)) w/3)
SCK2d (++
E Pa/2
(Y,) (1 1 3 (.)) w/3)
sup (5.5
(26)
Here the last line follows from Lemma 12.
Step III. Controlling expectations via Chebyshev and Li-Yau estimates.
To bound
E Elog 1 sup
note that
(V,).
1 pidu / du du=11
E
Pe (Y.)
Pt
By Chebyshev's inequality, we have
(Y)) 21 5x-, 1=0,
and then
log'
1
21 Se e 5x A
(Y,)
23

---

Published as a conference paper at ICLR 2026
Integrate with respect to A,
1
E log" + Pr(()) SC.
We then apply Li-Yau's estimate (Lemma 8) combined with Lemma 9, Lemma 11 to obtain sup Po/2 sup log (s x. y) S dllog K's, where the first inequality follows from Ps/2 being the convolution of po with (s/22, x. v). These together shows
sup Ps/2 V Cdlog + CKs)
CdK
E log
(Y,)
8
Here used log < c 1/4 and 2 2626.
Step IV. Controlling exit probability via stopping time. remains to bound the probability PYp(Y,, Yo,) w/3). This would follow from stopping time argument. We claim that given t --
t S h, we have
P(p Y,.) > w/3) = exp
S
w
(27)
ty
where the last inequality follows from (3). Plug this back into the desired conclusion of the lemma proved.
We now prove (27). Let a be the largest < tx such that p(Y,, Y..) > w/3. We have P(p(Y,, Y.,) w/3) P( 2 t).
In the interval [o, tu]. Y stays in the geodesic ball By, (w/3). and follows the SDE (2). In normal coordinates, this can be spelled explicitly:
dy = Va lo + A4 o dw?
1 v (Y)) + (0,
dt + Agaw?
(28)
2
where A is the square root of the matrix representing the coefficients of the Laplace- Laplace-Belrammi Beltrami oper- ator
det (.(V det 0g)
can be checked with the help of Lemma 5 that ||A - Ill I CdKw?, and ||,A4|| < CdK omit the computation that has similar pattern as many of the previous arguments. Further- more, we have log loog).(Y)) x (8-1 + K) log supy2a by the same argument via Lemma 12 as pt before. This time combine the uniform bound provided by Lemma 10 with Lemma 8 clude log sup $ (8-1+ K + do d)² Diam(MM² $ (8-1 K K d)² K22 by Assumption 1.
This shows
sup
logpell & (8-1 + d doogd)3K K?
(29)
Therefore, in view of Lemma to convert the above bound to normal coordinates, and together with the aforementioned bound for A and 0,A,
that the drift term up to time will not exceed
sup/ log Pe (Y+) (0, (7)4))177111 (tx o) - c(8-+ K + d log d) K2 a)) S T2
(30)
where the last inequality used (3). On the other hand, the bound on A implies that the quadratic variation of the martingale part does not exceed
|"
AC A2dt = 2(tx IoI.
By Burkholder- Drsiiiderr--vvi---ee Davis-Gundy inequality (Revuz & Yor, 2013), the tail of f Agaw; is 0(1)- subgaussian, thus we have
P a 2
|1 Agaw?l|
-olw 2 - t))²
P
< exp
12
l
In view of (3), combine this with (30) and (28), we have proved (27) claimed.
24

---

Published as a conference paper at ICLR 2026
Lemma 18. Under the same assumptions as in Theorem and assuming (3) without loss of gener- aliry, the discretization error obeys the following upper bound:
Cd K
E Vogg. - Pi U.YY (Y))
83 2
where C is a universal constant.
Proof. This follows directly from Lemma 17.
F BROWNIAN MOTION SIMULATION ERROR
In this section, we handle the Brownian motion simulation error using the machinery of Minakshisundaram-Poee parametrix. complete introduction to this heavy machinery would quire establish whole system of notations and lemmas in geometric analysis, which is unduly burdensome. We instead refer the interested reader to Berline al. (2003) for a comprehensive treatment, and point to results there whenever needed.
F.1 OVERVIEW
Our aim is to prove the following lemma. Lemma 19. Under the same assumptions in Theorem and assuming (3) without loss of gener- ality, we have
TV(p g6) S KT poly(d, k,
To better explain the idea of the proof, we ignore the rejection sampling procedure in the construction of K2 temporarily. Our starting point is the observation that by Makkerr-Plance equation, K2 is the heat kernel associated to the Euclidean Laplacian with drift St.YY,. in normal coordinates. On the other hand, Kat is also a heat kernel with the same drift, but associated to the manifold Laplace- Beltrami operator. The following lemma shows that the two solutions coincide up to first order in time, at least in polynomially small neighborhood of initial point and in polynomially short time.
Lemma 20. Let (1, x. y) be the (generalized) heat kernel for the operator H A + (Pawy V). Define the Euclidean density
(- 1|u 1i--Frnn Sen (x))(1
(e(l;z)= = (2it)d/2 exp
u T,M,
2t
and let (11,, y) be the density of the forh-forward by . of (:sP((;;))))) with respect to the volume
measure, where Yu the cutoff function defined (4). Then there exists polynomial poly K)
with universally constant coefficients, such that for all O < < poly(d.K,6 and all ply. Yr,) S
5/12 we have
F Yte.
S poly K, )).
(1, Yr,,
With Lemma 20 in hand, it is tempting to calculate the KL error with the following heuristic:
KaK || pk KR) S -E / pHt Yous .) log pH Yeas")
p(h, Yes.
(1, Yes. 2
$ E F
2hh, Y(..))
S poly K, 2n2
where we ignore the fact that Lemma 20 holds only in small neighborhood: the first line is post- processing incquality, and the second line stems from the fact that for two distributions p, q, we have
f plog q
-= = plog 1 9 p
p
25

---

Published as a conference paper at ICLR 2026
c1a- (P-p)
9 -
P c
2
p
p²
1 1 11-1
p
given - is sufficiently small, where the last line follows from p _ 1. From this, conclude that the accumulated error along N steps is bounded by poly K)h2N N poly (d, K)hT,
and the desired bound follows from Pinsker's inequality.
Apart from Lemma 20, above computation is the essence of this proof. The rest this section is mainly devoted to proving Lemma 20, and then formalizing the above computation by handling exceptional events of exiting the polynomially small neighborhood.
F.2 PROOF OF LEMMA 20: A PARAMETRIX ESTIMATE
We begin the proof of Lemma 20. For simplicity, denote by o the normal coordinate representa- tion of St. (Y,,). Naturally, our initial test solution the drifted heat kernel, as simulated by our discretized process:
(- [[ll utll?) 21
(r(u)= (2nt)a/22 exp
Rd
Before we compare this with the manifold heat kernel, there one subtlety we need to keep in mind. The density for is with respect to the Lebesgue measure on T,M, not with respect to the volume on M. We compute and define the corresponding density on M as follows:
|det d log, y1
(t,, (logg v) ()(A(z v), A(x,)) =
VE Bu(w).
det g(y)
Here quantities are computed in normal coordinates. The factor A(x, v) is known as the van Vleck-Morette determinant. From Lemma and Lemma 6, we know that
I Az 22 if p(x,))
c Kd
(31)
We consider the generalized Laplacian
H = -ーAA AM + (Dun V).
As in Lemma 20, denote by F the heat kernel of H at time to --- We also propose an approximation FH by
d y)
(1,,,,y) = G(I, exp(p((, v) (ax,))
G(I,,,y) = (2mt)d/2 exp
2t.
where for any two point v, E M, letting : [0, 1] + constant-speed geodesic connecting x to y, we define
sxx, v) = (( (7()).
ds 7(8) d>(s)) ds.
The auxiliary function V bridges FH and T in the following sense. On the one hand, relate + and v with the following lemma:
Lemma 21. There exists a polynomial poly(d, K)) with universally constant coefficients such that the following holds. For any By, (r), we have
S poly(d,K.5-T] and for all x, y
poly(d,K,3-T)
(t, x. 8) V (t, v)
polyydd, K, (r3 +t).
26

---

Published as a conference paper at ICLR 2026
On the other hand, we have the following asymptotic expansion:
(t, 2t)
Σ (z) - 0+,
where u, are smooth functions that can be computed explicitly via recursive formula (Berline al., 2003). We will not need the formula here, but instead require in and the remainder terms to be bounded properly. Such bounds have been well-established, which we wrap up into the following lemma. Recall the cutoff function Tue with radius w defined in (4). It is clear we can replace w with any > 0 to define cutoff of radius t.
Lemma 22 (adapted from Berline et (2003)). Fix a positive S po/yoodydd K). There exists smooth function in (I, y) on M x M such that
| o0 yalll & poly(d, K K,8-),, and all 0 < S 1/ (oly K), y E (2))), we have
((0, - H) In. (2(( (t, v) + (z, (r,(t,,,)) tr,
(32)
where
1,,4,, v, v) < poly (d, [ G(I,,,,), it
(33a)
(t, x, S poly(d, K, 10:56-1)) 1 p(2,))5 G(t, y).
(33b)
Proof. The inequality on un follows from Theorem 2.26 in Berline et al. (2003), with H the same our H and therefore F (in V). Note that all the coefficients in H are bounded in C2 by poly(d, K) + I,,,, (l))(M)), which is further bounded by poly (d, K, 1-) as we will show momentarily. In fact, by Lemma 5 and Assumption (A3), we have t1, [IC t,, Xi, 11~0 V. Mi, (i V Il poly(d, K) w ?1++ log Pt, (Y,, Il),
and then we control V Pt, (recall that w poly (d, K) by definition).
Pt.. (l,,). < poly (d, K, 11) via (29), yielding the claimed bound
We proceed to prove (32). We follow the proof of Theorem 2.29, item (iii) in Berline et al. (2003), and choose the cutoff function e there be 17 as defined in (4), and with the differential operator B defined in Berline et al. (2003), we have
((0. H) [V, x. (1 tun (2,, v))I
poly(d) G(I,,, 1. G((,,,y) 1(B (2) v)
Vy, related terms
Here By can be viewed as a coordinate-transformees version of H. applied to the variable y (precise definition can be found in the reference). and u2 is the second order term in the expansion. Similar to the argument we used to bound ||u| in virtue of Theorem 2.26 in Berline al. (2003), we have
1 (.)())) I poly(d, K, ))).
The claimed bound follows from combining the above inequalities.
We now state the Volterra series representation of heat kernel. Lemma 23 (Volterra series, Theorem 2.23 in Berline al. (2003)). Fix a > Let v (1, 77.
tm (2,
xt) = (a, - H) v (t, x x,
Define the time-space convolution operator * as
+ (),,))) (-] f(I- - s, 2)g(s, z, (dd(dsss..
Then we have
20
=H=
>
1-1) v
where r =
k=1
on any domain such that the series the right hand side converges absolutely uniformly.
27

---

Published as a conference paper at ICLR 2026
Lemma 24 (Iterative bounds for Volterra series). There exists polynomial poly(d, d, K, 55.8-1 with universally constant coefficients such that the following holds. Assume 0 < S poly(d, 1))-)) take 1 = 4/5/12 d in the definition of V. Then we have, for all p(x, y) & +5/12 (Ad), that
11 * r1 (1, S poly(d, K, K, x,y).
Proof. Recall Lemma 22, and denote by P the polynomial factor poly(d, K) therein. Denote by AA the dilated standard simplex
A+1
AA {(s).- Sk-1) S; 2 0, 8; = }}. A
Fix some h 2 Set 20 x and 2k = y: we have
V r1th (1,2,
d I
I th ph
ds
L
(11))) + 1s; ) (d))- (IR G( 2-1,
M
i=l
(34)
We split the integral in (34) into a local part and an outlier part. Define small "local" region
x( S 26, p(z=-1,2) 2.) S 1/2, 1,",, ,k]. R = {= ; =k-1)
We further define
[s) = Sfc II G 21))))) (de).
k
Tout (8) = f
II 1=1
G( 2-11 2,) 1 (1+ (lte)))))))))) (01))
1
187
It is clear that
1 Jea--1 (Moc(s) + Lour (s))
M + rrih (t, S shppp
(35)
We will establish bounds for hos and lout respectively.
Bounding the local integral. For ease of understanding, we begin by computing the first integral in Jie with respect to 21. Extracting the factors containing 21, we need to calculate
Jndensssnn (2T3s) dd/2 (2 d/22
1
p(2x211)
p(21.22)3
exp
m(dz1).
(36)
231
282
To proceed, we will invoke Lemma 13. Let x, be minimizer of
p(x,2)
222,22)
V(z)
p(x,, 22
$1 + 82) 82(s1 + S2) s2)
By Lemma 13, V(2.) 0. Since V(z) > whenever p(x, p(ar, 22), know that 2,
B. (2.). Moreover, Lemma 13 and Lemma 5 together imply
V(z) 2 (1 Cd)
2( V2 E Bz(4)).
$182
Here, the second inequality follows from strong convexity given by Lemma 13 and comparison of geometric distance and Euclidean distance in normal coordinates fueled by Lemma 5. Denote for the moment that
0
Caft R2,
28

---

Published as a conference paper at ICLR 2026
Plugging this back into (36), we obtain
J (0(2).))221) (2TT8) )u/2 (2TT2) dd/2 exp
2(2)122))) (l(2))
281
252
- (2T8)) 2/2 (2 dd/2 exp
V(zi)
((((2))
2(81 + 82)
< [ootaiossns (2 - s))d/2 (2TT)d122 exp
8(81 + 82 82) 21, 2.)
p(a, 22)
su(dz1)
28182
2(81 82)
((0,,22 22)2
+ 82)) / (2T (2Trz) dd/2 exp
2 0(s1
I 4ex
22(11 22(11
dz
28182
4(2 + 52)) G(s1 + s2, 22)
1
1
2T
8182
d/2
(2mr)) d1/2 (2T82))a/2
0(8) + 82)
S 40 G((ss + $2, x, 22),
where the third-to-lass line follows from change of variable to normal coordinates at 2, and from using (31) to bound the determinant; the penultimate line follows from Gaussian integration. Now, note that
0 g-4/2_ = (1 - cC 4/2 & cxp(2Cd ( 2,
give S 13 Putting these pieces together, we proved
1
p(x,
p p(21, 22)² u(dz)) S 8G(s1 + 82, 22).
(p(sua)sa) (2ms)) 4/2 (2m82)4/22 exp
281
282
Iterate the above argument for the integration over 22. -n- 1 to obtain
Σ S/ 2,
(sc(s) S g+ G
8 ,g(t,,z))).
(37)
Bounding the outlier integral. Next, we show how to control Lout. We first write
1,
G(s, 211 (2Ts;)))))
exp
(38)
28;
We claim that for any E Rº, we have
T II exp (-2 k )(
11 pp(s na)s (1+ 2s,(d (x,(d+1)
1 - p(s 122)21/2) S eexp 18;
16td
(39)
The claim is proved at the end of this proof. It is tempting to plug this back into (38), and argue that when is polynomially small, the integrand in Tout becomes exponentially small. However, this would not work since it does not resolve the singular factors I[= 8, d/2 in the integrand. For this purpose, we need the following crucial "freezing" trick, which follows trivially from 1 1+d a+1.
II 2(1,, 24 1 (1+ 1 1s,
(1 20 ( T.
The idea is to keep the Gaussian behavior resolve the s; d/2 factors, and only single out a very small proportion to demonstrate exponential smallness. Plug (39) into the above identity to obtain
(s) S exp (-
22
2
(11td)
[ +d +d= 21 23) p 11 (dk))) (1=1
29

---

Published as a conference paper at ICLR 2026
Integrate successively for each variable, convert to normal coordinates, and apply Lemma 5, Eqn. (31), and Gaussian integration as we did in bounding 4oc. we obtain
Jee [[G((Ldd)) 2-11 ((
1 (d2)))) - (1+d 4/2 - (ct)
(ct)-a/2
Therefore
12 16ld
lour()) S 32t (ct) d/2 exp
S 32 32 G(t,
(40)
where in the last inequality we used the assumption p( y) < 15/12 = (/d))) and t s
1/poly(d, K), that exp(-1//C0 (32td)) S exp(-p(x; y) (21)) and exp(-12/(322))
exp(-fdt-110) 1/6 (2mt)-d/2.
Putting things together.
We plug the bounds (37) and (40) into (35) to obtain
/ Jeat- (8+ +22 G(t 2, y)ds s
40
1 * *
x1. s th pk
(k - 1)! (40))
The desired conclusion Lemma 24 follows from the above inequality by summing over k and taking / S 1/ poly(d, K).
Proof of Claim (39). For 2 E R let
J {i ((2)-1, > 1/23-
By definition of R², either is nonempty, or there is in such that p(x, 2io) that
21. For i E J, we note
(_
(-
8s,d (1 1
1
62
exp
S exp
S exp
(1itdd]
(41)
28,d
18;
where the last inequality follows from 8, S 1, 415/12d, and that /. < 1/ poly(d, K). When J is nonempty, we readily deduce (39) as all the other factors are S 1.
When J is empty, let be such that p(x, zio) > 2c. We apply Cauchy- Sehoh--SchwarZ to obtain
p(5)-,
Σ
Ep(z)1122)
Σ p(2i-1)
Ei 8i
li=l
Then, by triangle inequalities, we have
ig
p 11 22
2t
(()-112)) 2 p(z0, zin) p(z,, zio) > 26, therefore
81
t
The desired claim (39) follows immediately, given that J is empty.
We now have all the ingredients prove Lemma 20.
Proof of Lemma 20. This follows immediately from Lemma 21, Lemma 23 and Lemma 24. Note that in applying Lemma 21, used < +5/12 thus 3 55/4 < t.
F.3 PROOF OF LEMMA 21
By definition, we can compute
@ x, v) V x, v)
exp (llgg y) - 101 v)) exp (-11
30

---

Published as a conference paper at ICLR 2026
Note that ||a|| < V (Y+,)] | by Lemma 5, which in turn is bounded by poly(d, 6 8-11 by
(29). When < voly(d,K.3-1) Therefore, it suffices to show
भा- we have exp(-Ilo||1/2) -1| z || s poly K, 6 8-1.
(logr y) v (2 3)1 S poly K, o
Recall the definition of th Note that since S poly(d.K,5-1j) when the polynomial poly(d, K, 8) sufficiently large, the geodesic 7 from x to y unique and inside By, (r). the normal coordinate on Y,, within radius r, the vector field 9 Yy.Y represented by the constant vector e. We also recognize that in normal coordinate, 7(1) - -(0)= logg y. We thus have
pasvn -1) / ( d ds (s)ds
(logg y). (xx,
I. (0n)(((5)) 80 a (s)ds
ds
- C L. [19(((s)) 1ll ||e|| d8 (s)
s <|||u|| p(x, / yk(sp(z,g))
S poly K, 0 ))
as desired, Here the penultimate line follows from Lemma 5.
F.4 PROOF OF LEMMA 19: HANDLING EXCEPTIONAL EVENTS
Proof of Lemma 19. Recall that pinx = PN PNrKN*K N- Ka and q0 = P Ry. We need to compare the kernel K with Kg, k 1 ..... N. To apply Lemma 20, we define two auxiliary kernels Kaux, K2 that are "localized" version of K and Kg. We show the auxiliary kernels close Kaux and Kg respectively in total variation, and establish bound KL(Ka.X K2)- Denote
= E M: plx, v) S 1//12), = M R.
Recall the notation in the proof of Lemma 24. To distinguish the kernels at different step, we denote 42 as the corresponding 4 at step k. Define Kank, KR by
KatM(r, dy) = KRtK (x, dyy R. (v) + K2tK (x, R2)
p(dy)
((R))
Jms Tu x, x,) ())(d)
Ka dy) = क (h, y) R. (v)
((dy))
(R:g)
By converting to normal coordinate and invoking Gaussian integration in the same way as in the proof of Lemma 24, we obtain
exp
S
Φx (()dd) & exp( 16/
(42)
h1/6
RR
When h < (/pyly K, 1), ) it is apparent (e.g., follows from Gromov's volume comparison theorem) that (R()) < 1/22, thus
2 ((RR) < H(M) 1.
We observe that Kg differs from K2 by a rejection sampling with radius 1/4, With the above bounds and the same Gaussian integration technique,
that the probability of rejection is bounded by
1
Prrreeetioon at step k) S exp(
S exp(
(43)
16h
16
31

---

Published as a conference paper at ICLR 2026
Summing up. We readily obtain
TV(Kg, Ks) exp
(44)
16h
On the other hand, by using the stopping time argument in the proof of (27), we have
K (ar R.) S exp( 16h-/6
(45)
Therefore, the following TV bound is obvious:
1
TV(K Karx) S exp
(46)
16A
Now we compute KL
KL(KA (Ka"((x,) (a )). By definition, we have
KL (L(K (x(x,))) (((z)))
/ R. Ka( dy) )
(Log Ket dy)
Ka dy) +
log
Kk (2., R.)
Ka(ar, RD)
JRg 4x (h, x, 2)
T2
We control the two terms separately.
Controlling Tj. We invoke Lemma 20 see
Kaux (ar, dy) Kx,2,, dy)
S poly(d, K, k,.-)).
Therefore, we use the elementary fact that log(1 1 + x) 2 x - 2x for x E [-1/2, to obtain
Kaux (x, dy)
KR (x, dy)
log
log
K2 (d)
Kaux (x, dy)
KR (z. dy)
Ka dy)
-1)
< 1
Katl (xx dy)
Kaux(I, dyy
K2 (xr dy)
<
poly(d, K 5 )0,
Kauk( (I, dy)
provided S (ol K, 11). Integrate with respect to Ka (z, dy) over y E R2 obtain
Ti < Ka (2, R2) Ka R₂) + poly (d, K, 8 1-
< 2 expp ( 16A + +oly (d, K, K,5-1)
I poly(d, K, 1,d. )
where the second line follows from (42) and (45), and the last line follows from h < 1/ (oly K, K,8-1) so that the exponential term sufficiently small.
Controlling Tz. This is strightforward given (45) and (42). We obtain in the same way above that
1
T2 S exp( 3// (, s/2 h²
Summarizing the above, we have shown that
KL (L(Ka"(<x,) -) (2x,)) K,
Accumulate the error over all N steps using post-processing inequality and apply Pinsker's inequal- ity, we obtain
TV aux (146) S Vpoly(d, K,81)) < VNT poly K,5
since J.N = T - 5 < T. claimed.
32

---

Published as a conference paper at ICLR 2026
G PROOF OF MAIN RESULTS
Proof of Lemma! I. This follows from combining Lemma 15 and Lemma 16.
Proof of Lemma 2. This follows from Lemma 18 and our choice of schedule IN = 7 < T.
Proof of Theorem I.
This follows from Lemma 1, Lemma and Lemma 19.
33


Reviews:
[{"review_id": "nsfQT3MYYS", "reviewer": "Reviewer_FiF9", "summary": "The paper studies Riemannian score-based diffusion models on a compact $d$-dimensional manifold $(M, g)$. Under bounded geometry (injectivity radius $\\rho > 0$; curvature/tensor bounds), an $L^2$-accurate score (Assumption 2) and a mild regularity bound on the learned score (Assumption 1. (A3)), it proves that a polynomially small step size suffices to achieve small total-variation (TV) error between the sampler's output law and an early-stopped target $p_{\\delta}$. The bound decomposes into (i) heat-flow mixing at rate $\\lambda_1$ (the spectral gap of $-\\Delta_M$), (ii) a term linear in the $L_2$ score error, and (iii) a discretization term scaling like $\\sqrt{hT}\\mathrm{poly}(d,K,\\rho^{-1},\\delta^{-1})$. Technically, the proof combines Li–Yau/Hamilton/Han–Zhang heat‑kernel/gradient estimates, a localization that freezes the drift in normal coordinates, and a Minakshisundaram–Pleijel parametrix control to bound the Brownian‑motion simulation error that arises on manifolds. The algorithm includes a \"reset to the uniform measure $\\mu$\" safeguard when a proposed step exits the injectivity radius.", "questions": "1. Reset term: Can you provide an explicit bound for\n$$\\sum^N_{k=1} \\mathrm{Pr} (\\| \\Delta_k \\| > h^{1/4}) $$ \nand show how it is absorbed into the final TV error? A short lemma that combines drift estimates from Lemma 11 with BDG would suffice.\n2. Generator sign: Which generator (and adjoint) is used in §E.2 for the parametrix? Please unify the sign and confirm that Lemma 19 and the Volterra‑series constants remain valid under the corrected operator.\n3. Geometric hypotheses: Please state (or the exact curvature condition you need) upfront, and propagate the dependence carefully through Lemmas 7/11 and Theorem 1.\n4. Heat‑kernel sup: Re‑derive the $ \\mathbb{E} [\\log^4(\\cdot)]$ control using a Li–Yau‑style $\\log \\sup H(t) \\lesssim \\frac{d}{2} \\log (1/t) + CKt$. Does the main scaling change? If not, please include the corrected calculation.\n5. Girsanov/Novikov: Please add a one‑page lemma verifying Novikov (or an equivalent condition) in your setting with cutoff. This would fully justify (5).\n6. (A3) implementability: How can one enforce (A3) in practice without access to $\\nabla \\log p_t$ ? If you intend (A3) to be a theoretical condition only, please say so and explain where a weaker data‑dependent clipping suffices.\n7. Schedule details: Please state the exact relationship among $h$, $\\omega$, $\\delta$ used to obtain the inequality $\\mathrm{Pr}(d > \\omega /3) \\leq \\omega^4$.\n8. Can you relate your work to https://proceedings.mlr.press/v251/sakamoto24a.html (https://openreview.net/forum?id=ahVFKFLYk2) ?", "limitations": "", "rating": 2, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- **Problem significance.** TV‑accurate sampling on manifolds with polynomial step size is a meaningful strengthening over prior Riemannian SGM bounds that scale poorly with $d$; cf. De Bortoli et al. 2022.\n- **Technique.** The parametrix approach to BM simulation error is original in this context and technically apt; it uses BGV-style local expansions and Volterra series to compare kernels; Berline-Getzler-Bergne (https://link.springer.com/book/10.1007/978-3-642-58088-8).\n- **Error decomposition.** Splitting the total error into mixing / score /discretization / BM simulation aligns with the Euclidean playbook but addresses the manifold-specific obstacle (non-Gaussian BM kernel).\n- **Assumptions.** Works under bounded geometry and does not require smooth/strictly positive data density; early stopping avoids small-time blow-ups.", "weaknesses": "1. **Reset-analysis mismatch.** Provide a uniform upper bound on\n$$\\sum^N_{k=1} \\mathrm{Pr} (\\| \\Delta_k \\| > h^{1/4}) $$ \nand propagate it into the final TV bound. One can combine drift bounds (from Li-Yau/Hamilton/Han-Zhang estimates) with a BDG-type tail on the martingale term to show a sub-Gaussian decay in $(\\omega^2 / (t_k - t))$; but the constants and dependence on $d$, $K$, $\\rho^{-1}$, and $\\delta^{-1}$ must be explicit.\n2. **Parametrix generator sign.** Unify the generator throughout §E.2 (either $-\\frac{1}{2}\\Delta_M + \\braket{S, \\nabla}$ or $+\\frac{1}{2}\\Delta_M + \\braket{S, \\nabla}$) ; then re-check Lemma 19 and the Volterra bounds against the exact BGV statements (e.g. Thm. 2.23, 2.26, 2.29).\n3. **Heat-kernel sup bound.** Replace $$\\sup \\log H(s/2, \\cdot,\\cdot) \\lesssim \\frac{d \\log d}{s} + Ks $$ by a standard Li–Yau form: $\\log \\sup_{x,y} H(t,x,y) \\lesssim \\frac{d}{2}\\log(1/t) + CKt$ (possibly after inserting volume factors via Bishop-Gromov/Günther). Re-derive the $\\mathbb{E}[\\log^4(\\cdot)]$ control under this bound.\n4. State geometric assumptions precisely. Add an explicit Ricci lower bound ($ \\mathrm{Ric} \\geq - K_{\\mathrm{Ric}} $) and track constant dependencies in Lemma 11 and downstream. This is standard in Hamilton's matrix Harnack and Han-Zhang's Hessian upper bounds.\n5. Girsanov conditions. Insert a lemma verifying Novikov/Kazamaki under (A3) and your cutoff—e.g., $ \\mathbb{E} \\exp( \\frac{1}{2} \\int \\| S \\|^2 dt ) < \\infty $ along the path—so that the KL steps in (5) are justified.\n6. Use the square‑root form where appropriate. In several places the proof appears to use a bound of the type $\\| \\nabla \\log p_t \\| \\lesssim (1/t + K) \\log (\\sup p_{t/2} / p_t )$, whereas standard Hamilton/Li-Yau give $\\| \\nabla \\log p_t \\|^2$ bounded by that expression. Ensure constants in, e.g., (11), remain valid with the square root in place. \n7. Assumption (A3) implementability. (A3) upper‑bounds $\\| s_{\\theta} (t_k , x) \\ |$ using the unknown true score $\\| \\nabla \\log p_{t_k} (x) \\|$. Clarify how clipping that depends only on $s_{\\theta}$ guarantees (A3) without access to the ground‑truth score; or restate (A3) as a condition on $s_{\\theta}$ alone and show it suffices for all places where (A3) is invoked.\n8. Make the $\\omega^4$ step explicit. When concluding $\\mathrm{Pr}(d(Y_y,Y_{t_k})) > \\omega /3 \\leq e^{-c\\omega^2 / (t_k - t)} \\leq \\omega^4$, spell out the exact regime relating $h$ and $\\omega$. This will help readers verify the scheduling choices that keep this probability polynomially small.\n9. Cite mixing to $\\mu$ precisely. The initialization/mixing bound uses the spectral‑gap‑driven decay of the heat semigroup—please cite Urakawa and state clearly how the \"$A$\" constant is controlled.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "k90HZPzl42", "reviewer": "Reviewer_5iY7", "summary": "The paper “Polynomial Convergence of Riemannian Diffusion Models” presents a significant advancement over the work of De Bortoli et al. (2022), which introduced generative diffusion processes on Riemannian manifolds. This framework, known as the Riemannian Score-Based Generative Model (RSGM), is motivated by the observation that many datasets, such as those arising in climate science and high-energy physics, naturally lie on Riemannian-like manifolds. De Bortoli et al. demonstrated that, in such settings, performing diffusion directly on the underlying manifold can substantially improve generative performance. However, their approach suffered from a severe limitation: achieving a small divergence between the model and target distributions required an exponentially small discretization step size, leading to prohibitive computational costs. The current paper overcomes this issue by proving that, under mild assumptions on the manifold structure, one can instead use a polynomially small step size to achieve comparable accuracy. This result dramatically improves the practical feasibility of RSGMs and paves the way for their broader application in real-world generative modeling tasks.", "questions": "I will now proceed to ask questions and raise issues that I would like the authors to address in order to both finalize my personal judgment over the paper and improve the manuscript. \n\n**Questions:**\n\n1. Why do I re-start from the uniform distribution if the diffusion step from $Y_k$ exits the injective radius ? What happens if one restarts from $Y_k$ and samples new noise? How often would the reset take place in practical applications? \n2. In the continuous-time variance-exploding prescription, a common choice for the prior distribution in backward diffusion is a wide Gaussian distribution, with zero mean and variance equal to the horizon $T$. Are there other prior distributions, alternative to the uniform one, that can reduce the initialization error? Can Urakawa (2006) be generalized to other initializations?\n3. Recent studies in this field, such as Ventura et al. (2025) \"Manifolds, Random Matrices and Spectral Gaps: The geometric phases of generative diffusion\", Wang et al. (2025) \"Diffusion Models Learn Low-Dimensional Distributions via Subspace Clustering\", Stanczuck et al. (2024) \"Diffusion Models Encode the Intrinsic Dimension of Data Manifolds\" show that, when data are supported by a low-dimensional manifold, the geometry of this manifold  progressively emerges in the backward diffusion in the Euclidean space, without the need of constraining the trajectory on the manifold itself. I would ask to the authors whether it would be possible to study the time dependence (that I guess translates into studying the $\\delta$ dependence) of the TV and see whether its upper bound becomes tighter at smaller times, when the geometry of the target manifold already emerges even when diffusing in the Euclidean space. \n\n**Readability issues:**\n\n1. [101] What does BM stand for?\n2. [157-158] Here and across all the paper there is a conflict between the dimension of the ambient Euclidean space and the dimension of the manifold used in the subsequent proofs. \n3. [270-271] You should probably state that M has dimension d, since it appears implicitly in A3.  \n4. Both in [306-307] and the very beginning in [088] the authors may introduce a bit better the role of the error constant 𝜀, as a quantity employed to tune the horizon time $T$ and from that the discretization stepsize. \n5. [395-396] May be useful to uniform the nomenclature for ℐ with respect to the previous expression in [394-395].\n\n**Typos**: \n1. [073] “exponential”\n2. [362] “time frozen”\n3. [324-325] \"poly-logarithmic\"\n\nI take advantage for thanking the authors for the nice job and I wait for their response.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The work has several notable strengths:\n\n1. The presentation is clear and pedagogical, making the material accessible even to readers without a strong background in differential geometry.\n\n2. The discussion of related work is thorough and well-integrated, situating the paper within a coherent and well-defined research landscape.\n\n3. The analytical results appear sound, and the way of deriving the bound over the total variation, as a measure of the distance between two distributions, could be employed in characterizing other contexts, such as classifier-free-guidance and other diffusion routines where an exact theory is lacking and sources of errors are similar.", "weaknesses": "The results are primarily analytical and the central goal of the work is to derive new performance bounds for RSGMs. I hence appreciate the total analytical vocation of the manuscript. On the other hand, one might think that the total absence of experimental validation could appear as a weakness of the work. De Bortoli et al. (2022) compare RSGM with other manifold-based diffusive routines on manifold supported datasets, like climate science spherical data. They quantify the quality of the performance in terms of the log-likelihood achieved by the model. I was wondering whether the authors of the current paper ever evaluated the idea of running similar experiments with both RSGM and other procedures (e.g. Mixture of Kent, Moser Flow etc.) and show that RSGM can reach a comparable degree of the performance still employing a polynomial number of steps. Such scaling could be analyzed as a function of relevant control parameters, including manifold size and early stopping time. Another possible, and possibly equivalent, experiment would instead imply verifying that Riemann supported real datasets verify the assumptions of the current analysis.  \n\nIn addition to addressing the questions outlined in the following section of this review, I would appreciate if the authors could argue around this point more in detail.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "FdRZVHgmtg", "reviewer": "Reviewer_qEnZ", "summary": "This paper extends the Riemannian diffusion model, studying the convergence behavior under the L2 accurate score estimate. The paper shows that a polynomially small stepsize is sufficient to guarantee small sampling error in TV distance with only mild curvature assumptions. For technical contribution, the proof uses Li-Yau gradient bounds for the heat kernel, a localization scheme for the drift, and Minakshisundaram–Pleijel parametrix for controlling deviations from the heat flow and its discretized proxy.", "questions": "Please address the questions raised in the weakness section", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The paper presents the first polynomial TV convergence results for Riemannian diffusion models to the best of my knowledge. While it extends the theoretical results in Riemannian diffusion models, the paper proves TV convergence instead of Wasserstein convergence under mild geometric conditions. \n\n- The paper provides analytical techniques to prove the polynomial convergence: Li-Yau estimate for the heat kernel, localization of the drift fields, and parametrix estimates. These tools may be meaningful for bridging diffusion theory and geometry and could be used for developing a sampling method for Riemannian diffusion models.", "weaknesses": "- I'm unsure what the impact of the paper is and how it contributes to the field. While the paper gives a theoretical guarantee for the Riemannian diffusion model, the empirical results from De Bortoli et al already validated that it worked in diverse settings. I agree that the paper provided new tools to use for future study in this field, but it is unclear how important it is to show polynomial TV convergence for the already working model. If the paper presented a recipe for stepsize (with empirical results) or an advanced sampling method, it should have been a strong signal for acceptance. I would appreciate it if the authors could give an explanation of the importance of the work.\n\n- While the theoretical contributions are important, I believe empirical results, at least on small toy settings, are necessary. I cannot say that the details of the proof are all correct due to an insufficient understanding of differential geometry. Empirical results should have complemented this, giving strong evidence.\n\n- Is the result only applicable to the Riemannian diffusion model from De Bortoli's design? For example, there are several designs of diffusion models, such as [Chin-Wei Huang et al., 2022], [Lou et al., 2023], Riemannian Flow Matching [Chen and Lipman, 2024], or Riemannian Diffusion Mixture [Jo et al., 2024], and I wonder if the convergence results or its argument could be applied to different designs.\n\n[Chin-Wei Huang et al., 2022] Riemannian Diffusion Models, NeurIPS 2022\n[Lou et al., 2023] Scaling Riemannian Diffusion Models, NeurIPS 2023\n[Chen and Lipman] Flow Matching on General Geometries, ICLR 2024\n[Jo et al., 2024] Generative Modeling on Manifolds Through Mixture of Riemannian Diffusion Processes, ICML 2024", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)