category: literaturenote citekey: chengeometricregularitydeterministicsampling2025 title: Geometric Regularity in Deterministic Sampling of Diffusion-based Generative Models authors: "Chen, Defang; Zhou, Zhenyu; Wang, Can; Lyu, Siwei" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.10177v1" zotero_key: 3GPW6HF6 zotero_storage: FT3BAMET collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Chen, Defang"

status: converted

Geometric Regularity in Deterministic Sampling of Diffusion-based Generative Models

Defang Chen defangch@buffalo.edu

Department of Computer Science and Engineering University at Buffalo, State University of New York Buffalo, NY 14260, USA

Zhenyu Zhou zhyzhou@zju.edu.cn

College of Computer Science Zhejiang University Hangzhou, 310027, China

Can Wang wcan@zju.edu.cn

College of Computer Science Zhejiang University Hangzhou, 310027, China

Siwei Lyu^∗

siweilyu@buffalo.edu

Department of Computer Science and Engineering University at Buffalo, State University of New York Buffalo, NY 14260, USA

Abstract

Diffusion-based generative models employ stochastic differential equations (SDEs) and their equivalent probability flow ordinary differential equations (ODEs) to establish a smooth transformation between complex high-dimensional data distributions and tractable prior distributions. In this paper, we reveal a striking geometric regularity in the deterministic sampling dynamics: each simulated sampling trajectory lies within an extremely lowdimensional subspace, and all trajectories exhibit an almost identical "boomerang" shape, regardless of the model architecture, applied conditions, or generated content. We characterize several intriguing properties of these trajectories, particularly under closed-form solutions based on kernel-estimated data modeling. We also demonstrate a practical application of the discovered trajectory regularity by proposing a dynamic programming-based scheme to better align the sampling time schedule with the underlying trajectory structure. This simple strategy requires minimal modification to existing ODE-based numerical solvers, incurs negligible computational overhead, and achieves superior image generation performance, especially in regions with only 5 ∼ 10 function evaluations.

Keywords: Diffusion-based generative models, sampling dynamics, trajectory regularity, low-dimensional structure, dynamic programming

1 Introduction

Diffusion-based generative models (Sohl-Dickstein et al., 2015; Song and Ermon, 2019; Ho et al., 2020; Song et al., 2021c; Karras et al., 2022; Chen et al., 2024) have gained significant attention and achieved remarkable results in image (Dhariwal and Nichol, 2021;

©XXXX Defang Chen, Zhenyu Zhou, Can Wang and Siwei Lyu.

License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/.

^∗. Corresponding Author.

Rombach et al., 2022), audio (Kong et al., 2021; Huang et al., 2023), video (Ho et al., 2022; Blattmann et al., 2023), and notably in text-to-image synthesis (Saharia et al., 2022; Ruiz et al., 2023; Podell et al., 2024; Esser et al., 2024). These models introduce noise into data through a forward process and subsequently generate data by sampling via a backward process. Both processes are characterized and modeled using stochastic differential equations (SDEs) (Song et al., 2021c). In diffusion-based generative models, the pivotal element is the score function, defined as the gradient of the log data density w.r.t. the input (Hyv¨arinen, 2005; Lyu, 2009; Raphan and Simoncelli, 2011; Vincent, 2011), irrespective of specific model configurations. Training such a model involves learning the score function, which can be equivalently achieved by training a noise-dependent denoising model to minimize the mean squared error in data reconstruction, using the data-noise pairings generated during the forward process (Karras et al., 2022; Chen et al., 2024). To synthesize new data, diffusionbased generative models solve the acquired score-based backward SDE through a numerical solver. Recent research has shown that the backward SDE can be effectively replaced by an equivalent probability flow ordinary differential equation (PF-ODE), preserving identical marginal distributions (Song et al., 2021c ,a; Lu et al., 2022a; Zhang and Chen, 2023; Zhou et al., 2024a). This deterministic ODE-based generation reduces the need for stochastic sampling to just the randomness in the initial sample selection, thereby simplifying and granting more control over the entire generative process (Song et al., 2021a,c). Under the PF-ODE formulation, starting from white Gaussian noise, the sampling trajectory is formed by running a numerical solver with discretized time steps. These steps collectively constitute the time schedule used in sampling.

Despite the impressive generative capabilities exhibited by diffusion-based models, many mathematical and statistical aspects of these models remain veiled in mystery, largely due to the complex nature of SDEs and neural network models, and the substantial data dimensions involved. In this paper, we reveal an intriguing regularity in the deterministic sampling dynamics of diffusion models, i.e., the tendency of sample paths to exhibit a regular "boomerang" shape, as illustrated in Figure 1. More precisely, we observe that each sampling trajectory barely strays from the displacement vector connecting its beginning and end points (Section 3.1), and the trajectory deviation can be effectively captured using two orthogonal bases (Section 3.2). Therefore, the sampling trajectory in the original high-dimensional data space can be faithfully approximated by its projection onto a three-dimensional subspace. These projected spatial curves are fully characterized by the Frenet-Serret formulas (Do Carmo, 2016) and consistently exhibits a strong geometric shape pattern, irrespective of initial random samples, applied control signals, or the corresponding real data samples (Figure 5, Section 3.3). This simple structure guarantees the common heuristic use of a shared time schedule for synthesizing different samples and large sampling steps without incurring much truncation error (Song et al., 2021a; Karras et al., 2022; Lu et al., 2022a), especially at the initial step (Dockhorn et al., 2022; Zhou et al., 2024a).

We hypothesize that the geometric regularity reflects some underlying structures of sampling trajectories. This work attempts to elucidate this phenomenon. We start by simplifying the ODE-based sampling trajectory (Section 2.4), which reveals an implicit denoising trajectory controlling the direction of the associated sampling trajectory and encapsulating its curvature information (Section 4.1). Built upon this insight, we show that the denoising trajectory admits a closed-form solution when considering kernel density estimates (KDEs)

Figure 1: A geometric picture of deterministic sampling dynamics in diffusion models. Each initial sample (from the noise distribution) starts from a big sphere and converges to the final sample (in the data manifold) along a regular sampling trajectory. The score direction points to the denoising output of the current position, and the denoising output forms an implicit denoising trajectory controlling the explicit sampling trajectory. Every sampling trajectory inherently lies in a low-dimensional subspace with almost the same shape.

of varying bandwidths to approximate the original data distribution from training samples (Section 4.2). This is analogous to the classical mean-shift algorithm (Fukunaga and Hostetler, 1975; Cheng, 1995; Comaniciu and Meer, 2002), albeit an important difference is that diffusion-based generative models use time-varying bandwidths. Though not tractable for the practical simulation of sampling trajectories, the KDE-based solution converges to the optimal solution based on the real data distribution asymptotically. Its closed form lends itself to theoretical analysis. We further characterize the deterministic sampling dynamics from the local behavior, i.e., stepwise rotation and monotone likelihood increase, and the global behavior, i.e., the linear-nonlinear-linear mode-seek path of approximately constant length, following naturally from this interpretation of the PF-ODE (Section 4.3). This trajectory regularity unifies prior observations and clarifies many existing heuristics to expedite diffusion sampling. Using the geometric shape regularity of sampling trajectories, we introduce an efficient and effective accelerated sampling approach based on dynamic programming to determine the optimal time schedule. Experimental results demonstrate that it significantly improves the performance of diffusion-based generative models in a few (≤ 10) function evaluations. Our main contributions are summarized as follows:1

We demonstrate and characterize a strong geometric regularity in deterministic sampling dynamics of diffusion-based generative models, i.e., every sampling trajectory follows a similar "boomerang" shape within an extremely low-dimensional subspace.
We attempt to explain this regularity through the closed-form solution of the denoising trajectory under KDE-based data modeling with time-varying bandwidths. Several derived properties offer insights into both the local and global trajectory structures.

^1. Part of this work was previously published at the forty-first International Conference on Machine Learning (ICML) (Chen et al., 2024).

• We develop a dynamic programming-based method that leverages the trajectory regularity to determine an optimal sampling time schedule. It introduces minimal computational overhead and improves image quality, particularly for few-step inference.

2 Preliminaries

2.1 Generative Modeling with Stochastic Differential Equations

For successful generative modeling, it is essential to connect the data distribution p^d with a manageable, non-informative noise distribution pn. Diffusion models achieve this objective by incrementally introducing white Gaussian noise into the data, effectively obliterating its structures, and subsequently reconstructing the synthesized data from noise samples via a series of denoising steps. A typical choice for pⁿ is an isotropic multivariate normal distribution with zero mean. The forward step can be modeled as a diffusion process {zt} for t ∈ [0, T] starting from the initial condition z⁰ ∼ pd, which corresponds to the solution of an Itˆo stochastic differential equation (SDE) (Song et al., 2021c; Oksendal, 2013)

$$d\mathbf{z}_t = \mathbf{f}(\mathbf{z}_t, t)dt + g(t)d\mathbf{w}_t, \quad \mathbf{f}(\cdot, t) : \mathbb{R}^d \to \mathbb{R}^d, \quad g(\cdot) : \mathbb{R} \to \mathbb{R},$$ (1)

where w^t denotes the Wiener process; f(·, t) is a vector-valued function referred to as drift coefficient and g(·) is a scalar function referred to as diffusion coefficient. The noise term in this case is independent of the state z^t (a.k.a. additive noise), and therefore the Itˆo and Stratonovich interpretations of the above SDE coincide (Stratonovich, 1968; S¨arkk¨a and Solin, 2019). A unique, strong solution of this SDE exists when the time-varying drift and diffusion coefficients are globally Lipschitz in both state and time (Oksendal, 2013). The temporal marginal distribution of z^t is denoted as pt(zt), with p0(z0) = pd(z0). By properly setting the coefficients and terminal time T, the data distribution p^d is smoothly transformed to the approximate noise distribution p^T (z^T ) ≈ pⁿ in a forward manner.

The solutions to Itˆo SDEs are always Markov processes, and they can be fully characterized by the transition kernel pst(z^t |zs) with 0 ≤ s < t ≤ T. This transition kernel becomes a Gaussian distribution when considering the linear SDE with an affine drift coefficient f(z^t , t) = f(t)z^t . In this case, data-noise pairs {(z0, zt)} with different levels of noise can be directly generated by ancestral sampling from the joint distribution pd(z0)p0t(z^t |z0) rather than pd(z0) Q^t ^k=1 p(k−1)^k (zk|zk−1), which largely simplifies the computation of forward process and ease the model training (Section 2.2), and therefore linear SDEs are widely used in practice.2 The transition kernel p0t(z^t |z0) in this context, which is derived with standard techniques (S¨arkk¨a and Solin, 2019; Karras et al., 2022), has the following form

$$p_{0t}(\mathbf{z}_t|\mathbf{z}_0) = \mathcal{N}\left(\mathbf{z}_t; s(t)\mathbf{z}_0, s^2(t)\sigma^2(t)\mathbf{I}\right), \tag{2}$$

or, z^t = s(t)z⁰ + [s(t)σ(t)] ϵ^t , where s(t) = exp(R ^t 0 ^f(ξ)dξ), ^σ(t) = qR ^t 0 [g(ξ)/s(ξ)]² dξ, and ϵ^t ∼ N (0, I). For notation simplicity, we hereafter denote them as s^t and σ^t , respectively. Then, we can rewrite the forward linear SDE (1) in terms of s^t and σ^t ,

$$d\mathbf{z}_t = \frac{d\log s_t}{dt}\mathbf{z}_t dt + s_t \sqrt{\frac{d\sigma_t^2}{dt}} d\mathbf{w}_t, \quad f(t) = \frac{d\log s_t}{dt}, \quad \text{and} \quad g(t) = s_t \sqrt{\frac{d\sigma_t^2}{dt}}.$$ (3)

^2. Some non-linear diffusion-based generative models also exist (Zhang and Chen, 2021; Chen et al., 2022; Liu et al., 2023a), but they are beyond the scope of this paper.

Table 1: We demonstrate that the original notations of VP-SDE, VE-SDE (Song et al., 2021c), Flow Matching (Lipman et al., 2023; Liu et al., 2023b; Albergo et al., 2023) can be recovered by properly setting the coefficients s^t and σ^t . As for the VP-SDE, β^t = −d log αt/dt, and α^t ∈ (0, 1] as a decreasing sequence with α⁰ = 1, α^T ≈ 0. As for the VE-SDE, σ^t is an increasing sequence with σ⁰ ≈ 0, σ^T ≫ 1. A typical flow matching-based instantiation defines the transition kernel directly without relying on the forward SDE.

	st	σt	Transition Kernel p0t(zt z0)	Forward SDE
VP-SDE	√ αt	p (1 − αt)/αt	√ N (zt ; αtz0,(1 − αt)I)	√ 1 dzt = − βtzt dt + βtdwt 2
VE-SDE	1	σt	; z0, σ2 N (zt t I)	dzt p 2 = dσ t /dt dwt
Flow Matching	1 − t/T	t/(T − t)	; (1 2I) N (zt − t/T)z0,(t/T)	—

Furthermore, following previous works (Kingma et al., 2021; Rombach et al., 2022), we define the signal-to-noise ratio (SNR) of the transition kernel (2) as SNR(t) = s 2 ^t /(s 2 t σ 2 t ) = 1/σ² t , which is a monotonically non-increasing function. A simple corollary is that any linear diffusion processes with the same σ^t exhibit identical signal-to-noise ratios at all time steps. In the literature, two specific forms of linear SDEs are widely used in large-scale diffusion models (Ramesh et al., 2022; Saharia et al., 2022; Rombach et al., 2022; Balaji et al., 2022; Peebles and Xie, 2023; Podell et al., 2024; Xie et al., 2025), namely, the variancepreserving (VP) SDE and the variance-exploding (VE) SDE (Song et al., 2021c; Karras et al., 2022). They correspond to the continuous versions of previously established models, i.e., DDPMs (Ho et al., 2020; Nichol and Dhariwal, 2021) and NCSNs (Song and Ermon, 2019, 2020), respectively. We summarize their notations, including recently proposed flow matching-based generative models (Lipman et al., 2023; Liu et al., 2023b; Albergo et al., 2023; Neklyudov et al., 2023; Heitz et al., 2023; Esser et al., 2024), in Table 1.

The reversal of the forward linear SDE as expressed in (3) is represented by another backward SDE, which facilitates the synthesis of data from noise through a backward sampling (Feller, 1949; Anderson, 1982). Based on the well-known Fokker-Planck-Kolmogorov (FPK) equation that describes the evolution of pt(zt) given the initial condition p0(z0) = pd(z0) (Oksendal, 2013), i.e.,

$$\frac{\partial p_t(\mathbf{z}_t)}{\partial t} = -\nabla \cdot \left[ p_t(\mathbf{z}_t) f(t) \mathbf{z}t - \frac{g^2(t)}{2} \nabla{\mathbf{z}_t} p_t(\mathbf{z}_t) \right], \tag{4}$$

it is straightforward to verify that a family of backward diffusion processes with varying η^t , as described by the following formula, all maintain the same temporal marginal distributions {pt(zt)} T ^t=0 as the forward SDE at each time throughout the diffusion process

$$d\mathbf{z}{t} = \left[ f(t)\mathbf{z}{t} - \frac{1 + \eta{t}^{2}}{2}g^{2}(t)\nabla{\mathbf{z}{t}}\log p{t}(\mathbf{z}{t}) \right] dt + \eta{t}g(t)d\mathbf{w},$$ (5)

where η^t controls the amount of stochasticity. Notably, there exists a particular deterministic process with the parameter η^t ≡ 0, and it is termed as probability flow ordinary differential equation (PF-ODE) in the literature (Song et al., 2021c; Karras et al., 2022). An equivalent differential equation-based generative modeling starts from a time-dependent vector field and then induces the associated probability path (Lipman et al., 2023; Liu et al., 2023b; Albergo et al., 2023). The deterministic nature of ODE offers several benefits in generative modeling, including efficient sampling, unique encoding, and meaningful latent manipulations (Song et al., 2021c,a; Chen et al., 2024). We thus choose this mathematical formula to analyze the sampling behavior of diffusion models throughout this paper.

2.2 Estimation of Score Function

Simulating the preceding PF-ODE requires having access to the score function ∇z^t log pt(zt) (Hyv¨arinen, 2005; Lyu, 2009), which can be estimated with nonparametric score matching, implicit score matching (Hyv¨arinen, 2005), sliced score matching (Song et al., 2019), or more typically, denoising score matching (DSM) (Vincent, 2011; Song and Ermon, 2019; Karras et al., 2022). The DSM objective function of training a score-estimation model sθ(z^t ;t) is

$$\mathcal{L}{\text{DSM}}\left(\boldsymbol{\theta}; \lambda(t)\right) := \int{0}^{T} \lambda(t) \mathbb{E}{\mathbf{z}{0} \sim p{d}} \mathbb{E}{\mathbf{z}{t} \sim p{0t}(\mathbf{z}{t}|\mathbf{z}{0})} | \overbrace{\nabla{\mathbf{z}{t}} \log p{\boldsymbol{\theta}}(\mathbf{z}{t}; t)}^{s{\boldsymbol{\theta}}(\mathbf{z}{t}; t)} - \nabla{\mathbf{z}{t}} \log p{0t}(\mathbf{z}{t}|\mathbf{z}{0}) |{2}^{2} dt. \tag{6}$$

The weighting function λ(t) across different noise levels reflects our preference for visual quality or density estimation during model training (Song et al., 2021b; Kim et al., 2022). The optimal estimator s ⋆ θ (z^t ;t) equals ∇z^t log pt(zt), and therefore we can use the converged score-estimation model as an effective proxy for the ground-truth score function.

Thanks to a profound connection between the score function and the posterior expectation from the perspective of empirical Bayes (Robbins, 1956; Morris, 1983; Efron, 2010, 2011; Raphan and Simoncelli, 2011; Saremi and Hyv¨arinen, 2019), we can also train a denoising autoencoder (DAE) (Vincent et al., 2008; Bengio et al., 2013b; Alain and Bengio, 2014) to estimate the conditional expectation E(z0|zt), and then convert it to the score function. We summarize this connection as the following lemma.

Lemma 1. Given the clean data z⁰ ∼ p^d and a transition kernel that adds Gaussian noise into the data p0t(z^t |z0) = N zt ; stz0, s² t σ 2 t I , we have the following connection between the score function and the posterior expectation

$$\nabla_{\mathbf{z}_t} \log p_t(\mathbf{z}_t) = \left[ s_t \sigma_t \right]^{-2} \left( s_t \mathbb{E}(\mathbf{z}_0 | \mathbf{z}_t) - \mathbf{z}_t \right), \tag{7}$$

or equivalently (due to the linearity of expectation),

$$\nabla_{\mathbf{z}_t} \log p_t(\mathbf{z}_t) = -\left[s_t \sigmat\right]^{-1} \mathbb{E}{p_{t0}(\mathbf{z}_0|\mathbf{z}_t)} \boldsymbol{\epsilon}_t, \qquad \boldsymbol{\epsilon}_t = \left[s_t \sigma_t\right]^{-1} (\mathbf{z}_t - s_t \mathbf{z}_0). \tag{8}$$

Therefore, we can train a data-prediction model rθ(z^t ;t) to approximate the posterior expectation E(z0|zt), or train a noise-prediction model ϵθ(z^t ;t) to approximate E zt−stz⁰ stσt |zt , and then substitute the score in (5) with the learned model for the diffusion sampling process (Section 2.3). The DAE objective function of training a data-prediction model rθ(z^t ;t) across different noise levels with the weighting function λ(t) is

$$\mathcal{L}{\text{DAE}}\left(\boldsymbol{\theta}; \lambda(t)\right) := \int{0}^{T} \lambda(t) \mathbb{E}{\mathbf{z}{0} \sim p{d}} \mathbb{E}{\mathbf{z}{t} \sim p{0t}(\mathbf{z}{t}|\mathbf{z}{0})} |r{\boldsymbol{\theta}}(\mathbf{z}{t}; t) - \mathbf{z}{0}|{2}^{2} dt. \tag{9}$$

The optimal estimator r ⋆ θ (z^t ;t), also known as Bayesian least squares estimator or minimum mean square error (MMSE) estimator, equals E(z0|zt). Therefore, we have ∇z^t log pt(zt) = [stσ^t −2 (str ⋆ θ (z^t ;t) − zt). In practice, it is assumed that this connection approximately holds for a converged model (rθ(z^t ;t) ≈ r ⋆ θ (z^t ;t) 3 ), and we can plug it into (5) with η^t ≡ 0 to derive the empirical PF-ODE for sampling as follows

$$\frac{d\mathbf{z}_t}{dt} = \frac{d\log s_t}{dt} \mathbf{z}_t - \frac{d\log \sigma_t}{dt} \left( st r{\boldsymbol{\theta}}(\mathbf{z}_t; t) - \mathbf{z}_t \right) = \frac{d\log s_t}{dt} \mathbf{z}_t + s_t \frac{d\sigmat}{dt} \epsilon{\boldsymbol{\theta}}(\mathbf{z}_t; t).$$ (10)

Both the data-prediction model rθ(z^t ;t) and the noise-prediction model ϵθ(z^t ;t) above are widely used in existing works (Ho et al., 2020; Song et al., 2021a; Bao et al., 2022; Karras et al., 2022; Lu et al., 2022b; Zhang and Chen, 2023; Zhou et al., 2024a; Chen et al., 2024).

2.3 Diffusion Sampling with Numerical Approximation

Given the empirical PF-ODE (10), we can synthesize novel samples by first drawing pure noises zˆt^N ∼ pⁿ as the initial condition, and then numerically solving this equation backward with N steps to obtain a sequence {zˆtⁿ } N ⁿ=0 with a certain time schedule Γ = {t⁰ = ϵ, · · · , t^N = T}. We adopt hat notations such as zˆtⁿ to denote the samples generated by numerical methods, which differs from the exact solutions denoted as ztⁿ . The final sample zˆt⁰ is considered to approximately follow the data distribution pd. We designate this sequence as a sampling trajectory generated by the diffusion model.

Generally, we have two formulas to calculate the exact solution from the current position ztn+1 to the next position ztⁿ (t⁰ ≤ tⁿ < tn+1 ≤ t^N ) in the ODE-based sampling to obtain the sampling trajectory from t^N to t0. One is the direct integral from tn+1 to tⁿ

$$\mathbf{z}_{tn} = \mathbf{z}{t{n+1}} + \int{t_{n+1}}^{t_n} \frac{\mathrm{d}\mathbf{z}t}{\mathrm{d}t} \mathrm{d}t = \mathbf{z}{t{n+1}} + \int{t_{n+1}}^{t_n} \left( \frac{\mathrm{d}\log s_t}{\mathrm{d}t} \mathbf{z}_t - \frac{1}{2} s_t^2 \frac{\mathrm{d}\sigmat^2}{\mathrm{d}t} \nabla{\mathbf{z}_t} \log p_t(\mathbf{z}_t) \right), \quad (11)$$

and another leverages the semi-linear structure in the PF-ODE to derive the following equation with the variant of constants formula (Lu et al., 2022a; Zhang and Chen, 2023)

$$\mathbf{z}_{tn} = \exp\left(\int{t_{n+1}}^{tn} f(t) dt\right) \mathbf{z}{t{n+1}} - \int{t_{n+1}}^{tn} \left(\exp\left(\int{t}^{tn} f(r) dr\right) \frac{g^2(t)}{2} \nabla{\mathbf{z}_t} \log p_t(\mathbf{z}_t)\right) dt$$

$$= \frac{s_{tn}}{s{t{n+1}}} \mathbf{z}{t{n+1}} - s{tn} \int{t_{n+1}}^{t_n} \left(s_t \sigma_t \sigmat' \nabla{\mathbf{z}_t} \log p_t(\mathbf{z}_t)\right) dt.$$ (12)

The above integrals, whether in (11) or (12), involving the score function parameterized by a neural network, are generally intractable. Therefore, deterministic sampling in diffusion models centers on approximating these integrals with numerical methods in each step. In practice, various sampling strategies inspired by classic numerical methods have been proposed to solve the backward PF-ODE (10), including the Euler method (Song et al., 2021a), Heun's method (Karras et al., 2022), Runge-Kutta method (Song et al., 2021c; Liu et al., 2022; Lu et al., 2022a), and linear multistep method (Liu et al., 2022; Lu et al., 2022b; Zhang and Chen, 2023; Zhao et al., 2023; Zhou et al., 2024a; Chen et al., 2024).

^3. We slightly abuse the notation and still denote the converged model as rθ(·;t) hereafter.

2.4 The Equivalence of Diffusion Models

We further demonstrate that various diffusion models modeled by linear SDEs are equivalent up to a scaling factor, provided they share the same SNR in their corresponding transition kernels (2). In particular, all other model types (e.g., the VP diffusion process) can be transformed into their VE counterparts via Lemma 2. We also establish the connections between their score functions and sampling behaviors.

Lemma 2. The linear diffusion process defined as (3) can be transformed into its VE counterpart with the change of variables x^t = zt/st, keeping the SNR unchanged.

Proofs are provided in Section A.2. Similarly, we provide the PF-ODE and its empirical version in terms of the x variable (or say, in the x-space) as follows

$$d\mathbf{x}{t} = -\sigma{t} \nabla{\mathbf{x}{t}} \log p{t}(\mathbf{x}{t}) d\sigma{t} = \frac{\mathbf{x}{t} - r{\boldsymbol{\theta}}(\mathbf{x}{t}; t)}{\sigma{t}} d\sigma{t} = \boldsymbol{\epsilon}{\boldsymbol{\theta}}(\mathbf{x}{t}; t) d\sigma_{t}, \tag{13}$$

with the score function ∇x^t log pt(xt) = st∇z^t log pt(zt), for t ∈ [0, T].

Corollary 1. With the same numerical method, the results obtained by solving (11) or (12) are not equal in the general cases. But they become exactly the same by first transforming the formulas into the x-space and then perform numerical approximation.

Corollary 2. With the same numerical method, the result obtained by solving (12) in the z-space is exactly the same as the result obtained by solving (11) or (12) in the x-space.

Because of the above analysis, we can safely remove the drift term in the forward SDE (3) by transforming them into the VE counterparts without changing the essential characteristics of the underlying diffusion model. In the following discussions, we merely focus on the mathematical properties and geometric behaviors of a standardized VE-SDE, i.e.,

$$d\mathbf{x}_t = \sqrt{d\sigma_t^2/dt} \, d\mathbf{w}_t, \quad \sigma_t : \mathbb{R} \to \mathbb{R}, \tag{14}$$

with a pre-defined increasing noise schedule σ^t . Lemma 2 guarantees the applicability of our conclusions to any other types of linear diffusion processes, including the typical flow matching-based models (Liu et al., 2023b; Lipman et al., 2023; Albergo et al., 2023). In this case, the sampling trajectory is denoted as {xˆtⁿ } N ⁿ=0 with the time schedule Γ = {t⁰ = ϵ, · · · , t^N = T} and the initial noise is denoted as xˆt^N ∼ pⁿ = N (0, σ² T I).

2.5 Conditional and Latent Diffusion Models

It is straightforward to extend the above framework of unconditional diffusion models into the conditional variants (Song et al., 2021c; Dhariwal and Nichol, 2021; Rombach et al., 2022). Given the class or text-based condition c, the modeled marginal distributions become pt(x^t |c) (or pt(z^t |c) in the z-space) instead of the original pt(xt), and the sampling process relies on the learned conditional score ∇x^t log pt(x^t |c) at each time. In general, discrete texts are first mapped into a continuous text embedding space (Nichol et al., 2022; Rombach et al., 2022; Saharia et al., 2022), which distinguishes text-conditional diffusion models from the diffusion models conditioned on discrete class labels. Another extension from the practical

Figure 2: Illustration of subspace projection techniques. The deterministic sampling trajectory begins with an initial noise xˆt^N and progresses to the synthesized data xˆt⁰ . (a) The trajectory deviation equals the reconstruction error when the d-dimensional point of the sampling trajectory is projected onto the displacement vector vtN→t⁰ := xˆt⁰ − xˆt^N . (b) We adopt vtN→t⁰ and several top principal components (PCs) from its (d − 1)-dimensional orthogonal complement to approximate the original d-dimensional sampling trajectory.

consideration is performing the diffusion process in a low-dimensional latent space rather than the original high-dimensional data space (Vahdat et al., 2021; Rombach et al., 2022). With the help of an autoencoder structure, latent diffusion models significantly reduce computational demand and scale up to high-resolution generation.

In the following empirical analysis (Section 3), we will demonstrate that strong trajectory regularity is widely present in unconditional, class-conditional, and text-conditional diffusion models. This observation motivates us to investigate the underlying mechanism behind (Section 4) and to develop an improved algorithm for sampling acceleration (Section 5).

3 Geometric Regularity in Deterministic Sampling Dynamics

As mentioned in Section 1, the sampling trajectories of diffusion-based generative models under the PF-ODE framework exhibit a certain regularity in their shapes, regardless of the specific content generated. To better demonstrate this concept, we undertake a series of empirical studies in this section, covering unconditional generation on CIFAR-10 (Krizhevsky and Hinton, 2009), class-conditional generation on ImageNet (Russakovsky et al., 2015), and text-conditional generation in the latent space with Stable Diffusion v1.5 (Rombach et al., 2022). The spatial resolutions used for these diffusion processes are 32 × 32, 64 × 64, 64×64, respectively. Given the complexity of visualizing the entire sampling trajectory and analyzing its geometric characteristics in the original high-dimensional space, we develop subspace projection techniques to better capture the intrinsic structure of diffusion models.

3.1 One-Dimensional Projection

We first examine the trajectory deviation from the straight line connecting the two endpoints, which serves to assess the linearity of the sampling trajectory. A sketch of this

Figure 3: Results of the 1-D trajectory projection. The sampling trajectory exhibits an extremely small trajectory deviation (red curve) compared to the final sample distance (blue curve) in the sampling process. Each trajectory is simulated with the Euler method and 100 number of function evaluations (NFE). The reported average and standard deviations are based on 5,000 randomly generated sampling trajectories, considering variations in initial noises, class labels, and text prompts.

computation is provided in Figure 2a. This approach allows us to align and collectively observe the general behaviors of all trajectories. Specifically, we denote the displacement vector between the two endpoints as vtN→t⁰ := xˆt⁰ − xˆt^N , and compute the trajectory deviation as the perpendicular Euclidean distance (L 2 ) from each intermediate sample xˆtⁿ to the vector vtN→t⁰ , i.e., dtd := q ∥vtn→t⁰ ² − v T tn→t⁰ · vtN→t⁰ /∥vtN→t⁰ ∥2 2 , Additionally, we calculate the L ² distance between each intermediate sample xˆtⁿ and the final sample xˆt⁰ , denoted as dfsd := ∥xˆtⁿ − xˆt⁰ ∥2, and refer to it as the final sample distance.

The empirical results of trajectory deviation dtd and final sample distance dfsd are depicted as the red curves and blue curves in Figure 3, respectively. Note that we use sampling time as the horizontal axis, which allows all sampling trajectories to be aligned and compared both within and across different time slices. From Figures 3a - 3b, we observe that the sampling trajectory's deviation gradually increases from t = 80 to approximately t = 10, then swiftly diminishes as it approaches the final samples. This pattern suggests that initially, each sample might be influenced by various modes, experiencing significant impact, but later becomes strongly guided by its specific mode after a certain turning point. This behavior supports the heuristic approach of arranging time intervals more densely near the minimum timestamp and sparsely towards the maximum one (Song et al., 2021a; Karras et al., 2022; Song et al., 2023; Chen et al., 2024). However, when we consider the ratio of the maximum deviation to the endpoint distance in Figure 3a - 3b, we find that the trajectory deviation is remarkably slight (e.g., 30/8800 ≈ 0.0034 for ImageNet), indicating a pronounced straightness. Additionally, the generated samples along the sampling trajectory tend to move monotonically from their initial points toward their final points (as illustrated by the blue curves). Similar results can be found for the text-conditional generation in the latent space, as shown in Figure 3c.

Figure 4: Visual comparison of trajectory reconstruction for (a) unconditional, (d) classconditional, and (g) text-conditional generation. The real sampling trajectories (top row) are reconstructed using vtN→t⁰ (1-D recon.), along with their top 1 or 2 principal components (2-D or 3-D recon.). To amplify visual differences, we present the denoising outputs of these trajectories. (b/e/h) The L ² distance between the real trajectory samples and their reconstructed counterparts is computed up to 5-D reconstruction. (c/f/i) The variance explained by the top k principal components is reported as the ratio of the sum of the top k eigenvalues to the sum of all eigenvalues.

The trajectory deviation also reflects the reconstruction error if we project all d-dimensional points of the sampling trajectory onto the displacement vector vtN→t⁰ . As demonstrated in Figure 4, the one-dimensional (1-D) approximation proves inadequate, leading to a significant deviation from the actual trajectory both in terms of visual comparison and quantitative results. These observations imply that while all trajectories share a similar macro-structure, the 1-D projection cannot accurately capture the full trajectory structure, probably due to the failure of modeling rotational properties. Therefore, we further develop a multi-dimensional subspace projection technique, as detailed below.

3.2 Multiple-Dimensional Projections

We then implement principal component analysis (PCA) on the orthogonal complement of the displacement vector vtN→t⁰ , which assists in assessing rotational properties of the sampling trajectory. A sketch of this computation is provided in Figure 2b. This (d − 1)-D orthogonal space relative to vtN→t⁰ is denoted as V = {u : u ^T vtN→t⁰ = 0, ∀ u ∈ R ^d}. We begin by projecting each d-D sampling trajectory into V, followed by conducting PCA.

As illustrated in Figure 4, the 2-D approximation using vtN→t⁰ and the first principal component markedly narrows the visual discrepancy with the real trajectory, thereby reducing the L 2 reconstruction error. This finding suggests that all points in each d-D sampling trajectory diverge slightly from a 2-D plane. Consequently, the tangent and normal vectors of the sampling trajectory can be effectively characterized in this manner. By incorporating an additional principal component, we enhance our ability to capture the torsion of the sampling trajectory, thereby increasing the total explained variance to approximately 85% (Figures 4c, 4f, 4i). This improvement allows for a more accurate approximation of the actual trajectory and further reduces the L 2 reconstruction error (Figures 4b, 4e, 4h). In practical terms, this level of approximation effectively captures all the visually pertinent information, with the deviation from the real trajectory being nearly indistinguishable (Figures 4a, 4d, 4g). Consequently, we can confidently utilize a 3-D subspace, formed by two principal components and the displacement vector vtN→t⁰ , to understand the geometric structure of high-dimensional sampling trajectories.

Expanding on this understanding, in Figure 5, we present a visualization of randomly selected sampling trajectories created by diffusion models under various generation settings. Note that the scale along the axis corresponding to xˆt⁰ − xˆt^N is orders of magnitude larger than those of the other two principal components. Since we focus on the geometric shape regularity of sampling trajectories rather than their absolute locations, we align all trajectories via orthogonal transformations to eliminate arbitrary orientation variations. These transformations, including rotations and reflections, are determined by solving the classic Orthogonal Procrustes Problem (Hurley and Cattell, 1962; Sch¨onemann, 1966; Golub and Van Loan, 2013), after which we visualize the calibrated trajectories. Specifically, we represent each projected sampling trajectory in the 3-D subspace as a matrix R ∈ R N×3 , where each row corresponds to a sample coordinate r ∈ R ³ at a particular time step, and N is the total number of generated samples in the trajectory. Then, we seek an orthogonal transformation O ∈ R 3×3 to minimize the Frobenius norm of the residual error matrix E when aligning one matrix R ∈ R ^N×³ with a reference projected trajectory matrix ^R^e ^∈ ^R ^N×³ of

Figure 5: We project 30 sampling trajectories generated by (a) unconditional, (b) classconditional, and (c) text-conditional diffusion models into the 3-D subspaces. Each trajectory is simulated with the Euler method and 100 number of function evaluations (NFE). These trajectories are first aligned to the direction of the displacement vector xˆt⁰ −xˆt^N (this direction is slightly different for each sample), and then projected to the top 2 principal components in the orthogonal space to xˆt⁰ − xˆt^N . See texts for more details.

the same dimensions. This optimization problem is formulated as

$$\mathbf{E} = \min{\mathbf{O}} |\widetilde{\mathbf{R}} - \mathbf{RO}|{\mathrm{F}}^{2}, \quad s.t. \quad \mathbf{O}^{\mathrm{T}}\mathbf{O} = \mathbf{I}_{3\times 3}.$$ (15)

The optimal solution can be easily derived by applying the method of Lagrange multipliers, yielding O^⋆ = UVT, where U and V are obtained via singular value decomposition (SVD) of ^RTR^e , i.e., ^RTR^e ⁼ UΣVT. As in Figure 3, we adopt sampling time as the horizontal axis to better observe the trajectory shapes by aligning different trajectories across time slices. Here, time is scaled by a factor of ^√ d to preserve relative magnitudes between axes. In other words, the visualization remains almost unchanged whether we use the displacement vector vtN→t⁰ or scaled time as the first axis. Besides, an orthogonal transformation ^O^e ^∈ ^R 2×2 that fixes the first axis suffices for alignment in this case.

As a result, the calibrated trajectories depicted in Figure 5 largely adhere to the straight line connecting their endpoints, corroborating the small trajectory deviation observed in our previous findings (Figure 3). Furthermore, Figure 5 accurately depicts the sampling trajectory's behavior, showing its gradual departure from the osculating plane during sampling. Interestingly, each trajectory consistently exhibits a simple, approximately linear-nonlinear structure. This reveals a strong regularity in all sampling trajectories, independent of the specific content generated and variations in initial noises, class labels, and text prompts.

3.3 Three-Dimensional Projection Revisited

Given the strong trajectory regularity of deterministic diffusion sampling manifested in the three-dimensional Euclidean space, as shown in Figure 5, we further resort to a differential geometry tool known as the Frenet–Serret formulas (Do Carmo, 2016) to precisely characterize geometric properties of the projected sampling trajectory.

We denote the projected sampling trajectory consisting of N discrete points in the 3-D subspace as r(ξ), where ξ = T − t with T as the terminal time and t as the sampling time (see Section 2 for detailed notations). Thanks to our proposed subspace projection techniques (Sections 3.1- 3.2), each projected sampling trajectory keeps starting from pure noise r(0) and ends at the synthesized data r(T). The arc-length of this spatial curve is denoted as s(ξ) = R ^ξ 0 ∥r ′ (u)∥du, which is a strictly monotone increasing function. We then use the arc-length s to parameterize the spatial curve, and define the tangent unit vector, normal unit vector, and binormal unit vector as T(s) := r ′ (s), N(s) := r ′′(s)/∥r ′′(s)∥, and B(s) := T(s) × N(s), respectively. These three unit vectors are interrelated, and their relationship is characterized by the well-known Frenet–Serret formulas listed below,

$$\frac{d\mathbf{T}(s)}{ds} = \kappa(s)\mathbf{N}(s), \quad \frac{d\mathbf{N}(s)}{ds} = -\kappa(s)\mathbf{T}(s) + \tau(s)\mathbf{B}(s), \quad \frac{d\mathbf{B}(s)}{ds} = -\tau(s)\mathbf{N}(s), \quad (16)$$

where the curvature κ(s) and the torsion τ (s) measure the deviation of a spatial curve r(s) to be a straight line and planar, respectively. In practice, we generally employ the following expressions for numerical calculation.4

$$\kappa(s) = \left[\mathbf{r}'(s) \times \mathbf{r}''(s)\right] / |\mathbf{r}'(s)|^3, \quad \tau(s) = \left[\left(\mathbf{r}'(s) \times \mathbf{r}''(s)\right) \cdot \mathbf{r}'''(s)\right] / |\mathbf{r}'(s) \times \mathbf{r}''(s)|^2. \quad (17)$$

Specifically, for each discrete point on the spatial curve, we employ its surrounding points within a given window size to estimate the first-, second-, and third-order derivatives based on the third-order Taylor expansion (Lewiner et al., 2005). With the help of least squares fitting using an appropriate window size, the torsion and curvature functions of each projected sampling trajectory can be reliably estimated with small reconstruction errors.

^4. These two equations also hold for the spatial curve r(ξ) parametrized by ξ.

(a) Unconditional generation (CIFAR-10), window size = 101. (b) Class-conditional generation (ImageNet), window size = 101. (c) Text-conditional generation (SDv1.5), window size = 151.

Figure 6: We compute the curvature and torsion functions of each 3-D spatial curve that approximates the original high-dimensional sampling trajectory. The Euler method with 1,000 NFEs is employed to faithfully simulate 30 sampling trajectories with (a) unconditional, (b) class-conditional, and (c) text-conditional diffusion models. The geometric properties of the curves are then estimated using least-squares fitting.

As shown in Figure 6, the projected sampling trajectory remains nearly straight for a large portion (≈ 80%) of the entire sampling process, with both curvature and torsion staying close to zero, aligning with the analysis in Sections 3.1- 3.2. Taking the class-conditional generation as an example (Figure 6b), the curvature and torsion gradually begin to increase when the arc-length s exceeds around 7,000. The peak torsion occurs at around 8,250, which corresponds to the sampling time of around four. Thereafter, the torsion decreases back to near zero. Note that Figure 6 serves as an important complement to Figure 5, helping to build an accurate intuition about the geometric behavior of deterministic sampling dynamics. The significant difference in axis magnitudes in Figure 5 may visually distort the true trajectory shape, whereas Figure 6 provides a more precise description of the rotational structure of the 3-D spatial curve throughout the sampling process. Furthermore, according to the fundamental theorem of curves in differential geometry (Do Carmo, 2016), the shape of any regular curve with non-zero curvature in 3-D space is completely determined by its curvature and torsion. The highly similar curvature and torsion evolution observed in Figure 6 provides strong evidence of trajectory regularity, suggesting that sampling trajectories are congruent, regardless of the specific diffusion model used or generation conditions.

4 Understanding the Deterministic Sampling Dynamics

In this section, we discuss several properties of deterministic sampling in diffusion models and attempt to explain the trajectory regularity observed in the previous section. We first point out that, in addition to the explicit sampling trajectory, there exists another important yet often overlooked trajectory that determines the position of each intermediate sampling point through a convex combination (Section 4.1). We then establish a connection between the sampling dynamics and the Gaussian kernel density estimate (KDE) approximation of data distribution (Section 4.2), which is the key to understanding the observed regularity. We further provide a detailed theoretical analysis of the sampling dynamics (Section 4.3), revealing stepwise rotation and monotone likelihood increase (local behavior), and characterizing the entire sampling trajectory as a linear-nonlinear-linear mode-seeking path of approximately constant length (global behavior).

4.1 Implicit Denoising Trajectory

Given a parametric diffusion model with the denoising output rθ(·; ·), the sampling trajectory is simulated by numerically solving the empirical PF-ODE (13). Meanwhile, an implicitly coupled sequence {rθ(xˆtⁿ , tn)} N ⁿ=1 is formed as a by-product. We designate this sequence, or simplified to {rθ(xˆtⁿ )} N ⁿ=1 if there is no ambiguity, as the denoising trajectory. We then rearrange the empirical PF-ODE (13) as rθ(x^t ;t) = x^t −σt(dxt/dσt), and take the derivative of both sides to obtain the differential equation of the denoising trajectory

$$dr_{\theta}(\mathbf{x}_t;t)/d\sigma_t = -\sigma_t \left[ d^2 \mathbf{x}_t / d\sigma_t^2 \right]. \tag{18}$$

This equation, although not directly applicable for simulation, reveals that the denoising trajectory encapsulates the curvature information of the associated sampling trajectory. The following proposition reveals how these two trajectories are inherently related.

Proposition 1. Given the probability flow ODE (13) and a current position xˆtn+1 , n ∈ [0, N − 1] in the sampling trajectory, the next position xˆtⁿ predicted by a k-order Taylor expansion with the time step size σtn+1 − σtⁿ equals

$$\hat{\mathbf{x}}_{tn} = \frac{\sigma{tn}}{\sigma{t{n+1}}} \hat{\mathbf{x}}{t{n+1}} + \frac{\sigma{t{n+1}} - \sigma{tn}}{\sigma{t{n+1}}} \mathcal{R}{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}), \tag{19}$$

which is a convex combination of xˆtn+1 and the generalized denoising output Rθ(xˆtn+1 ),

$$\mathcal{R}{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}) = r{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}) - \sum{i=2}^{k} \frac{1}{i!} \frac{\mathrm{d}^{(i)} \mathbf{x}{t}}{\mathrm{d}\sigma{t}^{(i)}} \Big|{\hat{\mathbf{x}}{t{n+1}}} \sigma{t{n+1}} (\sigma{t{n}} - \sigma{t{n+1}})^{i-1}.$$ (20)

We have Rθ(xˆtn+1 ) = rθ(xˆtn+1 ) for the Euler method (k = 1), and Rθ(xˆtn+1 ) = rθ(xˆtn+1 )+ σtn−σtn+1 drθ(xˆtn+1 ) dσt for second-order numerical methods (k = 2).

Corollary 3. The denoising output rθ(xˆtn+1 ) reflects the prediction made by a single Euler step from xˆtn+1 with the time step size σtn+1 .

Corollary 4. Each second-order ODE-based accelerated sampling method corresponds to a specific first-order finite difference of drθ(xˆtn+1 )/dσt.

The ratio σtⁿ /σtn+1 in (19) reflects a preference for maintaining the current position rather than transitioning to the generalized denoising output at tn+1. Since the denoising output starts from a spurious mode and converges toward a true mode over time, a reasonable strategy is to decrease this weight progressively during the sampling process. From this perspective, different time-schedule functions designed for diffusion sampling, such as uniform, quadratic, and polynomial schedules (Song et al., 2021a; Lu et al., 2022a; Karras et al., 2022; Chen et al., 2024), essentially represent various weighting functions. This interpretation also motivates the direct search for proper weights using a parametric neural network to improve visual quality (Kingma et al., 2021; Chen et al., 2023a).

We primarily focus on the Euler method to simplify subsequent discussions, though these insights can be readily extended to higher-order methods. The trajectory behavior in the continuous-time scenario is similarly discernible through examining the sampling process with an infinitesimally small Euler step.

4.2 Closed-form Solutions under Kernel Density Estimate

In this section, we show that by using a Gaussian kernel density estimate (KDE) with the training data samples, both the sampling trajectory and the denoising trajectory admit closed-form solutions. This helps establish a relationship between deterministic sampling dynamics and the classic mean-shift algorithm (Fukunaga and Hostetler, 1975; Cheng, 1995; Comaniciu and Meer, 2002), and reveals the memorization effect of diffusion models.

In particular, we focus on the practical setting where we have a dataset D := {yⁱ ∈ R ^d}i∈I containing |I| i.i.d. data points drawn from the unknown data distribution pd. The marginal density at each time step of the forward diffusion process (14) becomes a Gaussian KDE with bandwidth σ 2 t , i.e., ˆpt(xt) = R p0t(x^t |y)ˆpd(y) = (1/|I|) P ⁱ∈I N xt ; yⁱ , σ² t I , where the empirical data distribution ˆpd(y) is a summation of multiple Dirac delta functions, i.e., ˆpd(y) = (1/|I|) P ⁱ∈I δ(∥y − yi∥). The forward process can be paraphrased as the expansion of magnitude and manifold, which means that samples escape from the original small-magnitude low-rank manifold and settle into a large-magnitude high-rank manifold. The squared magnitude of the noisy sample is expected to be larger than that of the original sample. As the dimension d → ∞, the sample magnitude expands with probability one and the isotropic Gaussian noise distributes as a uniform distribution on the sphere (Vershynin, 2018; Chen et al., 2023a). The backward process has an opposite trend due to marginal preservation. As discussed in Section 2.2, once a diffusion model has converged to the optimum, i.e., ∀t, rθ(x^t ;t) → r ⋆ θ (x^t ;t) = E(x0|xt), it captures the score ∇x^t log pt(xt) at different noise levels. The exact formula for the optimal denoising output is provided below

Proposition 2. The optimal denoising output of training a denoising autoencoder with the empirical data distribution is a convex combination of original data points, where each weight is calculated based on the time-scaled and normalized Euclidean distance between the input x^t and yⁱ belonging to the dataset D, i.e.,

$$= \sum{i} \frac{\exp\left(-|\mathbf{x}{t} - \mathbf{y}{i}|{2}^{2}/2\sigma{t}^{2}\right)}{\sum{j} \exp\left(-|\mathbf{x}{t} - \mathbf{y}{j}|{2}^{2}/2\sigma{t}^{2}\right)} \mathbf{y}_{i}.$$ (21)

The above equation appears to be highly similar to the iterative formula used in mean shift (Fukunaga and Hostetler, 1975; Cheng, 1995; Comaniciu and Meer, 2002; Yamasaki and Tanaka, 2020). Mean shift is a non-parametric algorithm designed to seek the modes of a density function, typically a KDE, via iterative gradient ascent with adaptive step sizes. Given a current position x, mean shift with a Gaussian kernel and bandwidth h iteratively adds the vector m(x) − x, which points in the direction of the maximum increase in the KDE: ph(x) = (1/|I|) P|I| ⁱ=1 N (x; yⁱ , h² I), to itself, i.e., x ← [m(x) − x] + x. The mean vector is

$$\mathbf{m}(\mathbf{x}, h) = \sum_{i} \frac{\exp\left(-|\mathbf{x} - \mathbf{y}_i|2^2 / 2h^2\right)}{\sum{j} \exp\left(-|\mathbf{x} - \mathbf{y}_j|_2^2 / 2h^2\right)} \mathbf{y}_i.$$ (22)

As a mode-seeking algorithm, mean shift has shown particular success in clustering (Cheng, 1995; Carreira-Perpin´an, 2015), image segmentation (Comaniciu and Meer, 1999, 2002), and video tracking (Comaniciu et al., 2000, 2003). By treating the bandwidth σ^t in (21) as the bandwidth h in (22), we build a connection between the optimal denoising output of a diffusion model and annealed mean shift under the KDE-based data modeling. Furthermore, the time-decreasing bandwidth (σ^t → 0 as t → 0) in (21) is strongly reminiscent of annealed mean shift, or multi-bandwidth mean shift (Shen et al., 2005), which was developed as a metaheuristic algorithm to escape local maxima, where classical mean shift is susceptible to stuck, by monotonically decreasing the bandwidth in iterations. Additionally, similar to (19), each Euler step in the optimal case equals a convex combination of the annealed mean vector and the current position, with the PF-ODE dxt/dσ^t = (x^t −r ⋆ θ (x^t ;t))/σ^t = ϵ ⋆ θ (x^t ;t).

The above analysis also implies that all ODE trajectories generated by an optimal diffusion model are uniquely determined and governed by a bandwidth-varying mean shift. In this case, the forward (encoding) process and backward (decoding) process depend only on the data distribution and the given noise distribution, regardless of the model architecture or optimization algorithm. This property was previously referred to as uniquely identifiable encoding and empirically verified in (Song et al., 2021c), while we point out its theoretical connection to a global kernel density mode-seeking algorithm (Shen et al., 2005), and thus reveal the asymptotic sampling behavior of diffusion models. Although the optimal diffusion models simply memorize the dataset and replay discrete data points in sampling, we argue that in practice, a slight score deviation from the optimum ensures the generative ability of diffusion models while greatly alleviating the mode collapse issue (see Section B.2).

4.3 Theoretical Analysis of the Trajectory Structure

Next, we further characterize the local behavior and global behavior of deterministic sampling dynamics based on the closed-form solutions in the KDE-based data modeling.

4.3.1 Local Properties

We first show that (1) the denoising output governs the rotation of the sampling trajectory, and (2) each sampling trajectory converges monotonically in terms of sample likelihood, with its coupled denoising trajectory consistently achieving higher likelihood.

Figure 7a illustrates two successive Euler steps according to Proposition 1. The direction of the sampling trajectory (depicted as a polygonal chain with consecutive blue vertices) is controlled by the denoising outputs, while the vertex positions depend on the time schedule. In the optimal case, the sampling path exhibits a similar structure with the PF-ODE dx^t = ϵ ⋆ θ (x^t ;t)dσ^t . This equation defines a special vector field characterized by an approximately constant magnitude ∥ϵ ⋆ θ (x^t ;t)∥² across all marginal distributions pt(xt).

Figure 7: (a) An illustration of two consecutive Euler steps, starting from a current sample xˆtn+1 . A single Euler step in the ODE-based sampling is a convex combination of the denoising output and the current position to determine the next position. Blue points form a piecewise linear sampling trajectory, while red points form the denoising trajectory governing the rotation direction. (b) We have three likelihood orders in the ODE-based diffusion sampling: (1) ph(rθ(xˆtⁿ )) ≥ ph(xˆtⁿ ), (2) ph(xˆtn−¹ ) ≥ ph(xˆtⁿ ), and (3) ph(r ⋆ θ (xˆtⁿ )) ≥ ph(xˆtⁿ ). Note that ph(r ⋆ θ (xˆtⁿ )) may not possess the highest likelihood within the sphere.

Proposition 3. The magnitude of ϵ ⋆ θ (x^t ;t) is approximately distributed around ^√ d, where d denotes the data dimension. Accordingly, the total length of the sampling trajectory is approximately σ^T √ d, where σ^T denotes the maximum noise level.

The above results also hold for practical diffusion models. The proofs and supporting empirical evidence are provided in Section A.5. We further deduce that the position of each intermediate point xˆtⁿ , n ∈ [1, N −1] in the sampling trajectory is primarily determined by the chosen time schedule, given that ∥rθ(xˆtn+1 )−xˆtⁿ ∥² = (σtⁿ /σtn+1 )∥rθ(xˆtn+1 )−xˆtn+1 ∥² = σtⁿ ∥ϵθ(xˆtn+1 )∥² ≈ σtⁿ ∥ϵθ(xˆtⁿ )∥² = ∥rθ(xˆtⁿ ) − xˆtⁿ ∥2. In this scenario, the denoising output rθ(xˆtn+1 ) appears to be oscillating toward rθ(xˆtⁿ ) around xˆtⁿ , akin to the motion of a simple gravity pendulum (Young et al., 1996). The pendulum length effectively shortens by the coefficient σtⁿ /σtn+1 in each sampling step, starting from an initial length of roughly σ^T √ d. This specific structure is shared across all trajectories. In practice, the magnitude of each oscillation is extremely small (≈ 0 ◦ ), and the entire sampling trajectory remains nearly confined to a two-dimensional plane. The minor deviation can be effectively represented using a small number of orthogonal bases, as discussed in Section 3.

Next, we characterize the likelihood orders in the deterministic sampling process. To simplify notations, we denote the deviation of denoising output from the optimal counterpart as d1(xˆtⁿ ) = ∥r ⋆ θ (xˆtⁿ ) − rθ(xˆtⁿ )∥2 and the distance between the optimal denoising output and the current position as d2(xˆtⁿ ) = ∥r ⋆ θ (xˆtⁿ ) − xˆtⁿ ∥² .

Proposition 4. In deterministic sampling with the Euler method, the sample likelihood is non-decreasing, i.e., ∀n ∈ [1, N], we have ph(xˆtn−¹ ) ≥ ph(xˆtⁿ ) and ph(rθ(xˆtⁿ )) ≥ ph(xˆtⁿ ) in terms of the Gaussian KDE ph(x) = (1/|I|) P ⁱ∈I N (x; yⁱ , h² I) with any positive bandwidth h, assuming that all samples in the sampling trajectory satisfy d1(xˆtⁿ ) ≤ d2(xˆtⁿ ).

A visual illustration is provided in Figure 7b. The assumption requires that our learned denoising output rθ(xˆtⁿ ) falls within a sphere centered at the optimal denoising output r ⋆ θ (xˆtⁿ ) with a radius of d2(xˆtⁿ ). This radius controls the maximum deviation of the learned

(a) Unconditional generation (CIFAR-10). Total NFE = 35. (b) Class-conditional generation (ImageNet). Total NFE = 79. (c) Text-conditional generation with SDv1.5. Total NFE = 100.

Figure 8: The comparison of visual quality (top is sampling trajectory, bottom is denoising trajectory) and Fr´echet Inception Distance (FID (Heusel et al., 2017), lower is better) w.r.t. the number of score function evaluations (NFEs). The denoising trajectory converges much faster than the sampling trajectory in terms of FID and visual quality.

denoising output and shrinks during the sampling process. In practice, this assumption is relatively easy to satisfy for a well-trained diffusion model. Therefore, each sampling trajectory monotonically converges (ph(xˆtn−¹ ) ≥ ph(xˆtⁿ )), and its coupled denoising trajectory converges even faster (ph(rθ(xˆtⁿ )) ≥ ph(xˆtⁿ )) in terms of the sample likelihood. Given an empirical data distribution, Proposition 4 applies to any marginal distribution of the forward SDE {pt(x)} T ^t=0, which are a spectrum of kernel density estimates with the varying positive bandwidth t. Furthermore, with an infinitesimal step size, Proposition 4 can be extended to a continuous-time version. Finally, the standard monotone convergence property of mean shift is recovered when considering diffusion models trained to optimal parameters.

Corollary 5. We have $$ph(\mathbf{m}(\hat{\mathbf{x}}{t_n})) \geq ph(\hat{\mathbf{x}}{tn})$$ , when $r{\theta}(\hat{\mathbf{x}}_{tn}) = r{\theta}^{\star}(\hat{\mathbf{x}}_{tn}) = \mathbf{m}(\hat{\mathbf{x}}{t_n})$ .

Besides the monotone increase in sample likelihood, a similar trend is observed in image quality (Figure 8). Both qualitative and quantitative results indicate that the denoising trajectory converges significantly faster than the sampling trajectory. This observation motivates the development of a new technique, termed "ODE-Jump", which enables a direct transition from any sample at any time step in the sampling trajectory to its corresponding point in the denoising trajectory, and returns the denoising output as the final synthetic result. Specifically, we modify the standard sampling sequence from xˆt^N → · · · → xˆtⁿ → · · · → xˆt⁰ to xˆt^N → · · · → xˆtⁿ → rθ(xˆtⁿ ), thereby reducing the total NFE from N to N − n + 1, assuming one NFE per sampling step. This technique is highly flexible and simple to implement. It only requires monitoring the visual quality of intermediate denoising samples to determine an appropriate time to terminate the remaining steps. As an example, consider the sampling process with SDv1.5 shown in Figure 8c. We perform a jump from NFE=80 of the sampling trajectory to NFE=81 of the denoising trajectory. In this case, we

Figure 9: Further analysis of deterministic sampling dynamics based on optimal denoising outputs, includes (a/d) the curvature and torsion functions; (b) an illustration of the bi-level convex combination used to infer the next position; (e) the evolution of convex combination coefficients; and (c/f) the corresponding Shannon entropy along the sampling trajectories.

achieve a substantial improvement in FID, while obtaining a visually comparable result to the final sample at NFE=100 with significantly fewer NFEs. Additional results are provided in Figure 10. Figure 8 also highlights the insensitivity of FID to subtle differences in image quality, a limitation also noted in previous work (Kirstain et al., 2023; Podell et al., 2024).

4.3.2 Global Properties

We then show that (1) the sampling trajectory acts as a linear-nonlinear-linear mode-seeking path, and (2) the trajectory statistics change dramatically in a short phase transition period.

In the optimal case, the denoising output, also referred to as the annealed mean vector, starts from a spurious mode (approximately the dataset mean), i.e., r ⋆ θ (x^t ;t) ≈ (1/|I|) P ⁱ∈I yⁱ when the bandwidth σ^t is sufficiently large. Meanwhile, the sampling trajectory is initially located in an approximately uni-modal Gaussian distribution with a linear score function:

The first linear stage: $$\nabla_{\mathbf{x}_t} \log p_t(\mathbf{x}t) = (r{\boldsymbol{\theta}}^{\star}(\mathbf{x}_t; t) - \mathbf{x}_t) / \sigma_t^2 \approx -\mathbf{x}_t / \sigma_t^2$$ . (23)

This approximation is valid for large t, because the norm of the dataset mean is negligible relative to the norm of x^t due to the concentration of measure (see Lemma 4). This supports the validity of some heuristic methods proposing to use the Gaussian analytic score to replace the learned score at the first step (Dockhorn et al., 2022; Zhou et al., 2024a; Wang and Vastola, 2024). As σ^t monotonically decreases during sampling, the number of modes in the Gaussian KDE ˆpt(xt) = (1/|I|) P ⁱ∈I N xt ; yⁱ , σ² t I increases (Silverman, 1981), and the underlying distribution surface gradually shifts from a simple Gaussian to a complex multi-modal form. In this medium σ^t stage, the score function appears highly datadependent and nonlinear, since multiple data points exert a non-negligible effect. Finally,

with a sufficiently small bandwidth σ^t , the sampling trajectory is attracted to a specific real-data mode, and the score function appears approximately linear again, i.e.,

The second linear stage: $$\nabla_{\mathbf{x}_t} \log p_t(\mathbf{x}t) = (r{\theta}^{\star}(\mathbf{x}_t; t) - \mathbf{x}_t) / \sigma_t^2 \approx (\mathbf{y}_k - \mathbf{x}_t) / \sigma_t^2,$$ (24)

where y^k denotes the nearest data point to x^t . In other words, the posterior distribution can be well approximated by a Gaussian distribution again. This global linear-nonlinear-linear behavior allows the sampling trajectory to locate a true mode under mild conditions, similar to annealed mean shift (Shen et al., 2005). Intriguingly, we can provide a guarantee that the total trajectory length is around σ^T √ d (see Proposition 3), which implies a common structural property shared across trajectories originating from different initial conditions.

A straightforward piece of quantitative evidence supporting the above analysis comes from the trajectory statistics computed from the generated sampling paths based on optimal denosing outputs, as discussed in Section 3.3. The statistics presented in Figures 9a and 9d exhibit a distinct three-stage pattern, with the first and second linear stages characterized by near-zero curvature and torsion. The critical transition points can be readily identified by setting a threshold. For example, the curvature can be considered sufficiently small (e.g., less than 1e −5 , which corresponds to sampling times of 14.80 or 3.44 on average) to define a linear trajectory. We next delve into a more detailed analysis of the phase transition between linear and nonlinear stages in the sampling dynamics. Specifically, given the current position xˆtn+1 and the corresponding optimal denoising output r ⋆ θ (xˆtn+1 ), the next position xˆtⁿ predicted by the Euler method according to (19) becomes

$$\hat{\mathbf{x}}{t{n}} = \frac{\sigma{t{n}}}{\sigma{t{n+1}}} \hat{\mathbf{x}}{t{n+1}} + \frac{\sigma{t{n+1}} - \sigma{t{n}}}{\sigma{t{n+1}}} r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n+1}}) = \sum{i} \underbrace{\sum{j} \exp\left(-|\hat{\mathbf{x}}{t{n+1}} - \mathbf{y}{i}|{2}^{2}/2\sigma{t}^{2}\right)}{\mathbf{u}{i}(\hat{\mathbf{x}}{t{n+1}})} \underbrace{\left(\frac{\sigma{t{n}}}{\sigma{t{n+1}}} \hat{\mathbf{x}}{t{n+1}} + \frac{\sigma{t{n+1}} - \sigma{t{n}}}{\sigma{t{n+1}}} \mathbf{y}{i}\right)}{\hat{\mathbf{y}}{i}(\hat{\mathbf{x}}{t_{n+1}})}. \tag{25}$$

This implies that xˆtⁿ lies within a convex hull, where each vertex yˆⁱ is formed by a convex combination of the data point yⁱ and the current position xˆtn+1 , with time-schedule dependent coefficients (σtⁿ /σtn+1 , 1 − σtⁿ /σtn+1 ), as illustrated in Figure 9b. In contrast, another convex combination, parameterized by coefficients ui(xtn+1 ), quantifies the relative influence of each individual data point and plays a critical role in determining the transition points within the linear-nonlinear-linear path. As shown in Figure 9e, where a logarithmic scale is used with a small bias term 1e ⁻¹⁰ added for numerical stability, the evolution of coefficients begins with an approximately uniform distribution around 1/50,000 ≈ 10−4.⁷ and this distribution gradually converges toward a specific data point. Note that the influence of different data points varies: some decrease steadily, while others increase initially before declining. This behavior suggests the existence of multiple hierarchical clusters within the dataset, potentially reflecting corase-to-fine semantic structures in representing learning (Bengio et al., 2013a). Furthermore, there exists a non-eligible period (roughly the final 5% of the trajectory in Figure 9e) during which a specific data point already dominates the trajectory's trend. This phenomenon becomes more pronounced in Figures 9c and 9f, where we introduce a Shannon entropy-based criterion, H(u) = − P ⁱ ui(xtn+1 ) log ui(xtn+1 ), along with its first-order temporal derivative ∂H(u)/∂t to quantify and visualize the evolving influence of data points throughout the sampling dynamics.

5 Geometry-Inspired Time Scheduling

In this section, as a simple illustration, we propose a new technique inspired by the geometric regularity of deterministic sampling in diffusion models to accelerate sampling and enhance sample quality. This technique is compatible with any numerical solver-based sampler and model architecture, easy to implement, and incurs negligible computational overhead.

5.1 Algorithm

A deterministic ODE-based numerical solver such as the Euler (Song et al., 2021c) or Runge-Kutta (Liu et al., 2022; Zhang and Chen, 2023) relies on a pre-defined time schedule Γ = {t⁰ = ϵ, · · · , t^N = T} in the sampling process. Typically, given the initial time t^N and the final time t0, the intermediate time steps from t¹ to tN−¹ are determined by heuristic strategies such as uniform, quadratic (Song et al., 2021a), log-SNR (Lu et al., 2022a,b), and polynomial functions (Karras et al., 2022; Song et al., 2023). In fact, the time schedule reflects our prior knowledge of the sampling trajectory shape. Under the constraint of the total number of score function evaluations (NFEs), an improved time schedule can reduce the local truncation error in each numerical step, and hopefully minimize the global truncation error. In this way, the sample quality generated by numerical methods could approach that of the exact solutions of the given empirical PF-ODE (13).

Our previous discussions in Section 3 identified each sampling trajectory as a simple lowdimensional "boomerang" curve. We thus leverage this geometric structure to re-allocate the intermediate timestamps according to the principle that assigning a larger time step size when the trajectory exhibits a relatively small curvature, while assigning a smaller time step size when the trajectory exhibits a relatively large curvature. Additionally, different trajectories share almost the same shape, regardless of the model architecture used or generation conditions, which helps us estimate the common structure of the sampling trajectory by using just a few "warmup" samples. We name our approach to achieve the above goal as geometry-inspired time scheduling (GITS) and elaborate the details as follows.

The allocation of the intermediate timestamps can be formulated as an integer programming problem and solved using standard dynamic programming (DP) to search for an optimal time schedule (Cormen et al., 2022).5 We first define a searching space denoted as Γg, which is a fine-grained grid including all possible intermediate timestamps. Then, we measure the trajectory curvature by the local truncation errors. More precisely, we define the cost from the current position xtⁱ to the next position xt^j as the difference between an Euler step and the ground-truth prediction, i.e., cti→t^j := D(xˆti→t^j , xti→t^j ), where tⁱ and t^j are two intermediate timestamps from Γ^g and tⁱ > t^j . According to the empirical PF-ODE (13), the ground-truth prediction is calculated as xti→t^j = xtⁱ + R ^t^j ti ϵθ(xt)σ ′ ^tdt, and the Euler prediction is calculated as xˆti→t^j = xtⁱ + (σt^j − σtⁱ )ϵθ(xˆtⁱ )σ ′ ti . The cost function D can be defined as the Euclidean distance in the original pixel space, or any other user-specified metric. Given all computed pairwise costs, which form a cost matrix, the problem reduces to a standard minimum-cost path problem and can be solved with dynamic programming. Since the global truncation error is not equal to the accumulation of local

^5. We also tried the Branch and Bound algorithm (Land and Doig, 1960) and obtained similar results. Nevertheless, alternative approaches exist for determining the time schedule, such as using a trainable neural network (Tong et al., 2025; Frankel et al., 2025), by leveraging our discovered trajectory regularity.

truncation errors at each step, we introduce a hyperparameter γ, analogous to the discount factor used in reinforcement learning (Sutton et al., 1998), to compensate for this effect.

Dynamic programming is a fundamental concept widely used in computer science and many other fields (Cormen et al., 2022). Watson et al. (2021) was the first one leveraging dynamic programming to re-allocate the time schedule in diffusion models. However, our motivation differs significantly from that of this previous work. Watson et al. (2021) exploited the fact that the evidence lower bound (ELBO) can be decomposed into separate KL terms and utilized DP to find the optimal discrete-time schedule that maximizes the training ELBO. However, this strategy was reported to worsen sample quality, as acknowledged by the authors. In contrast, we first discovered a strong trajectory regularity shared by all sampling trajectories, and then used several "warmup" samples to estimate the trajectory curvature to determine a more effective time schedule for the sampling of diffusion models.

5.2 Experimental Results

We adhere to the setup and experimental designs of the EDM framework (Karras et al., 2022; Song et al., 2023), with ^f(t) = 0, ^g(t) = ^√ 2t, and σ^t = t. Under this parameterization, the forward VE-SDE is expressed as dx^t = √ 2t dw^t , while the corresponding empirical PF-ODE is formulated as dxt/dt = (x^t −r^θ (x^t ;t))/t. The temporal domain is segmented using a polynomial function tⁿ = (t 1/ρ ⁰ + n N (t 1/ρ ^N − t 1/ρ 0 ))^ρ , where t⁰ = 0.002, t^N = 80, n ∈ [0, N], and ρ = 7. We initiate the dynamic programming experiments with 256 "warmup" samples randomly selected from Gaussian noise to create a more refined grid, and then calculate the associated cost matrix. The ground-truth predictions are generated by iPNDM (Zhang and Chen, 2023), which employs a fourth-order multistep Runge-Kutta method with a lowerorder warming start, using the polynomial time schedule with 60 NFEs. This yields a grid size of |Γg| = 61. The default classifier-free guidance scale of 7.5 is used for Stable Diffusion (SDv1.5). We follow the standard FID and CLIP Score evaluation protocol for SDv1.5, using the reference statistics and 30k sampled captions from the MS-COCO validation set (Lin et al., 2014). For other datasets, we compute FID based on 50k generated samples (Heusel et al., 2017). All reported results for evaluated methods are obtained based on our developed open-source toolbox: https://github.com/zju-pi/diff-sampler.

Image generation. As shown in Tables 2-3, our simple time re-allocation strategy based on iPNDM (Zhang and Chen, 2023) consistently outperforms all existing ODE-based accelerated sampling methods with a significant margin, especially in the few NFE cases. In particular, all time schedules in these datasets are searched based on the Euler method, i.e., DDIM (Song et al., 2021a), but they are directly applicable to high-order methods such as iPNDM (Zhang and Chen, 2023). The trajectory regularity we uncovered guarantees that the schedule determined through 256 "warmup" samples is effective across all generated content. Furthermore, the experimental results suggest that identifying this trajectory regularity enhances our understanding of the mechanisms of diffusion models. This understanding opens avenues for developing tailored time schedules for more efficient diffusion sampling. Note that we did not adopt the analytical first step (AFS) that replaces the first numerical step with an analytical Gaussian score to save one NFE, proposed in (Dock horn et al., 2022) and later used in (Zhou et al., 2024a), as we found AFS is particularly effective only for datasets with low-resolution images. DPM-Solver-2 (Lu et al., 2022a) and

Table 2: Sample quality comparison in terms of Fr´echet Inception Distance (FID (Heusel et al., 2017), lower is better) on four datasets (resolutions ranging from 32×32 to 256×256). †: Results reported by authors. More results are provided in Table 9.

METHOD 5 CIFAR-10 32×32 (Krizhevsky and Hinton, 2009) DDIM (Song et al., 2021a) 49.66 DDIM + GITS (ours) 28.05 DPM-Solver-2 (Lu et al., 2022a) - DPM-Solver++(3M) (Lu et al., 2022b) 24.97 DEIS-tAB3 (Zhang and Chen, 2023) 14.39 UniPC (Zhao et al., 2023) 23.98 AMED-Solver (Zhou et al., 2024a) - AMED-Plugin (Zhou et al., 2024a) -	6 35.62 21.04 60.00 11.99 9.40 11.14 7.04	8 22.32 13.30 10.30 4.54 5.55 3.99	10 15.69 10.37 5.01 3.00






			4.09
			2.89
		5.56	4.14
	6.67	3.34	2.48
iPNDM (Zhang and Chen, 2023) 13.59	7.05	3.69	2.77
iPNDM + GITS (ours) 8.38	4.88	3.24	2.49
FFHQ 64×64 (Karras et al., 2019)
DDIM (Song et al., 2021a) 43.93	35.22	24.39	18.37
DDIM + GITS (ours) 29.80	23.67	16.60	13.06
DPM-Solver-2 (Lu et al., 2022a) -	83.17	22.84	9.46
DPM-Solver++(3M) (Lu et al., 2022b) 22.51	13.74	6.04	4.12
DEIS-tAB3 (Zhang and Chen, 2023) 17.36	12.25	7.59	5.56
UniPC (Zhao et al., 2023) 21.40	12.85	5.50	3.84
AMED-Solver (Zhou et al., 2024a) -	10.28	6.90	5.49
AMED-Plugin (Zhou et al., 2024a) -	9.54	5.28	3.66
iPNDM (Zhang and Chen, 2023) 17.17	10.03	5.52	3.98
iPNDM + GITS (ours) 11.22	7.00	4.52	3.62
ImageNet 64×64 (Russakovsky et al., 2015)
DDIM (Song et al., 2021a) 43.81	34.03	22.59	16.72
DDIM + GITS (ours) 24.92	19.54	13.79	10.83
DPM-Solver-2 (Lu et al., 2022a) -	44.83	12.42	6.84
DPM-Solver++(3M) (Lu et al., 2022b) 25.49	15.06	7.84	5.67
DEIS-tAB3 (Zhang and Chen, 2023) 14.75	12.57	6.84	5.34
UniPC (Zhao et al., 2023) 24.36	14.30	7.52	5.53
RES(M)† (Zhang et al., 2023) 25.10	14.32	7.44	5.12
AMED-Solver (Zhou et al., 2024a) -	10.63	7.71	6.06
AMED-Plugin (Zhou et al., 2024a) -	12.05	7.03	5.01
iPNDM (Zhang and Chen, 2023) 18.99	12.92	7.20	5.11
iPNDM + GITS (ours) 10.79	8.43	5.82	4.48
LSUN Bedroom 256×256 (Yu et al., 2015) (pixel-space)
DDIM (Song et al., 2021a) 34.34	25.25	15.71	11.42
DDIM + GITS (ours) 22.04	16.54	11.20	9.04
DPM-Solver-2 (Lu et al., 2022a) -	80.59	23.26	9.61
DPM-Solver++(3M) (Lu et al., 2022b) 23.15	12.28	7.44	5.71
UniPC (Zhao et al., 2023) 23.34	11.71	7.53	5.75
AMED-Solver (Zhou et al., 2024a) -	12.75	6.95	5.38
AMED-Plugin (Zhou et al., 2024a) -	11.58	7.48	5.70
iPNDM (Zhang and Chen, 2023) 26.65	20.73	11.78	5.57
iPNDM + GITS (ours) 15.85	10.41	7.31	5.28

AMED-Solver/Plugin (Zhou et al., 2024a) are thus not applicable with NFE = 5 (marked as "-") in Table 2. Ablation studies on AFS are provided in Table 9.

Table 3: Image generation results using Stable Diffusion v1.5 (two NFEs per sampling step).

| METHOD | Step | | | | | | | |-----------------------------------------------------------------------|----------------|----------------|----------------|----------------|--|--|--| | | 5 | 6 | 7 | 8 | | | | | DPM-Solver++(2M) (Lu et al., 2022b)
DPM-Solver++(2M) + GITS (ours) | 16.80
15.53 | 15.43
13.18 | 14.88
12.32 | 14.65
12.17 | | | |

Figure 10: Top: Visual comparison of samples generated by SDv1.5 using 10-step DPM-Solver++(2M) under various time schedules: (a) uniform, (b) AYS-optimized (Sabour et al., 2024), and (c) GITS-optimized. (d) Results from GITS+Jump, which further reduces the number of steps by 30%, are also presented for comparison. Bottom: FID and CLIP Scores (underlined) for GITS along the trajectories are reported.

A concurrent work named AYS was recently proposed to optimize time schedules for sampling by minimizing the mismatch between the true backward SDE and its linear ap-

TIME SCHEDULE			NFE
	5	6	8	10
DDIM-uniform	36.98	28.22	19.60	15.45
DDIM-logsnr	53.53	38.20	24.06	16.43
DDIM-polynomial	49.66	35.62	22.32	15.69
DDIM + GITS (ours)	28.05	21.04	13.30	10.09
iPNDM-uniform	17.34	9.75	7.56	7.35
iPNDM-logsnr	19.87	10.68	4.74	2.94
iPNDM-polynomial	13.59	7.05	3.69	2.77
iPNDM + GITS (ours)	8.38	4.88	3.24	2.49

Table 4: The comparison of FID results on CIFAR-10 across different time schedules.

proximation, utilizing tools from stochastic calculus (Sabour et al., 2024). In contrast, our GITS exploits the strong trajectory regularity inherent in diffusion models and yields time schedules via dynamic programming, requiring only a small number of "warmup" samples. Our method also gets rid of the time-consuming Monte-Carlo computation in AYS (Sabour et al., 2024) and therefore is several orders of magnitude faster. In Figure 10, we compare image samples generated under different time schedules, using the publicly released colab code and its default setting.6 The text prompts used are "a photo of an astronaut riding a horse on mars" (1st row); "a whimsical underwater world inhabited by colorful sea creatures and coral reefs" (2nd row); "a digital illustration of the Babel tower, 4k detailed, trending in artstation, fantasy vivid colors" (3rd row). The evaluated FID results for each schedule are 14.28 (uniform), 12.48 (AYS), and 12.01 (GITS). Besides, building on the significantly faster convergence of the denoising trajectory compared to the sampling trajectory, as discussed in Section 4.3.1, we propose "GITS-Jump" to further reduce sampling cost by 30% (from 10-step to 7-step), almost without degradation in image quality.

Time schedule comparison. From Table 4, we can see that time schedules considerably affect the image generation performance. Compared with existing handcrafted schedules, the schedule we found better fits the underlying trajectory structure in the sampling of diffusion models and achieves smaller truncation errors with improved sample quality.

Running time. Our strategy is highly efficient and incurs a very low computational cost, without requiring access to the real dataset. The procedure starts with a small number of initial "warmup" samples, followed by executing the given ODE-solver with both fine-grained and coarse-grained steps to construct the cost matrix for dynamic programming. Such a computation is performed only once per dataset, and it yields optimal time schedules for different NFE budgets simultaneously, thanks to the optimal substructure property (Cormen et al., 2022). As reported in Table 5, the entire algorithm takes less than or approximately one minute on datasets such as CIFAR-10, FFHQ, and ImageNet 64 ×64, and around 10 to 15 minutes for larger datasets such as LSUN Bedroom and LAION (Stable Diffusion), when evaluated on an NVIDIA A100 GPU.

Ablation studies. We provide ablation studies on the number of "warmup" sample sizes and the grid size used for generating the fine-grained sampling trajectory in Tables 6 and 7, respectively. The default experiments are conducted using iPNDM+GITS with the

^6. https://research.nvidia.com/labs/toronto-ai/AlignYourSteps/

Table 5: Time (in seconds) used at different stages of GITS. "warmup" samples are generated using 60 NFE, and the NFE budget for dynamic programming is set to 10.

DATASET	sample	cost	dynamic	total
	generation	matrix	programming	time (s)
CIFAR-10 32 × 32	27.47	5.29	0.015	32.78
FFHQ 64 × 64	51.90	10.88	0.016	62.79
ImageNet 64 × 64	71.77	13.28	0.016	85.07
LSUN Bedroom	517.63	122.13	0.015	639.78
LAION (SDv1.5)	877.62	24.00	0.016	901.62

Table 6: Ablation study on the grid size of the dynamic programming-based time scheduling.

| GRID SIZE | | NFE BUDGET | | | | | | | | | |--------------|-------|------------|------|------|------|------|------|--|--|--| | | 4 | 5 | 6 | 7 | 8 | 9 | 10 | | | | | 11 | 20.88 | 10.15 | 5.11 | 4.63 | 3.16 | 2.78 | 2.77 | | | | | 21 | 16.22 | 9.87 | 4.83 | 3.76 | 3.39 | 3.20 | 2.81 | | | | | 41 | 15.34 | 9.34 | 4.83 | 5.54 | 3.01 | 2.66 | 2.53 | | | | | 61 (default) | 15.10 | 8.38 | 4.88 | 5.11 | 3.24 | 2.70 | 2.49 | | | | | 81 | 15.74 | 8.57 | 5.09 | 5.38 | 3.10 | 2.93 | 2.38 | | | | | 101 | 15.03 | 8.72 | 5.02 | 5.19 | 3.12 | 2.81 | 2.41 | | | | | iPNDM | 24.82 | 13.59 | 7.05 | 5.08 | 3.69 | 3.17 | 2.77 | | | |

Table 7: Ablation study on the "warmup" sample size. † indicates that a unique time schedule is searched for each of the 50k generated samples. This special case is more timeconsuming while achieving similar results, owing to the strong trajectory regularity.

| NFE | | SAMPLE SIZE | | | | | | | | | | |-----|------|-------------|-----------|-----------|-----------|-----------|-----------|-----------|--|--|--| | | 1† | 16 | 64 | 128 | 256 | 512 | 1024 | 2048 | | | | | 5 | 9.25 | 9.55±0.75 | 9.57±0.97 | 9.21±0.44 | 8.84±0.30 | 8.81±0.04 | 8.89±0.11 | 8.88±0.12 | | | | | 6 | 5.12 | 5.36±0.61 | 5.16±0.28 | 4.99±0.18 | 5.03±0.25 | 5.20±0.27 | 5.01±0.19 | 4.92±0.08 | | | | | 8 | 3.13 | 3.25±0.13 | 3.22±0.08 | 3.28±0.10 | 3.27±0.11 | 3.30±0.11 | 3.29±0.08 | 3.33±0.10 | | | | | 10 | 2.41 | 2.46±0.11 | 2.46±0.05 | 2.45±0.05 | 2.46±0.04 | 2.45±0.04 | 2.44±0.05 | 2.44±0.05 | | | |

coefficient γ = 1.15 on CIFAR-10. We also provide a sensitivity analysis of the coefficient in Table 9. It is shown that the number of "warmup" samples is not a critical hyperparameter, but reducing it generally increases the variance, as shown in Table 7. Due to subtle differences among sampling trajectories (see Figure 5), we recommend utilizing a reasonable number of "warmup" samples to determine the optimal time schedule, such that this time schedule works well for all the generated samples.

6 Related Work and Discussions

The popular variance-exploding (VE) SDEs (Song and Ermon, 2019; Song et al., 2021c) are taken as our main examples for analysis, which are equivalent to their variance-preserving (VP) counterparts according to Itˆo's lemma (see Lemma 2 and Appendix A.2). The equivalence has been established in their corresponding PF-ODE (rather than SDE) forms by using the change-of-variable formula (see Proposition 1 of (Song et al., 2021a) and Proposition 3 of (Zhang and Chen, 2023)). Karras et al. (2022) also presented a set of steps to express different specific models in a common framework. Instead of training a noise-conditional score model with denoising score matching (Vincent, 2011; Song and Ermon, 2019; Song et al., 2021b) or training a noise-prediction model to estimate the added noise in each step (Ho et al., 2020; Song et al., 2021a; Nichol and Dhariwal, 2021; Vahdat et al., 2021; Bao et al., 2022), we follow (Kingma et al., 2021; Karras et al., 2022) and train a denoising model that predicts the reconstructed data from its corrupted version. With the help of simplified empirical PF-ODE (13), we could characterize an implicit denoising trajectory, draw inspiration from classical non-parametric mean shift (Fukunaga and Hostetler, 1975; Cheng, 1995; Comaniciu and Meer, 2002), and derive various trajectory properties.

Denoising trajectory has been observed since the renaissance of diffusion models (see Figure 6 of (Ho et al., 2020)) and later in Figure 3 of (Kwon et al., 2023), but has not been further investigated, perhaps due to the indirect model parameterization. Karras et al. (2022) first stated that the denoising output reflects the tangent of the sampling trajectory, as shown in our Corollary 3, but neither characterizes the denoising trajectory in differential equations nor discusses how it controls the sampling trajectory. In fact, Karras et al. (2022) mentioned this property to argue that the sampling trajectory of (13) is approximately linear due to the slow change in denoising output, and verified this in a 1D toy example. In contrast, we demonstrate the equivalence of linear diffusion models and provide an in-depth analysis of high-dimensional sampling trajectories with real data, highlighting their intrinsically low dimensionality and pronounced shape regularity.

We then describe a potential application to accelerate the sampling process. Different from most existing methods focusing on developing improved ODE-solvers (Song et al., 2021a; Karras et al., 2022; Liu et al., 2022; Lu et al., 2022a; Zhang and Chen, 2023; Zhao et al., 2023; Zhou et al., 2024a) while selecting time schedules through a handcrafted or trial-and-error way, we leverage the trajectory regularity in deterministic sampling dynamics to more effectively allocate discretized time steps. Our method is extremely faster than those distillation-based accelerated sampling methods (Luhman and Luhman, 2021; Sali mans and Ho, 2022; Zheng et al., 2023; Song et al., 2023; Kim et al., 2024; Zhou et al., 2024b) by several orders of magnitude. While Watson et al. (2021) was the first to employ dynamic programming for optimizing time schedules based on the decomposable nature of the ELBO objective function, their method was shown to degrade sample quality. Moreover, although various theoretical studies have explored the problems in convergence analysis and score estimation of diffusion models, none of them focus on the trajectory-level properties (De Bortoli, 2022; Pidstrigach, 2022; Lee et al., 2023; Chen et al., 2023b,c).

7 Conclusion

We reveal that a strong trajectory regularity consistently emerges in the deterministic sampling dynamics of diffusion-based generative models, regardless of the model implementation or generated content. This regularity is explained by characterizing and analyzing the implicit denoising trajectory, particularly its behavior under kernel density estimate-based data modeling. These insights into the underlying trajectory structure lead to an accelerated sampling method that enhances image synthesis quality with negligible computational overhead. We hope that the empirical and theoretical findings presented in this paper contribute to a deeper understanding of diffusion models and inspire further research into more efficient training paradigms and faster sampling algorithms.

Future works. We aim to explore deeper geometric regularities in sampling trajectories, characterize more precise structural patterns, and identify new applications inspired by these insights. Several promising directions are outlined below:

The geometric regularity of sampling trajectories analyzed in this paper may have potential connections to the behavior of random walk paths simulated in the forward process of diffusion models. In the limit of infinite dimensions and trajectory length, random walk-based trajectories exhibit similarly intriguing low-dimensional structures, with the explained variance taking an analytic form. Furthermore, their projections onto PCA subspaces follow Lissajous curves (Moore et al., 2018; Antognini and Sohl-Dickstein, 2018). Extending existing theoretical results from the forward diffusion process to the backward sampling process remains an open problem.
A distinct three-stage pattern emerges in the sampling dynamics when using optimal score functions. Concurrently with our earlier manuscripts (Chen et al., 2023a, 2024), Biroli et al. (2024) introduced concepts from statistical physics, such as symmetry breaking and phase transitions, to characterize sampling dynamics. They provided analytic solutions for critical points in a simplified setting (two well-separated Gaussian mixture classes), and discussed the trade-off between generalization and memorization. It is particularly intriguing to bridge these theoretical insights with realistic diffusion models, especially incorporating conditional signals into the framework.
In practice, sampling in diffusion-based generative models is typically performed using general-purpose numerical solvers, sometimes augmented with learned solver coefficients or sampling schedules in a data-driven way (Tong et al., 2025; Frankel et al., 2025). Our findings reveal that each integral curve of the gradient field defined by a diffusion model lies within an extremely low-dimensional subspace embedded in the high-dimensional data space, with a regular trajectory shape shared across all initial conditions. While we present a preliminary attempt to exploit this structure, further investigation in this direction holds great promise.

Acknowledgments and Disclosure of Funding

This work is supported by the University at Buffalo Office of Vice President for Research and Economic Development. We don't have competing interests to declare.

Appendix A. Proofs

In this section, we provide detailed proofs of the lemmas, propositions, and theorems presented in the main content.

A.1 Proof of Lemma 1

$$\nabla_{\mathbf{z}_t} \log p_t(\mathbf{z}_t) = [s_t \sigma_t]^{-2} (s_t \mathbb{E}(\mathbf{z}_0 | \mathbf{z}_t) - \mathbf{z}_t), \qquad (26)$$

or equivalently (due to the linearity of expectation),

Proof We take derivative of pt(zt) = R pd(z0)p0t(z^t |z0)dz⁰ with respect to z^t

$$\nabla{\mathbf{z}{t}} p{t}(\mathbf{z}{t}) = \int \frac{(s{t}\mathbf{z}{0} - \mathbf{z}{t})}{s{t}^{2} \sigma{t}^{2}} p{d}(\mathbf{z}{0}) p{0t}(\mathbf{z}{t} | \mathbf{z}{0}) d\mathbf{z}_{0}$$

$$s{t}^{2} \sigma{t}^{2} \nabla{\mathbf{z}{t}} p{t}(\mathbf{z}{t}) = \int s{t} \mathbf{z}{0} p{d}(\mathbf{z}{0}) p{0t}(\mathbf{z}{t} | \mathbf{z}{0}) d\mathbf{z}{0} - \mathbf{z}{t} p{t}(\mathbf{z}_{t})$$

$$s{t}^{2} \sigma{t}^{2} \frac{\nabla{\mathbf{z}{t}} p{t}(\mathbf{z}{t})}{p{t}(\mathbf{z}{t})} = s{t} \int \mathbf{z}{0} p{t}(\mathbf{z}{0} | \mathbf{z}{t}) d\mathbf{z}{0} - \mathbf{z}_{t}$$

$$\nabla{\mathbf{z}{t}} \log p{t}(\mathbf{z}{t}) = [s{t} \sigma{t}]^{-2} (s{t} \mathbb{E}(\mathbf{z}{0} | \mathbf{z}{t}) - \mathbf{z}{t}).$$

$$(28)$$

We further have

$$\nabla{\mathbf{z}{t}} \log p{t}(\mathbf{z}{t}) = [s{t}\sigma{t}]^{-2} \left( s{t}\mathbb{E}(\mathbf{z}{0}|\mathbf{z}{t}) - \mathbf{z}{t} \right) = [s{t}\sigma{t}]^{-1} \mathbb{E}\left( \frac{s{t}\mathbf{z}{0} - \mathbf{z}{t}}{s{t}\sigma{t}} | \mathbf{z}{t} \right)$$

$$= -[s{t}\sigma{t}]^{-1} \mathbb{E}{p{t0}(\mathbf{z}{0}|\mathbf{z}{t})} \boldsymbol{\epsilon}_{t}, \tag{29}$$

due to the linearity of expectation, where z^t = stz⁰ + stσtϵ^t , ϵ^t ∼ N (0, I) according to the transition kernel (2).

A.2 Proof of Lemma 2

Lemma 2. The linear diffusion process defined as (3) can be transformed into its VE counterpart with the change of variables x^t = zt/st, keeping the SNR unchanged.

Proof We adopt the change-of-variables formula x^t = ϕ(t, zt) = zt/s^t with ϕ : [0, T]×R ⁿ → R n , and denote the i-th dimension of z^t , x^t and w^t as z^t [i], x^t [i], and w^t [i] respectively; ϕ = [ϕ1, · · · , ϕⁱ , · · · , ϕn] ^T with a twice differentiable scalar function ϕi(t, z) = z/s^t of two real variables t and z. Since each dimension of z^t is independent, we can apply Itˆo's lemma (Oksendal, 2013) to each dimension with ϕi(t, z^t [i]) separately. We have

$$\frac{\partial \phi_i}{\partial t} = -\frac{z}{s_t^2} \frac{\mathrm{d}s_t}{\mathrm{d}t}, \quad \frac{\partial \phi_i}{\partial z} = \frac{1}{s_t}, \quad \frac{\partial^2 \phi_i}{\partial z^2} = 0, \quad \mathrm{d}\mathbf{z}_t[i] = \frac{\mathrm{d}\log s_t}{\mathrm{d}t} \mathbf{z}_t[i] \mathrm{d}t + s_t \sqrt{\frac{\mathrm{d}\sigma_t^2}{\mathrm{d}t}} \mathrm{d}\mathbf{w}_t[i], \quad (30)$$

then

$$d\phi{i}(t, \mathbf{z}{t}[i]) = \left(\frac{\partial \phi{i}}{\partial t} + f(t)\mathbf{z}{t}[i]\frac{\partial \phi{i}}{\partial z} + \frac{g^{2}(t)}{2}\frac{\partial^{2}\phi{i}}{\partial z^{2}}\right)dt + g(t)\frac{\partial \phi{i}}{\partial z}d\mathbf{w}{t}[i]$$

$$= \left(\frac{\partial \phi{i}}{\partial t} + \frac{g^{2}(t)}{2}\frac{\partial^{2}\phi{i}}{\partial z^{2}}\right)dt + \frac{\partial \phi{i}}{\partial z}d\mathbf{z}{t}[i]$$

$$d\mathbf{x}{t}[i] = -\frac{\mathbf{z}{t}[i]}{s{t}}\frac{\mathrm{d}\log s{t}}{\mathrm{d}t}dt + \frac{1}{s{t}}\left(\frac{\mathrm{d}\log s{t}}{\mathrm{d}t}\mathbf{z}{t}[i]dt + s{t}\sqrt{\frac{\mathrm{d}\sigma{t}^{2}}{\mathrm{d}t}}dw{t}\right)$$

$$d\mathbf{x}{t}[i] = \sqrt{\mathrm{d}\sigma{t}^{2}/\mathrm{d}t}\ d\mathbf{w}{t}[i] \quad \Rightarrow \quad d\mathbf{x}{t} = \sqrt{\mathrm{d}\sigma{t}^{2}/\mathrm{d}t}\ d\mathbf{w}{t},$$

$$(31)$$

with the initial condition x⁰ = z⁰ ∼ pd. Since σ^t in the above VE-SDE (x-space) is exactly the same as that used in the original SDE (z-space, (3)), the SNR remains unchanged. Similarly, we have the score function (t ∈ [0, T])

$$\nabla{\mathbf{x}{t}} \log p{t}(\mathbf{x}{t}) = \nabla{\mathbf{x}{t}} \log \int \mathcal{N}\left(\mathbf{x}{t}; \mathbf{x}{0}, \sigma{t}^{2} \mathbf{I}\right) p{d}(\mathbf{x}{0}) d\mathbf{x}{0}$$

$$= \nabla{\mathbf{z}{t}/s{t}} \log \int \mathcal{N}\left(\mathbf{z}{t}/s{t}; \mathbf{x}{0}, \sigma{t}^{2} \mathbf{I}\right) p{d}(\mathbf{x}{0}) d\mathbf{x}{0}$$

$$= s{t} \nabla{\mathbf{z}{t}} \log \int s{t}^{d} \mathcal{N}\left(\mathbf{z}{t}; s{t} \mathbf{z}{0}, s{t}^{2} \sigma{t}^{2} \mathbf{I}\right) p{d}(\mathbf{z}{0}) d\mathbf{z}{0} = s{t} \nabla{\mathbf{z}{t}} \log p{t}(\mathbf{z}_{t}).$$

$$(32)$$

Corollary 2. With the same numerical method, the result obtained by solving (12) in the z-space is exactly the same as the result obtained by solving (11) or (12) in the x-space.

Proof Given the sample ztⁿ obtained by solving the equation (12) in z-space starting from ztn+1 , we demonstrate that ztⁿ /stⁿ is exactly equal to sampling with the equation (11) in x-space starting from xtn+1 = ztn+1 /stn+1 to xtⁿ . We have

$$\mathbf{z}_{tn} = s{tn} \left( \frac{\mathbf{z}{t{n+1}}}{s{t{n+1}}} - \int{t_{n+1}}^{t_n} s_t \sigma_t \sigmat' \nabla{\mathbf{z}_t} \log p_t(\mathbf{z}_t) dt \right)$$

$$= s_{tn} \left( \mathbf{x}{t{n+1}} + \int{t_{n+1}}^{t_n} -\sigmat \nabla{\mathbf{x}_t} \log p_t(\mathbf{x}_t) \sigmat' dt \right) = s{tn} \left( \mathbf{x}{t{n+1}} + \int{t_{n+1}}^{t_n} \frac{d\mathbf{x}t}{dt} dt \right) = s{tn} \mathbf{x}{t_n}.$$ (33)

A.3 Proof of Proposition 1

All following proofs are conducted in the context of a VE-SDE dx^t = √ 2t dw^t , i.e., σ^t = t for notation simplicity, and the sampling trajectory always starts from xˆt^N ∼ N (0, T² I) and ends at xˆt⁰ . The PF-ODEs of the sampling trajectory and denoising trajectory are dxt ^d^t = ϵθ(x^t ;t) = ^xt−rθ(xt;t) , and ^drθ(xt;t) ^d^t = −t d ²x^t dt ² , respectively.

$$\hat{\mathbf{x}}_{t_n} = \frac{tn}{t{n+1}} \hat{\mathbf{x}}{t{n+1}} + \frac{t_{n+1} - tn}{t{n+1}} \mathcal{R}{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t_{n+1}}). \tag{34}$$

which is a convex combination of xˆtn+1 and the generalized denoising output Rθ(xˆtn+1 ),

$$\mathcal{R}{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}) = r{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}) - \sum{i=2}^{k} \frac{1}{i!} \frac{\mathbf{d}^{(i)} \mathbf{x}{t}}{\mathbf{d}t^{(i)}} \Big|{\hat{\mathbf{x}}{t{n+1}}} t{n+1} (t{n} - t{n+1})^{i-1}.$$ (35)

We have Rθ(xˆtn+1 ) = rθ(xˆtn+1 ) for the Euler method (k = 1), and Rθ(xˆtn+1 ) = rθ(xˆtn+1 )+ tn−tn+1 drθ(xˆtn+1 ) dt for second-order methods (k = 2).

Proof The k-order Taylor expansion at xˆtn+1 is

$$\hat{\mathbf{x}}{t{n}} = \sum{i=0}^{k} \frac{1}{i!} \frac{\mathbf{d}^{(i)} \mathbf{x}{t}}{\mathbf{d}t^{(i)}} \Big|{\hat{\mathbf{x}}{t{n+1}}} (t{n} - t_{n+1})^{i}$$

$$= \hat{\mathbf{x}}{t{n+1}} + (t{n} - t{n+1}) \frac{\mathbf{d} \mathbf{x}{t}}{\mathbf{d}t} \Big|{\hat{\mathbf{x}}{t{n+1}}} + \sum{i=2}^{k} \frac{1}{i!} \frac{\mathbf{d}^{(i)} \mathbf{x}{t}}{\mathbf{d}t^{(i)}} \Big|{\hat{\mathbf{x}}{t{n+1}}} (t{n} - t_{n+1})^{i}$$

$$= \hat{\mathbf{x}}{t{n+1}} + \frac{t{n} - t{n+1}}{t{n+1}} \left( \hat{\mathbf{x}}{t{n+1}} - r{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}) \right) + \sum{i=2}^{k} \frac{1}{i!} \frac{\mathbf{d}^{(i)} \mathbf{x}{t}}{\mathbf{d}t^{(i)}} \Big|{\hat{\mathbf{x}}{t{n+1}}} (t{n} - t_{n+1})^{i}$$

$$= \frac{t{n}}{t{n+1}} \hat{\mathbf{x}}{t{n+1}} + \frac{t{n+1} - t{n}}{t{n+1}} \mathcal{R}{\boldsymbol{\theta}}(\hat{\mathbf{x}}{t{n+1}}), \tag{36}$$

where Rθ(xˆtn+1 ) = rθ(xˆtn+1 ) − P^k i=2 i! (i)x^t dt xˆtn+1 tn+1(tⁿ − tn+1) i−1 . As for the first-order approximation (k = 1), we have Rθ(xˆtn+1 ) = rθ(xˆtn+1 ). As for the second-order approximation (k = 2), we have

$$\mathcal{R}{\theta}(\hat{\mathbf{x}}{t{n+1}}) = r{\theta}(\hat{\mathbf{x}}{t{n+1}}) - \frac{t_{n+1}(tn - t{n+1})}{2} \frac{\mathrm{d}^2 \mathbf{x}t}{\mathrm{d}t^2} \Big|{\hat{\mathbf{x}}{t{n+1}}} = r{\theta}(\hat{\mathbf{x}}{t_{n+1}}) + \frac{tn - t{n+1}}{2} \frac{\mathrm{d}r{\theta}(\hat{\mathbf{x}}{t_{n+1}})}{\mathrm{d}t}.$$ (37)

Corollary 3. The denoising output rθ(xˆtn+1 ) reflects the prediction made by a single Euler step from xˆtn+1 with the time step size tn+1.

Proof The prediction of such an Euler step equals to

$$\hat{\mathbf{x}}{t{n+1}} + (0 - t{n+1}) \left( \hat{\mathbf{x}}{t{n+1}} - r{\boldsymbol{\theta}} (\hat{\mathbf{x}}{t{n+1}}) \right) / t{n+1} = r{\boldsymbol{\theta}} (\hat{\mathbf{x}}{t{n+1}}).$$ (38)

A.4 Proof of Proposition 2

Lemma 3. The optimal estimator r ⋆ θ (x^t ;t), also known as the Bayesian least squares estimator, of the minimization of denoising autoencoder (DAE) objective is E (x0|xt).

$$\mathcal{L}{\text{DAE}} = \mathbb{E}{\mathbf{x}_0 \sim p_d(\mathbf{x}0)} \mathbb{E}{\mathbf{x}t \sim p{0t}(\mathbf{x}_t | \mathbf{x}0)} | r{\boldsymbol{\theta}}(\mathbf{x}_t; t) - \mathbf{x}_0 |_2^2 = \int p_t(\mathbf{x}t) p{t0}(\mathbf{x}_0 | \mathbf{x}t) | r{\boldsymbol{\theta}}(\mathbf{x}_t; t) - \mathbf{x}_0 |_2^2.$$ (39)

Proof The solution is easily obtained by setting the derivative of LDAE equal to zero. For each noisy sample x^t , we have

$$\nabla{r{\boldsymbol{\theta}}(\mathbf{x}{t};t)} \mathcal{L}{\text{DAE}} = \mathbf{0}$$

$$\int p{t0}(\mathbf{x}{0}|\mathbf{x}{t}) \left(r{\boldsymbol{\theta}}^{\star}(\mathbf{x}{t};t) - \mathbf{x}{0}\right) d\mathbf{x}_{0} = \mathbf{0}$$

$$\int p{t0}(\mathbf{x}{0}|\mathbf{x}{t}) r{\boldsymbol{\theta}}^{\star}(\mathbf{x}{t};t) d\mathbf{x}{0} = \int p{t0}(\mathbf{x}{0}|\mathbf{x}{t}) \mathbf{x}{0} d\mathbf{x}_{0}$$

$$r{\boldsymbol{\theta}}^{\star}(\mathbf{x}{t};t) = \mathbb{E}\left(\mathbf{x}{0}|\mathbf{x}{t}\right).$$

$$(40)$$

$$r{\boldsymbol{\theta}}^{\star}(\mathbf{x}{t};t) = \min{r{\boldsymbol{\theta}}} \mathbb{E}{\mathbf{y} \sim \hat{p}{d}} \mathbb{E}{\mathbf{x}{t} \sim p{0t}(\mathbf{x}{t}|\mathbf{y})} |r{\boldsymbol{\theta}}(\mathbf{x}{t};t) - \mathbf{y}|{2}^{2} = \sum{i} \frac{\exp\left(-|\mathbf{x}{t} - \mathbf{y}{i}|{2}^{2}/2\sigma{t}^{2}\right)}{\sum{j} \exp\left(-|\mathbf{x}{t} - \mathbf{y}{j}|{2}^{2}/2\sigma{t}^{2}\right)} \mathbf{y}{i}.$$

$$(41)$$

Proof Based on Lemmas 1 and 3, and the Gaussian kernel density estimate (KDE) ˆpt(xt) = R p0t(x^t |y)ˆpd(y) = ¹ |I| P ⁱ N xt ; yⁱ , σ² t I , the optimal denoising output is

$$r{\theta}^{\star}(\mathbf{x}{t};t) = \mathbb{E}\left(\mathbf{y}|\mathbf{x}{t}\right) = \mathbf{x}{t} + \sigma{t}^{2}\nabla{\mathbf{x}{t}}\log\hat{p}{t}(\mathbf{x}_{t})$$

$$= \mathbf{x}{t} + \sigma{t}^{2}\sum{i}\frac{\nabla{\mathbf{x}{t}}\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{i},\sigma{t}^{2}\mathbf{I}\right)}{\sum{j}\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{j},\sigma{t}^{2}\mathbf{I}\right)}$$

$$= \mathbf{x}{t} + \sigma{t}^{2}\sum{i}\frac{\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{i},\sigma{t}^{2}\mathbf{I}\right)}{\sum{j}\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{j},\sigma{t}^{2}\mathbf{I}\right)}\left(\frac{\mathbf{y}{i} - \mathbf{x}{t}}{\sigma_{t}^{2}}\right)$$

$$= \mathbf{x}{t} + \sum{i}\frac{\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{i},\sigma{t}^{2}\mathbf{I}\right)}{\sum{j}\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{j},\sigma{t}^{2}\mathbf{I}\right)}\left(\mathbf{y}{i} - \mathbf{x}_{t}\right)$$

$$= \sum{i}\frac{\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{i},\sigma{t}^{2}\mathbf{I}\right)}{\sum{j}\mathcal{N}\left(\mathbf{x}{t};\mathbf{y}{j},\sigma{t}^{2}\mathbf{I}\right)}\mathbf{y}_{i}$$

$$= \sum{i}\frac{\exp\left(-|\mathbf{x}{t} - \mathbf{y}{i}|{2}^{2}/2\sigma{t}^{2}\right)}{\sum{j}\exp\left(-|\mathbf{x}{t} - \mathbf{y}{j}|{2}^{2}/2\sigma{t}^{2}\right)}\mathbf{y}_{i}.$$

$$(42)$$

A.5 Proof of Proposition 3

Lemma 4 (see Section 3.1 in Vershynin (2018)). Given a high-dimensional isotropic Gaussian noise z ∼ N (0; σ 2 Id), σ > 0, we have E ∥z∥ ² = σ ²d, and with high probability, z stays within a "thin spherical shell": ∥z∥ = σ √ d ± O(1).

Proof We denote zⁱ as the i-th dimension of random variable z, then the expectation and variance is E [zⁱ ] = 0, V [zⁱ ] = σ 2 , respectively. The fourth central moment is E z 4 i = 3σ 4 . Additionally,

$$\mathbb{E}\left[\mathbf{z}{i}^{2}\right] = \mathbb{V}\left[\mathbf{z}{i}\right] + \mathbb{E}\left[\mathbf{z}{i}\right]^{2} = \sigma^{2}, \qquad \mathbb{E}\left[|\mathbf{z}|^{2}\right] = \mathbb{E}\left[\sum{i=1}^{d} \mathbf{z}{i}^{2}\right] = \sum{i=1}^{d} \mathbb{E}\left[\mathbf{z}_{i}^{2}\right] = \sigma^{2}d,$$

$$\mathbb{V}\left[|\mathbf{z}|^{2}\right] = \mathbb{E}\left[|\mathbf{z}|^{4}\right] - \left(\mathbb{E}\left[|\mathbf{z}|^{2}\right]\right)^{2} = 2d\sigma^{4},$$ (43)

Then, we have

$$\mathbb{E}\left[|\mathbf{x} + \mathbf{z}|^2 - |\mathbf{x}|^2\right] = \mathbb{E}\left[|\mathbf{z}|^2 + 2\mathbf{x}^T\mathbf{z}\right] = \mathbb{E}\left[|\mathbf{z}|^2\right] = \sigma^2 d. \tag{44}$$

Furthermore, the standard deviation of ∥z∥ 2 is σ 2 √ 2d, which means

$$|\mathbf{z}|^2 = \sigma^2 d \pm \sigma^2 \sqrt{2d} = \sigma^2 d \pm O\left(\sqrt{d}\right), \quad |\mathbf{z}| = \sigma\sqrt{d} \pm O\left(1\right),$$ (45)

holds with high probability.

Proof We next provide a sketch of proof. Suppose the data distribution lies in a smooth real low-dimensional manifold with the intrinsic dimension as m. According to the Whitney embedding theorem Whitney (1936), it can be smoothly embedded in a real 2m Euclidean space. We then decompose each ϵ ⋆ ^θ ^∈ ^R ^d vector as ϵ ⋆ θ,∥ and ϵ ⋆ θ,⊥ , which are parallel and perpendicular to the 2m Euclidean space, respectively. Therefore, we have ∥ϵ ⋆ θ ∥² = ∥ϵ ⋆ ^θ,^∥ + ϵ ⋆ θ,⊥ ∥² ≥ ∥ϵ ⋆ θ,⊥ ∥² ≈ √ d − 2m.

We provide a upper bound for the ∥ϵ ⋆ θ ∥² below

$$\mathbb{E}{p{t}(\mathbf{x}{t})} | \boldsymbol{\epsilon}{\boldsymbol{\theta}}^{\star} |{2} = \mathbb{E}{p{t}(\mathbf{x}{t})} | \frac{\mathbf{x}{t} - r{\boldsymbol{\theta}}^{\star}(\mathbf{x}{t})}{\sigma{t}} |{2} = \mathbb{E}{p{t}(\mathbf{x}{t})} | \frac{\mathbf{x}{t} - \mathbb{E}(\mathbf{x}{0}|\mathbf{x}{t})}{\sigma{t}} |{2} \ = \mathbb{E}{p{t}(\mathbf{x}{t})} | \mathbb{E}(\frac{\mathbf{x}{t} - \mathbf{x}{0}}{\sigma{t}} | \mathbf{x}{t}) |{2} = \mathbb{E}{p{t}(\mathbf{x}{t})} | \mathbb{E}{p{t0}(\mathbf{x}{0}|\mathbf{x}{t})} \boldsymbol{\epsilon} |{2} \ \leq \mathbb{E}{p{t}(\mathbf{x}{t})} \mathbb{E}{p{t0}(\mathbf{x}{0}|\mathbf{x}{t})} | \boldsymbol{\epsilon} |{2} = \mathbb{E}{p{0}(\mathbf{x}{0})} \mathbb{E}{p{0t}(\mathbf{x}{t}|\mathbf{x}{0})} | \boldsymbol{\epsilon} |_{2} \ \approx \sqrt{d} \qquad \text{(concentration of measure, Lemma 4)}.$$

Figure 11: The optimal noise prediction satisfies ∥ϵ ⋆ θ ∥² ≈ d throughout the entire sampling process, as guaranteed by theoretical results. The actual noise prediction √ ∥ϵθ∥² also remains d for most time steps, but exhibits a noticeable shrinkage in the final stage, when the time step approaches zero. This norm shrinkage almost does not affect the trajectory length, as the discretized time steps in the final stage are extremely small.

Table 8: PCA ratios computed by directly performing dimension reduction of the original sampling trajectory. This differs from PCA ratios in the main context that measures the reconstruction of the orthogonal complement of the displacement vector of each trajectory.

Dataset	Top-1	Top-2	Top-3
CIFAR-10 (32 × 32)	99.99729±0.00135%	99.99998±0.00018%	99.999994±0.000027%
FFHQ (64 × 64)	99.99843±0.00095%	99.99999±0.00016%	99.999994±0.000024%
ImageNet (64 × 64)	99.99805±0.00121%	99.99998±0.00019%	99.999994±0.000033%
SDv1.5, MS-COCO (64 × 64)	99.94101±0.01289%	99.99432±0.00150%	99.999195±0.000475%
LSUN Bedroom (256 × 256)	99.99953±0.00083%	99.99999±0.00013%	99.999994±0.000025%

Additionally, the variance of ∥ϵ ⋆ θ ∥² is relatively small.

$$\operatorname{Var}{p{t}(\mathbf{x}{t})} | \boldsymbol{\epsilon}{\boldsymbol{\theta}}^{\star} |{2} = \operatorname{Var}{p{t}(\mathbf{x}{t})} | \mathbb{E}{p{t0}(\mathbf{x}{0}|\mathbf{x}{t})} \boldsymbol{\epsilon} |{2} = \mathbb{E}{p{t}(\mathbf{x}{t})} | \mathbb{E}{p{t0}(\mathbf{x}{0}|\mathbf{x}{t})} \boldsymbol{\epsilon} |{2}^{2} - \left[ \mathbb{E}{p{t}(\mathbf{x}{t})} | \mathbb{E}{p{t0}(\mathbf{x}{0}|\mathbf{x}{t})} \boldsymbol{\epsilon} |_{2}^{2} \right]^{2}$$

$$\leq \mathbb{E}{p{t}(\mathbf{x}{t})} \mathbb{E}{p{t0}(\mathbf{x}{0}|\mathbf{x}{t})} | \boldsymbol{\epsilon} |{2}^{2} - (d - 2m) = \mathbb{E}{p{0}(\mathbf{x}{0})} \mathbb{E}{p{0t}(\mathbf{x}{0}|\mathbf{x}{t})} | \boldsymbol{\epsilon} |{2}^{2} - (d - 2m)$$

$$= d - (d - 2m)$$

$$= 2m$$

$$(47)$$

Therefore, the standard deviation of ∥ϵ ⋆ θ ^∥² is upper bounded by ^√ 2m. Since d ≫ m, we can conclude that in the optimal case, the magnitude of vector field is approximately constant, i.e., ∥ϵ ⋆ θ ∥² ≈ √ d. The total sampling trajectory length is PN−¹ ⁿ=0 (σtn+1 −σtⁿ )∥ϵ ⋆ θ (xtn+1 )∥² ≈ σT √ d. Empirical verification is provided in Figure 11.

A.6 Proof of Proposition 4

Proposition 4. In deterministic sampling with the Euler method, the sample likelihood is non-decreasing, i.e., ∀n ∈ [1, N], we have ph(xˆtn−¹ ) ≥ ph(xˆtⁿ ) and ph(rθ(xˆtⁿ )) ≥ ph(xˆtⁿ ) in terms of the Gaussian KDE ph(x) = (1/|I|) P ⁱ∈I N (x; yⁱ , h² I) with any positive bandwidth h, assuming all samples in the sampling trajectory satisfy d1(xˆtⁿ ) ≤ d2(xˆtⁿ ).

Proof We first prove that given a random vector v falling within a sphere centered at the optimal denoising output r ⋆ θ (xˆtⁿ ) with a radius of ∥r ⋆ θ (xˆtⁿ ) − xˆtⁿ ∥² , i.e., ∥r ⋆ θ (xˆtⁿ ) − xˆtⁿ ∥² ≥ ∥v∥² , the sample likelihood is non-decreasing from xˆtⁿ to r ⋆ θ (xˆtⁿ )+v, i.e., ph(r ⋆ θ (xˆtⁿ )+v) ≥ ph(xˆtⁿ ). Then, we provide two settings for v to finish the proof.

The increase of sample likelihood from xˆtⁿ to r ⋆ θ (xˆtⁿ ) + v in terms of ph(x) is

$$p{h}(r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v}) - p{h}(\hat{\mathbf{x}}{t{n}}) = \frac{1}{|\mathcal{I}|} \sum{i} \left[ \mathcal{N} \left( r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v}; \mathbf{y}{i}, h^{2} \mathbf{I} \right) - \mathcal{N} \left( \hat{\mathbf{x}}{t{n}}; \mathbf{y}_{i}, h^{2} \mathbf{I} \right) \right]$$

$$\stackrel{(i)}{\geq} \frac{1}{2h^{2}|\mathcal{I}|} \sum{i} \mathcal{N} \left( \hat{\mathbf{x}}{t{n}}; \mathbf{y}{i}, h^{2} \mathbf{I} \right) \left[ |\hat{\mathbf{x}}{t{n}} - \mathbf{y}{i}|{2}^{2} - |r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v} - \mathbf{y}{i}|_{2}^{2} \right]$$

$$= \frac{1}{2h^{2}|\mathcal{I}|} \sum{i} \mathcal{N} \left( \hat{\mathbf{x}}{t{n}}; \mathbf{y}{i}, h^{2} \mathbf{I} \right) \left[ |\hat{\mathbf{x}}{t{n}}|{2}^{2} - 2\hat{\mathbf{x}}{t{n}}^{T} \mathbf{y}{i} - |r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v}|{2}^{2} + 2 \left( r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v} \right)^{T} \mathbf{y}{i} \right]$$

$$\stackrel{(ii)}{=} \frac{1}{2h^{2}|\mathcal{I}|} \sum{i} \mathcal{N} \left( \hat{\mathbf{x}}{t{n}}; \mathbf{y}{i}, h^{2} \mathbf{I} \right) \left[ |\hat{\mathbf{x}}{t{n}}|{2}^{2} - 2\hat{\mathbf{x}}{t{n}}^{T} r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) - |r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v}|{2}^{2} + 2 \left( r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + \mathbf{v} \right)^{T} r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) \right]$$

$$= \frac{1}{2h^{2}|\mathcal{I}|} \sum{i} \mathcal{N} \left( \hat{\mathbf{x}}{t{n}}; \mathbf{y}{i}, h^{2} \mathbf{I} \right) \left[ |\hat{\mathbf{x}}{t{n}}|{2}^{2} - 2\hat{\mathbf{x}}{t{n}}^{T} r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) + |r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}})|{2}^{2} - |\mathbf{v}|_{2}^{2} \right]$$

$$= \frac{1}{2h^{2}|\mathcal{I}|} \sum{i} \mathcal{N} \left( \hat{\mathbf{x}}{t{n}}; \mathbf{y}{i}, h^{2} \mathbf{I} \right) \left[ |r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{t{n}}) - \hat{\mathbf{x}}{t{n}}|{2}^{2} - |\mathbf{v}|_{2}^{2} \right] \geq 0,$$

$$(48)$$

where (i) uses the definition of convex function f(x2) ≥ f(x1) + f ′ (x1)(x² − x1) with f(x) = exp − ∥x∥ , x¹ = (xˆtⁿ − yi) /h and x² = (r ⋆ θ (xˆtⁿ ) + v − yi) /h; (ii) uses the relationship between two consecutive steps xˆtⁿ and r ⋆ θ (xˆtⁿ ) in mean shift, i.e.,

$$r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}{tn}) = \mathbf{m}(\hat{\mathbf{x}}{tn}) = \sum{i} \frac{\exp\left(-|\hat{\mathbf{x}}_{t_n} - \mathbf{y}_i|2^2 / 2h^2\right)}{\sum{j} \exp\left(-|\hat{\mathbf{x}}_{t_n} - \mathbf{y}_j|_2^2 / 2h^2\right)} \mathbf{y}_i, \tag{49}$$

which implies the following equation also holds

$$\sum{i} \mathcal{N}\left(\hat{\mathbf{x}}{t_n}; \mathbf{y}_i, h^2 \mathbf{I}\right) \mathbf{x}i = \sum{i} \mathcal{N}\left(\hat{\mathbf{x}}_{t_n}; \mathbf{y}i, h^2 \mathbf{I}\right) r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}_{t_n}).$$ (50)

Since ∥r ⋆ θ (xˆtⁿ ) − xˆtⁿ ∥² ≥ ∥v∥² , or equivalently, ∥r ⋆ θ (xˆtⁿ ) − xˆtⁿ ∥ 2 ² ≥ ∥v∥ 2 2 , we conclude that the sample likelihood monotonically increases from xˆtⁿ to r ⋆ θ (xˆtⁿ )+v unless xˆtⁿ = r ⋆ θ (xˆtⁿ )+ v, in terms of the kernel density estimate ph(x) = ¹ |I| P ⁱ N (x; yⁱ , h² I) with any positive bandwidth h. We next provide two settings for v, which trivially satisfy the condition ∥r ⋆ θ (xˆtⁿ ) − xˆtⁿ ∥² ≥ ∥v∥² , and have the following corollaries:

• $$ph(r{\theta}(\hat{\mathbf{x}}_{t_n})) \ge ph(\hat{\mathbf{x}}{tn})$$ , when $\mathbf{v} = r{\theta}(\hat{\mathbf{x}}_{tn}) - r{\theta}^{\star}(\hat{\mathbf{x}}_{t_n})$ .

• $$ph(\hat{\mathbf{x}}{t_{n-1}}) \ge ph(\hat{\mathbf{x}}{tn})$$ , when $\mathbf{v} = r{\boldsymbol{\theta}}(\hat{\mathbf{x}}_{tn}) - r{\boldsymbol{\theta}}^{\star}(\hat{\mathbf{x}}_{tn}) + \frac{t{n-1}}{tn} (\hat{\mathbf{x}}{tn} - r{\boldsymbol{\theta}}(\hat{\mathbf{x}}_{t_n}))$ .

Figure 12: Trajectory deviation (1-D projection) on unconditional generation.

Figure 13: Trajectory reconstruction, visualization and statistics on FFHQ.

Appendix B. Additional Results

B.1 Visualization of Sampling Trajectories

Figures 12a and 12b provide more experiments about the 1-D trajectory projection on FFHQ and LSUN Bedroom. Figure 13 provide more results about Multi-D projections on FFHQ. Figure 14 visualizes more generated samples on three datasets.

B.2 Diagnosis of Score Deviation

In this section, we simulate four new trajectories based on the optimal denoising output r ⋆ θ (·) to monitor the score deviation from the optimum. We denote the optimal sampling

Figure 14: The visual comparison of samples generated by DDIM and DDIM + GITS (1st row: CIFAR-10, 2nd row: ImageNet 64 × 64, 3rd row: LSUN Bedroom).

trajectory as {xˆ ⋆ tn N ⁿ=0, where we generate samples as the standard sampling trajectory {xˆtⁿ } N ⁿ=0 with the same time schedule Γ = {t⁰ = ϵ, · · · , t^N = T}, but adopt optimal denoising output r ⋆ θ (·) rather than denoising output rθ(·) for score estimation. The other three trajectories are simulated by tracking the (optimal) denoising output of each sample in {xˆ ⋆ ^t } or {xˆt}, and designated as {rθ(xˆ ⋆ t )}, {r ⋆ θ (xˆ ⋆ t )}, {r ⋆ θ (xˆt)}. According to (19) and t⁰ = 0, we have xˆ ⋆ ^t⁰ = r ⋆ θ (xˆ ⋆ t1 ), and similarly, xˆt⁰ = rθ(xˆt¹ ). As t → 0, r ⋆ θ (xˆ ⋆ t ) and r ⋆ θ (xˆt) serve as the approximate nearest neighbors of xˆ ⋆ ^t and xˆ^t to the real data, respectively.

We calculate the deviation of denoising output to quantify the score deviation across all time steps using the L ² distance, though they should differ by a factor t 2 , and have the following observation: The learned score is well-matched to the optimal score in the largenoise region, otherwise, they may diverge or almost coincide depending on different regions. In fact, our learned score has to moderately diverge from the optimum to guarantee the generative ability. Otherwise, the ODE-based sampling reduces to an approximate (singlestep) annealed mean shift for global mode-seeking (see Section 4.2), and simply replays the dataset. As shown in Figure 15, the nearest sample of xˆ ⋆ t0 to the real data is almost the same as itself, which indicates the optimal sampling trajectory has a very limited ability to synthesize novel samples. Empirically, score deviation in a small region is sufficient to bring forth a decent generative ability.

From the comparison of {rθ(xˆ ⋆ t )}, {r ⋆ θ (xˆ ⋆ t )} sequences in Figures 15 and 16, we can clearly see that along the optimal sampling trajectory, the deviation between the learned denoising output rθ(·) and its optimal counterpart r ⋆ θ (·) behaves differently in three successive regions: The deviation starts off as almost negligible (about 10 < t ≤ 80), gradually increases (about 3 < t ≤ 10), and then drops down to a low level once again (about 0 ≤ t ≤ 3). This phenomenon was also validated by a recent work Xu et al. (2023) with a different perspective. We further observe that along the sampling trajectory, this

Figure 15: Top: We visualize a forward diffusion process of a randomly-selected image to obtain its encoding xˆt^N (first row) and simulate multiple trajectories starting from this encoding (other rows). Bottom: The k-nearest neighbors (k=5) of xˆt⁰ and xˆ ⋆ t0 to real samples in the dataset.

Figure 16: The deviation (Euclidean distance) of outputs from their corresponding optima.

phenomenon disappears and the score deviation keeps increasing (see {rθ(xˆt)}, {r ⋆ θ (xˆt)} sequences in Figures 15 and 16). Additionally, samples in the latter half of {r ⋆ θ (xˆt)} appear almost the same as the nearest sample of xˆt⁰ to the real data, as shown in Figure 15. This indicates that our score-based model strives to explore novel regions, and synthetic samples in the sampling trajectory are quickly attracted to a real-data mode but do not fall into it.

B.3 More Results

Time schedule. The uniform schedule is commonly used with the DDPM (Ho et al., 2020) backbone. Following EDMs (Karras et al., 2022), we rewrite this schedule from its original range [ϵs, 1] to [t0, t^N ], where ϵ^s = 0.001, t⁰ = 0.002 and t^N = 80. We first uniformly sample τⁿ (n ∈ [0, N]) from [ϵs, 1] and then calculate tⁿ by tⁿ = q exp( ¹ βdτⁿ ² + βminτn) − 1, where β^d = 2 ϵs−1 log(1+t⁰ ϵs − log(1 + t^N 2 ), βmin = log(1 + t^N 2 ) − 1 βd.

Table 9: Sample quality in terms of FID (Heusel et al., 2017) on four datasets (resolutions ranging from 32 × 32 to 256 × 256). †: After obtaining the DP schedule, we could further optimize the first time step with AFS, using the same "warmup" samples. The default setting in our main submission does not use AFS and keeps the coefficient in dynamic programming as 1.1 for LSUN Bedroom and 1.15 otherwise. Although the performance can be further improved by carefully tuning the coefficient and using AFS as shown below.

METHOD	Coeff				NFE
		3	4	5	6	7	8	9	10
CIFAR-10 32×32 (Krizhevsky and Hinton, 2009)
DDIM (Song et al., 2021a)	-	93.36	66.76	49.66	35.62	27.93	22.32	18.43	15.69
DDIM + GITS	1.10	88.68	46.88	32.50	22.04	16.76	13.93	11.57	10.09
DDIM + GITS (default)	1.15	79.67	43.07	28.05	21.04	16.35	13.30	11.62	10.37
DDIM + GITS	1.20	77.22	43.16	29.06	22.69	18.91	14.22	12.03	11.38
iPNDM (Zhang and Chen, 2023)	-	47.98	24.82	13.59	7.05	5.08	3.69	3.17	2.77
iPNDM + GITS	1.10	51.31	17.19	12.90	5.98	6.62	4.36	3.59	3.14
iPNDM + GITS (default)	1.15	43.89	15.10	8.38	4.88	5.11	3.24	2.70	2.49
iPNDM + GITS	1.20	42.06	15.85	9.33	7.13	5.95	3.28	2.81	2.71
iPNDM + GITS	1.10	34.22	11.99	12.44	6.08	6.20	3.53	3.48	2.91
iPNDM + GITS	1.15	29.63	11.23	8.08	4.86	4.46	2.92	2.46	2.27
iPNDM + GITS	1.20	25.98	10.11	6.77	4.29	3.43	2.70	2.42	2.28
FFHQ 64×64 (Karras et al., 2019)
DDIM (Song et al., 2021a)	-	78.21	57.48	43.93	35.22	28.86	24.39	21.01	18.37
DDIM + GITS	1.10	62.70	43.12	31.01	24.62	20.35	17.19	14.71	13.01
DDIM + GITS (default)	1.15	60.84	40.81	29.80	23.67	19.41	16.60	14.46	13.06
DDIM + GITS	1.20	59.64	40.56	30.29	23.88	20.07	17.36	15.40	14.05
iPNDM (Zhang and Chen, 2023)	-	45.98	28.29	17.17	10.03	7.79	5.52	4.58	3.98
iPNDM + GITS	1.10	34.82	18.75	13.07	7.79	8.30	4.76	5.36	3.47
iPNDM + GITS (default)	1.15	33.09	17.04	11.22	7.00	6.72	4.52	4.33	3.62
iPNDM + GITS	1.20	31.70	16.87	10.83	7.10	6.37	5.78	4.81	4.39
iPNDM + GITS	1.10	33.19	19.88	12.90	8.29	7.50	4.26	4.95	3.13
iPNDM + GITS	1.15	30.39	15.78	10.15	6.86	5.97	4.09	3.76	3.24
iPNDM + GITS	1.20	26.41	13.59	8.85	6.39	5.36	4.91	3.89	3.51
ImageNet 64×64 (Russakovsky et al., 2015)
DDIM (Song et al., 2021a)	-	82.96	58.43	43.81	34.03	27.46	22.59	19.27	16.72
DDIM + GITS	1.10	60.11	36.23	27.31	20.82	16.41	14.16	11.95	10.84
DDIM + GITS (default)	1.15	57.06	35.07	24.92	19.54	16.01	13.79	12.17	10.83
DDIM + GITS	1.20	54.24	34.27	24.67	19.46	16.66	14.15	13.41	11.87
iPNDM (Zhang and Chen, 2023)	-	58.53	33.79	18.99	12.92	9.17	7.20	5.91	5.11
iPNDM + GITS	1.10	36.18	19.64	13.18	9.58	7.68	6.44	5.24	4.59
iPNDM + GITS (default)	1.15	34.47	18.95	10.79	8.43	6.83	5.82	4.96	4.48
iPNDM + GITS	1.20	32.70	18.59	11.04	9.23	7.18	6.20	5.50	5.08
iPNDM + GITS	1.10	31.50	21.50	13.73	10.74	7.99	6.88	5.29	4.64
iPNDM + GITS	1.15	28.01	18.28	10.28	8.68	6.76	5.90	4.81	4.40
iPNDM + GITS	1.20	26.41	16.41	9.85	8.39	6.44	5.64	4.79	4.47
LSUN Bedroom 256×256 (Yu et al., 2015) (pixel-space)
DDIM (Song et al., 2021a)	-	86.13	54.45	34.34	25.25	19.49	15.71	13.26	11.42
DDIM + GITS	1.05	81.77	36.89	27.46	18.78	13.60	12.23	10.29	8.77
DDIM + GITS (default)	1.10	61.85	35.12	22.04	16.54	13.58	11.20	9.82	9.04
DDIM + GITS	1.15	60.11	31.02	23.65	17.18	13.42	12.61	10.89	10.57
iPNDM (Zhang and Chen, 2023)	-	80.99	43.90	26.65	20.73	13.80	11.78	8.38	5.57
iPNDM + GITS	1.05	59.02	24.71	19.08	12.77	8.19	6.67	5.58	4.83
iPNDM + GITS (default)	1.10	45.75	22.98	15.85	10.41	8.63	7.31	6.01	5.28
iPNDM + GITS	1.15	44.78	21.67	17.29	11.52	9.59	8.82	7.22	5.97

The logSNR time schedule is proposed for fast sampling in DPM-Solver (Lu et al., 2022a). We first uniformly sample λⁿ (n ∈ [0, N]) from [λmin, λmax] where λmin = − log t^N and λmax = − log t0. The logSNR schedule is given by tⁿ = e ⁻λⁿ .

Table 10: Comparison of various time schedules on CIFAR-10.

NFE	TIME SCHEDULE	FID
Uniform
3	[80.0000, 6.9503, 1.2867, 0.0020]	50.44
4	[80.0000, 11.7343, 2.8237, 0.8565, 0.0020]	18.73
5	[80.0000, 16.5063, 4.7464, 1.7541, 0.6502, 0.0020]	17.34
6	[80.0000, 20.9656, 6.9503, 2.8237, 1.2867, 0.5272, 0.0020]	9.75
7	[80.0000, 25.0154, 9.3124, 4.0679, 2.0043, 1.0249, 0.4447, 0.0020]	12.50
8	[80.0000, 28.6496, 11.7343, 5.4561, 2.8237, 1.5621, 0.8565, 0.3852, 0.0020]	7.56
9	[80.0000, 31.8981, 14.1472, 6.9503, 3.7419, 2.1599, 1.2867, 0.7382, 0.3401, 0.0020]	10.60
10	[80.0000, 34.8018, 16.5063, 8.5141, 4.7464, 2.8237, 1.7541, 1.0985, 0.6502, 0.3047, 0.0020]	7.35
LogSNR
3	[80.0000, 2.3392, 0.0684, 0.0020]	88.38
4	[80.0000, 5.6569, 0.4000, 0.0283, 0.0020]	35.59
5	[80.0000, 9.6090, 1.1542, 0.1386, 0.0167, 0.0020]	19.87
6	[80.0000, 13.6798, 2.3392, 0.4000, 0.0684, 0.0117, 0.0020]	10.68
7	[80.0000, 17.6057, 3.8745, 0.8527, 0.1876, 0.0413, 0.0091, 0.0020]	6.56
8	[80.0000, 21.2732, 5.6569, 1.5042, 0.4000, 0.1064, 0.0283, 0.0075, 0.0020]	4.74
9	[80.0000, 24.6462, 7.5929, 2.3392, 0.7207, 0.2220, 0.0684, 0.0211, 0.0065, 0.0020]	3.53
10	[80.0000, 27.7258, 9.6090, 3.3302, 1.1542, 0.4000, 0.1386, 0.0480, 0.0167, 0.0058, 0.0020]	2.94
	Polynomial (ρ = 7)
3	[80.0000, 9.7232, 0.4700, 0.0020]	47.98
4	[80.0000, 17.5278, 2.5152, 0.1698, 0.0020]	24.82
5	[80.0000, 24.4083, 5.8389, 0.9654, 0.0851, 0.0020]	13.59
6	[80.0000, 30.1833, 9.7232, 2.5152, 0.4700, 0.0515, 0.0020]	7.05
7	[80.0000, 34.9922, 13.6986, 4.6371, 1.2866, 0.2675, 0.0352, 0.0020]	5.08
8	[80.0000, 39.0167, 17.5278, 7.1005, 2.5152, 0.7434, 0.1698, 0.0261, 0.0020]	3.69
9	[80.0000, 42.4152, 21.1087, 9.7232, 4.0661, 1.5017, 0.4700, 0.1166, 0.0204, 0.0020]	3.17
10	[80.0000, 45.3137, 24.4083,12.3816, 5.8389, 2.5152, 0.9654, 0.3183, 0.0851, 0.0167, 0.0020]	2.77
	GITS (ours)
3	[80.0000, 3.8811, 0.9654, 0.0020]	43.89
4	[80.0000, 5.8389, 1.8543, 0.4700, 0.0020]	15.10
5	[80.0000, 6.6563, 2.1632, 0.8119, 0.2107, 0.0020]	8.38
6	[80.0000, 10.9836, 3.8811, 1.5840, 0.5666, 0.1698, 0.0020]	4.88
7	[80.0000, 12.3816, 3.8811, 1.5840, 0.5666, 0.1698, 0.0395, 0.0020]	3.76
8	[80.0000, 10.9836, 3.8811, 1.8543, 0.9654, 0.4700, 0.2107, 0.0665, 0.0020]	3.24
9	[80.0000, 12.3816, 4.4590, 2.1632, 1.1431, 0.5666, 0.2597, 0.1079, 0.0300, 0.0020]	2.70
10	[80.0000, 12.3816, 4.4590, 2.1632, 1.1431, 0.5666, 0.3183, 0.1698, 0.0665, 0.0225, 0.0020]	2.49

The polynomial time schedule tⁿ = (t 1/ρ ⁰ + n N (t 1/ρ ^N − t 1/ρ 0 ))^ρ is proposed in EDM (Karras et al., 2022), where t⁰ = 0.002, t^N = 80, n ∈ [0, N], and ρ = 7.

The optimized time schedules for SDv1.5 in Figure 10 include

AYS (Sabour et al., 2024): [999, 850, 736, 645, 545, 455, 343, 233, 124, 24, 0].
GITS: [999, 783, 632, 483, 350, 233, 133, 67, 33, 17, 0].

References

G. Alain and Y. Bengio. What regularized auto-encoders learn from the data-generating distribution. Journal of Machine Learning Research, 15(1):3563–3593, 2014.
M. S. Albergo, N. M. Boffi, and E. Vanden-Eijnden. Stochastic interpolants: A unifying framework for flows and diffusions. arXiv preprint arXiv:2303.08797, 2023.
B. D. Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3):313–326, 1982.
J. Antognini and J. Sohl-Dickstein. Pca of high dimensional random walks with comparison to neural network training. In Advances in Neural Information Processing Systems, volume 31, 2018.
Y. Balaji, S. Nah, X. Huang, A. Vahdat, J. Song, K. Kreis, M. Aittala, T. Aila, S. Laine, B. Catanzaro, et al. ediffi: Text-to-image diffusion models with an ensemble of expert denoisers. arXiv preprint arXiv:2211.01324, 2022.
F. Bao, C. Li, J. Zhu, and B. Zhang. Analytic-dpm: an analytic estimate of the optimal reverse variance in diffusion probabilistic models. In International Conference on Learning Representations, 2022.
Y. Bengio, A. Courville, and P. Vincent. Representation learning: A review and new perspectives. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(8): 1798–1828, 2013a.
Y. Bengio, L. Yao, G. Alain, and P. Vincent. Generalized denoising auto-encoders as generative models. In Advances in Neural Information Processing Systems, pages 899– 907, 2013b.
G. Biroli, T. Bonnaire, V. De Bortoli, and M. M´ezard. Dynamical regimes of diffusion models. Nature Communications, 15(1):9957, 2024.
A. Blattmann, R. Rombach, H. Ling, T. Dockhorn, S. W. Kim, S. Fidler, and K. Kreis. Align your latents: High-resolution video synthesis with latent diffusion models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 22563– 22575, 2023.
M. A. Carreira-Perpin´an. A review of mean-shift algorithms for clustering. arXiv preprint arXiv:1503.00687, 2015.
D. Chen, Z. Zhou, J.-P. Mei, C. Shen, C. Chen, and C. Wang. A geometric perspective on diffusion models. arXiv preprint arXiv:2305.19947, 2023a.
D. Chen, Z. Zhou, C. Wang, C. Shen, and S. Lyu. On the trajectory regularity of ode-based diffusion sampling. In International Conference on Machine Learning, pages 7905–7934, 2024.
M. Chen, K. Huang, T. Zhao, and M. Wang. Score approximation, estimation and distribution recovery of diffusion models on low-dimensional data. In International Conference on Machine Learning, pages 4672–4712, 2023b.
S. Chen, G. Daras, and A. Dimakis. Restoration-degradation beyond linear diffusions: A non-asymptotic analysis for ddim-type samplers. In International Conference on Machine Learning, pages 4462–4484, 2023c.
T. Chen, G.-H. Liu, and E. Theodorou. Likelihood training of schr¨odinger bridge using forward-backward sdes theory. In International Conference on Learning Representations, 2022.
Y. Cheng. Mean shift, mode seeking, and clustering. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(8):790–799, 1995.
D. Comaniciu and P. Meer. Mean shift analysis and applications. In Proceedings of the International Conference on Computer Vision, pages 1197–1203, 1999.
D. Comaniciu and P. Meer. Mean shift: A robust approach toward feature space analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(5):603–619, 2002.
D. Comaniciu, V. Ramesh, and P. Meer. Real-time tracking of non-rigid objects using mean shift. In Proceedings IEEE Conference on Computer Vision and Pattern Recognition, pages 142–149, 2000.
D. Comaniciu, V. Ramesh, and P. Meer. Kernel-based object tracking. IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(5):564–577, 2003.
T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to algorithms. MIT press, 2022.
V. De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022.
P. Dhariwal and A. Nichol. Diffusion models beat gans on image synthesis. In Advances in Neural Information Processing Systems, pages 8780–8794, 2021.
M. P. Do Carmo. Differential geometry of curves and surfaces: revised and updated second edition. Courier Dover Publications, 2016.
T. Dockhorn, A. Vahdat, and K. Kreis. Genie: Higher-order denoising diffusion solvers. In Advances in Neural Information Processing Systems, pages 30150–30166, 2022.
B. Efron. Large-scale inference: empirical Bayes methods for estimation, testing, and prediction. Cambridge University Press, 2010.
B. Efron. Tweedie's formula and selection bias. Journal of the American Statistical Association, 106(496):1602–1614, 2011.
P. Esser, S. Kulal, A. Blattmann, R. Entezari, J. M¨uller, H. Saini, Y. Levi, D. Lorenz, A. Sauer, F. Boesel, et al. Scaling rectified flow transformers for high-resolution image synthesis. In International Conference on Machine Learning, 2024.
W. Feller. On the theory of stochastic processes, with particular reference to applications. In Proceedings of the First Berkeley Symposium on Mathematical Statistics and Probability, pages 403–432, 1949.
E. Frankel, S. Chen, J. Li, P. W. Koh, L. J. Ratliff, and S. Oh. S4s: Solving for a diffusion model solver. In International Conference on Machine Learning, 2025.
K. Fukunaga and L. Hostetler. The estimation of the gradient of a density function, with applications in pattern recognition. IEEE Transactions on Information Theory, 21(1): 32–40, 1975.
G. H. Golub and C. F. Van Loan. Matrix computations, volume 3. The Johns Hopkins University press, 2013.
E. Heitz, L. Belcour, and T. Chambon. Iterative α-(de) blending: A minimalist deterministic diffusion model. In ACM SIGGRAPH 2023 Conference Proceedings, pages 1–8, 2023.
M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, and S. Hochreiter. GANs trained by a two time-scale update rule converge to a local Nash equilibrium. In Advances in Neural Information Processing Systems, pages 6626–6637, 2017.
J. Ho, A. Jain, and P. Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, pages 6840–6851, 2020.
J. Ho, T. Salimans, A. A. Gritsenko, W. Chan, M. Norouzi, and D. J. Fleet. Video diffusion models. In Advances in Neural Information Processing Systems, pages 8633–8646, 2022.
R. Huang, J. Huang, D. Yang, Y. Ren, L. Liu, M. Li, Z. Ye, J. Liu, X. Yin, and Z. Zhao. Make-an-audio: Text-to-audio generation with prompt-enhanced diffusion models. In International Conference on Machine Learning, pages 13916–13932. PMLR, 2023.
J. R. Hurley and R. B. Cattell. The procrustes program: Producing direct rotation to test a hypothesized factor structure. Behavioral science, 7(2):258, 1962.
A. Hyv¨arinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6:695–709, 2005.
T. Karras, S. Laine, and T. Aila. A style-based generator architecture for generative adversarial networks. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4401–4410, 2019.
T. Karras, M. Aittala, T. Aila, and S. Laine. Elucidating the design space of diffusionbased generative models. In Advances in Neural Information Processing Systems, pages 26565–26577, 2022.
D. Kim, S. Shin, K. Song, W. Kang, and I.-C. Moon. Soft truncation: A universal training technique of score-based diffusion model for high precision score estimation. In International Conference on Machine Learning, pages 11201–11228, 2022.
D. Kim, C.-H. Lai, W.-H. Liao, N. Murata, Y. Takida, T. Uesaka, Y. He, Y. Mitsufuji, and S. Ermon. Consistency trajectory models: Learning probability flow ode trajectory of diffusion. In International Conference on Learning Representations, 2024.
D. P. Kingma, T. Salimans, B. Poole, and J. Ho. Variational diffusion models. In Advances in Neural Information Processing Systems, pages 21696–21707, 2021.
Y. Kirstain, A. Polyak, U. Singer, S. Matiana, J. Penna, and O. Levy. Pick-a-pic: An open dataset of user preferences for text-to-image generation. In Advances in Neural Information Processing Systems, volume 36, pages 36652–36663, 2023.
Z. Kong, W. Ping, J. Huang, K. Zhao, and B. Catanzaro. Diffwave: A versatile diffusion model for audio synthesis. In International Conference on Learning Representations, 2021.
A. Krizhevsky and G. Hinton. Learning multiple layers of features from tiny images. Technical Report, 2009.
M. Kwon, J. Jeong, and Y. Uh. Diffusion models already have a semantic latent space. In International Conference on Learning Representations, 2023.
A. Land and A. Doig. An automatic method of solving discrete programming problems. Econometrica, 28(3):497–520, 1960.
H. Lee, J. Lu, and Y. Tan. Convergence of score-based generative modeling for general data distributions. In International Conference on Algorithmic Learning Theory, pages 946–985, 2023.
T. Lewiner, J. D. Gomes Jr, H. Lopes, and M. Craizer. Curvature and torsion estimators based on parametric curve fitting. Computers & Graphics, 29(5):641–655, 2005.
T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Doll´ar, and C. L. Zitnick. Microsoft coco: Common objects in context. In Proceedings of the European Conference on Computer Vision, pages 740–755, 2014.
Y. Lipman, R. T. Chen, H. Ben-Hamu, M. Nickel, and M. Le. Flow matching for generative modeling. In International Conference on Learning Representations, 2023.
G.-H. Liu, A. Vahdat, D.-A. Huang, E. A. Theodorou, W. Nie, and A. Anandkumar. I2sb: image-to-image schr¨odinger bridge. In International Conference on Machine Learning, pages 22042–22062, 2023a.
L. Liu, Y. Ren, Z. Lin, and Z. Zhao. Pseudo numerical methods for diffusion models on manifolds. In International Conference on Learning Representations, 2022.
X. Liu, C. Gong, and Q. Liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. In International Conference on Learning Representations, 2023b.
C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, and J. Zhu. Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps. In Advances in Neural Information Processing Systems, pages 5775–5787, 2022a.
C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, and J. Zhu. Dpm-solver++: Fast solver for guided sampling of diffusion probabilistic models. arXiv preprint arXiv:2211.01095, 2022b.
E. Luhman and T. Luhman. Knowledge distillation in iterative generative models for improved sampling speed. arXiv preprint arXiv:2101.02388, 2021.
S. Lyu. Interpretation and generalization of score matching. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence, pages 359–366, 2009.
J. Moore, H. Ahmed, and R. Antia. High dimensional random walks can appear low dimensional: Application to influenza h3n2 evolution. Journal of theoretical biology, 447:56–64, 2018.
C. N. Morris. Parametric empirical bayes inference: theory and applications. Journal of the American Statistical Association, 78(381):47–55, 1983.
K. Neklyudov, R. Brekelmans, D. Severo, and A. Makhzani. Action matching: Learning stochastic dynamics from samples. In International Conference on Machine Learning, pages 25858–25889. PMLR, 2023.
A. Q. Nichol and P. Dhariwal. Improved denoising diffusion probabilistic models. In International Conference on Machine Learning, pages 8162–8171, 2021.
A. Q. Nichol, P. Dhariwal, A. Ramesh, P. Shyam, P. Mishkin, B. Mcgrew, I. Sutskever, and M. Chen. Glide: Towards photorealistic image generation and editing with text-guided diffusion models. In International Conference on Machine Learning, pages 16784–16804, 2022.
B. Oksendal. Stochastic differential equations: an introduction with applications. Springer Science & Business Media, 2013.
W. Peebles and S. Xie. Scalable diffusion models with transformers. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 4195–4205, 2023.
J. Pidstrigach. Score-based generative models detect manifolds. In Advances in Neural Information Processing Systems, pages 35852–35865, 2022.
D. Podell, Z. English, K. Lacey, A. Blattmann, T. Dockhorn, J. M¨uller, J. Penna, and R. Rombach. Sdxl: Improving latent diffusion models for high-resolution image synthesis. In International Conference on Learning Representations, 2024.
A. Ramesh, P. Dhariwal, A. Nichol, C. Chu, and M. Chen. Hierarchical text-conditional image generation with clip latents. arXiv preprint arXiv:2204.06125, 2022.
M. Raphan and E. P. Simoncelli. Least squares estimation without priors or supervision. Neural Computation, 23(2):374–420, 2011.
H. E. Robbins. An empirical bayes approach to statistics. In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1956.
R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer. High-resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684–10695, 2022.
N. Ruiz, Y. Li, V. Jampani, Y. Pritch, M. Rubinstein, and K. Aberman. Dreambooth: Fine tuning text-to-image diffusion models for subject-driven generation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 22500– 22510, 2023.
O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. S. Bernstein, A. C. Berg, and F. Li. Imagenet large scale visual recognition challenge. International Journal of Computer Vision, 115(3):211–252, 2015.
A. Sabour, S. Fidler, and K. Kreis. Align your steps: Optimizing sampling schedules in diffusion models. In International Conference on Machine Learning, 2024.
C. Saharia, W. Chan, S. Saxena, L. Li, J. Whang, E. L. Denton, K. Ghasemipour, R. Gontijo Lopes, B. Karagol Ayan, T. Salimans, et al. Photorealistic text-to-image diffusion models with deep language understanding. In Advances in Neural Information Processing Systems, pages 36479–36494, 2022.
T. Salimans and J. Ho. Progressive distillation for fast sampling of diffusion models. In International Conference on Learning Representations, 2022.
S. Saremi and A. Hyv¨arinen. Neural empirical bayes. Journal of Machine Learning Research, 20:181:1–181:23, 2019.
S. S¨arkk¨a and A. Solin. Applied stochastic differential equations, volume 10. Cambridge University Press, 2019.
P. H. Sch¨onemann. A generalized solution of the orthogonal procrustes problem. Psychometrika, 31:1–10, 1966.
C. Shen, M. J. Brooks, and A. van den Hengel. Fast global kernel density mode seeking with application to localisation and tracking. In Proceedings of the International Conference on Computer Vision, pages 1516–1523, 2005.
B. W. Silverman. Using kernel density estimates to investigate multimodality. Journal of the Royal Statistical Society: Series B (Methodological), 43(1):97–99, 1981.
J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference on Machine Learning, pages 2256–2265, 2015.
J. Song, C. Meng, and S. Ermon. Denoising diffusion implicit models. In International Conference on Learning Representations, 2021a.
Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. In Advances in Neural Information Processing Systems, pages 11895–11907, 2019.
Y. Song and S. Ermon. Improved techniques for training score-based generative models. In Advances in Neural Information Processing Systems, pages 12438–12448, 2020.
Y. Song, S. Garg, J. Shi, and S. Ermon. Sliced score matching: A scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence, pages 574–584. PMLR, 2019.
Y. Song, C. Durkan, I. Murray, and S. Ermon. Maximum likelihood training of scorebased diffusion models. In Advances in Neural Information Processing Systems, pages 1415–1428, 2021b.
Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations, 2021c.
Y. Song, P. Dhariwal, M. Chen, and I. Sutskever. Consistency models. In International Conference on Machine learning, pages 32211–32252, 2023.
R. Stratonovich. Conditional Markov processes and their application to the theory of optimal control. American Elsevier Publishing, 1968.
R. S. Sutton, A. G. Barto, et al. Reinforcement learning: An introduction, volume 1. MIT press Cambridge, 1998.
V. Tong, T.-D. Hoang, A. Liu, G. V. d. Broeck, and M. Niepert. Learning to discretize denoising diffusion odes. In International Conference on Learning Representations, 2025.
A. Vahdat, K. Kreis, and J. Kautz. Score-based generative modeling in latent space. In Advances in Neural Information Processing Systems, pages 11287–11302, 2021.
R. Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018.
P. Vincent. A connection between score matching and denoising autoencoders. Neural Computation, 23(7):1661–1674, 2011.
P. Vincent, H. Larochelle, Y. Bengio, and P.-A. Manzagol. Extracting and composing robust features with denoising autoencoders. In International Conference on Machine learning, pages 1096–1103, 2008.
B. Wang and J. Vastola. The unreasonable effectiveness of gaussian score approximation for diffusion models and its applications. Transactions on Machine Learning Research, 2024.
D. Watson, J. Ho, M. Norouzi, and W. Chan. Learning to efficiently sample from diffusion probabilistic models. arXiv preprint arXiv:2106.03802, 2021.
H. Whitney. Differentiable manifolds. Annals of Mathematics, pages 645–680, 1936.
E. Xie, J. Chen, J. Chen, H. Cai, H. Tang, Y. Lin, Z. Zhang, M. Li, L. Zhu, Y. Lu, et al. Sana: Efficient high-resolution image synthesis with linear diffusion transformers. In International Conference on Learning Representations, 2025.
Y. Xu, S. Tong, and T. S. Jaakkola. Stable target field for reduced variance score estimation in diffusion models. In International Conference on Learning Representations, 2023.
R. Yamasaki and T. Tanaka. Properties of mean shift. IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(9):2273–2286, 2020.
H. D. Young, R. A. Freedman, T. Sandin, and A. L. Ford. University physics, volume 9. Addison-wesley Reading, MA, 1996.
F. Yu, A. Seff, Y. Zhang, S. Song, T. Funkhouser, and J. Xiao. Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop. arXiv preprint arXiv:1506.03365, 2015.
Q. Zhang and Y. Chen. Diffusion normalizing flow. Advances in Neural Information Processing Systems, pages 16280–16291, 2021.
Q. Zhang and Y. Chen. Fast sampling of diffusion models with exponential integrator. In International Conference on Learning Representations, 2023.
Q. Zhang, J. Song, and Y. Chen. Improved order analysis and design of exponential integrator for diffusion models sampling. arXiv preprint arXiv:2308.02157, 2023.
W. Zhao, L. Bai, Y. Rao, J. Zhou, and J. Lu. Unipc: A unified predictor-corrector framework for fast sampling of diffusion models. In Advances in Neural Information Processing Systems, pages 49842–49869, 2023.
H. Zheng, W. Nie, A. Vahdat, K. Azizzadenesheli, and A. Anandkumar. Fast sampling of diffusion models via operator learning. In International Conference on Machine Learning, pages 42390–42402, 2023.
Z. Zhou, D. Chen, C. Wang, and C. Chen. Fast ode-based sampling for diffusion models in around 5 steps. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7777–7786, 2024a.
Z. Zhou, D. Chen, C. Wang, C. Chen, and S. Lyu. Simple and fast distillation of diffusion models. In Advances in Neural Information Processing Systems, volume 37, pages 40831– 40860, 2024b.

chengeometricregularitydeterministicsampling2025.md 191 KB Historique Raw