Volterra Generative Models Yusen Jia
The Hong Kong University of Science and Technology (Guangzhou)
Bingyan Han
arXiv:2606.18071v1 [cs.LG] 16 Jun 2026
The Hong Kong University of Science and Technology (Guangzhou)
Abstract Score-based diffusion models typically use Brownian perturbations, which provide tractable reverse-time dynamics but impose memoryless noising. We introduce Volterra generative models, a continuous-time score-based framework whose forward process injects path-dependent noise through fractional kernels. To handle the non-Markovian and non-semimartingale dynamics, we construct finite-dimensional Markovian lifts using Gaussian quadrature in both regimes and a hybrid finitedifference exponential approximation in the smooth regime. We prove squared error bounds, derive an augmented linear-Gaussian forward process, and show that the learning can remain data-dimensional by considering residual states and analytic auxiliary Gaussian scores. We also identify covariance and reverse-time degeneracies caused by shared Brownian factors and signed smooth-regime weights. The degeneracy motivates stabilized conditioning and, for stiff larger lifts, a Gaussianbridge reconstruction sampler. Experiments on MNIST and CIFAR-10 show that persistent fractional perturbations with small Markovian lifts can improve scorebased generation on MNIST and provide a promising extension to natural images, while the bridge sampler provides a stability mechanism for larger lifts.
1
Introduction
Diffusion generative models have become a central methodology for learning high-dimensional data distributions. Their basic principle is to construct a forward noising process that gradually transforms data into a tractable reference distribution, and then to learn the reverse transformation from noise back to data. Early diffusion probabilistic models realized this idea through discrete nonequilibrium noising chains (Sohl-Dickstein et al., 2015; Ho et al., 2020), while score-based models learned gradients of perturbed data densities and used them for sampling (Song & Ermon, 2019). The continuous-time formulation of Song et al. (2020) unifies these views through stochastic differential equations (SDEs): a prescribed forward diffusion transports the data distribution to noise, and the reverse-time SDE is obtained by replacing an unknown drift term with the time-dependent score. This framework provides a flexible interface between stochastic analysis, denoising score matching (Hyvärinen & Dayan, 2005; Vincent, 2011), numerical SDE solvers, and neural score networks. Subsequent work has further clarified the design space of continuous-time diffusion models, including choices of noise schedules and samplers (Karras et al., 2022). A complementary line of work improves the diffusion models through scheduled sampling, distillation, and efficient architectures; examples include multi-step denoising scheduled sampling (Ren et al., 2024), spatial fitting-error reduction distillation (Zhou et al., 2024), and physics-inspired generative models (Xu et al., 2023). 1
Despite their empirical success, standard diffusion models usually rely on Brownian-driven Markovian perturbations. This choice is analytically convenient: Brownian increments are independent, the forward process is Markovian, and the conditional transition laws are often explicit Gaussians. However, the same structure also restricts the temporal dependence of the noising mechanism. The perturbation at a given time depends on the past only through the current state, and the injected noise has no memory. This is potentially limiting when one wants the forward process to encode persistent correlations, rough fluctuations, or multi-scale memory. Recent work has therefore begun to explore alternatives to Brownian noising, including Lévy-driven score models (Yoon et al., 2023) and fractional-noise-based generative models (Nobis et al., 2024). These directions suggest that the stochastic process used for noising is not merely an implementation detail; it is a modeling choice that can affect sample quality and numerical behavior. This paper develops a Volterra-based generative model in which the Brownian perturbation is replaced by a path-dependent stochastic convolution. The forward process takes the form Z t
Xt = ρ(t)X0 +
G(t − s)b(s) ds +
0
Z t
G(t − s)Σ(s) dBs ,
0
where the deterministic attenuation ρ(t)X0 controls the decay of the initial signal and the kernel G determines the memory structure of the injected perturbation. Our main examples are fractional kernels 1 tH− 2 GH (t) = , H ∈ (0, 1). Γ(H + 12 ) Stochastic Volterra equations with such kernels have been studied extensively, especially in connection with rough volatility and Markovian approximation methods (Carmona & Coutin, 1998; Gatheral et al., 2018; Alfonsi & Kebaier, 2024; Bayer & Breneis, 2023). The main difficulty is that Volterra processes are non-Markovian and non-semimartingale in general. Consequently, the usual reverse-time SDE machinery for Markov diffusions cannot be applied directly. The core idea of this work is to approximate the Volterra kernel by a finite exponential mixture and to lift the resulting process to a finite-dimensional Markovian system. For H < 12 , the fractional kernel is completely monotone and admits a positive Laplace representation, so Gaussian quadrature yields a positive exponential approximation. For H > 12 , complete monotonicity fails. We instead use a representation involving te−γt , apply Gaussian quadrature to obtain an intermediate approximation, and then convert each basis function into exponentials through a hybrid finite-difference rule. This produces a signed exponential mixture but keeps all final mean-reversion rates nonnegative. The lifted model introduces new challenges that are specific to the generative setting. All Markovian factors are driven by the same Brownian motion, so the auxiliary covariance matrix can be ill-conditioned even when the exponential weights are positive. In the smooth regime H > 21 , the finite-difference construction creates nearby signed pairs of exponential rates, which generates additional near-null covariance directions. Moreover, the signed weights cancel at the origin, eliminating the direct diffusion loading of the primary state and making a naive reverse Euler sampler unstable. These effects are structural rather than cosmetic: they determine how the score should be parameterized and how reverse-time sampling should be performed. Our main contributions are as follows. 1. We develop Gaussian-quadrature-based Markovian approximations for fractional Volterra kernels in both Hurst regimes. For H < 12 , we use the Laplace representation of the frac2
tional kernel and obtain a positive finite-exponential approximation. For H > 21 , where the kernel is no longer completely monotone, we use the te−γt representation and introduce a hybrid finite-difference construction that yields a signed finite-exponential approximation with nonnegative rates. We also prove non-asymptotic error bounds under generic quadrature tuning parameters, making the approximation theory compatible with numerical choices. 2. We derive the augmented reverse-time dynamics and clarify the regime-dependent sampling behavior of the Markovian lift. The augmented forward process is linear Gaussian, which allows us to learn only a data-dimensional residual score and compute the auxiliary Gaussian score analytically. In the smooth regime, the primary coordinate has no direct Brownian forcing or direct score correction. For small lifts, this degeneracy is not numerically harmful because the auxiliary coupling remains non-stiff, and ordinary explicit reverse discretization is stable. For larger lifts, however, the finite-difference rates can span several orders of magnitude, making explicit discretization unstable. We therefore introduce a Gaussianbridge reconstruction sampler as a stabilizer for stiff smooth-regime lifts. 3. We provide numerical analysis and MNIST experiments demonstrating the practical effect of the proposed Volterra noising mechanism. First, we demonstrate how the Hurst parameter and the lift size affect the final rates and weights, covariance singularity, and reverse-time stability. On MNIST, our method achieves its best performance with FID 0.52 when H = 0.9 and N = 2, outperforming the Brownian SDE baseline, GFDM, and representative MNIST benchmarks. We further provide a preliminary CIFAR-10 experiment under the same smalllift persistent regime, obtaining FID around 9.5, which suggests that the Volterra noising mechanism extends beyond grayscale digit generation. The closest works to ours are Yoon et al. (2023) and Nobis et al. (2024), which also modify the driving noise in score-based generative models. Yoon et al. (2023) replace Brownian motion by an isotropic α-stable Lévy process, thereby introducing heavy-tailed jumps while preserving independent increments. By contrast, our Volterra perturbation is Gaussian conditional on the data but has memory through a fractional convolution kernel. Hence, the main difficulty is non-Markovian temporal dependence rather than jump-driven heavy tails. Nobis et al. (2024) is more closely related, since they also use a finite-dimensional Markovian approximation of fractional Brownian noise and an augmented score decomposition. However, the approximation and sampling mechanisms are different. Nobis et al. (2024) approximate fractional Brownian motion by a linear combination of Ornstein–Uhlenbeck (OU) processes with geometrically spaced rates and coefficients chosen through an L2 (P) projection criterion. In contrast, we approximate the Volterra kernel itself by Gaussian quadrature, which gives explicit error bounds. The smooth regime also suffers from the zero diffusion issue. For small lifts, this can still be handled by explicit reverse discretization. For larger smooth-regime lifts, the finite-difference rates become stiff, and we introduce a Gaussian-bridge reconstruction sampler as a stabilization mechanism. The rest of the paper is organized as follows. Section 2 reviews continuous-time score-based generative modeling and introduces the stochastic Volterra forward process. Section 3 develops the Markovian approximation theory and the Gaussian quadrature constructions for H < 12 and H > 12 . Section 4 derives the augmented forward process, the augmented score-matching objective, and the reverse-time samplers. Section 5 presents the numerical implementation, MNIST generation results, and a preliminary CIFAR-10 study. The appendices contain the proofs and extra details on the algorithm design. 3
2
Formulation
This section introduces the continuous-time forward process used throughout the paper. We first recall the standard score-based diffusion formulation and fix our reverse-time convention. Then, we introduce the stochastic Volterra process, which replaces Brownian perturbations by pathdependent noise with fractional kernels. 2.1
Score-based Generative Models
Score-based generative modeling constructs a continuous path of probability distributions that transports the data distribution to a tractable reference distribution. In its standard continuoustime formulation, the forward noising process is an RD -valued Itô diffusion X0 ∼ p0 ,
dXt = f (Xt , t) dt + g(t) dBt ,
t ∈ [0, T ],
(2.1)
where B is a standard D-dimensional Brownian motion, f : RD × [0, T ] → RD is the drift, and g : [0, T ] → R+ is a scalar diffusion coefficient. We denote by pt the density of Xt . Throughout the paper, reverse time is parametrized by an increasing clock. More precisely, we define Xt := XT −t , pt := pT −t , g(t) := g(T − t), f (x, t) := −f (x, T − t), (2.2) for t ∈ [0, T ]. Thus X0 = XT is the terminal noised variable and XT = X0 has the data distribution. In this convention, reverse time evolves with dt > 0. Under standard regularity assumptions, the time reversal of (2.1) is again an Itô diffusion (Stratonovich, 1965; Anderson, 1982; Föllmer, 2006). With the convention (2.2), it is given by h
i
dXt = f (Xt , t) + g 2 (t) ∇x log pt (Xt ) dt + g(t) dBt ,
t ∈ [0, T ],
(2.3)
where B is a Brownian motion in the reversed filtration. Equivalently, ∇x log pt (x) = ∇x log pT −t (x). The unknown object in (2.3) is the time-dependent score function ∇x log pt (x). Learning this score is the central problem in score-based generative modeling. The reason score matching is tractable is that one usually chooses the forward process so that the conditional law of Xt given X0 is explicit, typically Gaussian. Let p0t (· | x0 ) denote this conditional density. Under mild regularity conditions, the marginal score can be expressed as a conditional expectation of the tractable conditional score. Indeed, since Z
pt (x) =
p0t (x | x0 )p0 (x0 ) dx0 ,
differentiating under the integral gives ∇x log pt (x) = E [∇x log p0t (x | X0 ) | Xt = x] .
(2.4)
This identity is the basis of denoising score matching (Song et al., 2020). It shows that the marginal score, which is unknown because it depends on the data distribution, can be learned by regressing against the conditional score, which is explicit for a suitably chosen forward process. 4
Accordingly, we train a time-dependent score model sθ : RD × [0, T ] → RD by minimizing h
i
L(θ) = Et EX0 EXt |X0 λ(t) ∥sθ (Xt , t) − ∇x log p0t (Xt | X0 )∥22 ,
(2.5)
where λ(t) > 0 is a weighting function. Here t is sampled from a prescribed distribution on [0, T ], usually the uniform distribution, X0 ∼ p0 is a data sample, and Xt | X0 ∼ p0t (· | X0 ) is generated by the forward noising process. For each fixed t, the L2 -optimal predictor in (2.5) is E [∇x log p0t (Xt | X0 ) | Xt = x] , and hence equals ∇x log pt (x) by (2.4). Once trained, sθ is substituted for the score in (2.3), and samples are generated by evolving the reverse dynamics from an approximate sample of pT back to time zero. Classical diffusion models rely on the Markovian and Brownian structure of (2.1). This structure is analytically convenient, but it also imposes a memoryless noising mechanism: increments of the driving noise are independent, and the perturbation at time t depends on the past only through the current state. This is restrictive when the desired forward perturbation is intended to encode temporal dependence, multi-scale memory, or fractional-type roughness. This motivates replacing Brownian noising by alternative choices. In this paper, we investigate a Volterra perturbation approach. 2.2
Stochastic Volterra Equations and Fractional Kernels
We now introduce the Volterra forward model used in this work. Let ρ(t) := e−κi⋆ t with constant κi⋆ > 0 chosen later. We consider the state-independent stochastic Volterra process Z t
Xt = ρ(t)X0 +
G(t − s)b(s) ds +
Z t
G(t − s)Σ(s) dBs ,
t ∈ [0, T ],
(2.6)
0
0
where X0 ∼ p0 , G : R+ → R ∪ {+∞} is a convolution kernel, b : [0, T ] → RD is a deterministic drift schedule, and Σ : [0, T ] → RD×D is a deterministic volatility schedule. The function ρ acts as a deterministic signal attenuation schedule that controls the decay of the initial state, whereas the Volterra kernel G controls the temporal dependence of the injected perturbation. In applications, these schedules are chosen so that the terminal distribution is close to a tractable reference distribution. Conditional on X0 , the process in (2.6) is Gaussian. Indeed, Xt | X0 ∼ N (ρ(t)X0 + mG (t), CG (t)) , where
Z t
mG (t) :=
Z t
G(t − s)b(s) ds,
CG (t) :=
0
G2 (t − s)Σ(s)Σ(s)⊤ ds.
(2.7)
(2.8)
0
Thus the conditional score appearing in the denoising objective remains explicit, as in the Markovian diffusion setting. The key difference is that the covariance is generated by a path-dependent convolution against the driving Brownian motion, rather than by instantaneous Brownian noise. The main kernels considered in this paper are fractional kernels. A canonical example is 1
GH (t) =
tH− 2
Γ H + 12
,
5
H ∈ (0, 1),
(2.9)
which is also related to Riemann–Liouville fractional Brownian motion (Mandelbrot & Van Ness, 1968). The Hurst parameter H controls the regularity of the perturbation. The case H = 12 corresponds to the Brownian scaling. When H < 21 , the kernel is singular near the origin and produces rougher perturbations. When H > 12 , the kernel is smoother near the origin and produces more persistent fractional-type dependence. This provides a principled way to interpolate between rough and smooth noising mechanisms. The additional modeling flexibility comes at a cost. The process (2.6) is generally non-Markovian and non-semimartingale, since the state at time t depends on the full past of the driving noise through the convolution kernel G. Consequently, the reverse-time SDE theory used for (2.1) does not apply directly to (2.6). The central approximation step in this paper is therefore to replace the Volterra kernel by a finite-dimensional exponential approximation. This produces a Markovian lift of the Volterra dynamics in an augmented state space.
3
Markovian Approximation
This section constructs finite-dimensional Markovian approximations of the Volterra forward process (2.6). The deterministic signal attenuation term ρ(t)X0 is kept unchanged; only the convolution kernel is approximated. We work on a fixed time horizon [0, T ] and assume that b ∈ L∞ ([0, T ]; RD ), Σ ∈ L∞ ([0, T ]; RD×D ), and the kernels under consideration belong to L2 ([0, T ]). 3.1
Stochastic Volterra Equations and Markovian Lifts
A classical route to finite-dimensional approximations of Volterra processes is through Markovian lifts, going back to Carmona & Coutin (1998). The construction starts from an integral representation of the kernel, Z ∞
ϕ(t, γ) µ(dγ),
G(t) =
t > 0,
(3.1)
0
where ϕ(t, γ) is a parametric family of functions whose finite linear combinations admit a Markovian realization. If µ is signed, the corresponding total-variation integral is required to be finite on the time interval under consideration. Discretizing the γ-integral in (3.1) then yields a finitedimensional approximation of the Volterra kernel. The canonical choice is ϕ(t, γ) = e−γt . By Bernstein’s theorem (Schilling et al., 2012, Theorem 1.4), G is completely monotone on (0, ∞) if and only if there exists a nonnegative Borel measure µ on [0, ∞) such that Z e−γt µ(dγ),
G(t) =
t > 0,
[0,∞)
with [0,∞) e−γt µ(dγ) < ∞ for all t > 0. Here, complete monotonicity means that G ∈ C ∞ (0, ∞) and (−1)n G(n) (t) ≥ 0 for all n ≥ 0 and t > 0. In this case, finite quadrature approximations of µ give finite exponential approximations of the kernel. The following proposition records the corresponding Markovian lift in the additive, state-independent setting of (2.6). It is the specialization of the Markovian-lift construction used for stochastic Volterra equations in Alfonsi & Kebaier (2024, Proposition 2.1) and Bayer & Breneis (2023, Proposition 1.1). R
Proposition 3.1 (Finite-dimensional Markovian lift). Let I be a finite index set and let b G(t) =
X
ωi e−λi t ,
λi ∈ R,
i∈I
6
ωi ∈ R.
(3.2)
For each i ∈ I, let Yi solve dYti = −λi Yti dt + b(t) dt + Σ(t) dBt ,
Y0i = 0D .
(3.3)
Then the process b t = ρ(t)X0 + X
X
ωi Yti
(3.4)
i∈I
satisfies Z t b t = ρ(t)X0 + X
b − s)b(s) ds + G(t
Z t
b − s)Σ(s) dBs . G(t
(3.5)
0
0
The proof follows by writing the integral form Yti =
Z t
e−λi (t−s) b(s) ds +
Z t
e−λi (t−s) Σ(s) dBs
0
0
and substituting it into (3.4). If all λi ≥ 0, the factors are OU-type mean-reverting factors. If some λi < 0, the factors remain well defined on the finite horizon [0, T ], but they are no longer mean reverting. b in The approximation error is controlled by the L2 error of the kernels. With X in (2.6) and X (3.5), for every t ∈ [0, T ], we have bt Xt − X
E
2
≤ CT
Z t
b G(s) − G(s)
2
ds,
(3.6)
0
where one may take CT = 2T ∥b∥2L∞ (0,T ) + 2∥Σ∥2L∞ (0,T ) . This follows from Cauchy’s inequality for the drift term and Itô’s isometry for the stochastic integral. Hence the construction of a Markovian approximation reduces to the construction of an accurate finite-exponential approximation of the kernel. 3.2
Gaussian Quadrature for H ∈ (0, 21 )
We first consider the rough regime. For H ∈ (0, 12 ), the normalized fractional kernel (2.9) admits the positive Laplace representation Z ∞
GH (t) = cH
1
e−γt γ −H− 2 dγ,
cH :=
0
1
. Γ(H + 12 )Γ( 21 − H)
(3.7)
This representation allows us to apply Gaussian quadrature with respect to the positive weight 1 wH (γ) = cH γ −H− 2 . Let N ∈ N denote the quadrature budget and let α, β, a, b > 0 be tuning parameters. Set A− :=
1 1 + 3 H 2 −H
!1/2
.
Choose integers m = mN and n = nN such that m≈
β √ N, A−
n≈
A− √ N, β
7
mn ≈ N,
where ≈ denotes the rounding convention used to make m and n integer-valued. Here, m is the level of the Gaussian quadrature rule applied to each individual subinterval. This means there are m nodes generated for each slice of the domain. Define √ α ξ0 := a exp − 3 N , ( 2 − H)A− !
ξn := b exp
α √ N , HA−
and
ξn j/n , j = 0, . . . , n. ξj := ξ0 ξ0 The quantities ξj are interval endpoints. On each interval [ξj , ξj+1 ], j = 0, . . . , n − 1, apply an m-point Gaussian quadrature rule with respect to wH . This produces local nodes and weights, which we relabel globally as (γi , ωi )mn i=1 . Finally, add the zero node
Z ξ0
γ0 := 0,
ω0 := cH 0
The resulting approximation is b H (t) = G
1
1 cH −H ξ02 . − H 2
γ −H− 2 dγ = 1
mn X
ωi e−γi t .
(3.8)
i=0
b H is completely monotone and can be inserted directly into All weights are nonnegative, so G Proposition 3.1 with nonnegative rates.
The following theorem gives a non-asymptotic error estimate for generic tuning parameters. Define the geometric ratio 1/n ξn . r−,N := ξ0 Theorem 3.2 (Generic rough-regime quadrature bound). Let H ∈ (0, 21 ), let X solve (2.6) with b be defined by (3.5) with G b=G b H . Then G = GH , and let X
E
bT XT − X
2
≤ CT c2H
2α √ T3 3−2H a exp − N A− ( 32 − H)2 3 −2H 2α √ + b exp − N 2H 2 A− 3 2H 2 5π T n r−,N − 1 4m+2 + . 12 H m2H 2
(3.9)
The bound is valid for arbitrary positive α, β, a, b. To turn it into a convergence statement, the middle term must also decay. For fixed a, b and the unrounded choices of m and n, one has r−,N → eαβ ; hence a simple sufficient condition for decay of the displayed Peano-kernel bound is αβ < log 3. The commonly used calibrated values of α, β, a, b may perform better in practice than this sufficient condition suggests, but then the above bound should be interpreted as a nonasymptotic upper bound rather than as an optimized convergence estimate. 3.3
Markovian Approximation for H ∈ ( 21 , 1)
For H ∈ ( 12 , 1), the normalized fractional kernel (2.9) is no longer completely monotone. Instead, for all t > 0, Z ∞ 1 1 GH (t) = cH te−γt γ 2 −H dγ, cH := . (3.10) 1 3 Γ(H + )Γ( 0 2 2 − H) 8
The additional factor t prevents a direct application of Proposition 3.1. We therefore first approximate the integral in (3.10) by Gaussian quadrature and then replace each basis function te−γt by a finite difference of exponentials. Set A+ :=
1 1 + 5 H 2 −H
!1/2
.
Let N ∈ N denote the quadrature budget, and let α, β, a, b > 0 be tuning parameters. Choose integers m = mN and n = nN such that m≈
β √ N, A+
n≈
A+ √ N, β
mn ≈ N,
where ≈ denotes the rounding convention used to make m and n integer-valued. Define √ α ξ0 := a exp − 5 N , ( 2 − H)A+ !
and
ξj := ξ0
ξn ξ0
ξn := b exp
α √ N , HA+
j/n
,
j = 0, . . . , n.
On each interval [ξj , ξj+1 ], j = 0, . . . , n−1, apply an m-point Gaussian quadrature rule with respect to the positive weight 1 wH (γ) = cH γ 2 −H . After relabelling the resulting mn positive nodes and weights as (γi , ωi )mn i=1 , add the zero node Z ξ0
γ0 := 0,
ω0 := cH 0
1
3 cH −H ξ02 . 2 −H
γ 2 −H dγ = 3
This gives the intermediate approximation GH (t) :=
mn X
ωi te−γi t .
(3.11)
i=0
We now convert (3.11) into a finite sum of exponentials. The zero node is treated by a one-sided difference, while the positive nodes are treated by a two-sided difference. Fix δ1 , δ2 > 0 and assume δ2 ≤ ξ0 .
(3.12)
Since the positive Gaussian nodes satisfy γi ≥ ξ0 for i ≥ 1, condition (3.12) ensures that all left-shifted rates γi − δ2 are nonnegative. Define e H (t) := ω0 G
mn 1 − e−δ1 t X e−(γi −δ2 )t − e−(γi +δ1 )t + ωi . δ1 δ1 + δ 2 i=1
e H is a finite linear combination of exponentials: Then G e H (t) = G
mn mn X ω0 −0·t ω0 −δ1 t X ωi ωi e − e + e−(γi −δ2 )t − e−(γi +δ1 )t . δ1 δ1 δ + δ δ + δ 1 2 1 2 i=1 i=1
9
(3.13)
Thus the approximation can be inserted into Proposition 3.1. The resulting exponential weights are signed, but all rates are nonnegative under (3.12). Since ξ0 decreases with N , δ2 needs to be small to satisfy (3.12). b denote the Volterra approximation obtained by replacing G by G e H in (2.6): Let X Z t b t = ρ(t)X0 + X
e H (t − s)b(s) ds + G
Z t
e H (t − s)Σ(s) dBs . G
(3.14)
0
0
Define the geometric ratio
r+,N := and set WN+ :=
ξn ξ0
mn X
cH ωi = 3 2 −H i=1
1/n
,
3 2
ξn
−H
3 −H 2
− ξ0
.
The equality follows because Gaussian quadrature integrates constants exactly on each middle interval. Theorem 3.3 (Generic smooth-regime quadrature bound). Let H ∈ ( 12 , 1), let X solve (2.6) with b be defined by (3.14). Suppose that the quadrature construction above is used G = GH , and let X with arbitrary α, β, a, b > 0 and that (3.12) holds. Define 3T 5 3 −2H 2α √ 2α √ + 2 5−2H QN,T := cH a exp − N + b exp − N A+ 4H 2 A+ 5( 25 − H)2 +
49π 3 T 2H 2 2−2H n m 24 H
and DN :=
r+,N − 1 2
4m+2
,
(3.15)
ω0 δ1 + WN+ max{δ1 , δ2 }. 2
(3.16)
Then the kernel approximation satisfies Z T
e H (t) GH (t) − G
0
2
dt ≤ 2Q+ N,T +
2T 5 2 DN . 5
(3.17)
Consequently,
E
bT XT − X
2
≤ CT
2T 2Q+ N,T + 5
5
! 2 DN
,
(3.18)
where CT is the constant in (3.6). √ Since ξn grows exponentially in N , WN+ may grow rapidly. Therefore, the convergence of the second term in DN needs δ1 and δ2 to decrease properly with N , such that DN → 0. The proof is given in Section B.2. The first term in (3.17) is the Gaussian quadrature error, and the second is the additional finite-difference error. The long Peano-kernel argument is needed only for the explicit expression of Q+ N,T ; the finite-difference part is elementary. If a, b are fixed and the unrounded choices of m, n are used, then r+,N → eαβ . Thus, conditions such as αβ < log 3 give simple sufficient conditions for the displayed Peano-kernel upper bound to decay, although empirically calibrated parameters may still perform well outside this sufficient regime. The Markovian approximation consists of coupled OU factors driven by common Brownian noise and therefore exhibits strong cross-factor dependence. It is an intrinsic consequence of the Volterra representation and plays a central role in both the augmented score-matching construction and the reverse-time sampling procedure developed later. 10
4
Generative Modeling
We build a score-based generative model from the finite-dimensional Markovian approximations developed in Section 3. To avoid confusion with the quadrature nodes and weights used in Section 3, we use a separate notation for the final exponential kernel entering the generative model. After all quadrature and, when H > 12 , finite-difference steps have been performed, write the approximation as X b G(t) = ψi e−κi t , κi ≥ 0, ψi ∈ R, (4.1) i∈I
b b H (t) in (3.8) with H ∈ (0, 1/2) and for where I is a finite index set. The notation G(t) stands for G e H (t) in (3.13) with H ∈ (1/2, 1). In the rough regime, the notation coincides with the quadrature G output: κi = γi and ψi = ωi . In the smooth regime, (κi , ψi )i∈I denotes the relabelled exponential terms obtained from the hybrid finite-difference approximation (3.13); these are not the original Gaussian quadrature nodes and weights. For example, the one-sided zero-node term contributes rates 0 and δ1 , while each positive quadrature node contributes the rates γi − δ2 and γi + δ1 .
4.1
Augmented Forward Processes
We specialize the Volterra forward process to the isotropic noising setting b(t) = µ(t)1D ,
Σ(t) = g(t)ID ,
where µ and g are deterministic scalar functions. Hence, all Markovian factors share the same D-dimensional Brownian motion. b In the remainder of the paper, the generative In Section 3, the approximation is denoted by X. model is defined using this Markovian approximation. To avoid excessive notation, we drop the hat and write X for the approximating forward process:
Xt = ρ(t)X0 +
X
ψi Yti ,
(4.2)
i∈I
where Yti satisfies h
i
dYti = −κi Yti + µ(t)1D dt + g(t) dBt ,
Y0i = 0D .
(4.3)
Recall that the attenuation schedule is specified as ρ(t) := e−κi⋆ t . The index i⋆ ∈ I is chosen such that the corresponding ψi⋆ ̸= 0 and κi⋆ > 0. Hence, any rate-zero term in (4.1), when present, is retained in the auxiliary state. The main purpose of considering the index i⋆ is to exclude Yti⋆ and keep the primary state Xt random. Indeed, it is direct to show that the primary process (4.2) satisfies
dXt = −κi⋆ Xt −
X
(κj − κi⋆ )ψj Ytj + µ(t)ψ̄ 1D dt + g(t)ψ̄ dBt ,
j∈J
where J := I \ {i⋆ },
ψ̄ :=
X i∈I
11
ψi .
X0 ∼ p0 ,
(4.4)
In contrast, including all factors Yti would make Xt deterministic if X0 and all factors are known. For later use, fix an ordering J = {j1 , . . . , jM }, where M = |J |, and define the stacked process
j1 jM jM j1 Zt := Xt,1 , Yt,1 , . . . , Yt,1 , . . . , Yt,D ∈ RD(M +1) . , . . . , Xt,D , Yt,D
Then Zt solves the linear SDE dZt = FZt dt + bZ (t) dt + LZ (t) dBt .
(4.5)
Here, the coefficients are
F = ID ⊗ A, with
⊤
q := ψ̄, 1, . . . , 1 and
bZ (t) = 1D ⊗ µ(t)q ,
−κi⋆ 0 A= .. .
−(κj1 − κi⋆ )ψj1 −κj1
0
LZ (t) = ID ⊗ g(t)q ,
∈ RM +1
· · · −(κjM − κi⋆ )ψjM 0 . .. .. . . −κjM
0
The matrix LZ (t) has size D(M +1)×D. Hence the augmented diffusion is rank-deficient whenever M ≥ 1, reflecting the fact that all factors for a given data coordinate are driven by the same Brownian component. For later use, we present the distributions of the primal variable Xt and the augmented state Zt . For each i ∈ I, define the scalar factor mean Z t
ui (t) :=
e−κi (t−s) µ(s) ds.
(4.6)
0
Then set νx (t) :=
X
ψi ui (t),
φ(t, s) :=
X
ψi e−κi (t−s) g(s).
(4.7)
i∈I
i∈I
Proposition 4.1 shows another representation of the primal state Xt . The proof is straightforward and thus omitted. Indeed, (4.8) follows by substituting the integral forms of the factors and exchanging the finite sum with the integrals. The conditional Gaussian law (4.9) then follows from Itô’s isometry. Proposition 4.1. Under (4.3)–(4.4), the primary state Xt admits the representation Z t
Xt = ρ(t)X0 + νx (t)1D +
φ(t, s) dBs .
(4.8)
0
Consequently,
σx2 (t) :=
Xt | X0 ∼ N ρ(t)X0 + νx (t)1D , σx2 (t)ID ,
Z t
φ(t, s)2 ds.
0
For the augmented state, we can solve the linear SDE (4.5) as Zt = eFt Z0 +
Z t
eF(t−s) bZ (s) ds +
0
Z t 0
12
eF(t−s) LZ (s) dBs .
(4.9)
Conditional on X0 , the first two terms are deterministic and the last term is an Itô integral with deterministic integrand. Hence, Zt | X0 is Gaussian. Its conditional mean and covariance are denoted by mZ (t; X0 ) := E[Zt | X0 ], CZ (t) := Cov(Zt | X0 ). They solve d mZ (t; X0 ) = FmZ (t; X0 ) + bZ (t), dt d CZ (t) = FCZ (t) + CZ (t)F⊤ + LZ (t)LZ (t)⊤ , dt
(4.10) (4.11)
with initial conditions mZ (0; X0 ) = (X0 , 0, . . . , 0) under the chosen stacking and CZ (0) = 0. These equations follow from taking conditional expectations and differentiating the covariance. Next, we present the conditional law of the primary state Xt given X0 , YtJ , where
YtJ := Ytj1 , . . . , YtjM ∈ RDM . For a single coordinate, define the M × M covariance matrix and the cross-covariance vector Z t
Cyy (t)ℓk :=
e−κjℓ (t−s) e−κjk (t−s) g(s)2 ds,
(4.12)
e−κjℓ (t−s) g(s)φ(t, s) ds,
(4.13)
0
Z t
cyx (t)ℓ := 0
for ℓ, k = 1, . . . , M . The explicit covariance entries (4.12) and (4.13) follow from the integral representations of the factors and from (4.8). Also, define uJ (t) := (uj1 (t), . . . , ujM (t))⊤ .
cxx (t) := σx2 (t), When Cyy (t) is nonsingular, define η(t) := Cyy (t)−1 cyx (t),
cx|y (t) := cxx (t) − cyx (t)⊤ Cyy (t)−1 cyx (t).
If Cyy (t) is singular, the exact degenerate Gaussian conditioning formula can be written with the Moore–Penrose pseudoinverse Cyy (t)† . In implementation, a truncated inverse may be used for numerical stability. See Section 5.2 for more details. The well-known Gaussian conditioning formula gives Xt
X0 , YtJ ∼ N mx|y (t), cx|y (t)ID ,
(4.14)
where mx|y (t) = ρ(t)X0 + νx (t)1D +
M X
ηℓ (t) Ytjℓ − ujℓ (t)1D .
(4.15)
ℓ=1
We assume cx|y (t) > 0 on the training time interval. This excludes the degenerate case in which the omitted anchor factor is determined by the retained factors at time t. 13
4.2
Augmented Score Matching
The augmented forward process Zt is Markovian, but its state dimension D(M + 1) is much larger than the data dimension D. Training a full score network on the entire augmented state would therefore be unnecessarily expensive. The conditional Gaussian structure in Section 4.1 allows us to reduce the learned part of the score to a data-dimensional network, while keeping the auxiliary contribution analytic. Define the residual variable ξt := Xt −
M X
(4.16)
(4.17)
ηℓ (t) Ytjℓ − ujℓ (t)1D .
ℓ=1
By (4.14), conditional on X0 and YtJ , we have
ξt ∼ N ρ(t)X0 + νx (t)1D , cx|y (t)ID . Hence, the conditional score with respect to the primary coordinate is ∇x log p0t Xt YtJ , X0 = −
ξt − ρ(t)X0 − νx (t)1D . cx|y (t)
(4.18)
Then we train a time-dependent data-dimensional score model sθ : RD × [0, T ] → RD on the residual variable ξt . The augmented denoising score-matching objective is
2
ξt − ρ(t)X0 − νx (t)1D Laug (θ) := Et E(X0 ,YJ )∼p0 ⊗qt EXt |X0 ,YJ λ(t) sθ (ξt , t) + . t t cx|y (t) 2
(4.19)
Here, qt denotes the forward-time marginal density of YtJ . Time t is sampled from the training time distribution and λ(t) > 0 is the loss weight. For each fixed t, the L2 -optimal predictor is "
#
ξt − ρ(t)X0 − νx (t)1D s∗ (ξ, t) = E − ξt = ξ , cx|y (t) which equals the marginal score ∇ξ log pξt (ξ) of the residual variable ξt . Proposition 4.2 shows the advantage of considering score matching for ξt . The decomposition (4.22) is computationally important. The network sθ is only data-dimensional. The remaining auxiliary score is available analytically from the Gaussian covariance matrix. Indeed, with the stacking convention used in Section 4.1,
∇y log qt (YtJ ) = − Cyy (t)† ⊗ ID
h
i
YtJ − uJ (t) ⊗ 1D ,
(4.20)
where Cyy (t)† denotes the Moore–Penrose pseudoinverse to cover the singular case. Proposition 4.2 (Augmented score decomposition). Assume that sθ is optimal for (4.19). Define the learned component of the augmented score by Sθ (Zt , t) := (sθ (ξt , t), −η1 (t)sθ (ξt , t), . . . , −ηM (t)sθ (ξt , t)) , 14
(4.21)
where the tuple is understood blockwise as the primary component followed by the M auxiliary vector components. Then
J ∇\ z log pt (Zt ) := Sθ (Zt , t) + 0D , ∇y log qt (Yt )
(4.22)
is the corresponding augmented score estimator. If the exact inverse in the Gaussian conditioning formula is used and the score model is exact, then (4.22) equals the score ∇z log pt (Zt ). 4.3
Reverse-Time Dynamics
Since the augmented process Z = (Zt )t∈[0,T ] is Markovian, its reverse dynamics are obtained from the standard reverse-SDE formula. We use the reverse-time convention from Section 2. Then, Zt := ZT −t ,
F := −F,
bZ (t) := −bZ (T − t),
LZ (t) := LZ (T − t),
where t ∈ [0, T ] denotes the increasing reverse-time clock. The reverse SDE is h
i
dZt = F Zt + bZ (t) + LZ (t)LZ (t)⊤ ∇z log pT −t (Zt ) dt + LZ (t) dBt .
(4.23)
The unknown score in (4.23) is replaced by the augmented score estimator in Proposition 4.2. This is the continuous-time target of the sampling procedure in both Hurst regimes. The finite-dimensional lift has different numerical structures in the two regimes. In the rough regime, the final weights are nonnegative and the primary coordinate has a nonzero diffusion loading. In the smooth regime, the hybrid finite-difference approximation gives a signed exponential P sum satisfying ψ̄ = i∈I ψi = 0. The primary coordinate then has no direct Brownian forcing and no direct score correction in (4.23). Therefore, we use different discretizations in the two cases. 4.3.1
Reverse-SDE discretization for H < 12
For H < 12 , the Laplace representation produces a Markovian approximation with ψ̄ > 0. Hence, both the primary diffusion loading g(t)ψ̄ dBt in (4.4) and the primary block of LZ L⊤ Z are nonzero. The reverse SDE supplies both stochastic forcing and score correction to the primary state. In this regime, we discretize (4.23) directly, for example by Euler–Maruyama or another standard SDE solver, using the score estimator from Proposition 4.2. This method is analogous to the standard score-based SDE samplers of Song et al. (2020) and to the augmented reverse-SDE sampler used in GFDM (Nobis et al., 2024). For stability, it is convenient to rescale the primary coordinate by ψ̄, so that its diffusion loading is comparable to that of the auxiliary factors. 4.3.2
Sampling in the smooth regime
b For H > 21 , the signed weights cancel at the origin: ψ̄ = i∈I ψi = G(0) = 0. Thus the primary coordinate has no direct Brownian forcing and no direct score correction in the augmented reverse SDE. This degeneracy is structural. Its numerical effect, however, depends on the lift size. For small lifts, such as N = 2 or N = 3 in our experiments, the auxiliary rates remain comparable and ordinary explicit reverse discretization is stable. For larger lifts, the finite-difference construction introduces rates spanning several orders of magnitude, and the deterministic coupling from the auxiliary variables to the primary state becomes stiff. In that stiff regime, we use the Gaussianbridge reconstruction sampler described in the Appendix D. P
15
Table 1: Weight sum ψ̄ under different quadrature budgets N . N
5
H
2
4
8
16
32
0.3 0.7 0.9
0.70 0.00 0.00
1.09 0.00 0.00
1.77 0.00 0.00
3.09 0.00 0.00
6.11 0.00 0.00
Numerical Study
In this section, we detail the numerical implementation of the proposed model. We first outline the core computational components, including the configuration of the quadrature rates and weights, the truncated pseudoinverse stabilization for the covariance matrix, the adaptive noise schedule, and the discrete-time sampling scheme. We subsequently validate the method on MNIST, where SVDM outperforms the evaluated baselines, and on CIFAR-10, where a preliminary small-lift experiment tests whether the same persistent Volterra perturbation can extend to natural images. Small persistent Volterra lifts can improve image generation quality, with the best MNIST results obtained at H = 0.9 and N = 2. Moreover, increasing the lift size is not automatically beneficial. Larger lifts can improve the kernel approximation but also introduce stiff rates, ill-conditioned auxiliary covariance matrices, and unstable reverse-time discretizations. The Gaussian-bridge sampler is introduced to stabilize these larger stiff lifts, while the best empirical performance is achieved in the non-stiff small-lift regime. 5.1
Factor Rates, Weights, and Stiffness
To understand the behavior of the factors Yi , this section presents the statistics of factor rates κi and weights ψi used in the augmented generative model. These quantities are not learned parameters. They are obtained from the quadrature construction in Section 3 and, in the smooth regime, together with the finite-difference steps. Table 1 reports the weight sum ψ̄ under a canonical diagnostic choice of quadrature parameters a = 1, b = 1, and α = 1.065. In this subsection, the sub-interval quadrature level is fixed to m = 1, which implicitly determines β as a function of H and N . In the rough regime H = 0.3, the quadrature weights are nonnegative and ψ̄ > 0. In the smooth regimes H = 0.7 and H = 0.9, the finite-difference construction produces signed exponential pairs that cancel at the origin, so ψ̄ = 0. Table 2 reports the largest rate κmax := maxi∈IN κi . Large values of κmax correspond √ to fast mean-reverting speed. The quadrature interval endpoint ξn grows exponentially in N , and the resulting quadrature nodes can span several orders of magnitude. The effect is especially pronounced in the rough regime, where the kernel is singular near the origin. The large rates in Table 2 create stiffness in the lifted system. A factor with rate κi evolves on the time scale κ−1 i . If a numerical step size is much larger than this time scale, an explicit integrator cannot resolve the factor accurately. This issue is common in Markovian lifts of fractional kernels because the approximation must represent behavior across a wide range of exponential decay scales. Hence, we should choose suitable sampling methods in the forward and reverse-time processes to tackle the instability caused by large mean-reverting rates. We present the implementation for 16
Table 2: Largest final rate κmax = maxi∈IN κi under different quadrature budgets N . N H
2
4
8
16
32
0.3 0.7 0.9
5.42 1.49 1.32
18.78 2.81 2.27
92.87 6.05 4.34
796.05 16.34 9.98
15343.90 62.26 30.50
simulating the augmented forward and reverse processes in Section 5.4 later. Before that, we address another numerical instability issue caused by possible singularity in factor covariance matrices. 5.2
Singularity in Factor Covariance
The factor covariance matrix Cyy (t) defined in Section 4.1 enters both the Gaussian conditioning coefficient η(t) and the auxiliary score in (4.20). Its non-singularity is therefore important for both training and reverse-time sampling. Since all auxiliary factors are driven by the same Brownian motion, Cyy (t) can be ill-conditioned in both regimes. Moreover, in the smooth regime, the finitedifference construction creates nearly paired exponential rates, which leads to small eigenvalues. We first give a simple calculation illustrating the source of this degeneracy. Consider one coordinate of two auxiliary factors with positive rates κi and κk . Suppressing the coordinate index, their covariance is Z Cov(Yti , Ytk ) =
t
e−(κi +κk )(t−s) g(s)2 ds,
0
which follows from the fact that the two factors are driven by the same Brownian component. If g(s) ≡ g is constant, then in the long run, we have i k Cov(Y∞ , Y∞ )=
g2 , κi + κk
i Var(Y∞ )=
g2 . 2κi
Hence, large rates produce small marginal variances. In addition, close rates lead to strong correlations. Indeed, the correlation is √ i ,Y k) 2 κi κk Cov(Y∞ ∞ q = . κi + κk i ) Var(Y k ) Var(Y∞ ∞ This is the geometric-to-arithmetic mean ratio of the two rates. It is close to one when κi and κk are close. Therefore, the conditioning of Cyy (t) is affected by both the scale of the rates and their spacing. Large rates generate small-variance factors, while nearby rates generate almost collinear factors. The H > 12 regime introduces a further source of ill-conditioning. Before relabelling the final exponential approximation, let (γk , ωk ) denote a positive Gaussian quadrature node and weight appearing in the intermediate representation. The finite-difference construction replaces this term by exponentials with the rate-weight pairs given by
(κk,− , ψk,− ) := γk − δ2 ,
ωk , δ1 + δ2
(κk,+ , ψk,+ ) := γk + δ1 , −
ωk . δ1 + δ2
Since κk,− and κk,+ are close, the corresponding factors are nearly collinear. This creates small eigenvalues in Cyy (t). In representative smooth-regime experiments, for example H = 0.7, the 17
smallest eigenvalues can be several orders of magnitude smaller than the largest ones, which makes raw inversion numerically unstable. To improve numerical stability, we use a truncated pseudoinverse for Cyy (t). Let Cyy (t) = V(t)Λ(t)V(t)⊤ be the spectral decomposition, with eigenvalues λ1 (t), . . . , λM (t). Denote λmax (t) := max λi (t). i
For a relative threshold τrel > 0, define τ (t) := τrel λmax (t). The truncated pseudoinverse is C†,τ yy (t) = V(t) diag
†,τ λ†,τ 1 (t), . . . , λM (t)
( ⊤
V(t) ,
λ†,τ i (t) =
λi (t)−1 , 0,
λi (t) > τ (t), λi (t) ≤ τ (t).
(5.1)
In the experiments, we use τrel = 10−3 . Thus, directions whose eigenvalue is below the numerical threshold are projected out rather than inverted. This is a rank truncation, instead of a ridge regularization. In the exact Gaussian conditioning formula, Cyy (t)−1 is used when the matrix is nonsingular, and the Moore–Penrose pseudoinverse can be used for a degenerate Gaussian law. In the numerical implementation, these inverses are replaced by C†,τ yy (t). Accordingly, the implemented regression coefficient is η τ (t) := C†,τ yy (t)cyx (t), and the auxiliary Gaussian score is evaluated in the stabilized form
∇y log qt (YtJ ) ≈ − C†,τ yy (t) ⊗ ID
h
i
YtJ − uJ (t) ⊗ 1D .
The truncation removes directions whose eigenvalue is below numerical resolution. This avoids amplifying tiny eigenvalues in the regression coefficients and in the auxiliary Gaussian score, at the cost of introducing a small stabilization bias in the discarded directions. 5.3
Noise Schedule
This section describes the choice of the deterministic schedules µ and g used in the augmented forward process. Design criteria. The schedules µ and g affect both the amount of injected noise and the deterministic drift of the augmented system. We choose them according to the following criteria. First, the schedule should be positive and simple to evaluate, since it is used repeatedly in the Gaussian covariance formulas. Second, the terminal primary innovation variance should be normalized across Hurst regimes and lift sizes. Otherwise, differences in sample quality or numerical stability may simply reflect different amounts of injected noise rather than properties of the Markovian lift. Third, the normalization should be based on the actual innovation variance vx (T ) := Var(XT − ρ(T )X0 | X0 ) , rather than on the direct diffusion loading g(t)ψ̄. This is essential in the smooth regime, where ψ̄ = 0. Fourth, the same construction should apply to both rough and smooth regimes without using a separate scale parameter for each H and N . Appendix C shows the design details and gives g(t) and µ(t) in (C.5). 18
N=4
N=8
N=10
1.0
1.0
0.8
0.8
0.6
0.6
Variance
Variance
N=2
0.4
0.2
0.0
0.0 0.2
0.4
Time
0.6
0.8
1.0
N=4
N=8
0.2
0.4
0.6
N=10
0.4
0.2
0.0
N=2
0.0
(a) H = 0.3 (rough regime).
Time
0.8
1.0
(b) H = 0.7 (smooth regime).
Figure 1: Primary innovation variance vx (t) under the noise schedule (C.5). The Hurst exponents are H ∈ {0.3, 0.7}, and the quadrature budgets are N ∈ {2, 4, 8, 10}. Numerical validation.
For D = 1, the primary innovation variance at time t is vx (t) =
Z th X 0
i2
ψi e−κi (t−s) g(s) ds.
(5.2)
i∈IN
For general D, the corresponding innovation covariance is vx (t)ID . Figure 1 reports vx (t) for H ∈ {0.3, 0.7} and N ∈ {2, 4, 8, 10}. As predicted by (C.5), all curves satisfy vx (T ) = 1. This confirms that the schedule removes the scale distortion introduced by different Markovian approximations. Once the terminal variance is fixed, the remaining differences in covariance conditioning and reverse-time stability can be attributed to the structure of the lifted system rather than to a trivial change in the amount of injected noise. 5.4
Sampling Schemes
The forward system (Xt , YtJ ) has explicit Gaussian marginal and conditional laws. Therefore, for each sampled time t, forward sampling is implemented directly by X0 ∼ p0 ,
YtJ ∼ qt ,
Xt | (X0 , YtJ ) ∼ N mx|y (t), cx|y (t)ID .
This direct Gaussian sampling avoids the stability restrictions that would arise from explicitly integrating high-rate auxiliary factors with an Euler scheme. Reverse-time generation follows the dynamics in Section 4.3: 1. When H < 21 , the final exponential weights are nonnegative and the primary diffusion loading g(t)ψ̄ is nonzero. In this case, we use a standard Euler–Maruyama discretization of the augmented reverse SDE, together with the augmented score estimator from Proposition 4.2. For numerical stability, the primary coordinate Xt can be rescaled by ψ̄, so that the primary and auxiliary diffusion loadings have comparable magnitude. 2. When H > 21 , the hybrid finite-difference construction gives ψ̄ = 0. The primary coordinate then has neither direct stochastic forcing nor direct score correction in a naive reverse Euler 19
step. In our experiments, the numerical treatment depends on the lift size. For N = 2 and N = 3, the auxiliary rates are comparable and the explicit reverse discretization remains stable; these are also the configurations that give the best FID scores. For larger lifts, such as N ≥ 4, the finite-difference rates become stiff, and we use the Gaussian-bridge reconstruction sampler in Appendix D. The bridge step should therefore be understood as a stabilizer for stiff smooth-regime lifts rather than as the default sampler for all H > 12 configurations. These numerical procedures are designed to stabilize sampling without changing the finite kernel approximation. Hence, the approximation error bounds in Section 3 remain applicable to the underlying continuous-time Markovian approximation. Numerical validation. We now compare the naive Euler reverse step with the bridge-based auxiliary update used in the sampler of Section 4.3.2. The test configuration is H = 0.7 and N = 8. The implemented sampler uses the full joint Gaussian bridge for YtJ . To visualize the stability mechanism, however, it is useful to look at the scalar marginal bridge coefficient of a single factor. For a single coordinate of a factor with rate κi , define Z t
vi (t) :=
e−2κi (t−u) g(u)2 du.
0
For a reverse step from forward time t to s < t, the scalar bridge coefficient is Ki (s, t) = e−κi (t−s)
vi (s) . vi (t)
This coefficient is the Gaussian regression coefficient in the conditional law of Ysi given Yti , where only a single coordinate is considered for simplicity. It determines how strongly the later-time factor value influences the earlier-time factor value in the bridge. Unlike the explicit reverse Euler amplification, it remains bounded for stiff factors. Figure 2 compares the Euler amplification factor with the bridge coefficient under the rates κi ∈ {42, 284}. The difference is substantial. The explicit Euler amplification grows rapidly and crosses the float32 overflow threshold for the largest rates after a small number of reverse steps. By contrast, the bridge coefficients remain bounded throughout the simulation. This illustrates why the bridge sampler is needed for larger smoothregime lifts. Combined with the DDIM-type reconstruction of the primary state, this yields a stable reverse sampler for stiff larger smooth-regime lifts. 5.5
Image Generation on MNIST
We evaluate the proposed stochastic Volterra diffusion model (SVDM) on MNIST. This experiment has two purposes. First, it compares the Volterra noising mechanism with Brownian and fractionaldiffusion baselines. Second, it examines how the Hurst parameter H and the Markovian lift size N affect generation quality. We use a conditional U-Net architecture following Ronneberger et al. (2015), with attention resolutions [4, 2], three residual blocks, and channel multipliers [1, 2, 2, 2, 2]. The learning rate of the Adam optimizer is controlled by the OneCycle schedule of Smith & Topin (2019). All models are trained for 50,000 iterations with batch size 1024 and maximum learning rate 10−4 . We do not use exponential moving average. 20
= 42
= 284
100
1.0
80
0.8
OU bridge coefficient
Euler amplification
= 42
60 40 20
= 284
0.6 0.4 0.2 0.0
0 0
200
400
600
Step index
800
0
1000
(a) Explicit Euler amplification.
200
400
600
Step index
800
1000
(b) Gaussian bridge coefficient.
Figure 2: Reverse-step stability with H = 0.7 and N = 8. The left panel shows log10 ((1+κi ∆t)step ) for explicit Euler amplification; the shaded region marks the float32 overflow range. The right panel shows the Gaussian bridge coefficient Ki (s, t).
To isolate the effect of the forward noising mechanism, the Brownian SDE baseline and SVDM use the same network architecture and training protocol. We consider H ∈ {0.3, 0.7, 0.9} and N ∈ {2, 3, 4, 5, 6}. The hyperparameters in Gaussian quadrature and finite difference are set as a = 4.108, b = 10.96, α = 1.801, β = 1.318, and δ1 = δ2 = 0.1. We also compare with GFDM (Nobis et al., 2024) under matched (H, N ) configurations. Unless otherwise stated, SVDM uses the noise schedule in Section 5.3. For H > 12 , the sampler is chosen according to the lift size: we use ordinary explicit reverse discretization for the non-stiff small-lift configurations N = 2, 3, and the Gaussian-bridge reconstruction sampler for larger stiff lifts. The best-performing configurations in Table 3 are small-lift runs and therefore use the explicit reverse discretization. It is important to note that the best configurations, N = 2 and N = 3, do not rely on the Gaussian-bridge sampler. In these small-lift regimes the auxiliary rates are sufficiently moderate that explicit reverse discretization is stable. The strong performance therefore comes from the persistent Volterra noising mechanism itself, rather than from the bridge stabilization. The bridge sampler becomes relevant for larger N , where the Markovian lift is closer to the fractional kernel but the reverse dynamics become stiff. In Table 3, the best result is obtained at H = 0.9, N = 2, with FID = 0.52 ± 0.0024. This is the lowest FID among all models and configurations reported in Table 3. It improves substantially over the Brownian SDE baseline, which has FID 7.49, and over the best matched GFDM configuration, whose FID is 7.34. The improvement remains strong at N = 3, where SVDM obtains FID values 0.68 for H = 0.7 and 0.67 for H = 0.9. The dependence on N is non-monotone. In the persistent regime, increasing N beyond 2 or 3 does not improve sample quality and eventually degrades performance. This agrees with the diagnostic results: a larger lift can improve the nominal kernel approximation, but it also worsens covariance conditioning and reverse-time stiffness. In the rough regime H = 0.3, SVDM is better than GFDM at N = 2, but both methods perform worse than in the persistent regime. This is consistent with the large-rate behavior of the rough-regime Markovian lift. 21
Table 3: MNIST generation quality measured by FID (↓). The Brownian SDE row is independent of H and N . Within each method–H row, the best value over N is shown in bold. The overall best value is achieved by SVDM with H = 0.9 and N = 2. Method
H
SDE
–
GFDM
0.3 0.7 0.9
11.82 ± 0.0412 8.43 ± 0.0409 8.44 ± 0.0383
14.95 ± 0.0391 8.23 ± 0.0389 7.74 ± 0.0412
SVDM
0.3 0.7 0.9
9.08 ± 0.0243 0.54 ± 0.0050 0.52 ± 0.0024
23.22 ± 0.0378 0.68 ± 0.0038 0.67 ± 0.0061
N=2
N=3
N=4
N=5
N=6
16.68 ± 0.0387 7.72 ± 0.0374 19.83 ± 0.0175
17.00 ± 0.0363 7.36 ± 0.0252 18.11 ± 0.0642
20.72 ± 0.0445 7.34 ± 0.0241 27.53 ± 0.0571
26.67 ± 0.0421 7.49 ± 0.0387 7.56 ± 0.0416
26.59 ± 0.0394 8.72 ± 0.0352 9.44 ± 0.0480
28.77 ± 0.0436 7.08 ± 0.0158 24.94 ± 0.0641
7.49 ± 0.0142
Table 4: Contextual comparison with representative published MNIST FID reports. Model SVDM, H = 0.9, N = 2 cU-Net with DDIM VPGA GANetic-loss DCGAN Brownian SDE baseline GFDM
Reported MNIST FID
Source
0.52 ± 0.0024 2.55 4.97 ± 0.07 5.96 7.49 ± 0.0142 7.34 ± 0.0241
This work Ordoñez et al. (2024) Zhang et al. (2020) Akhmedova & Körber (2025) This work This work, following Nobis et al. (2024)
Moreover, Table 4 compares the best SVDM result with representative published MNIST FID values. These external numbers are not protocol-matched, because FID implementations, sample counts, architectures, and conditioning mechanisms differ across papers. Nevertheless, they indicate that the FID 0.52 obtained by SVDM is highly competitive among reported results. Figure 3 in Appendix E shows class-conditional MNIST samples generated by SVDM at H = 0.9 for N ∈ {2, 4, 6}. This Hurst parameter gives the best quantitative result in Table 3. At N = 2, the generated digits are sharp and visually stable across classes. At N = 4, most digits remain recognizable, but stroke thickness and local shape deformation become less consistent. At N = 6, the degradation is more pronounced. This qualitative trend matches the FID scores and supports the conclusion that, in the persistent regime, a small Markovian lift is sufficient and often preferable. 5.6
Image Generation on CIFAR-10
We also conduct an experiment on CIFAR-10 to test whether the proposed Volterra noising mechanism extends beyond grayscale digit generation. Following the MNIST findings, we use a persistent small-lift configuration with H = 0.9 and N = 2. The Gaussian quadrature and finite-difference parameters are the same as before. Figure 4 in Appendix E shows 600 generated samples arranged in a 30 × 20 grid. The samples cover multiple CIFAR-10 categories, including animals, vehicles, ships, and aircraft. The FID is approximately 9.5 under this configuration. This result should be interpreted as a preliminary validation rather than a fully optimized CIFAR-10 benchmark. It nevertheless supports the main empirical observation from MNIST: persistent Volterra perturbations can be effective with a small 22
Table 5: Contextual comparison of CIFAR-10 FID values. The numbers are not strictly protocolmatched: architectures, training budgets, conditioning, EMA usage, sample counts, and sampling steps differ across papers. The SVDM result is a preliminary single-configuration experiment. Model
CIFAR-10 FID (↓)
Source
9.5 8.99 4.85 5.20 3.17 2.20 1.97 1.91 3.86 5.31
This work Nobis et al. (2024) Nobis et al. (2024) Nobis et al. (2024) Ho et al. (2020) Song et al. (2020) Karras et al. (2022) Xu et al. (2023) Ren et al. (2024) Zhou et al. (2024)
SVDM, H = 0.9, N = 2 GFDM, FVP H = 0.9, K = 2 Brownian VP retrained baseline Brownian VE retrained baseline DDPM Score-SDE / NCSN++ EDM, unconditional PFGM++, unconditional MDSS with DDIM, 100 steps SFERD, one-step distillation
Markovian lift, while avoiding the stiffness and covariance-conditioning issues that arise for larger N. Table 5 provides a contextual comparison with representative CIFAR-10 generative results. The comparison is not protocol-matched, and the SVDM experiment is not yet optimized for CIFAR-10. Nevertheless, it supports the qualitative conclusion from MNIST: persistent Volterra perturbations can be useful with a small Markovian lift.
6
Conclusion
This paper introduced Volterra generative models, a score-based generative framework in which the standard Brownian forward perturbation is replaced by a stochastic Volterra perturbation. The Hurst parameter H provides a principled way to control the temporal structure of the noising mechanism, ranging from rough perturbations when H < 12 to persistent perturbations when H > 12 . To make the resulting non-Markovian and non-semimartingale forward process tractable, we developed finite-dimensional Markovian approximations of the fractional kernel. On MNIST, persistent Volterra perturbations with small lift sizes substantially improve sample quality over the Brownian baseline and GFDM in the tested configurations. A preliminary CIFAR-10 experiment further suggests that the same small-lift persistent regime can be extended to natural-image generation, although the CIFAR-10 result is not yet optimized. Several limitations remain. The smooth-regime approximation still relies on finite differences, which introduce signed weights and cancellation. Future work should study better-conditioned exponential approximations, quantify the sampling bias induced by covariance regularization and bridge-based updates for stiff larger lifts, adaptively select H, the lift size, and the anchor factor, and investigate whether a meaningful limiting reverse-time object exists as the number of Markovian factors tends to infinity. Acknowledgments Bingyan Han is partially supported by The Hong Kong University of Science and Technology (Guangzhou) Start-up Fund G0101000197. This work is also supported by the HPC AI Intelligent 23
Computing Platform of The Hong Kong University of Science and Technology (Guangzhou). We thank Nobis et al. (2024) for making their source code publicly available, which facilitated our numerical implementation and comparison.
References Shakhnaz Akhmedova and Nils Körber. Ganetic loss for generative adversarial networks with a focus on medical applications. Neural Computing and Applications, 37(23):19177–19197, 2025. Aurélien Alfonsi and Ahmed Kebaier. Approximation of stochastic Volterra equations with kernels of completely monotone type. Mathematics of Computation, 93(346):643–677, 2024. Brian DO Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3):313–326, 1982. Christian Bayer and Simon Breneis. Markovian approximations of stochastic Volterra equations with the fractional kernel. Quantitative Finance, 23(1):53–70, 2023. H Brass and G Hämmerlin. Numerical integration IV (oberwolfach, 1992), 1993. Helmut Brass and Knut Petras. Quadrature Theory: The Theory of Numerical Integration on a Compact Interval. Number 178 in Mathematical Surveys and Monographs. American Mathematical Society, 2011. Philippe Carmona and Laure Coutin. Fractional Brownian motion and the Markov property. Electronic Communications in Probability, 3:95–107, 1998. Hans Föllmer. Time reversal on Wiener space. In Stochastic Processes—Mathematics and Physics: Proceedings of the 1st BiBoS-Symposium held in Bielefeld, West Germany, September 10–15, 1984, pp. 119–129. Springer, 2006. Jim Gatheral, Thibault Jaisson, and Mathieu Rosenbaum. Volatility is rough. Quantitative Finance, 18(6):933–949, 2018. Gene H. Golub and John H. Welsch. Calculation of gauss quadrature rules. Mathematics of Computation, 23(106):221–230, 1969. doi: 10.1090/S0025-5718-69-99647-1. Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840–6851, 2020. Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005. Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusionbased generative models. In Advances in Neural Information Processing Systems, volume 35, pp. 26565–26577, 2022. Derrick H Lehmer. On the maxima and minima of Bernoulli polynomials. The American Mathematical Monthly, 47(8):533–538, 1940. Benoit B Mandelbrot and John W Van Ness. Fractional Brownian motions, fractional noises and applications. SIAM Review, 10(4):422–437, 1968. 24
Gabriel Nobis, Maximilian Springenberg, Marco Aversa, Michael Detzel, Rembert Daems, Roderick Murray-Smith, Shinichi Nakajima, Sebastian Lapuschkin, Stefano Ermon, Tolga Birdal, et al. Generative fractional diffusion models. Advances in Neural Information Processing Systems, 37: 25469–25509, 2024. Sergio Calvo Ordoñez, Chun-Wun Cheng, Jiahao Huang, Lipei Zhang, Guang Yang, Carola-Bibiane Schönlieb, and Angelica I Aviles-Rivero. The missing U for efficient diffusion models. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. URL https://openreview.net/forum? id=Y4YWzBiTEV. Zhiyao Ren, Yibing Zhan, Liang Ding, Gaoang Wang, Chaoyue Wang, Zhongyi Fan, and Dacheng Tao. Multi-step denoising scheduled sampling: Towards alleviating exposure bias for diffusion models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pp. 4667– 4675, 2024. doi: 10.1609/aaai.v38i5.28267. Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computerassisted intervention, pp. 234–241. Springer, 2015. René L. Schilling, Renming Song, and Zoran Vondraček. Bernstein Functions: Theory and Applications, volume 37 of De Gruyter Studies in Mathematics. Walter de Gruyter, Berlin, 2 edition, 2012. Leslie N Smith and Nicholay Topin. Super-convergence: Very fast training of neural networks using large learning rates. In Artificial intelligence and machine learning for multi-domain operations applications, volume 11006, pp. 369–386. SPIE, 2019. Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 2256–2265. PMLR, 2015. Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021, 2021. Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In Advances in Neural Information Processing Systems, volume 32, 2019. Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456, 2020. Ruslan Leont’evich Stratonovich. Conditional Markov processes. In Non-linear transformations of stochastic processes, pp. 427–453. Elsevier, 1965. Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Computation, 23(7):1661–1674, 2011. Yilun Xu, Ziming Liu, Yonglong Tian, Shangyuan Tong, Max Tegmark, and Tommi Jaakkola. PFGM++: Unlocking the potential of physics-inspired generative models. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 38566–38591. PMLR, 2023. 25
Eunbi B. Yoon, Keehun Park, Sungwoong Kim, and Sungbin Lim. Score-based generative models with Lévy processes. In Advances in Neural Information Processing Systems, volume 36, pp. 40694–40707, 2023. Zhisheng Zhang, Yuting Xie, and Lin Yang. Perceptual generative autoencoders. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 11298–11306. PMLR, 2020. Shengzhe Zhou, Zejian Li, Shengyuan Zhang, Lefan Hou, Changyuan Yang, Guang Yang, Zhiyuan Yang, and Lingyun Sun. Reducing spatial fitting error in distillation of denoising diffusion models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pp. 7686–7694, 2024. doi: 10.1609/aaai.v38i7.28602.
A
Gaussian Quadrature
This appendix recalls the Gaussian quadrature facts used in the smooth-regime error analysis. We rely on Gaussian quadrature and Peano kernels from Brass & Petras (2011, Chapters 4 and 6). The specific Peano-kernel bound used below follows Brass & Hämmerlin (1993, Theorem 2), and the Bernoulli function maximum is taken from Lehmer (1940, Theorem 1). Let [u, v] be a finite non-degenerate interval and let w : [u, v] → (0, ∞) be continuous. The Gaussian quadrature rule of level m for the weight w consists of nodes x1 , . . . , xm ∈ (u, v) and positive weights η1 , . . . , ηm such that Z v
p(x)w(x) dx = u
m X
ηk p(xk )
k=1
for every polynomial p of degree at most 2m − 1. Equivalently, the nodes are the roots of the degree-m orthogonal polynomial associated with the inner product ⟨f, g⟩w :=
Z v
f (x)g(x)w(x) dx. u
The existence of the orthogonal polynomial and the location of its roots follow from Brass & Petras (2011, Theorems A.1–A.2); the Gaussian exactness and positivity are stated in Brass & Petras (2011, Theorem 6.1.2). The corresponding weights are uniquely determined and may be represented by the standard Gaussian-weight formulas in Brass & Petras (2011, Theorem 6.1.3). In numerical implementations, the nodes and weights are typically computed by the Golub–Welsch algorithm (Golub & Welsch, 1969). We use the Peano representation for quadrature remainders from Brass & Petras (2011, Theorem 4.2.5). Since an m-point Gaussian quadrature rule is exact for polynomials of degree at most 2m − 1, its remainder admits the representation Z v
f (x)w(x) dx −
u
m X
Z v
ηk f (xk ) =
k=1
f (2m) (x)K2m (x) dx,
(A.1)
u
for f ∈ C 2m ([u, v]), where K2m is the Peano kernel associated with the quadrature remainder. The following bound for this Peano kernel is a standard Gaussian-quadrature estimate: sup |K2m (x)| ≤ x∈[u,v]
(2π)2m (2m)!
v−u 2
2m
26
sup |B2m (x)| sup w(x), x∈[−1,1]
x∈[u,v]
(A.2)
where B2m is the Bernoulli function (Brass & Hämmerlin, 1993, Theorem 2). We use the Fourierseries convention ∞ X cos 2πℓx − πs 2 Bs (x) = −2 . (2πℓ)s ℓ=1 For even s, Lehmer (1940, Theorem 1) gives sup |Bs (x)| = sup |Bs (x)| = x∈[−1,1]
x∈[0,1]
2ζ(s) . (2π)s
In particular, 2ζ(2m) π2 ≤ , (2π)2m 3(2π)2m
sup |B2m (x)| = x∈[−1,1]
because Riemann zeta function ζ(2m) ≤ ζ(2) = π 2 /6 for m ≥ 1.
B
Proofs of Results
B.1
Proof of Theorem 3.2
b H . Decompose the γProof. By (3.6), it suffices to bound the L2 (0, T ) error between GH and G integral into (0, ξ0 ), [ξ0 , ξn ], and (ξn , ∞), and write the corresponding errors as EL , EM , and ER . Since (x + y + z)2 ≤ 3(x2 + y 2 + z 2 ), Z T
b H (t)|2 dt ≤ 3 |GH (t) − G
0
Z T
|EL (t)|2 + |EM (t)|2 + |ER (t)|2 dt.
0
For the left tail, the choice of ω0 gives Z ξ0
EL (t) = cH
1
(e−tγ − 1)γ −H− 2 dγ.
0
Using |e−x − 1| ≤ x for x ≥ 0, cH t 23 −H ξ0 . 2 −H
|EL (t)| ≤ 3 Therefore, Z T
|EL (t)|2 dt ≤
0
2α √ c2H T 3 a3−2H exp − N . 3 2 A− 3( 2 − H)
For the right tail
Z ∞
ER (t) = cH
1
e−γt γ −H− 2 dγ,
ξn
we regard its squared L2 norm as a double integral. Then, Fubini’s theorem and the inequality √ γ + η ≥ 2 γη imply Z T 0
c2 c2 2α √ |ER (t)| dt ≤ H2 ξn−2H = H2 b−2H exp − N . 2H 2H A−
2
For the middle interval, we use the same argument as in Bayer & Breneis (2023, Lemmas 2.8–2.9), with the actual geometric ratio 1/n ξn r−,N := ξ0 27
in place of eαβ . More precisely, on each interval [ξj , ξj+1 ], the single-interval Gaussian quadrature estimate gives !2m+1 !1/2 5π 3 cH ξj+1 H− 12 |Ej (t)| ≤ . −1 t 18 22m+1 mH ξj Since the grid is geometric, ξj+1 /ξj = r−,N for all j = 0, . . . , n − 1. Hence 5π 3 |EM (t)| ≤ n 18
!1/2
1 cH tH− 2 (r−,N − 1)2m+1 . 2m+1 H 2 m
Squaring and integrating over [0, T ] yields Z T
|EM (t)|2 dt ≤
0
5π 3 c2H T 2H n2 36 H m2H
r−,N − 1 2
4m+2
.
Combining the three estimates and applying (3.6) proves (3.9). B.2
Proof of Theorem 3.3
Lemma B.1 (Single-interval quadrature error). Let H ∈ ( 12 , 1) and let (γk , ωk )m k=1 be the level1 m Gaussian quadrature nodes and weights on [u, v] with respect to the weight wH (γ) = cH γ 2 −H . Then, for all t > 0, Z v
cH
1
te−tγ γ 2 −H dγ −
u
m X
ωk te−γk t ≤
k=1
7π 3/2 cH m1−H H− 1 t 2 6 22m+1
2m+1 v −1 . u
Proof. Let f (γ) = te−tγ . Then f (2m) (γ) = t2m+1 e−tγ . Applying (A.1)–(A.2) with wH and using 1 that γ 2 −H is decreasing on [u, v], we obtain Z v
cH
1
te−tγ γ 2 −H dγ −
m X
u
ωk te−γk t
k=1 2 1 π cH t2m+1 e−tu u 2 −H (v − u)2m+1 . ≤ 3 · 22m (2m)!
Writing (v − u)2m+1 = u2m+1 (v/u − 1)2m+1 and using −x
e
≤
η η
e
x−η ,
x > 0,
η > 0,
with η = 2m + 23 − H and x = tu, gives 2m+ 32 −H
t2m+1 e−tu u Stirling’s lower bound (2m)! > 1+
√
3 2 −H
2m
2m + 23 − H e
≤
!2m+ 3 −H 2
1
tH− 2 .
4πm(2m/e)2m , together with !2m
≤e
3 −H 2
3 H 2+ − 2m m
,
for m ≥ 1 and H ∈ ( 12 , 1), yields the stated constant. 28
3 −H 2
7 ≤ , 2
Lemma B.2 (Middle-interval quadrature error). Let r+,N = (ξn /ξ0 )1/n . In the setting of Theorem 3.3, Z ξn
Z T
cH
1
te−tγ γ 2 −H dγ −
ξ0
0
mn X
2
ωi te−γi t dt ≤
i=1
49π 3 c2H T 2H 2 2−2H n m 72 H
r+,N − 1 2
4m+2
.
Proof. On each interval [ξj , ξj+1 ], Lemma B.1 applies with v/u = ξj+1 /ξj = r+,N . Summing the pointwise bounds over j = 0, . . . , n − 1 gives Z ξn
te
cH
−tγ
γ
1 −H 2
dγ −
ξ0
mn X
ωi te−γi t ≤
i=1
7π 3/2 cH nm1−H H− 1 t 2 (r+,N − 1)2m+1 . 6 22m+1
Squaring and integrating over t ∈ [0, T ] yields the result, since Z T
t2H−1 dt =
0
T 2H . 2H
Lemma B.3 (Global quadrature error). In the setting of Theorem 3.3, Z T 0
|GH (t) − GH (t)|2 dt ≤ Q+ N,T ,
where Q+ N,T is defined in (3.15). Proof. Decompose the γ-integral into (0, ξ0 ), [ξ0 , ξn ], and (ξn , ∞), and write the corresponding errors as EL (t), EM (t), and ER (t). Then
|GH (t) − GH (t)|2 ≤ 3 |EL (t)|2 + |EM (t)|2 + |ER (t)|2 . For the left tail, the zero node gives Z ξ0
EL (t) = cH
1
t(e−tγ − 1)γ 2 −H dγ.
0
Using |e−x − 1| ≤ x, cH t2 25 −H . |EL (t)| ≤ 5 ξ0 2 −H Hence
Z T 0
c2 T 5 2α √ |EL (t)| dt ≤ 5 H a5−2H exp − N . 2 A+ 5( 2 − H)
2
For the right tail, Fubini’s theorem gives Z T 0
Since
|ER (t)|
2
dt = c2H Z ∞ 0
Z ∞ Z ∞ "Z T
# 2 −t(γ+η)
t e ξn
ξn
t2 e−t(γ+η) dt =
1
0
2 1 ≤ , 3 (γ + η) 4(γη)3/2 29
1
dt γ 2 −H η 2 −H dγ dη.
we obtain
Z T
|ER (t)|2 dt ≤
0
c2H 4
Z ∞
γ −H−1 dγ
2
c2H −2H 2α √ b exp − N . 2 4H A+
=
ξn
The middle term is controlled by Lemma B.2. Multiplying the three tail and middle estimates by the factor 3 from the inequality at the start of the proof gives (3.15). Proof of Theorem 3.3. By Lemma B.3, Z T 0
|GH (t) − GH (t)|2 dt ≤ Q+ N,T .
It remains to control the finite-difference error. For the zero node, t−
δ1 t2 1 − e−δ1 t ≤ , δ1 2
t ≥ 0.
For each positive node, the mean-value theorem gives a point γi′ ∈ (γi − δ2 , γi + δ1 ) such that e−(γi −δ2 )t − e−(γi +δ1 )t ′ = te−γi t . δ1 + δ 2 Condition (3.12) ensures that γi′ ≥ 0. Therefore ′
te−γi t − te−γi t ≤ t2 max{δ1 , δ2 }. Using the nonnegativity of the quadrature weights, e H (t)| ≤ |GH (t) − G
ω0 δ1 2 t + WN+ t2 max{δ1 , δ2 } = DN t2 . 2
Consequently, Z T 0
e H (t)|2 dt ≤ D2 |GH (t) − G N
Z T
t4 dt =
0
T5 2 D . 5 N
Combining this estimate with Z T
e H (t)|2 dt ≤ 2 |GH (t) − G
0
Z T
2
|GH (t) − GH (t)| dt + 2
0
Z T
e H (t)|2 dt |GH (t) − G
0
proves (3.17). The process-level bound (3.18) then follows directly from (3.6). B.3
Proof of Proposition 4.2
Proof. We prove the claim for the exact Gaussian conditioning coefficients, while the truncated pseudoinverse used in the implementation does not yield an exact score identity. Define Rt := ξt − ρ(t)X0 − νx (t)1D . By the definition of ξt in (4.16), Rt = Xt − ρ(t)X0 − νx (t)1D −
M X ℓ=1
30
ηℓ (t) Ytjℓ − ujℓ (t)1D .
Conditional on X0 , the vector
Rt , YtJ − uJ (t) ⊗ 1D
is Gaussian, since the augmented forward process is linear with deterministic coefficients. For each data coordinate, the covariance between the centered auxiliary vector and the residual innovation equals cyx (t) − Cyy (t)η(t) = 0, by the definition η(t) = Cyy (t)−1 cyx (t). Different data coordinates are driven by independent Brownian components, so the same orthogonality holds blockwise for the D-dimensional variables. Hence Rt is uncorrelated with YtJ conditional on X0 . Since the pair is jointly Gaussian, Rt is independent of YtJ conditional on X0 . Moreover, X0 is independent of the Brownian motion that generates both Rt and YtJ . Therefore ξt = ρ(t)X0 + νx (t)1D + Rt is independent of YtJ marginally. Now consider the affine change of variables (x, y) 7−→ (ξ, y),
ξ =x−
M X
ηℓ (t) yjℓ − ujℓ (t)1D ,
ℓ=1
where
y = y j1 , . . . , y jM
uses the same block ordering as YtJ . This map is invertible and has unit Jacobian determinant. Since ξt and YtJ are independent, the joint density of (Xt , YtJ ) can be written as pt (x, y) = pξt (ξ) qt (y), where pξt denotes the marginal density of ξt , and qt denotes the forward-time marginal density of YtJ . Differentiating this factorization gives ∇x log pt (x, y) = ∇ξ log pξt (ξ), because ∂ξ/∂x = ID . For the ℓ-th auxiliary block, ∂ξ = −ηℓ (t)ID , ∂yjℓ and hence the chain rule gives ∇yjℓ log pt (x, y) = −ηℓ (t)∇ξ log pξt (ξ) + ∇yjℓ log qt (y),
ℓ = 1, . . . , M.
The augmented denoising score-matching objective (4.19) trains sθ to approximate the marginal residual score ∇ξ log pξt . Substituting this approximation into the preceding identities yields the score decomposition ∇z log pt (Zt ) = Sθ (Zt , t) + 0D , ∇y log qt (YtJ ) , when the score model is exact. This is precisely (4.22).
31
C
Noise Schedule We first set the scale of the raw schedule. Let
Intrinsic raw scale.
κfloor := min{κi : κi > 0, i ∈ I}. Define the aggregate forcing scale as Ψ := |ψi⋆ | +
κi⋆ |ψj |. max{κj , κfloor } j∈J X
(C.1)
It can be understood as follows. First, the floor prevents division by zero. Second, for a positiverate factor, the deterministic part of the dynamics suggests the quasi-stationary approximation Ytj ≈ µ(t) κj 1D . It motivates the use of κi⋆ /κj in (C.1). Third, in the smooth regime, the weights ψj are signed and may cancel. Since such cancellation does not reflect the magnitude of intermediate deterministic coupling, we use absolute values in (C.1). The intrinsic scale is then given by β ⋆ :=
κi⋆ . Ψ
(C.2)
This balances the anchor mean-reversion scale κi⋆ against the aggregate forcing scale Ψ. Raw schedule and diffusion normalization. raw linear schedule ⋆ hraw min := 0.05 S β ,
⋆ hraw max := 2 S β ,
Given a strength parameter S > 0, define the
raw raw hraw (t) := hraw min + (hmax − hmin )
t . T
(C.3)
The constants 0.05 and 2 determine the shape of the linear profile, while Sβ ⋆ determines its raw p magnitude. We set graw (t) := hraw (t). The terminal primary innovation variance induced by graw is vraw (T ) :=
Z T "X 0
#2
ψi e−κi (T −r) graw (r)
dr.
(C.4)
i∈I
Assuming vraw (T ) > 0, we normalize the diffusion coefficient g(t) and set µ(t) as follows: graw (t) g(t) := p , vraw (T )
µ(t) := hraw (t).
(C.5)
If both µ and g 2 were rescaled by the same terminal-variance normalizer, then the global factor Sβ ⋆ would cancel from the final schedule. Hence, only g is normalized. The constants S and β ⋆ control the deterministic drift scale and the mean offset, while the stochastic innovation is normalized to a common terminal variance. Indeed, by construction, Z T "X 0
i∈I
#2 −κi (T −r)
ψi e
g(r)
1 dr = vraw (T )
Z T "X 0
#2
ψi e
−κi (T −r)
graw (r)
dr = 1.
i∈I
It remains valid even when the direct diffusion coefficient of the primary state vanishes as ψ̄ = 0. 32
D
Gaussian-bridge Reconstruction Sampler
When H > 1/2 and the finite-difference lift is stiff, we update the auxiliary variables YJ by a Gaussian bridge step and then reconstruct the primary variable X from the conditional Gaussian representation. To introduce the sampling scheme of YJ , we consider the Y -component of the augmented reverse SDE. Define
Ly (t) := 1M ⊗ g(t)ID ,
Ky := diag(κj1 , . . . , κjM ), so that
Ly (t)Ly (t)⊤ = g(t)2 1M 1⊤ M ⊗ ID .
The forward dynamics of YtJ can be written as dYtJ = ay (t, YtJ ) dt + Ly (t) dBt ,
ay (t, y) := −(Ky ⊗ ID )y + µ(t)1M ⊗ 1D .
(D.1)
Then the Y -component of the reverse SDE is J
J
h
J
dYt = −ay (T − t, Yt ) + Ly (T − t)Ly (T − t)⊤ ∇y log pT −t Xt , Yt
i
dt
(D.2)
+ Ly (T − t) dBt . In the same spirit of Proposition 4.2, we separate the drift by Bayes’ rule. Recall that qt denotes the forward-time marginal density of YtJ . Since pt (x, y) = qt (y) pt (x | y), we have ∇y log pt (x, y) = ∇y log qt (y) + ∇y log pt (x | y).
(D.3)
Substituting (D.3) into (D.2) gives a natural decomposition of the drift, and the full dynamics are J
J
h
J
dYt = −ay (T − t, Yt ) + Ly (T − t)Ly (T − t)⊤ ∇y log qT −t Yt
i
dt (D.4)
+ Ly (T − t) dBt J
+ Ly (T − t)Ly (T − t)⊤ ∇y log pT −t Xt | Yt dt. The first two lines are exactly the reverse SDE of the marginal auxiliary process. Since YtJ is linear Gaussian, its reverse transition from forward time t to s < t can be handled by the Gaussian bridge. More explicitly, define the cross-time covariance Z s
Cyy (s, t)ℓk :=
e−κjℓ (s−u) e−κjk (t−u) g(u)2 du,
s ≤ t.
0
The Gaussian bridge of the marginal YJ is
y YsJ YtJ = yt ∼ N ms|t , Cys|t ⊗ ID ,
(D.5)
where Ks,t := Cyy (s, t)Cyy (t, t)† , y ms|t := uJ (s) ⊗ 1D + (Ks,t ⊗ ID ) [yt − uJ (t) ⊗ 1D ] ,
Cys|t := Cyy (s, s) − Cyy (s, t)Cyy (t, t)† Cyy (t, s).
(D.6)
When Cyy (t, t) is nonsingular, the Moore–Penrose inverse † is the usual inverse and the bridge is exact. In degenerate Gaussian cases, the Moore–Penrose inverse gives the exact conditional 33
Gaussian law on the support. In the numerical implementation, we replace it by the truncated pseudoinverse, which gives a stabilized approximation of the exact bridge. Next, we note that the last line in (D.4) is not included in the Gaussian bridge. By the augmented score decomposition in Proposition 4.2, the conditional auxiliary score is approximated as follows: ∇y log pt (Xt | YtJ ) ≈ (−η1 (t)st , . . . , −ηM (t)st ) = −η(t) ⊗ st ,
st := sθ (ξt , t).
Since Ly (t)Ly (t)⊤ = g(t)2 (1M 1⊤ M ) ⊗ ID , we have
h
i
Ly (t)Ly (t)⊤ ∇y log pt (Xt | YtJ ) ≈ −g(t)2 (1M 1⊤ M ) ⊗ ID [η(t) ⊗ st ] = −g(t)2 1M 1⊤ M η(t) ⊗ st .
Equivalently, each auxiliary block receives the same correction 2
−g(t)
M X
!
ηk (t) st .
k=1
Therefore, after drawing the Gaussian bridge proposal (D.5), a first-order Euler step gives Ysjℓ ← Ysjℓ − g(t)2
M X
!
ηk (t) st (t − s),
ℓ = 1, . . . , M.
(D.7)
k=1
Finally, we reconstruct the primary variable X as follows. The conditional Gaussian representation in (4.17) implies that ε ∼ N (0, ID ),
ξt = ρ(t)X0 + νx (t)1D + σx|y (t)ε,
(D.8)
2 (t) = c 2 where σx|y x|y (t). Hence, the L -optimal predictor satisfies
"
#
ξt − ρ(t)X0 − νx (t)1D 1 s∗ (ξ, t) = E − ξt = ξ = − E[ε | ξt = ξ]. cx|y (t) σx|y (t) Since we train the network sθ (ξt , t) to estimate s∗ (ξ, t), an estimator of E[ε | ξt = ξ] is given by εθ (t) := −σx|y (t)sθ (ξt , t). Together with the representation (D.8), we obtain a Tweedie-type estimate of E[X0 | ξt ] given by b 0 below: X b0 = X
i i 1 h 1 h ξt − νx (t)1D − σx|y (t)εθ (t) = ξt − νx (t)1D + cx|y (t)sθ (ξt , t) . ρ(t) ρ(t)
After updating the auxiliary variables YJ to time s < t, we define the primary state Xs by b 0 + νx (s)1D + Xs = ρ(s)X
M X
ηℓ (s) Ysjℓ − ujℓ (s)1D + σx|y (s)εθ (t).
ℓ=1
34
This update should be interpreted as a DDIM-type plug-in update rather than as an exact conditional simulation. The Gaussian conditioning formula gives the law of Xs conditional on the true clean variable X0 and the auxiliary state YsJ . In the sampler, however, X0 is replaced by the b 0 , and the normalized residual noise is replaced by the score-based estimate denoising estimate X εθ (t) obtained at the previous time t. Thus, the reconstruction defines a deterministic coupling between the states at times t and s, in a similar spirit of the DDIM method in Song et al. (2021). It avoids the zero diffusion difficulty in X caused by ψ̄ = 0.
E
Generated Figure Samples
(a) H = 0.9, N = 2, FID = 0.52.
(b) H = 0.9, N = 4, FID = 7.56.
(c) H = 0.9, N = 6, FID = 24.94.
Figure 3: Class-conditional MNIST samples generated by SVDM in the best-performing Hurst regime H = 0.9. Each panel shows two rows containing the digits 0–9.
35
Figure 4: CIFAR-10 samples generated by SVDM with H = 0.9 and N = 2, displayed in a 30 × 20 grid. The FID is approximately 9.5.
36