Implicit Regularisation in Diffusion Models: An Algorithm-Dependent Generalisation Analysis
Authors: Tyler Farghly
Patrick Rebeschini
George Deligiannidis
Arnaud Doucet
Venue: ICLR 2026 Poster (2026)
Keywords: diffusion models
generalization
stability
algorithm dependent
learning theory
Primary area: learning theory
TL;DR: We propose a framework for analysing the generalisation properties of diffusion models. We use it to identify and analyse three key sources of implicit regularisation: early stopping, coarse sampler discretisation, and optimisation with SGD.

Abstract: The success of denoising diffusion models raises important questions regarding their generalisation behaviour, particularly in high-dimensional settings. Notably, it has been shown that when training and sampling are performed perfectly, these models memorise training data—implying that some form of regularisation is essential for generalisation. Existing theoretical analyses primarily rely on algorithm-independent techniques such as uniform convergence, heavily utilising model structure to obtain generalisation bounds. In this work, we instead leverage the algorithmic aspects that promote generalisation in diffusion models, developing a general theory of algorithm-dependent generalisation for this setting. Borrowing from the framework of algorithmic stability, we introduce the notion of score stability, which quantifies the sensitivity of score-matching algorithms to dataset perturbations. We derive generalisation bounds in terms of score stability, and apply our framework to several fundamental learning settings, identifying sources of regularisation. In particular, we consider denoising score matching with early stopping (denoising regularisation), sampler-wide coarse discretisation (sampler regularisation), and optimising with SGD (optimisation regularisation). By grounding our analysis in algorithmic properties rather than model structure, we identify multiple sources of implicit regularisation unique to diffusion models that have so far been overlooked in the literature.

Published as a conference paper at ICLR 2026
IMPLICIT REGULARISATION IN DIFFUSION MODELS: AN ALGORITHM-DEPENDENT GENERALISATION ANALYSIS
Tyler Farghly, Patrick Rebeschini, George Deligiannidis Arnaud Doucet Department of Statistics
University of Oxford {farghly, rebeschini, del deligian, ox
ABSTRACT
The success of denoising diffusion models raises important questions regarding their generalisation behaviour, particularly in high-dimensional settings. Notably, it has been shown that when training and sampling are performed perfectly, these models memorise training data- aata-implying that some form of regularisation es- sential for generalisation. Existing theoretical analyses primarily rely algorithm- independent techniques such as uniform convergence, heavily utilising model structure to obtain generalisation bounds. In this work, instead leverage the algorithmic aspects that promote generalisation in diffusion models, developing a general theory of sggorthmm-depender generalisation for this setting. Borrowing from the framework of algorithmic stability, we introduce the notion of score stability, which quantifies the sensitivity of core-mmaching algorithms dataset perturbations. We derive generalisation bounds in terms of score stability, and apply our framework to several fundamental learning settings, identifying sources of regularisation. In particular, we consider denoising score matching with early stopping (denoising regularisation), sampler-wide coarse discretisation (sampler regularisation) and optimising with SGD (optimisation regularisation). By ground- ing our analysis algorithmic properties rather than model structure, we identify multiple sources of implicit regularisation unique to diffusion models that have far been overlooked the literature.
INTRODUCTION
Diffusion models DSoht-Dickstein et al. 2015 Ho et al. 2020 Song et al. 2021 are class of generative models that have achieved state-of the--e-rrtttttt performance across image, audio, video, and protein synthesis tasks (Rombach et al. 2022: Saharia et al. 2022 Ramesh et al. 2022 Watson al. 2023 Esser al., 2024). Their ability to generate high-quality samples from complex, high- dimensional distributions with limited data motivates the need for a theoretical understanding of the mechanisms underpinning their strong generalisation capabilities.
The goal of diffusion models is to generate new synthetic samples from data distribution Vdata using /
a finite set of data points { j Central to the methodology a unique approach to generating data, formulating it as the iterative transformation of noise into data, or equivalently, the reversal of diffusion process (Song et al.) 202T). This diffusion process, called the forward process, is defined by the stochastic differential equation (SDE),
Xo Vdatta for 2 where W. denotes the Brownian motion in Rd and T >> 0 is the terminal time. It can then be shown that the time-reversal of this process, Y, == Xx admits a weak formulation as a solution to the SDE,
d.X, = -ax; dt + V2W,,
t E [0.T].
(1)
dY, aY,dt 2V
+ /2dW,, Yo PT.
E (O,
(2)
where P; denotes the marginal density of X, (Haussmann & Pardoux 1986). Therefore, simulating samples from Vittt Pn can be achieved by solving the diffusion process in (2). which requires

---

Published as a conference paper at ICLR 2026
Closest training data
C- 0.000 1- 3,437 5,302 1-5.997
L-8.702 t- 9.3774
t- t-9.743
1-9.9 1-9.998 <- 9-9999
X
X
X
X
X
Figure I: Samples generated using the empirical score function on CIFAR-10 compared to the closest image in the dataset, illustrating memorisation of the training data.
an approximation of the score function, V log Pt. This achieved by fitting a time-dependent deep neural network to minimise weighted -d-distance called the (population) score matching loss:
fam (s T) Ex, (xx(l)((X,,t))
lo (X.)) (dt),
where r is probability measure over (0, T| that determines the weighting of the timepoints. Since V log P. is unknown, typically the (population) denoising score matching loss fasm. which differs from &m (s) only by constant, is used instead and is then approximated using the dataset, forming the empirical denoising score matching loss ldsm (see equations (6) and The score network s(x, t) is trained on this objective using standard stochastic optimisation methods relying mini- batching. Once an approximation is obtained, samples are generated by numerically solving the reverse-time SDE, Both score matching and sampling introduce distinct challenges and design choices that impact the quality of model output (Karras et al. 2022).
Score matching presents key difference from standard supervised learning. In the space of all L score functions, the empirical objective ldsm possesses a unique minimiser empirical score function-as a result of the integration over X, Xo (see Lemma] D}. This contrasts with traditional supervised learning, where the empirical risk minimisation problem can have infinitely many solutions (e.g., in overparameterised regression) and often requires regularisation to be well-posed. As shown in Figure [ sampling with this empirical score leads to exact recovery of the training data (Pidstrigach) 2022). This behaviour is distinct from 'benign overfitting', a phenomenon from the deep learning literature where interpolating the data does not necessarily prevent generalisation (Bartlett al. 2021 Zhang et al. 2021). This divergence suggests that existing theory may be insufficient to explain the success of diffusion models, highlighting the need for new frameworks tailored to this setting.
Recently, there has been a drive towards developing theory for better understanding the unique structure of diffusion models. The most developed subset of this work focuses on connecting sample quality score matching by deriving upper bounds distribution error (c.g. KL divergence, total variation, or Wasserstein distance) between model samples and the data distribution, controlling it by the population score matching loss (De Bortoli et al.
2021 De Bortoli 2022; Lee al. 2022
Chen et al. 2023 Benton et al.) 2024 Potaptchik et al. convergence bounds, typically take the form,
2024
These results, often referred to as
Distribution error S fm +
where A is the discretisation eπrоr of the sampling scheme, which can be made small with sufficiently fine discretisation. However, since am is not computable, these bounds say little about performance under empirical guarantees- theat is, their generalisation properties. One line of work, initiated by
Oko et al. (2023) and extended (Azangulov et al. 2024 TTang & Yang 2024), applies classical uniform convergence theory to bound the generalisation gap from the decomposition.
lam (s) E.m (s) + lam (5) la
generalisation gap
where &am denotes the empirical counterpart to &m These results rely on covering number bounds for specific classes of neural networks and, while informative, they are limited to carefully chosen
2

---

Published as a conference paper at ICLR 2026
model classes and do not account for algorithmic properties. An alternative approach by De Bortoli (2022) uses a decomposition of the Wasserstein distance that leverages convergence properties of the empirical measure. Though more model-agnostic, this method overlooks how diffusion models generate novel data. Both lines of work are fundamentally algorihhn-independent in that they lack any utilisation of the algorithmic aspects that uniquely define diffusion models. Recent efforts aim to incorporate algorithmic effects by restricting the problem. For instance, Shah et al. (2023); Chen et al. (2024) consider Gaussian mixture targets, while Li et al. (2023); Yang (2022) study random feature models. These settings allow for finer analysis of the role of the score matching algorithm, but remain limited in scope, leaving open the challenge of developing a more general eegorithm-dependee theory of generalisation in diffusion models.
As noted earlier, if the empirical score matching loss was completely minimised and sampling was performed perfectly, the diffusion model would simply return training data, failing to generalise. Therefore, the observed success of diffusion models producing novel data implies that, in practice, they either avoid completely minimising &mm or must avoid perfectly sampling. This suggests that (implicit) regularisation in the score matching sampling algorithm is crucial for generalisation, making algorithmic considerations essential for understanding diffusion models.
1.1 OUR CONTRIBUTIONS
We introduce score stability, general, algorithn-dependen framework for analysing diffusion model generalisation based on the classical approach of algorithmic stability. This framework quantifies an algorithm'ss 's dependence on individual training examples, from which we derive expected generalisation gap bounds for score matching losses. Using the score stability framework, we then analyse several examples of score matching algorithms, identifying three distinct sources of implicit regularisation in diffusion model training and sampling: noising, sampler, and -idiiuuatioo-nndueedd regularisation.
Denoising regularisation. To begin with, we consider the empirical risk minimisation algorithm (ERM) that minimises &da over a hypothesis class H. Through score stability analysis, we reveal regularisation source within this objective when early stopping of the forward process is used standard practice in the diffusion model literature. Utilising properties of the noising forward process, we obtain generalisation gap bounds with near-linear rate, -d" (A(e- N2 + minz ( )c/2 for any 1, where > 0 is the early stopping time and d' is the dimension of the data support.
Sampler regularisation. We then apply this analysis to discrete-time sampling algorithms, deriving statistical guarantees for the expected KL divergence between the true data distribution and samples generated by the diffusion model. The bound we derive is formed of two stages: we obtain generic rates e-1/2 N ming f )o/d but when N and miny l.mm are sufficiently small relative to c, we obtain bounds with rates c c-d"/Accd" N2 + ming ) that are faster in N and minH l To derive this bound, discretisation. We find that by increasing discretisation coarseness, we can improve the generalisation gap bound at the expense of worsening the discretisation error term.
utilise regularisation brought about by the coarseness of the
Optimisation regularisation. Finally, we consider the role of the optimisation scheme, analysing stochastic gradient descent (SGD) with gradient clipping and weight decay. On the model class, we assume only structural assumptions typical in the optimisation literature, including non-global Lipschitz and smoothness assumptions. While this initially yields bounds that grow with the number of iterations, we more closely inspect the impact of the nigh-variance gradient estimator used in diffusion training. We show this gradient noise induces a contractive behaviour in the training dynamics, which we harness to obtain stability bounds that do not grow with the number of iterations (Proposition [4), showing that the noisy dynamics enable tighter generalisation guarantees.
BACKGROUND
Suppose that the data distribution Vaata is on R2 and we are provided a finite dataset of samples S= {x1,.... ..., Σn' which assumee are sampled independently and identically (i.id.) from Vaata- As discussed in the introduction, diffusion models are formed two distinct stages. The first stage, score
matching, consists of learning approximation to the score function V log Pt using the dataset S. In
3

---

Published as a conference paper at ICLR 2026
this work, we take score function be any function belonging the set (R x 10, 1); """, the set of Borel measurable functions of the form I x [0.T] is taken to be any mapping of the form Agm (UN= (R*) yoM) x Ω + H. where H. is a measurable subset of L (IR x [0, T]: 1R)) Here, Ω is the event space belonging to probability space (S2, F,,
score
R. Then, a score matching algorithm
The second stage diffusion models, sampling, consists generating samples with the learned score function. We take sampling algorithm to be mapping of the form, Asamp H + P( where ) denotes the set of Borel measures on F. Typically, sampling is performed using approximation to the reverse process given in (2) replacing V log Pt with the learned score function s(;,t) and replacing its initial distribution, pr with data-independent prior, Pprior NO, o In case of a > O, we choose o2 prioo = ax so that the prior coincides with the stationary distribution of the forward process, and when a O we simply set of = 2T. With this, we arrive at the SDE,
dÝ, a 2s(Y,, t)dt + VadW Yo Pptior-
(5)
Thus, sample generated by sampling from Yr, or more commonly, the process is terminated early, sampling from Yr_, for some small e 0. Therefore, diffusion models are density estimation algorithms formed from the composition Asamp Asm-
Denoising score matching and overfitting, As stated in the introduction, computing Sun requires access to the population score function, V log So instead, the (population) denoising score matching loss is used in its place:
dam (s; T) E Exo- If [xxxxoll(((xx,1 Pl ((t))
(6)
which differs from lam (8) only by constant Cam, (see Lemma 16 whilst being easier to approximate without access to V log (Hyvärinen, 2005). Since Piju is Gaussian kernel, its score is given by,
v log log)))())lz)
24
He =
of
a (1-1(1 --)) if a otherwise.
(7)
0f
2t,
In practice, the objective in 66 is further approximated via Monte Carlo estimation using the dataset which leads to the empirical denoising score matching loss,
lasam (s; S,T)
Σ [x.xxxll(((xe.11 N
o X, T(dt).
In the following lemma, we highlight the important property that this can equivalently be defined as the denoising score matching objective for the process X, which evolves as in [ but with the initial distribution given by the empirical distribution, Xo ~ 1 E 5n, (dx).
Lemma 1. The objective la (s; is identical, up to constant, to the objective
/ [x)]=(8,,)) - (ogg((()))111)))((()
la (s;
(9)
where De is the marginal density of Xc. Therefore, any minimiser of S,T)):; L x [0, T); IR)) is identical V Vlg a.e. for any E supp(T).
See Appendix A.2 for the proof. This lemma shows that, unlike in traditional supervised learning problems, the empirical objective here admits single unique minimiser, the empirical score function, V log f.. The nature of this score function and the samples it generates has been the focus several recent studies, notably (Pidstrigach 2022) which shows that with perfect sampling, any score function sufficiently close to log P. recovers the training data.
Other notation. When the score matching algorithm Agm is random, use Asm (S) shorthand for the random score function (x, t. + Agm (S. (x,, t). Given two random score functions s, s', we let T(s, s') denote the of all couplings of these random functions (Appendix A. for details).
4

---

Published as a conference paper at ICLR 2026
3 SCORE STABILITY AND GENERALISATION
Algorithmic stability is classical technique in learning theory used to understand the generalisation properties of variety of important learning algorithms (Kearns & Ron 1999| Devroye & Wagner 1979; Bousquet & Elisseeff 2002) Hardt et al.] 2016). While there are various formulations, they all share the common aim of connecting properties of learning algorithm to its robustness under changes in the dataset. Its use has primarily been focused around regression and classification problems- this section, we propose a notion of stability that applies specifically diffusion models.
We introduce the notion of score stability which quantifies how sensitive score matching algo- rithm Asm individual changes in the dataset. We do this by defining the adjacent dataset S == ..., (ii Ii--, H. 1i+1, .... EN' where i ~ Vdatas independent from S, and then measuring the similarity between the score functions $ (5) and s (S')
Definition 2. A score matching algorithm Aam is score stable with constant Estab > O for any i [N] it holds that,
I E||s(X,,)) (' 2. (B,
Es,2
inf
I Esabbr
where r; = [(A (S), As
Since Agm may be random, we define score stability terms of the best-case coupling of the random score functions 8, 8 We recall that [(::)) denotes the set of couplings between two random score functions, and when it is not random, is given by the singleton I; = {(Asm Asm (S²))) In the following theorem, generalisation gap by the score stability constant.
connect score stability generalisation by controlling the expected
Theorem 3. Suppose that the score matching algorithm Asm score stable with constant eaab Then, with 8 = Aum (S), it holds that
[[Caam [G:T]] 1/2 E[[asam S,T) S Estab.
(10)
Furthermore, it holds that
E Cl (8; 11---1165 (3; < 2 smub IE Edsm (5; {,,)) + =
(11)
With Theorem 3 we obtain that the generalisation gap for both the denoising score matching loss and the score matching loss decays at the same rate as score stability. We can further simplify the bound for the score matching loss using the fact that lassm and l are identical up to a constant, to obtain,
E[lam (8: [B] S IE [eam (B] + Estah C1/2 Estab
One should expect that if the score matching algorithm is effective, both å and 3 converge to the ground truth as N grows, and thus stab should decrease to O. Ascertaining the rate at which N decreases requires stability as an algorithm-dependent approach. This contrasts with uniform learning, which utilises control over the hypothesis class, providing worst-case bound that is independent from the algorithm.
analysis of the algorithm at hand, hence the categorisation of algorithmic
In the following sections, we apply the framework of score stability to some common learning settings for diffusion models. We derive estimates of the score stability constant for these algorithms and identify features that promote generalisation.
EMPIRICAL SCORE MATCHING AND IMPLICIT REGULARISATION
We begin our examples by considering the score matching algorithm that minimises the empirical denoising score matching loss. Given hypothesis class H. C (RR x [O,T]; Rd), we define this algorithm by.
Aemm (S) = argmin (s; S,
While this algorithm is not often used in practice, it is the natural analogue to empirical risk minimi- sation from traditional supervised learning and thus serves as canonical example. We consider the setting of the manifold hypothesis where the data distribution is supported submanifold of IRd
5

---

Published as a conference paper at ICLR 2026
Assumption 4. Suppose that Vdata is supported smooth submanifold IR that has dimension d' and reach Treach > O. Furthermore, its density on the submanifold, Pv. satisfies cy := inf > 0.
The reach describes the maximum distance where the projection to the manifold is uniquely defined and therefore, it quantifies the maximum curvature of the manifold. We refer to Appendix [A.3 for the full definition. Several recent works have considered the assumption that Vdata lies on a submanifold of RRd. These works argue that d can often be far smaller than d and so dependence with respect to d" over d is favourable (De Bortoli 2022 Pidstrigach, 2022: Loaiza- Ganem et al. 2024 Potaptehik (et al.) 20224; Huang et al. 2024). The assumption that the density is bounded from below has also appeared in several of these works (Potaptchik et al. 2024, Huang et al. 2024). We also make the following assumption about the class of score networks.
Assumption 5. Suppose there exists Dz 2 0 such that for any s, E H. it holds that
[s( - s' (.,,t)ll== S Du/o?, for all E supp(r).
Under these assumptions, we obtain the following estimate for the stability constant.
Proposition 6. Suppose that assumptions 5 hold and that e := inf supp(r) E (0, (rach) then for any E (0, 1) and sufficiently large N. the score matching algorithm Aem is score stable with.
Stab 2 C(CCsm 2 E(esmn ($))) s
C =
D V
An interesting feature of Proposition 6 is that we obtain generalisation bounds for the empirical risk minimisation algorithm under relatively light assumptions on the structure of the hypothesis class without any additional regularisation. This contrasts with algorithmic stability in the setting of traditional supervised learning, where ERM is stable only when restricting the hypothesis class or with the use of explicit regularisation (Zhang et al., 202T; Bousquet & Elisseeff, 2002). Here, we show that the denoising score matching loss possesses the unique property that it is stable without the need for additional regularisation, suggesting that the denoising score matching loss possesses a form of implicit regularisation.
Cam = O(d' c Since )) the bound only depends on d' and not d, this suggests that diffusion models are automatically The bound also heavily depends on with it being smaller for larger e and growing exponentially fast as approaches zero, indicating that the natural regularisation present in the score matching objective is more prevalent at larger noise scales. The requirement to have e >> O is closely related to the technique of early stopping which is frequently used in the diffusion model literature (Song & Kingma 2021 Karras et al. 2022). This is where the backwards process Y, is terminated early by some small amount of time to avoid irregularity issues of the score function when close to convergence. Other theoretical works have also identified the importance of early stopping in the generalisation properties of diffusion models (Oko et al., 2023 Azangulov et al.) 2024).
When d' > 4, for sufficiently small, we have that C = O(c,''
Proof summary We now provide a brief summary of the proof of Proposition 6 The first step of the proof technique utilises a fundamental property of the empirical denoising score matching objective, & (s; S, r): that is strongly convex in in a data-dependent weighted T.2-space. Strong convexity is often used in algorithmic stability analyses, especially in deriving stability bounds for linear models- we borrow a similar approach, but we analyse the stability the algorithm in function space. With this, we arrive at the following inequality (see Lemma
19
J [E())((,,,)) s' (X,, (dd) I B[Lam ($)] + Eue (C1/22 + Estab),
(12)
NV
where Estah is the (yet-to-be bounded) score stability constant of Aorm-
The second step of the proof technique utilises a characteristic property of the heat kernel that smooths out functions. In particular, we utilise the celebrated Harnack inequality of Wang (1977) that captures this property by showing that for any positive measurable ő : Rdd - R_, I, y E Rd holds that
Bio(X))X x E(((X))
2/2 exp
|a- 1))p
for any > 0.p > 1. Utilising this bound, we convert the upper bound in 12 to a bound on the stability constant. The full proof can be found in Appendix C

---

Published as a conference paper at ICLR 2026
5 STOCHASTIC SAMPLING AND SCORE STABILITY
In practice, the backwards process in cannot be sampled exactly, so we instead rely on approx- imations based on numerical integration schemes, this section, we investigate how algorithmic stability interacts with these sampling schemes. We consider the Euler-Manyyama-yya sampling scheme proposed in (Benton et al., 2024; Potaptchik et al. 2024) which discretises at the timesteps
(t=) where f = T (1 + K) choose K
for large k. The quantity K. > O is chosen freely and we log (1 + k)) so that t.k = c (see Appendix D for details). By sampling its terminating iterate ÿK, we obtain a sampling algorithm, Aoms that maps score function to the distribution law which approximates the distribution law(Y.).
In the previous section, we identified that early stopping of the backwards process benefits gen- eralisation. In the present section, we will consider how coarseness of the discretisation scheme produces similar benefits. It often the case that the score function is trained only at those time steps considered by the sampler, i.e. using the time- weighteigg (dt() = K Σx= 5r-t, (dt) (oo et al. 2020). As result, the effective stopping time of the algorithm can be much larger than the early stopping time, e. In the following proposition, we demonstrate how this benefits generalisation.
Proposition 7. Consider the setting of Proposition O with = 1 and set T = F. then for sufficiently large N, & S 1/ /4 and any e E (0, 1). we have that for qK Aem Aerm
EID(Pe S B( + B (1 + k))dd
+ C.m
(1 + +x) x=)=d
Be
2d' + x=1 + kdd log(e log(e-1)2 + de -27
where B
Cam Cam N-2 + 1 Cam.n = inf i (h; fx))
The second and third terms of the bound in Proposition 7 are due to the score stability of the ERM algorithm and decay increases. The fourth term of the bound captures the discretisation error and therefore increases with s. What this result captures is that there is trade-off between sampler accuracy and generalisation that is managed by the discretisation of the diffusion model. In the following corollary, this trade-off is optimised.
Corollary 8. Consider the setting of Proposition 7 then for any c (0, and sufficiently small c, there exists such that with qK- = Aem Aerm (S)
B C BrG
Br S log(e-1) -1)2,
[[D(pellak)] S
log(e- -1) 12 + (C + d²)) log(c - + de-2T otherwise.
The primary strength of this result over (Oko et al. 2023 Azangulov et al. 2024 is that we assume little about the hypothesis class. Their results require carefully constrained network architectures and specific early stopping time to control complexity. In contrast, our result holds for any sufficiently small early stopping time, relying instead on a carefully chosen discretisation scheme, which is usually tuned in practice (Karras al., 2022 Williams et al. (024). The main drawback is that our general approach does not exploit the model class to adapt to smoothness properties of the underlying measure, which we leave for future work.
In Figure we provide a toy experiment that reproduces the implicit regularisation brought about by coarse discretisation and carly stopping. Specifically, we track the population KL divergence of a 1-dimensional Gaussian diffusion model guided by an empirical score function evaluated on small dataset (N=40). We observe distinct U-shaped trajectories across the scaling parameter K, the number of discretisation steps, and the early stopping time e. The pronounced minima in these curves demonstrate that restricting the sampling process promotes generalisation by acting as form of effective regularisation.
STOCHASTIC OPTIMISATION AND IMPLICIT REGULARISATION
To learn the score function, it common to choose it from a parametric hypothesis class (so R" } (c.g. a deep neural network) by minimising larm via stochastic optimisation (Karras et al. 2024). In this section, we consider the score stability of this setting, focusing on stochastic gradient descent
7

---

Published as a conference paper at ICLR 2026
wwwmm
0.22
0.22
diseagasis
KL. divergrsiee
0.20
0.12
0.00
0.18
0.18
0.05
0.18
0.16
0.05
211
Discrctianiion stops
Eary stogroing time,
Figure 2: Impact of discretisation and early stopping on population KL divergence for a diffusion model using the empirical score function. Early stopping and coarse discretisation provide form of implicit regularisation. Evaluated on a I-dimensional Gaussian target with dataset of size 40.
(SGD) with gradient clipping and weight decay. We consider the standard gradient estimator: given the mini-batch ( of size NB & N we define the random estimator,
G(O, (z) = No
Wt,,, Vellso(Xisio) (2.2) log Pt.,, (X
(13)
where we define the random variables X1, ue.,, 4 + ot.,, ER. ti,s ~ w, (0,-i((t)) (dt), ) Ev,5 ~ N (0. Ia). The additional variance introduced by the random variables Ei, and ti. leads tor with significantly higher variance than in standard supervised learning. This presents several challenges during training, and various strategies have been proposed mitigate this issue (Karras eet al. 2024 Song Kingma| 2021). 1). For example, the weighting function w : [0, Ti + R. can be tuned reduce variance (Karras et 2022) the number of resamples P N can be increased. We consider the following iterative scheme, defined for a given weight decay constant À > 0 and clipping value C' >
gradient estima-
0K+1 (1 A)0x 7/ Clipo (Gx (Ox, (zi))en,,))
(14)
where 7/s and Ba C [N] the learning rates and mini-batch indices for each iteration k E N and we define the clipping operator (v) = (1 [Cl||e] l)) Both gradient clipping and weight decay are widely used in diffusion model training and are typically motivated by their stabilising effect on optimisation, minimising the impact of the high variance of the gradient estimator (Song al. 2021| Ho et al., 2020). Throughout this section, we take the mini-batch B& to be i.i.d. and uniformly sampled from [N] without replacement. For the sake of simplicity, we suppose that the iterative scheme is terminated after K iterations, where K fixed and independent of the data.
6.1 STABILITY OF SGD WITH WEIGHT DECAY AND CLIPPING
In our analysis, avoid restricting the score network to specific parametric class and instead make structural assumptions based on its smoothness properties. We recall that function is Lipschitz with constant L Z 0 if it is differentiable and its directional derivatives are uniformly bounded by L.
Assumption 9 (Smoothness of the score network). There exists L : Rd x (0,T] + R4 and M Rdd x (0,T) + R, such that for almost all E iRd,t t (0, TT, so(x, is Lipschitz and smooth (gradient-Lipsee Lipschitz) in E R with constants L(x, t) and (x, t). respectively. Furthermore, there exists constants L, M 2 such that for any 2 E supp(vaat))
[X 1) = 2) (dtt < It" JE[M(X, = M
The use of Lipschitz and smoothness assumptions is commonplace in the analysis of optimisation schemes (Nesterov 2018 Hardt et al. 2016). However, the assumption differs slightly from the usual that we only require these properties to hold almost everywhere with respect to the input distribution and allow the Lipschitz and smoothness constants vary with the input, provided their square averages remain bounded. This relaxation enables us to accommodate common models that would otherwise violate global smoothness assumptions, such as ReLU networks.
Assumption 10. Suppose there exists B, > O such that for any 0 E R", it holds that ld (80; (x}, < /3//ot/ for each x E supp(vaataa, supp(r).
(15)
8

---

Published as a conference paper at ICLR 2026
This property requires that the supported score functions made of denoising functions that concentrated on a compact set. To highlight that this can be achieved quite easily, we note that with the naive estimate s(z,
- (15) is satisfied with B [[]]X
In the following proposition we demonstrate score stability bounds in the case that the step size is decaying with rate of 1/k.
Proposition 11. Consider the score matching algorithm Agm : - sex for some fixed K N where (0x & is as given in (14). Suppose that assumptions 9 Jand hold and The S nikk f/ all k < K, for some 9) E (0, )) Then, we obtain that Asm score stable with constant,
stab 2
C x
VR
30
I
C
GU41
NBK no+1
80
(1t)
VI
N
where R² = 11|10. (MB + -X) V and C, = f fo, a *dtt).
Since the score matching algorithm is random, to control the stability constant we construct a coupling of the random score functions (S) and (S*) through a coupling of the optimisation trajectories associated with training on S versus S The construction of this coupling is such that the trajectories are identical for a large portion of the train-time. The proof also heavily relies on the
hence why we consider the setting of Nn < N. The proof-technique is a modification of methodology developed by Hardt al. (2016). Several recent works have explored
stochastic
the influence of first-order optimisation methods on generalisation (Neu et al. 2021 Pensia et al. 2018
Hardt et al. 2016 Clerico et al.) 2022: Dupuis al. 2025).
6.2 UTILISING NOISE IN THE GRADIENT ESTIMATOR
The primary drawback of Proposition II is that the bound grows with the number of iterations. This is particularly problematic since diffusion models often require numerous steps due the high-variance gradient estimator. In this section, we improve this dependence by explicitly leveraging the noise in the gradient estimator. The idea that stochasticity regularisation has motivated the development of numerous learning algorithms and theoretical works in recent years (Srivastava et al., 2014 Bishop 1995 Mou al., 2018 Pensia et al. 2018). Here, we investigate how the noise intrinsic to the gradient estimator for Edsm can play a similar role in promoting generalisation in diffusion models.
optimisation can act form of implicit
To incorporate the effects of the gradient noise, we consider a simplified model in which the noise from the stochastic gradient estimator is approximated with second-order Gaussian approximation:
0k+1 = (1 - --)0 - (G:))( Bg, s] + ry Cov (C:)pc(G:)((:: Bg, 8) 1/2 EA. (16)
where ε E IRd is a standard Gaussian and we use Gg == Gu B&). This approximation can be justified by observing that the inner summation in (13) is over conditionally i.i.d. variables, once conditioned on 0, B and S. Therefore, the gradient estimator G becomes approximately Gaussian as P grows large. For this analysis, we assume the following lower bound on the gradient noise.
Assumption 12. There exists positive semi-definite matrix Σ E Baxn such that for any x E supp(w) and 0 E R". Cove (Clipo - x13)|Xo- = (
Furthermore, the eigenvalues E, (Ae) possess the spectral gap Agap minx X > 0.
We use the matrix E to dictate the geometry on which we perform our analysis. In particular, we consider the weighted norm ||0|| := "' e where Σ is the pscudoinverse matrix.
Assumption 13. For almost all x Id t E (0, so(æ, t) is Lipschitz and smooth (gradient- Lipschitz) 0 E R with respect to the seminorm || l1st and with constants L(x,, 1) and M M(x, t), respectively. Furthermore, there exists constants L, M 2 O such that for any E supp(vasta)
= =) (d)) = JEIM(X x] r(d)) & M..

---

Published as a conference paper at ICLR 2026
By requiring that the Lipschitz and smoothness properties hold with respect to I|s effectively require that the gradient estimator adds noise in all directions aside from those that do not change the function (e.g, along symmetries in the parameter space). With this, arrive at our time-convergent score stability bound for SGD.
Proposition 14. Consider the score matching algorithm A&m + sex for some fixed K E N
where (0k)x & is as given in (16). Suppose that assumptions 70 12 and 13 [3 hold, then there exists some
fj > O such that, if supk 7k I 7) we obtain that Agm is score stable with constant
L'C N
(CPN&C2
K-1
@PNBC22
)
Stab
PNE
exp
1/min X2
min
MR, exp
Mmin X2
k=0
for some constant, 2(M, L, ARAPI Amax (E), Be, C+) > 0 and 7mmin
ming T/k-
In this bound, we recover the 1 score stability bounds from Proposition 11 while also introducing VN the property that the bound does not grow endlessly with the number of iterations. This property is obtained using the noise in the gradient estimator and is not possible without additional noise. Our proof methodology builds techniques developed in the literature analysing stochastic gradient Langevin dynamics, a modification of SGD that applies isotropic Gaussian noise at each step. In particular, we draw on the reflection coupling method of Farghly & Rebeschini (2021), which constructs coupled trajectories of the optimisation iterated that contract in expectation under a suitably defined metric (using the technical results of Eberle Majkal Maji) setting, this contraction arises naturally from the noise inherent to the gradient estimator for the denoising score matching objective. As P increases or r) decreases, the long-term stability bound increases exponentially fast as a result of the benefit of the noise weakening. Through this analysis, we identify the generalisation benefit of a property unique to diffusion models and how they interact with SGD.
al.] (2020)). In our
7 CONCLUSION AND FUTURE WORK
In this paper, we propose general algorithn-dependent framework for analysing the generalisation capabilities of diffusion models. We introduce score stability, which quantifies an algorithm's ' sensitivity to the dataset, and use it derive expected generalisation gap bounds. Applying this framework to several common algorithms, we derive closed-form bounds and identify several previously overlooked sources of implicit regularisation in diffusion models. First, our analysis of empirical risk minimisation finds that the denoising score matching objective provides inherent stability guarantees without further regularisation (denoising regularisation). We then analyse how score stability interacts with discrete-time samplers, identifying that coarse discretisation can improve generalisation guarantees (sampler regularisation). Finally, we consider stochastic optimisation schemes for score matching, obtaining stability guarantees (optimisation regularisation).
This work opens up several directions for future work. This includes identifying further relationships between score stability and generalisation by developing high probability bounds, bounds on notions of memorisation or privacy. The analysis of empirical score matching could be tightened by utilising more properties of the data distribution model class, such as smoothness. The analysis of sampling could be taken further by considering and comparing different sampling algorithms under the score stability framework (e.g. by considering the probability flow ODE).
ACKNOWLEDGEMENTS
Tyler Farghly was supported by Engineering and Physical Sciences Research Council (EPSRC) [grant number EP/T 1/1] and by the DeepMind scholarship. Patrick Rebeschini was funded by UK Research and Innovation (UKRI) under the UK government's Horizon Europe funding guarantee Igrant number EP/Y028333/1), George Deligiannidis acknowledges support from EPSRC [grant number EP/Y018273/11 The authors would like to thank Michael Hutchinson, Valentin De Bortoli, Peter Potaptchik, Sam Howard, Iskander Azangulov and Christopher J. Williams for valuable comments and stimulating discussions. We would like to give special thanks to Sam Howard for assisting in creating Figure] [
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Eddic Aamari. Convergence Rates for Geometric Inference. PhD thesis, Université Paris-Saclay, September 2017.
Eddie Aamari, Jisu Kim, Frédéric Chazal, Bertrand Michel, Alessandro Rinaldo, and Larry Wasser- man. Estimating the reach of manifold. Electronic Journal of Statistics, 13()):135--139, 2019.
Amit Attia and Tomer Koren. Uniform stability for first-order empirical risk minimization. In Po-Ling Loh and Maxim Raginsky (eds.). Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pp. 3313-3332. PMLR, 2022.
Iskander Azangulov, George Deligiannidis, and Judith Rousseau. Convergence of diffusion models under the manifold hypothesis in nigh-dimensions. arXiv preprint arXiv:2409. 2024.
Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and Geometry of Markov Diffusion Operators. Springer International Publishing, 2014.
Peter L. Bartlett, Andrea Montanari, and Alexander Rakhlin. Deep learning: statistical viewpoint. Acta Numerica, 30:87-201, 2021.
Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis. Nearly d-linear con- vergence bounds for diffusion models via stochastic localization. In International Conference on Learning Representations, 2024.
Chris Bishop. Training with noise is equivalent to tikhonov regularization. Neural Computation, (1):108-166, 1995.
Olivier Bousquet and André Elisseeff. Stability and generalization. The Journal of Machine Learning Research, 2:499-522 526, 2002.
Zachary Charles and Dimitris Papailiopoulos. Stability and generalization of learning algorithms that converge global optima. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 745-754. PMLR, 2018.
Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In International Conference on Learning Representations, 2023.
Sitan Chen, Vasilis Kontonis, and Kulin Shah. Learning general Gaussian mixtures with efficient score matching. arXiv preprint arXiv:2404.1 18893, 2024.
Eugenio Clerico, Tyler Farghly, George Deligiannidis, Benjamin Guedj, and Arnaud Doucet. Gen- eralisation under gradient descent via deterministic PAC-bayes. arXiv [stat.ML], September 2022.
Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022.
Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet. Diffusion Schrödinger bridge with applications to score-based generative modeling. In Advances in Neural Information Processing Systems, 2021.
Devroye and Wagner. Distribution-free performance bounds for potential function rules. IEEE Transactions on Information Theory, 25(5):601-604, September 1979.
Benjamin Dupuis, Maxime Haddouche, George Deligiannidis, and Umut Simsekli. Understanding the generalization enror of markov algorithms through poissonization. arXiv [stat.MI.], February 2025.
Andreas Eberle. Reflection couplings and contraction rates for diffusions. Probability Theory and Related Fields, 166(3):851-88) December 2016.

---

Published as a conference paper at ICLR 2026
Andreas Eberle and Mateusz B Majka. Quantitative contraction rates for Markov chains general state spaces. Electronic Journal of Probability, 24:1-36, January 2019.
Patrick Esser, Sumith Kulal, Andreas Blattmann, Rahim Entezari. Jonas Müller, Harry Saini, Yam
Levi, Dominik Lorenz, Axel Sauer, Frederic Boesel, et al. Scaling rectified flow transformers for high-resolution image synthesis. In Forty-furst international conference on machine learning, 2024.
Tyler Farghly and Patrick Rebeschini. Time- ieiee-indeeendent generalization bounds for SGLD in non- convex settings. In Advances in Neural Information Processing Systems, 2021.
Moritz Hardt, Ben Recht, and Yoram Singer. Train faster, generalize better: Stability of stochastic gradient descent. In International Conference on Machine Learning, pp. 1225-1234, 2016.
G Haussmann and Pardoux. Time reversal diffusions. The Annals of Probability, 14(4): 1408, 1986.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, 2020.
Zhihan Huang, Yuting Wei, and Yuxin Chen. Denoising diffusion probabilistic models are optimally adaptive unknown low dimensionality. arXiv [cs.LG], October 2024.
Aapo Hyvärinen. Estimation of Non-Normalized statistical models by score matching. Journal of Machine Learning Research, 6(24):695-709, 2005.
Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion- based generative models. In S. Koyejo and S. Mohamed and A. Agarwal and D. Belgrave and K. Cho and A. Oh (ed.), Advances in Neural Information Processing Systems, 2022.
Tero Karras, Miika Aittala, Jaakko Lehtinen, Janne Hellsten, Timo Aila, and Samuli Laine. Analyzing and improving the training dynamics of diffusion models. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 24174-24184, June 2024.
M Kearns and D Ron. Algorithmic stability and sanity-check bounds for leave-one-out cross- validation. Neural Computation, 1176):1427-1453, 1999.
Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence for score-based generative modeling with polynomial complexity. Advances in Neural Information Processing Systems, 35:22870-22882. 2022.
Puheng Li, Zhong Li, Huishuai Zhang, and Jiang Bian. On the generalization properties of diffusion models. In Advances in Neural Information Processing Systems, 2023.
Gabriel
Brendan Leigh Ross, Rasa Hosseinzadeh, Anthony Caterini, and Jesse C Cresswell. Deep generative models through the lens of the manifold hypothesis: A survey and new connections. arXiv [cs, April 2024.
Mateusz B Majka, Aleksandar Mijatovic, and Lukasz Szpruch. Nonasymptotic bounds for sampling algorithms without log.ccncyvity. Ann. Appl. Probab., 30(4):1534-155 1581, August 2020.
Wenlong Mou, Liwei Wang, Xiyu Zhai, and Kai Zheng. Generalization bounds of sgld for non-convex learning: Two theoretical viewpoints. In Bubeck, Sébastien and Perchet, Vianney and Rigollet, Philippe (ed.), Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pp. 605 -388.68.77 PMLR, 2018
Yurii Nesterov, Lectures on convex optimization. Springer optimization and its applications. Springer International Publishing, Cham, Switzerland, 2 edition, December 2018.
Gergely Neu, Gintare Karolina Dziugaite, Mahdi Haghifam, and Daniel M Roy. coammation-theoretti generalization bounds for stochastic gradient descent. In Proceedings of Thirty Fourth Conference Learning Theory, volume 134 of Proceedings of Machine Learning Research, pp. 3526-3545, PMLR, 2021.
Kazusato Oko, Shunta Akiyama, and Taiji Suzuki. Diffusion models are minimax optimal distribution estimators. In International Conference on Machine Learning, 2023.
12

---

Published as a conference paper at ICLR 2026
Ankit Pensia, Varun Jog. and Po-Ling Loh. Generalization enror bounds for noisy. iterative algorithms. In 2018 IEEE International Symposium on Information Theory (ISIT), pp. 546-550, June 2018
Jakiw Pidstrigach. Score-based generative models detect manifolds. In Advances in Neural Informa- fion Processing Systems, 2022.
Peter Potaptchik, Iskander Azangulov, and George Deligiannidis. Linear convergence of diffusion models under the manifold hypothesis. arXiv preprint arXiv:2410.09046. 2024.
Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen. Hierarchical text- conditional image generation with clip latents. arXiv preprint arXiv:2204.06125, 2022.
Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Bjorn Ommer. High- resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF Confer- ence on Computer Vision and Pattern Recognition (CVPR). IEEE, June 2022.
Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily L Denton, Kamyar Ghasemipour, Raphael Gontijo Lopes, Burcu Karagol Ayan, Tim Salimans, Jonathan Ho, David Fleet, and Mohammad Norouzi. Photorealistic Text-to-lmage diffusion models with deep language understanding. In Advances in Neural Information Processing Systems, 2022.
Kulin Shah, Sitan Chen, and Adam Klivans. Learning mixtures of Gaussians using the DDPM objective. In Advances in Neural Information Processing Systems, 2023.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. International Conference on Machine Learning. 2015.
Yang Song and Diederik P Kingma. How to train your energy-based models. arXiv preprint arXiv:2101. 03288, 2021.
Yang Song, Jascha Sohl- Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Boore-Based generative modeling through stochastic differential equations. In Interrationnal Conference on Learning Representations, 2021.
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15556):1999-1958 2014.
Rong Tang and Yun Yang. Adaptivity of diffusion models to manifold structures. In International Conference on Artificial Intelligence and Statistics, 2024.
Simon Vary, David Martínez-Rubio and Patrick Rebeschini. Black-box uniform stability for non- euclidean empirical risk minimization. arXiv [cs.LG], December 2024.
Martin Wainwright. High--Diensiooal Statistics: A Non-Asymptotic Viewpoint. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2019.
Feng- Yu Wang. Logarithmic sobolev inequalities on noncompact riemannian manifolds. Probability Theory and Related Fields, 10943)441--424 November 1997.
Joseph L Watson, David Juergens, Nathaniel R Bennett, Brian L Trippe, Jason Yim, Helen E Eisenach, Woody Ahern, Andrew Borst, Robert J Ragotte, Lukas F Milles, et al. De novo design of protein structure and function with RFdiffusion. Nature, 620(7976):1089 1100, 2023.
Christopher Williams, Andrew Campbell, Amaud Doucet, and Saifuddin Syed. Score-optimal diffusion schedules. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
Hongkang Yang. A mathematical framework for learning probability distributions. arXiv preprint arXiv:2212. 221121222.1... 2022.
Chiyuan Zhang, Samy Bengio, Moritz Hardt, Google Brain, Benjamin Recht, and Oriol Vinyals. Understanding deep learning (still) requires rethinking generalization. Communications of the ACM, 6473):007-155, 2021.
13

---

Published as a conference paper at ICLR 2026
A FURTHER BACKGROUND
We begin with some further details on notation and lemmas used throughout this work and provide proofs for the lemmas in Section 2
A.1 RANDOM SCORE MATCHING ALGORITHMS
We begin with some additional details on how random score matching algorithms are defined in this work. Recalling the probability space (St, F,),, define the set of random score functions,
S S R x 0, T] x s(..,w) E L (IR" [0,T]:R*)}
For any random score matching algorithm Asm : (UN=1 (R))ON) Ω + ( [0, T) R2), use Asm (S) as shorthand for the random score function (w, I, 1) + Asm (S, (, () belonging to S.
Given random functions 8, 8', let [(s, s') denote the of all couplings of these functions which we define as,
8's, = (5, 8') E x S 8
where s denotes the fact that for any bounded measurable test function Φ: (0,""x0,T]::R* + R. it holds that,
f ( = ((5(:,,))dP.
A.2 PRELIMINARY LEMMAS
For the score matching loss bound, we begin with the fact that the score matching loss is equivalent to the denoising score matching loss up an added constant (Song et al. (2021); Hyvärinen (2005).
Lemma 15. For any t O. y Rd we have
E X X, = V (v)
1 & X, = S, V log
(17)
oi
of
Proof. We begin by showing that the conditional score is an unbiased estimate of log For any Rd 0, we have
/ v log (oo (ojo(y)
E (x, = x]
logPro(Z(1)) Pelo(=/y)oo()) dy /
(z(x)
Po Po(v) dy- (x(x)
Vp(o(zly)
Therefore, using the exchangeability of gradients and integrals (note that Pejo is (?)), we arrive at
Vp((x)
E log (x,)) = x]
(t(r)
Vogg
(18)
Alternatively, using 7 obtain that the left-hand side takes the form,
V log (X = [x X, 2
completing the proof of the first equality in (17). For the second equality concerning the empirical score function, the proof follows similarly once the empirical measure & ΣNn 8,, is used place
of Vdata-
14

---

Published as a conference paper at ICLR 2026
Lemma 16. For any integrable score function 8, it holds that
Laam [dam(s;7 (8; = lem (3; + Cam:
where, given 8* (x, == V (o)), we define
Cam := =
2 d E|TH (oxx (dsm (
(19)
Proof. Let s be any score function. Using the equality in [18 we obtain the following bias-variance decomposition of (dmm(s) T);
(8; = E [|s(Xe. (I((r..))VVopppoo0
Vlogpro
= | 12 [|s(X || T(L)) + (
- V (xo8)((X)))])((())
= km (s:
+
E Cr (Vlog X.) T(dt).
Once we note that,
[uXo x:
Tr Cov (Vlog:1(X)/XX))
Tr Cov
Xx
22 a
Cor(XolXe).
we obtain the bound lasm (s; (8; + Cam from the statement. To derive the equality that (am (sn[s) and so we obtain (asm (s*;)) + Cam-
Cam Ldsm (s*; ;)),
Similarly, there is an equivalence between the empirical forms of the denoising score matching loss and the score matching loss,
Edsm (s; (s) E (s; (s; S, + Cam.
(20)
where
Cam f n T Cov (dt)) Ldssu SS of
(21)
and 6* (Ix, () V log This follows immediately from the above proof once the empirical
measure
Lim1 82, is considered in place of Vdatt- With this, we can now prove Lemma from
N
Section 2
Lemma
The objective
is identical, up constant, to the objective
(s; = [[I,s((r,,)) (ogg((())111)))((())
(22)
where Pe is the marginal density of Xt. Therefore, any minimiser of (: on L (RRd x [0, 1Rd) is identical V log Dr a.e. for any t supp(T).
Proof. The proof follows nearly immediately from 20 Since Prjo is 8 V log Prjo is measurable and thus its empirical average V log Pr must be also. Therefore, the score function s* (x, V log de (x) satisfies 8* E L (R. x [0, Rd) as well as,
Lam S,
Now let s E (R x J0, TO,T); Rdd) be any minimiser of (=: Through the equivalence of lasm and & up to constant, it follows that must also be minimiser of 4am (,; and, due to the existence of s*, must satisfy lam S, 7) = 0 also. Given that whenever > Prjo has full support, it follows that s(,()) s* s* almost everywhere, 7-almost surely.
15

---

Published as a conference paper at ICLR 2026
A.3 MANIFOLDS
We also introduce some basic properties of smooth manifolds, primarily referencing Aamari et al. (2019). We define the manifold reach and include known property of this quantity.
Definition 17. The reach of set A C Rd, is defined by inf d(p, M Med(AA where define the set,
Med(A) {a Rd 3p.q E A.t. p p q, 1p =11 11g =11}.
Lemma 18. Suppose that the measure μ is supported manifold M with reach TM 0 and dimension d', Then, for any r S TM, we have
(x),
inf Pp
B.(a)
where Py. denotes the density of y with respect the volume measure M.
For the proof of this lemma, we refer to the proof of Proposition 4.3 in Aamari et al. (2019) or Lemma III.23 Aamari (2017).
B PROOFS FOR THE GENERALISATION GAP BOUNDS
We now provide the proof of Theorem that bounds the generalisation gap under score stability guarantees.
Theorem 3, Suppose that the score matching algorithm Asm is score stable with constant e Then, with As)SS), holds that
[ (8; []: (8:7)) E[Rdsm (3; (,) S Estab.
(23)
Furthermore, it holds that
E Blesu (5; [S: B[lam (8; S8; $ 2 stab TE [[Edsmm (s: (.,, ($:S, +
(24)
Proof. Setting & Aum (S) and s Aam (S*), we use the property that (8, x) and (st, .x) are distributed identically to obtain that,
[ffdem (8; r] B[Easm (8; (f),
aasm (8'; (are):))
Ex, [xx - V log (x,/ xs, ((d()
Therefore, it follows from the triangle inequality in norm that,
N
1/2
B[easm (// (8
E
N
Note that if the algorithm A&m is stochastic, the right-hand side would hold regardless of how a|S, a and 8*| "S were coupled. Therefore, the most efficient coupling can be chosen, leading to the bound,
[eass (22 E.ladm (8; S))
(25)
N
1/2
< E
inf (e,
ElIS(X,, - 8'xx. SI(()
N
< Ssaab
(26)
completing the proof of the bound in
23
16

---

Published as a conference paper at ICLR 2026
To obtain the bound in 24 we use Lemma 16 derive
E[fam (8;
E fam (8; S, r)] El (8; - fasm (8; S, (stm log Pr; S)
Edsm (V log Pt;
(27)
Since ldsm(t) T) is unbiased estimator of lam (:T,, we have that
aasm V log 7) = B[fdsm logge; Pt; S, 7)] 2 E (Cttm (V log Pe; 7)],
(28)
where the inequality follows from the fact that V log P. minimises lasm- Furthermore, using (26). we deduce the bound,
[E[[asm ($; T) Lass (5; S, T)]
E[eaam (3; T)]//2 + Ellaam (s; (3))/)) B[eaam (S; Ellasm (,;
22leasm (8; +) + Estab Extab
28 B( l (8; S, estab
(29)
Thus, substituting
28 and
29 into
recovers the bound (24) in the statement.
We obtain upper bounds relying on the fact that the constant separating the score matching loss from the denoising score matching loss is larger on average in the empirical case. One could obtain lower bounds through our techniques, but this would require an analysis of the rate of convergence of this constant, which is beyond the scope of this paper.
C PROOFS FOR STABILITY OF EMPIRICAL DENOISING SCORE MATCHING
In this section, we provide the proof for Theorem where the algorithm that minimises ldam (c; S, over some class of score functions H. is shown be score stable.
C.1 ON-AVERAGE STABILITY OF THE ERM ALGORITHM
We begin with an important lemma that shows that under minimal assumptions, 8 Aerm (S) and
8 Aerm (S*) are close L space, averaged over the full dataset. The first half of this proof utilises the fact that fasm strongly convex weighted T,2 space, exploiting a well-known relationship between strong-convexity and algorithmic stability (e.g. see (Bousquet & Elisseeff 2002 Charles & Papailiopoulos 2018: Vary et al. 2024 Attia Koren 2022 ).
Lemma 19. Suppose that Aerm score stable with constant Estab, then for any E [N]. we obtain,
[[ 1 - t)1 pa (dy) r(dt) < 8E (8)] +
/
8 Esual (C1/2 + Estab)
EE
(30)
N
where 8 Aerm (S), ASrm
Proof. Choose E [N] and let = Aerm (S), ASrm(S') so that & E argming (,T,,), argming (::Si,7) The proof begins with the following simple expression, which holds for all
j E [N]:
/ (z) s(y, s(..t). logp)p()12,)) ()o((yyaz,
- [18((u,)) log (()o(dy)x;) [],(y,t) Vlogpio(vlZE))1 p.jo (dyl2))
/115/1.00 s(v. 1) ( Pro
17

---

Published as a conference paper at ICLR 2026
By averaging over [M] and integrating with respect to 7, we arrive at the upper bound,
2 N
| (#(u.t) $(v. 1),8 (v.
Σ
V log log Peo(dylzr)) r(dt)
JE[N]
Edsm (8):S,T) - ldsm (8; S, + / [15((v,)) s(v, t)11 (dyy r(dt)
2 / J [1=((,,) () (ldd)) (c
(31)
where the inequality follows from the fact that lasm S,T) S lasm (s; S, for any score function E H. Additionally, the left-hand side is upper bounded using the Cauchy- Schhwrr inequality to obtain.
2
Σ
(8..) 8(v, 1) Vlog (noo)
N
PES
2
Σ/ /
- ((0 po
N
RES
2
/ / (6((v. - (v,,1),8(,,,1) (1)))))())))))))
+
N
2
/ / (6 - (6(()tt)SS(utt))
- v (
N
S 20 (8'; ,7 ( (18((1,)) -i,. 1)1 (dd)
1/2
2 la (8': (2.). 7)2/2
( (l)t $(y. 1) p:jo (dyl) r(dt)
1/2
N
2 Raam (8'; /
1/2
[18 S(y, 1)1 p:o(dyl)) T(d()
(32)
N
where (dy) = N Enss Peo (dyl)) Combining the expressions in expectation, we derive the bound,
31 and
32 and taking the
E E[ / (g) 8(y, f)12 (dy) r(at)
< 21 (s'; St, /2 / [18((y.)) - s(y, ((dy)) r(dt)
1/2
1/2
B E dam (8'; (1,)) ( /2 / |1st -8(, t)1 Pr pr r(dt)
1/2
2 Bl Bldaam (8); (3), 1/28 / (1,((t,t) - +(1, (ojo(dy))) 7 (dtt)
21 (3; S, 12 / |1st (v. $(v, t)]? (c(d)) r(dt)
1/2
(Bleaam (8; S, 7))///2 + (B: 1122)
+ N Sstab
where we recall that estab is the stability constant for Aerm- Here, we have used the fact that (8, S) has the same law as (8*, St) and also Bllddm (8); (i), r)] Bflassm (8; S, r)] and ElLdam (8%; (xre) 7)]
18

---

Published as a conference paper at ICLR 2026
Bledsm (8; 7).. By solving the quadratic equation, deduce that the above inequality implies,
(v. s(y. t) pr(dy) T(dt)
2
Bl (s: S. + [[esmm (8; S.T] N Establ (E( (8; r)] Elldsm (s; (8;
< 4B| (8; S,
N uab (E(easm (8; r)} + ElLatm (8; S, 2)
We simplify the above expression further using Theorem R
Using the stability assumption, follows
from (23) that Ef ( ($)) & B[fasm (8)]2/2 + Estab. Furthermore, from Lemma 16, we have
E Eledsm ($)] = Elean (8) + [[[am]
[B am (s)] + Cam.
where we recall the definitions of Cam and Cam from (1 and 19) and recall that B[Cam] Cam from (28). Thus, from Young's inequality, we obtain the bound
E
(1, - -8(,, 1)1 -B(y,1)1l (o(dd))(d))
+
stab (2E (8)) + 2C + Estab)
NN
S 8E| Rum ($)
+
Estab (Estab/ /N N + 2C + Estab)
N
< 8E (5)
Estab (C1/2 Estab).
C.2 PROOF OF PROPOSITION
To obtain the stability bound in Proposition 6
we convert the result in Lemma 19 which is a bound in T.2 to the bound in T.2 (P (ā)) required of score stability. For this, we rely two further process, captured by the Harnack incquality of Wang (1997) (see Theorem 5.6.1 in (Bakry al. (014)).
lemmas, the first of which a fundamental property of the Ornstein-Uhlenber
Lemma 20 (Wang's 's Harnack inequality). For each positive measurable function & Rdd - R. every t >>> 0.p > and every x, E i it holds that,
[ 11 y 2(p 1 1)0?
Elo(X.) Xo x1 < E(o(X,)'IXo // exp
This result describes the stability of the diffusion semigroup under changes its initial position and shows that as & grows, the distribution of X, depends less on Xo. The second lemma, for which provide a proof, controls the empirical measure,
1
((zx)
8z, 82, (dr),
N
on balls around training examples. Lemma 21. Suppose that Assumption is satisfied, then for any i E [N], r E (0, Treach and any decreasing function o (0, (0) > R2,, we have the bound
E (l s 2(N exp ( -N + 6(cor""
whenever N (Ac71--d) v 2, where cu inf Pe-
19

---

Published as a conference paper at ICLR 2026
Proof. We rewrite the object (z:)) as empirical average of Bernoulli random variables,
N
(x,)
Σ N
1
1 N
+
Σ Lx,e N
1
When conditioned on 24, the random variables (sa)), are independently and identically distributed Bernoulli random variables with probability μ Vdata (z+)) To utilise concentration of the empirical process, we first rewrite the probability
((D(B,(a.)) - w/2 = (P S - 1
where SN-1 Eij (x,) Since μ S have Np/2 (N 1)p/2. 2. and hence
P(D(B,(t.)) S 4/2 x:) < IP SN S (N - 11)
Therefore, by the multiplicative Chernoff bound for binomial random variables (e.g. Lemma 2.2 in Wainwright 2019)), we obtain,
P(SN-I s (NN-1)//2r; S exp
(N 11)/88
I exp (NNp/16).
for all N 2 2.
Next, using the trivial bound ((Br(I:)) 2 N and the fact that is decreasing, we have
E [(o(D(B((x))) m: = (B.
(0(B. ( zi
+ +E[o(D(B,, (x:)) (10BB zi
s o(N O(N (IBBB)() S 11/2
s 0(p/2)+ O( exp ( (NNp/16)
To control use Lemma 18 which asserts that μ 2 card" and using again that @ is decreasing, we obtain
Substituting this into the above bound
B (((()(B(())))) mx] s a(cur (2) +DN=)N exp (- C., /16).
Taking the expectation with respect to x; yields the desired bound.
This now brings us to the proof of the proposition, which we first restate.
Proposition 6 Suppose that assumptions 4 5 and that inf supp(T) E (0, Theach): then for any c E (U, 1) and sufficienlly large N. the score matching algorithu Aer is score stable with.
stab E C(CC N 2 (( (3)))
D2 or coo
Proof. We use the shorthand & (s) = & (s; S, S,,),,,,,,)).),, la (s) Ra (s; S), (s) ( (5; for the sake of brevity. We start from Lemma 19 which provides bound on the difference between s and & in (p:)) We use it to develop a bound in L (+)))) as required by score stability. In particular, we define the quantity,
IE [
113
[1'( s(v, t) (nu(dy) T(dt)
so that, due to the fact that Aerm invariant with respect to dataset permutations, must be score stable with constant (we have that < from Assumption 5p. Therefore, from Lemma T9 have
/ 118 (v. s(y, t) p. (dy) rddt) $ S ( (8)] + N e + e).
20

---

Published as a conference paper at ICLR 2026
We proceed using Lemma 20 with (())
(v, -8(y, t)1? to obtain that for any j E [N]
p > 1,
/ (15((1,)) - s(y,t)]l t)l (ejo(yy)az)
(ʃ (f11st (v. - s(v. 1) p prlo(dy))))
1/8
[n [ 1: - 23g 2311 112
exp
2(p 1)02
Setting B = B2(22)) O S. we can average over the above bound to obtain,
/ (1,) s(v: 11 ((yo(yy)
1/8
RR
$
B
(1²( s(y,tt)) pro(dylzz)
exp
2(p - 1)o7
ZEB
1/p
1
|18 (v. 3(v, 1)1 p(jo(dy)z)))
exp
2(p 1)0
REB
We can further bound this using the 0 bound in Assumption 5
f |18t (u, 1) s(y, 1) Peoo(dy(z;)
(1B--1100 1/0 (/15((0.t) -|( $11, ba(du)) 1/r exp 2(
RRR2
(1)o/
(Dzl (D2/00)30-11000000 1/1 (f -) S(y, 1)1
1
n²R2 R2
pe(dy)
exp
2(p 1)07
Integrating with respect to r and taking the expectation, we obtain,
2² < (DH/o? 2/ ((B) w/2
1/9
/ (1) t)12 (c(d))
1/0
RR
exp
2(p - 1002
1/4 [8P(2am (5) + N $ SD /0 RR ((B) 4/p)
8
eCCl + +
1/2
RR.
exp
2(p
1)o?
where we define = (1 - 1/p)-1 > 1 and assume that is sufficiently small, so that o² s DH- Using Young's inequality, follows that for any X >>>
e² of و6
D E
(B)
4/p //] exp
que R² 2(p 1)0²
+ X"
8 (SE(Lmm (3)
ECCl/ +e))
N
Setting = SA8 /pN, we can rearrange this to obtain the quadratic inequality,
4/v D2 Di
914 R²
(1 k)= - 11/ Ke <
4/1)
+ NN (8)) 2(p - 1)o?
exp
Npx
029
Requiring that / 1 1/2, we solve the quadratic to obtain the inequality,
2
4 S Cumk? +
8 Nps
4/p
D2
B
((B)
v/rl
exp
quº R22
+ NR
(33)
o.4
2(p
1)o?
Next, we optimise R by setting R = a. and apply Lemma whenever a. < Treach we obtain.
21 with ((r) = --/ to obtain that
coNo
2 9/p
E (BB
4/p)
< N1/r exp
16
C,oo
2
9/v
2
cyort
21

---

Published as a conference paper at ICLR 2026
for all N sufficiently large (depending only cv. oe, p and q). Returning to 33 it follows from the above that
6² S
(
16
q/p
2D 0/4
2q
Cam
H
exp
+ N.E (8)
(34)
4
p-
We now choose & by optimising the second and third terms of this bound, by which we arrive at the choice,
2D
2g
16
9/p
Npccod" for some 7 > O. For any set of constants, we can take N sufficiently large, to ensure that x S 1/2. Substituting this into (34 we arrive at the bound,
J2pNN
exp
p
2
S Cam (Np)
2D
H
2/9
exp
4
16
2/v
7214
2/9
Cood
2D
1/4
2
16
1/p
1/D
1
H
exp
coo
q
+
p/2//
Bleam (3))
Applying Young's incquality to the three terms on the right-hand side, we may choose y > 0 so that,
8²
2D H
16
2D
16
4
exp
p-1
2
Goor
Cam
N
exp
pp-1
2
Currr
BI (8)
12
16
4
2D H
16
Com-N + Ellem (3)]
exp
cyood
p
Optimising p, we obtain,
2D
16
5
log(a -1)1/2 log(a-1)1/2 -2 =))
0
Good
exp
2/2
where,
22
16
H
Cam-N-2 + Eleand (amm(3)) [3)]
of
coord
from which the bound in the statement follows.
D PROOFS FOR SAMPLING AND SCORE STABILITY
In this section, we provide details for the discretisation scheme considered in Section 5 and give the proof for Proposition 7 and Corollary 8 In the work of Potaptchik al. (2024), they consider the following discretisation scheme, based on the scheme of (Benton et al. 2024):
ūk+1 = Mi,1 1, 2h + 1
2
OT OT-
s(Ok. T - +x)) + 012+1-4
Sk. OT-4
k (0, K -1,
-
where CE -N(O.L)) La) and we recall that the timesteps
(tx (+x) are given by,
kk.
if < T-1
t6
T
(1+
-
if T-1 < k < K,
where L -= T-1 0, K L log(c 1) log(1 + +n] and > 0, is chosen freely. We
recall the following result from Potaptchik et al. (0024).
Lemma 22. Suppose that 1 and Assumption 4 holds with diam supp(vaata) $ 1. Then, it holds that,
D(pel ( S lam (s; + ((prlpps) + AR.K: A = k + d² (K- (log(( -1) ) + sup log(p.) )) where we define the measure,
K-1
+(dt) = K
Σ k=0
5 (dt).
22

---

Published as a conference paper at ICLR 2026
D.1 COARSE DISCRETISATION AND REGULARISATION
Fix > 0 and suppose that k. is such that log(e / [og(1 + k) is an integer. Set K = log(e ()/ [o (1+ + k) that, according the discretisation scheme,
tk - (1 + k)) -oofe /) log(1+a) =T
Proof of Proposition Let $ Aerm (S). We begin with Lemma 22 which provides the bound, (pD (mem x Bllsm (8; S,+) D(prllpee) + A.,K-
For sufficiently small we have the bound,
log(1 x) log(e (log(e sup Iog(p=)])
Ax, = K d' n
S <<1 + r)dd 1o
Using Theorem 3 we obtain that if the algorithm is estab score stable, we have
El (3; f)] S [l (,) + 1saan (8; S, f)) + stab
x ElC (5: f)] + stah sm + =
Using Proposition we obtain that with i, Aerm score stable, with constant,
Suh 2 C(CCsm N-2 (Bimms())) (8))
S S o' T- ( (c) c dt --x-x1 N=m Can Bllam (8))
Now by definition, we have that
T so if we take sufficiently small so that 6(1 + x)) 1 we also have 0² ((14m)) 2 6(1 + x)) and thus obtain,
= +
stab o 1 12-11111211111 R) 1r 12 (c c 1z) -=) 12 C El (8))
S c -d (1) d' (c Nam Ellau (8))
where in the last inequality, that c(1 + S 1/2.
We now proceed by proving Corollary g in which the bound in Proposition 7 is optimised.
Proof of Corollary g Let 7, denote the weak limit of the measure Tx as K - 0+. Since supp(f.) S [c. T] and 6 > O, we know that inf & (1 < co. From this, we deduce that lim B, < 00.
With this there exists n² which is the smallest quantity satisfying.
B 1.
+ x* x2))dd log(e 1)2
In the case that B.- > log(e-)) )), we have that
1 (1 (1+xt)
+
Bx-
(1
+
+ + * + k")dd log(e-1)
Cam
F1 R
log(E -1)===)++ + (C_-1 + d' B3 क
< B |(
log + (G + d" ) B गम log( -1)2
Plus, if Bx- < log(e--) and therefore x* 1, then there exists such that,
B Cam
Bu Cam
B3/2 + x:) +
(1 (1+ x(1 +
log( (oog(--11) S B1/2
+ de
Combining these leads to the bound in the statement.
23

---

Published as a conference paper at ICLR 2026
E PROOFS FOR STABILITY OF SGD
In this section, we analyse the stochastic optimisation scheme in (14, deriving the score stability bounds given in Proposition 11 We begin with basic lemma that follows from weight decay and gradient clipping.
Lemma 23. Suppose that Teke x-1 for all lo E N, then for any K E N, it holds that
Ce
||0x1 - A
v |10oll--
Proof.
We begin with the bound,
10++1 ||
(1 m/ 0> |1 + mh (ll (aa)ee)))))
I (1 mh 10 + M.C.
By comparison, this leads to the bound
K-1
K-1
K-1
|10x11
S
C
7
II (1
n2 + II (1 A1)llll
k=0
K
h
R-1
CC
Th exp
ΣΣn 1=0
1/
+ eep
A
Σ Thk 11001
k=0
i=0
K-1 Σ
K-1 Σ
K-1 Σ 7 k=0
C
Th +
À max UK
?k exp
>
A
7
+ exp
|10o1l
Since the sum forms a left Riemann sum, approximating an integral of an increasing function, we can upper bound it by the integral over (xp(At). Furthermore, we have that À max; 7/s to the bound,
S 1, which leads
K-1
exp(At)dt + exp ( (-x
K-1
|0 |1 S Ceexp
-
Σ
1
Σ m
|1001
i=0
i=0
Ce
(-
K-1
Σ
mm))
(-xΣ m/6
K-1 110.11
exp
-exp
==
Ce V |10|||-
We now ready prove Proposition
Proposition 11 Consider the score matching algorithu Aam S + s0g for some fixed K E N where (0x) k as given in [4]. Suppose that assumptions Aand 10hold and qa S iik all k < K. for some 7 E (0, Then, we obtain that Agm is score stable with constant,
stob ( VR
(VVR)
1+
I.²
NBK nUFI
(1t)
17
N
where R = E| v ( + +I- -) and C, = f fo, T(dt).
Proof. Since the stochastic mini-batch scheme, and therefore the resulting score matching algorithm, is symmetric dataset permutations, consider stability under changes in the th entry of the dataset, without loss of generality. Let 0g be the process given in (4), using the dataset S and let 0k be the same process using S instead of S:
0 = (1 n)>)0 7h Clipo (Ō,, (ā))
0 = 0o.
where ă; x; for + N, HN = By having the processes share the same mini-batch indices Bş and gradient approximation G2 (i.e. sharing the same random time variables f.,j and noise E. we couple the processes 0g and 0z.
24

---

Published as a conference paper at ICLR 2026
We proceed by first controlling the stability of the gradient estimator, computing the bound, (() Gk( (#)))BB.)l
1 NBP Σ
Whi, |vo (Xa,,, t2.) (so, (Xis, "
V log Pt. (X
IE IEBL j=1
Vsg. (Xtnss tiss) (T (s0. (Kun)t Po. (N,,,, |z:))\||
(IV (Xen tis) - Vsgz (Xer, t [.l]ll[ser (Xes, tigs)
I NpP
Whit
iE BN 1=1
log ol o(X (xxi|| + V| (Xu,, (a.j) ( |l|so (Xus, (is)) Sb, (Xesys des)ll)
Whi (M(r,,,, ( 10 [z))]|
NpP
iEBN j=1
+ ((xennr4es) 100 0zll-
We control the expectation of this by first noting that,
E Wes,, (M t.,)llse (Xe,, :t,s) V Flog Pogg |x,)i| + L(Xes, (ao)) t2)))) 0k 02 S.
$ (ʃ 1)|20 =x)))===)) 1/2 (feammoon (so tzit, de) (dt)
1/2
1) - xz]T(d()
S MBB.
where we define the quantity Cr := f S,, (ddt). From this, it follows that
E || (x)))EB.) Gx((k,, (2:) 0x. Bk, 5, (MB.C//2 + +T.') 110 0xlll
Furthermore, we can control the difference between Gu(ōk, (x;)seB)) and Gx (i:)sEB)) )r using the fact that they are identical whenever N & Bg. Thus, obtaining,
E [ (x)))EB.) - Clipo (2:) )) 01 0z. S. S.
s |IG(On (xi)ie - G(On. ( 0., S, + +E[ICCipo (I)))))))))) Clipo
Ox. og. S.
$ (MB.C + L²) 1|0₂ 0all -C No
where we have used the fact that P(N E BB) Ne.
Na.
Thus, using 14 we obtain that for any ko S k.
E B[10kt1 0x+1l| 0kp: BRo S.a
(1 + m/ MMBC + +T²² (2=--)) B[100= 0z
0Ko
BRg S.
+ 2m Nn
N
< (1 + m) B 110 0 0e1l 0ro S 2m C
No N
where v C./2 +I² X. Thus, by comparison, we obtain,
K-1
K-1
K-1
NB N
|0x - 0ko ONy S < Σ
2m
II j=+
(1+ nyu) + 100,,
0mall II
(1 + no)..
1=ho
g =
25

---

Published as a conference paper at ICLR 2026
From this obtain the following:
NB
K-1
K-1
E B|10K 0x 10
:
20
Σ Th exp
Σ
N
immo
j=it1
2
K
K
6v
N
km
CNB
K
fju
Nu
kig
where we use the fact that 2j=it1 3 S log k) - (i).
By the law of total probability, we have
0|| - 0xx1000
= =B[10x - 0x11 = 01o [OOk = -ho (0)) + E1| - 0x # 0ko 0o)P(Ono # 0ka 10))
CNB K
- v 11601)
Ce
KoNB N
No
KD
+
where in the second inequality, we use Lemma
23
Thus, optimising ko leads to the bound,
0k 0x 1110 - (-) (1 (1+ (1+
Ce
NB K ..
1/cv)
|1661l)
Finally, we obtain score stability using the fact that
f
E| Ellssx () (Xe, 101111,, ā, S]r(d))
- L2110k - 0x11
Ce
$ 2E
I
[16o1|
110x - 0x1]
A
Ce
C
NB (1+1/co)
$ T.2
v
Keot1 N
A
where R² = 2
F WASSERSTEIN CONTRACTIONS
In this section, we derive the Wasserstein contraction result used in the proof of Proposition 14 We begin with the more abstract problem of deriving Wasserstein contractions for a discrete time diffusion process with anisotropic non-constant volatility. We consider stochastic processes given by the discrete-time update,
1k (1 1) r + 1)b( + (ara))
(35)
uk+1 = (1 - 17 + ) + v (nōō(y))
(36)
for some b, b. R. - I. o, R Rx where & N I,),, and we show that the laws of JA and l/k contract in Wasserstein distance. We borrow the strategy developed by Eberle (016) and extended in (Eberle & Majka| 2019 Majka et al.) (2220), constructing coupling and a metric for which exponential contractions of the coupling can be obtained. However, these works are restricted to the setting of isotropic noise with constant volatility (i.e. (xx) = cI,) and so some careful modification to the strategy is required. In particular, we analyse this process with respect to the pseudoinverse
seminorm || ||a- given by ||=1|2 = x- G where G denotes the of the matrix G. Furthermore, we allow for 3y and y/s have different bias and volatility terms and so controlling for this will also require some modifications to the proof technique.
To define coupling we first suppose that there exists symmetric positive semi-definite matrix d
G E R such that 2 ,y))²²² > G for all x E R". We define the residual volatility,
Ea(x) = (x 1)2/2
Ec(1) (ō())² 1/2/2
26

---

Published as a conference paper at ICLR 2026
which are well as real PSD square roots. Since Z' and Z independent, the noise [x)()) Z' + /nj G1/2 has covariance n)E G) ( recovering the correct law. We construct the coupling the matched setting where both processes use the same drift b and volatility No: the mismatched case will be handled by a perturbation argument in Section F22.4 Starting from z, E I", we define the deterministic drift step and the residual noise step,
i = + ;= ž+ =ー Z
j (1 n>jy ryb(y): ý = (v) Z',
where Z' ~ N(O, We then define the synchronous coupled processes,
X' = 2 vT
Y 9 + /
with Z ~ NO0 I.). We also consider the reflection coupling.
Y = / + /2 (1=
with E (6 00/1 - |
(37)
which has the noise act the mirrored direction. We combine these couplings to arrive at the final coupling (X', Y'):
X
if (x ('), (e, 2)12 m² /n/nn and S o (X') / (X', (ee, 2)12 m² /n and f S ml
Y
Y
i >
(38)
Y
otherwise,
for some fixed m O.
We assume the following regularity properties.
Assumption 24. Suppose that b is bounded, satisfying == suPaER= [Mx]] < 00 and we have the Lipschitz property. [1b(x) - b(y)la+ S Lo|l yllo+ and [[a(() Ea (y)ll S La || - Bll. for all x, IR" for some Lo. Lo 2 0.
We also allow for + b and a + o. making the following assumption.
Assumption 25. Suppose that b, b satisfy [lb(x) b(x)) s B,, [[Ea(x) Ex())llnna I Bo for all x E Ihn and for some Bo, B. 2 O.
We define the objects.
R= | y ylo+,
-= l Vllo. -= -
R' 1X''X''lca''
We wish to show that R' contracts expectation, i.c. is less than R on average. We modify the metric to guarantee this is possible. We define the function,
1 - e-a),, f(r) =
if < r2 1 - e- + 2r2 e ar - otherwise,
where a 6Lgr1/cu, F1 4(1 + Lo
= r1 No and co. To are defined below. The analysis in Majka al.
coupling and the strategy for proving contractions is closely based (2020) and for the sake of comparison, we rely on similar notation. We will also heavily borrow properties of the function f that proven in their work.
By allowing Eg to be non-constant (i.e. when a is non-constant), we run in to additional complica- tions that are controlled by making the following assumption about the scale of Lg. Assumption 26. Suppose that the following three inequalities hold:
(6 (4a)))= VT(1 e-ars )
n - 1,2 - 1)² 32
I < A/8n,
for some universal constant Kg-
Under these assumptions, obtain exponential contractions.
27

---

Published as a conference paper at ICLR 2026
Proposition 27. Suppose that assumptions 24 25 and 26 hold and m /nuu2. Suppose further that the mismatch is sufficiently small:
cpp
BB + vmn Bo -
(39)
256
Then for any 1l S y/ and x, y E R it holds that
BU((IX' - -Yll++)] (1-
1 16 ((R) + 2inBo + vm Bo);
where Y = =Y + w is the corrected output defined Section F22.4
X
Le 3
r2
X
-60r//00
3Lr1
c:= min
16
[1 ar
16'
200
/6 700
16
20p (ollll//))))))
3/2)A²
4B²
(ogo(22))
To == min
AL X 2L
4 1
22
2304L 132
4(A + La + 1) A
for some universal co and 2(Lb A)+ 4n 1/2 v 2(n 1).
F.1 THE COUPLING
Before provide the proof of Proposition 27 we provide explanation of how the coupling is arrived at. We begin by discussing the one- mineemnnsioonal coupling of the Gaussian distribution that the construction is ultimately based on. Consider the following coupling (t, 17) and ((s 7) for
t, s R: with 2 ~ N(0,
+t=+ vTz
(40)
if (t. (t'), //1zl m, and It s Ti if > 2n (t) (t'), ( Tpp m, and s S1
m2.
(41)
8 + VnZ,
otherwise.
This coupling has the following property given in lemmas 3.1 and 3.2 of Majka et al. (0020)
Lemma 28. For the coupling defined in 40 and (41) we have
EEIt l'll I s/,
and i S 4m we have
E B[(1t' 11-1t-s))11e-111111 Co ( - s))
1
I, (0. r + vT), (r ,.r), otherwise,
if vTI
where
for some universal constant co > 0.
Thus, through the second bound, we have control of the probability that |t' s'I contracts below s|. The coupling proposed in (38) is multivariate analogue of this that also accounts for the diffusion coefficient G1/2, Let the vector e Rn be defined in (37). then we obtain that,
(e, = Gta) + VT((G" Z), (e, G G Ø) + (2)
Therefore, G'), are a synchronous coupling of (((e,G+a)). and (lee, G (),)). 1). Furthermore, we have
(e, () + ( = (e, + /2( e, 2/7l,, G ()) e =c, -V((()// Z),
(G
28

---

Published as a conference paper at ICLR 2026
and so (e, G X'), (e, G is the oinedimnnsi reflection coupling. Finally we obtain,
Pggn (X')
PR (X')
1
1 27
IıTz (G ill +
VTTz /
exp
2ny
1
11 +1 ! (
exp
21
1 (0 (g - -))2 / 3))(((/2)t).2))
exp
2m
(VTle, (c. (5- +))) + (TT(G\//)t 2)12)
= exp
2n
2m
01) (VV(()
60. (G Z))
ii ((e
( (f,
From this, we deduce that (e, G x'), (e, G+Y") are coupled as in 40 follows by setting
41
The equivalence
= Ge,
s' s' (e,ee)
(42)
(e, Gie, 8).
s= = G) z= (()
(43)
Through this equivalence, we can extend the previous lemma to obtain the following result about the high dimensional coupling.
Lemma 29. For the coupling defined in 38 we obtain that for S 4m?, we have the following:
E[R' f,,
E (R' f) Mer z²|
1 Co min (ainnvi,f)
2
where co and L, is as in Lemma 28
Proof. Let {er be basis of R= with respect to the inner product (- with e₁ Then, have that
(R')²
(e, G (X'-X
==1
= 11' 8/2 + Σ G (X'--X'))
(44)
i=2
where t', s' are as defined in 42p. For any i * 1, we can use that 1 to obtain that
(e,, G G
(e., G g) + eT, + 2 Z)
(e. G () + 2//T2.
From this, we obtain that,
(e,,
= (e, G - 0)) = 0.
This also holds for the synchronous coupling and hence we obtain (e,, G (X''''''''''' (X'-Y')) = 00. Combined with (44), we obtain that R' = Il' 3'). Similarly can shown that == |1 - s| and thus, from Lemma 28 the statement of the lemma follows.
F.2 PROOF OF PROPOSITION 27
Since the coupling is defined using the matched dynamics (b, Eg) for both processes, we first carry the full coupling analysis in this setting, establishing the contraction [R') (1 ff/8). where R'
Y'Ilo+ + is the matched coupling output. The mismatched case is then handled in Section F F2.4 by perturbation argument.
X'
29

---

Published as a conference paper at ICLR 2026
F.2.1 CONTRACTIONS IN f
When R is large, we can rely on contractive properties following from the weight decay. For this, obtain the following.
Lemma 30. Suppose that Assumption 24 holds. Then whenever R 2 4B/X, we have, (1 m//4 + njo)
(45)
where we define,
AL |el|c |ell|c- -/8)///4
- 2/8)
Furthermore, if < 4B/X, X, we have that,
i P (1 + nLo + na)R.
(46)
Proof. From the triangle inequality and the Lipschitz property of b, we obtain
i S - n)>) æ ullo+ + n)lb(x) byy)la+
< (1 + vLLs R)
Next we bound using the decomposition,
S - 5 + V - Eo(1))
- ععا (z) - Ec(y) on lla
(1 77(4, A) + VnLo 12 ala )R
where in the final bound, we utilise the Lipschitz property of Eg. We can refine this bound, using the orthogonal decomposition,
z'= z0= To + +I--vu vz,
From this, we obtain,
2² = || + vT(io(z) - Ev)))/'12
< || ÿ + Vive"" (Co(()) Σc(r)) + ml[(1 vo (Zo(()) Eo(1))2'1a+
< (f+ v (Eg(z) En 1--) + - Ea (1)
$ + 27 2
u0 ue (Eo(z) - En()))2'lla + (Ea((( - Ec(u)) 1'
In the case that Ex(x) # [v), we can simplify this using the notation Ç == |1(Ec(z)
Ea(y))ull (Ec(z) - E'()))Z' we obtain
2 < 12 +
< (11 + (L 41/ 12, R².
From this, we obtain that
1 + ry( A) +
ALa
||ll| -118/) + (1.1171111
-1/8)
1 + TT + 17
//
(1 nLo n,)R,
which gives (46). Alternatively, we can use the fact that |lblaa+ & B to obtain, (1-1)))
In partieular, if R 2 4B/A, obtain F < (1 12 From this, we obtain the bound in (45)
f < (1 7
We will also need the following properties of the function f-
30

---

Published as a conference paper at ICLR 2026
Lemma 31. The function satisfies the property that for all r 2 r
(()
Proof.
For 2 r2, this is Lemma 4.3 of Majka al. (2020). For ri Sr< since is concave on [0,
[()) - ((11(1)))) (r) Ar e a
The ratio Ar///4 f(r)
4 is decreasing in r, so its minimum on [ri, is at
F2, where
1-e-ar
it equals arze-ar 1-e-ary
A From the definition of we have S 2(1-6 "rz] To = orse ar2
4e-aray A
4.
2(1-e-ar3]
To so
aree
À 2 2c, which gives f(r) f((1 - m1/4)r)) 2 (rn 2 cnf(r).
1-e-arg
Lemma 32. For any 1,8 2 O, we have
f(t) f(s) S (rz aara V) 8) 1))t s/.
We record two further properties of f that will be used repeatedly. Lemma 33. The function satisfies:
(i) f is concave on r2) with f(0) = 0 and fr() c arr for T2. In particular, f S 1 on r2].
(ii) For any d' 2 1 and r 2 0 with d'r S r2, we have f(cr) - c((())
(iii) For any 0 with r S r2, we have ear f(r).
(iv) For any r > we have r S 4e s(r).
(v) For any o, 2 0, we have f(a + b) S f(a) + b e c-ar (ab + )
Proof. Part (i) follows directly from the definition. For (ii), since f is concave on [0, with f(0) = 0, the function r >> ((r)rr is non-increasing on r2]. Therefore f(dr)/ S f(r)/r, giving f(c'r) c' f(r). For (iii), we have f(r) (1 e arr) )/a 2 re re-ar, where the inequality follows from 1 e- 2 ue for all u 2 0.
For (iv), the case r r2 follows from (iii) since car S car For T2, the explicit form
gives f(r) 2 2ra
<022 (r2 - rτ). Since r > T2, we have 72 - 1 2 r²/2 whenever r 2 T2 72
giving f(r) 2 ear,, and hence r S f(r)) < 4e f(r) (using r I 2 fr)
4r2
which rearranges to (4ee (()) < 12/rr < 1). For 12 < < r2 we use s(r) 2 (12)
(1 eL-e-ara//e / 2 r2e are,,, so ne f(r) 2 4r2 2r.
For (v), we consider three cases. a b < r2: since f < 1 on [0, r2). we have f(a ((a+b)) f(a)= S b.
and the remaining terms are non-negative. If a > T2: both a and a a+ b lie in the quadratic region,
-22 so f(a + b) f(a)= 2r2 ((a + b)2 a²)
-= (2ab + 12)
ar2 (ab + 102), and b 2 r2
2r2
gives the result. If a S T2 a b. we split f(a + b) fa)
(f(rz)
s(a)) [ʃ(a + b)
2ra (2ab + t²), where in the last step we used (a + 6)2 - r3 < (a + b)² d² 2ab 62 (since r2).
f(22) < (r2 a) + 2r. ((a - r S b+ + -
-r2
We also recall the following result from Majka et al. 2020) (Lemma 4.4 therein), which captures the net contraction arising from the interplay between the drift expansion and the coupling.
Lemma 34. The function f, satisfies the property that for all i [0, ri),
1 25' c n + 7 min( /7, () /sup ("()) S (f).
31

---

Published as a conference paper at ICLR 2026
F.2.2 MOMENT BOUNDS ON
We now establish that B[aa]] exponentially small, a fact that will be used in several places. Recall that a a1 + where
4La |ell< ) +
02 = (T. 2
Z 1ltt 8
a1 =
and Ç is a one-dimensional standard Gaussian, Z' ~ (1,
For ai, set 3 4La /lll/ VT. Since Ç NO0, 1), a standard Gaussian tail bound yields
Bla $ X2 exp ( (-oy g2
X2
(47)
for some universal constant C1 O.
concentra- For ag, since || has the law of weighted x² random variable, the Laurent-Massart concenttt tion inequality gives
Ad Adeap (G) 21
E[aš] S X2 exp (-02
(48)
L 7
for some universal C2 >
Combining via Minkowski's inequality (ll 1lzz triangle inequality), we obtain
E|a2]1/2 s B[a + [[a3)
(49)
Assumption 26 ensures that the right-hand side is small needed. In particular, for the remainder of the proof, we assume that Assumption 26 is strong enough to guarantee:
E|a2]1//5 32
P(a > 0) - 128(c + Ls + 1)ea
(50)
Both conditions are satisfied since B[a and Ia(a > 0) are bounded by the exponentially small quantities in (47) (49,, and Assumption 26 controls these exponents.
We also establish the following fact, which underpins the small-R analysis.
Lemma 35. Suppose R TI. If 0. then f $ r1.
Proof. If R 4B/X, then by 46 with 0, we have (1+mLs) < (1 (1+1)o))-44)) = ri where the last equality is the definition of r1. If 4B/A & R < Fi, then by (45) with a 0, we have A < (1 7//4)R < <1.
We also establish a lower bound on that will be needed for the coupling analysis.
Lemma 36. Suppose R Tı and a O. Then, provided тo S 4[++1s+]] and 110 S 2 (both of which are ensured by the minimum in the definition of no), we have i 2 R/4.
Proof. By the reverse triangle inequality,
$ ス i (T)Ea(() Eg (tt.)llloo++(111111
2 i - RnLoR R la
On laz- -0}, have 12122 < 2/8. so VTLLllZ'lla S vn/8. Since 7h S 2/A, have
Vn2/85 $188$1/1 1/2.
Case R 4B/A: We F > (> LR)R (from the reverse triangle inequality applied P). Therefore.
i 2 (1 7 L.) - nA/8) R.
Since 71 < [[A+a++]) we have ( + Lb) < 1 Combined with V 1A/88 < 1/2, the coefficient is bounded below by 1/4 1 1/2 1/4.
32

---

Published as a conference paper at ICLR 2026
Case 4B/X < R m We use F 2 (1 1 n//2)R (from R 2 4B/X). Then,
A (1 1A/2 vn\/8)R.
Since 1J0 2/A, we have 1//2 and v/2/8 - 1/2. Furthermore, 7/ S T(- I
gives 1 Σ 1/8. Hence the coefficient at least 1/8 - 1/2 = 3/8 2 1/4. Both cases yield f 2 R/4.
F.2.3 COUPLING CONTRACTION
We now combine the ingredients above to prove the contraction
[[)())) S (1 nc/8) f(R),
(51)
splitting into two regimes.
Case 1: R 2 r1. Since ri 4B/A, Lemma 30 gives f m//4 + n,) Applying Lemma with 1. = f and 8 = (1 - -A/4) followed by Lemma 31
32
f(f) - (1 (c) + (rz ar (1+ 1
na R + 1) naR.
For R 2 r2, the explicit form of f on [r2, oo) gives R² S 4 (R) and R S re f(R), so
ar erz' + na) na S 8ear2(1 + na)na
For rı R r2, since f'is is concave on r2) with f(0) 0, Lemma 33
lii) gives R S eaR caRs(()
S
ears (R), from which the same bound follows (with a larger constant). In either case,
f(f) S (1 - (c) Se (1 + na)).
(52)
Taking expectations of the error term, we expand [E[( + ma)aa = E[a] + r)B[a²]. By Jensen's inequality, [[[o] < R[a2]1/2 Since 2//2]1/2 < 1 (which follows from (50) ). we also have E[a?) s
1/2 Therefore, [E( + na)aa I (1 + Ea & 2F/2211/2 giving
E[f(F)) s (1 - nc) ( + 16e F) ((R) & (1 - nc/))
where the last inequality the smallness condition (50)
remains to pass from to It' We claim that (('( z'| S f(f). Indeed:
If 2 F1. the synchronous coupling is used, so R' =.
If f < r1. then R' < f /To I and f is concave on [0. r2). Since ZRR = i
(Lemma 29h, Jensen's inequality gives [[ < (E z")) f(f).
Therefore, (R)( B[((S) < (1 (I).
Case 2: R < T1. We condition on the event {a = 0} and its complement.
we have 3 < r1. The coupling then yields H' < To 12, so both and R' lie in the concave region of f. We combine the drift and coupling contributions using Lemma 34
On fa = 0}: By Lemma 35
For the drift from R A: since f is concave on [0,
f(f) S(R) S (F) - R-)-
When f < R, have f(f) f(R) < 0 (f) When i > R. use i - R
NLLR < L26 (from 46 with 0) together with fR( =c S e-at ai af(- S f (f
to obtain
f(A) S(R) < (
anLor nb < (F)
(53)
For the coupling from f to Rt: since < 1ı and ET = f. the same Taylor expansion argument
in Section of Majka et al. 22200 gives
[[((R)) f(f) < sup ln(0) 0o min( (nn, oelr
(54)
33

---

Published as a conference paper at ICLR 2026
Adding (53) and 54 and applying Lemma
-] - - -enf(f).
(55)
By Lemma 36
R/4, so Lemma (3 (applied with c' 4 and R/4) gives f(R)
4f(R/4) 4 (f). Therefore,
= f(R) f()) o1 5(B) (1 (1- {(t).
(56)
On {a > 0}: We use crude bound. The coupling satisfies ([)( S {()) in all cases (either by concavity when S ri, or by the synchronous coupling when F r1). From Lemma 30 f f + R. Applying Lemma 33 v) with R. and b 7L4 + a.
E
1 La + R2
f(A) S f(R) + MLLB + +a)
T2
Since R <r r2, have p-are 12
R se S R, and the above simplifies to
f(f) S f(R) + 3n + aa
using ri La + a) < 1 (which holds with room to spare on the event {a > 0f for no sufficiently small, since a has sub-exponential tails). By Lemma (3)ii,, I ealt f(R) S e f(R). Therefore,
E[f(R')z') s (( + 3n(L, + a)ea))
(57)
Combining:
Taking total expectations.
[[f(R') = E 1
$ (1 en O)
(, + are [as0]
(1 4 f(R) ((R) ((R)
01
c7
P(a 0) 3ne t E (Ls + a)
4
We bound the bracketed term, which we denote by B. By
Schwarz.
Ejalasul $ E(2 > 0) 1/2 = 1/2
Therefore,
B ? P(a > + (L, > 0) + 122
We bound each term using the conditions in (50)
For the first term,
T P(a > 0)
P((a
Cn
4
4
128(c + Lo + 11eani
512"
For the second term,
Lc
3cm
am LoP(a 0) < am Lb
128(c + Ly + 1)
128(c Lu + 1)
128
where the last step uses c+L++I 1. For the third term,
3cm
" 2/22/1/2 < 3nears
32eara
32
using T2. Summing, we obtain
Cn
3crj
3cm
c71(1 + 12 + 48)
61en
CT
<
<
512
128
32
512
512
8
where the last inequality follows from 61 64. Therefore,
Elf(H') < (1 (1- en +
61en
c7 f(R) = (1 67cm f(R) - (1- f(R). ((R).
4
512
512
This gives the stated bound (51)
34

---

Published as a conference paper at ICLR 2026
F.2.4 EXTENSION TO THE MISMATCHED CASE
We now treat the general case o b, + ō. The coupling (X', Y") is constructed using the matched dynamics (b, Eg), so Y' has the law of a single step of the y-process with drift b and residual volatility Eg. To correct for the mismatch, we define
= Y' w
n(b(() b(y)) vn(E_o()) Eo(y))
Since w is a deterministic shift (given Z') that adjusts the mean and residual covariance, Y has the correct mismatched law: conditionally on Z', we have Y ~ Ni nG) where Im (1 11)y v)b()) + vi (v) Z
- -Yllo+ R' + |wlla+. Applying Lemma 33 v) with Rt' and By the triangle inequality, ||X'
b = w
f(I - Ylo- (( f(R)) + |wlla+ +///)l + e ears
R [ + ||00 (
r2
Taking the conditional expectation over and using that ||w]|a is Z'-measurable and E R 2z (Lemma 29
B Yllo)) ([f(R') z4) + || + e
wl + 2 [111ll ).
(58)
T2
Taking full expectations and using the matched contraction
(51)
E )} (1- f(8) +[)))
-ara
[B] 1 =Bllllla.... (59)
72
We now bound the cross term [[i]| Since S (1 + Ths + njx) R (Lemma 30
(1+ + |+ |
By Ceuchy-Schware Efa [[][[ < ||a2)1/22/ [|= )1/2_ Since E[a2]1/2 is exponentially small by (50) and 1/ S 172 for S r2 (with similar control for R > r2), this second term absorbed into the quadratic error. For the leading term, we use Lemma 33 iv): R < ar f(R) for all R > 0. Therefore.
4(1 TTT)))
_ f(R) B[[l[[c]]:
8
(1 (1+nyoo R1o <
f (R))B[lwlal))
r2
T2
ri
using 1 + 11 < 2 (for nhu sufficiently small) and r2 2 r1. By the assumption in 339 and the estimate |]][w||a+] < nB. + Mm B., we obtain,
[[ll+++ < - (1 + Vpm Br
8 crjy
CT)
256
32
Substituting back into 59 and collecting the quadratic error terms,
E|f(IIX' Ylaat) ) < (1 cn
17 ()) f(R) + 32
C. f E[[[wlo++] + -e e - E[]lw]]
T2
for some universal constant > 0, yielding a contraction at rate c/16, We also obtain the second moment bound,
E
w 1/ < 2 D3 + 2n 12 27y² D + 2mm
using Ell 1 = tr( < n. Substituting, we arrive at
|f[(I' - - ]
(7)
f(R) + nB + vnn B, +
3
arz (1² BJ + mmn²)
2r2
<
CT)
[(R) BBn vmm B.),
16
completing the proof of Proposition 27
35

---

Published as a conference paper at ICLR 2026
G PROOFS FOR THE STABILITY OF THE NOISY GRADIENT ESTIMATOR
For simplicity, we will prove for the case of fixed step-size ? = but is easily extendable to the full setting. Using the Wasserstein contraction obtained in the previous section, will now prove Proposition [4
Proposition 14
Consider the score matching algorithm Aum - s0 for some fixed K e N
where (0k) is as given in (16). Suppose that assumptions 10 12 and 13 hold, then there exists some nj > such that, if supx Th S T) we obtain that Asm score stable with constant
I'C N
(ePN&C'
(GPNDC2
min Σ Ti exp
stab
V
BA
exp
1mmin X2
(mminx X2
Imii
for some constant, 2(M, L, Agap: Amax (E), Be, C.) > O and 1mmin ming 7]--
The proof applies Proposition 27 to the process in (16). As in the proof of Proposition I1 we obtain stability estimates by analysing the trajectories 9g and 0₂ trained on S and SN with coupled minibatch indices.
Assign variables of Proposition 27 the dynamics of 0% and 0g take the form 35
Given set of minibatch indices B C [N] with |B] = NB.
36 with
b(0) = E Clipp (G(0, (sx))i 0, B.S.
b(0) = Clipo (G) (E:)EB) B,S
a(0)² = (0,
ā(0)² = Vم (0, B).
Using the shorthand us.y (0) = Wi(),)) Vollso (X.,)) t(as))
V log (t 10 (2,)) 1 12, we bound
the covariance from below:
Eg(0, B)
Σ
Cov iEB
> PNB Σ.
PN
Therefore a0)2 (0, N Σ, and similarly for ō. We set
G= 2PN E,
so that (0)2 > G. The residual volatilities are then
Eg (0) (o(0) =VT(ss(B))) 2PN Σ) 1/2
and similarly for Eg The weighted norm satisfies |10lc+ S 2PPNB/ B/ (mdepp) 11011
Lipschitz constants. Since Clipo is 1-Lipschitz, Jensen's 's inequality gives [lb(0) b(0') Ila I E [[] (zx)ee) - G(B', (0 B.S). Expanding via the chain rule,
U (0) v., = 2we,, (V Vasgr) T (sa - V log Pog + (Vosor) T(so ( - sw))
where all functions are evaluated at (Xa,,, tr,)). Taking G -ooms and applying the triangle inequality. the Lipschitz constant L(z, of so in 0 and the smoothness constant (z, from Assumption 13 give
I| (0) (0))la=
< 2wn., ((X),,,t t( v log Pt.., (X1,) xx)ll L(X t())) 110 l'llo.
Taking expectations and applying Cauchy-Schward each term, using TE < BlAAIsS
M. and Assumption 10 to bound Eillsp V log lrol/ z;1 /2 < Bz/o?.
((]))),(() (0'))la s
B: M oī
+ +I:) 00 - 0 l'
36

---

Published as a conference paper at ICLR 2026
Averaging over E B and [P] and integrating against using C, f 01 ",,0"1100 T(dt),
E (Ta))i))) G(Ø', (x))EB)))) 22MM,B.C. +T)1l - 8 ll
giving L, 2(M, LT). Furthermore, due gradient clipping, [[b()]la++ v
/2PNB/ B (1p) B.
For the Lipschitz constant of Eg. we first bound |E, (0( n||=s(0, B)
E0( B+ Writing (0) = Clipo(G(0, (ra)sen)) and μ((0) == b(0), we decompose
Es( Es(O', B) = E (Z()) - (M()))))))((()0))))))) (m(0))(0)) -
For the first term, we write 2(0)2(0)* Z(0") Z(0))= = =(((0) Z(0)) + 2(0')
Z (6"))" and use (() (| ((0) C by the definition of Clipe (1). giving
[())) - Z(6)) Z(0)) op.o. 2CC ll+)]
- For the second term, the same bound gives |m()))1000~000))00)))00000)))))) 1op. I Since Clipo is 1-L.ipshhtzz both terms are controlled by G()) G(f')lla-]- and combining.
[Es (0 B)llop.G- & &CLo 8 l'laa.
(60)
To convert this to Lipschitz bound on Eg itself, we apply the following matrix square root perturbation bound.
Lemma 37. Let A, B, We E R positive semi-definite mutrices such that A, B b and both A and B commute with W. 6 is the smallest nonzero eigenvalue of W, then
1|A B Bllop.
|| - B2/221 1lp.2 S
2v0
Proof. Because A and B commute with W. their principal square roots A1/2 and Y B1/2 also commute with W. and in particular X, Y and X - Y all preserve the range of W. Restricted to (W), the identity A- B X (X -Y)+(X - Y)Y is a Sylvester equation for X Y.. Because (ang satisfies UT Ae X Range(W) satisfy Amin (X (Range(W)) > V5, and similarly for Y. Applying the standard Sylvester equation operator norm bound via the integral representation Y Jo e X XA B)e ds restricted Range(W),
A > W, every unit vector
Ww 6, so the eigenvalues of
|| Bll Bloop.W
AA lon.w
X - Yloop. S
Amin [ange(W)) + Amin (Y Range(W))
2v8
To bypass the fact that Eg may not commute with G, we evaluate the matrix square root perturbation in the standard unweighted operator norm, which corresponds setting W scalar multiple of the identity matrix restricted to the range G.
- E, the minimum non-zero eigenvalue of (o) on Because (o) = (),s(0, - G > ZPNB
the range of Σ is strictly bounded below by 5 2PNE. Thus, setting W öî restricted to this range perfectly satisfies the conditions of Lemma matrix, trivially commutes with A and B.
ndve
37 A, B > W, and because W is scaled identity
From Lemma 37 we obtain,
El, B) Es(O', B) Ilp.ss
[[E((0) Ec(('))l <
2
To convert this back the G+-norm, we the standard norm equivalence || S Mx(G and ||M[ S /ICC)] where the condition number is x(G) Amaz
Ago
37

---

Published as a conference paper at ICLR 2026
Chaining these bounds with (60) yields:
Amax [2GG0) (0((0)llop.d))
||x((0) Ea(O')llop.a+ S
Agap
Amax (1) n||_s((, - Es( B)llop.hi
Agap
2
Amax (E) n1, B) Es(0', B) B)olop.c+
Agap
2
naw
2PN
Amax () (8C 10 - 0'|lc+)
2PN
Aap
T7
Amax (
2nPN
L
110
8'11G+
Apap
App
We therefore obtain L. = LC sufficiently small.
Amaz
2rPNe App
Since L3 x 7, Assumption 26 is satisfied for
Aep
Mismatch bounds. Since ã; = for all i N, the gradient estimators are identical whenever N & B, justifying the indicators 1 NeB throughout. On the event {N E B}, since Clipc is 1 1-Lipschitz and only the N-th summand of G differs between S and SN
2C
[bb(0) - b(0)1 < E[||G- Gs* [|] =
E 1l (8:2a) - IN (6: xN) S NB
NB
where ON (0;(0;)) 1 プド= "N (0; 2) and the last step uses 1 (0; I C (since each || ill bounded by C. before clipping, by Assumptions 10 and Amax (G+) /2PNB/ (1))p))
13
Converting to the -+-morm via
2P
[b6(0) b(0)lo+ s $20
LNEB Bo.
B
For the residual volatility mismatch, we use that =(2 Egn (0, B)). which is nonzero only when N E B. On this event, the same covariance perturbation argument as above, with Zg := Clipo(G(B,, (xx)xeB)) Zgv == Clipo (εx)ee)) and Zs | s C by definition of Clipc. gives
8C2
(s B) Esx (0, B)llop S 4C Zgx [1] S
1
NB
where the second step with = / (2PNB).
7Zg Zg S ||Gs GsN || s 2C//NB. Applying Lemma 37 before
70/2ss(0, B3 VsN (0, B)llop.w
Amax (G+)
- (0000 S
2/ 2//pp (2PNg)
Amin (G+)
4n02
12PN
Amax G++
< NB
INEB
Amin (G+)
PrpP
Amax
= 20
LNEB =; B
N NB
Agap
To verify the condition 33). we require BB
/nn B₁ S cnri empi/256. We have
P
Amaxe (2) Agap
nB mm Br < V NBA
20 402 Vmmi
1NEB.
38

---

Published as a conference paper at ICLR 2026
Meanwhile, 4(1 + MuLL)B/X ~ (NgD)())ap)) emr PNgC/>gap Since the
right-hand side grows PNB/>up Agap while the mismatch terms grow vn(P (N the condition is satisfied for 7) sufficiently small.
Applying Proposition 27
With assumptions verified, Proposition 27 gives, writing d(0,0) =
s(1 0-)))-), )
B[d(0x+1, (0(0+++,0+++)000.. Bz Bx] - (1
Ca d(Ox, 0k) 2(nBo + vnn Bo).
(61)
16
Since B, and Bg are nonzero only on {N E Bxh, which has probability NB/N, taking full expectations gives
20 N
2np
Amax (2))
Bl + vm B.) S
1 + 2C 21m0
NAAepp
Agap
where the second inequality holds for v sufficiently small. Taking full expectations in (61).
E|d(0x+1, 0 11+1)] - (1
Cy Bld(Ok. 0+]+
4C N
2np NgAuap
+ 20
2mm Vm
Amax Σ)
(62)
16
Agap
Iterating (62) with d(Oo, 0o) 0 (since 0g section,
0o) and applying the formula from the start of this
K-1
E[d(OK, [x(0x,0x)] <
4C
2np
1 2C 20
Amax
2
Cn
(63)
N
V NadSp
Agap
16
k=U
4C
2P
Amax (I)
16
- N
1
min
nk,,
(64)
NBm
Xgap
Since f(r)
"02 2 for all 2 have |10 1'1 S 2rzeara 2ra
d(A, 8'). Since Assumption 13
states that 8g is Lipschitz in with respect to || || and ||. 113
7
|| 112
2PN,
/ b sex (X,,t) 50 (Xe, 1)1 Xo ā, (dtt)
(65)
< [11[ - -0x1
nL²
E E[11xx - 0x+1a=]
2PN
nT. nT.2 azear
Eld(Ox, 0k)1
PNB
C
2m
Xmax (2) Axap
16 min nk,
< 4L T2e
ara
1 2C Vmm
-
(66)
N
PN3
A
B gap
When is sufficiently small, we obtain yo 2 A-1 and therefore
PN 2
L LS + T. 1/2 1
X2 X2
and since L and grow 7 0+, also have rāc 2 exp(-6Lri/co) Substituting these estimates into (66) and absorbing are into the S yields the bound in the statement Proposition
14
39


Reviews:
[{"review_id": "5296ZpabH2", "reviewer": "Reviewer_bmqA", "summary": "This paper studies generalisation in denoising diffusion models from an algorithm-dependent perspective. It introduces score stability, a measure of how sensitive score-matching algorithms are to changes in the training data, and uses it to derive generalisation gap bounds. The authors identify three sources of implicit regularisation: (i) early stopping in the denoising process, (ii) coarse discretisation in the sampling algorithm, and (iii) optimisation noise from SGD with gradient clipping and weight decay. They show that each of these sources improves generalisation in different ways, such as limiting memorisation, controlling discretisation error, or inducing contractive training dynamics. By grounding the analysis in algorithmic properties rather than model structure, the paper provides a unifying theoretical framework for understanding why diffusion models generalise despite high dimensionality. Overall, it links practical training heuristics to formal generalisation guarantees, revealing overlooked mechanisms that naturally regularise diffusion models.", "questions": "**Q1.** How can score stability be practically estimated for large-scale diffusion models?\n\n**Q2.** Can the framework be extended to conditional or latent diffusion models? \n\n**Q3.** How sensitive are the generalisation bounds to different optimisation schemes?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "**1. Introduction of an algorithm-dependent generalization framework:** The paper introduces score stability, a novel, algorithm-dependent framework for analyzing generalization in diffusion models. Unlike previous works that rely on algorithm-independent uniform convergence bounds, this approach explicitly quantifies how training and sampling algorithms affect generalization. By linking the sensitivity of the learned score function to single-sample perturbations with the expected generalization gap, the framework provides theoretical insight into why diffusion models generalize despite high dimensionality.  \n\n**2. Identification of multiple sources of implicit regularization:** The paper systematically analyses three distinct sources of implicit regularization in diffusion model training: denoising regularisation via early stopping, sampler regularisation through coarse discretization, and optimization regularization induced by stochastic gradient dynamics. By grounding these effects in the algorithmic framework of score stability, the work reveals mechanisms that were previously overlooked in the literature.", "weaknesses": "**1. Difficulty in computing score stability:** While the paper introduces the notion of score stability and shows its theoretical connection to generalization, it does not provide a clear or practical method for computing or estimating $\\epsilon_{stab}$ for a given algorithm and dataset. In practice, evaluating score stability requires measuring the sensitivity of the learned score function to changes in individual data points, which may involve retraining or coupling multiple stochastic outputs—a process that can be computationally expensive, especially for large datasets and high-dimensional diffusion models.\n\n**2. The loss function for diffusion models has no integration over time:** The analysis in the paper treats the diffusion model loss as time-agnostic, meaning it does not explicitly integrate or weigh contributions across different noise levels or time steps. In practice, the difficulty of learning the score function varies with time: early steps with high noise are often easier to model, whereas later steps with low noise are more sensitive and crucial for high-quality generation. Ignoring this temporal structure may oversimplify the analysis and obscure how generalization and stability depend on the progression of the diffusion process.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "A6id0esgQM", "reviewer": "Reviewer_HHPm", "summary": "Diffusion models generalize their training data, and how this happens must depend on some combination of (i) the training objective, (ii) optimization (e.g., using SGD), and (iii) the details of sampling. Prior work has attempted to bound generalization error given various assumptions (e.g., by assuming the 'manifold hypothesis'), but these bounds might neglect various details important in practice. The authors attempt to rectify this by introducing the notion of \"score stability\", and show that this notion allows them to derive interesting bounds on generalization error.", "questions": "1. Can the authors show via some numerical experiments that their bounds are tight or at least qualitatively interesting?\n\n2. Can the authors come up with a toy setting where things can be explicitly computed, which may help show that their derived bounds are interesting?", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is relatively clearly written and the math claims appear to be rigorous (although I did not check them in detail). The notion of score stability is interesting and appears to be a useful conceptual/technical contribution. In the main text, the authors strike a good balance between presenting their results and presenting the intuition for them, without overwhelming the reader with technical detail.", "weaknesses": "I understand that the point of the paper is to prove bounds in the tradition of recent formal (mathematical) work on diffusion models, but I am still left wondering about the extent to which these bounds are tight and/or interesting in practice. There are no experiments in the paper; would it be possible to conduct any experiments to show that the bounds are interesting, or that they are at least qualitatively consistent with what one sees in experiments? Relatedly, figures (either showing schematics, or experiments) would be helpful in getting the reader to follow the arguments. \n\nAlso relatedly, are there any interesting (even toy) settings where the reader can analytically see that these bounds are interesting, because various things can be explicitly computed?\n\nThere is some interesting work on how mislearning in diffusion models supports generalization. Three examples I can think of are the Kadkhodaie et al. 2024 geometry-adaptive harmonic representations paper (https://arxiv.org/abs/2310.02557), the Kamb and Ganguli 2025 combinatorial creativity paper (https://arxiv.org/abs/2412.20292), and the Vastola 2025 generalization through variance paper (https://arxiv.org/abs/2504.12532). The authors don't have to cite these if they don't find them interesting/appropriate; I just mean to say that there's other work (not in the same formal mathematical tradition the authors follow) on the basic question they're interested in. This work may be useful in framing their contribution.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xJ0C2Tev8e", "reviewer": "Reviewer_dqKw", "summary": "This work studies the implicit regularization of denoising score matching, brought by the early stopping, timestep discretization and optimization algorithm. The analysis is based on a novel score stability framework, which converts the generalization bound into the sum of empirical score matching error plus a stability term.", "questions": "1. Can you provide some simulation to validate the theory? For example, just considering learning the score of a simple Gaussian distribution? Can you show by early stopping or using a coarser discretization, the generalization (KL divergence) is improved? Currently, the bound is not  intuitive for interpretation without solid example.\n\n2. How can your theory explain diffusion model can generalize despite the empirical optimum can only memorize?\n\n3. Although the theoretical result is rigorous, in practice, early stopping and discretization are not the key factors of generalization. What do you think about this? What are other factors that enable generalization?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The proposed score stability framework is sound and novel.\n\n2. The analysis is comprehensive, covering a broad range of aspects, including discretization, early stopping and optimization algorithm.", "weaknesses": "1. No empirical verification of the theorem is provided. \n\n2. it is unclear to me how the results explain how diffusion model generalizes despite the empirical optimum can only memorize.\n\n3. Empirically, early stopping and discretization are no the key to generalization.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)