Score-Based Generative Modeling through Stochastic Differential Equations

Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, Ben Poole

Introduction

Two successful classes of probabilistic generative models involve sequentially corrupting training data with slowly increasing noise, and then learning to reverse this corruption in order to form a generative model of the data. Score matching with Langevin dynamics (SMLD) (Song & Ermon, 2019) estimates the score (i.e., the gradient of the log probability density with respect to data) at each noise scale, and then uses Langevin dynamics to sample from a sequence of decreasing noise scales during generation. Denoising diffusion probabilistic modeling (DDPM) (Sohl-Dickstein et al., 2015; Ho et al., 2020) trains a sequence of probabilistic models to reverse each step of the noise corruption, using knowledge of the functional form of the reverse distributions to make training tractable. For continuous state spaces, the DDPM training objective implicitly computes scores at each noise scale. We therefore refer to these two model classes together as score-based generative models.

Score-based generative models, and related techniques (Bordes et al., 2017; Goyal et al., 2017; Du & Mordatch, 2019), have proven effective at generation of images (Song & Ermon, 2019; 2020; Ho et al., 2020), audio (Chen et al., 2020; Kong et al., 2020), graphs (Niu et al., 2020), and shapes (Cai et al., 2020). To enable new sampling methods and further extend the capabilities of score-based generative models, we propose a unified framework that generalizes previous approaches through the lens of stochastic differential equations (SDEs).

Specifically, instead of perturbing data with a finite number of noise distributions, we consider a continuum of distributions that evolve over time according to a diffusion process. This process progressively diffuses a data point into random noise, and is given by a prescribed SDE that does not depend on the data and has no trainable parameters. By reversing this process, we can smoothly mold random noise into data for sample generation. Crucially, this reverse process satisfies a reverse-time SDE (Anderson, 1982), which can be derived from the forward SDE given the score of the marginal probability densities as a function of time. We can therefore approximate the reverse-time SDE by training a time-dependent neural network to estimate the scores, and then produce samples using numerical SDE solvers. Our key idea is summarized in Fig. 1.

Our proposed framework has several theoretical and practical contributions:

Flexible sampling and likelihood computation: We can employ any general-purpose SDE solver to integrate the reverse-time SDE for sampling. In addition, we propose two special methods not viable for general SDEs: (i) Predictor-Corrector (PC) samplers that combine numerical SDE solvers with score-based MCMC approaches, such as Langevin MCMC (Parisi, 1981) and HMC (Neal et al., 2011); and (ii) deterministic samplers based on the probability flow ordinary differential equation (ODE). The former unifies and improves over existing sampling methods for score-based models. The latter allows for fast adaptive sampling via black-box ODE solvers, flexible data manipulation via latent codes, a uniquely identifiable encoding, and notably, exact likelihood computation.

Controllable generation: We can modulate the generation process by conditioning on information not available during training, because the conditional reverse-time SDE can be efficiently estimated from unconditional scores. This enables applications such as class-conditional generation, image inpainting, colorization and other inverse problems, all achievable using a single unconditional score-based model without re-training.

Unified framework: Our framework provides a unified way to explore and tune various SDEs for improving score-based generative models. The methods of SMLD and DDPM can be amalgamated into our framework as discretizations of two separate SDEs. Although DDPM (Ho et al., 2020) was recently reported to achieve higher sample quality than SMLD (Song & Ermon, 2019; 2020), we show that with better architectures and new sampling algorithms allowed by our framework, the latter can catch up—it achieves new state-of-the-art Inception score (9.89) and FID score (2.20) on CIFAR-10, as well as high-fidelity generation of 1024×10241024\times 1024 images for the first time from a score-based model. In addition, we propose a new SDE under our framework that achieves a likelihood value of 2.99 bits/dim on uniformly dequantized CIFAR-10 images, setting a new record on this task.

Background

Given sufficient data and model capacity, the optimal score-based model sθ(x,σ)\mathbf{s}_{{\boldsymbol{\theta}}^{*}}(\mathbf{x},\sigma) matches xlogpσ(x)\nabla_{\mathbf{x}}\log p_{\sigma}(\mathbf{x}) almost everywhere for σ{σi}i=1N\sigma\in\{\sigma_{i}\}_{i=1}^{N}. For sampling, Song & Ermon (2019) run MM steps of Langevin MCMC to get a sample for each pσi(x)p_{\sigma_{i}}(\mathbf{x}) sequentially:

2 Denoising diffusion probabilistic models (DDPM)

After solving Eq. 3 to get the optimal model sθ(x,i)\mathbf{s}_{{\boldsymbol{\theta}}^{*}}(\mathbf{x},i), samples can be generated by starting from xNN(0,I)\mathbf{x}_{N}\sim\mathcal{N}(\mathbf{0},\mathbf{I}) and following the estimated reverse Markov chain as below

Score-based generative modeling with SDEs

Perturbing data with multiple noise scales is key to the success of previous methods. We propose to generalize this idea further to an infinite number of noise scales, such that perturbed data distributions evolve according to an SDE as the noise intensifies. An overview of our framework is given in Fig. 2.

Our goal is to construct a diffusion process {x(t)}t=0T\{\mathbf{x}(t)\}_{t=0}^{T} indexed by a continuous time variable t[0,T]t\in[0,T], such that x(0)p0\mathbf{x}(0)\sim p_{0}, for which we have a dataset of i.i.d. samples, and x(T)pT\mathbf{x}(T)\sim p_{T}, for which we have a tractable form to generate samples efficiently. In other words, p0p_{0} is the data distribution and pTp_{T} is the prior distribution. This diffusion process can be modeled as the solution to an Itô SDE:

Typically, pTp_{T} is an unstructured prior distribution that contains no information of p0p_{0}, such as a Gaussian distribution with fixed mean and variance. There are various ways of designing the SDE in Eq. 5 such that it diffuses the data distribution into a fixed prior distribution. We provide several examples later in Section 3.4 that are derived from continuous generalizations of SMLD and DDPM.

2 Generating samples by reversing the SDE

By starting from samples of x(T)pT\mathbf{x}(T)\sim p_{T} and reversing the process, we can obtain samples x(0)p0\mathbf{x}(0)\sim p_{0}. A remarkable result from Anderson (1982) states that the reverse of a diffusion process is also a diffusion process, running backwards in time and given by the reverse-time SDE:

3 Estimating scores for the SDE

The score of a distribution can be estimated by training a score-based model on samples with score matching (Hyvärinen, 2005; Song et al., 2019a). To estimate xlogpt(x)\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}), we can train a time-dependent score-based model sθ(x,t)\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{x},t) via a continuous generalization to Eqs. 1 and 3:

We typically need to know the transition kernel p0t(x(t)x(0))p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)) to efficiently solve Eq. 7. When f(,t)\mathbf{f}(\cdot,t) is affine, the transition kernel is always a Gaussian distribution, where the mean and variance are often known in closed-forms and can be obtained with standard techniques (see Section 5.5 in Särkkä & Solin (2019)). For more general SDEs, we may solve Kolmogorov’s forward equation (Øksendal, 2003) to obtain p0t(x(t)x(0))p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)). Alternatively, we can simulate the SDE to sample from p0t(x(t)x(0))p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)) and replace denoising score matching in Eq. 7 with sliced score matching for model training, which bypasses the computation of x(t)logp0t(x(t)x(0))\nabla_{\mathbf{x}(t)}\log p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)) (see Appendix A).

4 Examples: VE, VP SDEs and beyond

The noise perturbations used in SMLD and DDPM can be regarded as discretizations of two different SDEs. Below we provide a brief discussion and relegate more details to Appendix B.

When using a total of NN noise scales, each perturbation kernel pσi(xx0)p_{\sigma_{i}}(\mathbf{x}\mid\mathbf{x}_{0}) of SMLD corresponds to the distribution of xi\mathbf{x}_{i} in the following Markov chain:

where zi1N(0,I)\mathbf{z}_{i-1}\sim\mathcal{N}(\mathbf{0},\mathbf{I}), and we have introduced σ0=0\sigma_{0}=0 to simplify the notation. In the limit of NN\to\infty, {σi}i=1N\{\sigma_{i}\}_{i=1}^{N} becomes a function σ(t)\sigma(t), zi\mathbf{z}_{i} becomes z(t)\mathbf{z}(t), and the Markov chain {xi}i=1N\{\mathbf{x}_{i}\}_{i=1}^{N} becomes a continuous stochastic process {x(t)}t=01\{\mathbf{x}(t)\}_{t=0}^{1}, where we have used a continuous time variable tt\in for indexing, rather than an integer ii. The process {x(t)}t=01\{\mathbf{x}(t)\}_{t=0}^{1} is given by the following SDE

Likewise for the perturbation kernels {pαi(xx0)}i=1N\{p_{\alpha_{i}}(\mathbf{x}\mid\mathbf{x}_{0})\}_{i=1}^{N} of DDPM, the discrete Markov chain is

As NN\to\infty, Eq. 10 converges to the following SDE,

Therefore, the noise perturbations used in SMLD and DDPM correspond to discretizations of SDEs Eqs. 9 and 11. Interestingly, the SDE of Eq. 9 always gives a process with exploding variance when tt\to\infty, whilst the SDE of Eq. 11 yields a process with a fixed variance of one when the initial distribution has unit variance (proof in Appendix B). Due to this difference, we hereafter refer to Eq. 9 as the Variance Exploding (VE) SDE, and Eq. 11 the Variance Preserving (VP) SDE.

Inspired by the VP SDE, we propose a new type of SDEs which perform particularly well on likelihoods (see Section 4.3), given by

When using the same β(t)\beta(t) and starting from the same initial distribution, the variance of the stochastic process induced by Eq. 12 is always bounded by the VP SDE at every intermediate time step (proof in Appendix B). For this reason, we name Eq. 12 the sub-VP SDE.

Since VE, VP and sub-VP SDEs all have affine drift coefficients, their perturbation kernels p0t(x(t)x(0))p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)) are all Gaussian and can be computed in closed-forms, as discussed in Section 3.3. This makes training with Eq. 7 particularly efficient.

Solving the reverse SDE

After training a time-dependent score-based model sθ\mathbf{s}_{\boldsymbol{\theta}}, we can use it to construct the reverse-time SDE and then simulate it with numerical approaches to generate samples from p0p_{0}.

Numerical solvers provide approximate trajectories from SDEs. Many general-purpose numerical methods exist for solving SDEs, such as Euler-Maruyama and stochastic Runge-Kutta methods (Kloeden & Platen, 2013), which correspond to different discretizations of the stochastic dynamics. We can apply any of them to the reverse-time SDE for sample generation.

Ancestral sampling, the sampling method of DDPM (Eq. 4), actually corresponds to one special discretization of the reverse-time VP SDE (Eq. 11) (see Appendix E). Deriving the ancestral sampling rules for new SDEs, however, can be non-trivial. To remedy this, we propose reverse diffusion samplers (details in Appendix E), which discretize the reverse-time SDE in the same way as the forward one, and thus can be readily derived given the forward discretization. As shown in Table 1, reverse diffusion samplers perform slightly better than ancestral sampling for both SMLD and DDPM models on CIFAR-10 (DDPM-type ancestral sampling is also applicable to SMLD models, see Appendix F.)

2 Predictor-corrector samplers

Unlike generic SDEs, we have additional information that can be used to improve solutions. Since we have a score-based model sθ(x,t)xlogpt(x)\mathbf{s}_{{\boldsymbol{\theta}}^{*}}(\mathbf{x},t)\approx\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}), we can employ score-based MCMC approaches, such as Langevin MCMC (Parisi, 1981; Grenander & Miller, 1994) or HMC (Neal et al., 2011) to sample from ptp_{t} directly, and correct the solution of a numerical SDE solver.

Specifically, at each time step, the numerical SDE solver first gives an estimate of the sample at the next time step, playing the role of a “predictor”. Then, the score-based MCMC approach corrects the marginal distribution of the estimated sample, playing the role of a “corrector”. The idea is analogous to Predictor-Corrector methods, a family of numerical continuation techniques for solving systems of equations (Allgower & Georg, 2012), and we similarly name our hybrid sampling algorithms Predictor-Corrector (PC) samplers. Please find pseudo-code and a complete description in Appendix G. PC samplers generalize the original sampling methods of SMLD and DDPM: the former uses an identity function as the predictor and annealed Langevin dynamics as the corrector, while the latter uses ancestral sampling as the predictor and identity as the corrector.

We test PC samplers on SMLD and DDPM models (see Algorithms 2 and 3 in Appendix G) trained with original discrete objectives given by Eqs. 1 and 3. This exhibits the compatibility of PC samplers to score-based models trained with a fixed number of noise scales. We summarize the performance of different samplers in Table 1, where probability flow is a predictor to be discussed in Section 4.3. Detailed experimental settings and additional results are given in Appendix G. We observe that our reverse diffusion sampler always outperform ancestral sampling, and corrector-only methods (C2000) perform worse than other competitors (P2000, PC1000) with the same computation (In fact, we need way more corrector steps per noise scale, and thus more computation, to match the performance of other samplers.) For all predictors, adding one corrector step for each predictor step (PC1000) doubles computation but always improves sample quality (against P1000). Moreover, it is typically better than doubling the number of predictor steps without adding a corrector (P2000), where we have to interpolate between noise scales in an ad hoc manner (detailed in Appendix G) for SMLD/DDPM models. In Fig. 9 (Appendix G), we additionally provide qualitative comparison for models trained with the continuous objective Eq. 7 on 256×256256\times 256 LSUN images and the VE SDE, where PC samplers clearly surpass predictor-only samplers under comparable computation, when using a proper number of corrector steps.

3 Probability flow and connection to neural ODEs

Score-based models enable another numerical method for solving the reverse-time SDE. For all diffusion processes, there exists a corresponding deterministic process whose trajectories share the same marginal probability densities {pt(x)}t=0T\{p_{t}(\mathbf{x})\}_{t=0}^{T} as the SDE. This deterministic process satisfies an ODE (more details in Section D.1):

which can be determined from the SDE once scores are known. We name the ODE in Eq. 13 the probability flow ODE. When the score function is approximated by the time-dependent score-based model, which is typically a neural network, this is an example of a neural ODE (Chen et al., 2018).

Exact likelihood computation Leveraging the connection to neural ODEs, we can compute the density defined by Eq. 13 via the instantaneous change of variables formula (Chen et al., 2018). This allows us to compute the exact likelihood on any input data (details in Section D.2). As an example, we report negative log-likelihoods (NLLs) measured in bits/dim on the CIFAR-10 dataset in Table 2. We compute log-likelihoods on uniformly dequantized data, and only compare to models evaluated in the same way (omitting models evaluated with variational dequantization (Ho et al., 2019) or discrete data), except for DDPM (LL/LsimpleL_{\text{simple}}) whose ELBO values (annotated with *) are reported on discrete data. Main results: (i) For the same DDPM model in Ho et al. (2020), we obtain better bits/dim than ELBO, since our likelihoods are exact; (ii) Using the same architecture, we trained another DDPM model with the continuous objective in Eq. 7 (i.e., DDPM cont.), which further improves the likelihood; (iii) With sub-VP SDEs, we always get higher likelihoods compared to VP SDEs; (iv) With improved architecture (i.e., DDPM++ cont., details in Section 4.4) and the sub-VP SDE, we can set a new record bits/dim of 2.99 on uniformly dequantized CIFAR-10 even without maximum likelihood training.

Manipulating latent representations By integrating Eq. 13, we can encode any datapoint x(0)\mathbf{x}(0) into a latent space x(T)\mathbf{x}(T). Decoding can be achieved by integrating a corresponding ODE for the reverse-time SDE. As is done with other invertible models such as neural ODEs and normalizing flows (Dinh et al., 2016; Kingma & Dhariwal, 2018), we can manipulate this latent representation for image editing, such as interpolation, and temperature scaling (see Fig. 3 and Section D.4).

Uniquely identifiable encoding Unlike most current invertible models, our encoding is uniquely identifiable, meaning that with sufficient training data, model capacity, and optimization accuracy, the encoding for an input is uniquely determined by the data distribution (Roeder et al., 2020). This is because our forward SDE, Eq. 5, has no trainable parameters, and its associated probability flow ODE, Eq. 13, provides the same trajectories given perfectly estimated scores. We provide additional empirical verification on this property in Section D.5.

Efficient sampling As with neural ODEs, we can sample x(0)p0\mathbf{x}(0)\sim p_{0} by solving Eq. 13 from different final conditions x(T)pT\mathbf{x}(T)\sim p_{T}. Using a fixed discretization strategy we can generate competitive samples, especially when used in conjuction with correctors (Table 1, “probability flow sampler”, details in Section D.3). Using a black-box ODE solver (Dormand & Prince, 1980) not only produces high quality samples (Table 2, details in Section D.4), but also allows us to explicitly trade-off accuracy for efficiency. With a larger error tolerance, the number of function evaluations can be reduced by over 90%90\% without affecting the visual quality of samples (Fig. 3).

4 Architecture improvements

We explore several new architecture designs for score-based models using both VE and VP SDEs (details in Appendix H), where we train models with the same discrete objectives as in SMLD/DDPM. We directly transfer the architectures for VP SDEs to sub-VP SDEs due to their similarity. Our optimal architecture for the VE SDE, named NCSN++, achieves an FID of 2.45 on CIFAR-10 with PC samplers, while our optimal architecture for the VP SDE, called DDPM++, achieves 2.78.

By switching to the continuous training objective in Eq. 7, and increasing the network depth, we can further improve sample quality for all models. The resulting architectures are denoted as NCSN++ cont. and DDPM++ cont. in Table 3 for VE and VP/sub-VP SDEs respectively. Results reported in Table 3 are for the checkpoint with the smallest FID over the course of training, where samples are generated with PC samplers. In contrast, FID scores and NLL values in Table 2 are reported for the last training checkpoint, and samples are obtained with black-box ODE solvers. As shown in Table 3, VE SDEs typically provide better sample quality than VP/sub-VP SDEs, but we also empirically observe that their likelihoods are worse than VP/sub-VP SDE counterparts. This indicates that practitioners likely need to experiment with different SDEs for varying domains and architectures.

Our best model for sample quality, NCSN++ cont. (deep, VE), doubles the network depth and sets new records for both inception score and FID on unconditional generation for CIFAR-10. Surprisingly, we can achieve better FID than the previous best conditional generative model without requiring labeled data. With all improvements together, we also obtain the first set of high-fidelity samples on CelebA-HQ 1024×10241024\times 1024 from score-based models (see Section H.3). Our best model for likelihoods, DDPM++ cont. (deep, sub-VP), similarly doubles the network depth and achieves a log-likelihood of 2.99 bits/dim with the continuous objective in Eq. 7. To our best knowledge, this is the highest likelihood on uniformly dequantized CIFAR-10.

Controllable generation

The continuous structure of our framework allows us to not only produce data samples from p0p_{0}, but also from p0(x(0)y)p_{0}(\mathbf{x}(0)\mid\mathbf{y}) if pt(yx(t))p_{t}(\mathbf{y}\mid\mathbf{x}(t)) is known. Given a forward SDE as in Eq. 5, we can sample from pt(x(t)y)p_{t}(\mathbf{x}(t)\mid\mathbf{y}) by starting from pT(x(T)y)p_{T}(\mathbf{x}(T)\mid\mathbf{y}) and solving a conditional reverse-time SDE:

In general, we can use Eq. 14 to solve a large family of inverse problems with score-based generative models, once given an estimate of the gradient of the forward process, xlogpt(yx(t))\nabla_{\mathbf{x}}\log p_{t}(\mathbf{y}\mid\mathbf{x}(t)). In some cases, it is possible to train a separate model to learn the forward process logpt(yx(t))\log p_{t}(\mathbf{y}\mid\mathbf{x}(t)) and compute its gradient. Otherwise, we may estimate the gradient with heuristics and domain knowledge. In Section I.4, we provide a broadly applicable method for obtaining such an estimate without the need of training auxiliary models.

We consider three applications of controllable generation with this approach: class-conditional generation, image imputation and colorization. When y\mathbf{y} represents class labels, we can train a time-dependent classifier pt(yx(t))p_{t}(\mathbf{y}\mid\mathbf{x}(t)) for class-conditional sampling. Since the forward SDE is tractable, we can easily create training data (x(t),y)(\mathbf{x}(t),\mathbf{y}) for the time-dependent classifier by first sampling (x(0),y)(\mathbf{x}(0),\mathbf{y}) from a dataset, and then sampling x(t)p0t(x(t)x(0))\mathbf{x}(t)\sim p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)). Afterwards, we may employ a mixture of cross-entropy losses over different time steps, like Eq. 7, to train the time-dependent classifier pt(yx(t))p_{t}(\mathbf{y}\mid\mathbf{x}(t)). We provide class-conditional CIFAR-10 samples in Fig. 4 (left), and relegate more details and results to Appendix I.

Imputation is a special case of conditional sampling. Suppose we have an incomplete data point y\mathbf{y} where only some subset, Ω(y)\Omega(\mathbf{y}) is known. Imputation amounts to sampling from p(x(0)Ω(y))p(\mathbf{x}(0)\mid\Omega(\mathbf{y})), which we can accomplish using an unconditional model (see Section I.2). Colorization is a special case of imputation, except that the known data dimensions are coupled. We can decouple these data dimensions with an orthogonal linear transformation, and perform imputation in the transformed space (details in Section I.3). Fig. 4 (right) shows results for inpainting and colorization achieved with unconditional time-dependent score-based models.

Conclusion

We presented a framework for score-based generative modeling based on SDEs. Our work enables a better understanding of existing approaches, new sampling algorithms, exact likelihood computation, uniquely identifiable encoding, latent code manipulation, and brings new conditional generation abilities to the family of score-based generative models.

While our proposed sampling approaches improve results and enable more efficient sampling, they remain slower at sampling than GANs (Goodfellow et al., 2014) on the same datasets. Identifying ways of combining the stable learning of score-based generative models with the fast sampling of implicit models like GANs remains an important research direction. Additionally, the breadth of samplers one can use when given access to score functions introduces a number of hyper-parameters. Future work would benefit from improved methods to automatically select and tune these hyper-parameters, as well as more extensive investigation on the merits and limitations of various samplers.

We would like to thank Nanxin Chen, Ruiqi Gao, Jonathan Ho, Kevin Murphy, Tim Salimans and Han Zhang for their insightful discussions during the course of this project. This research was partially supported by NSF (#1651565, #1522054, #1733686), ONR (N000141912145), AFOSR (FA95501910024), and TensorFlow Research Cloud. Yang Song was partially supported by the Apple PhD Fellowship in AI/ML.

References

Appendix

We include several appendices with additional details, derivations, and results. Our framework allows general SDEs with matrix-valued diffusion coefficients that depend on the state, for which we provide a detailed discussion in Appendix A. We give a full derivation of VE, VP and sub-VP SDEs in Appendix B, and discuss how to use them from a practitioner’s perspective in Appendix C. We elaborate on the probability flow formulation of our framework in Appendix D, including a derivation of the probability flow ODE (Section D.1), exact likelihood computation (Section D.2), probability flow sampling with a fixed discretization strategy (Section D.3), sampling with black-box ODE solvers (Section D.4), and experimental verification on uniquely identifiable encoding (Section D.5). We give a full description of the reverse diffusion sampler in Appendix E, the DDPM-type ancestral sampler for SMLD models in Appendix F, and Predictor-Corrector samplers in Appendix G. We explain our model architectures and detailed experimental settings in Appendix H, with 1024×10241024\times 1024 CelebA-HQ samples therein. Finally, we detail on the algorithms for controllable generation in Appendix I, and include extended results for class-conditional generation (Section I.1), image inpainting (Section I.2), colorization (Section I.3), and a strategy for solving general inverse problems (Section I.4).

Appendix A The framework for more general SDEs

In the main text, we introduced our framework based on a simplified SDE Eq. 5 where the diffusion coefficient is independent of x(t)\mathbf{x}(t). It turns out that our framework can be extended to hold for more general diffusion coefficients. We can consider SDEs in the following form:

According to (Anderson, 1982), the reverse-time SDE is given by (cf., Eq. 6)

where we define F(x):=(f1(x),f2(x),,fd(x))T\nabla\cdot\mathbf{F}(\mathbf{x}):=(\nabla\cdot\mathbf{f}^{1}(\mathbf{x}),\nabla\cdot\mathbf{f}^{2}(\mathbf{x}),\cdots,\nabla\cdot\mathbf{f}^{d}(\mathbf{x}))^{\mkern-1.5mu\mathsf{T}} for a matrix-valued function F(x):=(f1(x),f2(x),,fd(x))T\mathbf{F}(\mathbf{x}):=(\mathbf{f}^{1}(\mathbf{x}),\mathbf{f}^{2}(\mathbf{x}),\cdots,\mathbf{f}^{d}(\mathbf{x}))^{\mkern-1.5mu\mathsf{T}} throughout the paper.

The probability flow ODE corresponding to Eq. 15 has the following form (cf., Eq. 13, see a detailed derivation in Section D.1):

Finally for conditional generation with the general SDE Eq. 15, we can solve the conditional reverse-time SDE below (cf., Eq. 14, details in Appendix I):

When the drift and diffusion coefficient of an SDE are not affine, it can be difficult to compute the transition kernel p0t(x(t)x(0))p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)) in closed form. This hinders the training of score-based models, because Eq. 7 requires knowing x(t)logp0t(x(t)x(0))\nabla_{\mathbf{x}(t)}\log p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)). To overcome this difficulty, we can replace denoising score matching in Eq. 7 with other efficient variants of score matching that do not require computing x(t)logp0t(x(t)x(0))\nabla_{\mathbf{x}(t)}\log p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)). For example, when using sliced score matching (Song et al., 2019a), our training objective Eq. 7 becomes

Appendix B VE, VP and sub-VP SDEs

Below we provide detailed derivations to show that the noise perturbations of SMLD and DDPM are discretizations of the Variance Exploding (VE) and Variance Preserving (VP) SDEs respectively. We additionally introduce sub-VP SDEs, a modification to VP SDEs that often achieves better performance in both sample quality and likelihoods.

First, when using a total of NN noise scales, each perturbation kernel pσi(xx0)p_{\sigma_{i}}(\mathbf{x}\mid\mathbf{x}_{0}) of SMLD can be derived from the following Markov chain:

where zi1N(0,I)\mathbf{z}_{i-1}\sim\mathcal{N}(\mathbf{0},\mathbf{I}), x0pdata\mathbf{x}_{0}\sim p_{\text{data}}, and we have introduced σ0=0\sigma_{0}=0 to simplify the notation. In the limit of NN\to\infty, the Markov chain {xi}i=1N\{\mathbf{x}_{i}\}_{i=1}^{N} becomes a continuous stochastic process {x(t)}t=01\{\mathbf{x}(t)\}_{t=0}^{1}, {σi}i=1N\{\sigma_{i}\}_{i=1}^{N} becomes a function σ(t)\sigma(t), and zi\mathbf{z}_{i} becomes z(t)\mathbf{z}(t), where we have used a continuous time variable tt\in for indexing, rather than an integer i{1,2,,N}i\in\{1,2,\cdots,N\}. Let x(iN)=xi\mathbf{x}\left(\frac{i}{N}\right)=\mathbf{x}_{i}, σ(iN)=σi\sigma\left(\frac{i}{N}\right)=\sigma_{i}, and z(iN)=zi\mathbf{z}\left(\frac{i}{N}\right)=\mathbf{z}_{i} for i=1,2,,Ni=1,2,\cdots,N. We can rewrite Eq. 20 as follows with Δt=1N\Delta t=\frac{1}{N} and t\in\big{\{}0,\frac{1}{N},\cdots,\frac{N-1}{N}\big{\}}:

where the approximate equality holds when Δt1\Delta t\ll 1. In the limit of Δt0\Delta t\rightarrow 0, this converges to

For the perturbation kernels {pαi(xx0)}i=1N\{p_{\alpha_{i}}(\mathbf{x}\mid\mathbf{x}_{0})\}_{i=1}^{N} used in DDPM, the discrete Markov chain is

where zi1N(0,I)\mathbf{z}_{i-1}\sim\mathcal{N}(\mathbf{0},\mathbf{I}). To obtain the limit of this Markov chain when NN\to\infty, we define an auxiliary set of noise scales {βˉi=Nβi}i=1N\{\bar{\beta}_{i}=N\beta_{i}\}_{i=1}^{N}, and re-write Eq. 22 as below

In the limit of NN\to\infty, {βˉi}i=1N\{\bar{\beta}_{i}\}_{i=1}^{N} becomes a function β(t)\beta(t) indexed by tt\in. Let β(iN)=βˉi\beta\left(\frac{i}{N}\right)=\bar{\beta}_{i}, x(iN)=xi\mathbf{x}(\frac{i}{N})=\mathbf{x}_{i}, z(iN)=zi\mathbf{z}(\frac{i}{N})=\mathbf{z}_{i}. We can rewrite the Markov chain Eq. 23 as the following with Δt=1N\Delta t=\frac{1}{N} and t{0,1,,N1N}t\in\{0,1,\cdots,\frac{N-1}{N}\}:

where the approximate equality holds when Δt1\Delta t\ll 1. Therefore, in the limit of Δt0\Delta t\to 0, Eq. 24 converges to the following VP SDE:

So far, we have demonstrated that the noise perturbations used in SMLD and DDPM correspond to discretizations of VE and VP SDEs respectively. The VE SDE always yields a process with exploding variance when tt\to\infty. In contrast, the VP SDE yields a process with bounded variance. In addition, the process has a constant unit variance for all t[0,)t\in[0,\infty) when p(x(0))p(\mathbf{x}(0)) has a unit variance. Since the VP SDE has affine drift and diffusion coefficients, we can use Eq. (5.51) in Särkkä & Solin (2019) to obtain an ODE that governs the evolution of variance

where ΣVP(t)Cov[x(t)]\mathbf{\Sigma}_{\text{VP}}(t)\coloneqq\operatorname{Cov}[\mathbf{x}(t)] for {x(t)}t=01\{\mathbf{x}(t)\}_{t=0}^{1} obeying a VP SDE. Solving this ODE, we obtain

from which it is clear that the variance ΣVP(t)\mathbf{\Sigma}_{\text{VP}}(t) is always bounded given ΣVP(0)\mathbf{\Sigma}_{\text{VP}}(0). Moreover, ΣVP(t)I\mathbf{\Sigma}_{\text{VP}}(t)\equiv\mathbf{I} if ΣVP(0)=I\mathbf{\Sigma}_{\text{VP}}(0)=\mathbf{I}. Due to this difference, we name Eq. 9 as the Variance Exploding (VE) SDE, and Eq. 11 the Variance Preserving (VP) SDE.

Inspired by the VP SDE, we propose a new SDE called the sub-VP SDE, namely

VE, VP and sub-VP SDEs all have affine drift coefficients. Therefore, their perturbation kernels p0t(x(t)x(0))p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)) are all Gaussian and can be computed with Eqs. (5.50) and (5.51) in Särkkä & Solin (2019):

As a result, all SDEs introduced here can be efficiently trained with the objective in Eq. 7.

Appendix C SDEs in the wild

Below we discuss concrete instantiations of VE and VP SDEs whose discretizations yield SMLD and DDPM models, and the specific sub-VP SDE used in our experiments. In SMLD, the noise scales {σi}i=1N\{\sigma_{i}\}_{i=1}^{N} is typically a geometric sequence where σmin\sigma_{\text{min}} is fixed to 0.010.01 and σmax\sigma_{\text{max}} is chosen according to Technique 1 in Song & Ermon (2020). Usually, SMLD models normalize image inputs to the range $.Since. Since\{\sigma_{i}\}_{i=1}^{N}isageometricsequence,wehaveis a geometric sequence, we have\sigma(\frac{i}{N})=\sigma_{i}=\sigma_{\text{min}}\left(\frac{\sigma_{\text{max}}}{\sigma_{\text{min}}}\right)^{\frac{i-1}{N-1}}forfori=1,2,\cdots,N.Inthelimitof. In the limit ofN\to\infty,wehave, we have\sigma(t)=\sigma_{\text{min}}\left(\frac{\sigma_{\text{max}}}{\sigma_{\text{min}}}\right)^{t}forfort\in(0,1]$. The corresponding VE SDE is

and the perturbation kernel can be derived via Eq. 29:

There is one subtlety when t=0t=0: by definition, σ(0)=σ0=0\sigma(0)=\sigma_{0}=0 (following the convention in Eq. 20), but σ(0+)limt0+σ(t)=σmin0\sigma(0^{+})\coloneqq\lim_{t\to 0^{+}}\sigma(t)=\sigma_{\text{min}}\neq 0. In other words, σ(t)\sigma(t) for SMLD is not differentiable since σ(0)σ(0+)\sigma(0)\neq\sigma(0^{+}), causing the VE SDE in Eq. 21 undefined for t=0t=0. In practice, we bypass this issue by always solving the SDE and its associated probability flow ODE in the range t[ϵ,1]t\in[\epsilon,1] for some small constant ϵ>0\epsilon>0, and we use ϵ=105\epsilon=10^{-5} in our VE SDE experiments.

For DDPM models, {βi}i=1N\{\beta_{i}\}_{i=1}^{N} is typically an arithmetic sequence where βi=βˉminN+i1N(N1)(βˉmaxβˉmin)\beta_{i}=\frac{\bar{\beta}_{\text{min}}}{N}+\frac{i-1}{N(N-1)}(\bar{\beta}_{\text{max}}-\bar{\beta}_{\text{min}}) for i=1,2,,Ni=1,2,\cdots,N. Therefore, β(t)=βˉmin+t(βˉmaxβˉmin)\beta(t)=\bar{\beta}_{\text{min}}+t(\bar{\beta}_{\text{max}}-\bar{\beta}_{\text{min}}) for tt\in in the limit of NN\to\infty. This corresponds to the following instantiation of the VP SDE:

where x(0)pdata(x)\mathbf{x}(0)\sim p_{\text{data}}(\mathbf{x}). In our experiments, we let βˉmin=0.1\bar{\beta}_{\text{min}}=0.1 and βˉmax=20\bar{\beta}_{\text{max}}=20 to match the settings in Ho et al. (2020). The perturbation kernel is given by

For DDPM, there is no discontinuity issue with the corresponding VP SDE; yet, there are numerical instability issues for training and sampling at t=0t=0, due to the vanishing variance of x(t)\mathbf{x}(t) as t0t\to 0. Therefore, same as the VE SDE, we restrict computation to t[ϵ,1]t\in[\epsilon,1] for a small ϵ>0\epsilon>0. For sampling, we choose ϵ=103\epsilon=10^{-3} so that the variance of x(ϵ)\mathbf{x}(\epsilon) in VP SDE matches the variance of x1\mathbf{x}_{1} in DDPM; for training and likelihood computation, we adopt ϵ=105\epsilon=10^{-5} which empirically gives better results.

As a sanity check for our SDE generalizations to SMLD and DDPM, we compare the perturbation kernels of SDEs and original discrete Markov chains in Fig. 5. The SMLD and DDPM models both use N=1000N=1000 noise scales. For SMLD, we only need to compare the variances of perturbation kernels since means are the same by definition. For DDPM, we compare the scaling factors of means and the variances. As demonstrated in Fig. 5, the discrete perturbation kernels of original SMLD and DDPM models align well with perturbation kernels derived from VE and VP SDEs.

For sub-VP SDEs, we use exactly the same β(t)\beta(t) as VP SDEs. This leads to the following perturbation kernel

We also restrict numerical computation to the same interval of [ϵ,1][\epsilon,1] as VP SDEs.

Empirically, we observe that smaller ϵ\epsilon generally yields better likelihood values for all SDEs. For sampling, it is important to use an appropriate ϵ\epsilon for better Inception scores and FIDs, although samples across different ϵ\epsilon look visually the same to human eyes.

Appendix D Probability flow ODE

The idea of probability flow ODE is inspired by Maoutsa et al. (2020), and one can find the derivation of a simplified case therein. Below we provide a derivation for the fully general ODE in Eq. 17. We consider the SDE in Eq. 15, which possesses the following form:

based on which we can continue the rewriting of Eq. 36 to obtain

same as the probability flow ODE given by Eq. 17. Therefore, we have shown that the probability flow ODE Eq. 17 induces the same marginal probability density pt(x)p_{t}(\mathbf{x}) as the SDE in Eq. 15.

D.2 Likelihood computation

The probability flow ODE in Eq. 17 has the following form when we replace the score xlogpt(x)\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}) with the time-dependent score-based model sθ(x,t)\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{x},t):

With the instantaneous change of variables formula (Chen et al., 2018), we can compute the log-likelihood of p0(x)p_{0}(\mathbf{x}) using

In our experiments, we use the RK45 ODE solver (Dormand & Prince, 1980) provided by scipy.integrate.solve_ivp in all cases. The bits/dim values in Table 2 are computed with atol=1e-5 and rtol=1e-5, same as Grathwohl et al. (2018). To give the likelihood results of our models in Table 2, we average the bits/dim obtained on the test dataset over five different runs with ϵ=105\epsilon=10^{-5} (see definition of ϵ\epsilon in Appendix C).

D.3 Probability flow sampling

where ziN(0,I)\mathbf{z}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{I}). We assume the discretization schedule of time is fixed beforehand, and thus we absorb the dependency on Δt\Delta t into the notations of fi\mathbf{f}_{i} and Gi\mathbf{G}_{i}. Using Eq. 17, we can obtain the following probability flow ODE:

We may employ any numerical method to integrate the probability flow ODE backwards in time for sample generation. In particular, we propose a discretization in a similar functional form to Eq. 41:

where the score-based model sθ(xi,i)\mathbf{s}_{{\boldsymbol{\theta}}^{*}}(\mathbf{x}_{i},i) is conditioned on the iteration number ii. This is a deterministic iteration rule. Unlike reverse diffusion samplers or ancestral sampling, there is no additional randomness once the initial sample xN\mathbf{x}_{N} is obtained from the prior distribution. When applied to SMLD models, we can get the following iteration rule for probability flow sampling:

D.4 Sampling with black-box ODE solvers

For producing figures in Fig. 3, we use a DDPM model trained on 256×256256\times 256 CelebA-HQ with the same settings in Ho et al. (2020). All FID scores of our models in Table 2 are computed on samples from the RK45 ODE solver implemented in scipy.integrate.solve_ivp with atol=1e-5 and rtol=1e-5. We use ϵ=105\epsilon=10^{-5} for VE SDEs and ϵ=103\epsilon=10^{-3} for VP SDEs (see also Appendix C).

Aside from the interpolation results in Fig. 3, we demonstrate more examples of latent space manipulation in Fig. 6, including interpolation and temperature scaling. The model tested here is a DDPM model trained with the same settings in Ho et al. (2020).

Although solvers for the probability flow ODE allow fast sampling, their samples typically have higher (worse) FID scores than those from SDE solvers if no corrector is used. We have this empirical observation for both the discretization strategy in Section D.3, and black-box ODE solvers introduced above. Moreover, the performance of probability flow ODE samplers depends on the choice of the SDE—their sample quality for VE SDEs is much worse than VP SDEs especially for high-dimensional data.

D.5 Uniquely identifiable encoding

As a sanity check, we train two models (denoted as “Model A” and “Model B”) with different architectures using the VE SDE on CIFAR-10. Here Model A is an NCSN++ model with 4 layers per resolution trained using the continuous objective in Eq. 7, and Model B is all the same except that it uses 8 layers per resolution. Model definitions are in Appendix H.

We report the latent codes obtained by Model A and Model B for a random CIFAR-10 image in Fig. 7. In Fig. 8, we show the dimension-wise differences and correlation coefficients between latent encodings on a total of 16 CIFAR-10 images. Our results demonstrate that for the same inputs, Model A and Model B provide encodings that are close in every dimension, despite having different model architectures and training runs.

Appendix E Reverse diffusion sampling

and suppose the following iteration rule is a discretization of it:

where ziN(0,I)\mathbf{z}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{I}). Here we assume the discretization schedule of time is fixed beforehand, and thus we can absorb it into the notations of fi\mathbf{f}_{i} and Gi\mathbf{G}_{i}.

Based on Eq. 45, we propose to discretize the reverse-time SDE

with a similar functional form, which gives the following iteration rule for i{0,1,,N1}i\in\{0,1,\cdots,N-1\}:

where our trained score-based model sθ(xi,i)\mathbf{s}_{{\boldsymbol{\theta}}^{*}}(\mathbf{x}_{i},i) is conditioned on iteration number ii.

When applying Eq. 46 to Eqs. 20 and 10, we obtain a new set of numerical solvers for the reverse-time VE and VP SDEs, resulting in sampling algorithms as shown in the “predictor” part of Algorithms 2 and 3. We name these sampling methods (that are based on the discretization strategy in Eq. 46) reverse diffusion samplers.

As expected, the ancestral sampling of DDPM (Ho et al., 2020) (Eq. 4) matches its reverse diffusion counterpart when βi0\beta_{i}\to 0 for all ii (which happens when Δt0\Delta t\to 0 since βi=βˉiΔt\beta_{i}=\bar{\beta}_{i}\Delta t, see Appendix B), because

Therefore, the original ancestral sampler of Eq. 4 is essentially a different discretization to the same reverse-time SDE. This unifies the sampling method in Ho et al. (2020) as a numerical solver to the reverse-time VP SDE in our continuous framework.

Appendix F Ancestral sampling for SMLD models

The ancestral sampling method for DDPM models can also be adapted to SMLD models. Consider a sequence of noise scales σ1<σ2<<σN\sigma_{1}<\sigma_{2}<\cdots<\sigma_{N} as in SMLD. By perturbing a data point x0\mathbf{x}_{0} with these noise scales sequentially, we obtain a Markov chain x0x1xN\mathbf{x}_{0}\to\mathbf{x}_{1}\to\cdots\to\mathbf{x}_{N}, where

Here we assume σ0=0\sigma_{0}=0 to simplify notations. Following Ho et al. (2020), we can compute

If we parameterize the reverse transition kernel as pθ(xi1xi)=N(xi1;μθ(xi,i),τi2I)p_{\boldsymbol{\theta}}(\mathbf{x}_{i-1}\mid\mathbf{x}_{i})=\mathcal{N}(\mathbf{x}_{i-1};{\boldsymbol{\mu}}_{\boldsymbol{\theta}}(\mathbf{x}_{i},i),\tau^{2}_{i}\mathbf{I}), then

where Lt1L_{t-1} is one representative term in the ELBO objective (see Eq. (8) in Ho et al. (2020)), CC is a constant that does not depend on θ{\boldsymbol{\theta}}, zN(0,I)\mathbf{z}\sim\mathcal{N}(\mathbf{0},\mathbf{I}), and xi(x0,z)=x0+σiz\mathbf{x}_{i}(\mathbf{x}_{0},\mathbf{z})=\mathbf{x}_{0}+\sigma_{i}\mathbf{z}. We can therefore parameterize μθ(xi,i){\boldsymbol{\mu}}_{\boldsymbol{\theta}}(\mathbf{x}_{i},i) via

where sθ(xi,i)\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{x}_{i},i) is to estimate z/σi\mathbf{z}/\sigma_{i}. As in Ho et al. (2020), we let τi=σi12(σi2σi12)σi2\tau_{i}=\sqrt{\frac{\sigma_{i-1}^{2}(\sigma_{i}^{2}-\sigma_{i-1}^{2})}{\sigma_{i}^{2}}}. Through ancestral sampling on i=1Npθ(xi1xi)\prod_{i=1}^{N}p_{\boldsymbol{\theta}}(\mathbf{x}_{i-1}\mid\mathbf{x}_{i}), we obtain the following iteration rule

where xNN(0,σN2I)\mathbf{x}_{N}\sim\mathcal{N}(\mathbf{0},\sigma_{N}^{2}\mathbf{I}), θ{\boldsymbol{\theta}}^{*} denotes the optimal parameter of sθ\mathbf{s}_{\boldsymbol{\theta}}, and ziN(0,I)\mathbf{z}_{i}\sim\mathcal{N}(\mathbf{0},\mathbf{I}). We call Eq. 47 the ancestral sampling method for SMLD models.

Appendix G Predictor-Corrector samplers

The predictor can be any numerical solver for the reverse-time SDE with a fixed discretization strategy. The corrector can be any score-based MCMC approach. In PC sampling, we alternate between the predictor and corrector, as described in Algorithm 1. For example, when using the reverse diffusion SDE solver (Appendix E) as the predictor, and annealed Langevin dynamics (Song & Ermon, 2019) as the corrector, we have Algorithms 2 and 3 for VE and VP SDEs respectively, where {ϵi}i=0N1\{\epsilon_{i}\}_{i=0}^{N-1} are step sizes for Langevin dynamics as specified below.

We take the schedule of annealed Langevin dynamics in Song & Ermon (2019), but re-frame it with slight modifications in order to get better interpretability and empirical performance. We provide the corrector algorithms in Algorithms 4 and 5 respectively, where we call rr the “signal-to-noise” ratio. We determine the step size ϵ\epsilon using the norm of the Gaussian noise z2\left\lVert\mathbf{z}\right\rVert_{2}, norm of the score-based model sθ2\left\lVert\mathbf{s}_{{\boldsymbol{\theta}}^{*}}\right\rVert_{2} and the signal-to-noise ratio rr. When sampling a large batch of samples together, we replace the norm 2\left\lVert\cdot\right\rVert_{2} with the average norm across the mini-batch. When the batch size is small, we suggest replacing z2\left\lVert\mathbf{z}\right\rVert_{2} with d\sqrt{d}, where dd is the dimensionality of z\mathbf{z}.

For both SMLD and DDPM models, the generated samples typically contain small noise that is hard to detect by humans. As noted by Jolicoeur-Martineau et al. (2020), FIDs can be significantly worse without removing this noise. This unfortunate sensitivity to noise is also part of the reason why NCSN models trained with SMLD has been performing worse than DDPM models in terms of FID, because the former does not use a denoising step at the end of sampling, while the latter does. In all experiments of this paper we ensure there is a single denoising step at the end of sampling, using Tweedie’s formula (Efron, 2011).

We use the same architecture in Ho et al. (2020) for our score-based models. For the VE SDE, we train a model with the original SMLD objective in Eq. 1; similarly for the VP SDE, we use the original DDPM objective in Eq. 3. We apply a total number of 1000 noise scales for training both models. For results in Fig. 9, we train an NCSN++ model (definition in Appendix H) on 256×256256\times 256 LSUN bedroom and church_outdoor (Yu et al., 2015) datasets with the VE SDE and our continuous objective Eq. 7. The batch size is fixed to 128 on CIFAR-10 and 64 on LSUN.

Models in this experiment are all trained with 1000 noise scales. To get results for P2000 (predictor-only sampler using 2000 steps) which requires 2000 noise scales, we need to interpolate between 1000 noise scales at test time. The specific architecture of the noise-conditional score-based model in Ho et al. (2020) uses sinusoidal positional embeddings for conditioning on integer time steps. This allows us to interpolate between noise scales at test time in an ad-hoc way (while it is hard to do so for other architectures like the one in Song & Ermon (2019)). Specifically, for SMLD models, we keep σmin\sigma_{\text{min}} and σmax\sigma_{\text{max}} fixed and double the number of time steps. For DDPM models, we halve βmin\beta_{\text{min}} and βmax\beta_{\text{max}} before doubling the number of time steps. Suppose {sθ(x,i)}i=0N1\{\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{x},i)\}_{i=0}^{N-1} is a score-based model trained on NN time steps, and let {sθ(x,i)}i=02N1\{\mathbf{s}_{\boldsymbol{\theta}}^{\prime}(\mathbf{x},i)\}_{i=0}^{2N-1} denote the corresponding interpolated score-based model at 2N2N time steps. We test two different interpolation strategies for time steps: linear interpolation where sθ(x,i)=sθ(x,i/2)\mathbf{s}_{\boldsymbol{\theta}}^{\prime}(\mathbf{x},i)=\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{x},i/2) and rounding interpolation where sθ(x,i)=sθ(x,i/2)\mathbf{s}_{\boldsymbol{\theta}}^{\prime}(\mathbf{x},i)=\mathbf{s}_{\boldsymbol{\theta}}(\mathbf{x},\lfloor i/2\rfloor). We provide results with linear interpolation in Table 1, and give results of rounding interpolation in Table 4. We observe that different interpolation methods result in performance differences but maintain the general trend of predictor-corrector methods performing on par or better than predictor-only or corrector-only samplers.

For Predictor-Corrector and corrector-only samplers on CIFAR-10, we search for the best signal-to-noise ratio (rr) over a grid that increments at 0.01. We report the best rr in Table 5. For LSUN bedroom/church_outdoor, we fix rr to 0.075. Unless otherwise noted, we use one corrector step per noise scale for all PC samplers. We use two corrector steps per noise scale for corrector-only samplers on CIFAR-10. For sample generation, the batch size is 1024 on CIFAR-10 and 8 on LSUN bedroom/church_outdoor.

Appendix H Architecture improvements

We explored several architecture designs to improve score-based models for both VE and VP SDEs. Our endeavor gives rise to new state-of-the-art sample quality on CIFAR-10, new state-of-the-art likelihood on uniformly dequantized CIFAR-10, and enables the first high-fidelity image samples of resolution 1024×10241024\times 1024 from score-based generative models. Code and checkpoints are open-sourced at https://github.com/yang-song/score_sde.

Unless otherwise noted, all models are trained for 1.3M iterations, and we save one checkpoint per 50k iterations. For VE SDEs, we consider two datasets: 32×3232\times 32 CIFAR-10 (Krizhevsky et al., 2009) and 64×6464\times 64 CelebA (Liu et al., 2015), pre-processed following Song & Ermon (2020). We compare different configurations based on their FID scores averaged over checkpoints after 0.5M iterations. For VP SDEs, we only consider the CIFAR-10 dataset to save computation, and compare models based on the average FID scores over checkpoints obtained between 0.25M and 0.5M iterations, because FIDs turn to increase after 0.5M iterations for VP SDEs.

All FIDs are computed on 50k samples with tensorflow_gan. For sampling, we use the PC sampler discretized at 1000 time steps. We choose reverse diffusion (see Appendix E) as the predictor. We use one corrector step per update of the predictor for VE SDEs with a signal-to-noise ratio of 0.16, but save the corrector step for VP SDEs since correctors there only give slightly better results but require double computation. We follow Ho et al. (2020) for optimization, including the learning rate, gradient clipping, and learning rate warm-up schedules. Unless otherwise noted, models are trained with the original discrete SMLD and DDPM objectives in Eqs. 1 and 3 and use a batch size of 128. The optimal architectures found under these settings are subsequently transferred to continuous objectives and deeper models. We also directly transfer the best architecture for VP SDEs to sub-VP SDEs, given the similarity of these two SDEs.

Our architecture is mostly based on Ho et al. (2020). We additionally introduce the following components to maximize the potential improvement of score-based models.

Upsampling and downsampling images with anti-aliasing based on Finite Impulse Response (FIR) (Zhang, 2019). We follow the same implementation and hyper-parameters in StyleGAN-2 (Karras et al., 2020b).

Rescaling all skip connections by \nicefrac12\nicefrac{{1}}{{\sqrt{2}}}. This has been demonstrated effective in several best-in-class GAN models, including ProgressiveGAN (Karras et al., 2018), StyleGAN (Karras et al., 2019) and StyleGAN-2 (Karras et al., 2020b).

Replacing the original residual blocks in DDPM with residual blocks from BigGAN (Brock et al., 2018).

Increasing the number of residual blocks per resolution from 22 to 44.

Incorporating progressive growing architectures. We consider two progressive architectures for input: “input skip” and “residual”, and two progressive architectures for output: “output skip” and “residual”. These progressive architectures are defined and implemented according to StyleGAN-2.

We also tested equalized learning rates, a trick used in very successful models like ProgressiveGAN (Karras et al., 2018) and StyleGAN (Karras et al., 2019). However, we found it harmful at an early stage of our experiments, and therefore decided not to explore more on it.

The exponential moving average (EMA) rate has a significant impact on performance. For models trained with VE perturbations, we notice that 0.999 works better than 0.9999, whereas for models trained with VP perturbations it is the opposite. We therefore use an EMA rate of 0.999 and 0.9999 for VE and VP models respectively.

H.2 Results on CIFAR-10

All architecture components introduced above can improve the performance of score-based models trained with VE SDEs, as shown in Fig. 10. The box plots demonstrate the importance of each component when other components can vary freely. On both CIFAR-10 and CelebA, the additional components that we explored always improve the performance on average for VE SDEs. For progressive growing, it is not clear which combination of configurations consistently performs the best, but the results are typically better than when no progressive growing architecture is used. Our best score-based model for VE SDEs 1) uses FIR upsampling/downsampling, 2) rescales skip connections, 3) employs BigGAN-type residual blocks, 4) uses 4 residual blocks per resolution instead of 2, and 5) uses “residual” for input and no progressive growing architecture for output. We name this model “NCSN++”, following the naming convention of previous SMLD models (Song & Ermon, 2019; 2020).

We followed a similar procedure to examine these architecture components for VP SDEs, except that we skipped experiments on CelebA due to limited computing resources. The NCSN++ architecture worked decently well for VP SDEs, ranked 4th place over all 144 possible configurations. The top configuration, however, has a slightly different structure, which uses no FIR upsampling/downsampling and no progressive growing architecture compared to NCSN++. We name this model “DDPM++”, following the naming convention of Ho et al. (2020).

The basic NCSN++ model with 4 residual blocks per resolution achieves an FID of 2.45 on CIFAR-10, whereas the basic DDPM++ model achieves an FID of 2.78. Here in order to match the convention used in Karras et al. (2018); Song & Ermon (2019) and Ho et al. (2020), we report the lowest FID value over the course of training, rather than the average FID value over checkpoints after 0.5M iterations (used for comparing different models of VE SDEs) or between 0.25M and 0.5M iterations (used for comparing VP SDE models) in our architecture exploration.

Switching from discrete training objectives to continuous ones in Eq. 7 further improves the FID values for all SDEs. To condition the NCSN++ model on continuous time variables, we change positional embeddings, the layers in Ho et al. (2020) for conditioning on discrete time steps, to random Fourier feature embeddings (Tancik et al., 2020). The scale parameter of these random Fourier feature embeddings is fixed to 16. We also reduce the number of training iterations to 0.95M to suppress overfitting. These changes improve the FID on CIFAR-10 from 2.45 to 2.38 for NCSN++ trained with the VE SDE, resulting in a model called “NCSN++ cont.”. In addition, we can further improve the FID from 2.38 to 2.20 by doubling the number of residual blocks per resolution for NCSN++ cont., resulting in the model denoted as “NCSN++ cont. (deep)”. All quantitative results are summarized in Table 3, and we provide random samples from our best model in Fig. 11.

Similarly, we can also condition the DDPM++ model on continuous time steps, resulting in a model “DDPM++ cont.”. When trained with the VP SDE, it improves the FID of 2.78 from DDPM++ to 2.55. When trained with the sub-VP SDE, it achieves an FID of 2.61. To get better performance, we used the Euler-Maruyama solver as the predictor for continuously-trained models, instead of the ancestral sampling predictor or the reverse diffusion predictor. This is because the discretization strategy of the original DDPM method does not match the variance of the continuous process well when t0t\to 0, which significantly hurts FID scores. As shown in Table 2, the likelihood values are 3.21 and 3.05 bits/dim for VP and sub-VP SDEs respectively. Doubling the depth, and trainin with 0.95M iterations, we can improve both FID and bits/dim for both VP and sub-VP SDEs, leading to a model “DDPM++ cont. (deep)”. Its FID score is 2.41, same for both VP and sub-VP SDEs. When trained with the sub-VP SDE, it can achieve a likelihood of 2.99 bits/dim. Here all likelihood values are reported for the last checkpoint during training.

H.3 High resolution images

Encouraged by the success of NCSN++ on CIFAR-10, we proceed to test it on 1024×10241024\times 1024 CelebA-HQ (Karras et al., 2018), a task that was previously only achievable by some GAN models and VQ-VAE-2 (Razavi et al., 2019). We used a batch size of 8, increased the EMA rate to 0.9999, and trained a model similar to NCSN++ with the continuous objective (Eq. 7) for around 2.4M iterations (please find the detailed architecture in our code release.) We use the PC sampler discretized at 2000 steps with the reverse diffusion predictor, one Langevin step per predictor update and a signal-to-noise ratio of 0.15. The scale parameter for the random Fourier feature embeddings is fixed to 16. We use the “input skip” progressive architecture for the input, and “output skip” progressive architecture for the output. We provide samples in Fig. 12. Although these samples are not perfect (e.g., there are visible flaws on facial symmetry), we believe these results are encouraging and can demonstrate the scalability of our approach. Future work on more effective architectures are likely to significantly advance the performance of score-based generative models on this task.

Appendix I Controllable generation

Consider a forward SDE with the following general form

and suppose the initial state distribution is p0(x(0)y)p_{0}(\mathbf{x}(0)\mid\mathbf{y}). The density at time tt is pt(x(t)y)p_{t}(\mathbf{x}(t)\mid\mathbf{y}) when conditioned on y\mathbf{y}. Therefore, using Anderson (1982), the reverse-time SDE is given by

Since pt(x(t)y)pt(x(t))p(yx(t))p_{t}(\mathbf{x}(t)\mid\mathbf{y})\propto p_{t}(\mathbf{x}(t))p(\mathbf{y}\mid\mathbf{x}(t)), the score xlogpt(x(t)y)\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}(t)\mid\mathbf{y}) can be computed easily by

This subsumes the conditional reverse-time SDE in Eq. 14 as a special case. All sampling methods we have discussed so far can be applied to the conditional reverse-time SDE for sample generation.

When y\mathbf{y} represents class labels, we can train a time-dependent classifier pt(yx(t))p_{t}(\mathbf{y}\mid\mathbf{x}(t)) for class-conditional sampling. Since the forward SDE is tractable, we can easily create a pair of training data (x(t),y)(\mathbf{x}(t),\mathbf{y}) by first sampling (x(0),y)(\mathbf{x}(0),\mathbf{y}) from a dataset and then obtaining x(t)p0t(x(t)x(0))\mathbf{x}(t)\sim p_{0t}(\mathbf{x}(t)\mid\mathbf{x}(0)). Afterwards, we may employ a mixture of cross-entropy losses over different time steps, like Eq. 7, to train the time-dependent classifier pt(yx(t))p_{t}(\mathbf{y}\mid\mathbf{x}(t)).

To test this idea, we trained a Wide ResNet (Zagoruyko & Komodakis, 2016) (Wide-ResNet-28-10) on CIFAR-10 with VE perturbations. The classifier is conditioned on logσi\log\sigma_{i} using random Fourier features (Tancik et al., 2020), and the training objective is a simple sum of cross-entropy losses sampled at different scales. We provide a plot to show the accuracy of this classifier over noise scales in Fig. 13. The score-based model is an unconditional NCSN++ (4 blocks/resolution) in Table 3, and we generate samples using the PC algorithm with 2000 discretization steps. The class-conditional samples are provided in Fig. 4, and an extended set of conditional samples is given in Fig. 13.

I.2 Imputation

Imputation is a special case of conditional sampling. Denote by Ω(x)\Omega(\mathbf{x}) and Ωˉ(x){\bar{\Omega}}(\mathbf{x}) the known and unknown dimensions of x\mathbf{x} respectively, and let fΩˉ(,t)\mathbf{f}_{{\bar{\Omega}}}(\cdot,t) and GΩˉ(,t)\mathbf{G}_{{\bar{\Omega}}}(\cdot,t) denote f(,t)\mathbf{f}(\cdot,t) and G(,t)\mathbf{G}(\cdot,t) restricted to the unknown dimensions. For VE/VP SDEs, the drift coefficient f(,t)\mathbf{f}(\cdot,t) is element-wise, and the diffusion coefficient G(,t)\mathbf{G}(\cdot,t) is diagonal. When f(,t)\mathbf{f}(\cdot,t) is element-wise, fΩˉ(,t)\mathbf{f}_{{\bar{\Omega}}}(\cdot,t) denotes the same element-wise function applied only to the unknown dimensions. When G(,t)\mathbf{G}(\cdot,t) is diagonal, GΩˉ(,t)\mathbf{G}_{{\bar{\Omega}}}(\cdot,t) denotes the sub-matrix restricted to unknown dimensions.

For imputation, our goal is to sample from p(Ωˉ(x(0))Ω(x(0))=y)p({\bar{\Omega}}(\mathbf{x}(0))\mid\Omega(\mathbf{x}(0))=\mathbf{y}). Define a new diffusion process z(t)=Ωˉ(x(t))\mathbf{z}(t)={\bar{\Omega}}(\mathbf{x}(t)), and note that the SDE for z(t)\mathbf{z}(t) can be written as

The reverse-time SDE, conditioned on Ω(x(0))=y\Omega(\mathbf{x}(0))=\mathbf{y}, is given by

Although pt(z(t)Ω(x(0))=y)p_{t}(\mathbf{z}(t)\mid\Omega(\mathbf{x}(0))=\mathbf{y}) is in general intractable, it can be approximated. Let AA denote the event Ω(x(0))=y\Omega(\mathbf{x}(0))=\mathbf{y}. We have

where Ω^(x(t))\hat{\Omega}(\mathbf{x}(t)) is a random sample from pt(Ω(x(t))A)p_{t}(\Omega(\mathbf{x}(t))\mid A), which is typically a tractable distribution. Therefore,

where [z(t);Ω^(x(t))][\mathbf{z}(t);\hat{\Omega}(\mathbf{x}(t))] denotes a vector u(t)\mathbf{u}(t) such that Ω(u(t))=Ω^(x(t))\Omega(\mathbf{u}(t))=\hat{\Omega}(\mathbf{x}(t)) and Ωˉ(u(t))=z(t)\bar{\Omega}(\mathbf{u}(t))=\mathbf{z}(t), and the identity holds because zlogpt([z(t);Ω^(x(t))])=zlogpt(z(t)Ω^(x(t)))+zlogp(Ω^(x(t)))=zlogpt(z(t)Ω^(x(t)))\nabla_{\mathbf{z}}\log p_{t}([\mathbf{z}(t);\hat{\Omega}(\mathbf{x}(t))])=\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z}(t)\mid\hat{\Omega}(\mathbf{x}(t)))+\nabla_{\mathbf{z}}\log p(\hat{\Omega}(\mathbf{x}(t)))=\nabla_{\mathbf{z}}\log p_{t}(\mathbf{z}(t)\mid\hat{\Omega}(\mathbf{x}(t))).

We provided an extended set of inpainting results in Figs. 14 and 15.

I.3 Colorization

Colorization is a special case of imputation, except that the known data dimensions are coupled. We can decouple these data dimensions by using an orthogonal linear transformation to map the gray-scale image to a separate channel in a different space, and then perform imputation to complete the other channels before transforming everything back to the original image space. The orthogonal matrix we used to decouple color channels is

Because the transformations are all orthogonal matrices, the standard Wiener process w(t)\mathbf{w}(t) will still be a standard Wiener process in the transformed space, allowing us to build an SDE and use the same imputation method in Section I.2. We provide an extended set of colorization results in Figs. 16 and 17.

I.4 Solving general inverse problems

Suppose we have two random variables x\mathbf{x} and y\mathbf{y}, and we know the forward process of generating y\mathbf{y} from x\mathbf{x}, given by p(yx)p(\mathbf{y}\mid\mathbf{x}). The inverse problem is to obtain x\mathbf{x} from y\mathbf{y}, that is, generating samples from p(xy)p(\mathbf{x}\mid\mathbf{y}). In principle, we can estimate the prior distribution p(x)p(\mathbf{x}) and obtain p(xy)p(\mathbf{x}\mid\mathbf{y}) using Bayes’ rule: p(xy)=p(x)p(yx)/p(y)p(\mathbf{x}\mid\mathbf{y})=p(\mathbf{x})p(\mathbf{y}\mid\mathbf{x})/p(\mathbf{y}). In practice, however, both estimating the prior and performing Bayesian inference are non-trivial.

Leveraging Eq. 48, score-based generative models provide one way to solve the inverse problem. Suppose we have a diffusion process {x(t)}t=0T\{\mathbf{x}(t)\}_{t=0}^{T} generated by perturbing x\mathbf{x} with an SDE, and a time-dependent score-based model sθ(x(t),t)\mathbf{s}_{{\boldsymbol{\theta}}*}(\mathbf{x}(t),t) trained to approximate xlogpt(x(t))\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}(t)). Once we have an estimate of xlogpt(x(t)y)\nabla_{\mathbf{x}}\log p_{t}(\mathbf{x}(t)\mid\mathbf{y}), we can simulate the reverse-time SDE in Eq. 48 to sample from p0(x(0)y)=p(xy)p_{0}(\mathbf{x}(0)\mid\mathbf{y})=p(\mathbf{x}\mid\mathbf{y}). To obtain this estimate, we first observe that

where y(t)\mathbf{y}(t) is defined via x(t)\mathbf{x}(t) and the forward process p(y(t)x(t))p(\mathbf{y}(t)\mid\mathbf{x}(t)). Now assume two conditions:

p(y(t)y)p(\mathbf{y}(t)\mid\mathbf{y}) is tractable. We can often derive this distribution from the interaction between the forward process and the SDE, like in the case of image imputation and colorization.

pt(x(t)y(t),y)pt(x(t)y(t))p_{t}(\mathbf{x}(t)\mid\mathbf{y}(t),\mathbf{y})\approx p_{t}(\mathbf{x}(t)\mid\mathbf{y}(t)). For small tt, y(t)\mathbf{y}(t) is almost the same as y\mathbf{y} so the approximation holds. For large tt, y\mathbf{y} becomes further away from x(t)\mathbf{x}(t) in the Markov chain, and thus have smaller impact on x(t)\mathbf{x}(t). Moreover, the approximation error for large tt matter less for the final sample, since it is used early in the sampling process.

where y^(t)\hat{\mathbf{y}}(t) is a sample from p(y(t)y)p(\mathbf{y}(t)\mid\mathbf{y}). Now we can plug Eq. 50 into Eq. 48 and solve the resulting reverse-time SDE to generate samples from p(xy)p(\mathbf{x}\mid\mathbf{y}).