High-accuracy and dimension-free sampling with diffusions

Summary

This paper studies the iteration complexity of diffusion-based samplers as a function of target accuracy $\varepsilon$. Prior analyses of probability-flow-ODE and DDPM-style samplers achieve iteration counts scaling polynomially in $1/\varepsilon$; the authors' contribution is the first diffusion-based sampler whose iteration count scales polylogarithmically in $1/\varepsilon$ while only assuming access to (Lipschitz, sub-exponentially accurate) score estimates. The target class is "bounded plus noise" — distributions $q = q_{\mathrm{pre}} \star \mathcal{N}(0,\sigma^2 I)$ with $q_{\mathrm{pre}}$ supported in a ball of radius $R$ — which naturally captures Gaussian mixtures and models what practical diffusion pipelines target through early stopping.

The conceptual move is to view the drift of the reverse ODE, along the trajectory, as a function of time that is well-approximated by a low-degree polynomial. The authors prove a dimension-free bound on the $k$-th time derivative of the score along the reverse ODE — depending only on $R/\sigma_t$ up to $(kp)^{1/2}$ moment terms — and feed this into the collocation / Picard-iteration ODE solver of Lee et al. (2018) applied on short windows. Because collocation converges exponentially per Picard step under a bounded polynomial basis, and low-degree Taylor approximation of the drift gives approximation error controlled by the derivative bound rather than by the discretization step size, the overall iteration budget is $\tilde O((R/\sigma)^2 \cdot \mathrm{polylog}(1/\varepsilon))$, and the dimension $d$ enters only through $R$ (and hidden log factors). A short underdamped-Langevin corrector step is then used to upgrade the $W_2$ guarantee of Theorem 3.7 to the TV guarantee of Corollary 3.9.

The main conceptual limitation is scope. The high-accuracy rate does not come "for free": it requires (a) sub-exponentially small pointwise score error rather than the standard $L^2$ assumption, (b) Lipschitz score estimates, and (c) the target being a genuine convolution with Gaussian noise ($\sigma > 0$), with the admissible tail scale $\varepsilon_{\mathrm{err}}$ tightening as $(R/\sigma)^{-2/3}$. The final TV guarantee is reached only after post-processing with underdamped Langevin, a step whose necessity is essentially unmotivated in the main body. And $R$ can scale with $d$ for many natural distributions, so "dimension-free" describes the parameterization rather than always the concrete regime. There are also no experiments and no code, so the polylog scaling and the sensitivity of the pipeline to each assumption are not empirically probed.

Strengths

• Conceptual contribution. First diffusion-based sampler with polylog($1/\varepsilon$) iteration complexity assuming only score access — a genuine acceleration over the polynomial-in-$1/\varepsilon$ regime that prior discretization analyses were stuck in, even under perfect scores.
• Structural result on time derivatives. Lemma 3.3's dimension-free bound on higher time derivatives of the reverse-ODE drift is novel and load-bearing. The $(74k/\sigma_t^2)^k (R/\sigma_t + (kp)^{1/2})^{2k}$ form is exactly what makes low-degree Taylor approximation of the drift viable.
• Robustness to perturbed initialization. The coupling machinery of Lemma 3.4 and Theorem 3.5 shows the derivative bound survives perturbation of the initial marginal in $W_2$, which is what lets errors reset across windows rather than compound exponentially.
• $W_2$-to-TV via short Langevin corrector. Reusing the short-time smoothing idea of Gupta et al. (2024a) and Chen et al. (2024c) to avoid Metropolis adjustment gives a TV guarantee at essentially no additional iteration cost, showing that the ODE-side of the argument is tight.
• Honest positioning against related work. The comparison to Huang et al. (2024a), Wainwright (2025), Li et al. (2025a), and higher-order solvers (Huang et al. 2025a,b; Li et al. 2025b) carefully identifies the axes along which each result is or is not comparable, including the key point that prior higher-order accelerations only give polynomial gains under the same access model.

Weaknesses

• Score-error assumption is materially stronger than standard. Assumption 2 requires sub-exponentially small pointwise error, not just an $L^2(q_t)$ bound. The paper flags this and lists it as open, but does not argue whether it is a proof artefact or fundamental to the polylog rate, nor whether trained neural scores plausibly satisfy it. Since the whole high-accuracy claim rests on this, more discussion is warranted.
• "Dimension-free" is a delicate framing. The iteration count depends on $R$ (Assumption 1's ball radius), and $R$ can be $\Theta(\sqrt d)$ or worse for many natural distributions. The Gaussian-mixture regime $R = \Theta(\sqrt{\log k})$, $\sigma = 1$ illustrates a strict improvement, but the abstract's headline reads more strongly than the actual guarantee, and Remark 3.8's hidden $\log d$ caveat is easy to miss.
• $\sigma > 0$ is fundamental and the trade-off is left implicit. The rate degrades as $(R/\sigma)^2$ and the admissible $\varepsilon_{\mathrm{err}}$ shrinks as $(R/\sigma)^{-2/3}$. The practical regime where diffusion models are used (early-stopped diffusion with small $\sigma$) is not the regime of the sharpest guarantee; the outlook mentions this but only briefly.
• Underdamped-Langevin corrector arrives with no signposting. Section 3.3 announces that the TV upgrade "has to run underdamped Langevin Monte Carlo" although nothing in Section 2 or Algorithm 2 mentions it. A reader following the body has no way to know the pipeline is a two-stage ODE + Langevin scheme rather than pure collocation until Appendix H.
• No experiments and no code. For a paper whose central quantitative claim is a scaling law in $1/\varepsilon$ and $R/\sigma$, there is no numerical illustration — not even on a toy Gaussian mixture — that the polylog rate manifests, that the constants absorbed by $\tilde O$ are practical, that Picard depth $m = O(\log(1/\varepsilon))$ matches theory, or that the corrector step behaves as promised.
• Lipschitz-score assumption without diagnostics. Assumption 3 is defended as "the true score is Lipschitz," but score estimators — particularly deep networks — can have effective Lipschitz constants much larger than $\nabla\log q_t$, and $\tilde L$ enters the Picard depth as $\ln(\tilde L)$; how loose $\tilde L$ can be before the polylog rate degrades is not discussed.
• Notational and bookkeeping issues in the appendix. Several concrete slips make the proof hard to check on a first read: Assumption 2's tail bound is quantified over the wrong variable ("for all $y \geq 0$" where $z$ appears), Definition 2.1 bounds an integral by "$\gamma_\phi t$" using an undefined $t$, the scaling of $X_t$ differs between Lemma 3.1 ($e^{-(T-t)} \bar X$, $\bar X \sim q_{\mathrm{pre}}$) and the appendix notation ($\sqrt{1-\sigma_t^2}\,\bar X$, $\bar X \sim q$), $F_t$ vs $F_t^\ast$ are used heavily but $F_t^\ast$ is never defined, and the failure-probability chain from Lemma E.6 ($1 - 2\varepsilon_1$) to Theorem 3.5 ($1 - \delta^2$) requires implicitly setting $\varepsilon_1 = \delta^2/2$ rather than the $\varepsilon_1 = \delta$ used elsewhere.
• Lipschitz constant $L$ in Algorithm 3 is circularly specified. The corrector step sets $\eta$ in terms of $L$ and bounds $L$ by $\mathrm{poly}(d/\eta)$, leaving the actual numerical setting to an unresolved fixed point. Also, "the Lipschitz constant of $q$" should presumably be of $\nabla\log q$.
• Aggressive characterization of prior work. The claim that Huang et al. (2025a,b) and Li et al. (2025b) cannot reach polylog($1/\varepsilon$) rests on "implicitly assum[ing] that the number of iterations is at least exponential in $K$" attributed to a discussion in Li et al. (2025b). This is a substantive comparative claim used to distinguish the contribution; readers cannot verify the "at best polynomial" conclusion from what is stated in-paper.
• Chebyshev-basis presentation. The rendered formula for $\tilde\phi_j(x)$ is unparseable, and the accompanying statement "$\tilde\phi_j$ is a degree-$j$ polynomial satisfying $\tilde\phi_j(c_i) = \mathbf{1}[i=j]$" is dimensionally inconsistent (a degree-$j$ polynomial cannot vanish at $D-1$ prescribed nodes unless $j \geq D-1$).

Reproducibility & code

• No implementation released. Algorithms 1 (PICARD), 2 (COLLOCATIONDIFFUSION), and 3 (CORRECTEDCOLLOCATIONDIFFUSION) are specified only in pseudocode; the accompanying code artefact is empty. For a paper whose headline is a scaling law in $1/\varepsilon$ and $R/\sigma$, even a small reference implementation on a synthetic Gaussian mixture would materially strengthen confidence that the polylog rate manifests and that the constants absorbed in $\tilde O$ are not enormous.
• Denoising schedule is not tabulated. Algorithm 2 takes as input a schedule $0 < t_1 < \dots < t_n = T$, but the schedule itself is derived only implicitly in the proof of Theorem G.1 (windows of length $\tilde O(\sigma_t^2/R^2)$). An implementer must reconstruct it from the appendix, which is a genuine barrier.
• Constants inside $\tilde O$ are not made numerical. The step size $h = \gamma \cdot k / ((1+(R/\sigma)^2)(\ln(1/\varepsilon_1)+\ln(d)))$ uses a "small enough constant $\gamma$", and the Chebyshev-basis bound rests on the constant $\leq 2000$ imported from Lee & Vempala (2017). None of these are consolidated into a table an implementer could plug in.
• Corrector-step parameters are circularly specified. The friction and step size for underdamped Langevin depend on $L$, an upper bound on the Lipschitz constant, which is itself only bounded as $L \leq \mathrm{poly}(d/\eta)$ with $\eta$ in turn depending on $L$.
• $\varepsilon_{\mathrm{err}}$ and $\tilde L$ are used as inputs but not diagnosed. Corollary 3.9 assumes $\varepsilon_{\mathrm{err}} \leq \tilde O(\varepsilon^{2/3} L^{1/2} d^{1/6}/(R/\sigma)^{2/3})$, but there is no procedure for estimating $\varepsilon_{\mathrm{err}}$ (or $\tilde L$) from a given score model. This is a real gap between the theoretical guarantee and any actual pipeline.
• Support radius $R$ is used as an input but not estimated. Since the iteration count scales as $(R/\sigma)^2$, $R$ is load-bearing, but in practice $R$ is not directly computable from samples for a real target.

Recommended Changes

Essential

• Release reference code and at least one numerical experiment. Implement Algorithms 2 and 3 for a Gaussian mixture (e.g. $R = \Theta(\sqrt{\log k})$, $\sigma = 1$) and empirically show iteration count scaling as polylog($1/\varepsilon$) and near-independence of $d$ at fixed $R/\sigma$. Addresses the No implementation released and No experiments and no code items.
• Introduce the underdamped-Langevin corrector in the main body. Add a paragraph in Section 2.4 (or the Algorithm 2 block) noting that the TV guarantee of Corollary 3.9 requires the corrector step defined in Algorithm 3. Addresses Underdamped-Langevin corrector arrives with no signposting.
• Fix the failure-probability accounting between Theorem 3.5 and Lemma E.6. Make the choice of $\varepsilon_1$ explicit throughout the induction so that the promised $1 - \delta^2$ bound in Theorem 3.5 actually follows from the $1 - 2\varepsilon_1$ conclusion of Lemma E.6. Addresses the corresponding item in Notational and bookkeeping issues.
• Correct Assumption 2's quantifier. Replace "for all $y \geq 0$" with "for all $z \geq 0$" to match the variable in the tail bound, and clarify whether $\varepsilon_{\mathrm{err}}$ is the $\psi_1$ Orlicz norm or the tail-scale parameter. Addresses Score-error assumption is materially stronger than standard and the bookkeeping item.
• Reconcile the two definitions of $X_t$. Lemma 3.1 uses $X_t = e^{-(T-t)} \bar X$ with $\bar X \sim q_{\mathrm{pre}}$, while Appendix B uses $X_t = \sqrt{1-\sigma_t^2}\,\bar X$ with $\bar X \sim q$. Add one line stating that these coincide under $\sigma_t^2 = 1-e^{-2(T-t)}$ (and match base distributions), or otherwise unify the notation.
• Define $F_t$ and $F_t^\ast$ explicitly. Appendix B introduces $F_t$ as the score estimate but never defines $F_t^\ast$; state which is which and how each relates to $\nabla\log q_t$ and to the ODE drift $y + \nabla\log q_t$.
• Fix the Lipschitz-constant specification in Algorithm 3. Break the circularity between $\eta$ and $L$ by stating either an a-priori upper bound on $L$ (e.g. via bounded second derivative of the log-density) or how to run the algorithm without knowing $L$ in advance. Also clarify whether $L$ refers to $q$ or to $\nabla\log q$.

Suggested

• Consolidate constants into a table. Add an "implementer's summary" listing the Chebyshev-basis constant ($\leq 2000$ from Lee & Vempala), the small constant $\gamma$ in the step-size formula, the Picard-depth prefactor $m$, and the admissible range of $\varepsilon_{\mathrm{err}}$ under Corollary 3.9. Addresses Constants inside $\tilde O$ are not made numerical.
• Make the denoising schedule explicit. Provide either the closed-form schedule or pseudocode for constructing it from $R, \sigma, \varepsilon$ alongside Algorithm 2. Addresses Denoising schedule is not tabulated.
• Motivate the sub-exponential score-error assumption. Include a short discussion in Section 2.1 explaining whether Assumption 2 is a proof artefact or fundamental, ideally connecting it to known concentration properties of trained score networks. Addresses Score-error assumption is materially stronger than standard.
• Discuss the practical role of $\tilde L$. Note how $\ln(\tilde L)$ enters $m$ in Theorem G.1 and comment on the regime where a trained score has looser Lipschitz constant than the true score. Addresses Lipschitz-score assumption without diagnostics.
• Correct Definition 2.1's variable typo. The bound reads "$\leq \gamma_\phi t$" where $t$ is undeclared; presumably this should be "$\leq \gamma_\phi x$" or a constant. Addresses the corresponding item in Notational and bookkeeping issues.
• Restore the Chebyshev-basis formula and reconcile the degree claim. Supply the intended closed form $\tilde\phi_j(x) = \sqrt{1-c_j^2}\,T_D(x)/(D(x-c_j))$ and either state that this is a common degree ($D-1$) polynomial family or clarify the upper-triangular structure. Addresses Chebyshev-basis presentation.
• Address the "34 new terms" and $(74k)^k$ bookkeeping. Enumerate the term types produced when differentiating Eq. (11), and rewrite the multiplicative accumulation as an explicit induction. Addresses the "34 new terms" and $(74k)^k$ items in Notational and bookkeeping issues.
• Rewrite the Section G intro to attribute the TV bound to Appendix H. Make clear that Theorem G.1 delivers $W_2$ and that the TV guarantee comes only after the corrector step. Addresses the corresponding item in Underdamped-Langevin corrector arrives with no signposting.
• Provide a self-contained sketch for the "$K$-dependence" claim. Add one line showing why setting $K = \Theta(\log(1/\varepsilon))$ in Huang et al. / Li et al. blows up the iteration bound to $\mathrm{poly}(1/\varepsilon)$. Addresses Aggressive characterization of prior work.
• Fix the citation duplication. "Li et al. (2024a); Wu et al. (2024); Li et al. (2024a)" repeats a reference; almost certainly one should be Li et al. (2023).
• Add a train.py-style diagnostic for $\varepsilon_{\mathrm{err}}$. Even a short recipe for empirically bounding the sub-exponential tail scale of a trained score would make the corrector-step tolerance concretely usable.