High-accuracy and dimension-free sampling with diffusions

Abstract: Diffusion models have shown remarkable empirical success in sampling from rich multi-modal distributions. Their inference relies on numerically solving a certain differential equation. This differential equation cannot be solved in closed form, and its resolution via discretization typically requires many small iterations to produce \emph{high-quality} samples. More precisely, prior works have shown that the iteration complexity of discretization methods for diffusion models scales polynomially in the ambient dimension and the inverse accuracy $1/\varepsilon$. In this work, we propose a new solver for diffusion models relying on a subtle interplay between low-degree approximation and the collocation method (Lee, Song, Vempala 2018), and we prove that its iteration complexity scales \emph{polylogarithmically} in $1/\varepsilon$, yielding the first ``high-accuracy'' guarantee for a diffusion-based sampler that only uses (approximate) access to the scores of the data distribution. In addition, our bound does not depend explicitly on the ambient dimension; more precisely, the dimension affects the complexity of our solver through the \emph{effective radius} of the support of the target distribution only.

• The paper presents the first dimension-free sampler achieving polylogarithmic iteration complexity in 1/ε for high-accuracy diffusion-based sampling.
• It leverages low-degree polynomial approximation along the ODE path via a collocation method to bypass conventional dimension-dependent discretization.
• The approach controls score estimation errors with sub-exponential tails, ensuring robustness even with complex multimodal distributions.

High-Accuracy, Dimension-Free Diffusion-Based Sampling

Introduction and Motivation

Diffusion models have established themselves as effective frameworks for sampling from complex, multimodal distributions, including high-dimensional, non-log-concave targets. Sampling with diffusion models proceeds by simulating a reverse-time SDE or an associated ODE using a learned score function, typically estimated from data. Practical implementations discretize these continuous processes, resulting in bias and requiring numerous steps for precise sampling. Prior theoretical analyses have demonstrated that the number of required discretization steps, or iteration complexity, typically scales polynomially with both the ambient dimension $d$ and the inverse target accuracy $1/\epsilon$.

The problem of high-accuracy sampling—where the goal is to achieve iteration complexity scaling only polylogarithmically with $1/\epsilon$—remains largely open for general (non-log-concave) distributions and under the standard model of score estimate access. In parallel, recent applications demand algorithms whose performance is not explicitly dictated by the ambient dimension, motivating dimension-free or geometry-adaptive analyses.

Main Contributions

This work presents the first dimension-free, high-accuracy sampler for a large class of distributions with the following properties:

• Polylogarithmic iteration complexity in $1/\epsilon$: The number of iterations required to obtain samples within $\epsilon$ total variation (or Wasserstein-2) distance of the target is $\tilde{O}(\mathrm{polylog} (1/\epsilon))$ (where $\tilde{O}$ suppresses polylogarithmic factors).
• No explicit dependence on ambient dimension $d$: Complexity depends only on the effective support radius of the distribution, specifically scaling with $(R/\sigma)^2$ for a target that is a convolution of a compactly supported distribution (disk of radius $R$) and isotropic Gaussian noise of variance $\sigma^2$.
• No reliance on MCMC accept/reject mechanisms: Unlike prior high-accuracy log-concave samplers (e.g., MALA), no Metropolis-Hastings correction or density ratio access is assumed; only (possibly approximate, sub-exponentially tailed) access to the score is required.

The core technical insight is the construction of an ODE solver that leverages low-degree polynomial approximation of the score function along the ODE path. By using a collocation (Picard iteration) scheme with polynomial basis functions, the scheme sidesteps the adverse dimension dependence traditionally introduced by naive discretization.

Conceptual and Technical Framework

Problem Setting

The target $q$ is assumed to be the convolution of a distribution $\nu$ (supported within a ball of radius $R$) and a Gaussian $\mathcal{N}(0, \sigma^2 I)$. This includes mixtures of well-separated isotropic Gaussians as a special and important case. Sampling from $q$ is cast as simulating its reverse-time diffusion, i.e., the probability flow ODE:

$\frac{dy_t}{dt} = y_t + \nabla \log q_t(y_t)$

where $q_t$ is the time-evolved version of $q$ under Brownian motion.

Avoiding Dimension-Dependent Discretization

Traditional Euler-type discretizations incur errors growing as $O(h^2 \sqrt{d})$, forcing step sizes $h = O(d^{-1/4})$ to control bias. This work circumvents the issue by proving that, for the target class considered, the vector field defined by $y_t + \nabla \log q_t(y_t)$ can be closely approximated by a low-degree polynomial in time along solution trajectories. This is made precise via bounds on high-order time derivatives of the score function, which are shown to be controlled in terms of $R$, $\sigma$, and $k$ (the degree).

Collocation and Exponential Convergence

Using a tailored version of the collocation method—specifically, Picard iteration with polynomial interpolation at Chebyshev nodes—the authors implement a solver that converges exponentially (in the number of iterations) to the true ODE solution over appropriately chosen time-windows. The time step size is determined by the effective support radius $R$, noise level $\sigma$, and polynomial approximation error—not by the dimension.

This guarantees that, with each window, the sampling trajectory remains in close Wasserstein proximity to the ideal reverse process. Chaining such high-accuracy, locally-exact solvers yields the overall $\tilde{O}((R/\sigma)^2 \polylog (1/\epsilon))$ complexity.

Robustness to Score Estimation

The main convergence result requires that the score estimates are not only Lipschitz but also have sub-exponential error tails (rather than merely $L_2$ accuracy as in standard analyses). Stronger error control is necessary since the method contracts error at an exponential rate, and polynomial tail error would otherwise dominate.

Key Theoretical Results and Claims

• Polylogarithmic Total Variation Convergence: There exists an algorithm that, given a target as above, outputs samples whose law is within $\epsilon$ total variation distance of $q$ using only $\tilde{O}((R/\sigma)^2 \polylog (1/\epsilon))$ evaluations of the score oracle (Theorem 1; Corollary).
• Independence from ambient dimension: The iteration complexity depends on $R/\sigma$ but not explicitly on $d$. For natural settings such as Gaussian mixtures with sub-logarithmic separation between centers, this yields much lower complexity compared to all previous diffusion-based methods.
• Robustness to initialization and error accumulation: Analysis accommodates non-ideal initialization (sampling noise) and demonstrates control on accumulation of local approximation errors across time windows, exploiting exponential contraction at each stage.
• Generality beyond Gaussian mixtures: Applies to all bounded-plus-noise targets, significantly generalizing previously dimension-free analyses that were specific to mixtures of isotropic Gaussians.

Numerical Rates and Scaling

The main theoretical guarantee asserts that, for error $\epsilon$ and $(R/\sigma)$ moderate, the sample complexity is $\tilde{O}((R/\sigma)^2 \polylog (1/\epsilon))$—an exponential improvement in $\epsilon$ dependence over previous analyses, which were at best polynomial. Notably, even under perfect score estimation, existing discretized samplers fall short of this scaling.

Relation to Prior and Contemporary Work

The theoretical literature on diffusion model convergence (e.g., [lee2023convergence], [chen2023sampling]) has focused on establishing polynomial rates in $d$ and $1/\epsilon$ under various smoothness, tail, or data geometry assumptions. Only very recently have a handful of works provided strong polynomial acceleration in $1/\epsilon$ via higher-order numerical solvers ([huang2025convergence], [wu2024stochastic], [li2025faster]), but with assumptions or preconditions (e.g., requiring access to log density ratios, or exponential in the polynomial order $K$) that fundamentally limit their practical high-accuracy performance.

The present work distinguishes itself not only by achieving exponential acceleration with only score access, but by analyzing a broad and practically motivated class of distributions, thereby situating itself as a canonical result for high-accuracy sampling in generative diffusion models.

Implications and Future Directions

Practical Implications

• Efficient high-fidelity generation: The method enables, in principle, sampling from complex high-dimensional structured distributions using orders of magnitude fewer score function queries at high precision.
• Applicability to challenging geometries: Situations where the data lie in intrinsically low-dimensional or bounded subspaces are covered, permitting practical sampling for distributions ill-suited to generic log-concave methods.
• Theoretical readiness for next-generation diffusion samplers: The advances here lay the groundwork for developing efficient, theoretically grounded ODE solvers for real-world diffusion-based generation tasks (e.g., in vision, speech, or scientific simulation).

Theoretical Implications

• Bridges gap with log-concave literature: The result mirrors properties familiar in log-concave sampling (e.g., dimension-free Metropolis-adjusted Langevin MCMC) but under the diffusion model and general multimodal structure.
• Connections to polynomial approximation in SDE/ODE solutions: The low-degree approximation property suggests deeper links between the pathwise regularity of diffusion/semi-group flows and efficient sample generation.
• Potential for further relaxation of assumptions: Open directions include weakening the strong subexponential tail requirement on the score error, and replacing compact support with sufficiently strong moment conditions.

Future Work

Key questions left open include:

• Can high-accuracy, dimension-free results be extended to more general or heavy-tailed distributions, or to settings with weaker error guarantees on the score function?
• Does similar exponential contraction hold for stochastic diffusions (e.g., DDPM) as opposed to deterministic ODE flows?
• What are the minimal sufficient access models (e.g., can log-density access be dispensed with for all high-accuracy, exponentially-accelerated samplers)?

Conclusion

This work provides the first theoretical construction of a diffusion-based sampler with polylogarithmic accuracy scaling and explicit independence from the ambient dimension, under score estimation models accessible with current techniques. The core strategy involves bounding high-order time derivatives of the score function, enabling low-degree polynomial pathwise approximation and exploitation of the collocation method for ODE resolution—an approach novel to the generative diffusion literature. The theoretical and practical promise of this strategy establish it as a standard reference point for future high-accuracy diffusion-based generative modeling research.

Citations:

• "High-accuracy and dimension-free sampling with diffusions" (2601.10708)
• Related works: [lee2023convergence], [chen2023sampling], [huang2025convergence], [li2025faster], [wu2024stochastic]