1 Introduction↩︎

Denoising Score Matching (DSM) and its variants [1]–[4] have driven the success of score-based diffusion models [5]–[8] by providing a scalable score-estimation objective. To approximate an unknown target distribution \(p^\star\), a score-based diffusion model transports it to an easy-to-sample prior distribution (typically, a standard Gaussian), from which samples are brought back, via a reverse diffusion, towards \(p^\star\). To implement this reverse diffusion, one must estimate the unknown scores \(\nabla \log p^\star_s\) of the time-\(s\) noised marginals of the target distribution along the reverse path (precise definitions in Section 2). DSM produces such estimators by minimizing the time-integrated \(L^2\) error from the true unknown score. The accuracy of this estimation critically determines the sampling quality of the diffusion model [9], as illustrated by the following fundamental upper bound (up to constants) on the KL divergence between the target \(p^\star\) and the learned distribution \(\hat{p}\): \[\label{eq:intro-upper-bound} \mathrm{KL}(p^\star\,\| \,\hat{p}) \;\lesssim\; \frac{1}{2}\int_{0}^T \mathbb{E}\!\left[\| \mathbf{e}_s\|^2_2\right]\,\mathrm{d}s,\qquad{(1)}\] where \(\mathbf{e}_s = \hat{\mathbf{s}}_{s}-\nabla \log p^\star_s\) is the time-\(s\) score estimation error, i.e. the difference between the estimated score \(\hat{\mathbf{s}}_{s}\) and the true score \(\nabla \log p^\star_s\). Inequality ?? underpins much of the current theory of diffusion models [10], from convergence analyses [9], [11]–[13], to minimax optimality results [14]–[17], and interpretations of training objectives in terms of sample quality [18], [19].

Inequality ?? ties the full \(L^2\) score error to the sampling quality of the diffusion model. However, we show that a learned diffusion model can incur arbitrarily large \(L^2\) score error while perfectly matching the target distribution (Theorem 1 (ii)): \[\label{eq:surprise} \mathrm{KL}(p^\star\,\|\, \hat{p}) = 0 \qquad \text{ while } \qquad \int_{0}^T \mathbb{E}\!\left[\|\mathbf{e}_s\|^2_2\right]\,\mathrm{d}s \gg 0.\tag{1}\]

More strongly, no monotone function of the \(L^2\) score error can uniformly lower-bound any meaningful divergence between \(\hat{p}\) and \(p^\star\) (Theorem 1 (iii)). The \(L^2\) error is therefore not the intrinsic quantity controlling marginal sampling quality, raising a natural question:

We identify the geometric reason behind observation 1 (discussed in Section 3): the marginal distributions of a diffusion process solve a Fokker-Planck equation [20], which depends on the score error \(\mathbf{e}_s\) only through a weighted divergence. Decomposing the error vector field into gradient and solenoidal components, we observe that this weighted divergence filters out the solenoidal component of the error field, making it structurally invisible to the marginal dynamics, while the gradient component affects marginals (see Figure 1). Therefore, a model with purely solenoidal errors satisfies 1 , as illustrated in the middle panel of Figure 1. Decomposing errors into gradient and solenoidal components (a Helmholtz–Hodge decomposition [21]) is standard in scientific machine learning [22]–[26]. However, this framework has not yet been used to characterize the geometry of diffusion models, as done in this work.

The discussion above suggests that, although DSM reduces the gradient component of the score error by minimizing the full \(L^2\) objective, this overall score matching error also penalizes solenoidal components that do not affect marginals. Therefore, it is not always a reliable indicator of sampling quality (as in the first two panels of Figure 1). Instead, considering only the \(L^2\) norm of the observable gradient component of errors (denoted \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\)), we derive a new bound on sampling quality (Theorem 2): \[\label{eq:kl-upper-new-intro} \mathrm{KL}(p^\star\,\|\,\hat{p}) \;\lesssim\; \frac{1}{2}\int_{t_0}^T \mathbb{E}\big[\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2\big]\mathrm{d}s.\qquad{(2)}\] This new bound is strictly tighter than the standard one in ?? whenever the solenoidal component of errors is nonzero, motivating a better diagnostic of learning quality based on the geometry of diffusion. We propose a tractable surrogate to estimate the gradient component of the score error.

1.0.0.1 Contributions.

We make three contributions.

• We prove that no monotone function of the full \(L^2\) error can uniformly lower-bound any meaningful distributional divergence between the learned and target models (Theorem 1). The geometric reason is that only the gradient component enters the marginal Fokker-Planck dynamics, while the solenoidal component is structurally invisible.

• Building on this geometry, we derive a new upper bound on the Kullback-Leibler divergence between the learned and target models (Theorem 2, ?? ), tightening the standard Girsanov bound [9], [12], [13]. We further show this new bound is recoverable from a Girsanov argument applied to a marginally-equivalent representative of the learned reverse process, identifying the looseness of the standard bound as the cost of operating on path-space rather than marginal-space dynamics.

• We give a tractable estimator for the norm of the gradient component of the score error via a dual Sobolev identity (Section 4.2), and demonstrate empirically on Fashion-MNIST and CIFAR-10 that it correlates substantially better with sample quality (FID) than the full \(L^2\) error.

1.0.0.2 Other related work.

Recent studies have found that standard denoising score matching can be inefficient in various settings. When the data concentrate near a low-dimensional manifold, score matching can waste capacity by forcing the score network to fit large directions normal to the data manifold [27], [28]. In another vein, a popular line of work modifies the relative importance of score errors across noise levels or frequency regimes through timestep weighting and training reweighting schemes [29]–[32]. These studies reflect the empirical fact that not all score errors are equally relevant for sample quality. However, all these approaches treat the full score estimation error as the central quantity to analyze, without examining its geometric structure.

1.0.0.3 Notation.

Bold symbols (e.g.\(\mathbf{f}, \mathbf{X}_t\), \(\mathbf{Y}_t\)) denote \(\mathbb{R}^d\)-valued random variables or processes. \(\mathcal{P}(\mathbb{R}^d)\) (resp. \(\mathcal{P}_2(\mathbb{R}^d)\)) denotes the space of Borel probability measures on \(\mathbb{R}^d\) (resp. with finite second moment). For \(\mu, \nu \in \mathcal{P}_2(\mathbb{R}^d)\), the Kullback-Leibler divergence between these two measures is defined as \(\mathrm{KL}(\mu \,\|\,\nu)=\int_{\mathbb{R}^d} \log(\tfrac{\mu}{\nu})\mu\mathrm{d}\mathbf{x}\) if \(\mu\) is absolutely continuous with respect to \(\nu\) and \(+\infty\) otherwise.

For a distribution \(\mu \in \mathcal{P}_2(\mathbb{R}^d)\), the space of square integrable vector fields relative to \(\mu\) is defined as \(L^2(\mu, \mathbb{R}^d):=\{\mathbf{v}: \mathbb{R}^d \to \mathbb{R}^d~:~\mathbb{E}_{\mu}[\|\mathbf{v}\|^2_2]<\infty\}\). It is a Hilbert space when equipped with its \(L^2\) scalar product defined as follows. For \(\mathbf{v}, \mathbf{w}\in L^2(\mu, \mathbb{R}^d)\), we have \(\langle \mathbf{v}\,,\mathbf{w}\rangle_{L^2(\mu, \mathbb{R}^d)} = \mathbb{E}_\mu[\mathbf{v}\cdot \mathbf{w}]\), where \(\cdot\) denotes the standard Euclidean scalar product. This geometry yields the \(L^2\)-norm, defined as \(\|\mathbf{v}\|_{L^2(\mu)} = (\mathbb{E}_{\mu}[\|\mathbf{v}\|^2_2])^{1/2}\).

The space of continuous functions from \([0, T]\) to \(\mathbb{R}^d\) is denoted by \(C([0,T],\mathbb{R}^d)\). For a stochastic process \((\mathbf{Z}_t)_{t\in[0, T]} \in C([0,T],\mathbb{R}^d)\), \(\mathrm{Law}((\mathbf{Z}_t)_t)\) denotes its path-space distribution and \(\mathrm{Law}(\mathbf{Z}_t)\) denotes the law of its marginal at time \(t\). The space of differentiable scalar-valued functions on \(\mathbb{R}^d\) with continuous gradient is denoted by \(C^1(\mathbb{R}^d)\), and \(C_c^\infty(\mathbb{R}^d)\) denotes the space of smooth (i.e. infinitely differentiable) scalar-valued functions with compact support on \(\mathbb{R}^d\).

2 Background↩︎

2.0.0.1 Forward diffusion.

A forward diffusion on \([0, T]\) (\(T > 0\)) is defined by the stochastic differential equation (SDE): \[\label{eq:forward-sde} \mathrm{d}\mathbf{X}_t = \boldsymbol{f}_t(\mathbf{X}_t)\,\mathrm{d}t + \sigma_t\,\mathrm{d}\mathbf{B}_t, \qquad \mathbf{X}_0 \sim p^\star, \qquad t \colon 0 \to T,\tag{2}\] where \(\mathbf{B}_t\) is a standard \(d\)-dimensional Brownian motion, \(\boldsymbol{f}_t : \mathbb{R}^d \to \mathbb{R}^d\) is the drift, and \(\sigma_t > 0\) is the scalar diffusion coefficient. Under standard regularity conditions [33], this equation admits a unique solution, and its marginals \(p^\star_t := \mathrm{Law}(\mathbf{X}_t)\) satisfy the forward Fokker–Planck equation [20], [34], [35]: \[\label{eq:fokker-planck-forward} \partial_t p^\star_t = -\nabla\!\cdot\!\big(\boldsymbol{f}_t\,p^\star_t\big) + \frac{\sigma_t^2}{2}\,\Delta p^\star_t, \qquad p^\star_0 = p^\star.\tag{3}\] Because the data distribution \(p^\star\) may be non-smooth or even singular, the score \(\nabla \log p_t^\star\) can become unstable as \(t\downarrow 0\) [36], [37]. Therefore, as is common in the literature [9], [13], [38], we work on a truncated interval \([t_0, T]\) with \(t_0 >0\), where the diffusion has already regularized the law: under suitable nondegeneracy and regularity assumptions [39], [40], \(p^\star_t\) admits a smooth density for every \(t>0\), so that the score is well-defined and better behaved away from the singular endpoint.

2.0.0.2 Reverse-time diffusion.

Under regularity conditions on the drift and diffusion coefficients (assumed throughout), the forward process admits a time reversal [41]–[44]. The reverse-time SDE, running backwards from \(T\) down to \(t_0\), is written in the standard backward-time convention \(s:T \downarrow t_0\) (so that \(\mathrm{d}s < 0\)): \[\label{eq:reverse-sde-true} \mathrm{d}\mathbf{Y}_s = [\boldsymbol{f}_s(\mathbf{Y}_s) - \sigma_s^2\,\nabla_{\mathbf{x}} \log p^\star_s(\mathbf{Y}_s)]\,\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \mathbf{Y}_T \sim p^\star_T, \qquad s \colon T \downarrow t_0,\tag{4}\] where \(\bar{\mathbf{B}}_s\) is a Brownian motion under the reverse filtration. The marginals of this reverse SDE 4 exactly match those of the forward so that \(\mathrm{Law}(\mathbf{Y}_s)=p^\star_s\) for all \(s\), and they satisfy the backward Fokker-Planck equation (once again, backward in time i.e. \(\mathrm{d}s < 0\)): \[\label{eq:FP-true-reverse} \partial_s p^\star_s = -\nabla\!\cdot\!\Big(\big[\boldsymbol{f}_s - \sigma_s^2\,\nabla_\mathbf{x}\log p^\star_s\big]\,p^\star_s\Big) - \frac{\sigma_s^2}{2}\,\Delta p^\star_s, \quad \qquad s \colon T \downarrow t_0,\tag{5}\]

2.0.0.3 Score estimation.

In practice, the true score \(\nabla_\mathbf{x}\log p^\star_s\) is unknown, and is estimated via score matching [1], [2] or denoising score matching (DSM) [3]. In modern diffusion and score-based generative modeling, DSM is used to learn a time-dependent score field \(\mathbf{s}_\theta : \mathbb{R}^d \times [t_0, T] \to \mathbb{R}^d\) (parametrized by \(\theta\)) of noise-perturbed marginals, by regressing toward the conditional score [4], [7], [8], \[\label{eq:dsm95objective} \mathcal{L}_{\mathrm{DSM}}(\theta) = \mathbb{E}_{s \sim \mathcal{U}[t_0, T]}\!\left[\lambda(s)\mathbb{E}_{\mathbf{x}_0} \mathbb{E}_{\mathbf{x}_s | \mathbf{x}_0} \,\big\|\mathbf{s}_\theta(\mathbf{x}_s, s) - \nabla_{\mathbf{x}_s}\log p^\star_s(\mathbf{x}_s \mid \mathbf{x}_0)\big\|^2_2 \right],\tag{6}\] where \(\lambda:[0,T] \to \mathbb{R}_{>0}\) is a positive time weighting. At the population level and under mild regularity, the minimizer of 6 approximates the true score field of the noised marginals.

2.0.0.4 Learned reverse process.

Plugging the learned score into 4 defines the learned reverse process \[\label{eq:reverse-sde-learned} \mathrm{d}\hat{\mathbf{Y}}_s = [\boldsymbol{f}_s(\hat{\mathbf{Y}}_s) - \sigma_s^2\,\mathbf{s}_\theta(\hat{\mathbf{Y}}_s, s)]\,\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \hat{\mathbf{Y}}_T \sim p^\star_T, \qquad s \colon T \to t_0.\tag{7}\] Defining \(\mathbf{e}_s(\mathbf{x}) := \mathbf{s}_\theta(\mathbf{x}, s) - \nabla_\mathbf{x}\log p^\star_s(\mathbf{x})\) as the vector field of score estimation error, the marginals \(\hat{p}_s := \mathrm{Law}(\hat{\mathbf{Y}}_s)\) of the reverse process 7 satisfy the backward Fokker-Planck equation: \[\label{eq:FP-learned-reverse} \partial_s \hat{p}_s = -\nabla\!\cdot\!\Big(\big(\boldsymbol{f}_s - \sigma_s^2\,\nabla_\mathbf{x}\log p^\star_s\big)\,\hat{p}_s\Big) - \frac{\sigma_s^2}{2}\,\Delta \hat{p}_s + \sigma_s^2\,\nabla\!\cdot\!\big(\mathbf{e}_s\,\hat{p}_s\big), \quad \hat{p}_T = p^\star_T, \quad s \colon T \to t_0.\tag{8}\]

The learned generative distribution is \(\hat{p}_{t_0} := \mathrm{Law}(\hat{\mathbf{Y}}_{t_0})\). Moreover, writing \(\nabla_\mathbf{x}\log p^\star_s(\mathbf{x}) = \mathbf{s}_\theta(\mathbf{x}, s) -\mathbf{e}_s(\mathbf{x})\) in the true reverse Fokker–Planck 5 yields, \[\label{eq:FP-true-reverse-v2} \partial_s p^\star_s = -\nabla\!\cdot\!\Big(\big[\boldsymbol{f}_s - \sigma_s^2\,\mathbf{s}_\theta(\cdot, s)\big]\,p^\star_s\Big) - \frac{\sigma_s^2}{2}\,\Delta p^\star_s - \sigma_s^2 \nabla\cdot(p^\star_s \mathbf{e}_s). \quad \qquad s \colon T \to t_0,\tag{9}\]

2.0.0.5 Assumptions.

We make three assumptions, which are standard in the literature on diffusion-model convergence [12], [15], [36], [45] and required in the proofs of Theorems 1 and 2.

1. Marginal regularity. For all \(s \in [t_0, T]\), the marginals \(p^\star_s\) and \(\hat{p}_s\) are strictly positive, \(C^1\) probability densities on \(\mathbb{R}^d\) with finite second moment.

2. Score-error integrability. For all \(s \in [t_0, T]\), \(\mathbf{e}_s \in L^2(p^\star_s; \mathbb{R}^d)\cap L^2(\hat{p}_s; \mathbb{R}^d)\).

3. Score-error differentiability. For all \(s \in [t_0, T]\), the error field \(\mathbf{e}_s\) is differentiable in the sense that the divergence \(\nabla \cdot (q (\mathbf{x})\mathbf{e}_s(\mathbf{x}))\) is well-defined on all bounded sets, for \(q \in \{p^\star_s, \hat{p}_s\}\).

2.0.0.6 Integration by parts.

The assumptions above ensure that for \(s\in [t_0, T]\), the following integration-by-parts identity holds for all test functions \(\varphi \in C_c^\infty(\mathbb{R}^d)\) and \(q \in \{p^\star_s, \hat{p}_s\}\): \[\int_{\mathbb{R}^d} (\nabla \varphi(\mathbf{x})\cdot \mathbf{e}_s(\mathbf{x})) \, q(\mathbf{x}) \mathrm{d}\mathbf{x}= - \int_{\mathbb{R}^d} \varphi(\mathbf{x})\, \nabla \cdot (q (\mathbf{x})\mathbf{e}_s(\mathbf{x})) \mathrm{d}\mathbf{x}. \label{eq:int95by95parts}\tag{10}\]

2.0.0.7 Weighted divergence-free fields.

For a density \(\mu\), a field \(\mathbf{w}\in L^2(\mu;\mathbb{R}^d)\) is called divergence-free if \(\nabla\!\cdot(\mu\mathbf{w})=0\) in the weak sense, namely \(\int _{\mathbb{R}^d}\nabla\varphi(\mathbf{x})\cdot \mathbf{w}(\mathbf{x})\,\mu(\mathbf{x})\,\mathrm{d}\mathbf{x}=0\) for all \(\varphi\in C_c^\infty(\mathbb{R}^d)\).

3 Geometric Structure of Score Errors: Observable and Invisible Components↩︎

The true densities \((p^\star_s)_{s\in[t_0,T]}\) evolve according to the backward Fokker-Planck given by 9 , in which the score estimation errors \((\mathbf{e}_s)_{s\in[t_0,T]}\) enter exclusively through the scalar weighted divergence terms \((\nabla \cdot (p^\star_s \mathbf{e}_s))_{s\in[t_0,T]}\). This induces an invariance at the level of marginals: two errors with the same weighted divergences generate the same marginal dynamics. In particular, adding any divergence-free perturbation \(\mathbf{w}_s \in L^2(p^\star_s; \mathbb{R}^d)\), i.e. any \(\mathbf{w}_s\) such that \(\nabla \cdot (p^\star_s \mathbf{w}_s)=0\), leaves the backward Fokker-Planck unchanged. We therefore refer to these locally mass-preserving perturbations as invisible: they may contribute to the ambient \(L^2\) score error, but they do not affect the reverse marginal dynamics.

For clarity, in this section we fix a time \(s \in [t_0, T]\) in the reverse diffusion and study how quantities evolve at that time.

3.0.0.1 Helmholtz-Hodge Decomposition.

To characterize these invisible perturbations, we observe that for any differentiable vector field \(\mathbf{w}_s\), the following are equivalent (by weak integration by parts): \[\nabla\!\cdot\!\big(p^\star_s\mathbf{w}_s\big) = 0 \iff \text{ for all } \varphi\in C_c^\infty(\mathbb{R}^d), \text{ we have } \langle \nabla \varphi, \mathbf{w}_s\rangle_{L^2(p^\star_s;\mathbb{R}^d)}=0 \iff \mathbf{w}_s \in \mathcal{G}_s^\perp, \label{eq:divfree95cond}\tag{11}\] where \(\mathcal{G}_s:= \mathrm{cl}({\{\nabla\varphi:\;\varphi\in C_c^\infty(\mathbb{R}^d)\}})\) and \(\mathrm{cl}\) denotes the closure with respect to \(L^2(p^\star_s;\mathbb{R}^d)\); the inner product notation in 11 is defined on page . Thus, the set of vector fields \(\mathbf{w}_s\) such that \(\nabla\!\cdot\!\big(\mathbf{w}_sp^\star_s\big) = 0\) is the orthogonal complement of gradient fields in the space \(L^2(p^\star_s;\mathbb{R}^d)\). We can decompose \(L^2(p^\star_s;\mathbb{R}^d)\) into a gradient subspace \(\mathcal{G}_s\) and its orthogonal complement \(\mathcal{G}_s^\perp\). Expressing the error using this orthogonal decomposition yields the weighted Helmholtz–Hodge decomposition [21]: \[\label{eq:HH-decomposition-errors} \mathbf{e}_s=\Pi_{\mathcal{G}_s}\mathbf{e}_s+\Pi_{\mathcal{G}_s^\perp}\mathbf{e}_s \in L^2(p^\star_s;\mathbb{R}^d)=\mathcal{G}_s\oplus \mathcal{G}_s^\perp\tag{12}\] where \(\Pi_{\mathcal{G}_s}\) (resp. \(\Pi_{\mathcal{G}_s^\perp}\)) denotes the orthogonal projection operator on \(\mathcal{G}_s\) (resp. on \(\mathcal{G}_s^\perp\)). The second component is naturally divergence-free and thus structurally invisible at the reference density since \(\nabla \!\cdot\! (p^\star_s \mathbf{e}_s) = \nabla\!\cdot\!\big(p^\star_s\Pi_{\mathcal{G}_s}\mathbf{e}_s\big)\). Plugging 12 in the backward Fokker-Planck 9 implies, \[\label{eq:FP-HH-learned-reverse} \partial_s p^\star_s = -\nabla\!\cdot\!\Big(\big[\boldsymbol{f}_s - \sigma_s^2\,\mathbf{s}_\theta(\cdot, s)\big]\,p^\star_s\Big) - \frac{\sigma_s^2}{2}\,\Delta p^\star_s - \sigma_s^2 \nabla\!\cdot\!\big(p^\star_s\Pi_{\mathcal{G}_s}\mathbf{e}_s\big), \qquad s: T \downarrow t_0.\tag{13}\] Intuitively, the component \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\) corresponds to a transport of mass across space, which modifies the density. In contrast, \(\Pi_{\mathcal{G}_s^\perp}\mathbf{e}_s\) preserves mass locally (e.g. rotation) without affecting the density. As a result, two error fields with identical \(L^2\)-norm can induce radically different effects on the generated distribution, as illustrated in Figure 1.

The following result summarizes the discussion above and its consequences (proof in Appendix 6.1).

Theorem 1 implies that the ambient \(L^2\) score error is not the intrinsic quantity controlling sampling quality. Indeed the theorem shows that purely solenoidal errors may have arbitrarily large \(L^2\) norm while leaving all marginals unchanged, and that no distributional discrepancy can be uniformly lower-bounded by the full score error alone. Notably, this implies that the classical bound \[\label{eq:kl-upper-old} \mathrm{KL}(p^\star_{t_0}\;\|\;\hat{p}_{t_0}) \;\le\; \frac{1}{2}\int_{t_0}^T \sigma_s^2 \|\mathbf{e}_s\|_{L^2(p^\star_s)}^2\mathrm{d}s\tag{15}\] may be loose whenever score errors have large invisible (solenoidal) components.

By contrast, Theorem 1 identifies the gradient projection of the error, \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\), as the component entering the Fokker-Planck dynamics at time \(s\). Empirically, Figure 2 shows that \(\|\Pi_{\mathcal{G}}\mathbf{e}\|^2\) is a better proxy of sampling quality than the full error \(\|\mathbf{e}\|^2\) as it correlates better with the Fréchet Inception Distance (FID) [46], widely used to measure sampling quality of generative models for images.

3.0.0.2 Improved upper bound.

Based on Theorem 1, we replace the ambient \(L^2\) geometry by an observable geometry, to obtain the following sharper endpoint upper bound (proof in Appendix 6.2).

Theorem 2 identifies the component of the score error that is intrinsic to the learned marginal dynamics: since the learned Fokker–Planck equation depends on \(\mathbf{e}_s\) through \(\nabla\cdot(p^\star_s\mathbf{e}_s)\), the natural projection is the \(L^2(p^\star_s)\)-projection onto gradient fields.

3.0.0.3 Remarks

1. Link with Girsanov-based approach. In contrast with the Girsanov-based approach, this new bound is a direct consequence of the corrected Fokker-Planck dynamics 13 (discussed in Appendix 6.2). Appendix 7 reconciles the two by deriving the same bound via Girsanov’s theorem applied to a marginally-equivalent representative of the learned process whose drift retains only the observable component. This identifies the standard bound’s looseness as the cost of operating on path-space rather than marginal-space dynamics.

2. Conservative parametrizations. If the learned score vector field is parameterized as \(\mathbf{s}_\theta(\mathbf{x}, s) = \nabla_\mathbf{x}E_\theta(\mathbf{x},s)\) for a scalar potential \(E_\theta\) (e.g. as in energy-based models [48]), then \(\mathbf{e}_s = \nabla(E_\theta - \log p^\star_s)\) is a pure gradient field without any solenoidal component. In that case, the standard Girsanov bound coincides with 16 . For unconstrained architectures (e.g.U-Nets [6], [49] or Transformers [50], [51]), which dominate practice, the gap between \(\|\mathbf{e}_s\|_{L^2(p^\star_s)}\) and \(\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}\) can be substantial and our new bound is tighter and correlates more strongly with the Fréchet Inception Distance, as shown in Figure 2.

3. Gain with new upper bound. The new upper bound 16 is expected to significantly improve on the standard one 15 when the gradient components of the errors have small norms compared to the solenoidal components. Figure 1 illustrates this: by artificially creating solenoidal errors, our upper bound is not affected while that given in 15 explodes.

Figure 3 illustrates this effect during training: after an initial phase in which the observable gradient component is reduced, the remaining score error is dominated by the solenoidal component, creating a persistent gap between the classical full-error bound and the observable bound.

4 From Observable Bounds to a Computable Diagnostic↩︎

This section introduces a practical procedure to estimate the gradient component of score estimation errors. Theorem 2 bounds the divergence between the target and learned data distributions using the \(L^2(p^\star_s)\)-gradient projection \(\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}\). Computing expectations with respect to \(p^\star_s\) is practical since one has access to samples from the target distribution during training. However, explicitly projecting errors over the space of gradient fields is computationally hard. In this paragraph, we present an approach that avoids explicit projections by characterizing the gradient component of \(\mathbf{e}_s\) through an \(H^{-1}\) dual variational identity.

4.1 A Variational Reformulation↩︎

The projection \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\) is encoded in the Fokker-Planck equation 13 and can be characterized without explicitly projecting. Indeed, at each time \(s\), the dynamics of the true densities depend on \(\mathbf{e}_s\) only through \(\nabla\!\cdot(p^\star_s\,\mathbf{e}_s)\). Thus any vector field \(\mathbf{v}_s\) such that \(\nabla\!\cdot(p^\star_s\,\mathbf{e}_s) = \nabla \cdot (p^\star_s\,\mathbf{v}_s)\) induces the same effect on marginals; we refer to such fields as admissible. In Section 3, we observed that \[\label{eq:formulation-errors-gs} \nabla\!\cdot(p^\star_s\mathbf{e}_s) = \nabla\!\cdot(p^\star_s\Pi_{\mathcal{G}_s}\mathbf{e}_s),\tag{17}\] making the gradient-projection of errors one such admissible vector field, i.e. \(\Pi_{\mathcal{G}_s}\mathbf{e}_s \in \{\mathbf{v}_s: \nabla\!\cdot(p^\star_s\mathbf{e}_s) = \nabla\!\cdot(p^\star_s \mathbf{v}_s)\}\). Moreover, since \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\) is the only observable part of \(\mathbf{e}_s\), it contains the minimal and irreducible energy of \(\mathbf{e}_s\) that is visible through \(\nabla\!\cdot(p^\star_s\mathbf{e}_s)\).

In other words, among all vector fields \(\mathbf{v}_s\) that induce the weighted divergence \(\nabla\cdot(p^\star_s\mathbf{v}_s)=\nabla\cdot(p^\star_s\mathbf{e}_s)\), the field \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\) is the unique minimum-energy representative in \(L^2(p^\star_s;\mathbb{R}^d)\). Notably, the right-hand term of ?? does not require computing a projection, and just requires minimizing \(L^2(p^\star_s)\)-norms. However, this right-hand term is still not satisfactory since it requires solving a Poisson equation \(\nabla\cdot(p^\star_s\mathbf{v}_s)=\nabla\cdot(p^\star_s\mathbf{e}_s)\) which is known to be computationally intractable in high dimensions due to the curse of dimensionality inherent in standard numerical PDE solvers [53], [54]. Fortunately, the optimization problem in ?? exactly corresponds to an \(H^{-1}\) norm [55], which enjoys a tractable dual characterization. The \(H^{-1}(p^\star_s)\)-norm of \(\nabla\cdot(p^\star_s\mathbf{e}_s)\) is defined as \[\label{def:Hminus-1} \|\nabla\cdot(p^\star_s\mathbf{e}_s)\|_{H^{-1}(p^\star_s)}^2 := \inf_{\mathbf{v}_s:\, \nabla\cdot(p^\star_s\mathbf{e}_s)=\nabla\!\cdot(p^\star_s \mathbf{v}_s)} \int_{\mathbb{R}^d}\|\mathbf{v}_s(\mathbf{x})\|^2_2\,p^\star_s(\mathbf{x})\,\mathrm{d}\mathbf{x},\tag{18}\] and can be expressed in its dual form as [55], \[\label{eq:def-Hminus1-dual} \|\nabla\cdot(p^\star_s\mathbf{e}_s)\|_{H^{-1}(p^\star_s)}^2 = \sup_{\varphi\in C_c^\infty(\mathbb{R}^d)} \left\{ - 2\int_{\mathbb{R}^d}\varphi\,\nabla\cdot(p^\star_s\mathbf{e}_s)\,\mathrm{d}\mathbf{x} - \int_{\mathbb{R}^d}\|\nabla\varphi\|^2_2\,p^\star_s\,\mathrm{d}\mathbf{x} \right\}.\tag{19}\] From ?? 18 , we have \[\label{eq:final-approach} \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2 = \sup_{\varphi\in C_c^\infty(\mathbb{R}^d)} \left\{ - 2\int_{\mathbb{R}^d}\varphi\,\nabla\cdot(p^\star_s\mathbf{e}_s)\,\mathrm{d}\mathbf{x} - \int_{\mathbb{R}^d}\|\nabla\varphi\|^2_2\,p^\star_s\,\mathrm{d}\mathbf{x} \right\}.\tag{20}\]

These \(H^{-1}\)-norms have a natural geometric interpretation in optimal transport theory (discussed in Appendix 8). Equation 20 is an unconstrained optimization problem over the set of smooth and scalar-valued test functions, and directly yields our method, described below, to estimate \(\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2\).

4.2 Estimation of the Observable Error↩︎

The variational characterization of \(\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2_{L^2(p^\star_s)}\) in 20 suggests a clear training procedure to compute this norm: aim to solve the optimization by learning a scalar-valued critic network \(\varphi_\psi\) with parameters \(\psi\), via gradient ascent on the objective function in 20 . In the remainder of this section, we adopt the standard framework of variance-preserving forward diffusion [8] (i.e. we fix \(\boldsymbol{f}_s (\mathbf{x})= -\frac{1}{2} \beta(s)\mathbf{x}\) and \(\sigma_s = \sqrt{\beta(s)}\) in 2 for some schedule \(\beta:[t_0, T] \to \mathbb{R}_{>0}\)), widely used in practical training of diffusion models. Applying integration by parts 10 to the first integral in 20 , we obtain \[\label{eq:Hminus195force95dual} \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2 = \sup_{\varphi\in C_c^\infty(\mathbb{R}^d)} \left\{ 2\,\mathbb{E}_{\mathbf{x}\sim p^\star_s}\!\big[\nabla\varphi(\mathbf{x})\cdot \mathbf{e}_s(\mathbf{x})\big] - \mathbb{E}_{\mathbf{x}\sim p^\star_s}\!\big[\|\nabla\varphi(\mathbf{x})\|^2_2\big] \right\},\tag{21}\] which only requires expectations over \(p^\star_s\), which are accessible during training by sampling noised data \(\mathbf{X}_s\sim p^\star_s\). We introduce a critic potential \(\varphi_\psi(\mathbf{x},s)\), parameterized by a scalar-valued neural network with parameters \(\psi\), and define \[\label{eq:J95theta95psi95pop} \mathcal{J}(\psi) := \mathbb{E}_{s\sim \nu}\, \mathbb{E}_{\mathbf{x}\sim p^\star_s} \Big[ 2\,\nabla_\mathbf{x}\varphi_\psi(\mathbf{x},s)\cdot \mathbf{e}_s(\mathbf{x}) - \|\nabla_\mathbf{x}\varphi_\psi(\mathbf{x},s)\|^2_2 \Big],\tag{22}\] where \(\nu\) is a training distribution over times \([t_0, T]\). By 21 and the universal approximation properties of neural networks [56], [57], at the population level we have \(\sup_{\psi}\,\mathcal{J}(\psi) = \mathbb{E}_{s\sim \nu}\,\big[\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2\big]\).

To avoid the dependence of 22 on the unknown true score \(\nabla_\mathbf{x}\log p^\star\) through \(\mathbf{e}_s=\mathbf{s}_\theta(\cdot, s) -\nabla_\mathbf{x}\log p^\star_s\), we exploit the DSM identity [3]: for \(\alpha(s) = \exp(-\tfrac{1}{2}\int_{0}^s \beta(u)\mathrm{d}u)\) and \(\sigma(s)^2 = 1 - \alpha(s)^2\), one has that \(\mathbf{X}_s=\alpha(s)\mathbf{X}_0+\sigma(s)\boldsymbol{\varepsilon}\) with \(\mathbf{X}_0\sim p^\star\), and \(\boldsymbol{\varepsilon}\sim\mathcal{N}(0,\mathbf{I}_d)\), and \[\label{eq:dsm95target95conditional} \nabla_\mathbf{x}\log p^\star_s(\mathbf{x}) = \mathbb{E}\big[\mathbf{s}_{\mathrm{target}}(\mathbf{X}_s,s)\,\big|\,\mathbf{X}_s=\mathbf{x}\big], \quad \text{ where } \quad \mathbf{s}_{\mathrm{target}}(\mathbf{X}_s,s) = -\frac{\mathbf{X}_s-\alpha(s)\mathbf{X}_0}{\sigma(s)^2}.\tag{23}\] Consequently, the residual \(\mathbf{r}_\theta(\mathbf{x},s) := \mathbf{s}_\theta(\mathbf{x},s)-\mathbf{s}_{\mathrm{target}}(\mathbf{x},s)\) is an unbiased proxy for the score error in the sense that \(\mathbb{E}[\mathbf{r}_\theta(\mathbf{X}_s,s)\mid\mathbf{X}_s=\mathbf{x}]=\mathbf{e}_s(\mathbf{x})\). Plugging this proxy in 22 yields, \[\label{eq:J95theta95psi95emp} \widehat{\mathcal{J}}(\psi) := \mathbb{E}_{s\sim \nu}\, \mathbb{E}_{\mathbf{X}_0,\boldsymbol{\varepsilon}} \Big[ 2\,\nabla_\mathbf{x}\varphi_\psi(\mathbf{X}_s,s)\cdot \mathbf{r}_\theta(\mathbf{X}_s,s) - \|\nabla_\mathbf{x}\varphi_\psi(\mathbf{X}_s,s)\|^2 \Big], \quad \mathbf{X}_s=\alpha(s)\mathbf{X}_0+\sigma(s)\boldsymbol{\varepsilon}.\tag{24}\] We then train our critic network \(\varphi_\psi\) to solve the following maximization problem, \[\label{eq:forcing95loss95final} \hat{\mathcal{L}}_{\mathrm{obs}} := \sup_{\psi}\,\widehat{\mathcal{J}}(\psi).\tag{25}\] Intuitively, optimizing 25 makes the critic search for the steepest gradient direction for the frozen score network (i.e., the direction that best exposes the gradient component of the score error). Algorithm 4 summarizes the resulting estimator. The score network is kept fixed throughout the procedure.

At convergence and at the population level, the optimal critic value estimates \(\mathbb{E}_{s\sim\nu}\left[\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2 \right]\).

4.2.0.1 Uses of this diagnostic.

We use \(\widehat{\mathcal{L}}_{\mathrm{obs}}\) as a diagnostic rather than as a training objective, decoupling critic estimation from score-network optimization. The diagnostic can be applied in two ways.

During training: one may periodically refit the critic \(\varphi_\psi\) against the frozen current score network \(\hat{\mathbf{s}}_\theta\), so that \(\widehat{\mathcal{L}}_{\mathrm{obs}}\) tracks the gradient component of the score error along the optimization trajectory.

Post hoc: once fully trained, the score network can be frozen and the critic optimized more thoroughly to estimate the upper bound in Theorem 2. This is precisely the procedure that produces the gradient-component curves of Figures 2, 3, 5 and 6: at each checkpoint we freeze the score network, fit a fresh critic via 24 , and report \(\widehat{\mathcal{L}}_{\mathrm{obs}}\) as our estimate of the observable error (full setup, hyperparameters, and a critic-suboptimality ablation in Appendices 9.1–9.3). A single evaluation of \(\widehat{\mathcal{L}}_{\mathrm{obs}}\) takes a few minutes on one GPU and adds roughly \(5\)–\(10\%\) overhead to standard DSM training, with cost driven by critic optimization, since the estimator only requires forward-noised samples \(\mathbf{X}_s = \alpha(s)\mathbf{X}_0 + \sigma(s)\boldsymbol{\varepsilon}\). Unlike sample-based metrics such as FID [46], the diagnostic \(\widehat{\mathcal{L}}_{\mathrm{obs}}\) targets a quantitative upper bound on \(\mathrm{KL}(p^\star_{t_0}\|\hat{p}_{t_0})\) rather than perceptual quality and avoids reverse-SDE sampling entirely, making it complementary to FID and applicable in settings where sampling is expensive or where domain-specific feature extractors are unavailable.

5 Discussion and Future Work↩︎

We showed that the score estimation error decomposes orthogonally into a gradient component, which affects the marginal dynamics of the learned reverse process, and a solenoidal component that is structurally invisible to the Fokker-Planck equation, and therefore does not affect the marginals (Theorem 1). The consequence of this filtering effect of the Fokker-Planck dynamics is an impossibility result, making the full \(L^2\) error inadequate for any lower-bound on any distributional divergence between the learned and target distributions. This geometric perspective allowed us to derive an upper bound on the KL divergence between the learned and target data distributions (Theorem 2). The bound considers only the observable gradient components of the score error, making it strictly tighter than the standard Girsanov-based bound whenever the solenoidal component is non-zero on a set of positive time measure. Identifying the looseness of the Girsanov bound as the cost of operating at a path-space level to bound a divergence between marginals, we showed that our new bound can be recovered by working on a marginally-equivalent representative of the learned reverse process (Appendix 7). Our improved bound motivated a diagnostic that empirically correlates better with sample quality than the ambient \(L^2\) score error (Figure 2).

5.0.0.1 Limitations and future work.

The bound 16 integrates the instantaneous gradient component of the errors along the path, while the sampling distribution \(\hat{p}_{t_0}\) is shaped by how these gradient errors accumulate at the endpoint. We expect the diffusive part of the reverse SDE to have a smoothing effect, with some gradient errors injected at time \(s\) dissipated before reaching \(\hat{p}_{t_0}\), and others amplified. Identifying which gradient components along the path affect the endpoint marginal is a fundamental open question, and a sharper notion of endpoint observability building on our geometric framework would refine the bound further.

Beyond this refinement, our decomposition admits a natural interpretation in Wasserstein-2 geometry: \(\Pi_{\mathcal{G}_s}\mathbf{e}_s\) is the component of the score error that lies in the tangent space at \(p^\star_s\) of the Otto manifold; the solenoidal part lies in the orthogonal complement and does not move mass. In addition to providing geometric intuition, this perspective opens connections to gradient flows in Wasserstein space and to geometry-aware regularizers. Finally, our estimator of the norm of the gradient component of errors is currently used only as a training diagnostic. It freezes the score network to track gradient components during training. This could be extended into a training objective, either in combination with the standard DSM loss, or to make the critic guide the score network towards observable directions that translate into better learning and better capacity. We note that this raises subtle min-max optimization questions, since errors in estimating the critic could contaminate the score-network training signal if not carefully addressed.

Acknowledgements↩︎

NBK is supported by a G-Research Trinity College Studentship, and RET is supported by the EPSRC Probabilistic AI Hub (EP/Y028783/1). RV was supported in part by an EPSRC Mathematical Sciences Small Grant.

6 Proofs of Main Results↩︎

6.1 Proof of Theorem 1↩︎

6.2 Proof of Theorem 2↩︎

We first state a useful and standard lemma.

The lemma is similar to [12], but we give a proof for completeness.

Finally, here is the final proof of Theorem 1.

7 Connection with Girsanov’s Theory↩︎

In this section, we formally develop the link between our Helmholtz-Hodge approach and its implications on the standard Girsanov-based approach to bound \(\mathrm{KL}(p^\star_{t_0}\|\hat{p}_{t_0})\). We prove that the upper bound in Theorem 2 can also be established (with an additional assumption) using Girsanov’s Theorem [33], [59], [60] by applying the Helmholtz-Hodge decomposition to the stochastic integral of errors and its quadratic variation.

7.1 Technical Background↩︎

We work on a complete filtered probability space \((\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \in [t_0, T]}, \mathbb{P})\) equipped with a \(\mathbb{P}\)-Brownian motion \((\mathbf{B}_t)_{t \in [t_0, T]}\). We denote by \((\bar\mathbf{B}_t)_{t \in [t_0, T]}\) a Brownian motion under the reverse filtration and recall the true reverse SDE, \[\label{eq:reverse-sde-true-girsanov-part} \mathrm{d}\mathbf{Y}_s = [\boldsymbol{f}_s(\mathbf{Y}_s) - \sigma_s^2\,\nabla_{\mathbf{x}} \log p^\star_s(\mathbf{Y}_s)]\,\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \mathbf{Y}_T \sim p^\star_T,\tag{41}\] and the learned reverse SDE, \[\label{eq:reverse-sde-learned-girsanov-part} \mathrm{d}\hat{\mathbf{Y}}_s = [\boldsymbol{f}_s(\hat{\mathbf{Y}}_s) - \sigma_s^2\,\nabla_\mathbf{x}\log p^\star_s(\hat{\mathbf{Y}}_s) - \sigma_s^2\mathbf{e}_s(\hat{\mathbf{Y}}_s)]\,\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \hat{\mathbf{Y}}_T \sim p^\star_T,\tag{42}\] where \(\mathbf{e}_s(\mathbf{x}) := \mathbf{s}_\theta(\mathbf{x}, s) - \nabla_\mathbf{x}\log p^\star_s(\mathbf{x})\). Under mild (and assumed) regularity assumptions both SDEs are well-posed in law, and we denote \((\mathbf{Y})_s\) (resp. \((\hat{\mathbf{Y}}_s)_s\)) the weak solution [33], [61] of 41 (resp. of 42 ). In the original probability space \((\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \in [t_0, T]}, \mathbb{P})\), these are two distinct processes (functions of \(\Omega \to C([t_0, T], \mathbb{R}^d)\)) under a same Brownian motion \(\bar{\mathbf{B}}\).

Based on these two weak solutions, we define \(\mathbb{P}^\star=\mathrm{Law}((\mathbf{Y}_s)_s)\) and \(\hat{\mathbb{P}}=\mathrm{Law}((\hat{\mathbf{Y}}_s)_s)\) to be the path laws of 41 and 42 . The path-space distributions \(\mathbb{P}^\star\) and \(\hat{\mathbb{P}}\) live in the canonical space of paths \(C([t_0, T], \mathbb{R}^d)\), and are two different measures on this space. We equip this space with the coordinate process \((\mathbf{X}_s = \omega(s))_{s\in[t_0,T]}\), which is the only process we define on \(C([t_0, T], \mathbb{R}^d)\) and which corresponds to the value of the path at time \(s\). On this new path-space, there is one single process \((\mathbf{X}_s)_{s\in[t_0,T]}\) and two distributions such that, under \(\mathbb{P}^\star\), \((\mathbf{X}_s)_{s\in[t_0,T]}\) satisfies, \[\mathrm{d}\mathbf{X}_s = [\boldsymbol{f}_s(\mathbf{X}_s) - \sigma_s^2\,\nabla_{\mathbf{x}} \log p^\star_s(\mathbf{X}_s)]\,\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}^{\mathbb{P}^\star}_s, \qquad \mathbf{X}_T \sim p^\star_T, \label{eq:X95Pstar}\tag{43}\] and under \(\hat{\mathbb{P}}\), \((\mathbf{X}_s)_{s\in[t_0,T]}\) satisfies, \[\mathrm{d}\mathbf{X}_s = [\boldsymbol{f}_s(\mathbf{X}_s) - \sigma_s^2\,\nabla_\mathbf{x}\log p^\star_s(\mathbf{X}_s) - \sigma_s^2\mathbf{e}_s(\mathbf{X}_s)]\,\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}^{\hat{\mathbb{P}}}_s, \qquad \mathbf{X}_T \sim p^\star_T, \label{eq:X95Phat}\tag{44}\] where \(\bar{\mathbf{B}}^{\mathbb{P}^\star}_s\) (resp. \(\bar{\mathbf{B}}^{\hat{\mathbb{P}}}_s\)) is a Brownian motion under \(\mathbb{P}^\star\) (resp. under \(\hat{\mathbb{P}}\)). As discussed below, under suitable assumptions we have \(\mathbb{P}^* \ll \hat{\mathbb{P}}\), with the Radon-Nikodym derivative \(\frac{\mathrm{d}\mathbb{P}^\star}{\mathrm{d}\hat{\mathbb{P}}}\) given by Girsanov’s theorem (see 47 ). Furthermore, the processes \(\bar{\mathbf{B}}^{\mathbb{P}^\star}_s\) and \(\bar{\mathbf{B}}^{\hat{\mathbb{P}}}_s\) are related as \[\label{eq:innovation-shift} \bar{\mathbf{B}}^{\mathbb{P}^\star}_t = \bar{\mathbf{B}}^{\hat{\mathbb{P}}}_t - \int_{T}^{t} \sigma_s\,\mathbf{e}_s\,\mathrm{d}s.\tag{45}\]

7.1.0.1 Girsanov density.

Beyond Assumptions [ass:minimal-assumptions-marginals], [ass:minimal-assumptions-errors], and [ass:assumption-errors-diff], Girsanov-based arguments further require the following assumption:

1. Define the stochastic integral and its quadratic variation (both random variables), \[\label{eq:MZ95def} M_s := \int_{t_0}^s \sigma_u\, \mathbf{e}_u(\mathbf{X}_u)\cdot\mathrm{d}\bar\mathbf{B}^{\hat{\mathbb{P}}}_u, \qquad \langle M\rangle_s := \int_{t_0}^s \sigma_u^2\,\|\mathbf{e}_u(\mathbf{X}_u)\|_2^2\,\mathrm{d}u,\tag{46}\] and the associated Doléans–Dade exponential \[Z_s := \exp\!\left(M_s - \tfrac{1}{2}\langle M\rangle_s\right).\] We assume that \((Z_s)_{s\in[t_0, T]}\) is a true \(\hat{\mathbb{P}}\)-martingale on \([t_0, T]\).

Assumption [ass:A4-girsanov] ensures that the Girsanov transformation yields a valid probability measure (i.e., \(\mathbb{E}_{\hat{\mathbb{P}}}[Z_T] = 1\)). While \(Z_s\) is guaranteed to be a local martingale by construction, different sufficient conditions can be found in the literature to guarantee the full martingale status (Novikov’s condition [62], Beneš’s condition [63]).

Under Assumptions [ass:minimal-assumptions-marginals], [ass:minimal-assumptions-errors], [ass:assumption-errors-diff], and [ass:A4-girsanov], Girsanov’s theorem [33], [59], [60] implies that the Radon–Nikodym derivative of \(\mathbb{P}^\star\) with respect to \(\hat{\mathbb{P}}\) is given by the Doléans–Dade exponential [64] \(Z_T\), \[\frac{\mathrm{d}\mathbb{P}^\star}{\mathrm{d}\hat{\mathbb{P}}} = Z_T = \exp\Big(\int_{t_0}^T \sigma_s \mathbf{e}_s\cdot\mathrm{d}\bar\mathbf{B}^{\hat{\mathbb{P}}}_s - \tfrac12\int_{t_0}^T \sigma_s^2\|\mathbf{e}_s\|_2^2\,\mathrm{d}s\Big). \label{eq:Girsanov95ratio}\tag{47}\]

For a similar application of Girsanov’s theorem to analyze the sampling quality diffusion models, we refer the reader to [9].

7.2 Standard Girsanov-Based Upper Bound↩︎

To capture errors along the entire true 41 and learned 42 reverse trajectories, we define \[\begin{align} \label{eq:pathwise-score-error-energy} \varepsilon_{\star}^2 := \int_{t_0}^T \sigma_s^2 \| \mathbf{e}_s\|^2_{L^2(p^\star_s)}\,\mathrm ds , \qquad \hat{\varepsilon}^2 := \int_{t_0}^T \sigma_s^2 \| \mathbf{e}_s\|^2_{L^2(\hat{p}_s)}\,\mathrm ds, \end{align}\tag{48}\] where \(p^\star_s\) and \(\hat{p}_s\) are the time-\(s\) marginals of \(\mathbb{P}\) and \(\hat{\mathbb{P}}\), respectively. We note that, Assumption [ass:minimal-assumptions-errors] implies that, \[\varepsilon_\star^2 = \mathbb{E}_{\mathbb{P}^\star}\left[\int_{t_0}^T \sigma_s^2 \|\mathbf{e}_s(\mathbf{X}_s)\|_2^2 \mathrm{d}s \right] < \infty.\]

We now formally state the standard upper bound, which Theorem 2 improves upon.

We provide a short proof for completeness, although the result is standard [9], [12], [13], [38]. This proof also clearly indicates the critical points where modifications are needed to obtain our new bound (Theorem 2).

7.3 A Helmholtz–Hodge View of Girsanov’s Theorem↩︎

In this section we explain how the projected bound of Theorem 2 can also be recovered from a Girsanov-type argument. The key point is that Girsanov’s theorem is a path-space statement, whereas the Helmholtz–Hodge decomposition identified in Section 3 is intrinsic to marginal dynamics. Consequently, the projected bound is obtained not by applying Girsanov directly to the original learned process, but by first replacing the learned process by a marginally equivalent representative whose drift error is the observable component of the score error.

Recall that the true reverse process has path law \(\mathbb{P}^\star\) and satisfies \[\label{eq:app-true-reverse-girsanov} \mathrm{d}\mathbf{Y}_s = \big[ \boldsymbol{f}_s(\mathbf{Y}_s) - \sigma_s^2\nabla \log p^\star_s(\mathbf{Y}_s) \big]\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \mathbf{Y}_T\sim p^\star_T, \qquad s:T\to t_0 .\tag{51}\] The learned reverse process has path law \(\hat{\mathbb{P}}\) and satisfies \[\label{eq:app-learned-reverse-girsanov} \mathrm{d}\hat{\mathbf{Y}}_s = \big[ \boldsymbol{f}_s(\hat{\mathbf{Y}}_s) - \sigma_s^2\mathbf{s}_\theta(\hat{\mathbf{Y}}_s,s) \big]\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \hat{\mathbf{Y}}_T\sim p^\star_T , \qquad s:T\to t_0 .\tag{52}\] Moreover, plugging \(\nabla \log p^\star_s = \mathbf{s}_\theta(\cdot, s) - \mathbf{e}_s\) in 51 yields, \[\label{eq:app-true-reverse-girsanov-v2} \mathrm{d}\mathbf{Y}_s = \big[ \boldsymbol{f}_s(\mathbf{Y}_s) - \sigma_s^2\mathbf{s}_\theta(\mathbf{Y}_s, s) + \sigma_s^2 \mathbf{e}_s(\mathbf{Y}_s) \big]\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \mathbf{Y}_T\sim p^\star_T, \qquad s:T\to t_0 .\tag{53}\] The drift discrepancy between 52 and 51 is given by weighted score estimation errors \(\sigma_s^2\mathbf{e}_s\). Recall from the proof of Proposition 4 that

\[\label{eq:app-classical-girsanov-kl} \mathrm{KL}(p^\star_{t_0}\,\|\,\hat{p}_{t_0}) \le \mathrm{KL}(\mathbb{P}^\star\,\|\,\hat{\mathbb{P}}) = \frac{1}{2} \int_{t_0}^T \sigma_s^2 \|\mathbf{e}_s\|_{L^2(p^\star_s)}^2 \,\mathrm{d}s .\tag{54}\] The looseness of this bound comes from the fact that it considers all \(\|\mathbf{e}_s\|_{L^2(p^\star_s)}^2\) errors, even if not every component of \(\mathbf{e}_s\) actually affects marginals.

7.3.0.1 Bridge between Girsanov’s approach and ours.

To bridge this gap, we decompose errors between its observable gradient component and its invisible solenoidal component, as explained in Section 3. Formally, at each time \(s\), \[\mathbf{e}_s = \Pi_{\mathcal{G}_s}\mathbf{e}_s + \Pi_{\mathcal{G}_s^\perp}\mathbf{e}_s,\] where we recall that \(\mathcal{G}_s:= \mathrm{cl}\{\nabla \varphi:\varphi\in C_c^\infty(\mathbb{R}^d)\}\), with the closure taken with respect to \(L^2(p^\star_s;\mathbb{R}^d)\). We now construct a marginally equivalent representative in which these invisible directions are removed. As discussed in Section 3, this decomposition implies \(\nabla\cdot(p^\star_s\mathbf{e}_s) = \nabla\cdot\!\big(p^\star_s\Pi_{\mathcal{G}_s}\mathbf{e}_s\big)\), which allows us to rewrite the true marginal dynamics associated with the reverse SDE 53 as in 13 , which we restate for convenience:

\[\label{eq:app-observable-FP} \partial_s p^\star_s = -\nabla\cdot \Big( [ \boldsymbol{f}_s - \sigma_s^2 \mathbf{s}_\theta(\cdot, s) ]p^\star_s \Big) - \frac{\sigma_s^2}{2}\Delta p^\star_s - \sigma_s^2\nabla\cdot \big( p^\star_s\Pi_{\mathcal{G}_s}\mathbf{e}_s \big).\tag{55}\]

Now introduce a process \(\tilde{\mathbf{Y}} = (\tilde{\mathbf{Y}}_s)_{s\in[t_0, T]}\) that is the solution to the following SDE \[\label{eq:app-observable-representative} \mathrm{d}\tilde{\mathbf{Y}}_s = \big[ \boldsymbol{f}_s(\tilde{\mathbf{Y}}_s) - \sigma_s^2\mathbf{s}_\theta(\tilde{\mathbf{Y}}_s, s) + \sigma_s^2\Pi_{\mathcal{G}_s}\mathbf{e}_s(\tilde{\mathbf{Y}}_s) \big]\mathrm{d}s + \sigma_s\,\mathrm{d}\bar{\mathbf{B}}_s, \qquad \tilde{\mathbf{Y}}_T\sim p^\star_T.\tag{56}\] Then the marginals \((\tilde{p}_s)_{s \in [t_0, T]}\), where \(\tilde{p}_s:=\mathrm{Law}(\tilde{\mathbf{Y}}_s)\), solve the Fokker-Planck equation 55 . We denote by \(\tilde{\mathbb{P}}\) the path-space law of 56 .

On the other hand, since the true marginals \(p^\star_s\) solve the Fokker-Planck associated with the reverse SDE 53 , we have \[\label{eq:fp-true-learned-reverse-v2} \partial_s p^\star_s = -\nabla\cdot \Big( [ \boldsymbol{f}_s - \sigma_s^2\mathbf{s}_\theta(\cdot, s) ]p^\star_s \Big) - \frac{\sigma_s^2}{2}\Delta p^\star_s - \sigma_s^2\nabla\cdot (p^\star_s\mathbf{e}_s).\tag{57}\] But, since \(\nabla\cdot(p^\star_s\mathbf{e}_s) = \nabla\cdot\!\big(p^\star_s\Pi_{\mathcal{G}_s}\mathbf{e}_s\big)\) and initial conditions match, the marginal dynamics 57 and 55 exactly match, i.e., \(\tilde{p_s} = p^\star_s\) for all \(s \in [t_0, T]\).

In particular, endpoint marginals match \(\tilde{p}_{t_0}=p^\star_{t_0}\).

We may now apply Girsanov’s theorem between \(\hat{\mathbb{P}}\) and the observable representative \(\tilde{\mathbb{P}}\) exactly as is done in Proposition 4, replacing \(\mathbb{P}^\star\) by \(\tilde{\mathbb{P}}\). The SDEs 52 and 56 have the same diffusion coefficient and the same terminal law \(p^\star_T\); their drift discrepancy is \(\sigma_s^2\Pi_{\mathcal{G}_s}\mathbf{e}_s\). Assuming [ass:minimal-assumptions-errors] and [ass:A4-girsanov], the corresponding Doléans–Dade exponential (obtained from 46 by replacing \(\mathbf{e}_s\) with \(\Pi_{\mathcal{G}_s} \mathbf{e}_s\) ) is a true martingale, Girsanov’s formula gives \[\label{eq:app-projected-path-kl} \mathrm{KL}(\tilde{\mathbb{P}}\,\|\,\hat{\mathbb{P}}) = \frac{1}{2} \int_{t_0}^T \sigma_s^2 \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(\tilde{p}_s)}^2 \,\mathrm{d}s .\tag{58}\] Using 58 , \(\tilde{p}_{t_0}=p^\star_{t_0}\) and data-processing, this becomes, \[\label{eq:app-projected-path-kl-phat} \mathrm{KL}(p^\star_{t_0}\;\|\;\hat{p}_{t_0}) = \mathrm{KL}(\tilde{p}_{t_0}\;\|\;\hat{p}_{t_0}) \le \mathrm{KL}(\tilde{\mathbb{P}}\;\|\;\hat{\mathbb{P}}) = \frac{1}{2} \int_{t_0}^T \sigma_s^2 \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2 \,\mathrm{d}s,\tag{59}\] which is exactly the bound of Theorem 2.

8 Optimal Transport Interpretation↩︎

The variational characterization to compute the gradient component of errors, given in Proposition 3 by, \[\label{eq:grad-error-minimizes-L2-v2} \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2 = \inf_{\mathbf{v}_s:\, \nabla\cdot(p^\star_s\mathbf{v}_s)=\nabla\cdot(p^\star_s\mathbf{e}_s)} \int_{\mathbb{R}^d} \|\mathbf{v}_s(\mathbf{x})\|_2^2\,p^\star_s(\mathbf{x})\,\mathrm{d}\mathbf{x}.\tag{60}\] has an interpretation in the dynamic formulation of optimal transport. Recall \(\mathcal{P}_2(\mathbb{R}^d)\) denotes the space of probability densities on \(\mathbb{R}^d\) with finite second moment. Let \(\mu\in\mathcal{P}_2(\mathbb{R}^d)\) be a smooth positive density, and consider a differentiable curve of probability densities \((\mu_\tau)_{\tau\in(-\epsilon,\epsilon)}\) such that \(\mu_0=\mu\). Its tangent vector at \(\tau=0\) is the signed density \[\zeta := \left.\partial_\tau \mu_\tau\right|_{\tau=0}, \qquad \int_{\mathbb{R}^d}\zeta(\mathbf{x})\,\mathrm{d}\mathbf{x}=0.\] In the Benamou–Brenier geometry [65], this tangent perturbation is represented by a velocity field \(\mathbf{v}\in L^2(\mu;\mathbb{R}^d)\) through the continuity equation [66], [67] \[\zeta = -\nabla\cdot(\mu\mathbf{v}).\] This representation is not unique: if \(\mathbf{w}\) satisfies \(\nabla\cdot(\mu\mathbf{w})=0\), then \(\mathbf{v}+\mathbf{w}\) represents the same tangent vector \(\zeta\). Wasserstein geometry removes this ambiguity by selecting the minimum-kinetic-energy representative, \[\inf_{\mathbf{v}:\,-\nabla\cdot(\mu\mathbf{v})=\zeta} \int_{\mathbb{R}^d}\|\mathbf{v}(\mathbf{x})\|_2^2\,\mu(\mathbf{x})\,\mathrm{d}\mathbf{x}.\] Equivalently, this minimum-energy representative belongs to the closure of gradient fields [67]–[69], \[T_\mu\mathcal{P}_2 = \mathrm{cl}(\{\nabla\varphi:\varphi\in C_c^\infty(\mathbb{R}^d)\}),\] where the closure \(\mathrm{cl}\) is taken in \(L^2(\mu;\mathbb{R}^d)\). The Otto metric [70] equips this tangent space with the norm obtained by pulling back the \(L^2(\mu)\) geometry through the minimal kinetic energy gradient representative, \[\|\zeta\|_{T_\mu\mathcal{P}_2}^2 = \inf_{\mathbf{v}:\,-\nabla\cdot(\mu\mathbf{v})=\zeta}\int_{\mathbb{R}^d}\|\mathbf{v}\|_2^2\,\mu\,\mathrm{d}\mathbf{x},\] the second equality being exactly the minimum-energy selection described above. This is the \(H^{-1}(\mu)\)-norm [55] of the signed density \(\zeta\). To recover the setting of Section 4, fix a diffusion time \(s \in [t_0,T]\), take \(\mu=p^\star_s\), and let \(\zeta_s := -\nabla\cdot(p^\star_s\mathbf{e}_s)\) be the tangent vector generated by the score error (in the Fokker–Planck equation 57 the error contributes \(-\sigma_s^2\,\nabla\cdot(p^\star_s\mathbf{e}_s)\) to \(\partial_sp^\star_s\)). Then, by the definition above together with Proposition 3, \[\|\zeta_s\|_{T_{p^\star_s}\mathcal{P}_2}^2 =\inf_{\mathbf{v}:\,\nabla\cdot(p^\star_s\mathbf{v})=\nabla\cdot(p^\star_s\mathbf{e}_s)} \int_{\mathbb{R}^d}\|\mathbf{v}\|_2^2\,p^\star_s\,\mathrm{d}\mathbf{x} =\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|_{L^2(p^\star_s)}^2 .\] Equivalently, this tangent norm is the \(H^{-1}(p^\star_s)\) norm \(\|\nabla\cdot(p^\star_s\mathbf{e}_s)\|_{H^{-1}(p^\star_s)}\) of 18 : the Otto metric tensor and the \(H^{-1}(p^\star_s)\) norm coincide. Thus, the target-weighted observable error is the squared Wasserstein tangent norm of the infinitesimal curve of marginals generated by the score error by the velocity field given by the score error at \(p^\star_s\). The solenoidal component of \(\mathbf{e}_s\) is invisible because it changes the velocity representation without changing the tangent vector \(\zeta_s\), i.e. without changing the infinitesimal motion of probability mass.

9 Experimental Setting↩︎

9.1 Experimental Setup Fashion-MNIST↩︎

We detail the experimental setup used to produce the Fashion-MNIST error-decomposition curves, KL upper-bound estimates, and Spearman correlations of Section 3. All image experiments (on CIFAR-10 and Fashion-MNIST) were run on a single NVIDIA RTX 6000 Blackwell GPU with 96GB memory; a full CIFAR-10 score network training run required approximately 10 hours of wall-clock time.

9.1.0.1 Training setup.

We train standard score-based diffusion models on Fashion-MNIST [47] using a variance-preserving SDE with linear schedule \(\beta(t) = \beta_{\min} + t(\beta_{\max} - \beta_{\min})\), \(\beta_{\min} = 0.1\), \(\beta_{\max} = 20\), time horizon \(T = 4\) and truncation \(t_0 = 5\times 10^{-2}\). The score network \(\mathbf{s}_\theta\) is a three-resolution-level U-Net with self-attention [49], [50] parametrized by a base channel multiplier \(c\) and a per-resolution depth \(n_{\mathrm{res}}\). We sweep three capacity tiers, tiny (\(c{=}32\), \(n_{\mathrm{res}}{=}1\), \(0.7\)M parameters), small (\(c{=}48\), \(n_{\mathrm{res}}{=}2\), \(2.0\)M) and full (\(c{=}64\), \(n_{\mathrm{res}}{=}2\), \(3.6\)M), and train each for \(400\) epochs with AdamW [71] (learning rate \(5\times 10^{-4}\) with linear warmup over \(2{,}000\) steps and cosine decay, weight decay \(10^{-4}\), batch size \(512\), gradient clipping at \(1.0\)). We use a log-uniform time-sampling distribution over \([t_0, T]\), an exponential moving average of the score-network weights with decay \(0.9999\) used for all evaluations, and five independent seeds per capacity. Permanent checkpoints are saved every \(20\) epochs, yielding \(20\) checkpoints per run.

9.1.0.2 Estimating the error decomposition.

At each saved checkpoint we evaluate the gradient and solenoidal components of the score error on the EMA score network via the computational procedure of Section 4.2. The full \(L^2\) error is computed directly from the unbiased residual \(\mathbf{r}_\theta(\mathbf{X}_s,s) = \hat{\mathbf{s}}_\theta(\mathbf{X}_s,s) - \mathbf{s}_{\mathrm{target}}(\mathbf{X}_s,s)\) defined in equation 23 , by Monte Carlo averaging over a fixed pool of \(6{,}400\) noised samples drawn from the Fashion-MNIST training set with \(t\) sampled from the same log-uniform distribution used during training (the same \((\mathbf{X}_0, t, \boldsymbol{\varepsilon}, \mathbf{X}_t)\) tuples are reused across all checkpoints of a run to remove sampling noise from the curves). The gradient component is estimated via the dual variational identity of equation 21 : at each evaluation, we instantiate a freshly initialized critic potential \(\varphi_\psi\) and train it for \(1{,}500\) steps to maximize the empirical objective \(\widehat{\mathcal{J}}(\psi)\) of equation 24 , using Adam [72] (learning rate \(10^{-3}\) with cosine decay to \(10^{-5}\), gradient clipping at \(1.0\), batch size \(256\)). The critic architecture mirrors the score network’s encoder–decoder structure but with reduced width (base channel multiplier \(24\), scalar output, single attention block at the bottleneck), and its gradient \(\nabla_{\mathbf{x}} \varphi_\psi\) is taken via automatic differentiation. The solenoidal component is obtained as \(\|\Pi_{\mathcal{G}_s^\perp}\mathbf{e}_s\|^2 = \max(\|\mathbf{e}_s\|^2 - \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2, 0)\), with the floor at zero handling the rare cases where finite-sample noise makes the dual lower bound exceed the empirical full norm. For the final-checkpoint tight bound used in Theorem 2, we re-estimate the gradient component with a \(3\times\) wider critic (base channel multiplier \(48\)) trained for \(5{,}000\) steps with three independent restarts, reporting the across-restart mean and standard deviation.

9.1.0.3 KL upper-bound estimates.

The right panel of Figure 5 reports both Girsanov-style and Helmholtz–Hodge bounds on \(\mathrm{KL}(p^\star_{t_0}\|\hat{p}_{t_0})\). The two bounds, equations 15 and 16 , are respectively the time integrals \(\tfrac{1}{2}\int_{t_0}^T \sigma_s^2 \,\|\mathbf{e}_s\|^2_{L^2(p^\star_s)}\,\mathrm{d}s\) and \(\tfrac{1}{2}\int_{t_0}^T \sigma_s^2 \,\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2_{L^2(p^\star_s)}\,\mathrm{d}s\). As in the CIFAR-10 setup (described in the following section), we use a common multiplicative constant \(K = (T - t_0)\bar{\sigma}^2 / 2\) where \(\bar{\sigma}^2\) is the schedule-averaged \(\sigma_s^2\) over \([t_0, T]\). The constant \(K\) is identical for both curves and does not depend on the score network, epoch, or seed, so the ratio of the two curves — and therefore the inobservability gap — is faithful to theoretical predictions independently of \(K\).

9.1.0.4 Sample-quality metric.

Sample quality is measured by feature-FID computed in the \(128\)-dimensional penultimate-layer feature space of a small CNN classifier [73] (\(> 92\%\) test accuracy on Fashion-MNIST), evaluated on \(5{,}000\) samples generated via a \(200\)-step predictor–corrector reverse SDE with the EMA score network. Because the diagnostic and FID are computed from independent quantities (training residual vs.generated samples), no leakage is possible between the two columns of Figure 2.

9.1.0.5 Spearman rank correlation.

Per-seed Spearman correlations are computed across the \(20\) checkpoints of a run between feature-FID and each of \(\|\mathbf{e}\|^2\), \(\|\Pi_\mathcal{G}\mathbf{e}\|^2\). Pooled correlations are computed across the resulting \(60\) (seed, epoch) pairs at each capacity. We report Spearman rather than Pearson because both error norms and FID range over multiple orders of magnitude during training, so a rank-based statistic is more meaningful than a linear one. Across all nine (capacity, seed) combinations the gradient-component correlation lies in \([0.93, 0.99]\), whereas the full-error correlation ranges as low as \(0.65\), indicating that \(\|\Pi_\mathcal{G}\mathbf{e}\|^2\) is not only a better but also a substantially more reliable proxy for sample quality than the full \(L^2\) score error.

9.1.0.6 Critic-suboptimality control.

The dual estimator of \(\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2\) is, at finite critic capacity, a lower bound. To verify that the apparent solenoidal plateau is not an artifact of critic underfitting, we ran an ablation on the small capacity tier in which the eval-time critic was trained for \(4{,}500\) steps (\(3\times\) longer) and with a \(2\times\) wider architecture (base channel multiplier \(48\)). The resulting gradient-component curve was within the inter-quartile band of the reported one, and the solenoidal plateau was unchanged at the displayed log-scale resolution. We therefore interpret the plateau as a genuine property of the score network rather than as a critic-capacity artifact.

9.2 Additional Results Fashion-MNIST↩︎

The same qualitative behavior holds on Fashion-MNIST across three capacity tiers (tiny, small, full): the gradient component of the score error decreases by more than three orders of magnitude during training while the solenoidal component remains almost flat and dominates the total error throughout (Figure 6, Appendix 9.1); the standard Girsanov bound therefore stays effectively constant, and is a loose proxy of learning quality while our bound from Theorem 2 tracks sample quality and continues to decrease (Figure 5). This matches the Spearman correlation results between the FID and the two error components in Figure 2. We deliberately let the model size as well as the dataset vary (see results on CIFAR-10 in Figure 3) to illustrate that the asymmetry between observable and invisible components is not an artifact of architecture, dataset, or model capacity, but a structural property of standard score matching predicted by Theorem 1.

The same separation between observable and invisible error components, and the same gap between the standard and observable KL bounds, hold on Fashion-MNIST across three capacity tiers (Appendix 9.1, Figures 6–5).

9.3 Experimental Setup CIFAR-10↩︎

We detail the experimental setting used to produce the error decomposition and KL upper-bound curves of Figure 3.

9.3.0.1 Training setup.

We train a standard score-based diffusion model on CIFAR-10 with the variance-preserving SDE of [8]. With our notations in Section 2, we thus chose \(\boldsymbol{f}_t(\mathbf{x}) = -\tfrac{1}{2}\beta(t)\mathbf{x}\) and \(\sigma_t = \sqrt{\beta(t)}\), with a linear schedule \(\beta(t) = \beta_{\min} + t(\beta_{\max} - \beta_{\min})\), \(\beta_{\min} = 0.1\), \(\beta_{\max} = 20\), and time horizon \(T=1\) with truncation at \(t_0 = 0.01\). We parametrize the score network \(\mathbf{s}_\theta\) via a three-resolution level U-Net [49] with self-attention [50] (approximately \(12.6\)M trainable parameters). Models are trained for \(350\) epochs with AdamW [71] (learning rate \(2\times 10^{-4}\) with linear warm-up over \(3000\) steps and cosine decay to \(10^{-5}\), weight decay \(0.01\), batch size \(128\), gradient clipping at \(1.0\)). We use Min-SNR-\(\gamma\) loss weighting [29] with \(\gamma = 5\), and maintain an exponential moving average of the score-network weights with decay \(0.9999\) that is used for all evaluations. Curves are reported as median and inter-quartile range over five independent seeds.

9.3.0.2 Estimating the error decomposition.

Every \(10\) epochs, we evaluate the gradient and solenoidal components of the score error on the EMA score network via the computational procedure described in Section 4.2. The full \(L^2\) error is computed directly from the unbiased residual \(\mathbf{r}_\theta(\mathbf{X}_s,s) = \hat{\mathbf{s}}_\theta(\mathbf{X}_s,s) - \mathbf{s}_{\mathrm{target}}(\mathbf{X}_s,s)\) defined in equation 23 , by Monte Carlo averaging over a fixed pool of \(\sim 6{,}400\) noised samples drawn from the CIFAR-10 training set with \(t \sim \mathcal{U}[t_0, T]\) (the same \((\mathbf{X}_0, t, \boldsymbol{\varepsilon}, \mathbf{X}_t)\) tuples are reused across all evaluations to remove sampling noise from the curves). The gradient component is estimated via the dual variational identity of equation 21 : at each evaluation, we instantiate a freshly initialized critic potential \(\varphi_\psi\) and train it for \(4{,}000\) steps to maximize the empirical objective \(\widehat{\mathcal{J}}(\psi)\) in equation 24 , using Adam [72] (learning rate \(10^{-3}\) with cosine decay to \(10^{-5}\), gradient clipping at \(1.0\)). The critic architecture mirrors the score network’s encoder-decoder structure but with reduced width (base channel multiplier \(96\), scalar output), and its gradient \(\nabla_\mathbf{x}\varphi_\psi\) is taken via automatic differentiation. To stabilize the estimator we average over three independent critic restarts per evaluation; the standard deviation of the gradient estimate across restarts is reported alongside the mean. The solenoidal component is then obtained as \(\|\Pi_{\mathcal{G}_s^\perp}\mathbf{e}_s\|^2 = \max(\|\mathbf{e}_s\|^2 - \|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2, 0)\), with the floor at zero handling the rare cases where finite-sample noise makes the dual lower bound exceed the empirical full norm.

9.3.0.3 KL upper-bound estimates.

The right panel of Figure 3 reports both Girsanov-style and Helmholtz–Hodge bounds on \(\mathrm{KL}(p^\star_{t_0}\|\hat{p}_{t_0})\). The two bounds, equations 15 and 16 , are respectively the time integrals \(\tfrac{1}{2}\int_{t_0}^T \sigma_s^2 \,\|\mathbf{e}_s\|^2_{L^2(p^\star_s)}\,\mathrm{d}s\) and \(\tfrac{1}{2}\int_{t_0}^T \sigma_s^2 \,\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2_{L^2(p^\star_s)}\,\mathrm{d}s\). For simplicity, we chose the same constant for these two bounds, corresponding to \(K = (T - t_0) \bar{\sigma}^2/2\) where \(\bar{\sigma}\) is the schedule-averaged \(\sigma_s^2\) over \([t_0, T]\). The constant \(K\) is identical for both curves and does not depend on the score network, the epoch, or the seed, making the ratio of the two curves at any epoch (and therefore the inobservability gap) faithful to theoretical predictions.

9.3.0.4 Critic-suboptimality control.

The dual estimator of \(\|\Pi_{\mathcal{G}_s}\mathbf{e}_s\|^2\) is, at finite critic capacity, a lower bound. To verify that the apparent solenoidal plateau is not an artifact of critic underfitting, we ran an ablation in which the eval-time critic was trained for \(12{,}000\) steps (3\(\times\) longer) and with a wider architecture (base channel multiplier \(192\) instead of \(96\)). The resulting gradient-component curve was within the interquartile band of the reported one, and the solenoidal plateau was unchanged at the displayed log-scale resolution. We therefore interpret the plateau as a genuine property of the score network rather than as a critic-capacity artifact.

9.4 Additional Experiments on CIFAR-10↩︎

References↩︎