Year: 2021

Abstract:
Score-based generative models (SGMs) have demonstrated remarkable synthesis quality. SGMs rely on a diffusion process that gradually perturbs the data towards a tractable distribution, while the generative model learns to denoise. The complexity of this denoising task is, apart from the data distribution itself, uniquely determined by the diffusion process. We argue that current SGMs employ overly simplistic diffusions, leading to unnecessarily complex denoising processes, which limit generative modeling performance. Based on connections to statistical mechanics, we propose a novel critically-damped Langevin diffusion (CLD) and show that CLD-based SGMs achieve superior performance. CLD can be interpreted as running a joint diffusion in an extended space, where the auxiliary variables can be considered “velocities” that are coupled to the data variables as in Hamiltonian dynamics. We derive a novel score matching objective for CLD and show that the model only needs to learn the score function of the conditional distribution of the velocity given data, an easier task than learning scores of the data directly. We also derive a new sampling scheme for efficient synthesis from CLD-based diffusion models. We find that CLD outperforms previous SGMs in synthesis quality for similar network architectures and sampling compute budgets. We show that our novel sampler for CLD significantly outperforms solvers such as Euler–Maruyama. Our framework provides new insights into score-based denoising diffusion models and can be readily used for high-resolution image synthesis.

heading: 1 INTRODUCTION
text: Score-based generative models (SGMs) and denoising diffusion probabilistic models have emerged as a promising class of generative models (Sohl-Dickstein et al., 2015; Song et al., 2021c;b; Vahdat et al., 2021; Kingma et al., 2021). SGMs offer high quality synthesis and sample diversity, do not require adversarial objectives, and have found applications in image (Ho et al., 2020; Nichol & Dhariwal, 2021; Dhariwal & Nichol, 2021; Ho et al., 2021), speech (Chen et al., 2021; Kong et al., 2021; Jeong et al., 2021), and music synthesis (Mittal et al., 2021), image editing (Meng et al., 2021; Sinha et al., 2021; Furusawa et al., 2021), super-resolution (Saharia et al., 2021; Li et al., 2021), image-to-image translation (Sasaki et al., 2021), and 3D shape generation (Luo & Hu, 2021; Zhou et al., 2021). SGMs use a diffusion process to gradually add noise to the data, transforming a complex data distribution to an analytically tractable prior distribution. A neural network is then utilized to learn the score function—the gradient of the log probability density—of the perturbed data. The learnt scores can be used to solve a stochastic differential equation (SDE) to synthesize new samples. This corresponds to an iterative denoising process, inverting the forward diffusion.
In the seminal work by Song et al. (2021c), it has been shown that the score function that needs to be learnt by the neural network is uniquely determined by the forward diffusion process. Consequently, the complexity of the learning problem depends, other than on the data itself, only on the diffusion. Hence, the diffusion process is the key component of SGMs that needs to be revisited to further improve SGMs, for example, in terms of synthesis quality or sampling speed.
Inspired by statistical mechanics (Tuckerman, 2010), we propose a novel forward diffusion process, the critically-damped Langevin diffusion (CLD). In CLD, the data variable, xt (time t along the diffusion), is augmented with an additional “velocity” variable vt and a diffusion process is run in the joint data-velocity space. Data and velocity are coupled to each other as in Hamiltonian dynamics, and noise is injected only into the velocity variable. Similarly as in Hamiltonian Monte
Carlo (Neal, 2011), the Hamiltonian component helps to efficiently traverse the joint data-velocity space and to transform the data distribution into the prior distribution more smoothly. We derive the corresponding score matching objective and show that for CLD the neural network is tasked with learning only the score of the conditional distribution of velocity given data∇vt log pt(vt|xt), which is arguably easier than learning the score of diffused data distribution directly. Using techniques from molecular dynamics (Bussi & Parrinello, 2007; Tuckerman, 2010; Leimkuhler & Matthews, 2013), we also derive a novel SDE integrator tailored to CLD’s reverse-time synthesis SDE.
We extensively validate CLD and the new SDE solver: (i) We show that the neural networks learnt in CLD-based SGMs are smoother than those of previous SGMs. (ii) On the CIFAR-10 image modeling benchmark, we demonstrate that CLD-based models outperform previous diffusion models in synthesis quality for similar network architectures and sampling compute budgets. We attribute these positive results to the Hamiltonian component in the diffusion and to CLD’s easier score function target, the velocity-data conditional distribution. (iii) We show that our novel sampling scheme for CLD significantly outperforms the popular Euler–Maruyama method. (iv) We perform ablations on various aspects of CLD and find that CLD does not have difficult-to-tune hyperparameters.
In summary, we make the following technical contributions: (i) We propose CLD, a novel diffusion process for SGMs. (ii) We derive a score matching objective for CLD, which requires only the conditional distribution of velocity given data. (iii) We propose a new type of denoising score matching ideally suited for scalable training of CLD-based SGMs. (iv) We derive a tailored SDE integrator that enables efficient sampling from CLD-based models. (v) Overall, we provide novel insights into SGMs and point out important new connections to statistical mechanics.

heading: 2 BACKGROUND
text: Consider a diffusion process ut ∈ Rd defined by the Itô SDE dut = f(ut, t) dt+G(ut, t) dwt, t ∈ [0, T ], (1)
with continuous time variable t ∈ [0, T ], standard Wiener process wt, drift coefficient f : Rd × [0, T ]→ Rd and diffusion coefficientG : Rd×[0, T ]→ Rd×d. Defining ūt := uT−t, a corresponding reverse-time diffusion process that inverts the above forward diffusion can be derived (Anderson, 1982; Haussmann & Pardoux, 1986; Song et al., 2021c) (with positive dt and t ∈ [0, T ]): dūt = [ −f(ūt, T − t) +G(ūt, T − t)G(ūt, T − t)>∇ūt log pT−t(ūt) ] dt+G(ūt, T − t)dwt, (2)
where∇ūt log pT−t(ūt) is the score function of the marginal distribution over ūt at time T − t. The reverse-time process can be used as a generative model. In particular, Song et al. (2021c) model data x, setting p(u0)=pdata(x). Currently used SDEs (Song et al., 2021c; Kim et al., 2021) have drift and diffusion coefficients of the simple form f(xt, t)=f(t)xt andG(xt, t)=g(t)Id. Generally, f and G are chosen such that the SDE’s marginal, equilibrium density is approximately Normal at time T , i.e., p(uT )≈N (0, Id). We can then initialize x0 based on a sample drawn from a complex data distribution, corresponding to a far-from-equilibrium state. While we let the state x0 relax towards equilibrium via the forward diffusion, we can learn a model sθ(xt, t) for the score ∇xt log pt(xt), which can be used for synthesis via the reverse-time SDE in Eq. (2). If f and G take the simple form from above, the denoising score matching (Vincent, 2011) objective for this task is:
min θ
Et∼[0,T ]Ex0∼p(x0)Ext∼pt(xt|x0) [ λ(t)‖sθ(xt, t)−∇xt log pt(xt|x0)‖22 ] (3)
If f and G are affine, the conditional distribution pt(xt|x0) is Normal and available analytically (Särkkä & Solin, 2019). Different λ(t) result in different trade-offs between synthesis quality and likelihood in the generative model defined by sθ(xt, t) (Song et al., 2021b; Vahdat et al., 2021).

heading: 3 CRITICALLY-DAMPED LANGEVIN DIFFUSION
text: We propose to augment the data xt ∈ Rd with auxiliary velocity1 variables vt ∈ Rd and utilize a diffusion process that is run in the joint xt-vt-space. With ut = (xt,vt)> ∈ R2d, we set
f(ut, t) :=
( 0 βM−1 −β −ΓβM−1 ) ⊗ Id ut , G(ut, t) := ( 0 0 0 √ 2Γβ ) ⊗ Id, (4)
where ⊗ denotes the Kronecker product. The coupled SDE that describes the diffusion process is( dxt dvt ) = ( M−1vt −xt ) βdt︸ ︷︷ ︸
Hamiltonian component=:H
+
( 0
−ΓM−1vt
) βdt+ ( 0√ 2Γβ ) dwt︸ ︷︷ ︸
Ornstein-Uhlenbeck process=:O
, (5)
which corresponds to Langevin dynamics in each dimension. That is, each xi is independently coupled to a velocity vi, which explains the blockwise structure of f and G. The mass M ∈ R+ is a hyperparameter that determines the coupling between the xt and vt variables; β ∈ R+ is a constant time rescaling chosen such that the diffusion converges to its equilibrium distribution within t ∈ [0, T ] (in practice, we set T=1) when initialized from a data-defined non-equilibrium state and is analogous to β(t) in previous diffusions (we could also use time-dependent β(t), but found constant β’s to work well, and therefore opted for simplicity); Γ ∈ R+ is a friction coefficient that determines the strength of the noise injection into the velocities. Notice that the SDE in Eq. (5) consists of two components. The H term represents a Hamiltonian component. Hamiltonian dynamics are frequently used in Markov chain Monte Carlo methods to accelerate sampling and efficiently explore complex probability distributions (Neal, 2011). The Hamiltonian component in our diffusion process plays a similar role and helps to quickly and smoothly converge the initial joint data-velocity distribution to the equilibrium, or prior (see Fig. 1). Furthermore, Hamiltonian dynamics on their own are trivially invertible (Tuckerman, 2010), which intuitively is also beneficial in our situation when using this diffusion for training SGMs. TheO term corresponds to an OrnsteinUhlenbeck process (Särkkä & Solin, 2019) in the velocity component, which injects noise such that the diffusion dynamics properly converge to equilibrium for any Γ>0. It can be shown that the equilibrium distribution of this diffusion is pEQ(u) = N (x; 0d, Id)N (v; 0d,MId) (see App. B.2). There is a crucial balance between M and Γ (McCall, 2010): For Γ2<4M (underdamped Langevin dynamics) the Hamiltonian component dominates, which implies oscillatory dynamics of xt and vt that slow down convergence to equilibrium. For Γ2>4M (overdamped Langevin dynamics) the Oterm dominates which also slows down convergence, since the accelerating effect by the Hamiltonian component is suppressed due to the strong noise injection. For Γ2=4M (critical damping), an ideal balance is achieved and convergence to pEQ(u) occurs as fast as possible in a smooth manner without oscillations (also see discussion in App. A.1) (McCall, 2010). Hence, we propose to set Γ2=4M and call the resulting diffusion critically-damped Langevin diffusion (CLD) (see Fig. 1).
Diffusions such as the VPSDE (Song et al., 2021c) correspond to overdamped Langevin dynamics with high friction coefficients Γ (see App. A.2). Furthermore, in previous works noise is injected directly into the data variables (pixels, for images). In CLD, only the velocity variables are subject to direct noise and the data is perturbed only indirectly due to the coupling between xt and vt.

heading: 3.1 SCORE MATCHING OBJECTIVE
text: Considering the appealing convergence properties of CLD, we propose to utilize CLD as forward diffusion process in SGMs. To this end, we initialize the joint p(u0)=p(x0) p(v0)=pdata(x0)N (v0; 0d, γMId) with hyperparameter γ<1 and let the distribution diffuse towards the tractable equilibrium—or prior—distribution pEQ(u). We can then
1We call the auxiliary variables velocities, as they play a similar role as velocities in physical systems. Formally, our velocity variables would rather correspond to physical momenta, but the term momentum is already widely used in machine learning and our mass M is unitless anyway.
learn the corresponding score functions and define CLD-based SGMs. Following a similar derivation as Song et al. (2021b), we obtain the score matching (SM) objective (see App. B.3):
min θ
Et∈[0,T ]Eut∼pt(ut) [ λ(t)‖sθ(ut, t)−∇vt log pt(ut)‖22 ] (6)
Notice that this objective requires only the velocity gradient of the log-density of the joint distribution, i.e., ∇vt log pt(ut). This is a direct consequence of injecting noise into the velocity variables only. Without loss of generality, pt(ut)=pt(xt,vt)=pt(vt|xt)pt(xt). Hence,
∇vt log pt(ut) = ∇vt [log pt(vt|xt) + log pt(xt)] = ∇vt log pt(vt|xt) (7)
This means that in CLD the neural network-defined score model sθ(ut, t) only needs to learn the score of the conditional distribution pt(vt|xt), an arguably easier task than learning the score of pt(xt), as in previous works, or of the joint pt(ut). This is the case, because our velocity distribution is initialized from a simple Normal distribution, such that pt(vt|xt) is closer to a Normal distribution for all t≥0 (and for any xt) than pt(xt) itself. This is most evident at t=0: The data and velocity distributions are independent at t=0 and the score of p0(v0|x0)=p0(v0) simply corresponds to the score of the Normal distribution p0(v0) from which the velocities are initialized, whereas the score of the data distribution p0(x0) is highly complex and can even be unbounded (Kim et al., 2021). We empirically verify the reduced complexity of the score of pt(vt|xt) in Fig. 2. We find that the score that needs to be learnt by the model is more similar to a score corresponding to a Normal distribution for CLD than for the VPSDE. We also measure the complexity of the neural networks that were learnt to model this score via the squared Frobenius norm of their Jacobians. We find that the CLD-based SGMs have significantly simpler and smoother neural networks than VPSDE-based SGMs for most t, in particular when leveraging a mixed score formulation (see next section).

heading: 3.2 SCALABLE TRAINING
text: A Practical Objective. We cannot train directly with Eq. (6), since we do not have access to the marginal distribution pt(ut). As presented in Sec. 2, we could employ denoising score matching (DSM) and instead sample u0,
and diffuse those samples, which would lead to a tractable objective. However, recall that in CLD the distribution at t=0 is the product of a complex data distribution and a Normal distribution over the initial velocity. Therefore, we propose a hybrid version of score matching (Hyvärinen, 2005) and denoising score matching (Vincent, 2011), which we call hybrid score matching (HSM). In HSM, we draw samples from p0(x0)=pdata(x0) as in DSM, but then diffuse those samples while marginalizing over the full initial velocity distribution p0(v0)=N (v; 0d, γMId) as in regular SM (HSM is discussed in detail in App. C). Since p0(v0) is Normal (and f and G affine), p(ut|x0) is also Normal and this remains tractable. We can write this HSM objective as:
min θ
Et∈[0,T ]Ex0∼p0(x0)Eut∼pt(ut|x0) [ λ(t)‖sθ(ut, t)−∇vt log pt(ut|x0)‖22 ] . (8)
In HSM, the expectation over p0(v0) is essentially solved analytically, while DSM would use a sample-based estimate. In particular, when drawing a sample u0 to diffuse in DSM, we are essentially placing an infinitely sharp Normal with unbounded score (Kim et al., 2021) at u0, which requires undesirable modifications or truncation tricks for stable training (Song et al., 2021c; Vahdat et al., 2021). Hence, with DSM we could lose some benefits of the CLD framework discussed in Sec. 3.1, whereas HSM is tailored to CLD and fundamentally avoids such unbounded scores. Closed form expressions for the perturbation kernel pt(ut|x0) are provided in App. B.1. Score Model Parametrization. (i) Ho et al. (2020) found that it can be beneficial to parameterize the score model to predict the noise that was used in the reparametrized sampling to generate
perturbed samples ut. For CLD, ut = µt(x0) + Lt 2d, where Σt = LtL>t is the Cholesky decomposition of pt(ut|x0)’s covariance matrix, 2d ∼ N ( 2d; 02d, I2d), and µt(x0) is pt(ut|x0)’s mean. Furthermore,∇vt log pt(ut|x0) = −`t d:2d, where d:2d denotes those d components of 2d that actually affect∇vt log pt(ut|x0) (since we take velocity gradients only, not all are relevant).
With Σt = ( Σxxt Σ xv t
Σxvt Σ vv t ) ︸ ︷︷ ︸
“per-dimension” covariance matrix
⊗ Id, we have `t := √
Σxxt Σxxt Σ vv t − (Σxvt )2 . (9)
(ii) Vahdat et al. (2021) showed that it can be beneficial to assume that the diffused marginal distribution is Normal at all times and parametrize the model with a Normal score and a residual “correction”. For CLD, the score is indeed Normal at t = 0 (due to the independently initialized x and v at t=0). Similarly, the target score is close to Normal for large t, as we approach the equilibrium.
Based on (i) and (ii), we parameterize sθ(ut, t) = −`tαθ(ut, t) with αθ(ut, t) = `−1t vt/Σvvt + α′θ(ut, t), where Σ vv t corresponds to the v-v component of the “per-dimension” covariance matrix of the Normal distribution pt(ut|x0 = 0d). In other words, we assumed p0(x0) = δ(x) when defining the analytic term of the score model. Formally, −v/Σvvt is the score of a Normal distribution with covariance Σ̂vvt Id. Following Vahdat et al. (2021), we refer to this parameterization as mixed score parameterization. Alternative model parameterizations are possible, but we leave their exploration to future work. With this definition, the HSM training objective becomes (details in App. B.3):
min θ
Et∈[0,T ]Ex0∼p0(x0)E 2d∼N ( 2d;02d,I2d) [ λ(t)`2t‖ d:2d − αθ(ut=µt(x0)+Lt 2d, t)‖22 ] , (10)
which corresponds to training the model to predict the noise only injected into the velocity during reparametrized sampling of ut, similar to noise prediction in Ho et al. (2020); Song et al. (2021c).
Objective Weightings. For λ(t) = Γβ, the objective corresponds to maximum likelihood learning (Song et al., 2021b) (see App. B.3). Analogously to prior work (Ho et al., 2020; Vahdat et al., 2021; Song et al., 2021b), an objective better suited for high quality image synthesis can be obtained by setting λ(t) = `−2t , which corresponds to “dropping the variance prefactor” ` 2 t .

heading: 3.3 SAMPLING FROM CLD-BASED SGMS
text: To sample from the CLD-based SGM we can either directly simulate the reverse-time diffusion process (Eq. (2)) or, alternatively, solve the corresponding probability flow ODE (Song et al., 2021c;b) (see App. B.5). To simulate the SDE of the reverse-time diffusion process, previous works often relied on Euler-Maruyama (EM) (Kloeden & Platen, 1992) and related methods (Ho et al., 2020; Song et al., 2021c; Jolicoeur-Martineau et al., 2021a). We derive a new solver, tailored to CLD-based models. Here, we provide the high-level ideas and derivations (see App. D for details).
Our generative SDE can be written as (with ūt = uT−t, x̄t = xT−t, v̄t = vT−t):( dx̄t dv̄t ) = ( −M−1v̄t x̄t ) βdt︸ ︷︷ ︸
AH
+
( 0
−ΓM−1v̄t
) βdt+ ( 0√
2Γβdwt ) ︸ ︷︷ ︸
AO
+
( 0
2Γ [ s(ūt, T − t) +M−1v̄t ])βdt︸ ︷︷ ︸ S
(11)
It consists of a Hamiltonian component AH , an Ornstein-Uhlenbeck process AO, and the score model term S. We could use EM to integrate this SDE; however, standard Euler methods are not well-suited for Hamiltonian dynamics (Leimkuhler & Reich, 2005; Neal, 2011). Furthermore, if S was 0, we could solve the SDE in closed form. This suggests the construction of a novel integrator.
We use the Fokker-Planck2 operator formalism (Tuckerman, 2010; Leimkuhler & Matthews, 2013; 2015). Using a similar notation as Leimkuhler & Matthews (2013), the Fokker-Planck equation corresponding to the generative SDE is ∂pt(ūt)/∂t=(L̂∗A+L̂∗S)pt(ūt), where L̂∗A and L̂∗S are the non-commuting Fokker-Planck operators corresponding to the A:=AH+AO and S terms, respectively. Expressions for L̂∗A and L̂∗S are in App. D. We can construct a formal, but intractable solution of the generative SDE as ūt = et(L̂ ∗ A+L̂∗S)ū0, where the operator et(L̂ ∗ A+L̂∗S) (known as the classical
2The Fokker-Planck operator is also known as Kolmogorov operator. If the underlying dynamics is fully Hamiltonian, it corresponds to the Liouville operator (Leimkuhler & Matthews, 2015; Tuckerman, 2010).
propagator in statistical physics) propagates states ū0 for time t according to the dynamics defined by the combined operators L̂∗A + L̂∗S . Although this operation is not analytically tractable, it can serve as starting point to derive a practical integrator. Using the symmetric Trotter theorem or Strang splitting formula as well as the Baker–Campbell–Hausdorff formula (Trotter, 1959; Strang, 1968; Tuckerman, 2010), it can be shown that:
et(L̂ ∗ A+L̂∗S) = lim
N→∞
[ e δt 2 L̂∗AeδtL̂ ∗ Se δt 2 L̂∗A ]N ≈ [ e δt 2 L̂∗AeδtL̂ ∗ Se δt 2 L̂∗A ]N +O(Nδt3) (12)
with time step δt := t/N , where N ∈ N+. The expression suggests that instead of directly evaluating the intractable et(L̂ ∗ A+L̂∗S), we can discretize the dynamics over t into N pieces of step size δt, such that we only need to apply the individual e δt 2 L̂∗A and eδtL̂ ∗ S many times one after another for small time steps δt. A finer discretization results in a smaller error (since N=t/δt the error effectively scales as O(δt2) for fixed t). Hence, this implies an integration method. Indeed, e δt2 L̂∗A ūt is available in closed form, as mentioned before. However, eδtL̂ ∗ S ūt is not. Therefore, we approximate this latter component of the integrator via a standard Euler step. Thus, the integrator formally has an error of the same order as standard EM methods. Nevertheless, as long as the dynamics is not dominated by the S component, our proposed integration scheme is expected to be more accurate than EM, since we split off the analytically tractable part and only use an Euler approximation for the S term. Recall that the model only needs to learn the score of the conditional distribution pt(vt|xt), which is close to Normal for much of the diffusion, in which case the S term will indeed be small. This suggests that the generative SDE dynamics are in fact dominated by AH and AO in practice. Note that only the propagator eδtL̂ ∗ S is computationally expensive, as it involves evaluating the neural network. We coin our novel SDE integrator for CLD-based SGMs Symmetric Splitting CLD Sampler (SSCS). A detailed derivation, analyses, and a formal algorithm are presented in App. D.

heading: 4 RELATED WORK
text: Relations to Statistical Mechanics and Molecular Dynamics. Learning a mapping between a simple, tractable and a complex distribution as in SGMs is inspired by annealed importance sampling (Neal, 2001) and the Jarzynski equality from non-equilibrium statistical mechanics (Jarzynski, 1997a;b; 2011). However, after Sohl-Dickstein et al. (2015), little attention has been given to the origins of SGMs in statistical mechanics. Intuitively, in SGMs the diffusion process is initialized in a non-equilibrium state u0 and we would like to bring the system to equilibrium, i.e., the tractable prior distribution, as quickly and as smoothly as possible to enable efficient denoising. This “equilibration problem” is a much-studied problem in statistical mechanics, particularly in molecular dynamics, where a molecular system is often simulated in thermodynamic equilibrium. Algorithms to quickly and smoothly bring a system to and maintain at equilibrium are known as thermostats. In fact, CLD is inspired by the Langevin thermostat (Bussi & Parrinello, 2007). In molecular dynamics, advanced thermostats are required in particular for “multiscale” systems that show complex behaviors over multiple time- and length-scales. Similar challenges are also encountered when modeling complex data, such as natural images. Hence, the vast literature on thermostats (Andersen, 1980; Nosé, 1984; Hoover, 1985; Martyna et al., 1992; Hünenberger, 2005; Bussi et al., 2007; Ceriotti et al., 2009; 2010; Tuckerman, 2010) may be valuable for the development of future SGMs. Also the framework for developing SSCS is borrowed from statistical mechanics. The same techniques have been used to derive molecular dynamics algorithms (Tuckerman et al., 1992; Bussi & Parrinello, 2007; Ceriotti et al., 2010; Leimkuhler & Matthews, 2013; 2015; Kreis et al., 2017).
Further Related Work. Generative modeling by learning stochastic processes has a long history (Movellan, 2008; Lyu, 2009; Sohl-Dickstein et al., 2011; Bengio et al., 2014; Alain et al., 2016; Goyal et al., 2017; Bordes et al., 2017; Song & Ermon, 2019; Ho et al., 2020). We build on Song et al. (2021c), which introduced the SDE framework for modern SGMs. Nachmani et al. (2021) recently introduced non-Gaussian diffusion processes with different noise distributions. However, the noise is still injected directly into the data, and no improved sampling schemes or training objectives are introduced. Vahdat et al. (2021) proposed LSGM, which is complementary to CLD: we improve the diffusion process itself, whereas LSGM “simplifies the data” by first embedding it into a smooth latent space. LSGM is an overall more complicated framework, as it is trained in two stages and relies on additional encoder and decoder networks. Recently, techniques to accelerate sampling from pre-trained SGMs have been proposed (San-Roman et al., 2021; Watson et al., 2021;
Kong & Ping, 2021; Song et al., 2021a). Importantly, these methods usually do not permit straightforward likelihood estimation (Song et al., 2021a). Therefore, it is unclear whether they suffer from decreased diversity in the few-step limit. Furthermore, they are usually developed for discrete-step diffusion models and not the continuous time framework, which we use.
A complementary work to CLD is “Gotta Go Fast” (GGF) (Jolicoeur-Martineau et al., 2021a), which introduces an adaptive SDE solver for SGMs, tuned towards image synthesis. GGF uses standard Euler-based methods under the hood (Kloeden & Platen, 1992; Roberts, 2012), in contrast to our SSCS that is derived from first principles. Furthermore, our SDE integrator for CLD does not make any data-specific assumptions and performs extremely well even without adaptive step sizes.
Some works study SGMs for maximum likelihood training (Song et al., 2021b; Kingma et al., 2021; Huang et al., 2021). Note that we did not focus on training our models towards high likelihood. Furthermore, Chen et al. (2020) and Huang et al. (2020) recently trained augmented Normalizing Flows, which have conceptual similarities with our velocity augmentation. Methods leveraging auxiliary variables similar to our velocities are also used in statistics—such as Hamiltonian Monte Carlo (Neal, 2011)—and have found applications, for instance, in Bayesian machine learning (Chen et al., 2014; Ding et al., 2014; Shang et al., 2015). As shown in Ma et al. (2019), our velocity is equivalent to momentum in gradient descent and related methods (Polyak, 1964; Kingma & Ba, 2015). Momentum accelerates optimization; the velocity in CLD accelerates mixing in the diffusion process. Lastly, our CLD method can be considered as a second-order Langevin algorithm, but even higher-order schemes are possible (Mou et al., 2021) and could potentially further improve SGMs.

heading: 5 EXPERIMENTS
text: Architectures. We focus on image synthesis and implement CLD-based SGMs using NCSN++ and DDPM++ (Song et al., 2021c) with 6 input channels (for velocity and data) instead of 3.
Relevant Hyperparameters. CLD’s hyperparameters are chosen as β=4, Γ=1 (or equivalently M−1=4) in all experiments. We set the variance scaling of the inital velocity distribution to γ=0.04 and use the proposed HSM objective with the weighting λ(t)=`−2t , which promotes image quality.
Sampling. We generate model samples via: (i) Probability flow using a Runge–Kutta 4(5) method; reverse-time generative SDE sampling using either (ii) EM or (iii) our SSCS. For methods without adaptive stepsize (EM and SSCS), we use evaluation times chosen according to a quadratic function, like previous work (Song et al., 2021a; Kong & Ping, 2021; Watson et al., 2021) (indicated by QS).
Evaluation. We measure image sample quality for CIFAR-10 via Fréchet inception distance (FID) with 50k samples (Heusel et al., 2017). We also evaluate an upper bound on the negative loglikelihood (NLL): − log p(x0)≤−Ev0∼p(v0) log pε(x0,v0)−H , where H is the entropy of p(v0) and log pε(x0,v0) is an unbiased estimate of log p(x0,v0) from the probability flow ODE (Grathwohl et al., 2019; Song et al., 2021c). As in Vahdat et al. (2021), the stochasticity of log pε(x,v) prevents us from performing importance weighted NLL estimation over the velocity distribution (Burda et al., 2015). We also record the number of function—neural network—evaluations (NFEs) during synthesis when comparing sampling methods. All implementation details in App. B.5 and E.

heading: 5.1 IMAGE GENERATION
text: Following Song et al. (2021c), we focus on the widely used CIFAR-10 unconditional image generation benchmark. Our CLD-based SGM achieves an FID of 2.25 based on the probability flow ODE and an FID of 2.23 via simulating the generative SDE (Tab. 1). The only models marginally outperforming CLD are LSGM (Vahdat et al., 2021) and NSCN++/VESDE with 2,000 step predictorcorrector sampling (Song et al., 2021b). However, LSGM uses a model with ≈475M parameters to achieve its high performance, while we obtain our numbers with a model of≈100M parameters. For a fairer comparison, we trained a smaller LSGM also with≈100M parameters, which is reported as “LSGM-100M” in Tab. 1 (details in App. E.2.7). Our model has a significantly lower FID score than “LSGM-100M”. In contrast to NSCN++/VESDE, we achieve extremely strong results with much fewer NFEs (for example see n∈{150, 275, 500} in Tab. 3 and also Tab. 2)—the VESDE performs poorly in this regime. We conclude that when comparing models with similar network capacity and under NFE budgets ≤500, our CLD-based SGM outperforms all published results in terms of FID. We attribute these positive results to our easier score matching task. Furthermore, our model reaches an NLL bound of 3.31, which is on par with recent works such as Nichol & Dhariwal (2021); Austin
Figure 3: CIFAR-10 samples.
Figure 4: CelebA-HQ-256 samples.
et al. (2021); Vahdat et al. (2021) and indicates that our model is not dropping modes. However, our bound is potentially quite loose (see discussion in App. B.5) and the true NLL might be significantly lower. We did not focus on training our models towards high likelihood.
To demonstrate that CLD is also suitable for high-resolution image synthesis, we additionally trained a CLD-SGM on CelebA-HQ-256, but without careful hyperparameter tuning due to limited compute resources. Model samples in Fig. 4 appear diverse and high-quality (additional samples in App. F).

heading: 5.2 SAMPLING SPEED AND SYNTHESIS QUALITY TRADE-OFFS
text: We analyze the sampling speed vs. synthesis quality trade-off for CLD-SGM and study SSCS’s performance under different NFE budgets in Tabs. 2 and 3. We compare to Song et al. (2021c) and use EM to solve the generative SDE for their VPSDE models and PC (reverse-diffusion + Langevin sampler) for VESDE models, which corresponds to their recommended settings3. We also compare to the GGF (Jolicoeur-Martineau et al., 2021a) solver for the generative SDE as well as probability flow sampling with a higher-order adaptive stepsize solver. Further, we compare to LSGM (Vahdat et al., 2021) (using our LSGM-100M model), which also relies on probability flow sampling. We do not compare to few-step sampling methods such as DDIM (Song et al., 2021a), since they do not permit likelihood estimation. This may imply that they trade off synthesis quality for diversity, which makes a fair comparison difficult. With one exception (VESDE with 2,000 NFE) our CLDSGM outperforms all baselines, both when comparing to adaptive and fixed-step size methods.
Several observations stand out: (i) As expected (Sec. 3.3), for CLD, SSCS significantly outperforms EM under limited NFE budgets. When using a fine discretization of the SDE (high NFE), the two perform similarly, which is also expected, as the errors of both methods will become negligible. (ii) In the adaptive solver setting, using a simpler ODE solver, we even outperform GGF, which is tuned towards image synthesis. (iii) Our CLD-SGM also outperforms the LSGM-100M model in terms of FID. It is worth noting, however, that LSGM was designed primarily for faster synthesis, which it achieves by modeling a smooth distribution in latent space instead of the more complex data distribution directly. This suggests that it would be promising to combine LSGM with CLD and train a CLD-based LSGM, to combine the strengths of the two approaches. It would also be
3colab.research.google.com/drive/1dRR_0gNRmfLtPavX2APzUggBuXyjWW55
Table 2: (right) Performance using adaptive stepsize solvers (ODE is based on probability flow, GGF simulates generative SDE). †: taken from Jolicoeur-Martineau et al. (2021a). LSGM corresponds to the small LSGM100M model for fair comparison (denoted as LSGM-100M in Tab. 1, details in App. E.2.7). Error tolerances were chosen to obtain approximately similar NFEs.
FID at n function evaluations ↓ Model Sampler n=50 n=150 n=275 n=500 n=1000 n=2000
CLD EM-QS 52.7 7.00 3.24 2.41 2.27 2.23 CLD SSCS-QS 20.5 3.07 2.38 2.25 2.30 2.29 VPSDE EM-QS 28.2 4.06 2.65 2.47 2.66 2.60
VESDE PC 460 216 11.2 3.75 2.43 2.23†
interesting to develop a more advanced, adaptive SDE solver that leverages SSCS as the backbone and, for example, potentially test our method within a framework like GGF. Our current SSCS only allows for fixed step sizes—nevertheless, it achieves excellent performance.
5.3 ABLATION STUDIES
We perform ablation studies to study CLD’s new hyperparameters (run with a smaller version of our CLD-SGM used above; App. E for details).
Mass Parameter: Tab. 4 shows results for a CLD-SGM trained with different M−1 (also recall that M−1 and Γ are tied together via Γ2 = 4M ; we are always in the critical-damping regime). Different mass values perform mostly similarly. Intuitively, training with smaller M−1 means that noise flows from the velocity variables vt into the data xt more slowly, which necessitates a larger time rescaling β. We found that simply tyingM−1 and β together via β=8 √ M works well and did not further fine-tune.
Initial Velocity Distribution: Tab. 5 shows results for a CLD-SGM trained with different initial velocity variance scalings γ. Varying γ similarly has only a small effect, but small γ seems slightly beneficial for FID, while the NLL bound suffers a bit. Due to our focus on synthesis quality as measued by FID, we opted for small γ. Intuitively, this means that the data at t=0 is “at rest”, and noise flows from the velocity into the data variables only slowly.
Mixed Score: Similar to previous work (Vahdat et al., 2021), we find training with the mixed score (MS) parametrization (Sec. 3.2) beneficial. With MS, we achieve an FID of 3.14, without only 3.56.
Hybrid Score Matching: We also tried training with regular DSM, instead of HSM. However, training often became unstable. As discussed in Sec. 3.2, this is likely because when using standard DSM our CLD would suffer from unbounded scores close to t=0, similar to previous SDEs (Kim et al., 2021). Consequently, we consider our novel HSM a crucial element for training CLD-SGMs.
We conclude that CLD does not come with difficult-to-tune hyperparameters. We expect our chosen hyperparameters to immediately translate to new tasks and models. In fact, we used the same M−1, γ, MS and HSM settings for CIFAR-10 and CelebA-HQ-256 experiments without fine-tuning.

heading: 6 CONCLUSIONS
text: We presented critically-damped Langevin diffusion, a novel diffusion process for training SGMs. CLD diffuses the data in a smoother, easier-to-denoise manner compared to previous SGMs, which results in smoother neural network-parametrized score functions, fast synthesis, and improved expressivity. Our experiments show that CLD outperforms previous SGMs on image synthesis for similar-capacity models and sampling compute budgets, while our novel SSCS is superior to EM in CLD-based SGMs. From a technical perspective, in addition to proposing CLD, we derive CLD’s score matching objective termed as HSM, a variant of denoising score matching suited for CLD, and we derive a tailored SDE integrator for CLD. Inspired by methods used in statistical mechanics, our work provides new insights into SGMs, and implies promising directions for future research.
We believe that CLD can potentially serve as the backbone diffusion process of next generation SGMs. Future work includes using CLD-based SGMs for generative modeling tasks beyond images, combining CLD with techniques for accelerated sampling from SGMs, adapting CLD-based SGMs towards maximum likelihood, and utilizing other thermostating methods from statistical mechanics.

heading: 7 ETHICS AND REPRODUCIBILITY
text: Our paper focuses on fundamental algorithmic advances to improve the generative modeling performance of SGMs. As such, the proposed CLD does not imply immediate ethical concerns. However, we validate CLD on image synthesis benchmarks. Generative modeling of images has promising applications, for example for digital content creation and artistic expression (Bailey, 2020), but can also be used for nefarious purposes (Vaccari & Chadwick, 2020; Mirsky & Lee, 2021; Nguyen et al., 2021). It is worth mentioning that compared to generative adversarial networks (Goodfellow et al., 2014), a very popular class of generative models, SGMs have the promise to model the data more faithfully, without dropping modes and introducing problematic biases. Generally, the ethical impact of our work depends on its application domain and the task at hand.
To aid reproducibility of the results and methods presented in our paper, we will make source code to reproduce all quantitative and qualitative results of the paper publicly available, including detailed instructions. Furthermore, all training details and hyperparameters are already in detail described in the Appendix, in particular in App. E.

heading: A LANGEVIN DYNAMICS
text: Here, we discuss different aspects of Langevin dynamics. Recall the Langevin dynamics, Eq. (5), from the main paper:(
dxt dvt
) = ( M−1vt −xt ) βdt︸ ︷︷ ︸
Hamiltonian component=:H
+
( 0
−ΓM−1vt
) βdt+ ( 0√ 2Γβ ) dwt︸ ︷︷ ︸
Ornstein-Uhlenbeck process=:O
. (13)

heading: A.1 DIFFERENT DAMPING RATIOS
text: As discussed in Sec. 3, Langevin dynamics can be run with different ratios between mass M and squared friction Γ2. To recap from the main paper:
(i) For Γ2 < 4M (underdamped Langevin dynamics), the Hamiltonian component dominates, which implies oscillatory dynamics of xt and vt that slow down convergence to equilibrium.
(ii) For Γ2 > 4M (overdamped Langevin dynamics), the O-term dominates which also slows down convergence, since the accelerating effect by the Hamiltonian component is suppressed due to the strong noise injection.
(iii) For Γ2 = 4M (critical-damping), an ideal balance is achieved and convergence to pEQ(u)occurs as fast as possible in a smooth manner without oscillations.
In Fig. 5, we visualize diffusion trajectories according to Langevin dynamics run in the different damping regimes. We observe that underdamped Langevin dynamics show undesired oscillatory behavior, while overdamped Langevin dynamics perform very inefficiently, too. Critical-damping achieves a good balance between the two and mixes and converges quickly. In fact, it can be shown to be optimal in terms of convergence; see, for example, McCall (2010).
Consequently, we propose to set Γ2 = Γ2critical := 4M in CLD.
A.2 VERY HIGH FRICTION LIMIT AND CONNECTION TO PREVIOUS SDES IN SGMS
Let us re-write the above Langevin dynamics and consider the more general case with timedependent β(t):
dxt = M −1vtβ(t)dt, (14) dvt = − xtβ(t)dt︸ ︷︷ ︸ (ii): potential term −ΓM−1vtβ(t)dt︸ ︷︷ ︸ (iii): friction term + √ 2Γβ(t)dwt︸ ︷︷ ︸ (iv): noise term . (15)
To solve this SDE, let us assume a simple Euler-based integration scheme, with the update equation for a single step at time t (this integration scheme would not be optimal, as discussed in Sec. 3.3., however, it would be accurate for sufficiently small time steps and we just need this to make the connection to previous works like the VPSDE):
xn+1 = xn + β(t)M −1vn+1δt, (16) vn+1 = vn︸︷︷︸ (i): current step velocity − β(t)xnδt︸ ︷︷ ︸ (ii): potential term −β(t)ΓM−1vnδt︸ ︷︷ ︸ (iii): friction term + √ 2β(t)ΓN (0, δt)︸ ︷︷ ︸ (iv): noise term , (17)
Now, let’s assume a friction coefficient Γ = Γmax := Mβ(t)δt . Since the time step δt is usually very small, this correspond to a very high friction. In fact, it can be considered the maximum friction limit, at which the friction is so large that the current step velocity (i) is completely cancelled out by the friction term (iii). We obtain:
xn+1 = xn + β(t)M −1vn+1δt (18) vn+1 = −β(t)xtδt+ √ 2 M
δt N (0, δt). (19)
Now the velocity update, Eq. (19), does not depend on the current step velocity on the right-handside anymore. Hence, we can insert Eq. (19) directly into Eq. (18) and obtain:
xn+1 = xn − β(t)2M−1xnδt2 + √
2β(t)2δtM−1N (0, δt) (20) = xn − β(t)2M−1xnδt2 + √ 2β(t)2δt2M−1N (0, 1). (21)
(22)
Re-defining δt′ := δt2 and β′(t) := β(t)2, we obtain xn+1 = xn − β′(t)M−1xnδt′ + √
2β′(t)δt′M−1N (0, 1), (23) (24)
which corresponds to the high-friction overdamped Langevin dynamics that are frequently run, for example, to train energy-based generative models (Du & Mordatch, 2019; Xiao et al., 2021).
Let’s further absorb the mass M−1 and the time step δt′ into the time rescaling, defining β̂(t) := 2β′(t)M−1δt′. We obtain:
xn+1 = xn − 1
2 β̂(t)xn +
√ β̂(t)N (0, 1) (25)
= (1− 1 2 β̂(t))xn +
√ β̂(t)N (0, 1) (26)
≈ √ 1− β̂(t)xn + √ β̂(t)N (0, 1), (27)
(28)
where the last approximation is true for sufficiently small β̂(t). However, this expression corresponds to
xn+1 ∼ N (xn+1; √
1− β̂(t)xn, β̂(t)) (29) which is exactly the transition kernel of the VPSDE’s Markov chain (Ho et al., 2020; Song et al., 2021c). We see that the VPSDE corresponds to the high-friction limit of a more general Langevin dynamics-based diffusion process of the form of Eq. (13).
If we assume a diffusion as above but with the potential term (ii) set to 0, we can similarly derive the VESDE Song et al. (2021c) as a high-friction limit of the corresponding diffusion. Generally, all previously used diffusions that inject noise directly into the data variables correspond to such high-friction diffusions.
In conclusion, we see that previous high-friction diffusions require an excessive amount of noise to be injected to bring the dynamics to the prior, which intuitively makes denoising harder. For our CLD in the critical damping regime we can run the diffusion for a much shorter time or, equivalently, can inject less noise to converge to the equilibrium, i.e., the prior.

heading: B CRITICALLY-DAMPED LANGEVIN DIFFUSION
text: Here, we present further details about our proposed critically-damped Langevin diffusion (CLD). We provide the derivations and formulas that were not presented in the main paper in the interest of brevity.

heading: B.1 PERTURBATION KERNEL
text: To recap from the main text, in this work we propose to augment the data xt ∈ Rd with auxiliary velocity variables vt ∈ Rd. We then run the following diffusion process in the joint xt-vt-space
dut := ( xt vt ) = f(ut, t) +G(ut, t)dwt (30)
f(ut, t) = f ⊗ Idut, f := ( 0 β(t)M−1 −β(t) −β(t)ΓM−1 ) , (31)
G(ut, t) = G⊗ Id, G := ( 0 0 0 √ 2Γβ(t) ) , (32)
where wt is a standard Wiener process in R2d and β : [0, T ]→ R+0 is a time rescaling 4. In particular, we consider the critically-damped Langevin diffusion which can be obtained by setting M = Γ2/4, resulting in the following drift kernel
fCLD(ut, t) = fCLD ⊗ Idut, fCLD := ( 0 4β(t)Γ−2 −β(t) −4β(t)Γ−1 ) . (33)
Since we only consider the critically-damped cased in this work, we redefine f := fCLD and f := fCLD for simplicity. Since our drift f and diffusion G coefficients are affine, ut is Normally distributed for all t ∈ [0, T ] if ut′ is Normally distributed for some t′ ∈ [0, T ] (Särkkä &
4For our experiments, we only used constant β; however, for generality, we present all derivations for time dependent β(t).
Solin, 2019). In particular, given that ut′ ∼ N (ut′ ;µt′ ,Σt′ = Σt′ ⊗ Id), we derive expressions for µt and Σt = Σt ⊗ Id, the mean and the covariance matrix of ut, respectively. Here Σt′ is a positive-definite diagonal 2-by-2 matrix and Σt is a symmetric positive-definite 2-by-2 matrix.
Following (Section 6.1, Särkkä & Solin, 2019), the mean and covariance matrix of ut follow the following respective ordinary differential equations (ODEs)
dµt dt
= f ⊗ Idµt, (34) dΣt dt = (f ⊗ Id)Σt + [(f ⊗ Id)Σt]> + ( GG> ) ⊗ ID. (35)
Notating µt′ = (xt′ ,vt′)>, the solutions to the above ODEs are
µt = ( 2B(t)Γ−1xt′ + 4B(t)Γ−2vt′ + xt′ −B(t)xt′ − 2B(t)Γ−1vt′ + vt′ ) e−2B(t)Γ −1 , (36)
and
Σt = Σ(t)⊗ Id, (37)
Σ(t) = ( Σxx(t) Σxv(t) Σxv(t) Σvv(t) ) e−4B(t)Γ −1 , (38)
Σxx(t) = Σxxt′ + e 4B(t)Γ−1 − 1 + 4B(t)Γ−1 (Σxxt′ − 1) + 4B2(t)Γ−2 (Σxxt′ − 2) + 16B(t)2Γ−4Σvvt′ ,
(39)
Σxv(t) = −B(t)Σxxt′ + 4B(t)Γ−2Σvvt′ − 2B2(t)Γ1 (Σxxt′ − 2)− 8B2(t)Γ−3Σvvt′ , (40) Σvv(t) = Γ 2
4
( e4B(t)Γ −1 − 1 ) + B(t)Γ + Σvvt′ ( 1 + 4B(t)2Γ−2 − 4B(t)Γ−1 ) + B(t)2 (Σxxt′ − 2) ,
(41) where B(t) = ∫ t t′ β(t̂) dt̂. For constant β(t) = β (as is used in all our experiments), we simply have B(t) = (t − t′)β. The correctness of the proposed mean and covariance matrix can be verified by simply plugging them back in their respective ODEs.
With the above derivations, we can find analytical expression for perturbation kernels p(ut | ·). For example, when conditioning on the the inital data and variance (used in denoising score matching (DSM)), the mean and covariance matrix of the perturbation kernel p(ut | u0) can be obtained by setting µ0 = (x0,v0)>, Σxx0 = 0, and Σ vv 0 = 0.
In our experiments, the initial velocity distribution is set to N (0d, γMId). Conditioning only on the initial data and marginalizing over the full initial velocity distribution (used in hybrid score matching (HSM), see Sec. C), the mean and covariance matrix of the perturbation kernel p(ut | x0) can be obtained by setting µ0 = (x0,0d)>, Σxx0 = 0, and Σ vv 0 = γM .

heading: B.2 CONVERGENCE AND EQUILIBRIUM
text: Our CLD-based training of SGMs—as well as denoising diffusion models more generally—relies on the fact that the diffusion converges towards an analytically tractable equilibrium distribution for sufficiently large t. In fact, from the above equations we can easily see that,
lim t→∞
Σxx(t) = 1, (42)
lim t→∞
Σxv(t) = 0, (43)
lim t→∞
Σvv(t) = Γ2
4 = M, (44)
lim t→∞ µt = 02d, (45)
which establishes pEQ(u) = N (x; 0d, Id)N (v; 0d,MId). Notice that our CLD is an instantiation of the more general Langevin dynamics defined by(
dxt dvt
) = ( M−1vt
∇xt log ppot(xt)
) βdt+ ( 0
−ΓM−1vt
) βdt+ ( 0√ 2Γβ ) dwt. (46)
which has the equilibrium distribution p̂EQ(u) = ppot(x)N (v; 0d,MId) (Leimkuhler & Matthews, 2015; Tuckerman, 2010). However, the perturbation kernel of this Langevin dynamics is not available analytically anymore for arbitrary ppot(x). In our case, however, we have the analytically tractable ppot(x) = N (x; 0d, Id). Note that this correspond to the classical “harmonic oscillator” problem from physics.

heading: B.3 CLD OBJECTIVE
text: To derive the objective for training CLD-based SGMs, we start with a derivation that targets maximum likelihood training in a similar fashion to Song et al. (2021b). Let p0 and q0 be two densities, then
DKL(p0 ‖ q0) = DKL(p0 ‖ q0)−DKL(pT ‖ qT ) +DKL(pT ‖ qT ) (47) = − ∫ T
0
∂DKL(pt ‖ qt) ∂t +DKL(pT ‖ qT ), (48)
where pt and qt are the marginal densities of p0 and q0, respectively, diffused by our criticallydamped Langevin diffusion. As has been shown in Song et al. (2021b), Eq .equation 48 can be written as a mixture of score matching losses. To this end, let us consider the Fokker–Planck equation associated with the the critically-damped Langevin diffusion: ∂pt(ut)
∂t = ∇ut ·
[ 1 2 ( GG> ) ⊗ Id∇upt(ut)− pt(ut)f ⊗ Idut ] (49)
= ∇ut · [hp(ut, t)pt(ut)] , hp(ut, t) := 12 ( GG> ) ⊗ Id∇u log pt(u)− f ⊗ Idut.
(50)
Similarly, we have ∂qt(ut)∂t = ∇ut · [hq(ut, t)qt(ut)]. Assuming log pt(ut) and log qt(ut) are smooth functions with at most polynomial growth at infinity, we have
lim ut→∞ hp(ut, t)pt(ut) = lim ut→∞ hq(ut, t)qt(ut) = 0. (51)
Using the above fact, we can compute the time-derivative of the Kullback–Leibler divergence between pt and qt as ∂DKL(pt ‖ qt)
∂t =
∂
∂t
∫ pt(ut) log pt(ut)
qt(ut) dut (52) = − ∫ pt(ut) [hp(ut, t)− hq(ut, t)]> [∇ut log pt(ut)−∇ut log qt(ut)] dut
(53)
= −1 2
∫ pt(ut) [∇ut log pt(ut)−∇ut log qt(ut)]> ( GG> ) ⊗ Id[∇ut log pt(ut)
−∇ut log qt(ut)] dut (54) = −β(t)Γ ∫ pt(ut)‖∇vt log pt(ut)−∇vt log qt(ut)‖22 dut. (55)
Notice that due to the form ofG, we now have only gradients with respect to the velocity component vt. Combining the above with Eq. (48), we have
DKL(p0 ‖ q0) = ΓEt∼U [0,1],ut∼pt(u) [ β(t)‖∇vt log pt(ut)−∇vt log qt(ut)‖22 ] +DKL(pT ‖ qT )
(56) ≈ Et∼U [0,1],ut∼pt(ut) [ λ(t)‖∇vt log pt(ut)−∇vt log qt(ut)‖22 ] , (57)
where λ(t) = Γβ(t). Note that the approximation holds if DKL(pT ‖ qT ) is insignificantly small. As shown in App. C, the above can be rewritten, up to irrelevant constant terms, as either of the following two objectives:
HSM : Et∼U [0,1],x0∼p0(x0),ut∼pt(ut|x0) [ λ(t)‖∇vt log pt(ut | x0)−∇vt log qt(ut)‖22 ] , (58)
DSM : Et∼U [0,1],u0∼p0(u0),ut∼pt(ut|u0) [ λ(t)‖∇vt log pt(ut | u0)−∇vt log qt(ut)‖22 ] . (59)
For both HSM and DSM, we have shown in App. B.1 that the perturbation kernels pt(ut | x0) and pt(ut | x0) are Normal distributions with the following structure of the covariance matrix:
Σt = Σt ⊗ Id, Σt = ( Σxxt Σ xv t
Σxvt Σ vv t
) . (60)
We can use this fact to compute the gradient∇utpt(ut | ·) ∇utpt(ut | ·) = −∇ut 12 (ut − µt)Σ−1t (ut − µt) (61)
= −Σ−1t (ut − µt) (62) = −L−>t L−1t (ut − µt) (63) = −L−>t 2d, (64)
where 2d ∼ N (0, I2d) and Σt = LtL>t is the Cholesky factorization of the covariance matrix Σt. Note that the structure of Σt implies that Lt = Lt ⊗ Id, where LtL>t is the Cholesky factorization of Σt, i.e,
Lt =
( Lxxt L xv t
Lxvt L vv t
) = √Σxxt 0Σxvt√ Σxxt √ Σxxt Σ vv t −(Σxvt )2 Σxxt  . (65) Furthermore, we have
L−>t = L −> t ⊗ Id (66)
= √Σxxt Σxvt√Σxxt 0 √ Σxxt Σ
vv t −(Σxvt )2 Σxxt −1 ⊗ Id (67) =
 1√Σxxt −Σxzt√Σxxt √Σxxt Σzzt −(Σxvt )2 0 √ Σxxt
Σxxt Σ vv t −(Σxvt )2 ⊗ Id. (68) Using the above, we can compute
∇vtpt(ut | ·) = [∇utpt(ut | ·)]d:2d (69) = [ −L−>t 2d ] d:2d (70)
= −`t d:2d, (71) where
`t :=
√ Σxxt
Σxxt Σ vv t − (Σxvt )2
, (72)
and d:2d denotes those (latter) d components of 2d that actually affect∇vt log pt(ut|·) . Note that `t depends on the conditioning in the perturbation kernel, and therefore `t is different for DSM, which is based on p(ut | u0), and HSM, which is based on p(ut | x0). Therefore, we will henceforth refer to `HSMt and ` DSM t if distinction of the two cases is necessary (otherwise we will simply refer to `t for both).
As discussed in Section 3.2, we model∇vt log qt(ut) as sθ(ut, t) = −`tαθ(ut, t). Plugging everything back in the maximum likelihood objective function, Eq. (58) and Eq. (59), we obtain
HSM : Et∼U [0,1],x0∼p0(x0),ut∼pt(ut|x0) [ λ(t) ( `HSMt )2 ‖ d:2d − αθ(ut, t)‖22] , (73) DSM : Et∼U [0,1],u0∼p0(u0),ut∼pt(ut|u0) [ λ(t) ( `DSMt
)2 ‖ d:2d − αθ(ut, t)‖22] , (74) where ut is sampled via reparameterization:
ut = µt +Lt = µt +
( Lxxt 0:d
Lxvt 0:d + L vv t d:2d
) . (75)
t is analytically bounded. Nevertheless,
numerical computation can be unstable (even when using double precision) in which case adding numerical stabilization of num = 1e−9 suffices to make HSM work.
Note again that Lt is different for HSM and DSM.
Analogously to prior work (Ho et al., 2020; Vahdat et al., 2021; Song et al., 2021b) an objective better suited for high quality image synthesis can be obtained by “dropping the variance prefactor”:
HSMreweighted :Et∼U [0,1],x0∼p0(x0),ut∼pt(ut|x0) [ λHSMreweighted(t) ( `HSMt )2 ‖ d:2d − αθ(ut, t)‖22] = (76)
Et∼U [0,1],x0∼p0(x0),ut∼pt(ut|x0) [ ‖ d:2d − αθ(ut, t)‖22 ] , (77)
DSMreweighted :Et∼U [0,1],u0∼p0(u0),ut∼pt(ut|u0) [ λHSMreweighted(t) ( `DSMt )2 ‖ d:2d − αθ(ut, t)‖22] = (78)
Et∼U [0,1],u0∼p0(u0),ut∼pt(ut|u0) [ ‖ d:2d − αθ(ut, t)‖22 ] , (79)
where we set λHSMreweighted(t) = ( `HSMt )−2 and λDSMreweighted(t) = ( `DSMt )−2 . In App C, we show that HSM is equal to DSM (up to addition of a constant).

heading: B.4 CLD-SPECIFIC IMPLEMENTATION DETAILS
text: Analytically, `HSMt is bounded (in particular, ` HSM 0 = 1/ √ γM ), whereas `DSMt is diverging for t→ 0. In practice, however, we found that computation of `HSMt can also be numerically unstable, even when using double precision. As is common practice for computing Cholesky decompositions, we add a numerical stabilization matrix numI2d to Σt before computing `t. In Figure 6, we visualize `HSMt and ` DSM t for different values of num using our main experimental setup of M = 0.25 and γ = 0.04. Note that very small numerical stabilization of num = 1e−9 in combination with the use of double precision makes HSM work in practice.

heading: B.5 LOWER BOUNDS AND PROBABILITY FLOW ODE
text: Given the score model sθ(ut, t), we can synthesize novel samples via simulating the reverse-time diffusion SDE, Eq. (2) in the main text. This can be achieved, for example, via our novel SSCS, Euler-Maruyama, or methods such as GGF (Jolicoeur-Martineau et al., 2021a). However, Song et al. (2021b;c) have shown that a corresponding ordinary differential equation can be defined that generates samples from the same distribution, in case sθ(ut, t) models the ground truth scores perfectly. This ODE is:
dūt = [ −f(ūt, T − t) + 1
2 G(ūt, T − t)G(ūt, T − t)>∇ūt log pT−t(ūt)
] dt (80)
This ODE is often referred to as the probability flow ODE. We can use it to generate novel data by sampling the prior and solving this ODE, like previous works (Song et al., 2021c). Note that in practice sθ(ut, t) won’t be a perfect model, though, such that the generative models defined by simulating the reverse-time SDE and the probability flow ODE are not exactly equivalent (Song et al., 2021b). Nevertheless, they are very closely connected and it has been shown that their performance is usually very similar or almost the same, when we have learnt a good sθ(ut, t). In addition to sampling the generative SDE in our paper, we also sample from our CLD-based SGMs via this probability flow approach.
With the definition of our CLD, the ODE becomes:( dx̄t dv̄t ) = ( −M−1v̄t x̄t ) βdt︸ ︷︷ ︸
A′H
+
( 0 Γs(ūt, T − t) ) βdt︸ ︷︷ ︸
S′
(81)
Notice the interesting form of this probability ODE for CLD: It corresponds to Hamiltonian dynamics (A′H ) plus the score function term S
′. Compared to the generative SDE (Sec. 3.3), the Ornstein-Uhlenbeck term disappears. Generally, symplectic integrators are best suited for integrating Hamiltonian systems (Neal, 2011; Tuckerman, 2010; Leimkuhler & Reich, 2005). However, our ODE is not perfectly Hamiltonian, due to the score term, and modern non-symplectic methods, such as the higher-order adaptive-step size Runge-Kutta 4(5) integrator (Dormand & Prince, 1980), which we use in practice to solve the probability flow ODE, can also accurately simulate Hamiltonian systems over limited time horizons.
Importantly, the ODE formulation also allows us to estimate the log-likelihood of given test data, as it essentially defines a continuous Normalizing flow (Chen et al., 2018; Grathwohl et al., 2019), that we can easily run in either direction. However, in CLD the input into this ODE is not just the data x0, but also the velocity variable v0. In this case, we can still calculate a lower bound on the log-likelihood:
log p(x0) = log (∫ p(x0,v0)dv0 ) (82)
= log (∫ p(v0) p(x0,v0)
p(v0) dv0
) (83)
≥ Ev0∼p(v0) [log p(x0,v0)− log p(v0)] (84) = Ev0∼p(v0) [log p(x0,v0)] +H(p(v0)) (85)
where H(p(v0)) denotes the entropy of p(v0) (we have H(p(v0)) = 12 log (2πeγM)). We can obtain a stochastic, but unbiased estimate of log p(x0,v0) ≈ log pε(x0,v0) via solving the probability flow ODE with initial conditions (x0,v0) and calculating a stochastic estimate of the logdeterminant of the Jacobian via Hutchinson’s trace estimator (and also calculating the probability of the output under the prior), as done in Normalizing flows (Chen et al., 2018; Grathwohl et al., 2019) and previous works on SGMs (Song et al., 2021c;b). In the main paper, we report the negative of Eq. (85) as our upper bound on the negative log-likelihood (NLL).
Note that this bound can be potentially quite loose. In principle, it would be desirable to perform an importance-weighted estimation of the log-likelihood, as in importance-weighted autoencoders (Burda et al., 2015), using multiple samples from the velocity distribution. However, this isn’t possible, as we only have access to a stochastic estimate log pε(x0,v0). The problems arising
from this are discussed in detail in Appendix F of Vahdat et al. (2021). We could consider training a velocity encoder network, somewhat similar to Chen et al. (2020), to improve our bound, but we leave this for future research.

heading: B.6 ON INTRODUCING A HAMILTONIAN COMPONENT INTO THE DIFFUSION
text: Here, we provide additional high-level intuitions and motivations about adding the Hamiltonian component to the diffusion process, as is done in our CLD.
Let us recall how the data distribution evolves in the forward diffusion process of SGMs: The role of the diffusion is to bring the initial non-equilibrium state quickly towards the equilibrium or prior distribution. Suppose for a moment, we could do so with “pure” Hamiltonian dynamics (no noise injection). In that case, we could generate data from the backward model without learning a score or neural network at all, because Hamiltonian dynamics is analytically invertible (flipping the sign of the velocity, we can just integrate backwards in reverse time direction). Now, this is not possible in practice, since Hamiltonian dynamics alone usually cannot convert the non-equilibrium distribution to the prior distribution. Nevertheless, Hamiltonian dynamics essentially achieves a certain amount of mixing on its own; moreover, since it is deterministic and analytically invertible, this mixing comes at no cost in the sense that we do not have to learn a complex score function to invert the Hamiltonian dynamics. Our thought experiment clearly shows that we should strive for a diffusion process that behaves as deterministically (meaning that deterministic implies easily invertible) as possible with as little noise injection as possible. And this is exactly what is achieved by adding the Hamiltonian component in the overall diffusion process. In fact, recall that it is the diffusion coefficient G of the forward SDE that ultimately scales the score function term of the backward generative SDE (and it is the score function that is hard to approximate with complex neural nets). Therefore, in other words, relying more on a deterministic Hamiltonian component for enhanced mixing (mixing just like in MCMC in that it brings us quickly towards the target distribution, in our case the prior) and less on pure noise injection will lead to a nicer generative SDE that relies less on a score function that requires complex and approximate neural network-based modeling, but more on a simple and analytical Hamiltonian component. Such a SDE could then be solved easier with an appropriate integrator (like our SSCS). In the end, we believe that this is the reason why our networks are “smoother” and why given the same network capacity and limited compute budgets we essentially outperform all previous results in the literature (on CIFAR-10).
We would also like to offer a second perspective, inspired by the MCMC literature. In MCMC, “mixing” helps to quickly traverse the high probability parts of the target distribution and, if an MCMC chain is initialized far from the high probability manifold, to quickly converge to this manifold. However, this is precisely the situation we are in with the forward diffusion process of SGMs: The system is initialized in a far-from-equilibrium state (the data distribution) and we need to traverse the space as efficiently as possible to converge to the equilibrium distribution, this is, the prior. Without efficient mixing, it takes longer to converge to the prior, which also implies a longer generation path in the reverse direction—which intuitively corresponds to a harder problem. Therefore, we believe that ideas from the MCMC literature that accelerate mixing and traversal of state space may be beneficial also for the diffusions in SGMs.

heading: C HSM: HYBRID SCORE MATCHING
text: We begin by recalling our maximum likelihood objective function from App. B.3: Et∼U [0,1] [ λ(t)Eut∼pt(u)[‖∇vt log pt(ut)− sθ(u, t)‖ 2 2] ] , (86)
where sθ(u, t) is our score model. In the following, we dissect the “score matching” part of the above objective:
LSM := Eut∼pt(ut) [ ‖∇vt log pt(ut)− sθ(u, t)‖22 ] (87)
= Eut∼pt(ut)‖sθ(ut, t)‖ 2 2 − 2S(θ) + C2. (88) where C2 := Eut∼pt(ut) [ ‖∇vt log pt(ut)‖22 ] and S(θ) is the cross term discussed below. Following Vincent (2011), we can rewrite LSM as an equivalent (up to addition of a time dependent
constant) denoising score matching objective LDSM: LDSM := Eu0∼p(u0),ut∼pt(ut|u0)‖∇vt log pt(ut | u0)− sθ(ut, t)‖ 2 2 (89)
= LSM + C3 − C2, (90) where C3 := Eu0∼p(u0),ut∼pt(ut|u0) [ ‖∇vt log pt(ut | u0)‖22 ] . Something that might not necessar-
ily be quite obvious is that there is no fundamental need to “denoise” with the distribution p(u0).
“Denoising” only with the data distribution p(x0) and marginalizing over the entire initial velocity distribution p(v0) results in
LHSM := Ex0∼p(x0),ut∼pt(ut|x0)‖∇vt log pt(ut | x0)− sθ(ut, t)‖ 2 2 (91)
= Eut∼pt(ut)‖sθ(ut, t)‖ 2 2 − 2Ex0∼p(x0),ut∼pt(ut|x0)[〈∇vt log pt(ut | x0), sθ(ut, t)〉] + C4,
(92) where C4 := Ex0∼p(x0),ut∼pt(ut|x0) [ ‖∇vt log pt(ut | x0)‖22 ] . In our case, this makes sense since p(v0) is Normal, and therefore (as shown in App B.1), the perturbation kernel pt(ut | x0) is still Normal.
In the following, only for completeness, we redo the derivation of Vincent (2011) and show that LSM is equivalent to LHSM (up to addition of a constant). Starting from S(θ), we have
S(θ) = Eut∼pt(ut) 〈∇vt log pt(ut), sθ(ut, t)〉 (93)
= ∫ ut pt(ut) 〈∇vt log pt(ut), sθ(ut, t)〉 dut (94)
= ∫ ut 〈∇vtpt(ut), sθ(ut, t)〉 dut (95)
= ∫ ut 〈 ∇vt ∫ x0 pt(ut | x0)p0(x0) dx0, sθ(ut, t) 〉 dut (96)
= ∫ ut 〈∫ x0 pt(ut | x0)p0(x0)∇vt log pt(ut | x0) dx0, sθ(ut, t) 〉 dut (97)
= ∫ ut ∫ x0 pt(ut | x0)p0(x0) 〈∇vt log pt(ut | x0), sθ(ut, t)〉 dx0 dut (98)
= Ex0∼p0(x0),ut∼p(ut|x0) [〈∇vt log pt(ut | x0), sθ(ut, t)〉] . (99) where 〈·, ·〉 donates the inner product (notation chosen to be consistent with Vincent (2011)). Hence, we have that
LHSM = LSM + C4 − C2. (100) This further implies that
LHSM = LDSM + C4 − C3. (101) Using the analysis from App B.1, we realize that C3 and C4 can be simplified to d ( `DSMt )2 and
d ( `HSMt )2 . Here, we used the fact that the expected squared norm of a multivariate standard Normal random variable is equal to its dimension, i.e., Eε∼N (0d,Id)‖ε‖ 2 2 = d. This analysis then simplifies Equation (101) to
LHSM = LDSM + d (( `DSMt )2 − (`HSMt )2) . (102) Using this relation, we can also find a connection between the maximum likelihood objective functions. In particular, we have
HSM = Et∼U [0,1] [λ(t)LHSM] (103) = Et∼U [0,1] [λ(t)LDSM] + dEt∼U [0,1] [ λ(t) (( `DSMt )2 − (`HSMt )2)] (104) = DSM + dΓEt∈U [0,1] [ β(t) (( `DSMt )2 − (`HSMt )2)] , (105)
where HSM and DSM are defined in Equation (73) and in Equation (74), respectively. Furthermore, we have the following relation for our reweighted objectives
HSMreweighted = Et∼U [0,1] [ λHSMreweighted(t)LHSM ] (106)
= Et∼U [0,1] [( `HSMt )−2 LHSM] (107) = Et∼U [0,1] [( `HSMt
)−2 LDSM]+ dEt∼U [0,1] [(`HSMt )−2 ((`DSMt )2 − (`HSMt )2)] (108)
= Et∼U [0,1] [( `HSMt )−2 LDSM]+ dEt∼U [0,1] [(`HSMt )−2 (`DSMt )2 − 1] . (109)

heading: C.1 GRADIENT VARIANCE REDUCTION VIA HSM
text: Above, we derived that LHSM = LDSM + C, so one might wonder why we advocate for HSM over DSM. In practice, we never solve expectations analytically but rather approximate them using Monte Carlo estimates. In the remainder of this section, we will show that in practice (Monte Carlo) gradients based on HSM have lower variance than those based on DSM.
From Eq. (103) and Eq. (104), we have
∇θHSM = Et∈U [0,1] [λ(t)∇θLHSM] (110) = Et∈U [0,1] [λ(t)∇θLDSM] , (111)
where θ are the learnable parameters of the neural network. Instead of comparing the above expectations directly, we instead compare λ(t)∇θLHSM with λ(t)∇θLDSM for t ∈ [0, 1] at discretized time values (as is done in practice). Replacing LHSM and LDSM with a single Monte Carlo estimate (as is used in practice), we have
λ(t)∇θLHSM ≈ λ(t)∇θsθ(ut, t)∇sθ(ut,t)‖∇vt log pt(ut | x0)− sθ(ut, t)‖ 2 2, (112) x0 ∼ p(x0),ut ∼ pt(ut | x0), (113) λ(t)∇θLDSM ≈ λ(t)∇θsθ(ut, t)∇sθ(ut,t)‖∇vt log pt(ut | u0)− sθ(ut, t)‖ 2 2, (114)
u0 ∼ p(u0),ut ∼ pt(ut | u0). (115) where we applied the chain-rule. In both Eq. (112) and Eq. (114), ut is sampled from the same distribution. Hence, λ(t)∇θsθ(ut, t) acts as a common scaling factor, with the variance difference between HSM and DSM originating from the squared norm term. Hence, we ignore λ(t)∇θsθ(ut, t) and only focus our analysis on the gradient of the norm terms, which we can further simplify:
(HSM) 1
2 ∇sθ(ut,t)‖∇vt log pt(ut | x0)− sθ(ut, t)‖ 2 2 = sθ(ut, t)− log pt(ut | x0) =: KHSM,
(116)
and
(DSM) 1
2 ∇sθ(ut,t)‖∇vt log pt(ut | x0)− sθ(ut, t)‖ 2 2 = sθ(ut, t)− log pt(ut | u0) =: KDSM.
(117)
We explore this difference in a realistic setup; in particular, we evaluate KHSM and KDSM for all data points in the CIFAR-10 training set. We choose sθ to be our trained ablation CLD model (with the standard setup of M−1 = β = 4, see Sec. E.2.1 for model details).
We then use these samples to estimate the covariance matrices CovHSM and CovDSM of the random variables KHSM and KDSM, respectively. As is common practice in statistics, we consider only the trace of the estimated covariance matrices 5. The trace of the covariance matrix (of a random variable) is also commonly referred to as the total variation (of a random variable).
We visualize our results in Fig. 7. For HSM, there is barely any visual difference in Tr(Cov) for γ = 0.04 and γ = 1. For DSM, both γ = 0.04 and γ = 1 result in very large Tr(Cov) values for
5Arguably, the most prominent algorithm that follows this practice is principal component analysis (PCA).
small t. For large t, Tr(Cov) is considerably smaller for γ = 0.04 than for γ = 1. However, in practice, we found that DSM is even unstable for small γ. Given this analysis, we believe this is due to the large gradient variance for small t. In conclusion, these results demonstrate a clear variance reduction by the HSM objective, in particular for large γ.
Note that re-introducing λ(t) weightings would allow us to scale the Tr(Cov) curves according to the “reweighted” objective or the maximum likelihood objective. However, we believe it is most instructive to directly analyze the gradient of the norm term itself.

heading: D SYMMETRIC SPLITTING CLD SAMPLER (SSCS)
text: In this section, we present a more complete derivation and analysis of our novel Symmetric Splitting CLD Sampler (SSCS).

heading: D.1 BACKGROUND
text: Our derivation is inspired by methods from the statistical mechanics and molecular dynamics literature. In particular, we are leveraging symmetric splitting techniques as well as (Fokker–Planck) operator concepts. The high-level idea of symmetric splitting as well as the operator formalism are well-explained in Tuckerman (2010), in particular in Section 3.10, including simple examples. Symmetric splitting methods for stochastic dynamics in particular are discussed in detail in Leimkuhler & Matthews (2015). We also recommend Leimkuhler & Matthews (2013), which discusses splitting methods for Langevin dynamics in a concise but insightful manner.

heading: D.2 DERIVATION AND ANALYSIS
text: Generative SDE. From Sec. 3.3, recall that our generative SDE can be written as (with ūt = uT−t, x̄t = xT−t, v̄t = vT−t):( dx̄t dv̄t ) = ( −M−1v̄t x̄t ) βdt︸ ︷︷ ︸
AH
+
( 0
−ΓM−1v̄t
) βdt+ ( 0√
2Γβdwt ) ︸ ︷︷ ︸
AO
+
( 0
2Γ [ s(ūt, T − t) +M−1v̄t ])βdt︸ ︷︷ ︸ S
Fokker–Planck Equation and Fokker–Planck Operators. The evolution of the probability distribution pT−t(ūt) is described by the general Fokker Planck equation (Särkkä & Solin, 2019):
∂pT−t(ūt) ∂t = − 2d∑ i=1 ∂ ∂ūi [µi(ūt, T − t)pT−t(ūt)] + 2d∑ i=1 2d∑ j=1 ∂2 ∂ūi∂ūj [Dij(ūt, T − t)pT−t(ūt)]
(118)
with µ(ūt, T − t) = ( −M−1v̄t
x̄t
) β + ( 0
−ΓM−1v̄t
) β + ( 0
2Γ [ s(ūt, T − t) +M−1v̄t ])β (119) D(ūt, T − t) = ( 0 0 0 2Γβ ) ⊗ Id (120)
For our SDE, we can write the Fokker Planck equation in short form as
∂pT−t(ūt) ∂t = (L̂∗A+L̂∗S)pT−t(ūt) (121)
with the Fokker–Planck operators (defined via their action on functions of the variables φ(ūt)):
L̂∗Aφ(ūt) := −βM−1v̄t∇x̄tφ(ūt) + βx̄t∇v̄tφ(ūt)− ΓM−1∇v̄t [v̄tφ(ūt)] + 2Γβ∆v̄tφ(ūt)
(122) L̂∗Sφ(ūt) := 2Γβ [ ∇v̄t ( [s(ūt, T − t) +M−1v̄t ] φ(ūt) ) (123)
∆v̄t := d∑ i=1 ( ∂2 ∂x̄2i + ∂2 ∂v̄2i ) (124)
We are providing these formulas for transparency and completeness. We do not directly leverage them. However, working with these operators can be convenient. In particular, the operators describe the time evolution of states ūt under the stochastic dynamics defined by the SDE. Given an initial state ū0, we can construct a formal solution to the generative SDE via (Tuckerman, 2010; Leimkuhler & Matthews, 2015):
ūt = e t(L̂∗A+L̂∗S)ū0 (125)
where the operator et(L̂ ∗ Aφ+L̂∗Sφ) is known as the classical propagator that propagates states ū0 for time t according to the dynamics defined by the combined Fokker–Planck operators L̂∗A + L̂∗S (to avoid confusion, note that in Eq. (125) the operator et(L̂ ∗ Aφ+L̂∗Sφ) is applied on ū0 in an elementwise or “vectorized” fashion on all elements of ū0 in parallel). The problem with that expression is that we cannot analytically evaluate it. However, we can leverage it to design an integration method.
Symmetric Splitting Integration. Using the symmetric Trotter theorem or Strang splitting formula as well as the Baker–Campbell–Hausdorff formula (Trotter, 1959; Strang, 1968; Tuckerman, 2010), it can be shown that:
et(L̂ ∗ A+L̂∗S) = lim
N→∞
[ e δt 2 L̂∗AeδtL̂ ∗ Se δt 2 L̂∗A ]N ≈ [ e δt 2 L̂∗AeδtL̂ ∗ Se δt 2 L̂∗A ]N +O(Nδt3) (126)
with time step δt := t/N , where N ∈ N+. The expression suggests that instead of directly evaluating the intractable et(L̂ ∗ A+L̂∗S), we can discretize the dynamics over t into N pieces of step size δt, such that we only need to apply the individual e δt 2 L̂∗A and eδtL̂ ∗ S many times one after another for small time steps δt. A finer discretization implies a smaller error (sinceN=t/δt the error effectively scales as O(δt2) for fixed t). Hence, this implies an integration method. The general idea of such splitting schemes is to split an initially intractable propagator into separate terms, each of which is analytically tractable. In that case, the overall integration error for many steps is only due to the splitting scheme error6, but not due to the evaluation of the individual updates. Such techniques are,
6In principle, the error of the splitting scheme is defined more specifically by the commutator of the noncommuting Fokker–Planck operators. See, for example Leimkuhler & Matthews (2013; 2015); Tuckerman (2010).
for example, popular in molecular dynamics to develop symplectic integrators as well as accurate samplers (Tuckerman et al., 1992; Tuckerman, 2010; Leimkuhler & Matthews, 2013; 2015; Bussi & Parrinello, 2007).
Analyzing the Splitting Terms. Next, we need to analyze the two individual terms:
(i) Let us first analyze e δt2 L̂ ∗ A ūt: This term describes the stochastic evolution of ūt under the dynamics of a SDE like Eq. (118), but with S set to zero. However, if S is set to zero, the remaining SDE has affine drift and diffusion coefficients. In that case, if the input is Normal (or a discrete state corresponding to a Normal with 0 variance) then the distribution is Normal at all times and we can calculate the evolution analytically. In particular, we can solve the differential equations for the mean µ̄ δt
2 and covariance Σ̄ δtn 2 of the Normal (these differential equation are, for example, pre-
sented in Särkkä & Solin (2019), also see Sec. B.1. The result below can be validated by inserting it into the differential equations—similar to Sec. B.1), and obtain
µ̄ δt 2 (ūt) =
( 2β δt2 Γ
−1x̄t − 4β δt2 Γ−2v̄t + x̄t β δt2 x̄t − 2β δt2 Γ−1v̄t + v̄t
) e−2β δt 2 Γ −1 , (127)
as well as
Σ̄ δt 2 = Σ̄
( δt
2
) ⊗ Id, (128)
Σ̄
( δt
2
) = ( Σ̄xx( δt2 ) Σ̄
xv( δt2 ) Σ̄xz( δt2 ) Σ̄ vv( δt2 )
) e−4β δt 2 Γ −1 , (129)
Σ̄xx ( δt
2
) = e4β δt 2 Γ −1 − 1− 4β δt
2 Γ−1 − 8
( β δt
2
)2 Γ−2, (130)
Σ̄xv ( δt
2
) = 4β δt
2 Γ−2, (131)
Σ̄vv ( δt
2
) = Γ2
4
( e4β δt 2 Γ −1 − 1 ) + β δt
2 Γ− 2
( β δt
2
)2 . (132)
Hence, we can write the action of the the propagator e δt 2 L̂∗A on a state ūt as:
e δt 2 L̂∗A ūt ∼ N (ūt+ δt2 ; µ̄ δt2 (ūt), Σ̄ δt2 ). (133)
(ii): Next, we need to analyze eδtL̂ ∗ S ūt. Unfortunately, we cannot calculate the action of the propagator eδtL̂ ∗ S on ūt analytically and we need to make an approximation. From Eq. (118), we can easily see that the propagator eδtL̂ ∗ S describes the evolution of the velocity component v̄t for time step δt under the ODE: dv̄t = 2βΓ [ s(ūt, T − t) +M−1v̄t ] dt (134)
We propose to simply solve this ODE for the step δt via a simple step of Euler’s method, resulting in:
eδtL̂ ∗ S ūt ≈ ūt + δt
( 0
2βΓ [ s(ūt, T − t) +M−1v̄t ])+O(δt2) (135) = eδtL̂ ∗ Euler S ūt +O(δt2), (136)
with the informal definition
eδtL̂ ∗ Euler S ūt := ūt + δt
( 0
2βΓ [ s(ūt, T − t) +M−1v̄t ] .
) (137)
Error Analysis. It is now instructive to study the overall error of our proposed integrator. With the additional Euler integration in one of the splitting terms, we have
et(L̂ ∗ A+L̂∗S) ≈ [ e δt 2 L̂∗A ( eδtL̂ ∗ Euler S +O(δt2) ) e δt 2 L̂∗A ]N +O(Nδt3) (138)
= [ e δt 2 L̂∗A ( eδtL̂ ∗ Euler S ) e δt 2 L̂∗A ]N +NO(δt2) (139)
= [ e δt 2 L̂∗A ( eδtL̂ ∗ Euler S ) e δt 2 L̂∗A ]N +O(δt), (140)
where we used N = tδt and only kept the highest error terms. We see that, just like Euler’s method, also our SSCS is a first-order integrator with local error ∼δt2 and global error ∼δt. This is of course expected, considering that we used an Euler step for the S term. Nevertheless, as long as the dynamics is not dominated by the S component, our proposed integration scheme is expected to still be more accurate than EM, since we split off the analytically tractable part and only use an Euler approximation for the S term.
To this end, recall that the model only needs to learn the score of the conditional distribution pt(vt|xt), which is close to Normal for much of the diffusion, in which case the S term will indeed be small. This suggests that the generative SDE dynamics are in fact dominated by AH and AO in practice. From another perspective, note that (recalling that sθ(ut, t) = −`tαθ(ut, t) with αθ(ut, t) = ` −1 t vt/Σ vv t + α ′ θ(ut, t) from Sec. 3.2):
s(ūt, T − t) +M−1v̄t = −`tα′θ(ut, t)− vt Σvvt +M−1v̄t (141)
= −`tα′θ(ut, t) + v̄t ( 1
M − 1
Σvvt
) . (142)
For large parts of the diffusion, Σvvt is indeed close to 1/M , such that the v̄t term is very small (this cancellation is the reason why we pulled the M−1v̄t term into the S component). In Sec. 3, we have also seen that our neural network component can be much smoother than that of previous SGMs. Overall, this suggests that the error of SSCS indeed might be smaller than the error we would obtain when applying a naive Euler–Maruyama integrator to the full generative SDE. Our positive experimental results in Sec. 5.2 validate that. Only in the limit for very small steps, both our SSCS and EM make only very small errors and are expected to perform equally well, which is exactly what we observe in our experiments. Our SSCS turns out to be well-suited for integrating the generative SDE with relatively few synthesis steps.
Note that in practice we do not use SSCS to integrate all the way from t = 0 to t = T , but only up to t = T − , and perform a denoising step, similar to previous works (Jolicoeur-Martineau et al., 2021a; Song et al., 2021c). It is worth noting that our SSCS scheme would also be applicable when we used time-dependent β(t), as in our more general derivation of the CLD perturbation kernel in App. B. However, since we only used constant β in the main paper, we also presented SSCS in that way.
A promising direction for future research would be to extend our SSCS towards adaptive step sizes and to use techniques to facilitate a higher-order integration, while still leveraging the unique advantages of SSCS.
SSCS Algorithm. Finally, we summarize SSCS in terms of a concise algorithm:
Algorithm 1 Symmetric Splitting CLD Sampler (SSCS) Input: Score function sθ(ūt, T − t), CLD parameters Γ, β, M = Γ2/4, number of sampling steps N , step sizes {δtn ≥ 0}N−1n=0 chosen such that := T − ∑N−1 n=0 δtn ≥ 0 (stepsizes can vary, for example in QS).
Output: Synthesized model sample x̄′N , along with a velocity sample v̄′N
x̄0 ∼ N (x̄0; 0d, Id), v̄0 ∼ N (v̄0; 0d,MId), ū0 = (x̄0, v̄0) . Draw initial prior samples from pEQ(u) t = 0 . Initialize time for n = 0 to N − 1 do
ūn+ 1 2 ∼ N (ūn+ 1 2 ; µ̄ δtn 2 (ūn), Σ̄ δtn 2 ) . First half-step: apply exp{ δtn 2 L̂∗A} on ūn
ū′ n+ 1 2 ← ūn+ 1 2 + δtn
( 0
2βΓ [ s(ūt, T − t) +M−1v̄t ]) . Full step: apply exp{δtnL̂∗S} on ūn+ 1 2
ūn+1 ∼ N (ūn+1; µ̄ δtn 2 (ū′ n+ 1 2 ), Σ̄ δtn 2 ) . Second half-step: apply exp{ δtn 2 L̂∗A} on ū′n+ 1 2
t← t+ δtn . Update time end for ū′N ← ūN − ( 0 βM−1 −β −ΓβM−1 ) ⊗ Id ūN + ( 0 2βΓs(ūt, ) ) . Denoising (x̄′N , v̄ ′ N ) = ū ′ N . Extract output data and velocity samples
Note that the algorithm uses the expressions in Eqs. (127) and (128) for µ̄t and Σ̄t. Furthermore, in practice in the denoising step at the end, we usually only update the x̄′N component of ū ′ N , since we are only interested in the data sample. This saves us the final neural network call during denoising, which only affects the v̄′N component. However, we did write the algorithm in the general way, which also allows to correctly generate the velocity sample v̄′N .
Also note that we could combine the second half-step from one iteration of SSCS with the first half-step from the next iteration of SSCS. This is commonly done in the Leapfrog integrator, which follows a similar structure as our SSCS. However, it is not important in our case, as the only computational costly operation is in the center full step, which involves the neural network evaluation. The first and last half-step come at virtually no computational cost.
E IMPLEMENTATION AND EXPERIMENT DETAILS

heading: E.1 SCORE AND JACOBIAN EXPERIMENT
text: In this section, we provide details for the experiments presented in Sec. 3.1. For both experiments, we consider a two-dimensional simple mixture of Normals of the form
pdata(x) = 9∑ k=1 1 9 p(k)(x), (143)
where p(k)(x) = N (x;µk; 0.042I2) and
µ1 = ( −a 0 ) , µ2 = ( −a/2 a/2 ) , µ3 = ( 0 a ) ,
µ4 = ( −a/2 −a/2 ) , µ5 = ( 0 0 ) , µ6 = ( a/2 a/2 ) ,
µ7 = ( 0 −a ) , µ8 = ( a/2 −a/2 ) , µ9 = ( a 0 ) ,
and a = 2− 1 2 . The choice of this data distribution is not arbitrary. In fact, mixture of Normal distributions are diffused by simply diffusing the components, i.e., setting p0(x0) = pdata(x), we have
pt(xt) = 9∑ k=1 1 9 p (k) t (xt), (144)
where p(k)t are the diffused components (analogously for CLD with velocity augmentation). This means that for both CLD as well as VPSDE Song et al. (2021c) we can diffuse pdata(x) analytically. This allows us to perform interesting analyses that would be possible when working, for example, with image data. We visualize the data distribution in Fig. 9.
Score experiment: We empirically verify the reduced complexity of the score of pt(vt|xt), which is learned in CLD, compared to the score of pt(xt), which is learned in VPSDE. To avoid scaling issues between VPSDE and CLD, we chose M = γ = 1 for CLD in this experiment; this results in an equilibrium distribution ofN (0, 1) (for both data and velocity component, which are independent at equilibrium), which is the same as the equilibrium distribution of VPSDE. We then measure the difference of the respective scores at time t and the equilibrium (or prior) scores, i.e. (recall that the score of a Normal distribution p(x) = N (x; 0, 1) is simply∇x log p(x) = −x),
ξVPSDE(t) := Ext∼p(xt)‖∇xt log pt(xt) + xt‖ 2 2, (145)
ξCLD(t) := Eut∼p(ut)‖∇vt log pt(vt | xt) + vt‖ 2 2. (146)
The expectations are approximated using 105 samples from p(xt) and p(ut) for VPSDE and CLD, respectively. As can be seen in Fig. 8a, ξCLD(t) is smaller than ξVPSDE(t) for all t ∈ [0, T ]. The difference is particularly striking for small time values t. Other previous SDEs, such as the VESDE, sub-VPSDE, etc., are expected to behave similarly. This result implies that the ground truth scores that need to be learnt in CLD are closer to Normal scores than the ground truth scores in previous SDEs like the VPSDE. Since the score of a Normal is very simple—and indeed directly leveraged in our mixed score formulation—we would intuitively expect that the CLD training task is easier.
Complexity experiment: Therefore, to understand the above observations in terms of learning neural networks, we train a small ResNet architecture (less than 100k parameters) for each of the following four setups: both CLD and VPSDE each with and without a mixed score parameterization. The mixed score of the VPSDE simply assumes a standard normal data distribution (which is also the equilibrium distribution of VPSDE) resulting in adding −xt to the score function. Formally, −xt is the score of a Normal distribution with unit variance. We train the models for 1M iterations using fresh data synthesized from pdata at a batch size of 512. The model and data distributions are visualized in Fig. 9. We measure the complexity of the trained neural networks using the squared Frobenius norm of the networks’ Jacobians. For CLD, we have
J CLDF (t) := Eut∼p(ut)‖∇utα′θ(ut)‖ 2 F (147)
Similarly, for VPSDE, we compute
J VPSDEF (t) := Ext∼p(xt)‖∇xtα′θ(xt)‖ 2 F (148)
For both CLD and VPSDE, expectations are again approximated using 105 samples. As can be seen in Fig. 8b the neural network complexity is lower for CLD compared to VPSDE. A mixed score formulation helps decreasing the neural network complexity for both CLD and VPSDE. This result implies that the arguably simpler training task in CLD indeed also translates to reduced model complexity in that the neural network is smoother as measured by JF(t). In large-scale experiments, this would mean that, given similar model capacity, a CLD-based SGM could potentially have a higher expressivity. Or, on the other hand, similar performance could be achieved with a smoother and potentially smaller model. Indeed these findings are in line with our strong results on the CIFAR10 benchmark.
E.2 IMAGE MODELING EXPERIMENTS
We perform image modeling experiments on CIFAR-10 as well as CelebA-HQ-256. We report FID scores on CIFAR-10 for our main model in for a sweep of different solvers; see Tab. 3 and
Tab. 2. We further present generated samples for both models in Sec. 5 using Euler–Maruyama with 2000 quadratic striding steps and Runge–Kutta 4(5) for CIFAR10 and CelebA-HQ-256, respectively. We present additional samples for variety of solver settings in App. F. All (average) NFEs for the Runge–Kutta solver are computed using a batch size of 128.

heading: E.2.1 TRAINING DETAILS AND MODEL ARCHITECTURES
text: Our models are based on the NCSN++ and the DDPM++ architectures from Song et al. (2021c). Importantly, we changed the number of input channels from three to six to facilitate the additional velocity variables. Note that the number of additional parameters due to this change is negligible.
For fair comparison, we train our models using the same t-sampling cutoff during training as is used for VESDE and VPSDE in Song et al. (2021c). Note, however, that this is not strictly necessary for CLD as we do not have any “blow-up” of the SDE due to unbounded scores as t→ 0. We summarize our three model architectures as well as our SDE and training setups in Tab. 6.

heading: E.2.2 CIFAR10 RESULTS FOR VESDE AND VPSDE
text: The results reported for VESDE and VPSDE using the GGF sampler are taken from JolicoeurMartineau et al. (2021a). All other results for VESDE and VPSDE are generated by using the provided PyTorch code as well as the provided checkpoints from Song et al. (2021c).

heading: E.2.3 QUADRATIC STRIDING
text: When we simulate our generative SDE numerically, using for example EM or our SSCS, we need to choose a time discretization. Given a certain NFE budget NNFE, how do we choose time step sizes?
The standard approach is to use an equidistant discretization, corresponding to a set of evaluation time steps {δti ≥ 0}NNFE−1i=0 with δti = 1NNFE ∀i ∈ [0, NNFE − 1]. However, prior work (Song et al., 2021a; Kong & Ping, 2021; Watson et al., 2021) has shown that it can be beneficial to focus function evaluation on times t close to the “data”. This is because the diffusion process distribution is most complex close to the data and almost perfectly Normal close to the prior. Among other techniques, these works used a useful heuristic, denoted as quadratic striding (QS), which discretizes the integration interval such that the evaluation times follow a quadratic schedule and the individual time steps follow a linear schedule. We also used this QS approach in our experiments.
We can formally define it as follows (assuming time interval [0.0, 1.0] here for simplicity): Denote the evaluation times as τi (including 0.0 and 1.0) and define:
τi = cτ i 2 ∀i ∈ [0, NNFE]. (149)
Hence, δti = τi − τi−1 = cτ (2i− 1) ∀i ∈ [1, NNFE], (150)
and cτ = 1N2NFE to ensure that τNFE = 1.0.
This describes the time steps as going from t = 0 to t = 1. During synthesis, however, we are going backwards. Hence, we can define our time steps as
δtj = cτ [2NNFE − 2j + 1] ∀j ∈ [1, NNFE], (151) where j now counts time steps in the other direction. Note that this can be easily adapted to general integration intervals [ , T ].

heading: E.2.4 DENOISING
text: As has been pointed out in Jolicoeur-Martineau et al. (2021b), samples that are generated with models similar to ours can contain noise that is hard to detect visually but worsens FID scores significantly.
Denoising Formulas. For fair comparison we use the same denoising setup for all experiments we conducted (including VESDE (PC/ODE) and VPSDE (EM/ODE)) except for LSGM 7. We simulate the underlying generative ODE/SDE until the time cutoff ε = 1e−3 and then take a single denoising step of the form
u0 = uε − f(uε, ε)ε+G(uε, ε)G(uε, ε)> (
0d sθ(uε, ε)
) ε. (152)
This denoising step can be considered as an Euler–Maruyama step without noise injection. For SDEs acting on data directly (VESDE, VPSDE, etc.) the corresponds denoising formula is
x0 = xε − f(xε, ε)ε+G(xε, ε)G(xε, ε)>sθ(xε, ε)ε. (153)
Influence of Denoising on Results. For SDEs acting in the data space directly, it has been reported that this denoising step is crucial to obtain good FID numbers Jolicoeur-Martineau et al. (2021b); Song et al. (2021c). When we simulate the generative probability flow ODE we found that denoising is important in order for the Runge–Kutta solver not to “blow-up” as t→ 0. On the other hand, when simulating CLD using our new SSCS solver, we found that denoising only slightly influences FID (see Tab 7). We believe that this might be because the neural network does not have any influence on the denoising step for CLD. Putting the drift and diffusion coefficients of CLD in the denoising formula in Eq. (152), we obtain
u0 = uε − ( 0 βM−1 −β −βΓM−1 ) ⊗ Iduεε+ ( 0d 2Γβ ) sθ(uε, ε)ε =⇒ x0 = xε − βM−1vεε.
(154)
This means the denoising for CLD results in u0 being a linear transformation of uε.
7Denosing has not been used in the original LSGM work (Vahdat et al., 2021) and is not needed in their case, since the output of the latent SGM lives in a smooth latent space and is further processed by a decoder.

heading: E.2.5 SOLVER TOLERANCES FOR RUNGE–KUTTA 4(5)
text: In Tab. 2, we report FID numbers for a Runge–Kutta 4(5) solvers (Dormand & Prince, 1980) as well as the “Gotta Go Fast” solver from Jolicoeur-Martineau et al. (2021a) (see their Table 1). For simulating CLD with Runge–Kutta 4(5) we chose the solver tolerances to hit certain regimes of NFEs to facilitate comparison with VPSDE and VESDE. We obtain a mean number of function evaluations of 312 and 137 using Runge–Kutta 4(5) solver tolerances of 1e−5 and 1e−3, respectively. For VESDE, VPSDE and LSGM we used 1e−5 as the ODE solver tolerance.

heading: E.2.6 ABLATIONS EXPERIMENTS
text: The model architecture used for all ablation experiments can be found in Tab. 6. As pointed out in Sec. 5 we found that the hyperparameters γ and M only have small effect on CIFAR-10 FID numbers. On the other hand, we found that the mixed score parameterization helps significantly in obtaining competitive FIDs.
E.2.7 LSGM-100M MODEL
Our CLD-based SGM has≈108M parameters, while the original CIFAR-10 Latent SGM from Vahdat et al. (2021), to which we compare in Tab. 1, uses≈476M parameters. To establish a fairer comparison between our CLD-based SGMs and LSGM (Vahdat et al., 2021), we train another smaller LSGM model with ≈109M parameters. To do this, we followed the exact setup of the “CIFAR10 (balanced)” model from LSGM (see Table 7 in Vahdat et al. (2021)), with a few minor modifications: We used a VAE backbone model with only 10 groups instead of 20, which corresponds to a reduction in parameters by a factor of 2 in the encoder and decoder networks. We also reduced the convolutional channels in the latent space SGM from 512 to 256 and reduced the number of the residual cells per scale from 8 to 4. With these modifications the resulting “LSGM-100M” uses only ≈109M parameters overall with approximately half of them in the encoder and decoder networks and the other half in the latent SGM. Other than these architecture modifications, our model is trained in exactly the same way as the bigger, original models in Vahdat et al. (2021).
For evaluation, we follow the recommended setting by Vahdat et al. (2021) and use the same RungeKutta 4(5) ODE solver with an error tolerance of 1e−5 to solve the probability flow ODE in LSGM’s latent space. LSGM-100M achieves an FID of 4.60, an NLL bound of 2.96 bpd, and requires on average 131 NFE for sampling new images. We report these results in Tabs. 1 and 2 in the main text.
Note that we also tried training a model following the “CIFAR10 (best FID)” setup, but found training to be unstable (however the orignal “CIFAR10 (best FID)” model from Vahdat et al. (2021) only performs marginally better in FID than their “CIFAR10 (balanced)” model anyway). Furthermore, we also tried training another small LSGM with a similar number of parameters but with more parameters in the latent SGM and less in the encoder and decoder networks, compared to the reported LSGM-100M. However, this model performed significantly worse.

heading: F ADDITIONAL EXPERIMENTS

heading: F.1 TOY EXPERIMENTS

heading: F.1.1 ANALYTICAL SAMPLING
text: In order to test combinations of diffusions and numerical samplers in isolation, we consider a dataset for which we know the ground truth score function (for all t) analytically. In particular, we use the mixture of Normals introduced in App. E.1; see Fig. 9a for a visualization of the data distribution. In Fig. 10, we show samples for VPSDE (Euler–Maruyama (EM) sampler) and CLD (EM and SSCS samplers). For quantitative comparison, we also compute negative log-likelihoods for the three combinations (which can be done easily due to our access to the ground truth distribution): as can be seen in Tab. 8, for each number of steps n ∈ {20, 50, 100, 200} CLD with SCSS outperforms both VPSDE and CLD with EM. As discussed in Sec. 3.3, we can see in Tab. 8 that EM is not wellsuited for CLD. This is true, in particular, when using a small number of synthesis steps n (function evaluations). In Fig. 10, we see that CLD with EM leads to sampling distributions which are too broad. These results are exactly in line with the “diverging” dynamics that is observed when solving Hamiltonian dynamics with a non-symplectic integrator, such as the standard Euler method (Neal, 2011). This problematic behavior of Euler-based techniques is more pronounced when using fewer steps with larger stepsizes, which is also what we observe in our experiments. These results further motivate the use of our novel SSCS, which addresses these challenges, for sampling from our CLDbased SGMs.

heading: F.1.2 MAXIMUM LIKELIHOOD TRAINING
text: For maximum likelihood training, models based on overdamped Langevin dynamics such as VPSDE need to learn an unbounded score for t→ 0. Our model, on the other hand, only ever needs to learn a bounded score even for t = 0. For our image data experiments, we use a reweighted objective function to improve visual quality of samples (as is general practice).
Here, we also study training towards maximum likelihood on toy dataset tasks. To explore this, we repeat the neural network complexity experiment from App. E.1 with maximum likelihood training (instead of the reweighted objective). Furthermore, we also train VPSDE-based and CLD-based SGMs on a challenging toy dataset and find that CLD significantly outperforms VPSDE. We leave the study of CLD with maximum likelihood training for high-dimensional (image) datasets to future work.
Complexity Experiment. The setup of this experiment is equivalent to the setup in App. E.1 up to the training objective: in this experiment we do maximum likelihood learning, i.e., we train CLD models with the objective from Eq. (8) with λ(t) = Γβ. 8 Furthermore, we test CLD in this setup for three different values of γ. The results of this experiment can be found in Fig. 11. For CLD, we find that larger values of γ generally lead to less complex networks, in particular for smaller times t. However, even for γ = 0.04 the learned neural network is still significantly smoother than the network learned for VPSDE when a mixed score parameterization is used. 9
8For the ML objective of VPSDE, we refer the reader to Song et al. (2021b). 9VPSDE with mixed score parameterization did not converge to the target distribution, and therefore is not
included in Fig. 11.
Challenging Toy Dataset. Using the same simple ResNet architecture (less than 100k parameters) from the above experiment, we trained a VPSDE-based as well as a CLD-based SGM to maximize the likelihood of a more challenging toy dataset (the dataset is essentially “multi-scale”, as it involves both large scale—the placement of the swiss rolls—and fine scale—the swiss rolls themselves— structure). Similar to the other toy datasets, the models are trained for 1M iterations using fresh data synthesized from the data distribution in each batch at a batch size of 512.
In Fig. 12, we compare samples of the models to the data distribution. Even with our simple model architecture, CLD is able to capture the multi-scale structure of the dataset: the five rolls are adequately resembled and only a few samples are in between modes. VPSDE, on the other hand, only captures the modes, but not the fine structure. Furthermore, VPSDE has the undesired behavior of “connecting” the modes.
Overall, we conclude that also in the maximum likelihood training setting CLD is a promising diffusion showing superior behavior compared to the VPSDE in our toy experiments.

heading: F.2 CIFAR-10 - EXTENDED RESULTS
text: In this section we provide additional results on CIFAR-10. For an extended version of Tab. 3, please see Tab. 9. Note that quadratic striding improves VPSDE and CLD for all settings but VPSDE with n = {1000, 2000}. For VESDE (using PC sampling), quadratic striding significantly worsens FID numbers. This is potentially due to variance of VESDE already following an exponential schedule (see Figure 5 in Song et al. (2021c)).
Furthermore, we present additional generated samples from our model: see Fig. 13 and Fig. 14 for samples from EM-QS with 2000 evaluations and SSCS with 150 evaluations, respectively.

heading: F.3 CELEBA-HQ-256 - EXTENDED RESULTS
text: In this section we provide additional qualitative results on CelebA-HQ-256. For high quality samples using our new SSCS solver see Fig. 15.
Samples generated with Runge–Kutta 4(5) at different solver tolerances can be found in Fig 16. We found that our model still generates very good samples even for solver error tolerance of 1e−3 using an average of 129 function evaluations.

Review summaries:
The paper proposes a novel class of continuous-time diffusion-based generative models. Specifically, the paper proposes forward diffusions that are simple forms of stochastic Hamiltonian dynamics, unlike the previous methods, where linear stochastic differential equations (linear SDE) diffuse the data distribution. Thus, the proposed forward diffusion transport a joint distribution of data $x_0$ and auxiliary random variables $v_0$ to a prior joint distribution. Then, the generative models are their reverse-time diffusion processes on the joint space whose initial condition is the prior joint distribution.  In the Hamiltonian dynamics-type forward diffusion, the auxiliary random variable can be interpreted as velocity (as in physics), while data corresponds to position. The forward diffusion updates data by the velocities without any noise, while the velocities are updated via linear SDEs, similarly to deterministic Hamiltonian dynamics. In particular, the paper shows that the data distribution will be transported to a prior distribution and that velocities will map to a prior distribution similar to previous approaches. The prior joint distributions are commonly chosen to simple distributions, such as standard Normal distributions.  The authors emphasize two interesting properties of the proposed method. First, the paper points out that the reverse-time SDEs should only contain $\nabla_{v_t} \log p_t(v_t| x_t)$, but not logarithmic gradient wrt $x_t$, which implies that the model complexity (the size of models) is almost similar to the previous approaches. Second, the authors shows that $p_t(v_t| x_t)$ is Gaussians for all $t \in [0, T]$, including $p_0(v_0|x_0) = p_0(v_0)$. This indicates that $\nabla_{v_t} \log p_t(v_t| x_t)$ is potentially bounded unlike previous diffusion-based generative models, where the scores closer to data are possibly unbounded. As a result, esp. with common practice to learn scores at all time $t$'s by a single neural networks, training the proposed model will be more stable.  In order to learn the $\nabla_{v_t} \log p_t(v_t| x_t)$, the paper first propose to use denoising score matching (DSM) following previous diffusion-based generative models. However, observing that training the proposed model via DSM is unstable, the paper proposes a modified objective called hybrid score matching (HSM). HSM can be obtained by marginalizing out $v_0$ from score matching loss similarly to deriving DSM; the marginalization is possible as the auxiliary random variable defined to follow a known distribution, such as Normal distribution.  In addition, as the reverse-time diffusion is also a Hamiltonian SDE, the paper proposes a new integrator for generations, benefitting from the symplectic structure of the Hamiltonian dynamics. Note that it has been well appreciated that the discretization methods that utilize symplectic structure are more accurate than the Euler-Maruyama method under similar computational costs. Thus, the proposed integrator will have a better quality of samples, esp. when the number of discretization is small.  In order to show the effectiveness of the proposed method, the paper runs three main experiments: First, with toy experiments, the authors demonstrate that learning $\nabla_{v_t} \log p_t(v_t| x_t)$ is less difficult in comparison to learning the scores wrt the data directly. Then, the authors evaluate the generation qualities of the proposed models on two image modeling benchmark datasets; CIFAR-10 and CelebA-256. Finally, the paper demonstrates the proposed integrator's effectiveness compared to the Euler-Maruyama method, including additional analysis ablating the effect of hyperparameters of the proposed integrator.

This paper proposes a critically-damped version of score-based diffusion models, by extending the inference process to an augmented state-space and diffusing the data coupled with an auxiliary velocity variable. The authors further proposed a hybrid score matching loss for training the reverse generative process, which provides an empirical advantage of learning a conditional score function that does not blow up (e.g. to infinity) and therefore stabilizes training. For sampling, the authors derived a new numerical integration method based on first principles of statistical mechanics and MD, which they found has a better convergence behavior in practice if the computing budget is limited.  

The authors propose critically-damped Langevin dynamics (CLD) for score-based generative modeling. This consists of a higher-order Langevin dynamics scheme with particle velocity and position coupled to each other, as in Hamiltonian dynamics. The Langevin dynamics is critical in the sense that it is neither over- nor under-damped. A corresponding score matching objective is derived as an objective, with proof given that it is simply necessary to approximate the score of the velocity given the position. Empirical evidence is provided that this score is easier to estimate on a synthetic example. As DSM is analytically intractable for the higher-order scheme, Hybrid Score Matching (HSM) is proposed and the integration integral to this objective is addressed with a new numerical integration scheme. This approximation scheme, called Symmetric Splitting CLD Sample (SSCS), decomposes the SDE to be integrated into a tractable expression and a (hopefully small) Euler-Maruyama integration for improved accuracy (although still first order) overall. Synthetic examples are used in both the main text and the supplementary material to motivate the theory. Benchmark image datasets exhibit exceptionally strong performance, with improved sample efficiency after training and robust hyperparameters.

In this paper, the authors introduce a novel approach for training a score-based generative model.  With prior score based networks, generation is performed by solving a stochastic differential equation (SDE) based on Langevin Dynamics, using an estimate of gradient of the log likelihood of the underlying signal distribution.  In the present paper, the authors present a novel forward process, where diffusion is run in a joint data-velocity space.  The noise term is only applied to the velocity component.  This reduces the learning problem to only needing to learn the score of the conditional distribution of velocity given data, which is easier than learning the score of the data distribution directly.  The paper shows that the novel scheme (CLD) yields higher quality for image generation when compared to prior models of similar capacity and number of neural network evaluations.