Score-Based Generative Modeling with Critically-Damped Langevin Diffusion
Authors: Tim Dockhorn
Arash Vahdat
Karsten Kreis
Venue: ICLR 2022 Spotlight (2022)
Keywords: Score-based generative modeling
denoising diffusion models
image synthesis

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. Project page and code: https://nv-tlabs.github.io/CLD-SGM.

Published as a conference paper at ICLR 2022
SCORE-BASED GENERATIVE MODELING CRITICALLY-DAMPA LANGEVIN DIFFUSION
WITH
Tim Dockhorn 123,
Arash Vahdat
Karsten Kreis!
NVVIDI
University Waterloo
Vector Institute (avehdat, kkreis (akeei>)owwe wwwideeeidddc.oooo.o.m
tim.
122772211200766884 ML] 25 Mar 2022
ABSTRACT
Score-based generative models (SGMs) have demonstrated remarkable synthe- sis quality. SGMs rely on a diffusion process that gradually perturbs the data towards a tractable distribution, while the generative model learns 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 em- play overly simplistic diffusions, leading to unnecessarily complex denoising pro- cesses, which limit generative modeling performance. Based on connections to statistical mechanics, we propose novel (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 auxil- iary variables can be considered "velocities" that are coupled to the data variables 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 CL.D significantly outperforms solvers such as Euler Fuler-Maruyama. Our framework provides new insights into score-based denois- ing diffusion models and can be readily used for high-resolution image synthesis. Project page and code: https: htt//vvtllllla tabs.. github. io/CLD DDDD- CLD-SGM.
Langevin diffusion
INTRODUCTION
Score-based generative models (SGMs) and denoising diffusion probabilistic models have emerged 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 al., 2021), and 3D shape generation (Luo & Hu, 2021; Zhou et al., 2021). SGMs use diffusion process to gradually add noise to the data, transforming a complex data distribution into 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 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.
Work done during internship NVIDIA.

---

Published as a conference paper at ICLR 2022
Denoising Modeling V., log p(vl(x))
v(vn)
Denoising
Velocity
Component V Marginal plv.)
Diffusion
p(xa)
Smooth Denoising
Data Component X2
NW
Marginal p(xx)
Diffission Joint Distribution p(x., v.)
Xt
Diflision
In critically damped Langevin diffusion, the data augmented with velocity Vr. A diffusion coupling x, and v, run in the joint data-velocity space (probabilities in red). Noise is injected only into Ve- This leads to smooth diffusion trajectories (green) for the data X2. Denoising only requires Vv, log p(ve(x)).
Figure 1:
Inspired by statistical mechanics (Tuckerman, 2010), we propose novel forward diffusion process, the critically-damped Langevin diffusion (CLD). In CLD, the data variable, Xt (time along the diffusion), augmented with an additional "velocity" variable v, and a diffusion process is run 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. As in Hamiltonian Monte Carlo (Duane al., 1987; Neal, 201 1), 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 the conditional distribution velocity given data Vv, log PPe(v (V₂ [x)), which is arguably easier than learning the score of diffused data directly. Using techniques from molecular dynamics (Bussi & Parrinello, 2007; Tuckerman, 2010; Leimkuhler & Matthews, 2013), also derive new SDE integrator tailored to CLD's reverse-time synthesis SDE.
We extensively validate CLD and the novel SDE solver: (i) We show that the neural networks learnt in CLD--ased SGMs are smoother than those of previous SGMs. (ii) On the CIFAR-10 image modeling benchmark, we demonstrate that CLDD-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 score of the velocity-data conditional distribution Vv. Po (v (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, novel diffusion process for SGMs. (ii) We derive score matching objective for CLD, which requires only the conditional distribution of velocity given data. (iii) We propose 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.
2 BACKGROUND
Consider diffusion process E IRd defined by the Itô SDE
du, f(u,, t) dt G(u,, dw,, t [0, T].
(1)
with continuous time variable [0, T]. standard Wiener process Wt. drift coefficient f: Rdd [0, 7] + Rd and diffusion coefficient G: I&d x [0, 7] Rdxd Defining u, == ur-1, correspond- ing -eveese-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 E T):
du. ff t) G(ü,, T- t)G T -)) Vn, log (ogpr-((u)) dt G(ū, T t)dw., (2)
where Va, log (ū:) is the score function of the marginal distribution over ū, 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 (u))=Pdsta(xx Currently used SDEs (Song et al., 2021c; Kim et al., 2021) have drift and diffusion coefficients of the simple form f(x,,)==f(t f and G are chosen such that the SDE's marginal, equilibrium density is approximately Normal at time T, i.e., p<u Iz). We can then initialize xo based sample drawn from complex
and G(x, t)) Generally,
2

---

Published as a conference paper at ICLR 2022
data distribution, corresponding to far-from- from-equilibrium state. While the state X relaxes towards equilibrium via the forward diffusion, we can learn a model sg(x,,)) t) for the score Vx, log Pe(r),
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, 201 1) objective for this task is:
min E 2100 1/ (e px ((t)]lso(x,, t) Vx, log (x//xx//)))1) xx)
(3)
If f and G are affine, the conditional distribution (x) is Normal and available analyti- cally (Särkkä & Solin, 2019). Different ((t) result in different trade-offs between synthesis quality and likelihood in the generative model defined by so (x,, t) (Song et al., 2021b; Vahdat et al., 2021).
3 CRITICALY-DAMPER LANGEVIN DIFFUSION
We propose to augment the data X₁ E Rd with auxiliary velocity! variables vi E IRd and utilize a diffusion process that is run in the joint Sp-ve- With u, (Xe, ve) E IRdd we set
f(u.,) () = ((-B O
3M-1 -13M-1
G(u,, = 6
2
ut
21
8 la.
(4)
where @ denotes the Kronecker product. The coupled SDE that describes the diffusion process is
dx:; dv:
M M- -X
(-r
0a -
0 (22r3)
3dt
Bdt +
dw:,
(5)
Hamiltonian component=:;
Ormstein- Ornseein-Uhienbee process=:O
which corresponds to Langevin dynamics in each dimension. That is, each x; is independently cou- pled to a velocity vi, which explains the blockwise structure of f and G. The mass M E R is a hyperparameter that determines the coupling between the x, and v. variables; stant time rescaling chosen such that the diffusion converges to its equilibrium distribution within equilibriuu state and analogous to B(t) in previous diffusions (we could also use time-dependent (()), but found constant 3's to work well, and therefore opted for simplicity); N E R + is a friction coefficient that deter- mines 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 ex- plore complex probability distributions (Neal, 2011). The Hamiltonian component 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. The O term corresponds to an Drnstein-Uhlenbeed process (Särkkä & Solin, 2019) in the velocity component, which injects noise such that the diffu- sion dynamics properly converge to equilibrium for any 1'>0. It can be shown that the equilibrium distribution of this diffusion is (u) (x; O2, Ia) Nv;; O2,, Ma) (see App. B.2).
RR+ is a con-
t [0, (in practice, we T=1) when initialized from a data-defined
There is a crucial balance between M and I' (McCall, 2010): For I2 (underdamped Langevin dynamics) the Hamiltonian component dominates, which implies oscillatory dynamics of x, and v, that slow down convergence to equilibrium. For 1²>4 (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. For 42M (critical damping). an ideal balance is achieved and convergence to PBo(u) occurs as fast as possible in a smooth manner without oscillations (also discussion in App. A. (McCall, 2010). Hence, we propose set 1 and call the resulting diffusion critically-damppe Langevin diffusion (CLD) (see Fig. 1).
Diffusions such as the VPSDE (Song et al., 2021c) correspond to overdamped Langevin dynamics with high friction coefficients I (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 X, and v..
We call the auxiliary variables velocities, as they play 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.
3

---

Published as a conference paper at ICLR 2022
3.1 SCORE MATCHING OBJECTIVE
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(u))=p(Xg) P(vo)=Paataxx)) (Vo; Od: (MTa) with hyperparameter ><1 and let the dis- tribution diffuse towards the tractable equilibrium- or prior distribution (Eo(u). We can then learn the corresponding score functions and define CLD-base SGMs. Following a similar derivation as Song et al. (2021b), we obtain the score matching (SM) objective (see App. B.3);
min (0, (u. (att)] t) Vv. log
Notice that this objective requires only the velocity gradient of the log-density of the joint distribu- tion, i.c. Vv, lo Pe (u,). This is a direct consequence of injecting noise into the velocity variables only. Without loss of generality, Pr (r(u))=pp =p: (X,, v.)=p((vv xx (xx)). Hence,
Vv, log pe = Vv. [v. + log Pt (x:) = Vv. log pe (Ve(2x))
This means that in CLD the neural network-defined score model so Uy t) only needs to learn the score of the condi- tional distribution pe >x(ve(x)) (xx), an arguably easier task than learning the score of (x:). as in previous works, or of the joint Pr (u:). This is the case, because our velocity dis- tribution is initialized from simple Normal distribution, such that Pe (Ve) (xx) is closer to a Normal distribution for
1200
1000
30*
Ve-De
800
10
400
all t20 (and for any x,) than (x.) itself. This is most ev-
ident at (=0: The data and velocity distributions are inde-
200-
0,50
pendent at t=0 and the score of Pv xa)= (Vo) simply corresponds to the score of the Normal distribution po (Vo) from which the velocities are initialized, whereas the score of the data distribution (Xo) is highly complex and can even be unbounded (Kim et al., 2021). We empirically verify the reduced complexity of the score of (v)( in Fig. 2. We find that the score that needs to be leartt 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 CLDD-based SGMs have significantly simpler and smoother neural networks than VPSDE-based SGMs for most I, in particular when lever- aging a mixed score formulation (see next section).
0
0.0
0.3
0.4
0.6
1.0
diffussoo
cimee
300
cLD SGM wioe Misel Skare
vissuk wost Mitend Hoors
250
virept SC3t w/o Sthoo Scare
200
10a
150
10
100-
10-4
0.0
0.2
0.4
0.6
0.8
1.0
ditfusson time
Figure 2: Top: Difference {(t) (via L2 norm) between score of diffused data and score of Normal distribution. Bottom: Frobenius norm of Jacobian Jr(t) of the neural network defining the score function for different The underlying data distri- bution is a mixture of Normals. Insets: Dif- ferent axes (see App. E.1 for detailed defi- nitions of E() and Jr)))).
3.2 SCALABLE TRAINING
Practical Objective. We cannot train directly with Eq. (6), since we do not have access to the marginal dis- tribution (r.))). As presented in See. 2, we could employ denoising score matching (DSM) and instead sample ug. 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 Normal distribution over the initial velocity. Therefore, we propose a hybrid version of score matching (Hyvärinen. 2005) and denoising score matching (Vincent, 2011), which we call hybrid score matching (HSM). In HSM, we draw samples from Po (x)))=PPaaaa (Xo) as in DSM, but then diffuse those samples while marginalizing over the full initial velocity distribution po (vo)=N((
(v; Og, 7M Ia) as in regular SM (HSM is discussed in detail in App. C). Since po (v)) is Normal (and ff and G affine), p(u, Xx) is also Normal and this remains tractable. We can write this HSM objective as:
(8) In HSM, the expectation over po Pvo)) is essentially solved analytically, while DSM would use sample-basee estimate. Hence, HSM reduces the variance of training objective gradients compared to pure DSM, which we validate in App. C.I. Furthermore, when drawing sample un to diffuse in
min ErE (o, .7 Ex (Xo) Eu (p(p/xx))))) [A(t)]sa(u, t)- Vv.
4

---

Published as a conference paper at ICLR 2022
DSM, we essentially placing an infinitely sharp Normal with unbounded score (Kim et al., 2021) at ug. which requires undesirable modifications or truncation tricks for stable training (Song et al., 2021c; Vahdat et al., 2021). Hence, with DSM 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 pa (u, (x)) are provided in App. B.1.
Score Model Parametrization. (i) Ho et al. (2020) found that it can be beneficial to parameter- ize the score model predict the noise that was used in the reparametrized sampling to generate perturbed samples ur- For CLD, u, (x(xx) + L, where E, L,L, is the Cholesky de- composition of pr(u, (u(o covariance matrix, Ezd (N(zzd; O2d: I24), and (x(x)) Pe u, (o)x mean. Furthermore, Vv, log Pt (u, (xo) - -4Ed: where Ed:2d denotes those d components of E2d that actually affect Vv. log pt (u,,/x)) Xo) (since we take velocity gradients only, not all are relevant).
FEF F
DF
With Σg =
Ia.
we have Ee == V DT (2
"per-dimension" covariance matrix
(ii) Vahdat et al. (2021) showed that can be beneficial to assume that the diffused marginal distri- bution is Normal at all times and parametrize the model with a Normal score and residual "correc- tion". For CLD, the score is indeed Normal at f (due to the independently initialized x and v at (-)). Similarly, the target score is close to Normal for large t, we approach the equilibrium.
Based on (i) and (ii), we parameterize so (u, t) -000 (u, t) with ag (u,, t) = = ve + dg (u. l). where En"h corresponds to the a-u component of the "per-dimension" covariance matrix of the Normal distribution pe Pe(u:|xo = 0₂)) In other words, we assumed po (Xo) = 8(x) when defining the analytic term of the score model. Formally, -V/It" is the score of Normal distribution with covariance Êtt" Ia. Following Vahdat al. (2021), refer this parameterization mixed score parameterization. Alternative model parameterizations are possible, but we leave their exploration future work. With this definition, the HSM training objective becomes (details in App. B.3):
min 17 (Bxx-xo(xx))) Besas-WeessDaa..aaa Bee
[A(t)& |1a:22) Ea:2d (x( t)
which corresponds to training the model to predict the noise only injected into the velocity during reparametrized sampling of u, similar to noise prediction in Ho et al. (2020); Song et al. (2021c).
Objective Weightings. For A(t) B, 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), objective better suited for high quality image synthesis can be obtained by ET which corresponds to "dropping the variance prefactor" E?.
setting A(1)
3.3 SAMPLING FROM CLD-BASED SGMS
To sample from the CLD-based SGM we can either directly simulate the reverse-time diffusion pro- cess (Eq. (2)) or, alternatielly. 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 re- lied on Euler-Maruyama (EM) (Kloeden & Platen, 1992) and related methods (Ho et al., 2020; Song et al., 2021c; Jolicoeur- Deliceeu--Martineu Martineau et al., a)) We derive a new solver, tailored to CL.D-based models. Here, provide the high-level ideas and derivations (see App. D for details).
Our generative SDE can be written as (with ũ, ur X = XT V VT-1):
(dx)) (v.))
Mv, Bdt
0a PM -TM-lv V.
Os Bdt + (verddw.)
0
(u, 1) M v. Bdt
consists of a Hamiltonian component Ans an Omstein-UUlenbeed process Ao, and the score model term S. We could use EM to integrate this SDE; however, standard Fuler methods are not well-suited for Hamiltonian dynamics (Leimkuhler & Reich, 2005; Neal, 2011). Furthermore, if S was O. we could solve the SDE in closed form. This suggests the construction of novel integrator.
We use the Fokker-Planck operator² formalism (Tuckerman, 2010; Leimkuhler & Matthews, 2013; 2015). Using& similar notation Leimkuhler & Matthews (2013), the Fokker-Planck equation
Thh Fokker-Planck operator is also known as Kolmogorov operator. If the underlying dynamics is fully Hamiltonian, it corresponds the Liouville operator (Leimkuhler & Matthews, 2015; Tuckerman, 2010).
5

---

Published as a conference paper at ICLR 2022
where L3 and Lss the corresponding to the generative SDE is ōp: ( non-commuting Fokker-Plancle operators corresponding to the A and S terms, respec- tively. Expressions for L and L can be found in App. D. We can construct a formal, but intractable solution of the generative SDE i, e((2+(2)) ūg where the operator c(((2+(3) (known as the classical propagator in statistical physics) propagates states un for time according to the dynamics defined by the combined operators م4 + L. 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 Maker-Campbell Haapbell-Hausdorff formula (Trotter, 1959; Strang, 1968: Tuckerman, 2010), it can be shown that:
(((a+(()== = lim Te 6r 22a] Re [e 0l 2 N-+00
+ ((N81))
(10)
for large N E N+ and time step 8t == 1/N. The expression suggests that instead of directly evaluating the intractable el() we can discretize the dynamics over t into N pieces of step size öt, such that only need to apply the individual e e and e c many times one after another for small steps öt. A finer discretization results in a smaller error (since N=t/St, the error effectively scales as O for fixed t). Hence, this implies an integration method. Indeed, 44 u, is available in closed form, as mentioned before; however, &&t u, 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 learn the score of the conditional distribution pc(v, [xx), 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 etts 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.
4 RELATED WORK
Relations to Statistical Mechanics and Molecular Dynamics. Learning a mapping between a simple, tractable and complex distribution as in SGMs inspired by annealed importance sam- pling (Neal, 2001) and the Jarzynski equality from non-equilibrium statistical mechanics (Jarzynski, 1997a;b; 2011: Bahri et al., 2020). However, after Sohl-Dickstein al. (2015), little attention has been given to the origins of SGMs in statistical mechanics. Intuitively, in SGMs the diffusion pro- cess initialized in a -quilibbii state Ug 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 denois- ing. This "equilibration problem" 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 also arise when modeling complex data, such as natural images. Hence, the vast literature on thermostats (Ander- sen, 1980; Nosé, 1984; Hoover, 1985; Martyna al., 1992; Hü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 tech- niques 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 a 2020). We build on Song al. (2021c). which introduced the SDE framework for modern SGMs. Nachmani et al. (2021) recently introduced Gaussinn 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

---

Published as a conference paper at ICLR 2022
the diffusion process itself, whereas LSGM "simplifies the data" by first embedding it into 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 log- likelihood estimation. Furthermore, they are originally not based on the continuous time framework, which we use, and have been developed primarily for discrete-step diffusion models.
A complementary work to CLD is "Gotta Go Fast" (GGF) Pellicoeur-aattined 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 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 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 -uchh 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 accelcrates 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.
5 EXPERIMENTS
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 3=4, r-=1 (or equivalently M- in all experiments. We set the variance scaling of the inital velocity distribution to -o-=00004 and use the proposed HSM objective with the weighting A(t)=E,² 2 which promotes image quality.
Sampling. We generate model samples via: (i) Probability flow using a Runge-KKtte 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 al., 2021) (indicated by QS).
Evaluation. We measure image sample quality for CIFAR-10 via Fréchet inception distance (FID) with 50k samples (Heusel al., 2017). We also evaluate an upper bound on the negative log- likelihood (NLL): - log
log Ps (Xo, vo)-H, where H is the entropy of p(vo) and log Pa (X Vo) is an unbiased estimate of log p(Xg. vo) from the probability flow ODE (Grath- wohl et al., 2019; Song et al., 2021c). As in Vahdat al. (2021), the stochasticity of legpe((, v) prevents us from performing importance weighted NLL estimation over the velocity distribution (Burda al., 2015). We also record the number of function-nnuurl network-evaluation (NFEs) during synthesis when comparing sampling methods. All implementation details in App. B.5 and E.
5.1 IMAGE GENERATION
Following Song et al. (2021c), we focus on the widely used CIFAR- unconditional image genera- tion 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 out- performing CLD are LSGM (Vahdat et al., 2021) and //CN+/VESDEE with 2,000 step predictor- corrector (PC) sampling (Song et al., 2021c). However, LSGM uses a model with =475M pa- rameters to achieve its high performance, while we obtain our numbers with a model of x100M parameters. For fairer comparison, we trained a smaller LSGM also with 1000M parameters, which is reported as "LSGMM better FID score than "LSGM-100M"' In contrast NSCN++/VESDE, we achieve extremely strong results with much fewer NFEs (for example, see ne (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
in Tab. I (details in App. E.2.7). Our model has a significantly
7

---

Published as a conference paper at ICLR 2022
Table I: Unconditional CIFAR-10 generative performance.
Class
Model
NLL
FIDJ
CLD.SGM (Proh. Flow) (ours)
$3.31
2.25
Score
CLD-SGM (SDE) (owns)
2:23
DDPM+ VPSDE (Prob. Flow) (Song et 2021c)
3.13
3.08
DDPM VPSID (SDE) (Song et al., 2021c)
2.41
DDPM+ subb-/P (Phob. Flow) (Song et 2021e)
2,99
2.92
DDPM+ sulh-VP (SDK) (Song ah, 2021e)
2.41
NCSN++ VESDE (SDE) (Song al., 2021c)
2.20
(Vabdar al., 2021)
53.43
LSGM--00MM (Vahdatted et al., 2021)
$2.96
4,60
DDPM (Ho al., 2020)
53.75
3.17
NCSN (Song Ermon, 2019)
Score
Adversarial DSM Mariiincnaa et al. 2021b)
6.10
Likelihood SDE (Song et al., 2021b)
2.84
2.87
DDIM steps) (Sonn
20211)
FasIDDPM (100 steps) (Kong Ping. 2021)
2.86
Improved DDPM (Nichol & Dhariwal, 2021)
3.37
2,90
VDM (Kingma al., 2021)
$2.49
7.41 (4.00)
Figure 3: CIFAR- 10 samples.
(Kim al., 2021)
3.04
D3PM (Austin et al., 2021)
$3.44
7.34
Gotin Go Fastt (Iilicoeur-Mattinsuu et al., 2021a)
2.44
DDPM Distillation
Lalunan, 2021)
SNGAN
SNGAN DGflow (Nisari al., 2021)
AutoGAN (Gong ot ol., 2019)
12.4
GANs
TranGGAN (Jian 2021)
9.26
StyleGAN2 w/o
(Karra
Sty 12 w/ ADA (Karms et al., 2020)
2.92
StyleGAN2 w/ Diffang (Zhao al, 2020)
5.79
DistAug
PixelCNN (Oard et al., 2016)
3.14
65.9
Glow (K (Kingma & Dharrwal, 2018)
3.35
48.9
Aut.-Reg.
Residual Flow (Chen
2019)
3.28
46.37
Flows,
NYAE (Valdat & Kautz, 2020)
2.91
23.5
YAES,
NCP-VAE (Aneja al., 2021)
24.08
EBMs
DC.VAE Parmar al. (2021)
17.90
IGEBM (1u & Mordacch, 2019)
40.6
VAEBM (Xiao et al., 2021)
12.2
Rocovery EBM (Gao al., 2021)
3,18
9.58
Figure 4: CelebA-HQ-2556 samples.
network capacity and under NFE budgets <500, our CLD-SGM outperforms all published results in terms of FID. We attribute these positive results to our easier matching task. Furthermore, model reaches an NLL bound of 3.31, which is on par with recent works such as Nichol & Dhariwal (2021); Austin et al. (2021): Vahdat 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 CLD-SGM on -Hobb-HQ--26. but without careful hyperparameter tuning due limited compute resources. Model samples in Fig. 4 appear diverse and high-quality (additional samples in App. F).
5.2 SAMPLING SPEED AND SYNTHESIS QUALITY TRADE-OFFS
We analyze the sampling speed synthesis quality trade-off for CLD-SGMs and study SSCS's performance under different NFE budgets (Tabs. 2 and 3). We compare to Song et al. (2021c) and use EM to solve the generative SDE for their VPSDE and PC (reverse-diffusion + Langevin sampler) for the VESDE model. We also compare to the GGF (Oolicoeur-Martineauu et al., 2021a) solver for the generative SDE well as probability flow ODE sampling with a higher-order adaptive step size solver. Further, we compare to LSGM (Vahdat et al., 2021) (using our 100M- which also uses probability flow sampling. With one exception (VESDE with 2,000 NFE) our CLD-SGM outperforms all baselines, both for adaptive and fixed-step size methods. More results in App. F.2.
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 simpler ODE solver, we even outperform GGF, which is tuned towards image synthesis. (iii) Our CLD-SGM also outperforms the LSGM-I00M 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, combining the strengths of the two approaches. It would also be interesting to develop more advanced, adaptive SDE solver that leverages SSCS as the backbone
8

---

Published as a conference paper at ICLR 2022
Table 2: (right) Performance using adaptive stepsize solvers (ODE is based on probability flow, GGF simulates generative SDE). taken from JolicoerrMMartnnea et al. (2021a). LSGM corresponds the small LSGM- 100M model for fair comparison (details App. E.2.7). Error tolerances were chosen to obtain similar NFEs.
Table 3: (bottom) Performance using non-adaptive stepsize solvers (for PC, QS performed poorly). t: 2.23 FID our evaluation. Song et al. (2021c) reports 2.20 FID. See Tab. 9 in App. F.2 for extended results.
Model
Solver
NFEs
FIDI
CLD
ODE
312
2.25
VPSDE
GGF
330
2.56
VESDE
GGF
488
2.99
FID at function evaluations
Model
Sampler
no=50
150
n=-25
n=500
n= 1000
R
2000
CLD
ODE
147
2.71
VPSDE ODE
141
2.76
CLD
EM-QS
52.7
7.00
3.24
2.41
2.27
2.23
VPSDE GGF
151
2.73
CLD
SSCS-QS
20.5
3.07
2.38
2.25
2.30
2.29
VESDE
ODE
182
7.63
VPSDE
EM-QS
28.2
4.06
2.65
2.47
2.66
2.60
VESDE
GGF
170
10.15'
VESDE
PC
460
216
11.2
3.75
2.43
2.23
LSGM
ODE
131
4.60
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.
Table 4: Mass hyperpa- rameter.
5.3 ABLATION STUDIES
M-1
NLL.
FIDI
We perform ablation studies study CLD's new hyperparameters (run with a smaller version of our CLD-SGM used above; App. E for details).
$3,30
3.23
$3.37
3.14
16
<3.26
3,16
Mass Parameter: Tab. 4 shows results for CLD-SGM trained with dif- ferent M-1 (also recall that M and I are tied together via n 4.; we always in the critical-damping regime). Different mass values perform mostly similarly. Intuitively, training with smaller M- means that noise flows from the velocity variables v, into the data x; more slowly, which necessitates a larger time rescaling 3. We found that simply tying M and 33 together via B=8/M M works well and did not further fine-tune.
Table 5; Initial velocity distribution width.
Y
NLLL
FIDU
0.04
$3.37
3.14
0.4
$3.15
3.21
$3.15
3.27
Initial Velocity Distribution: Tab. shows results for CLD-SGM trained with different initial velocity variance scalings y. Varying y similarly has only small effect, but small seems slightly beneficial for FID, while the NLL bound suffers bit. Due our focus on synthesis quality as measued by FID, we opted for small y. 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 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 See. 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 crucial element for training CLD-SGMs.
We conclude that CLD does not come with difficult-to-tung hyperparameters. We expect our chosen hyperparameters to immediately translate to new tasks and models. In fact, we used the same M- yi MS and HSM settings for CIFAR-10 and SelebA-HQ-255 experiments without fine-tuning.
CONCLUSIONS
We presented critically-dmmped 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-paramettin pressivity. 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--bsed SGMs. From a technical perspective, in addition to proposing CLD, we derive CLD'ss 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.
score functions, fast synthesis, and improved ex-
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.

---

Published as a conference paper at ICLR 2022
7 ETHICS AND REPRODUCIBILITY
Our paper focuses on fundamental algorithmic advances to improve the generative modeling perfor- mance 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). I). It is worth mentioning that compared to generative adversarial networks (Goodfellow et al., 2014), very popular class of generative models, SGMs have the promise model the data more faithfully, without dropping modes and introducing problematic biases. Generally, the ethical impact of our work depends its application domain and the task hand.
To aid reproduciiility of the results and methods presented in our paper, we made source code to reproduce the main results of the paper publicly available, including detailed instructions; see project page https:: nv- t abs git wit.ub..ttti.. 1o/ CLD-SGM and the code repository https: //github com/ nv- tlabs/ / CCDDSGMM, Furthermore, all training details and hyperparameters are already in detail described in the Appendix, in particular in App. E.
REFERENCES
Guillaume Alain, Yoshua Bengio, Li Yao, Jason Yosinski, Éric Thibodeau-Laufed Saizheng Zhang, and Pascal Vincent. GSNs: generative stochastic networks. Information and Inference: Journal of the IMA, 5(2):210-22) 249, 03 2016. ISSN 2049-8766
Hans C. Andersen. Molecular dynamics simulations at constant pressure and/or temperature. The Journal of Chemical Physics, 72(4):2384-2393, 1980.
Brian DO Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Ap- plications, 12(3):313-326, 1982.
Jyoti Aneja, Alexander Schwing, Jan Kautz, and Arash Vahdat. A Contrastive Learning Ap- proach for Training Variational Autoencoder Priors. In Neural Information Processing Systems (NeurIPS), 2021.
Abdul Fatir Ansari, Ming Liang Ang. and Harold Soh. Refining Deep Generative Models via Dis- criminator Gradient Flow. In International Conference on Learning Representations, 2021.
Jacob Austin, Daniel Johnson, Jonathan Ho, Danny Tarlow, and Rianne van den Berg. Structured Denoising Diffusion Models in Discrete State-Spaces. In Neural Information Processing Systems (NeurIPS), 2021.
Yasaman Bahri, Jonathan Kadmon, Jeffrey Pennington, Sam S. Schoenholz, Jascha Sohl- Dickkktiint,, and Surya Ganguli. Statistical Mechanics of Deep Learning. Annual Review of Condensed Matter Physics, 11:0015228, 2020.
Bailey, The tools of generative art, from flash to neural networks, Art in America, 2020.
Yoshua Bengio, Eric Laufer, Guillaume Alain, and Jason Yosinski. Deep Generative Stochastic Net- works Trainable by Backprop. In Proceedings of the 31st International Conference Machine Learning, 2014.
Florian Bordes, Sina Honari, and Pascal Vincent. Learning Generate Samples from Noise through Infusion Training. 5th International Conference on Learning Representations, ICLR, 2017.
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. Importance Weighted Autoencoders. arXiv: 150.00550 00, 2015.
Giovanni Bussi and Michele Parrinello, Accurate sampling using Langevin dynamics. Phys. Rev. E, 75:056707, 2007.
Giovanni Bussi, Davide Donadio, and Michele Parrinello. Canonical sampling through velocity rescaling. The Journal of Chemical Physies, 126(1):014101, 2007.
10

---

Published as a conference paper at ICLR 2022
Michele Ceriotti, Giovanni Bussi, and Michele Parrinello. Langevin Equation with Colored Noise for Constant- Temperaturr Molecular Dynamics Simulations. Physical Review Letters, 102: 020601, 2009.
Michele Ceriotti, Michele Parrinello, Thomas E. Markland, and David E. Manolopoulos. Efficient stochastic thermostatting of path integral molecular dynamics. The Journal of Chemical Physics, 133(12):124104, 2010.
Jianfei Chen, Cheng Lu, Biqi Chenli, Jun Zhu, and Tian Tian. VFlow: More Expressive Generative Flows with Variational Data Augmentation. International Conference on Machine Learning, pp. 1660-1669. PMLR, 2020.
Nanxin Chen, Yu Zhang, Heiga Zen, Ron Weiss, Mohammad Norouzi, and William Chan. Wave- Grad: Estimating Gradients for Waveform Generation. In International Conference on Learning Representations, 2021.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud, Neural Ordinary Dif- ferential Equations. Advances in Neural Information Processing Systems, 2018.
Ricky T. Q. Chen, Jens Behrmann, David Duvenaud, and Momn-Heriik Jacobsen. Residual Flows for Invertible Generative Modeling. In Advances in Neural Information Processing Systems, 2019.
Tianqi Chen, Emily Fox, and Carlos Guestrin. Stochastic Gradient Hamiltonian Monte Carlo. In Proceedings of the 31st International Conference on Machine Learning, 2014.
Prafulla Dhariwal and Alex Nichol. Diffusion Models Beat GANs on Image Synthesis. In Neural Information Processing Systems (NeurIPS), 2021.
Nan Ding, Youhan Fang, Ryan Babbush, Changyou Chen, Robert D) Skeel, and Hartmut Neven. Bayesian Sampling Using Stochastic Gradient Thermostats. In Advances in Neural Information Processing Systems, 2014.
J. R. Dormand and P. Prince. A family of embedded Runge- Kutta formulae. Journal of Compu- tational and Applied Mathematics, 6(1):19-26, 1980.
Yilun Du and Igor Mordatch. Implicit Generation and Modeling with Energy Based Models. In Advances in Neural Information Processing Systems, pp. 3608-3618, 2019.
Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics Letters B, 195(2):216-22) 1987.
Chie Furusawa, Shinya Kitaoka, Michael Li, and Yuri Odagiri. Generative Probabilistic Image Colorization. arXiv:2109.14516: 2021.
Ruiqi Gao, Yang Song, Ben Poole, Ying Nian Wu, and Diederik P Kingma. Learning Energy- Basedd Models by Diffusion Recovery Likelihood. In International Conference on Learning Represen- tations, 2021.
Xinyu Gong, Shiyu Chang, Yifan Jiang, and Zhangyang Wang- AutoGAN: Neural Architecture Search for Generative Adversarial Networks. Proceedings of the IEEF/CVF International Conference Computer Vision, pp. 3224-3234 2234, 2019.
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley. Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative Adversarial Nets. Advances in neural informa- fion processing systems, 27, 2014.
Anirudh Goyal, Nan Rosemary Ke, Surya Ganguli, and Yoshua Bengio. Variational Walkback: Learning a Transition Operator Stochastic Recurrent Net. In Proceedings of the 31st Inter- national Conference on Neural Information Processing Systems, 2017.
Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud. FFJORD: Free-form Continuous Dynamics for Scalable Reversible Generative Models. Inter- national Conference Learning Representations, 2019.

---

Published as a conference paper at ICLR 2022
Ulrich G Haussmann and Etienne Pardoux. Time Reversal of Diffusions. The Annals of Probability, pp. -185,1205, 1986.
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. GANs Trained by a Two Time-Scale Update Rule Converge to a Local Nash Equilibrium. In I. Guyon, U. V. Luxburg, S. Bengio. H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising Diffusion Probabilistic Models. In Advances in Neural Information Processing Systems, 2020.
Jonathan Ho, Chitwan Saharia, William Chan, David Fleet, Mohammad Norouzi, and Tim Sali- mans. Cascaded Diffusion Models for High Fidelity Image Generation. arXiv:2106.15282, 2021.
William G. Hoover. Canonical dynamies: Equilibrium phase-space distributions. Physical Review A, 16695-667, 1985.
Chin-We Wei Huang, Laurent Dinh, and Aaron Courville. Augmented Normalizing Flows: Bridging the Gap Between Generative Flows and Latent Variable Models. arXiv:2002.07101, 2020.
Chin-Wei Huang, Jae Hyun Lim, and Aaron Courville. A Variational Perspective on Diffusion- Based Generative Models and Score Matching. Neural Information Processing Systems (NeurIPS), 2021.
Philippe H. Hünenberger. Thermostat Algorithms for Molecular Dynamics Simulations, volume 173 of Advanced Computer Simulation. Advances in Polymer Science. Springer, Berlin, Heidelberg, 2005.
Aapo Hyvärinen. Estimation of Non-Normalized Statistical Models by Score Matching. Journal of Machine Learning Research, 6:695-709, 2005. ISSN 1532-4435.
Christopher Jarzynski. Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach. Physical Review E, 56:5018-500 5035, 1997a.
Christopher Jarzynski. Nonequilibrium Equality for Free Energy Differences. Physical Review Letters, 78:2690-2693, 2693, 1997b.
Christopher Jarzynski. Equalities and Inequalities: Irreversibility and the Second Law of Thermo- dynamics at the Nanoscale. Anmual Review of Condensed Matter Physics, 2(1):329-351, 2011.
Myeonghun Jeong, Hyeongju Kim, Sung Jun Cheon, Byoung Jin Choi, and Nam Soo Kim. Diff- TTS: Denoising Diffusion Model for Text-to-Speech. arXiv preprint arXiv:2 10Xiv:2104.01409 2021.
Yifan Jiang, Shiyu Chang, and Zhangyang Wang. TransGAN: Two Pure Transformers Can Make One Strong GAN, and That Can Scale Up. arXiv:2102.00074, 2021.
Alexia Jolicoeur-Martnnaau, Ke Li, Rémi Piché-Taillefee Tal Kachman, and Ioannis Mitliagkas. Gotta Go Fast When Generating Data with Score-Based Models. arXiv:2/0S.140800 2021a.
Alexia BaiiceerrMMartnnau Rémi Piché-Taillefee loannis Mitliagkas, and Remi Tachet des Combes. Adversarial score matching and improved sampling for image generation. Inter- national Conference on Learning Representations, 2021b.
Heewoo Jun, Rewon Child, Mark Chen, John Schulman, Aditya Ramesh, Alec Radford, and Ilya Sutskever. Distribution Augmentation for Generative Modeling. In International Conference on Machine Learning, pp. 5006-5019, PMLR, 2020.
Tero Karras, Miika Aittala, Janne Hellsten, Samuli Laine, Jaakko Lehtinen, and Timo Aila. Training Generative Adversarial Networks with Limited Data. In Neural Information Processing Systems (NeurIPS), 2020.
Dongjun Kim, Seungjae Shin, Kyungwoo Song, Wanmo Kang, and II-Chul Moon. Score Matching Model for Unbounded Data Score. arXiv:2/06.05527, 2021.
12

---

Published as a conference paper at ICLR 2022
Diederik P Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In International Conference on Learning Representations, 2015.
Diederik P Kingma, Tim Salimans, Ben Poole, and Jonathan Ho. Variational Diffusion Models. arXiv:2107.00630, 2021.
Durk P Kingma and Prafulla Dhariwal. Glow: Generative Flow with Invertible 1x1 Convolutions. In Advances in neural information processing systems, pp. 10215-10224, 2018.
Peter E. Kloeden and Eckhard Platen. Numerical Solution of Stochastic Differential Equations.
Springer, Berlin, 1992.
Zhifeng Kong and Wei Ping. On Fast arXiv:2106.00 2021.
Sampling of Diffusion Probabilistic Models.
Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro, Dave A Versatile Diffusion Model for Audio Synthesis. In International Conference on Learning Representations, 2021.
Karsten Kreis, Kurt Kremer, Raffaello Potestio, and Mark E. Tuckerman. From classical to quantum and back: Hamiltonian adaptive resolution path integral, ring polymer, and centroid molecular dynamics. The Journal of Chemical Physics, 147(24):244104, 2017.
Benedict Leimkuhler and Charles Matthews. Rational Construction of Stochastic Numerical Meth- ods for Molecular Sampling. Applied Mathematics Research eXpress, 2013(1):34 3::4-56, 2013.
Benedict Leimkuhler and Charles Matthews. Molecular Dynamics: With Deterministic and Stochas- lic Numerical Methods. Interdisciplinary Applied Mathematics. Springer, 2015.
Benedict Leimkuhler and Sebastian Reich. Simulating Hamiltonian Dynamics. Cambridge Mono- graphs on Applied and Computational Mathematics. Cambridge University Press, 2005.
Haoying Li, Yifan Yang, Meng Chang, Huajun Feng. Zhihai Xu, Qi Li, and Yueting Chen. SRDiff:
Single Image Super- Resoluteso with Diffusion Probabilistic Models. 204v:2/044149/// 2021.
Eric Luhman and Troy Luhman. Knowledge Distillation in Iterative Generative Models for Im- proved Sampling Speed. hit::2001.02386 2021.
Shitong Luo and Wei Hu. Diffusion Probabilistic Models for 3D Point Cloud Generation. In Pro- ceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2021.
Siwei Lyu. Interpretation and Generalization of Score Matching. In Proceedings of the Twenty-Fifte Conference on Uncertainty in Artificial Intelligence, UAI '09, pp. 359-366, Arlington, Virginia, USA, 2009. AUAI Press.
Yi-An Ma, Niladri Chatterji, Xiang Cheng, Nicolas Flammarion, Peter Bartlett, and Michael I Jor- dan. There an Analog of Nesterov Acceleration for MCMC? arXiv:/902.0099: 2019.
Glenn Martyna, Michael L. Klein, and Mark Tuckerman. Hosé-Hoover chains: The canonical ensemble via continuous dynamics. The Journal of Chemical Physics, 97(4):2635-243, 1992.
Martin W McCall. Classical Mechanics: From Newton to Einstein: A Modern Introduction, 2nd Edition. Wiley, Hoboken, N.J., 2010.
Chenlin Meng, Yang Song, Jiaming Song, Jiajun Wu, Jun- Yan Zhu, and Stefano Ermon. SDEdit: Image Synthesis and Editing with Stochastic Differential Equations. arXiv:2108.0107., 2021.
Yisroel Mirsky and Wenke Lee. The Creation and Detection of Deepfakes: A Survey. ACM Comput. Surv., 54(1), 2021.
Gautam Mittal, Jesse Engel, Curtis Hawthorne, and Ian Simon. Symbolic music generation with diffusion models. In Proceedings of the 22nd International Society for Music Information Re- trieval Conference, 2021. URL https: : archives. smir. ir. .et iittiismiiii /22///ppexx/ 000058. .df.
13

---

Published as a conference paper at ICLR 2022
Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral Normalization for Generative Adversarial Networks. In International Conference on Learning Representations (ICLR), 2018.
Wenlong Mou, Yi-An Ma, Martin J. Wainwright, Peter L. Bartlett, and Michael Jordan. High- Order Langevin Diffusion Yields an Accelerated MCMC Algorithm. Journal of Machine Learn- ing Research, 22(42):1-41, 2021.
Javier R. Movellan. Contrastive Divergence in Gaussian Diffusions. Neural Computation, 20(9): 2238-2252, 2008.
Eliya Nachmani, Robin San Roman, and Lior Wolf. Non Gaussian Denoising Diffusion Models. arXiv:2106.07582, 2021.
Radford M. Neal. Annealed importance sampling. Statistics and Computing, 2001.
Radford M. Neal. MCMC Using Hamiltonian Dynamics. Handbook of Markov Chain Monte Carlo, 14:113-16,, 2011.
Thanh Thi Nguyen, Quoc Viet Hung Nguyen, Cuong M. Nguyen, Dung Nguyen, Duc Thanh Nguyen, and Saeid Nahavandi. Deep Learning for Deepfakes Creation and Detection: Sur- vey. arXiv: 11573, 2021.
Alexander Quinn Nichol and Prafulla Dhariwal. Improved Denoising Diffusion Probabilistic Mod- cls. In International Conference on Machine Learning, 2021.
Shuichi Nosé. A unified formulation of the constant temperature molecular dynamics methods. The Journal of Chemical Physics, 8 11(1)5 1984.
Aaron van den Oord, Nal Kalchbrenner, and Koray Kavukcuoglu. Pixel Recurrent Neural Networks. International Conference on Machine Learning, 2016.
Gaurav Parmar, Dacheng Li, Kwonjoon Lee, and Zhuowen Tu. Dual Contradistinetive Generative Autoencoder. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pp. 823-832, 2021.
B. T. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computa- tional Mathematics and Mathematical Physics, 4(5):1-17, 1964. ISSN 0041-5553.
A. Roberts. Modify the Improved Fuler scheme to integrate stochastic differential equations. arXiv: 2012.
Chitwan Saharia, Jonathan Ho, William Chan, Tim Salimans, David J Fleet, and Mohammad Norouzi. Image Super-Resolution via Iterative Refinement. arXiv:2104.07633 2021.
Robin San-Roman, Eliya Nachmani, and Lior Wolf. Noise Estimation for Generative Diffusion Models. arXiv:2104. 026000 2021.
Simo Särkkä and Aro Solin. Applied Stochastic Differential Equations, volume 10. Cambridge University Press, 2019.
Hiroshi Sasaki, Chris G. Willcocks, and Toby P. Breckon. UNIT-DDPM: UNpaired Image Transla- tion with Denoising Diffusion Probabilistic Models. arXiv:2104.05338, 2021.
Xiaocheng Shang. Zhanxing Zhu, Benedict Leimkuhler, and Amos J Storkey. Covariance- Controlled Adaptive Langevin Thermostat for Large-Scale Bayesian Sampling. In Advances in Neural Information Processing Systems, 2015.
Abhishek Sinha, Jiaming Song. Chenlin Meng, and Stefano Ermon. D2C: Ilffusion-Denoising Models for Few-shot Conditional Generation. NXXi::/06.068/// 2021.
Jascha Sohl-Dickstein Peter Battaglino, and Michael R. DeWeese. Minimum Probability Flow Learning. In Proceedings of the 28th International Conference on International Conference on Machine Learning, 2011
14

---

Published as a conference paper at ICLR 2022
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep Unsuper- vised Learning using Nonequilibrium Thermodynamics. In International Conference on Machine Learning. 2015.
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising Diffusion Implicit Models. In Inter- national Conference on Learning Representations, 2021a.
Yang Song and Stefano Ermon. Generative Modeling by Estimating Gradients of the Data Distribu- tion. In Proceedings of the 33rd Annual Conference on Neural Information Processing Systems, 2019.
Yang Song, Conor Durkan, lain Murray, and Stefano Ermon. Maximum Likelihood Training of Score-Based Diffusion Models. In Neural Information Processing Systems (NeurIPS), 2021b.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-Based Generative Modeling through Stochastic Differential Equations. In Interna- tional Conference on Learning Representations, 2021c.
Gilbert Strang. On the Construction and Comparison of Difference Schemes. SIAM Journal on Numerical Analysis, 5(3):506-317, 1968.
H. F. Trotter. On the Product of Semi- -roips of Operators. Proceedings of the American Mathe- matical Society, 10:545-551, 1959.
M. Tuckerman, The Journal of Chemical Physics, 1990)119,- 1992.
Berne, and G. J. Martyna. Reversible multiple time scale molecular dynamics.
Mark E. Tuckerman. Statistical Mechanics: Theory and Molecular Simulation. Oxford University Press, New York, 2010.
Cristian Vaccari and Andrew Chadwick. Deepfakes and Disinformation: Exploring the Impact of Synthetic Political Video Deception, Uncertainty, and Trust in News. Social Media Society, 6(1):2056305 120566000000000000 20903408, 2020.
Arash Vahdat and Jan Kautz. NVAE: A Deep Hicrarchical Variational Autoencoder. In Neural Information Processing Systems (NeurIPS). 2020.
Arash Vahdat, Karsten Kreis, and Jan Kautz. Score-based Generative Modeling in Latent Space. In Neural Information Processing Systems (NeurIPS), 2021.
Pascal Vincent. A Connection Between Score Matching and Denoising Autoencoders. Neural Computation, 23(7):1661-1674, 2011.
Daniel Watson, Jonathan Ho. Mohammad Norouzi, and William Chan. Learning to Efficiently
Sample from Diffusion Probabilistic Models. arXiv:2106.0388 2021.
Zhisheng Xiao, Karsten Kreis, Jan Kautz, and Arash Vahdat. VAEBM: A Symbiosis between Vari- ational Autoencoders and Energy-based Models. In International Conference on Learning Rep-
resentations, 2021.
Richard Zhang. Making Convolutional Networks Invarannv Again. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 7324-7334. PMLR, 09-15 Jun 2019.
Shengyu Zhao, Zhijian Liu, Lin, Jun- Van Zhu, and Song Han. Differentiable Augmentation for Data-Efficient GAN Training. Advances in Neural Information Processing Systems, 33, 2020.
Lingi Zhou, Yilun Du, and Jiajun Wu. 3D Shape Generation and Completion through Point- Voxel Diffusion. In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), 2021.
15

---

Published as a conference paper at ICLR 2022
CONTENTS
Introduction
Background
3
Crimieally-Dmmped Langevin Diffusion
3.1
Score Matching Objective
3.2 Scalable Training
3.3 Sampling from CLD-based SGMs
4
Related Work
5
Experiments
5.1
Image Generation
5.2
Sampling Speed and Synthesis Quality Trade-Offss
8
5.3
Ablation Studies
6
Conclusions
7 Ethics and Reproducibility
10
References
10
A Langevin Dynamies
18
A.1 Different Damping Ratios
18
A.2 Very High Friction Limit and Connections previous SDEs in SGMs
18
B Criticaly-Daapped Langevin Diffusion
20
B.1 Perturbation Kernel
20
B.2
Convergence and Equilibrium
21
B.3
CLD Objective
................................
22
B.4
CLD-specific Implementation Details
24
B.5
Lower Bounds and Probability Flow ODE
24
B.6 On Introducing a Hamiltonian Component into the Diffusion
................................
26
HSM: Hybrid Score Matching
26
C.1
Gradient Variance Reduction via HSM
28
Symmetric Splitting CLD Sampler (SSCS)
29
D.I Background
29
Derivation and Analysis
Implementation and Experiment Details
34
16

---

Published as a conference paper at ICLR 2022
E.1
Score and Jacobian Experiments
34
E.2
Image Modeling Experiments
35
E.2.1
Training Details and Model Architectures
36
E.2.2
CIFAR-I0 Results VESDE and VPSDE
................................
36
E.2.3
Quadratic Striding
37
E.2.4
Denoising
................................
37
E.2.5
Solver Error Tolerances for Runge Kutta 4(5)
Rutte-Kutta
38
E.2.6
Ablation Experiments
38
E.2.7
LSGM-100M Model
38
Additional Experiments
39
F.1
Toy Experiments
39
F.1.1
Analytical Sampling
F.1.2
Maximum Likelihood Training
39
F2
CIFAR-10 Extended Results
F3
CelebA-HQ-256 - Extended Results
45
Proofs of Perturbation Kernels
50
G.1 Forward Diffusion...
50
G.1.1 Proof of Correctness of the Mean
50
G.1.2 Proof of Correctness of the Covariance
51
G.2 Analytical Splitting Term of SSCS
52
G.2.1 Proof of Correctness of the Mean
53
G.2.2 Proof of Correctness of the Covariance
53
17

---

Published as a conference paper at ICLR 2022
A LANGEVIN DYNAMICS
Here, discuss different aspects of Langevin dynamics. Recall the Langevin dynamics, Eq. (5), from the main paper:
dx: dv:
1 -X,
M
Oa + -IM 'v, Bet (V2T3)
0
3dt
dw.
(11)
Hamiltonian component==:
Ornstein- OrastcinUhllebbeck Uhlenbeck process=:O
A. DIFFERENT DAMPING RATIOS
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 x, and v, that slow down convergence equilibrium.
(ii) For r > 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- cintiial-damming)..... ideal balance is achieved and convergence to Peo ( occcurs quickly in a smooth manner without oscillations.
In Fig. 5, we visualize diffusion trajectories according Langevin dynamics in the different damping regimes. We observe that underdamped Langevin dynamics show undesired oscillatory behavior, while overdamped Langevin dynamics perform very inefficientl,, too. Critical-damping achieves a good balance between the two and mixes and converges quickly. In fact, it can be shown be optimal terms of convergence; see, for example, McCall (2010).
Consequently, we propose to set 1 r
critical := 4M in CLD.
A.2 VERY HIGH FRICTION LIMIT AND CONNECTIONS TO PREVIOUS SDES IN SGMS
Let us re-write the above Langevin dynamics and consider the more general case with time- dependent B()):
dx; = M v,)dd,,
(12)
dv: (t)
TM V
B())d
(13)
potential
(ii)): Friction term
(iv):
To solve this SDE, let assume simple Fuler-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):
X x, + 3(1M 1n
(14)
Vn+1 =
((t)x,, 5t
B() 'v S (3(()) StI,).
(15)
(i) current stop velocity
(8)) patential term
(iii): friction term
(iv): noise term
Now, let us assume a friction coefficient -= max a(t) Since the time step öt 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:
M
Xn+1 = =X, 3(1)M S6
(16)
Vn+1 = - + N(Oa, St1a).
(17)
18

---

Published as a conference paper at ICLR 2022
a
r = 0.1 I critical
b
r
critical
p(xo)
Denoising
((2))
p (x0)
Denoising
p(x1)
W
Diflusion
Diffusion
p(vo)
Denoising
v(v1)
(vo)
Denoising
pl
Diffusion
Diffusion
r 10T critical
(d) r
r max
((xo)
Denoising
x1
p p(xo)
Denoising
p {1
nw
Diffusion
Diffusion
p(vo)
Denoising
p(v1)
p(vo)
Denoising
p (((1)
Diffusion
Diffusion
Figure 5: Langevin dynamics in different damping regimes. Each pair of visualizations corresponds to the (coupled) evolution of data and velocities We show the marginal (red) probabilities and the projections of the (green) trajectories. The probabilities always correspond the optimal setting I' = Pemical (recall that Pemical = 2/M and T max- = M/(B(t)ōtt; see Sec. A.2). The trajectories correspond to different Langevin trajectories run in the different regimes with indicated friction enefficients Т. We see in (b). that for critical damping the trajectories quickly explore the space and converge according the distribution indicated by the underlying probability. In the damped regime (a), even though trajectories mix quickly we observe undesired oscillatory behavior. For over-damped Langevin dynamics, (c) and (d), the x, trajectories mix and converge only very slowly. Note that the visualized diffusion uses different hyperparameters compared the diffusion shown in Fig. 1 in the main text: Here, have chosen much larger 9, such that also the slow overdamped Langevin dynamics trajectories shown here mix little bit over the visualized diffusion time (while the probability distribution and the trajectories critical damping converge almost instantly).
Now velocity update, Eq. (17), does not depend on the current step velocity on the right-hand- side anymore. Hence, we can insert Eq. (17) directly into Eq. (16) and obtain:
X,+1 X, (()) xx, xt, + V28())5tM--X(00 8tI.) Xn. -(112 -01))2--xx,0 xx,0t + V20()))0MM- Ia).
(18)
Re-defining 5t' 512 and B²(1) = 3() we obtain
X X2 -3' -8 + (2 Ia),
(19)
which corresponds to the high-friction overdamped Langevin dynamics that are frequently run, for example, train energy- based generative models (Du & Mordatch, 2019; Xiao et al., 2021). Let's further absorb the mass M-1 and the time step St' into the time rescaling, defining B(1)
19

---

Published as a conference paper at ICLR 2022
2,3 ( We obtain:
1
3()xx + V B(eEN(O... 1a)
Xn+1 = Xn
2
= (1
1
3()))xx, + VBUNN Ia)
(20)
2
R3 1 - -(1)x,, + B())( Ia),
where the last approximation is true for sufficiently small B(I). However, this expression corre- sponds to
Xa+1 M(x++1: / B(t)x,, B(t)
(21)
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. (11).
we assume diffusion as above but with the potential term (ii) set to O, we can similarly derive the VESDE Song et al. (2021c) as 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, inject less noise converge the equilibrium, i.e., the prior.
B CRITIALLY-DAMPED LANGEVIN DIFFUSION
Here, present further details about proposed aitically-dampeds Langevin diffusion (CLD). We provide the derivations and formulas that were not presented the main paper in the interest of brevity.
B.1 PERTURBATION KERNEL.
To recap from the main text. in this work we propose to augment the data x, E Rd with auxiliary ve- locity variables V, E Rdd We then run the following diffusion process in the joint Xy-Vee
[dx; (dv,
du,
((,,, C(uu,()dww
(22)
0
8(t)M-1
f(u, (J(t) @ Iaju,, ʃ(1) Btt)
-((0)L-11)
(23)
- @ la, G(I) ( (2
0
(24)
where w, a standard Wiener process 2dd and 8: [0, T] - RS is time rescaling. 3 In particular, consider the pitticall-damppe Langevin diffusion which can be obtained by setting M r2/4, resulting in the following drift kernel
(u, t) (f) & Ia)u, fCLD())= (-9()
43(1)r-2) -43(1)1-1)
(25)
Since we only consider the critieally-damped case in this work, redefine f fCLD and f == FCLD for simplicity. Since our drift f and diffusion G' coefficients are affine, u, is Nor- mally distributed for all E [0, 7] if us is Normally distributed at 0 (Särkkä & Solin, 2019). In particular, given that ug ~ N (uo: No, En
Eo 0 Ia), where En = diag(2 2g") is positive
For our experiments, only used constant 8; however, for generality, present all derivations for time-dependent (t))
20

---

Published as a conference paper at ICLR 2022
semi-definite diagonal 2-by-2 matrix (we restrict our derivation to diagonal covariance matrices t == O for simplicity, since in our situation velocity and data are generally independent at = 0), we
derive expressions for μe and Σ,, the mean and the covariance matrix of u, respectively.
Following Särkkä & Solin (2019) (Section 6.1), the mean and covariance matrix of u, obey the following respective ordinary differential equations (ODEs)
dup dt (J(t) @ La)te
(26)
dEc
(f(t) @ LJ), + ((F(1) @ ()) + (G(()(()))
@ la-
(27)
dt
Notating 140- (xo: vo) the solutions to the above ODEs are
F4 [2 + 4B() vv + Xo
e 2B(1)
(28)
-B(()xx 2 Vo
and
Σ; = Σg ② Id.
(29)
Σ,
2x D
4B(1J)
(30)
Fize 1
F
F e 1B((T - + 44B( (-0 1) + 43 (I)T~ (=6" 2) + 16B(() -
v
(31)
27 B(1)2% + 4B((t)1- -B (1)1 ( - 2) (B2
(32)
E. _ (613(10)) -1) + B(t)I+++++ + 28" (1 +1+ 2
+ B(t) (2% 2)
(33)
where ((1) - 1o (i) dř. For constant ((1) = (as is used in all our experiments), we simply have B(t) t3. The correctness of the proposed mean and covariance matrix can be verified by simply plugging them back into their respective ODEs; see App. G.I.
With the above derivations, we can find analytical expressions for the perturbation kernel p(u,) I-). For example, when conditioning on initial data and velocity samples Xo and vo (as in denoising score matching (DSM)), the mean and covariance matrix of the perturbation kernel p(u,/up) can be obtained by setting μ4 (Xo, vo) Eg = O, and 28 = 0.
In our experiments, the initial velocity distribution is set to (O,,,, Ia). Conditioning only on initial data samples Xo and marginalizing over the full initial velocity distribution (as in our hybrid score matching (HSM), see Sec. C), the mean and covariance matrix of the perturbation kernel p(u (xo) can be obtained by setting p40 (xo, 0a) 20 0, and M y.
B.2 CONVERGENCE AND EQUILIBRIUM
Our CLD-based training of SGMs-as well as denoising diffusion models more generally-relies on the fact that the diffusion converges towards analytically tractable equilibrium distribution for sufficiently large t. In fact, from the above equations we can easily see that,
lim D =1 1, t
(34)
lim (+000
= O, DEV
(35)
lim yow 1 t-y00
M, lim pt t-000
(36)
O2d,
(37)
which establishes (u) = N/ 0,, Nv; 0g, ML.).
Notice that our CLD is an instantiation of the more general Langevin dynamics defined by
(xx dv.
(xx, log (p))(xx)) ( Bdtt + (-
M-lv:
0g (-TM- Bdt +
0
=
dw,-
(38)
21

---

Published as a conference paper at ICLR 2022
which has the equilibrium distribution (eo(u) ( (x) (v: O2, MTa) (Leimkuhler & Matthews, 2015; Tuckerman, 2010). However, the perturbation kernel of this Langevin dynamics is not avail- able analytically anymore for arbitrary Ppar (x). In our case, however, tractable Prou (x) = (x; Oa, Ia). Note that this corresponds to the classical "harmonic oscillator" problem from physics.
have the analytically
B.3 CLD OBJECTIVE
To derive the objective for training CLD-based SGMs, we start with derivation that targets maxi- mum likelihood training in similar fashion to Song et al. (2021b). Let po and be two densities, then
(RI(PO |l go) = DKT. (po go) - (PT (r) + (RL (r)
- I 0DkL (pe gc) dt + (KL(PT || (T),
(39)
Öt
where Pt and lt are the marginal densities of po and go, respectively, diffused by critically- damped Langevin diffusion. As has been shown in Song et al. (2021b), Eq. (39) can be written as a mixture (over t) of score matching losses. To this end, let us consider the Fokker-Planck equation associated with the critically-damped Langevin diffusion:
(p:(ur) Vu, 1 (G())(()) @ Ia) Vu.P((u)) ((() @ La)u.] Øt
Vu. [up (u))l t) (G(t)G(t)) (G())))))) Ia) Vu, logpe(u.) (f(t) Ia) u
(40)
Similarly, we have (u.) Vu, (ug () (u:)). Assuming log pr and log (u,) are at smooth functions with at most polynomial growth at infinity, we have
lim (u,. (u,) = lim h t) (u,) = 0.
(41)
Using the above fact, we can compute the time-derivative of the Kullback Leibler divergence be- tween pt and qe
0 (Pr |l 4t)
/ pi (u() log Pr( du, ge (u:)
Öt
Ot.
pe (wi) (up( t) (u, Vu, logp:(u:) - Vu. log (u:)] du,
pe(ue) (Vu. log) V., log (G(I)((IT 6 Ia) (Vu, log).(u.)
Vu. logg:((;)) du,
-(()T / (e((u))) logp((u.) Vv, log qe (u:))2 du,-
(42)
Notice that due the form of G(!), we now have only gradients with respect the velocity com- ponent Vy. Combining the above with Eq. (39), we have
(KL(Po |l go) E..=44(..T).uumpppii
TB(() V lo:):(u:) V., log (u, )12
+ DKL(pr || gr)
EE..uq0..TT [B([) lo:):(u:) V., log qr uu ||
(43)
Note that the approximation holds if pr is sufficiently "close" dr- We obtain general objective function by replacing 1'3(t) with an arbitrary function A(t), i.c.,
E2 [(t)] Vo. Vv, log (u)))
(44)
22

---

Published as a conference paper at ICLR 2022
As shown in App. C, the above be rewritten, up to irrelevant constant terms, either of the following two objectives:
HSM(A(E))
xx) [(t))lVv, log xo) Vv, log q: (u
(45)
DSM(X(I) Be-ero.,T], (o0po)))) (oppuu//u
[(())VV log uo) Vv, log qe (u
(46)
For both HSM and DSM, we have shown in App. B.1 that the perturbation kernels (u, xo) and
Pu( un) are Normal distributions with the following structure of the covariance matrix:
Eer
Σ. E, La
Σ.
(47)
Dan
We can use this fact to compute the gradient Vn, log Pr (e 1))
Vn, n Pa
-V 1 (u: me) Dr - 14)
=-E, 22, 104-10 1e)
(48)
--,," 1 pk)
-, -,- E2d.
where E2d (N 22)) and E, = L,L, is the Cholesky factorization of the covariance matrix E, Note that the structure Σ, implies that L; L, @ Id. where L,L, the Cholesky factorization of Σ,, i.e
Lxx D
L L
L
(49)
SYE
Furthermore, we have
LJT=
Li
= Li T 2
la
Mer
X Ia (
0
(50)
E
2 Ia-
Using the above, we can compute
V log
.) V., (u; laaad
= -- TT d:2d
(51)
= -2E
where
(52)
EF (2
and 6d:2d denotes those (latter) d components of e2d that actually affect v., log logp:(u.l.).
Note that & depends on the conditioning in the perturbation kernel, and therefore l is different for DSM, which is based on p(u, un). and HSM, which based p(u, xo). Therefore, will henceforth refer to PHSM and εDSM if distinction of the two cases is necessary (otherwise we will simply refer to for both).
As discussed in Section 3.2, we model Vv, log qe (u:) as so so t) = -lap (u, t). Plugging every- thing back into our objective functions, Eq. (45) and Eq. (46), we obtain
HSM(A(!)) (1 uu lxo)
X(E) (CHIMM)² 2 || co(u,, 1)
(53)
((t) (EDSM)² 2 ||Ea:ed ao(u,
DSM(A(t))
(54)
23

---

Published as a conference paper at ICLR 2022
140
HSM w/
num le 12
HSM w/ Caum 1c 9
120
HSM W/
6num le 6
DSM w/ 6num le - 12
100
DSM W/ Cuum = 1c 9
DSM W/ Cum
6
80
60
40
20
10
10
10
10
10
10
Figure 6: Comparison of P (in green) and EDSM (in orange) for our main hyperparameter setting with M = 0.25 and y = 0.04. In contrast to 2DSM, fHSM analytically bounded. Nevertheless, numerical computation
be unstable (even when using double precision) in which case adding a numerical stabilization of Enum = 10-9 to the covariance matrix before computing & suffices to make HSM work (see App. B.4).
where u, is sampled via reparameterization:
1 cg:d L. Eo.d + L Ea:2d)
u μt L,E Mt
(55)
Note again that L, 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":
HSM (A(() (RNTMM) -2)
||Ea::2 uo 4)
1
(56)
DSM (>(t) (eDSM) -2)
||ea:2d ao(u, (57)
B.4 CLD-SPECIFIC IMPLEMENTATION DETAILS
Analytically, #HSM is bounded (in particular, e 1/V^M), whereas &DSM is diverging for t 0. In practice, however, we found that computation of #HSM also be numerically unstable, even when using double precision. As is common practice for computing Cholesky decompositions, add numerical stabilization matrix Cnum Iaa to E, before computing l, In Fig. we visualize t PHSM and @DSM for different values of mum using our main experimental setup of M 0.25 and
7 0.04 (also, recall that practice we have
1). Note that very small numerical stabilization of 6mm 10 in combination with the use of double precision makes HSM work in practice.
B.5 LOWER BOUNDS AND PROBABILITY FLOW ODE
Given the score model so (u, t), we can synthesize novel samples via simulating the reverse-ting diffusion SDE Eq. (2) in the main text. This be achieved, for example, via novel SSCS,
24

---

Published as a conference paper at ICLR 2022
Euler-Mann Maruyama, methods such as GGF (2021b;c) have shown that a corresponding ordinary differential equation can be defined that gener- ates samples from the same distribution, in case so (u, t) models the ground truth scores perfectly. This ODE is:
et al., 2021a). However, Song al.
dů, = ( (I),, Vn, log logpr-r(u)) di (58)
This ODE is often referred as the probability flow ODE. We can use generate novel data by sampling the prior and solving this ODE, like previous works (Song et al., 2021c). Note that practice (o(u,,)) won't be 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 perfor- mance is usually very similar almost the same, when we have learnt good sæ so(u, t). In addition to sampling the generative SDE in our paper, we also sample from our CL.D-based SGMs via this probability flow approach.
With the definition of our CLD, the ODE becomes:
dxi dv.,
M 've
0g 1) M-lv.))
3dt +
[s
3dt
(59)
X
A/r
Notice the interesting form of this probability flow ODE for CLD: It corresponds to Hamiltonian dynamics (A'r) plus the score function term S'. Compared to the generative SDE (Sec. 3.3), the Constein-Uhlenbect term disappears. Generally, symplectic integrators are best suited for integrat- ing Hamiltonian systems (Neal, 2011; Tuckerman, 2010; Leimkuhler & Reich, 2005). However, ODE is not perfectly Hamiltonian, due the score term, and modern non-symplectic methods, such as the higher-order udaptive-step size Runge- Kutta 4(5) ODE integrator (Dormand & Prince, 1980), which we use practice 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 liglliillhood 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, CLD the input into this ODE is not just the data Xo. but also the velocity variable In this case, we can still calculate lower bound on the log-likelihood:
logp(xo) = log (frizo, volavo)
p(vto) p(xo, vo) p(vo)
= log
dvo
(60)
Evor~p(vo) p(xpp( vo) log p(vo) -vorpp(vo [logp(xo, vo)l (p(vo))
where H (p(vo) denotes the entropy of p(vo) (we have (I((vv)) 1/2 log (2TeeMM) We can obtain a stochastic, but unbiased estimate of p(xo vo) R log Pe (xo. vo) via solving the prob- ability flow ODE with initial conditions (Xu, Vu) and calculating stochastic estimate of the log-
determinant of the Jacobian via Hutchinson's 's trace estimator (and also calculating the probability of the output under the prior), as done Normalizing flows (Chen et al.. 2018; Grathwohl al., 2019) and previous works on SGMs (Song et al., 2021c;b). In the main paper, we report the negative of Eq. (60) as our upper bound on the negative ligliikld (NILL).
Note that this bound can be potentially quite loose. In principle, it would be desirable to per- form an importance-wiighter estimation of the log-likelihood, as in importance-weightee autoen- coders (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 Pe (Xp, vo). 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.
25

---

Published as a conference paper at ICLR 2022
B.6 ON INTRODUCING HAMILTONIAN COMPONENT INTO THE DIFFUSION
Here, we provide additional high-level intuitions and motivations about adding the Hamiltonian component to the diffusion process, 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 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 (flippngg 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 Iearn a complex score function to invert the Hamiltonian dynamics. Our thought experiment shows that we should strive for 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 setwork-based modeling, but more on a simple and analytical Hamiltonian component. Such an SDE could then be solved easier with 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 essentially outperform all previous results in the literature (on CIFAR-10).
We would also like to offer a second perspective, inspired by the Markov chain Monte Carlo (MCMC) literature. In MCMC, "mixing" helps 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-equilibriug data distribution) and we need traverse the space as efficiently as possible 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-wwe 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. In fact, leveraging Hamiltonian dynamics to accelerate sampling is popular in the MCMC field (Neal, 2011). Note that this line of reasoning extends to thermostating techniques from statistical mechan- ics and molecular dynamics, which essentially tackle similar problems like MCMC methods from the statistics literature (see discussion in Sec. 4).
state (the
C HSM: HYBRID SCORE MATCHING
We begin by recalling our objective function from App. B.3 (Eq. (44)):
Em-ujo.27 [(()) Vr, log P: (u.) se (u, f 2)]3]]
(61)
where se (u, f) is our score model. In the following, we dissect the "score matching" part of the above objective:
LsM == Eu [IV., log so(u,, 1)
(62)
= Eurp (w,) |150 (u,, 1) - 25(0) + C2(().
where C2(() = Bu) [V. log (r:))(u))) [2]] and S()) is the cross term discussed below. Fol- lowing Vincent (2011), we can rewrite Lsm as an equivalent (up addition of time-dependent
26

---

Published as a conference paper at ICLR 2022
constant) denoising score matching objective LDSM:
DDMM == _uo (m) vv. po un) - so(u,, t)ll
(63)
= LsM + Ct(t) C2(t),
where (t) = E _un (on V log uo)l Something that might not neces-
sarily be quite obvious is that there is no fundamental need to "denoise" with the distribution p(uo) (this is, use samples from the joint Xo-Vo distribution p(un), perturb them, and learn the score for denoising).
Instead, we can "denoise" only with the data distribution p(xo) and marginalize over the entire initial velocity distribution p( vo), which results in
LHSM = Exo www.p (u.|xo) V log xo) - (u, 1)l
Bu (u.) (ou 1) 2F xixx (v,, log peruee (u, xo); so 4)}] + Ca(l), (64)
x 11 and (+,-) donates the inner product (notation chosen to consistent with Vincent (2011)). In our case, this makes sense since p( p(vo) is Normal, and therefore (as shown in App B.1), the perturbation kernel pe(ue xx) is still Normal.
where Ca(t) = Exx (ux (ur V Po
In the following, for completeness, we redo the derivation of Vincent (2011) and show that LsM is equivalent to LHSM (up to addition of constant). Starting from S(0), we have
S(0) (u. vVv, log pogg:(uu) (ux), so 4))
pr( vVv. logp: (u:). so( ()) du,
[n (Vv.Pa(u)) so du:
Vv. Pe(ue (x)))u(XX) dxw, so(u..t) du,
(65)
Pe(u: xx)Pa V., log (u xx) dx so(u,,)) t) du.
X0
pr(u: (xo) (V., logp: (ue xo), so (ur. ()) dxo du,
uu J80
E Pxxo=po(xx (V log pi (n, x so (u, ().
Hence, we have that
LUSM = LsM + C4(1) C2(1).
(66)
This further implies that
LHSM LDSM + -1(1)- C3(1).
-
(67)
Using the analysis from App B.I, we realize that C3 and C4 can be simplified to d (2PSM) and (MHSM)2 respectively. Here, we used the fact that the expected squared norm of a multivariate standard Normal random variable is equal to its dimension, i.e., Ee |lel2 d. This analysis then simplifies Eq. (67) to
CHSM = LDSM + d ( (
(68)
Using this relation, we can also find connection between our CLD objective functions from App. B.3. In particular, we have HSM()(()) Ex-ecjo.21 [A(I)
= E,- CDsMl + + d' [A(t) ((eDSM) -(2,SM))
(69)
DSM(A(!) + E'Eare [>(1) ( (cHIS))))
27

---

Published as a conference paper at ICLR 2022
C.I GRADIENT VARIANCE REDUCTION VIA HSM
Above, we derived that LHSM - DSM + const, so one might wonder why we advocate for HSM over DSM. As discussed in Sec. 3.2, one advantage of HSM is that avoids unbounded scores at - O. However, there is second advantage: 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. (69), we have
V HHHS(((()) = E [A(1)Vo CHSM)
= EE-utOT
(70)
= VoDSM(((!)),
where A are the learnable parameters of the neural network. Instead of comparing the above expec-
tations directly, we instead compare Vt)V&LHSM with Vet)V for t E [0, 1] (we use T
in all experiments) at discretized time values (as is done in practice). Replacing LUSM and LDSM with a single Monte Carlo estimate (as is used in practice), we have
X()VV R X(L)V V so(ue.)) Vv, log Xo) so(
(71)
p(x)). Uy (u, xo),
VODDS
V s)) vv, logee(u, un) se t)l2
(72)
Uo p(u), u1 ~ (u, uo):
where we applied the chain-rule. Note that in Eq. (71) and Eq. (72), u, is sampled from the same distribution. Hence. X(t) (use((, t) t) acts as common scaling factor, with the variance difference between HSM and DSM originating from the squared norm term. Hence, we ignore A(t) Vasg (u, t) and only focus our analysis on the gradient of the norm terms, which we can further simplify:
2 (n log xb) so 1)1 = so (u, t) -v, log pr Xo) =: KHSM, (73)
and
V ()( Vv, le (u)) so so 1)12 = so(u, t) Vv, log Uo) =: KDSM. (74)
evaluate KHSM and KDSM for all data points in the CIFAR-10 training set. We choose so to be our trained ablation CLD model (with the standard setup of M- B = see Sec. E.2.1 for model details). We then these samples to compute the empirical covariance matrices COVHSM and COVDSM of the random variables KuSM and KDSM. respectively.
We explore this difference in realistic setup: in particular,
As is common practice in statistics, we consider only the trace of the estimated covariance matrices." 4 The trace of the covariance matrix (of random variable) is also commonly referred to 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 y 1. For DSM, both y 0.04 and = result in very large Tr(Cov) values for small For large t, Tr( Tr(Cos (Cov) practice, we found that DSM is even unstable for small 7. Given this analysis, we believe this is due to the large gradient variance for small t. conclusion, these results demonstrate a clear variance reduction by the HSM objective, in particular for large 7. Ultimately, this is expected: In HSM, we are effectively integrating out the initial velocity distribution when estimating gradients, while in DSM we use noisy samples for the initial velocity.
considerably smaller for 7 0.04 than for 7 = 1. However, in
Note that re-introducing X(t) weightings would allow us to scale the Tr(Cov) curves according to the "reweighted" objective the maximum likelihood objective. However, we believe it is most instructive to directly analyze the gradient of the relevant norm term itself.
Arguablyy the most prominent algorithm that follows this practice is principal component analysis (PCA).
28

---

Published as a conference paper at ICLR 2022
7
107
107
7
HSM
10
10
10a
Tr(Cov)
10
10-4
10
10a
0.0
0.2
0.4
0.8
1.0
diffusion time
Figure 7: Traces of the estimated covariance matrices.
D SYMMETRIC SPLITTING CLD SAMPLER (SSCS)
In this section, present a more complete derivation and analysis of our novel Symmetric Splitting CLD Sampler (SSCS).
D.1 BACKGROUND
Our derivation is inspired by methods from the statistical mechanics and molecular dynamics liter- ature. In particular, are leveraging symmetric splitting techniques as well (Fokker-Planck) operator concepts. The high-level idea symmetric splitting as well as the operator formalism are well-explained in Tuckerman (2010), in particular in their Section 3.10, which includes simple amples. 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.
D.2 DERIVATION AND ANALYSIS
Generative SDE. From Sec. 3.3, recall that our generative SDE can be written as (with i, ur-1- x, = XT, Vi VT-)):
dxx
M 'vi
Oa
02 ( t) M -M
dv.
Xe
Bdt
-1'M
V. (2
Adt +
2r
Bdt.
Ag
Ao
(75)
Fokker-Planck Equation and Fokker-Planck Operators. The evolution of the probability distri- bution ( is described by the general Fokker-Planck equation (Särkkä & Solin, 2019):
(ur)
2d
2d
2d
82
11t ( T - (p))=((u)))
+
Duy T -)pr--((x))l
01
Dū, Oū,
1=1
(76)
29

---

Published as a conference paper at ICLR 2022
with
v
Og B
(-TM-19 ) -IMM
),)) 8+ (2r
2r [s(ū, T - + M vvi) 3, (77)
0d
p(ū T -t)=
X
D(ū, T -
@ Ia-
(78)
r3
For our SDE, can write the Fokker-Planck equation in short form
(p.) (Ć) + (2))+r-(uu
(79)
0t
with the Mikkerr-Plancle operators (defined via their action on functions of the variables o(ū))):
(ū)((.) v, o(ū)) 3x, Vo, ((u)) --M- [v.ó(ū,)] r,A., (ūx),
(80)
(30(ū,) == -V [(s(ū,, T 1) M p(ū)))
(81)
d
0²
8²
4v. = Σ 1=11
(82)
dz?
We 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 de- scribe the time evolution of states u, under the stochastic dynamies defined by the SDE. Given an initial state ūg we can construct formal solution to the generative SDE via (Tuckerman, 2010; Leimkuhler Matthews, 2015):
ū, = c(( No:
(83)
where the operator ef is known as the classical propagator that propagates states ün for time according to the dynamics defined by the combined avoid confusion, note that in Eq. (83) the operator (²()) is applied üg in element-wise or "vectorized" fashion on all elements of úg in parallel). The problem with that expression is that cannot analytically evaluate it. However, we can leverage
operators E + C (to
design an integration method.
Symmetric Splitting Integration. Using the symmetric Trotter theorem or Strang splitting formula as well the Baker- CaakeerCambelll--e Hausdorff formula (Trotter, 1959; Strang, 1968: Tuckerman, 2010), it can be shown that:
(C)+ lim [e ot a 221
Rs 2 e 22 2 + O(N81))
(84)
N
for large N E N+ and time step St = 1/N. The expression suggests that instead of directly evaluating the intractable ((27+ we can discretize the dynamics over / into N pieces of step size ôt, such that we only need to apply the individual e E1 and 5tt2 many times one after another for small time steps St. A finer discretization implies a smaller error (since N-t/o öt the error effectively scales ( 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 splicting scheme error, but not due to the evaluation of the individual updates. Such techniques are, 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 4 A'A ú,: This term describes the stochastic evolution of ü, under the dy- namics of an SDE like Eq. (75), but with S set to zero. However, if S is set to zero, the remaining
"I principle. the error of the splitting scheme is defined more specifically by the commutator of the non- commuting Fokker-Planck operators. See, for example Leimkuhler & Matthews (2013; 2015); Tuckerman (2010).
30

---

Published as a conference paper at ICLR 2022
SDE has affine drift and diffusion coefficients. In that case, if the input is Normal (or discrete state corresponding to a Normal with O variance) then the distribution is Normal at all times and we solve the differential equations for the
can calculate the evolution analytically. In particular, mean j de and covariance E a: of the Normal (see Sec. B) and obtain
233 r -43 f r 2v, + X
(85)
3 x1
x
2 -
Vv V/
as well as
Σ 2 = Eos & La,
(86)
2
FEB
Eg
(87)
D 2
51 r-
FE
48 r -
43
8
3
r-2
(88)
2
2
81
8
r-1
(89)
2
2
Eue
T
13 r 1
-1)
+ 3
5t
(
8t
(90)
The correctness of the proposed mean and covariance matrix can be verified by simply plugging them back in their respective ODEs; see App. G.2.
Now, we can write the action of the the propagator e 4a on a state u, as:
e Lin ur N( z a (ūr). Eg).
(91)
(ii): Next, we need analyze ct ú, Unfortunately, we cannot calculate the action of the prop- agator ettl3 on in analytically and we need to make an approximation. From Eq. (75), we can easily see that the propagator étt describes the evolution of the velocity component V; for time step dt under the ODE (this can be easily seen by noticing that the S term in Eq. (75) only acts on the velocity component of the joint state úz):
dv. 2.1 S,, T - t) + M- dt.
(92)
We propose to simply solve this ODE for the step öt via a simple step of Euler's method, resulting in:
stC's ū. R ūr + öt
2 (u T t) + M 0(882)
O2
(93)
5ihs Bake
U, + O
with the informal definition
8t Eo Fater
ū, + 5t 231 (ū, T t) + M -
02
(94)
Error Analysis. It now instructive to study the overall error of our proposed integrator. With the additional Euler integration in one of the splitting terms, we have
((5)+(3)
Rv
e
2
( Bate
+ 066t))
e E. N + +NNSt)
$522 E2
$ 8-
N
+ NO(812)
(95)
Te 4 E
24 24
N O(t)).
where we used N dr and only kept the dominating error terms of lowest order in St. We see that, just like Fuler's method, also our SSCS is first-order integrator with local error ~812 and global
31

---

Published as a conference paper at ICLR 2022
(a) Euler-Maruyama
Buler-based Approximation
uv 5t ( (2))
Ug
11g
13g2
U1
Euler-based
(b) Symmetric Splitting CLD Sampler
Approximation
Exact
uN 6 6
uo
Analytical Step
Exact Analytical Step
14
12
i1
Figure 8: Conceptual visualization of our new SSCS sampler and comparison to Fuler- MaleerMaruyamaa (for image synthesis): (a) EM-based sampling, in each integration step the entire SDE is inte- grated using Euler-based approximation. This can be formally expressed solving the full-step propagator exp 5t(2% + [g)] via Fuler-based approximation for N small steps of size St (see red steps; for simplicity, this visualization assumes constant öt). (b): In contrast, our SSCS the propagator is partitioned into analytically tractable component exp bt 2 (blue) and the score model term exp 5t (brown). Only the latter requires numerical approximation, which results in an overall more accurate integration scheme.
enror ~ôt, which can be also seen from the last two lines of Eq. (95). This is expected, considering that we used an Euler step for the term. Nevertheless, as long as the dynamics not dominated by the S component, our proposed integration scheme is still 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.
To this end, recall that the model only needs to learn the score of the conditional distribution Pr(v)(x)) which is close to Normal for much of the diffusion. in which the S term will in- deed be small. This suggests that the generative SDE dynamics are in fact dominated by Arr and Ao in practice. From another perspective, note that (recalling that so (ue, -l,cro (u,, () with (o, -, vi/E a's (u,, from Sec. 3.2):
s s(ū, T t) M- = -,06
Vz Dey
+ M -
(96)
= --,066uu, v.
M
For large parts of the diffusion, D/ is indeed close to M, such that the v, term is very small (this cancellation is the reason why we pulled the M- term into the component). In Sec. 3, we have also seen that our neural network component a' (u, t) can be much smoother than that of previous SGMs. Overall, this suggests that the error of SSCS indeed might be smaller than the error would obtain when applying naive Euler-Maruyana. 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 SSCS and EM make only very small errors and are expected 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 of CLD-SGMs with relatively few synthesis steps.
32

---

Published as a conference paper at ICLR 2022
Note that error analysis of stochastic differential equation solvers is usually performed in terms of weak and strong convergence (Kloeden & Platen, 1992). Due to the use of Euler's s method for the S component, argued above, expect our SSCS to formally have the same weak and strong vergence properties like EM, this is, weak convergence of order 1 and strong convergence of order 1 well, since the noise is additive in our case (and assuming appropriate smoothness conditions for the drift and diffusion coefficients; furthermore, without additive noise, we would have strong convergence of order 0.5). We leave more detailed analysis to future work.
In practice, we do not use SSCS to integrate all the way from (=0 to t=T,, but only up to (-7 and perform denoising step, similar to previous works (Jolicoeur- Martinnea et al., 2021a; Song et al.. 2021c). It is worth noting that our SSCS scheme would also be applicable when we used dime-dependent ((t), as in our more general derivation of the CLD perturbation kernel in App. B. However, since we only used constant ,3 in the main paper, we also presented SSCS in that way.
promising direction for future work would be to extend SSCS to adaptive step sizes and to use techniques to facilitate higher-order integration, while still leveraging the advantages of SSCS.
SSCS Algorithm. Finally, we summarize SSCS terms of a concise algorithm:
Algorithm Symmetric Splitting CLD Sampler (SSCS)
Input: Score function se (ū, T t), CLD parameters I, B, M r /1, number of sampling steps N, step sizes (8t, 2 N= chosen such that := T - Ex=o 8t,, 2 0 (stepsizes can vary, for example in QS). Output: Synthesized model sample X'y, along with a velocity sample v'y.
Xu N(xo; Od, Ia), Va ~ NiVoo Oa, Mla). uo = (Xo, Vo) t= 0 for = to N - do
b Draw initial prior samples from pro (u) b. Initiaaize time
Bn+ ~ -N(ū fig Paa
(ūn).
Eag 7
> First half-step: apply exp( 6tu L')
02
a' - 4444 8bn
P Full step: apply exp exp( on u
2 [s T t) + M vv.]
ün+1 NN(nnt am (0',
E Emy
>> Second half-step: apply exp -= L%) on u' Update time
t t 6tn
end for
U'y + üN (-)
Os -T8M-1) Ia uN +e
3M-1
D Denoising >> Extract output data and velocity samples
2.
(XW, v'x) = u'y
Note that the algorithm uses the expressions in Eqs. (85) and (86) for ea, and Σ, Furthermore, in practice in the denoising step at the end, we usually only update the Xxy component of u'y since we are only interested in the data sample. This saves us the final neural network call during denoising, which only affects the V'y component (also see App. E.2.4). However, we wrote the algorithm in the general way, which also allows to correctly generate the velocity sample v'y- In Fig. 8, we show conceptual visualization of our SSCS and contrast it to EM.
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 (Leimkuhler & Reich, 2005; Tuckerman, 2010; Neal, 2011; Leimkuhler & Matthews, 2015), which follows a similar structure our SSCS. However, it is not important in case, as the only computationally costly operation is in the center full step, which involves the neural network evaluation. The first and last half-steps come virtually computational cost.
33

---

Published as a conference paper at ICLR 2022
CLD SGM wy Mixeed Seare
CLD SCM wyj Mixud Store
1000
10a
CLD
VI5IB
VPSOE SCM Mixed Store
VPSDD eG34
800
((
600
10-3
1072
400
10-4
101aa
0.00
0.25
diffasion time
0,50
0.75
1.00
50
10
10-2
10a
didfusion time
dilfusion time
dilfusion time
(a) Difference E(L) (via L2 norm) between score of (b) Frobenius norm of Jacobian Jr(() of the neural network defining the function for different
diffused data and score of Normal distribution.
Figure 9: Toy experiments for mixture of Normals dataset.
E IMPLEMENTATION AND EXPERIMENT DETAILS
E. SCORE AND JACOBIAN EXPERIMENTS
In this section, we provide details for the experiments presented in Sec. 3.1. For both experiments, we consider wwo-dimensional simple mixture of Normals of the form
Pdata (x) =
g
(x),
(97)
where (x) -N(x; pn; 0.04 and
a/2
11
112
a/2
P3
a/2
a/2
#11
0
Ns
0
16
0122
a/2
a
#7
1aa/22
1g
and a 2-2. 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 Po(xo) Paata(x), we have
(x():) = Σ Di (x:),
(98)
k=1
where p( 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) with analytical access to the diffused marginal (x(x,) ps (u.). This allows us to perform interesting analyses that would be impossible when working, for example, with image data. We visualize the data distribution in Fig. 10.
Score experiment: We empirically verify the reduced complexity of the score of pf (v(llx)). which learned in CLD, compared the score of pr (xx), which is learned in VPSDE. To avoid scaling issues between VPSDE and CLD, we chose M = 7 = 1 for CLD in this experiment: this results an equilibrium distribution of (0 I2) (for both data and velocity components, which are inde- pendent at equilibrium), which is the same as the equilibrium distribution of the VPSDE. We then
TThe Leapfrog integrator corresponds the velocity Verlet integrator in molecular dynamies.
34

---

Published as a conference paper at ICLR 2022
(a) Data Pdata
(b) CLD w/ MS
(c) CLD w/ MS (d) VPSDE w/ MS (e) VPSDE w/o MS
Figure 10: Mixture of Normals: data and trained models (samples).
measure the difference of the respective scores at time / and the equilibrium (or prior) scores, i.e. (recall that the score of a Normal distribution p(x) = (02,,, I2) is simply Vx ((x) = -x).
<VPSDE (t) = Ex, V x. logPi(x;) + x 1
(99)
CLDD (1) = Eu Puu) Vv. log xx) + vellz.
(100)
The expectations are approximated using 10t samples from p(x)) and p(u,) for VPSDE and CLD, respectively. As be seen in Fig. 9a, (CCD(() (1) smaller than <VPSDE(() (l) for all / [0, T]. The difference is particularly striking for small time values l. Other previous SDEs, such as the VESDE, APSDE, etc., are expected 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 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 small ResNet architecture (less than 100k parameters) for each of the following four setups: both CLD and VPSDE each with and without mixed score parameterization. The mixed score of the VPSDE simply assumes standard Normal data distribution (which is also the equilibrium distribution of VPSDE) resulting in adding X to the score function. Formally, -x; is the score of Normal distribution with unit variance.
We train the models for 1M iterations using fresh data synthesized from Pdara at a batch size of 512. The model and data distributions are visualized in Fig. 10. We see that all models have learnt good representations of the data. We measure the complexity of the trained neural networks using the squared Frobenius norm of the networks' Jacobians. For CLD, we have
J CCLD (D) = Eurpp(u)) Vn, do(u))))
(101)
Similarly, for the VPSDE we compute
J (YPSDE(t) = Expx))x) V (x,))1
(102)
For both CLD and VPSDE, expectations are again approximated using 100h samples. As can be in Fig. 9b the neural network complexity is significantly lower for CLD compared to VPSDE. A mixed score formulation further 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 smoother as measured by Jt(t). In large- scale experiments, this would mean that, given similar model capacity, CLD-based SGM could potentially have 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 CIFAR-10 benchmark.
E.2 IMAGE MODELING EXPERIMENTS
We perform image modeling experiments CIFAR-10 10 well as CelebA-HQ-256 We report FID scores on CIFAR-10 for our main model for various different solvers; see Tab. 3 and Tab. 2. We further present generated samples for both models in Sec. 5 using Maler- with 2000
35

---

Published as a conference paper at ICLR 2022
Hyperparameter
CIFARIO (Main)
CelebA (Qualitative)
CIFARIO (Ablation)
Model
EMA rate
0.9999
0.9999
0.9999
# of ResBlooks per Resolution Normalization
8
Group Normalization
Group Normalization
Group Normalization
Scaling by a Nonlincarity
x
X
Swish 16
Swish 16
Swish
Attention resolution
16
Embedding type Progressive
Fourier None Residual Sum
Positional None None
Positional
None
Progressive input
None
Progressive combine
N/A
N/A
Finite Impulse Response (Zhang, 2019)
V
X
X
# of parameters
108M
R 68.
2 39M
Training
# of iterations
800k
320k
IM 100k
# of warmup iterations
100k
100k
Optimizer
Adam /
Adam
Adam
Mixed precision
Learning rate Gradient norm clipping
2 10
10
2 10
1.0
1.0
1.0
Dropout
0.1
0.1
0.1
Batch size per GPU # GPUs
8
16
16
t-sampling cutoff during training
10-3
10-5
10-5
SDE
.M
0.25
0,25
varies
Y
0.04
0.04
varies
8
4
4
varies
Cnum
10% 9
10
10 9
Table 6: Model architectures as well as SDE and training setups for our experiments on CIFAR-10 and CelebA-HQ-256.
quadratic striding steps and Runge-Kutta 4(5) for CIFARIO and CelebA-HQ-256,, respectively. We present additional samples for various solver settings in App. F. All (average) NFEs for the Runge- Kutta solver computed using batch size of 128.
E.2.1 TRAINING DETAILS AND MODEL ARCHITECTURES
Our models are based on the NCSN+ and the DDPM++ architectures from Song et al. (2021c). Importantly, we changed the number input channels from three six facilitate the additional velocity variables. Note that the number of additional neural network parameters due to this change negligible.
For fair a comparison, we train our models using the same t-sampling cutoff during training as is used for VESDE and VPSDE in Song al. (2021c). Note, however, that this not strictly necessary for CLD as we do not have any "blow-up" of the SDE due to unbounded scores as /. -> (also see Fig. 18 and Fig. 19).
We summarize three model architectures as well as SDE and training setups in Tab. 6.
E.2.2 CIFAR-10 RESULTS FOR VESDE AND VPSDE
The results reported for VESDE and VPSDE using the GGF sampler are taken from Jolicoeur- Martineau al. (2021a). All other results for VESDE and VPSDE generated using the provided PyTorch code as well as the provided checkpoints from Song et al. (2021c).7 We used EM and PC
httpps://ee thu.. wwwphh.w.....aacccccccc.eeecaaae
36

---

Published as a conference paper at ICLR 2022
to sample from the VPSDE and VESDE models, respectively (see Sec. 5.2), since these choices correspond to their recommended settings."
Furthermore, in App. F.2 also used DDIM (Song et al., 2021a) to sample the VPSDE. DDIM's update rule is
X
0t-1 X +
s0 t)] -0_100(5((xx
(103)
at
where a2 eep(p00 B(t) dt of of - exp(-5o a(t) dt). and ((t) = 0.1 + 19.9t.
E.2.3 QUADRATIC STRIDING
When we simulate our generative SDE numerically, using for example EM or our SSCS, we need to choose a time discretization. Given 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 {ōt, 2 0);-0 NNRE-1 with 86; Nsr Vi [0, NNFF 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 evaluations (neural network calls) 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 useful heuristic, denoted quadratic striding (QS), which discretizes the integration interval such that the evaluation times follow quadratic schedule and the individual time steps follow a linear schedule. We also used this QS approach in our experiments.
We formally define it as follows (assuming time interval j0.0, 1.0] here for simplicity): Denote the evaluation times as T, (including 0.0 and 1.0) and define:
T i
Vi [0, NNFE]-
(104)
Hence,
81, = T, T:-1
c, c2i
Vi E NNEE),
(105)
and c = N
to ensure that TNFE = 1.0.
This describes the time steps going from t t = 1. During synthesis, however, going backwards. Hence, we can define our time steps as
81j C, [2NN 2j 1]
Vj E [1, NNFF)
(106)
where j now counts time steps in the other direction. Note that this can be easily adapted to general integration intervals [e, T].
E.2.4 DENOISING
As has been pointed out in Jolicoeur-Martinaau et al. (2021b), samples that are generated with models similar ours can contain noise that is hard to detect visually but worsens FID scores significantly.
Denoising Formulas. For a fair comparison we use the same denoising setup for all experiments we conducted (including VESDE (PC/ODE) and VPSDE (EM/ODE)) except for LSGM. We sim- ulate the underlying generative ODE/SDE until the time cutoff = 10-3 and then take a single denoising step of the form
O₂ ()(
Un = Ug - (f + (ea)++GGGie G((,,,,)
(107)
This denoising step can be considered as an Huler-Maruyama step without noise injection. For SDEs acting on data directly (VESDE, VPSDE, etc.) the corresponding denoising formula is
X X - 8f (x,, + + eG G((,, so (x(x);)).
(108)
'https: / / col ab. researcc googl e. ccm Denoising has not heen used in the original LSGM work (Vahdat al., 2021) and is not needed in their case, since the output of the latent SGM lives in smooth latent space and is further processed by decoder.
dr dr ve/ ddRR_ [dRR_ogRNN km wt..tttttP2.... wtw..ognn......aaa22.a2 uu
37

---

Published as a conference paper at ICLR 2022
Table 7: Influence of denoising step on FID scores (using our main CIFAR-IO model).
FID
function evaluations
Sampler
Denoising
n==00
n==50
SSCS
/ 81.1
2.30
SSCS-QS /
20.5
2.25
SSCS X
78.9
2.32
SSCS-QS
X
28.5
2.3
Influence of Denoising on Results. For SDEs acting in the data space directly, has been reported that this denoising step is crucial to obtain good FID scores Jolicocur-Martinenl et al. (2021b); Song al. (2021c). When we simulate the generative probability flow ODE we found that denoising important in order for the Lunge-Kutta solver not to "blow-up" as 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. More specifically, the neural network only denoises the velocity component. However, diffusion coefficients of CLD in the denoising formula in Eq. (107), we obtain
primarily interested in the data component. Putting the drift and
Ug u, -8 ((9B
3M-1 1 -TM
Og 3r
e))
@ Iá
Ur
Xo -33M-
- Ve-
(109)
E.2.5 SOLVER ERROR TOLERANCES FOR RUNGE-KUTTA 4(5)
In Tab. 2, we report FID scores for Runge-Kutta 4(5) solver (Dormand & Prince, 1980) well as the "Gotta Go Fast" solver from Dolicoeur-Martnne et al. (2021a 1a) (see their Table 1). For simulating CLD with Runge-Kutta 4(5) chose the solver error tolerances hit certain regimes of NFEs to facilitate comparisons with VPSDE and VESDE. We obtain a mean number of function evaluations of 312 and 137 using Kutte-Kutta 4(5) solver error tolerances of 10 and 10-3 respectively. For VESDE, VPSDE and LSGM we used 10-5 as the ODE solver error tolerance, following the recommended default setups (Song et al., 2021c; Vahdat et al., 2021). These values are used for both relative and absolute error tolerances.
E.2.6 ABLATION EXPERIMENTS
The model architecture used for all ablation experiments can be found in Tab. 6. As pointed out in Sec. we found that the hyperparameters y and M only have small effects on CIFAR- 10 FID scores. 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 1088M parameters, while the original CIFAR-10 Latent SGM from Vahdat et al. (2021), to which we compare in Tab. 1, uses x476M parameters. To establish a fairer comparison between our CLD-based SGMs and LSGM (Vahdat al., 2021), we train an- other smaller LSGM model with 110M parameters. To do this, we followed the exact setup of the "CIFAR-10 (balanced)" model from LSGM (see Table minor modifications: We used 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 the residual cells per scale from 8 to 4. With these modifications the resulting "LSGM- 100M" uses only x109M 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, model is trained in exactly the same way as the bigger, original models in Vahdat et al. (2021).
Vahdat et al. (2021)), with a few
38

---

Published as a conference paper at ICLR 2022
Performance (measured negative log lilkliiiioodd using analytical scores for non adaa-tdtiie stepsize solvers for varying numbers of synthesis steps n (function evaluations).
Table 8:
p(r) at r function evaluations
Model
Sampler
n==0
n=
n=100
n=200
CLD CLD
EM
60.6
9.71
0.72
-1.04
SSCS
10.5
1.55
-1.25
-1.54
VPSDE
EM
14.2
4.68
-0.35
-1.11
For evaluation, we follow the recommended setting by Vahdat al. (2021) and use the same Runge- Kutta 4(5) ODE solver with an error tolerance of 10-8 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. and 2 in the main text.
Note that we also tried training model following the "CIFAR-10 10 (best FID)" setup, but found train- ing to be unstable (however, the orignal "CIFAR-1O (best FID)" model from Vahdat al. (2021) only performs marginally better in FID than their "CIFAR-10 (balanced)" model anyway). Fur- thermore, 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.
F ADDITIONAL EXPERIMENTS
F.1 TOY EXPERIMENTS
F.1.1 ANALYTICAL SAMPLING
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. particular, use the mixture of Normals introduced in App. E.I; see Fig. 10a for visualization of the data distribution. In Fig. 11. we show samples for VPSDE Maruyaaaa (EM) sampler) and CLD (EM and SSCS samplers). For quantitative comparison, we also compute negative ligeliiels for the three combinations (which can done easily due to our access to the ground truth distribution): as can be seen in Tab. 8, for each number of steps E (20, 50, 100, 200) CLD with SCSS outperforms both VPSDE and CLD with EM. discussed in Sec. 3.3, we can see in Tab. 8 that EM is not well- suited for CLD. This is true, in particular, when using small number of synthesis steps (function evaluations). In Fig. 11, 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-symples sympleeleecic 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 our novel SSCS, which addresses these challenges, for sampling from our CLD- based SGMs.
1.1 MAXIMUM LIKELIHOOD TRAINING
For maximum likelihood training, models based on overdamped Langevin dynamics such as VPSDE need to learn an unbounded score for - O. Our model, on the other hand, only ever needs to learn a bounded score even for 0. For our image data experiments, we use 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. with maximum likelihood training (instead of the reweighted objective). Furthermore, we also train VPSDE- and CLD-based SGMs 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.
39

---

Published as a conference paper at ICLR 2022
(a) VPSDE+ EM, 20
(b) CM,++M,,M 20 20
(c) CLD SSCS, n 20
(d) VPSDE + EM, n = 50
(c) CLD + EM, n = 50
(f) CLD+ + SSCS, n = 50
(g) VPSDE+ EM, n = 100
(h) CLD + EM, n 100
(i) CLD + SSCS, n = 100
() VPSDE+ EM, 200
(k) CLD EM. 200
(1) CLD SSCS, 200
Figure 11: Mixture of Normals: numerical simulation with analytical score function for different diffusions (VPSDE with EM vs. CLD with EM/SSCS) and number of synthesis steps n. A visual- ization of the data distribution can be found in Fig. 10a.
40

---

Published as a conference paper at ICLR 2022
VPSDE SGM w/o Mixed Score
CLD SGM w/ Mixed Score and 7
CLD SGM w/o Mixed Score and 7
300
CLD SGM w/ Mixed Score and 7 =
CLD SGM w/o Mixed Score and 7 0.4
CLD SGM w/ Mixed Scure and y
CLD SGM w/o Mixed Score and 7 0.04
250
10a
200
JF(t)
10
150
Jr(t)
10 10-2
10 10-4
100
50
10
10
10
diffusion time t
0
0.0
0.2
0.4
0.6
0.8
1.0
diffusion time t
Figure 12: Frobenius norm JF of the neural network defining the score function for different t.
Complexity Experiment. The setup of this experiment is equivalent the setup in App. E. 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 A(l) r 10 Furthermore, we test CLD in this setup for three different values of 7. The results of this experiment can be found in Fig. 12. For CLD, we find that larger values of > generally lead less complex networks, in particular for smaller times t. However, even for 7 = 0.04 the learned neural network is still significantly smoother than the network learned for the VPSDE when mixed score parameterization is used." !!
Challenging Toy Dataset. Using the same simple ResNet architecture (less than 100k parameters) from the above experiment, we trained a VPSDE-based the likelihood of more challenging toy dataset (the dataset is essentially "multi-scale" as it involves both large scale-the the placement of the swiss rolls-and and fine scale- tte swiss rolls themselves- structure). Similar to the other toy datasets, the models trained for 1M iterations using fresh data synthesized from the data distribution in each batch at batch size of 512.
well CLD--baed SGM maximize
In Fig. 13, 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 ade- quately resembled and only few samples are in between modes. VPSDE, the other hand, only captures the main 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 the VPSDE in our toy experiments.
For the MI. objective the VPSDE, we refer the reader to Song et al. (2021b). The VPSDE-hased model with mixed score parameterization did not converge the target distribution, and therefore not included Fig. 12.
41

---

Published as a conference paper at ICLR 2022
(a) Data
VPSDE
(c) CLD
Figure 13: Data distribution and model samples for multi-scale toy experiment.
Table 9: Performance using non-adaptive step size solvers. Extended version of Tab. 3.
FID function evaluations +
Model
Sampler
n==0
n=150
n=575
n=500
==1000
n=2000
CLD
EM
143
31.5
10.9
3.96
2.50
2.27
CLD
EM-QS
52.7
7.00
3.24
2.41
2.27
2.23
CLD
SSCS
81.1
10.5
2.86
2.30
2.32
2.29
CLD
SSCS-QS
20.5
3.07
2.38
2.25
2.30
2.29
VPSDE
EM
92.0
30.3
13.1
4.42
2.46
2.43
VPSDE
EM-QS
28.2
4.06
2.65
2.47
2.66
2.60
VPSDE DDIM
6.04
4.04
3.53
3.26
3.09
3,01
VPSDE
DDIM-QS
3.78
3.15
3.05
2.99
2.96
2.95
VESDE
460
216
11.2
3.75
2,43
2.23
VESDE
PC-QS
461
388
155
5.47
11.4
11.2
F.2 CIFAR-10 EXTENDED RESULTS
In this section, provide additional results the CIFAR-10 image modeling benchmark.
An extended version of Tab. 3 (sampling the generative SDE with different fixed-step size solvers for different compute budgets) including additional baselines can be found in Tab. 9, Note that time stepping with quadratic striding (QS) improves sampling from VPSDE- and CLD-based models for all settings except for the combination of VPSDE and EM sampling in the setting For the VESDE (using PC sampling), QS significantly worsens FID scores. The reason for this could be that the variance of the VESDE already follows exponential schedule (see Fig. 5 in Song al. (2021c)). We additionally present results for the VPSDE using the DDIM (Denoising Diffusion Implicit Models) sampler (Song et al., 2021a). As observed by Song et al. (2021a), QS also helps for DDIM. Importantly. for any n 2 150, our CLD with our novel SSCS (and QS) even outperforms 50. DDIM performs better. It needs to be mentioned, however, that the DDIM sampler was specifically designed for few-step sampling, whereas our CLD with SSCS is derived in a general fashion without this particular regime in mind. In particular, DDIM sampling non-Markovian sampling method and it is not clear how to calculate the log- likelihood of hold-out validation data under this non-Markoviat synthesis approach. Nevertheless, it would be interesting to also explore non-Markovian DDIM-inspire techniques for CLD-SGMs to further improve sampling speed in CLD-SGMs.
{1000, 2000).
DDIM. Only for
can be interpreted
Note that DDIM results shown in Tab. are better than those presented in Song et al. (2021a) itself, because we are relying on the DDPM+ model trained in Song et al. (2021c), whereas Song et al. (2021a) uses the DDPM model from Ho al. (2020).
Finally, we present additional generated samples from our CLD-SGM model: see Fig. 14 and Fig. 15 for samples from EM-QS with 2000 evaluations and SSCS-QS with 150 evaluations, respectively.
42

---

Published as a conference paper at ICLR 2022
Figure 14: Additional samples using EM-QS with 2000 function evaluations. This setup gave us our best FID score of 2.23.
43

---

Published as a conference paper at ICLR 2022
Figure 15: Additional samples using SSCS-QS. This setup resulted in an FID score 3.07 using only 150 function evaluations.
44

---

Published as a conference paper at ICLR 2022
F.3 CELEBA-HO HQ-256 EXTENDED RESULTS
In this section, we provide additional qualitative results on CelebA-HQ-256. For high quality sam- ples using our new SSCS solver Fig. 16.
Samples generated with an adaptive step size Runge-Kutta 4(5) solver at different solver tolerances can be found in Fig. 17. We found that our model still generates very good samples even for solver error tolerance of 10 using an average of 129 neural network evaluations.
Lastly, we show "generation paths" of samples from our CelebA-HQ-255 model: see Fig. 18 and Fig. 19 for samples from the probability flow ODE and the generative SDE, respectively. We visu- alize the continuous generation paths time steps. Interestingly, we can see that the velocity variables "encode" the data at intermediate t. On the other hand, at time = 1.0, by construction, both data and velocity are distributed according to the "equilibrium distribution" of the diffusion, namely, (u) Nxx: Oa, Ia) (v; O2, MTa)- Furthermore, as + 0 the data variables approximately converge to the data distribution, while the velocity variables approximately converge to another Normal distribution (v; Oa, Y Ia) (with y 0.04 in our experiments).
snapshots of data and velocity variables at eight different
Recall that for CLD, the neural network approximates the score Vv, log (v(v, We believe that the generation paths are further evidence that CLD-SGMs need to learn simpler models: for fixed t the velocity variable V, appears to be "noisy" version of the data x,, and therefore we believe Pe(V.(x,) to be relatively smooth and simple when compared to the marginal (x+)..
Finally, note that in Figs. 18 and 19, when visualizing the velocity variables, used colorization scheme that corresponds exactly to the inverse of the color scheme used for visualizing the images themselves. Alternatively, we could also interpret this in such way that we are not actually visu- alizing velocities, but negative velocities with flipped signs. When using this inverse colorization scheme for the velocities, we see that at intermediate where the velocities encode the data, the color values visualizing image data and velocities are, apart from the additional noise in the veloc- ities, similar (i.e. the velocities appear as noisy versions of the actual images). This implies that image pixel values x, translate into corresponding negative velocities back towards the mean of the equilibrium distribution. This consequence of the Hamiltonian coupling between the data and velocity variables. In other words, it is result of the negative sign front of x, the H term in Eq, (5) (and analogously for the reverse-time generative SDE), Also the visualizations
that pull the pixel values
project page (https: //tp///vtaaaaaaaa//e gi thub 1o/ CCM)
45

---

Published as a conference paper at ICLR 2022
IE
U
r
Figure 16: Samples generated by our model on the CelebA-HQ-256 dataset using our SSCS solver.
46

---

Published as a conference paper at ICLR 2022
(a) ODE solver error tolerance 10 9; 273 average NFE.
(b) ODE solver error tolerance 10-4; 190 average NFE.
(c) ODE solver error tolerance 10 3, 129 average NFE.
I
(d) ODE solver tolerance 10-2; 99.4 average NFE.
Figure 17: Samples generated by our model on the HelebA-HQ-255 256 dataset using Runge-Kutte 4(5) adaptive ODE solver solve the probability flow ODE. We show the effect of the ODE solver emor tolerance on the quality of samples ((a), (b), (c) and (d) were generated using the same prior samples). Little visual differences can seen between 10-5 and 10-1, Low frequency artifacts can be observed at 10-3 Deterioration starts to set in at 10-2
47

---

Published as a conference paper at ICLR 2022
Figure 18: Generation paths samples from our CelebA-HQ-2566 model (Runge-Kutta 4(5) solver;
mean NFE: 288). Odd and even rows visualize data and velocity variables, respectively. The eight columns correspond to times t e { {1.0, 0.5, 0.3, 0.2, 0.1, 10-2, 2 10-3, 3 10 53 (from left to right). The velocity distribution converges
Normal (different variances) for both t 0 and t 1. See App. F.3 for visualization details and discussion.
48

---

Published as a conference paper at ICLR 2022
Figure 19: Generation paths of samples from our CelebA-HQ-256 model (SSCS-QS QS using only 150 steps). Odd and even rows visualize data and velocity variables, respectively. The eight columns correspond to times E {1.0, 0.5, 0.3, 0.2, 0.1, 10 1003 10 1-0-. $j (from left to right). The velocity distribution converges to Normal (different variances) for both 0 and + 1. See App. E3 for visualization details and discussion.
49

---

Published as a conference paper at ICLR 2022
G PROOFS OF PERTURBATION KERNELS
In this section, we prove the correctness of the perturbation kernels of the forward diffusion (App. B.I) well as for the analytical splitting term in our SSCS (App. D.2). All derivations are presented for general time-dependent ((t).
G.1 FORWARD DIFFUSION
We have the following ODEs describing the evolution of the mean and the covariance matrix
dyp (f(t) @ La)pe
(110)
dt
d
(f(t) @ La)=; + [(f(t) @ 1)) + (G(()G @ la.
(111)
dt
where
0
48(t)r-2
f(t) =
(112)
-(t)
-3((t)r-1)
G(t) = 0
0
0 V2r
(113)
In App. B.1, we claim the following solutions:
Ht = Cipr
(114)
P
A =
(115)
C. = e -2B(()
(116)
u = (B(()I- -1 X0 + 4B((L)I 3vo + Xo:
(117)
μe = - B(t)xo -B((tr- Vv + Vo.
(118)
and
E. = 2i @ la:
(119)
Σ. = D.E.,
(120)
Exz
Σ. =
(121)
D F
D. = e 48(())
(122)
E EeF AB(() - 1 4B()) - 1)) - 2) 16
(123)
EFV = -B(t)2 + 4B(t)I- -2 -26- - 81)
(124)
D.
12 1 (1 -1) B(t)I - Eg" (1 11+483 41+ -0B())B-))) (t) (E
2).
(125)
where B(() f/ a(6) di and no = =[xo, vol T as well as D5 and 2 are initial conditions.
G.1.1 PROOF OF CORRECTNESS OF THE MEAN
Plugging the claimed solution (Eqs.
(114)-(118) back into the ODE (Eq, 110), we obtain
dCr.
dip C. Ct(f(t) La)
pr: dt dt
(126)
The above can be decomposed into two equations:
dud
-23(t)r- kr +
dt 4,(t)r- k
(127)
duy
-23(/)r-1 1 + di
3(0 43(1))-111) p
(128)
50

---

Published as a conference paper at ICLR 2022
where we used the fact that dC, -2
Eq. (127): Plugging the claimed solution into Eq. (127), obtain:
-28(t)1 [28(1) 1o Vottoo
[200 48 Pvo]
2 [BB(t)xi 2B(e) -
(129)
Eq. (128): After simplification, plugging in the claimed solution into Eq. (128), we obtain:
" - + [=000))6-2- -B(1) 25( xo +
(130)
This completes the proof of the correctness of the mean.
G.1.2 PROOF CORRECTNESS OF THE COVARIANCE
Plugging the claimed solution (Eqs. (119)-(125) back in the ODE (Eq. (111)), we obtain
dE dt
dD: E
Ia = Di(f()) @ (2, @ Ia) + D. [((()) @ (a) & Ta] + (G())G(()) @ Ia.
D.
+
dt
(131)
Noting that
(S(1) 8 L2) @ Ia) = (f((),,) la
(->,- 41
41 25
41 2
(132)
= B(t)
x la.
E
Ex - 41
and
G(t)G (t) (0 2
(133)
well the fact that dD
--((tr Di, decompose Eq. (131) into three equations:
-
-38(()
83(()T
2500 -
(134)
dt
" " + ev -48(t)
3(t)
- - 4T -p-1111 DEV + 4T
-
(135)
dt
ww
43(()r 1 dt
3(t)
22x 8I 81 + 2,)
(136)
Eq. (134): Plugging the claimed solution into Eq. (134), we obtain
48(()1 4B
4444(()) (26 43+ (t) (2 26B (1)
2
+ 48(0)) (0())=-LABI0)) ABOTT 43(()) (20 ~ +8.( 2Ex - -))+323())((1 Sie
= -3(()r-2 (BB))5 -
2B (())= (25 -2) 2)-8B2 vo
(137)
Eq. (135): After simplification, plugging the claimed solution into Eq. (135), we obtain
-00EE 48(t)P=2 48(())((t)) ( (25 -1)33( E
=-B(t) S
LBBrr ABror - -1+4 (S8* - +41 (1)) (25 2)+168"
+33(t) T2 (B()) 1- +B())
(
- + Σ." (144844)13411883 +1²
(S
2]
(138)
51

---

Published as a conference paper at ICLR 2022
Eq. (136): After simplification, plugging the claimed solution into Eq. (136), obtain
18(0r AA -1) + +B Ed
+ 1 18(1) 18101 +3(1IT + DNd (L11B 5 (2) - 2)
AB(()) + 23(1)B(() (E6 -2] = -2)((() [B(t) +4B(A) (t)r-) 2)) 20
(139)
+2T, (B
This completes the proof of the correctness of the covariance.
G.2 ANALYTICAL SPLITTING TERM OF SSCS
We have the following ODEs describing the evolution of the mean and the covariance matrix
dit
(J(T - @ Ia)
(140)
dt
dEc
dl (f(T -) @ 1.) + (U(T- I.). + (G (T- (T--)) @ la. (141)
where
0 -43(T - 1)
f(T -) t)= = +B((T ()
(142)
-43(T - 1)r-1
G(T -) = (0 V2T t)
0
(143)
These ODEs are very similar to the ODEs of the forward diffusion in App. G.1. the only difference being flipped signs in the off-diagonal terms of S(T - 1) (highlighted in red).
In App. D.2, we claim the following solutions
A Cefae
(144)
(145)
C. = -8(())-
(146)
i = 2B(() (xx; -x:()B()) Re'
(147)
3. +B(()xx -B((t) 'v. Ve.
(148)
and
Σ. = Σ₁ @ lai
(149)
Σ, DYEr,
(150)
FEE
(151)
D. 4B(1JT
(152)
FE = /1()T-1 - 4B(t) - (t)1
(153)
Zvv = (1)
(154)
I T (11(1)T+ -1) -1) + B((T - 2B)
(155)
where B(t)= = 3T 2) di and jay = [xx,, Vel is an initial condition. Differences of the above solution to the solutions of the forward diffusion are again highlighted in red. Note that by con- struction the initial covariance for the analytical splitting term of SSCS is the zero matrix, i.c. Saz Eee Dg 0, since always initialize from an "updated sample", which itself does not have any uncertainty. Also note that in this derivation we use general initial t' (whereas in App. G.I set l' for simplicity).
52

---

Published as a conference paper at ICLR 2022
G.2.1 PROOF OF CORRECTNESS OF THE MEAN
Plugging the claimed solution (Eqs.
(144)-(148) into the ODE (Eq. (140)), we obtain
ji dC,
dīt C. = C,U((T - 1) & Ia)ja.
(156)
dt di
The above can be decomposed into two equations:
-2 -)r- fu: + du:
-3 -t] 2
(157)
dt
dit
-3 - -)T i dl
3(T -t)Øū 43(T 1111-11'11" p.;
(158)
where we used the fact that dC. 2T 1)--
dt
Eq. (157): Plugging the claimed solution into Eq. (157), we obtain:
-33(T 1)r'' 2 X -B( Vedk x 28(T f)1 N,' 48(T t)1 ve
(159)
--13(T-t)T" ( 28( vol
Eq. (128): After simplification, plugging the claimed solution into Eq. (128), we obtain:
2 :)r- [B(A)xF- ( Vv
vel BT L) 2B( X 43(()L Vec
S(T +)x/ -2( 1)T
(160)
This completes the proof of the correctness of the mean.
G.2.2 PROOF OF CORRECTNESS OF THE COVARIANCE
Plugging the claimed solution (Eqs. (149)-(155) into the ODE (Eq. (141)), we obtain
ds. dD x de
D.(f @ L)) La) + D. [If La) @ Ia)]
(GG") @ (161)
D.+
dt
with f = -[[T - t) and = G(T t).
Noting that
(f(T - ( @ la) = (f(T - t)Ei)) @ Ia
= B -t) 5ix -4 - -5
41 -4 25
(162)
8 la
4
5
40
5w0
and
0 - -1) -1) 3T 2 :)):
G(T - -)G (T-
(163)
as well as the fact
dD: d5
- 43 -) Dr, can decompose Eq. (161) into three equations:
di
-43(T t)r- 15,0 +
-8 -)
(164)
dt
-43(T
E. 15
5:00 3(T 1) 1Sa² 41 11500
-2
(165)
dl
dove
-43(T - - Σ.. +
3(T - t) |2S: 8I 8 Soo
+ (T -
(166)
di
Eq. (164): Plugging the claimed solution into Eq. (164), we obtain
4B(T -) 4
-
48(T- -
-48(T--)) - - -1630-0--01--
(167)
--8(T -|
-AB2
53

---

Published as a conference paper at ICLR 2022
Eq. (165): After simplification, plugging the claimed solution into Eq. (165), obtain
-88(T-t)B(t)1
3(' ( 1-4B(t1 1-AB(0T-8BB"(11 -8B"
(168)
3 - [ ABIOTT 1) -) (1)
Eq. (166): After simplification, plugging the claimed solution into Eq. (166), obtain
4 1)
[ the 4B(ETT
tBP tjT -88(T - tt)(())
(169)
-3 18 (t) +2 t)e 1801511
This completes the proof of the correctness of the covariance.
To connect back the SSCS presented in App. D.2, recall that in practice we use constant | (and T = 1) and that we solve for small time steps of size 2 expressions presented App. D.2.
5t such that (t()) = B which leads to the
54


Reviews:
[{"review_id": "lZim7VZM-MH", "reviewer": "Reviewer_AJW8", "summary": "", "questions": "", "limitations": "", "rating": 10, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Strengths**\n\n- Section 1 and 4 provide a very thorough overview of recent work, with a fantastic visualisation of the method.\n- I found the arguments and derivations throughout Section 3 compelling, with the appropriate mathematics reserved for the appendix. The synthetic experiment demonstrating the ease of numerical approximation is also clearly explained and useful.\n- Particular praise should be afforded to the consistent use of appropriate references in Section 3, placing the derivations in context, as well as the connections drawn to the high-friction limit in Section A.2.\n- Model evaluation is in-line with literature and on-par with the very best performing models. Up to extreme levels of compute, the method is essentially SoTA. Given the higher-order nature of the sampling, evaluation of NFEs is also well-motivated and is a clear improvement for constant compute. Finally, the robust hyperparameters are of practical relevance.\n\n**Weaknesses**\n\nThere are very few weaknesses beyond typos and personal preferences for presentation.\n\nMy only comment would be the absence of references to the line of work on higher-order MCMC by Michael I Jordan and others which I believe to be pertinent:\n\n[1] Ma, Yi-An, et al. \"Is there an analog of Nesterov acceleration for MCMC?.\" arXiv preprint arXiv:1902.00996 (2019).\n\n[2] Mou, Wenlong, et al. \"High-order Langevin diffusion yields an accelerated MCMC algorithm.\" arXiv preprint arXiv:1908.10859 (2019).\n\nMinor comments:\n- The reference for Song et al.'s SDE paper is dated as ICLR 2020, rather than 2021.\n- Penultimate para of sec 1: \"diffusionand\" -> \"diffusion and\"\n- Section B.2: \"SGMs\" does not need to be repeated in sentence 1.\n- One of the $\\Sigma_{t}^{xv}$ matrices in Equation (9) should $\\Sigma_{t}^{vt}$? This notation is used at the end of Section B.2, so consistency would make things a little smoother.\n- Sec 5.1: \"We are significantly outperforming this model\" could be rephrased.\n- Sec 5.1 \"outperforms\" -> \"outperform\"\n- Sec 5.2, para 2: \"as backbone\" -> \"as a/the backbone\"\n- Tab. 2 caption: \"denotes\" -> \"denoted\"\n- A.1: \"very inefficient too\" -> \"very inefficiently too\"\n- Sec B2, penultimate para. and Equation (67): What is $\\Sigma_{t}^{zz}$? I assume this is a typo?", "summary_of_the_paper": "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.", "summary_of_the_review": "This paper is exceptionally well put together, in terms of derivations, presentation, and experimentation. The narrative is well-written, well-motivated, and logical throughout. To the best of my abilities, the proofs are entirely correct. The only way I can see to improve the experiments is to offer the authors more compute or data.\n\nAs a result, I believe this paper to be of the highest quality and recommend that it is given particular attention at the conference.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "i_p47Jp0wCP", "reviewer": "Reviewer_tK3A", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper claims the following contributions.\n1. First, the paper proposes a novel class of continuous-time diffusion-based generative models, called critically-damped Langevine diffusion (CLD).\n2. Second, the authors show that a score matching objective for the proposed method requires only $\\nabla_{v_t} \\log p_t(v_t | x_t)$, not $\\nabla_{(x_t, v_t)} p_t(x_t, v_t)$.\n3. Third, for the proposed model, the paper proposes a variant of denoising score matching, called hybrid score matching (HSM), observing that the training of CLD-based models is unstable by using DSM.\n4. Fourth, the authors propose a leapfrog-style integrator designed for CLD-based models, benefitting from Hamiltonian dynamics' symplectic structure.\n5. Finally, the paper claims to provide novel insights into continuous-time diffusion-based generative models and their connections to statistic mechanics.\n\n\n**Strengths of the paper**\n\nIn general, the paper's contributions wrt the novelty are clear for several reasons:\n- First, the paper provides novel diffusion-based models inspired by stochastic Hamiltonian dynamics. Moreover, the paper provides sufficient proofs to support that the proposed methods are well defined. I find that the proposed method will be very interesting contribution to the generative model community.\n- Second, the authors demonstrate the scalability of the proposed methods by (1) showing that the increment of the model complexity is marginal thanks to reverse-time SDEs only requires log-density gradient of conditional distributions of velocities given data, and (2) proposing a modified DSM loss suited for CLD-based models.\n    \nIn addition, I found that the paper has a well-organized structure so that it is clear to understand the proposed and other practical techniques to improve training.\n\n**Weaknesses of the paper**  \nThe following aspects of the paper can be improved.\n\n- First, the discussion about the proposed integrator can be improved. The authors emphasize that the proposed leapfrog-style integrator for CLD-based models can generate better quality samples even when the number of discretization steps is small compared to the Euler-Maruyama method. However, it isn't easy to assess the statement based on the results described in Tables 2 and 3. In the experiments, the results are dependent on both the integrators and the qualities of learned models. To resolve this issue, one can add toy experiments to solidify the argument, where ground truth reverse-time SDEs are (approximately) accessible; thus, the resulting analysis can highlight the difference of the integrators.\n\n- Second, the paper claims that the proposed method only needs to learn smoother signals, and thus the models can perform better. However, in the current submission, the discussion about this seems limited.  \n    \n  While the authors have already discussed the property by using toy experiments (such as in Fig. 2), the benefits of the proposed model seem marginal in high-dimensional datasets.  \n    \n  I assume that dropping weights to improve the generation qualities has diluted the benefits of the proposed method. Thus I believe that discussion about the proposed methods wrt maximum likelihood training can strengthen the significance of the proposed method. In my understanding, dropping weights for training diffusion-based generative models prevents modeling unbounded scores at $t$'s close to the data, where the weights are extremely high. As a result, the previous methods don't learn unbounded scores in practice, which can be a potential reason for the marginal benefits of the proposed method in high-dimensional settings. Therefore, the proposed method may benefit from not dropping weights, possibly resulting in high-quality samples being achieved with maximum likelihoods training. I anticipate that the comparison can be more observable if the models are trained in real space after logit transformation + uniform dequantization.  \n  \n  Moreover, it can be helpful to add discussions about the trade-off between the variance of the DSM/HSM objective estimates and the variance of the $p_0(v_0)$. For example, I imagine that the high variance of $p_0(v_0)$ results in the high variance in $p_t(v_t|x_t)$, implying that the models only need to learn smoother scores. Again, I assume this aspect can be more observable if the models are trained to maximize likelihoods.\n    \n- Third, the motivation of Hamiltonian dynamics-style generative models can be improved. Aside from the applicability of the leapfrog-style integrator for the proposed method, the authors claim that the velocity in CLDs accelerates mixing in the diffusion process (for example, in Sec 4). However, it is unclear how important is this \"mixing\" in the context of diffusion-based generative models.  \n    In standard MCMC methods, the mixing properties are important since such methods aim at matching the distribution of recorded states (from a generated chain) to target distributions. On the other hand, the diffusion-based generative models rely on the distributions of collection of chains, esp. the chains' last states.  \nKnowing that the previous diffusion-based generative models are already expressive and capable of generating diverse samples, it is difficult to motivate the proposed method while its novelty.   \n    I believe that improving the discussion about the proposed integrators will resolve the issue, regardless of the proposed method's connection to statistic mechanics.\n    \n- Finally, the discussion about the motivation for HSM can be improved. In Sec 5.3, the authors state that for the proposed method, training with regular DSM is unstable. While the authors state some reasoning about the instability, but the analysis about it is limited. I believe that the single sample Monte Carlo estimate for the regular DSM objective can have a high variance for the proposed method, and thus I believe that the marginalizing out $v_0$ in HSM has reduced the variance. I consider that additional toy experiments that discuss the difference between DSM vs. HSM wrt their variance will further motivate the importance of HSM for CLD-based models.\n", "summary_of_the_paper": "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.\n\nIn 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.\n\nThe 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.\n\nIn 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.\n\nIn 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.\n\nIn 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.", "summary_of_the_review": "In general, the paper's contributions wrt the novelty are clear, and the proposed methods are well-defined. In addition, I found that the paper has a well-organized structure so that it is clear to understand the proposed and other practical techniques to improve training. Thus, I'm inclined to accept the paper. However, I found that several aspects of the submission can be improved. I'm inclined to improve my evaluation if the aforementioned weak points are well addressed.  \n\n  \n===== POST-REBUTTAL COMMENTS ======== \n\nThe rebuttal had addressed most of my concerns. Consequently, I raised my score from 6 to 8, and I also raised the \"Empirical Novelty And Significance\" score from 2 to 3.\n\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "fa8cgh27ZYv", "reviewer": "Reviewer_F5kX", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths of the paper:\n- The paper uses tools from physics to design a novel score based generative model, allowing it to leverage existing insights in that field\n- The paper results in training of generative networks that outperform existing approaches when making sure to balance compute budgets and model size\n- The paper derive a SDE integrator that allows for efficient sampling from the proposed class of models.\n\n\nMy only question to the authors is: Do they envision that a class of methods of this sort are also possible?  For example, would the method work if they had an acceleration term in which the noise was injected?\n", "summary_of_the_paper": "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.  \n", "summary_of_the_review": "I recommend the paper be accepted.  It presents a potentially significant improvement to the training of score-based generative models which is well grounded in theory from physics and demonstrates strong empirical performance for sample generation when compared to baselines of similar complexity and compute budgets. \n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "bHKuKHnTrxz", "reviewer": "Reviewer_95SQ", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This is a paper well-done. Diffusion models require designing (or sometimes training) an inference process that brings the data distribution closer to some tractable prior. Working on an augmented state-space to introduce a momentum (velocity) variable to simulate damped Langevin dynamics is very a natural extension that complements well the current score-based diffusion-based models literature. The authors then went on exploring several benefits of such an extended framework, including a smoother (conditional) score function to approximate, a novel smoother (score-matching) loss function to optimize, and a specialized numerical integrator for efficient sampling. I believe the novelties in all different aspects introduced in this work could largely benefit the community.  \n\nOverall, this is a great paper that is easy to follow and has multiple neat ideas. Strengths of the paper include \n* Novelty: although it is a straightforward extension, the authors explored different aspects of the central idea of augmentation/acceleration and proposed a new loss that is smooth and a numerical method that works better for certain settings.\n* Clarity: the paper is well written, ideas presented with clear theoretical development. Most of the time the design choices are explained, leaving room for future research on some relevant topics (such as optimizing for NLL, and adaptive solvers).  \n* Impact: the problem the proposed method solves is very specialized, particularly for SGDMs. The idea of state-space augmentation could potentially mitigate the exploding score of DSM-type training and can be a principled way of designing an SGDM that naturally admits faster numerical simulation.\n  \nI do not have any serious criticisms, but I do have the following (more detailed) questions / comments:\n1. There is a paragraph in Section 3 describing the trade-off between being overdamped and underdamped. The explanation is quite intuitive, but is there any more quantitative evidence for opting for critical damping (either in theory or empirically)? I suppose an optimal ratio between $\\Gamma$ and $M$ will be problem-specific. Would such generality be useful in practice?\n2. A missing point about using damped Langevin dynamics is the potential speedup in convergence to the prior distribution (A la Ma et al. 2019, Proposition 1). Diffusing along with the momentum variable could allow for significant acceleration, which might reduce the integration time T needed for sampling as well. \n3. For HSM, is marginalizing out $v_0$ the key to smoothing out the loss function (compared to DSM)? If so, is this a form of Blackwellization, and does the gradient estimate always have a smaller variance compared to DSM?\n4. In 3.3, Leimkuhler & Reich, 2005 was cited for Euler’s methods not being suitable for Hamiltonian dynamics. Is there any intuition behind this? Can the standard midpoint integration scheme be employed which is well suited for Hamiltonian systems?\n5. For SSCS, can the second half step of the current iteration be merged with the first half step of the next iteration, similarly to leap-frog integration. \n6. Experiment-wise, are the numbers presented in the tables with one random seed (for experiment and for evaluation)? Just reading the number, I wouldn’t say the difference between EM and SSCS for large n in Table 3 is negligible, but perhaps that’s more due to the stochasticity of the nature of FID? \n7. About $p(v_0)$, the reason for worse NLL for small $\\gamma$ could be due to the fact that the upper bound requires subtracting the entropy of the augmented distribution (which will be smaller). It’d be worth noting that this could be potentially optimized for NLL. Aside from that, smaller $\\gamma$ would also mean the magnitude for the score at $t\\rightarrow0$ is also larger? \n\nMa et al 2019: Is There an Analog of Nesterov Acceleration for MCMC?\n", "summary_of_the_paper": "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. \n", "summary_of_the_review": "This is a timely paper that explores acceleration methods for diffusion-based models. The authors also explored different aspects / benefits of the proposed framework which I think will benefit the community of researchers working on SGDMs, including the smoothness of loss function and a new numerical integration method. \n\nMy recommendation is on the lower band of 8: accept, good paper. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper develops a diffusion-process based generative model that perturbs the data using a critically damped Langevin diffusion. The diffusion is set up through an auxiliary velocity term like in Hamiltonian dynamics. The idea is that picking a process that diffuses faster will lead to better results.The paper then constructs a new score matching objective adapted to this diffusion, along with a sampling scheme for critically damped Langevin score based generative models. The idea of a faster diffusion to make generative models is a good one. The paper is a solid accept.  

Reviewer tK3A was lukewarm as evidenced by their original 2 for empirical novelty that moved to a 3. From my look, it felt like a straightforward application of ideas in one domain, sampling, to another, generative modeling. It's a good paper, but it does not stand out relative to other accepts.
Decision: Accept (Spotlight)