Convergence for score-based generative modeling with polynomial complexity

Introduction

A key task in machine learning is to learn a probability distribution from data, in a way that allows efficient generation of additional samples from the learned distribution. Score-based generative modeling (SGM) is one empirically successful approach that implicitly learns the probability distribution by learning how to transform white noise into the data distribution, and gives state-of-the-art performance for generating images and audio [SE19, Dat+19, Gra+19, SE20, Son+20a, Men+21, Son+21a, Son+21, Jin+22]. It also yields a conditional generation process for inverse problems [DN21]. The basic idea behind score-based generative modeling is to first estimate the score function from data [Son+20] and then to sample the distribution based on the learned score function. Other approaches for generative modeling include generative adversarial networks (GANs) [Goo+14, ACB17], normalizing flows [DSB16], variational autoencoders [KW19], and energy-based models [ZML16]. While score-based generative modeling has achieved great success, its theoretical analysis is still lacking and is the focus of our work.

The score function of a distribution $P$ with density $p$ is defined as the gradient of the log-pdf, $\nabla\ln p$. Its significance arises from the fact that knowing the score function allows running a variety of sampling algorithms, based on discretizations of stochastic differential equations (SDE’s), to sample from $p$. SGM consists of two steps: first, learning an estimate of the score function for a sequence of “noisy” versions of the data distribution $P_{\textup{data}}$, and second, using the score function in lieu of the gradient of the log-pdf in the chosen sampling algorithm. We now describe each of these steps more precisely.

First, a method of adding noise to the data distribution is fixed; this takes the form of evolving a (forward) stochastic differential equation (SDE) starting from the data distribution. We fix a sequence of noise levels $\sigma_{1}<\cdots<\sigma_{N}$. For $\sigma\in\{\sigma_{1},\ldots,\sigma_{N}\}$, let the resulting distributions be $P_{\sigma^{2}}$ and the distributions conditional on the starting data point be $P_{\sigma^{2}}(\cdot|x)$. Typically, $\sigma_{1}$ is chosen so that $P_{\sigma_{1}^{2}}\approx P_{\textup{data}}$ and $P_{\sigma_{N}^{2}}$ is close to some “prior” distribution that is easy to sample from, such as $N(0,\sigma_{N}^{2}I_{d})$. While the score $\nabla\ln p_{\sigma^{2}}$ cannot be estimated directly, it turns out that a de-noising objective that is equivalent to the score-matching objective can be calculated [SE19]. This de-noising objective can be estimated from samples $(X,\widetilde{X})$ where $\widetilde{X}\sim P_{\sigma^{2}}(\cdot|x)$. The objective is represented and optimized within an expressive function class, typically neural networks, to obtain a $L^{2}$-estimate of the score, that is, $s_{\theta}(x,\sigma^{2})$ such that

The reason we estimate the score function $\nabla\ln p_{\sigma^{2}}$ is that there are a variety of sampling algorithms—based on simulating SDE’s—that can sample from $p$ given access to $\nabla\ln p$, including Langevin Monte Carlo and Hamiltonian Monte Carlo. The second step is then to use the estimated score function $s_{\theta}(x,t)$ in lieu of the exact gradient in the sampling algorithm to successively obtain samples from $p_{\sigma_{N}^{2}},\ldots,p_{\sigma_{1}^{2}}$. This sequence interpolates smoothly between the prior distribution (e.g., $N(0,\sigma_{N}^{2}I_{d})$) and the data distribution $P_{\textup{data}}$; such an “annealing” or “homotopy” method is required in practice to generate good samples [Son+20a].

Examples of SGM’s.

There have been several instantiations of this general approach. [SE19] add gaussian noise to the data and then use Langevin diffusion at a discrete set of noise levels $\sigma_{N}>\cdots>\sigma_{1}$ as the sampling algorithm. [Son+20a] take the continuous perspective and consider a more general framework, where the forward process can be any reasonable SDE. Then a natural reverse SDE evolves the final distribution $p_{\sigma_{N}^{2}}$ back to the data distribution; this process can be simulated with the estimated score. They consider methods based on two different SDE’s: score-matching Langevin diffusion (SMLD) based on adding Gaussian noise and denosing diffusion probabilistic models (DDPM) [Soh+15, HJA20], based on the Ornstein-Uhlenbeck process. Note that a difference with MCMC-based methods is that these SDE’s are evolved for a fixed amount of time, rather than until convergence. However, they can be combined with MCMC-based methods such as Langevin diffusion in the predictor-corrector approach for improved convergence. [DVK21] include Hamiltonian dynamics: they augment the state space with a velocity variable and consider a critically-damped version of the Ornstein-Uhlenbeck process. Finally, we note the work of [De ̵+21], who introduce the Diffusion Schrödinger Bridge method to learn a diffusion that more quickly transforms the prior into the data distribution.

We will give a general analysis framework for SGM’s that applies to the algorithms in both [SE19] and [Son+20a].

2 Prior work and challenges for theory

Although the literature on convergence for Langevin Monte Carlo [DM17, CB18, Che+18, Dal17, DK19, MMS20, EHZ21] and related sampling algorithms is extensive, prior works mainly consider the case of exact or stochastic gradients. In contrast, by the structure of the loss function (1), the score function learned in SGM is only accurate in $L^{2}(p)$. This poses a significant challenge for analysis, as the stationary distribution of Langevin diffusion with $L^{2}(p)$-accurate gradient can be arbitrarily far from $p$ (see Appendix D). Hence, any analysis must be utilizing the short/medium-term convergence, while overcoming the potential issue of long-term behavior of convergence to an incorrect distribution.

[BMR20] give the first theoretical analysis of SGM, and in particular, Langevin Monte Carlo with $L^{2}(p)$-accurate gradients. First, they show using uniform generalization bounds that optimizing the de-noising autoencoder (DAE) objective does in fact give a $L^{2}(p)$-accurate score function, with sample complexity depending on the complexity of the function class. They analyze convergence of LMC in Wasserstein distance. However, the error they obtain (Theorem 13) only decreases as $\varepsilon^{1/d}$ where $\varepsilon$ is the accuracy of the score estimate—so it suffers from the curse of dimensionality—and increases exponentially in the time that the process is run, the dimension, and the smoothness of the distribution, as in ODE/SDE discretization arguments that do not depend on contractivity.

[De ̵+21] give an analysis for [Son+20a] in TV distance that requires a $L^{\infty}$-accurate score function and depends exponentially on the amount of time the reverse SDE is run. Although exponential dependence is bad in general, it is mollified using their Diffusion Schrödinger Bridge (DSB) approach, as it allows running for a shorter, fixed amount of time, before the forward SDE converges to the prior distribution. However, this supposes that a good solution can be found for the DSB problem, and theoretical guarantees may be difficult to obtain.

We overcome the challenges of analysis with a $L^{2}(p)$-accurate gradient, and give the first analysis with only polynomial dependence on running time, dimension, and smoothness of the distribution, with rates that are a fixed power of $\varepsilon$. Our convergence result is in TV distance. We assume only smoothness conditions and a bounded log-Sobolev constant of the data distribution, a weaker condition than the dissipativity condition required by [BMR20]. We introduce a general framework for analysis of sampling algorithms given $L^{2}$-accurate gradients (score function) based on constructing a “bad set” with small measure and showing convergence of the discretized process conditioned on not hitting the bad set. We use our framework to give an end-to-end analysis for both the algorithms in [SE19] and [Son+20a], and illuminate the relative performance of different methods in practice.

3 Notation and organization

In Section 2 we explain our main results for Langevin Monte Carlo with $L^{2}(p)$-accurate score estimate and use it to derive convergence bounds for the annealed LMC method of [SE19]. In Section 3, we give our main results for the predictor-corrector algorithms of [Son+20a] based on simulating reverse SDE’s. Our proofs are based on a common framework which we introduce in Section 4. Full proofs are in the appendix.

Results for Langevin dynamics with estimated score

We make the following assumptions on the density $p$ and the score estimate $s$, which we will use throughout this paper.

$\ln p$ is $C^{1}$ and $L$-smooth, that is, $\nabla\ln p$ is $L$-Lipschitz. We assume $L\geq 1$.

$p$ satisfies a log-Sobolev inequality with constant $C_{\textup{LS}}$. We assume $C_{\textup{LS}}\geq 1$.

We note that the uniform Lipschitzness assumption (1) helps ensure a unique strong solution to the Langevin diffusion, as in [BMR20]. One special case where one can prove Lipschitzness for all $t$ is when $p_{0}$ is strongly log-concave [Lee+21, Lemma 28]. Although satisfying a log-Sobolev inequality (3) is a significant assumption, it is standard for analysis of Langevin Monte Carlo [VW19]. It is much weaker than assumptions in previous works [BMR20], including log-concave distributions and distributions satisfying strong dissipativity, and is stable under bounded perturbations. See Section E.1 for background on functional inequalities.

$s$ is a $C^{1}$ function that is $L_{s}$-Lipschitz. We assume $L_{s}\geq 1$.

The error in the score estimate is bounded in $L^{2}$:

Our first main result gives an error bound between the sampled distribution and $p$, assuming $L^{2}$-accurate score function estimate.

then running (LMC-SE) with score estimate $s$, step size h=\Theta\Bigl{(}\frac{\varepsilon_{\chi}^{2}}{dL^{2}C_{\textup{LS}}}\Bigr{)}, and time T=\Theta\Bigl{(}C_{\textup{LS}}\ln\bigl{(}\frac{2K_{\chi}}{\varepsilon_{\chi}^{2}}\bigr{)}\Bigr{)} results in a distribution $p_{T}$ such that $p_{T}$ is $\varepsilon_{\operatorname{TV}}$-far in TV distance from a distribution $\overline{p}_{T}$, where $\overline{p}_{T}$ satisfies $\chi^{2}(\overline{p}_{T}||p)\leq\varepsilon_{\chi}^{2}.$ In particular, taking $\varepsilon_{\chi}=\varepsilon_{\operatorname{TV}}$, we have the error guarantee that $\operatorname{TV}(p_{T},p)\leq 2\varepsilon_{\operatorname{TV}}$.

Note that the error bound is only achieved when running LMC for a moderate time; this is consistent with the fact that the stationary distribution of LMC with a $L^{2}$-score estimate can be arbitrarily far from $p$. Note also that we need a warm start in $\chi^{2}$-divergence: to obtain fixed errors $\varepsilon_{\operatorname{TV}},\varepsilon_{\chi}$, the required accuracy for the score estimate is inversely proportional to $K_{\chi}$. Intuitively, we must suffer from such a dependence because if the starting distribution is very far away, then there is no guarantee that $\|\nabla\ln p(x_{t})-s(x_{t})\|^{2}$ is small on average during the sampling algorithm. Finally, although we can state a result purely in terms of TV distance, we need this more precise formulation to prove a result for annealed Langevin dynamics.

2 Annealed Langevin dynamics with estimated score

In light of the warm start requirement in Theorem 2.1, we typically cannot directly sample from $p_{\textup{data}}$ or its approximation. Hence, [SE19] proposed using annealed Langevin dynamics: consider a sequence of noise levels $\sigma_{N}>\cdots>\sigma_{1}\approx 0$ giving rise to a sequence of distributions $p_{\sigma_{N}^{2}},\ldots,p_{\sigma_{1}^{2}}\approx p_{\textup{data}}$, where $p_{\sigma^{2}}=p*\varphi_{\sigma^{2}}$, $\varphi_{\sigma^{2}}$ being the density of $N(0,\sigma^{2}I_{d})$. For large enough $\sigma_{N}$, $\varphi_{\sigma_{N}^{2}}\approx p_{\sigma_{N}^{2}}$ provides a warm start to $p_{\sigma_{N}^{2}}$. We then successively run LMC using score estimates for $p_{\sigma_{k}^{2}}$, with the approximate sample for $p_{\sigma_{k}^{2}}$ giving a warm start for $p_{\sigma_{k-1}^{2}}$. We obtain the following algorithm and error estimate.

then $x^{(1)}$ is a sample from a distribution $q$ such that $\operatorname{TV}(q,p_{\sigma_{1}^{2}})\leq\varepsilon_{\operatorname{TV}}$.

Note that we assume a score estimate with error $\varepsilon$ at all noise scales; this corresponds to using an objective function that is a maximum of the score-matching objective over all noise levels, rather than an average over all noise levels as more commonly used in practice. However, these two losses are at most a factor of $M$ apart.

The proof shows that the noise levels $\sigma_{k}$ can be chosen as a geometric sequence, which matches the choice used in practice [SE20]. The additional dependence on $d$ and $\varepsilon_{\operatorname{TV}}$ in Theorem 2.2 compared to Theorem 2.1 comes from requiring a sequence of $\widetilde{O}(\sqrt{d})$ noise levels and an additional factor in $\chi^{2}$-divergence we suffer at the beginning of each level $m$. In the next section, we will find that using a reverse SDE to evolve the samples between the noise levels—called a predictor step—will improve the rate and time complexity.

Results for reverse SDE’s with estimated score

To improve the empirical performance of score-based generative modeling, [Son+20a] consider a general framework where noise is injected into a data distribution $p_{\textup{data}}$ via a forward SDE,

where $\widetilde{x}_{0}\sim\widetilde{p}_{0}:=p_{\textup{data}}$. Let $\widetilde{p}_{t}$ denote the distribution of $\widetilde{x}_{t}$ ($\widetilde{p}_{t}$ is used instead of $p_{t}$ to distinguish with the Gaussian-convolved distribution used in Annealed Langevin dynamics as in §2.2). Remarkably, $\widetilde{x}_{t}$ also satisfies a reverse-time SDE,

The case where $f\equiv 0$ and $g\equiv 1$ recovers the simple case of convolving with a Gaussian as used in §2.2; note, however that the reverse-time SDE differs from Langevin diffusion in having a larger (and time-varying) drift relative to the diffusion. [Son+20a] highlight the following two special cases. We will focus on DDPM while noting that our analysis applies more generically.

Score-matching Langevin diffusion: $f\equiv 0$. In this case, $\widetilde{p}_{t}=\widetilde{p}_{0}*\varphi_{\int_{0}^{t}g(s)^{2}\,ds}$, so [Son+20a] call this a variance-exploding (VE) SDE. As is common for annealing-based algorithms, [SE19, Son+20a] suggest choosing an exponential schedule, so that $g(t)=ab^{t}$ for constants $a,b$. We take $p_{\textup{prior}}=N(0,\int_{0}^{T}g(s)^{2}\,ds\cdot I_{d})$.

Denoising diffusion probabilistic modeling: $f(x,t)=-\frac{1}{2}g(t)^{2}x$. This is an Ornstein-Uhlenbeck process with time rescaling, $\widetilde{p}_{t}=M_{-\frac{1}{2}\int_{0}^{t}g(s)^{2}\,ds\sharp}\widetilde{p}_{0}*\varphi_{1-e^{-\int_{0}^{t}g(s)^{2}\,ds}}$, where $M_{\alpha}(x)=\alpha x$. [Son+20a] call this a variance-preserving (VP) SDE, as the variance converges towards $I_{d}$. Because it displays exponential convergence towards $N(0,I_{d})$, it can be run for a smaller amount of normalized time $\int_{0}^{t}g(s)^{2}\,ds$. [Son+20a] suggest the choice $g(t)=\sqrt{b+\alpha t}$. We take $p_{\textup{prior}}=N(0,(1-e^{-\int_{0}^{t}g(s)^{2}\,ds})I_{d})\approx N(0,I_{d})$.

To obtain an algorithm, we consider the following discretization and approximation of (4); note that in all cases of interest the integrals can be analytically evaluated. We reverse time so that $t$ corresponds to $T-t$ of the forward process. As we are free to rescale time in the SDE, we assume without loss of generality that the step sizes are constant. The predictor step is

running (4) starting from $p_{\textup{prior}}$ for time $T=\Theta\left({\ln(C_{\textup{LS}}d)\vee C_{\textup{LS}}\ln\left({\frac{1}{\varepsilon_{\operatorname{TV}}}}\right)}\right)$ and step size $h=\Theta\left({\frac{\varepsilon_{\operatorname{TV}}^{2}}{C_{\textup{LS}}(C_{\textup{LS}}+d)(L\vee L_{s})^{2}}}\right)$ results in a distribution $q_{T}$ so that $\operatorname{TV}(q_{T},p_{\textup{data}})\leq\varepsilon_{\operatorname{TV}}$.

A more precise statement of the Theorem can be found in the Appendix. Although we state our theorem for DDPM, we describe in Appendix C how it can be adapted to other SDE’s like SMLD and the sub-VP SDE; the primary SDE-dependent bound we need is a bound on $\nabla\ln\frac{\widetilde{p}_{t}}{\widetilde{p}_{t+h}}$. Because the predictor is tracking a changing distribution $p_{t}$, we incur more error terms and worse dependence on parameters ($C_{\textup{LS}},L$) than in LMC (Theorem 2.1). Motivated by this, we intersperse the predictor steps with LMC steps—called corrector steps in this context—to give additional time for the process to mix, resulting in improved dependence on parameters.

Keep the setup of Theorem 3.1. Then for $\varepsilon_{\operatorname{TV}}^{3}=O\left({\frac{1}{(1+L_{s}/L)^{2}(1+C_{\textup{LS}}/d)(\ln(C_{\textup{LS}}d)\vee C_{\textup{LS}})}}\right)$, if

then Algorithm 2 with appropriate choices of $T=\Theta\left({\ln(C_{\textup{LS}}d)\vee C_{\textup{LS}}\log\left({\frac{1}{\varepsilon_{\operatorname{TV}}}}\right)}\right)$, $N_{m}$, corrector step sizes $h_{m}$ and predictor step size $h$, produces a sample from a distribution $q_{T}$ such that $\operatorname{TV}(q_{T},p_{\textup{data}})<\varepsilon_{\operatorname{TV}}$.

The assumption on $\varepsilon_{\operatorname{TV}}$ is for convenience in stating our bound. In comparison to using the predictor step alone (Theorem 3.1), note that in the bound on $\varepsilon$, we obtain the improved rate of the corrector step as in Theorem 2.1; this is because the predictor step only needs to track the actual distribution in $\chi^{2}$-divergence with error $O(1)$, and the final corrector steps are responsible for decreasing the error to $\varepsilon_{\operatorname{TV}}$. In comparison to the Annealed Langevin sampler (Algorithm 1, Theorem 2.2), which can be viewed as using the corrector step alone, adding a predictor step provides a better warm start for the distribution at the next smaller noise level, resulting in better dependence on parameters. Thus the predictor-corrector algorithm combines the strengths of the predictor and corrector steps. For real-world data, it can be challenging to estimate TV-distance between distributions given only samples, and hence difficult to check consistency with empirical observations. However, our claim that using a corrector can improve the convergence rate of DDPM/SMLD is consistent with the simulation results in Section 4.2 of [Son+20a].

Theoretical framework and proof sketches

The main idea of our analysis framework is to convert a $L^{2}$ error guarantee to a $L^{\infty}$ error guarantee by excluding a bad set, formalized in the following theorem.

If $Z_{k}\in B_{k}^{c}$ for all $0\leq k\leq n-1$, then $Z_{n}=\overline{Z}_{n}$. (For $n=0$, this says $Z_{0}=\overline{Z}_{0}$.)

$\chi^{2}(\overline{q}_{n}||p_{n})\leq D_{n}^{2}$.

For our setting, we will take the “bad sets” $B_{n}$ to be the set of $x$ where $\left\|{s_{\theta}(x)-\nabla\ln p}\right\|$ is large, $q_{n}$ to be the discretized process with estimated score, and $\overline{q}_{n}$ to be the discretized process with estimated score except in $B_{n}$ where the error is large. Because $\overline{q}_{n}$ uses an $L^{\infty}$-accurate score estimate, we can use existing techniques for analyzing Langevin Monte Carlo [VW19, EHZ21, Che+21] to bound $\chi^{2}(\overline{q}_{n}||p_{n})$.

First note that if some $Z_{k}\in B_{k}$ for $0\leq k\leq n-1$, then for the smallest such $k$, we have $\overline{Z}_{k}=Z_{k}\in B_{k}$; the same is true if $\overline{Z}_{k}\in B_{k}$ for some $0\leq k\leq n-1$. We then bound using condition 1 and Cauchy-Schwarz:

The second inequality then follows from the triangle inequality and Cauchy-Schwarz:

It now remains to give $\chi^{2}$ convergence bounds under $L^{\infty}$-accurate score estimate. The following theorem may be of independent interest.

Following [Che+21], we prove this by first defining a continuous-time interpolation $q_{t}$ of the discrete process, and then deriving a differential inequality for $\chi^{2}(q_{t}||p)$ using the log-Sobolev inequality for $p$. Compared to [Che+21], we incur an extra error term arising from the inaccurate gradient.

This allows us to sketch the proof of Theorem 2.1; a complete proof is in Section B.

We first define the bad set where the error in the score estimate is large,

for some $\varepsilon_{1}$ to be chosen. Then by Chebyshev’s inequality, $P(B)\leq\left({\frac{\varepsilon}{\varepsilon_{1}}}\right)^{2}=:\delta$. Let $\overline{q}_{nh}$ be the discretized process, but where the score estimate is set to be equal to $\nabla\ln p$ on $B$; note it agrees with $q_{nh}$ as long as it has not hit $B$. Because $\overline{q}_{nh}$ uses a score estimate that has $L^{\infty}$-error $\varepsilon_{1}$, Theorem 4.2 gives a bound for $\chi^{2}(\overline{q}_{Nh}||p)$. Then Theorem 4.1 gives

The theorem then follows from choosing parameters so that $\chi^{2}(\overline{q}_{T}||p)\leq\varepsilon_{\chi}^{2}$ and $\operatorname{TV}(q_{T},\overline{q}_{T})\leq\varepsilon_{\operatorname{TV}}$. ∎

We remark that the main inefficiency in the proof comes from the use of Chebyshev’s inequality, and a $L^{p}$ bound on the error for $p>2$ will improve the bound.

Choosing the sequence $\sigma_{1}<\cdots<\sigma_{M}$ to be geometric with ratio $1+\frac{1}{\sqrt{d}}$ ensures that the $\chi^{2}$-divergence between successive distributions $p_{\sigma_{m}^{2}}$ is $O(1)$. Then, choosing $\sigma_{M}^{2}=\Omega(C_{\textup{LS}}d)$ ensures we have a warm start for the highest noise level: $\chi^{2}(p_{\textup{prior}}||p_{\sigma_{M}^{2}})=O(1)$. This uses $O\left({\sqrt{d}\log\left({\frac{dC_{\textup{LS}}}{\sigma_{\min}^{2}}}\right)}\right)$ noise levels. Chebyshev’s inequality can be used to show that the distribution of the final sample $x^{(m)}$ for $p_{\sigma_{m}^{2}}$ is $O(\varepsilon_{\operatorname{TV}}/M)$ close to a distribution that is $O(M/\varepsilon_{\operatorname{TV}})$ in $\chi^{2}$-divergence from $p_{\sigma_{m+1}^{2}}$. This gives the warm start parameter $K_{\chi}=(M/\varepsilon_{\operatorname{TV}})^{1/2}$; substituting into Theorem 2.1 then gives the required bound for $\varepsilon$. Note that the TV errors accrued from each level add to $O(\varepsilon_{\operatorname{TV}})$. ∎

To analyze the predictor-based algorithms, we also first prove convergence bounds under $L^{\infty}$-accurate score estimate.

Consider DDPM with $g\equiv 1$, $T\geq 1\vee\ln(C_{\textup{LS}}d)$, and $h=O\left({\frac{1}{C_{\textup{LS}}(d+C_{\textup{LS}})(L\vee L_{s})^{2}}}\right)$. (Recall that $p_{kh}$ and $q_{kh}$ are the $k$-th iterate of LMC with step size $h$ and true/estimated score respectively.) Then

and if $\varepsilon_{1}<\frac{1}{128C_{\textup{LS}}}$,

Moreover, for $q_{0}=p_{\textup{prior}}$, $\chi^{2}(q_{0}||p_{0})\leq e^{-T/2}C_{\textup{LS}}d$.

We give a more precise statement in Section C. Note that unlike the case for LMC as in Theorem 4.2, the base density $p_{t}$ is also evolving in time, which produces additional error terms and necessitates a more involved analysis. The additional error terms can be bounded using the Donsker-Varadhan variational principle, concentration for distributions satisfying LSI, and error bounds between $p_{t}$ and $p_{t+h}$ for small $h$.

Here, we only state the result about DDPM, which has better bounds than SMLD (when $g\equiv 1$) because both the forward and backwards processes exhibit better mixing properties: the warm start improves exponentially rather than inversely with $T$, and the log-Sobolev constant is uniformly bounded by that of $p_{\textup{data}}$ rather than increasing. However, the analysis in Section C can be directly applied to SMLD and other models as well. We also note there is a sense in which DDPM and SMLD are equivalent under a rescaling in time and space (see discussion in Section C.2).

Note that the choice of $h$ is necessary for exponential decay of error; as if $h$ is not small enough, we would get an exponential growing instead of decaying factor in the one-step error (See Section C for details). Such an $h$ may however still be a suitable choice when used in conjunction with a corrector step. Moreover, as $\varepsilon_{1}\to 0$, with appropriate choice of $T$ and $h$, $q_{Nh}$ and $p_{Nh}$ can be made arbitrarily close.

Theorem 3.1 now follows from the $L^{\infty}$ result (Theorem 4.3) in the same way that Theorem 2.1 follows from Theorem 4.2.

To prove Theorem 3.2, it suffices to run the corrector steps only at the lowest noise level, that is, set $N_{m}=0$ for $1\leq m<T/h$, although we note that interleaving the predictor and corrector steps does empirically help with mixing. The proof follows from using the predictor and the corrector theorems in series: first apply Theorem 3.1 with $\varepsilon_{\chi}=O(1)$ to show that the predictor results a warm start $p_{\textup{data}}$, then use Theorem 2.1 to show the corrector reduces the error to the desired $\varepsilon_{\operatorname{TV}}$.

Conclusion

We introduced a general framework to analyze SDE-based sampling algorithms given a $L^{2}$-error score estimate, and used it to obtain the first convergence bounds for several score-based generative models with polynomial complexity in all parameters. Our analysis can potentially be adapted to other SDE’s and sampling algorithms beyond Langevin Monte Carlo. There is also room for improving our analysis to better use smoothing properties of the SDE’s and compare different choices of the diffusion speed $g$.

We present several interesting further directions to explore. In addition to extending the analysis to other SGM’s and comparing their theoretical performance (relative to each other as well as other approaches to generative modeling), we propose the following.

Our assumption of a bounded log-Sobolev constant essentially limits the analysis to distributions that are close to unimodal. However, SGM’s are empirically successful at modeling multimodal distributions [SE19], and in fact perform better with multimodal distributions than other approaches such as GAN’s. Can we analyze the convergence for simple multimodal distributions, such as a mixture of distributions each with bounded log-Sobolev constant? Positive results on sampling from multimodal distributions such as [GLR18] suggest this is possible, as the sequence of noised distributions is natural for annealing and tempering methods (see [GLR18, Remark 7.2]).

Weakening conditions on the score estimate.

The assumption that we have a score estimate that is $O(1)$-accurate in $L^{2}$, although weaker than the usual assumptions for theoretical analysis, is in fact still a strong condition in practice that seems unlikely to be satisfied (and difficult to check) when learning complex distributions such as distributions of images. What would a reasonable weaker condition be, and in what sense can we still obtain reasonable samples?

Guarantees for learning the score function.

Our analysis assumes a $L^{2}$-estimate of the score function is given, but the question remains of when we can find such an estimate. What natural conditions on distributions allow their score functions to be learned by a neural network? Various works have considered the representability of data distributions by diffusion-like processes [TR19], but the questions of optimization and generalization appear more challenging.

Acknowledgements

We thank Andrej Risteski for helpful conversations. This work was done in part while HL was visiting the Simons Institute for the Theory of Computing. The work was supported in part by National Science Foundation via awards DMS-2012286 and CCF-1934964 (Duke Tripods).

References

Appendix A Computations

We start the proofs by collecting some preliminary results. In the following, we will consider the SDE

and the interpolation of the discretization of an approximation

when $t\geq t_{-}$. Let $P_{t}$ and $Q_{t}$ denote the law of $x_{t}$ and $z_{t}$, respectively. We will take $t_{-}=kh$ and $t\in[kh,(k+1)h)$. We will assume that $f,\widehat{f},G$ are continuous and the functions $f(\cdot,t)$, $\widehat{f}(\cdot,t)$ are uniformly Lipschitz for each $t\in[kh,(k+1)h]$.

In this section, we will make some computations that will be used in both Sections B and C. First, we derive how the density evolves in time.

Let $Q_{t}$ denote the law of the interpolated process (8). Then

Let $q_{t|t_{-}}$ denote the distribution of $z_{t}$ conditioned on $z_{t_{-}}$. Then the Fokker-Planck equation gives

Taking expectation with respect to $z_{t_{-}}$ we get

Note that for fixed $z$, $\int q_{t_{-}|t}(y|z)dy=1$. Hence

We now use Lemma A.1 to compute how the $\chi^{2}$-divergence between the approximate and exact densities changes. The following generalizes the calculation of [EHZ21] in the case where $x_{t}$ is a non-stationary stochastic process. For simplicity of notation, from now on, we wil consider the case $G(t)$ being a scalar.

Let $P_{t}$ and $Q_{t}$ be the laws of (7) and (8) for $G(t)=g(t)I_{d}$. Then

For the second term, using integration by parts,

Finally, we will make good use of the following lemma to bound the second term in Lemma A.2.

Appendix B Analysis for LMC

then running (LMC-SE) with score estimate $s$ and step size $h=\frac{\varepsilon_{\chi}^{2}}{2720dL^{2}C_{\textup{LS}}}$ for any time $T\in[T_{\min},C_{T}T_{\min}]$, where $T_{\min}=4C_{\textup{LS}}\ln\left({\frac{2K_{\chi}}{\varepsilon_{\chi}^{2}}}\right)$, results in a distribution $p_{T}$ such that $p_{T}$ is $\varepsilon_{\operatorname{TV}}$-far in TV distance from a distribution $\overline{p}_{T}$, where $\overline{p}_{T}$ satisfies $\chi^{2}(\overline{p}_{T}||p)\leq\varepsilon_{\chi}^{2}.$ In particular, taking $\varepsilon_{\chi}=\varepsilon_{\operatorname{TV}}$, we have the error guarantee that $\operatorname{TV}(p_{T},p)=2\varepsilon_{\operatorname{TV}}$.

The main difficulty is that the stationary distribution of LMC using the score estimate may be arbitrarily far from $p$, even if the $L^{2}$ error of the score estimate is bounded. (See Section D.) Thus, a long-time convergence result does not hold, and an upper bound on $T$ is required, as in the theorem statement.

We instead proceed by showing that conditioned on not hitting a bad set, if we run LMC using $s$, the $\chi^{2}$-divergence to the stationary distribution will decrease. This means that the closeness of the overall distribution (in TV distance, say) will decrease in the short term, despite it will increase in the long term, as the probability of hitting the bad set increases. This does not contradict the fact that the stationary distribution is different from $p$. By running for a moderate amount of time (just enough for mixing), we can ensure that the probability of hitting the bad set is small, so that the resulting distribution is close to $p$. Note that we state the theorem with a $C_{T}$ parameter to allow a range of times that we can run LMC for.

More precisely, we prove Theorem B.1 in two steps.

Defining a bad set and bounding the hitting time (Section B.2).

The idea is now to reduce to the case of $L^{\infty}$ error by defining the “bad set” $B$ to be the set where $\left\|{s-\nabla f}\right\|\geq\varepsilon_{1}$, where $\varepsilon\ll\varepsilon_{1}\ll 1$. This set has small measure by Chebyshev’s inequality. Away from the bad set, Theorem 4.2 applies; it then suffices to bound the probability of hitting $B$. Technically, we define a coupling with a hypothetical process where the $L^{\infty}$ error is always bounded, and note that the processes disagree exactly when it hits $B$; this is the source of the TV error.

We consider the probability of being in $B$ at times $0,h,2h,\ldots$. we note that Theorem B.1 bounds the $\chi^{2}$-divergence of this hypothetical process $X_{t}$ at time $t$ to $p$. If the distribution were actually $p$, then the probability $X_{t}\in B^{\prime}$ is exactly $p(B^{\prime})$; we expect the probability to be small even if the distribution is close to $p$. Indeed, by Cauchy-Schwarz, we can bound the probability $X\in B$ in terms of $P(B)$ and $\chi^{2}(q_{t}||p)$; this bound is given in Theorem 4.1. Note that the eventual bound depends on $\chi^{2}(q_{t}||p)$, so we have to assume a warm start, that is, a reasonable bound on $\chi^{2}(q_{0}||p)$.

The following gives a long-time convergence bound for LMC with inaccurate gradient, with error bounded in $L^{\infty}$; this may be of independent interest.

See 4.2 Following [Che+21], convergence in Rényi divergence can also be derived; we only consider $\chi^{2}$-divergence because we will need a warm start in $\chi^{2}$-divergence for our application. Note that by letting $N\to\infty$ and $h\to 0$, we obtain the following.

Keep the assumptions in Theorem 4.2. The stationary distribution $q$ of Langevin diffusion with score estimate $s$ satisfies

We follow the proof of [Che+21, Theorem 4], except that we work with the $\chi^{2}$ divergence directly, rather than the Rényi divergence, and have an extra term from the inaccurate gradient (17). Given $t\geq 0$, let $t_{-}=h\left\lfloor\frac{t}{h}\right\rfloor$. Define the interpolated process by

and let $q_{t}$ denote the distribution of $X_{t}$ at time $t$, when $X_{0}\sim q_{0}$.

where $A_{1},A_{2},A_{3}$ are obtained by substituting in the 3 terms in (13), and given in (15), (16), and (17). Let $V(x)=-\ln p(x)$. We consider each term in turn.

when $h^{2}\leq\frac{1}{288L^{2}}$. By [Che+21, p. 15]

when $h\leq\frac{1}{1536L^{2}C_{\textup{LS}}}$. Finally,

Combining (12), (14), (15), (16), and (17) gives

if $h\leq\left({\frac{1}{12\cdot 72dL^{3}C_{\textup{LS}}}}\right)^{1/2}\wedge\frac{1}{12\cdot 336dC_{\textup{LS}}}$ and $\varepsilon_{1}\leq\left({\frac{1}{48C_{\textup{LS}}}}\right)^{1/2}$. Then for $t\in[kh,(k+1)h)$,

using $h\leq\frac{1}{12\sqrt{2}L}$. Unfolding the recurrence and summing the geometric series gives

when $h\leq\frac{1}{1360dL^{2}C_{\textup{LS}}}$ and $\varepsilon_{1}^{2}\leq\frac{1}{40C_{\textup{LS}}}$. We can check that the given condition on $h$ and the fact that $LC_{\textup{LS}}\geq 1$ (Lemma E.5) imply all the required inequalities on $h$. ∎

B.2 Proof of Theorem B.1

We first define the bad set where the error in the score estimate is large,

Given $t\geq 0$, let $t_{-}=h\left\lfloor\frac{t}{h}\right\rfloor$. Given a bad set $B$, define the interpolated process by

In other words, run LMC using the score estimate as long as the point is in the good set at the previous discretization step, and otherwise use the actual gradient $\nabla\ln p$. Let $\overline{q}_{t}$ denote the distribution of $\overline{z}_{t}$ when $\overline{z}_{0}\sim q_{0}$; note that $q_{nh}$ is the distribution resulting from running LMC with estimate $b$ for $n$ steps and step size $h$. Note that this auxiliary process is defined only for purposes of analysis; it cannot be used for practical algorithm as we do not have access to $\nabla f$.

We can couple this process with LMC using $s$ so that as long as $X_{t}$ does not hit $B$, the processes agree, thus satisfying condition 1 of Theorem 4.1.

For this to be bounded by $\varepsilon_{\chi}^{2}$, it suffices for the terms to be bounded by $\frac{\varepsilon_{\chi}^{2}}{2},\frac{\varepsilon_{\chi}^{2}}{4},\frac{\varepsilon_{\chi}^{2}}{4}$; this is implied by

(We choose $h$ so that the condition in Theorem 4.2 is satisfied; note $\varepsilon_{\chi}\leq 1$.) By Theorem 4.1,

In order for this to be $\leq\varepsilon_{\operatorname{TV}}$, it suffices for

Supposing that we run for time $T$ where $T_{\min}\leq T\leq C_{T}T_{\min}$, we have that $N=\frac{T}{h}\leq\frac{C_{T}T_{\min}}{h}$. Thus it suffices for

B.3 Proof of Theorem 2.2

We restate the theorem for convenience. See 2.2

Choose the sequence $\sigma_{\min}^{2}=\sigma_{1}^{2}<\cdots<\sigma_{M}^{2}$ to be geometric with ratio $1+\Theta\left({\frac{1}{\sqrt{d}}}\right)$. Note that

For $\sigma_{2}^{2}=(1+\varepsilon)\sigma_{1}^{2}$, this equals $(1+\varepsilon)^{-d/2}(1-\varepsilon)^{-d/2}=(1-\varepsilon^{2})^{-d/2}-1$. For $\varepsilon=\Theta\left({\frac{1}{\sqrt{d}}}\right)$, this is $d\cdot O\left({\frac{1}{d}}\right)=O(1)$. Hence, the $\chi^{2}$-divergence between successive distributions $p_{\sigma_{m}^{2}}$ is $O(1)$. Choosing $\sigma_{M}^{2}=\Omega(d(M_{1}+C_{\textup{LS}}))$ ensures we have a warm start for the highest noise level by Lemma E.9: $\chi^{2}(p_{\textup{prior}}||p_{\sigma_{M}^{2}})=O(1)$. This uses $O\left({\sqrt{d}\log\left({\frac{dC_{\textup{LS}}}{\sigma_{\min}^{2}}}\right)}\right)$ noise levels.

Write $p_{m}=p_{\sigma_{m}^{2}}$ for short. Let $q_{m}$ be the distribution of the final sample $x^{(m)}$. We show by downwards induction on $m$ that there is $\overline{q}_{m}$ such that

For $m=M$, this follows from the assumption on $\varepsilon$ and Theorem 2.1 with $K_{\chi}=O(1)$ (given by the warm start).

Fix $m<M$ and suppose it holds for $m+1$. We use the closeness between $q_{m+1}$ and $p_{m+1}$ combined with $\chi^{2}(p_{m+1}||p_{m})=O(1)$ to obtain compute how close $q_{m+1}$ and $p_{m}$ are. Because the triangle inequality does not hold for $\chi^{2}$, we will incur an extra TV error.

Let $\overline{q}_{m,m+1}$ be the distribution of the final sample if $x^{(m+1)}_{0}\sim\overline{q}_{m}$. We have $\operatorname{TV}(q_{m+1},\overline{q}_{m,m+1})\leq\operatorname{TV}(q_{m},\overline{q}_{m})\leq\frac{(M+1)-m}{M+1}\varepsilon_{\operatorname{TV}}$.

By Markov’s inequality, when $\chi^{2}(p_{m+1}||p_{m})\leq 1$,

so $\overline{q}_{m+1,m}\leq 2\overline{q}_{m+1}$ and

Let $\overline{q}_{m+1,m}^{\prime}$ be the distribution of $x_{N_{m}}^{(m)}$ when $x_{0}^{(m)}\sim\overline{q}_{m+1,m}$. Then by assumption on $\varepsilon$ (3) and Theorem 2.1 (with $K_{\chi}={4\sqrt{\frac{M+1}{\varepsilon_{\operatorname{TV}}}}}$, $\varepsilon_{\chi}=\frac{\varepsilon_{\operatorname{TV}}}{4(M+1)}$, and $\varepsilon_{\operatorname{TV}}\leftarrow\frac{\varepsilon_{\operatorname{TV}}}{2(M+1)}$), there is $\overline{q}_{m}$ such that $\operatorname{TV}(\overline{q}_{m,m+1}^{\prime},\overline{q}_{m})\leq\frac{\varepsilon_{\operatorname{TV}}}{2(M+1)}$ and $\chi^{2}(\overline{q}_{m}||p_{m})\leq{\frac{\varepsilon_{\operatorname{TV}}}{4(M+1)}}$. It remains to bound

where we use (19) in the last line. This finishes the induction step.

Finally, the theorem follows by taking $m=1$ and noting

Appendix C Analysis for SGM based on reverse SDE’s

In this section, we analyze score-based generative models based on reverse SDE’s. In Section C.2, we prove convergence of the predictor algorithm under $L^{\infty}$-accurate score estimate (Theorem 4.3, restated as C.1) using lemmas proved in Section C.3, C.4, C.5, and C.6. In Section C.7, we prove convergence of the predictor algorithm under $L^{2}$-accurate score estimate (Theorem 3.1, restated as C.16). In Section C.8, we prove convergence of the predictor-corrector algorithm (Theorem 3.2).

With a change of variable in (4), we define the sampling process $x_{t}$ on $[0,T]$ by

where $h=T/N$ is the step size and $\eta_{k}$ is a sequence of independent Gaussian random vectors. As we run (20) from to $N$ with $h$ small enough, we should expect that the distribution of $z_{T}$ is close to that of $x_{T}$. However, in both SMLD or DDPM models, for fixed $z_{k}$, the integration

can be exactly computed, as can the diffusion term. Therefore, we can consider the following process $z_{t}$ as an “interpolation” of (20):

Note that by running this process instead, we can reduce the discretization error. Now if we denote the distribution of $z_{t}$ by $q_{t}$, with $q_{0}\approx p_{0}$, we can expect that $q_{T}$ is close to $p_{T}$. Here the estimated score $s$ satisfies for all $x$

Observe that in either SMLD or DDPM, the function $g(t)^{2}$ is Lipschitz on $[0,T]$. So in the following sections, we will assume that $g(t)^{2}$ is $L_{g}$-Lipschitz on $[0,T]$.

C.2 Predictor

In this section, we present the main result (Theorem C.1) on the one-step error of the predictor in $\chi^{2}$-divergence, which can be obtained by directly applying the Gronwall’s inequality to the differential inequality derived in Lemma C.3. Note that Theorem C.1 is a more precise version of Theorem 4.3; see the remark following the theorem.

With the setting in Section C.1, assume $g$ is non-decreasing and let

where $C_{t}$ is the log-Sobolev constant of $p_{t}$, bounded in Lemma E.7. Suppose that $\nabla\ln p_{t}$ is $L$-Lipschitz for all $t\in[kh,(k+1)h]$, $s(\cdot,kh)$ is $L_{s}$-Lipschitz, $L,L_{s}\geq 1$, and $\varepsilon_{kh}$ is such that (22) holds. Then

are defined in (24), (30), (33), (35) and (36), respectively.

The theorem follows from applying Gronwall’s inequality to the result of Lemma C.3. ∎

Remark. Note that in DDPM, $E=O(L_{s}^{2}+L^{2}d)$. Therefore, when $g\equiv 1$, $C_{t,kh}=O(\varepsilon_{1}^{2}+(L_{s}^{2}+L^{2}d)h)$, where we denote the upper bound of $\varepsilon_{kh}$ for all $k\in\{0,...,N\}$ by $\varepsilon_{1}$. Using the bound on the log-Sobolev constant (Lemma E.7) and second moment (Lemma E.8) for DDPM, we note that the restriction on $h$ for all steps is implied by

with appropriate constants. Then we can conclude the first inequality in Theorem 4.3 by combining Theorem C.1 and Lemma E.7 and the second inequality from unfolding the first one and evaluating the geometric series. Likewise, we have the following analogue for SMLD, for which we omit the proof.

and letting $t=T-Nh$, if $\varepsilon_{1}<\frac{1}{128C_{T}}$,

Moreover, for $q_{0}=p_{\textup{prior}}$, $q_{0}=\varphi_{T}$, $\chi^{2}(q_{0}||p_{0})\leq\frac{C_{\textup{LS}}d}{T}$.

Remark. We note that in a sense SMLD and DDPM are equivalent, as we can get from one to the other by rescaling in time and space. First we recall that, as discussed in Section 3, all the SMLD models are equivalent under rescaling in time. Therefore we can assume $g(t)=e^{t/2}$ and consider the forward SDE for SMLD

where $w_{t}$ is a standard Brownian Motion. Now let $y_{t}=e^{-t/2}x_{t}$; then

which is exactly DDPM with $g(t)=1$. Note that Theorem C.2 uses a different parameterization for SMLD and the resulting complexity is slightly worse.

C.3 Differential Inequality

Now we prove a differential inequality involving $\chi^{2}(q_{t}||p_{t})$. As in [Che+21], the key difficulty is to bound the discretization error. We decompose it into two error terms and bound them in Lemma C.4 and Lemma C.5 separately.

Let $(q_{t})_{0\leq t\leq T}$ denote the law of the interpolation (21). With the setting in Lemma C.1, we have for $t\in[kh,(k+1)h]$,

where $C_{t}$ is the LSI constant of $p_{t}$, $\varepsilon_{kh}$ is the $L^{\infty}$-score estimation error at time $kh$ and $E$ is defined in (24).

Using the fact that $p_{t}$ satisfies a log-Sobolev inequality with constant $C_{t}$,

In the setting of Lemma C.3, we have the following bound for term $A$:

In SMLD, $f(x,t)=0$ and hence $A=0$; while in DDPM, $f(x,t)=-\frac{1}{2}g(t)^{2}x$. Therefore, by Lemma A.3,

In the setting of Lemma C.3, we have the following bound for term $B$:

Now we bound these error terms separately. For $B_{1}$, by the Lipschitz assumption, we have by Lemma A.3, for a constant $C_{2}>0$ to be chosen later,

For $B_{2}$, recalling the assumption that $\left\|{s(x,T-kh)-\nabla\ln p_{kh}(x)}\right\|\leq\varepsilon_{kh}$ for all $x$, we have by Lemma A.3

Now for the last error term $B_{3}$, we have by Lemma A.3 that

where $C_{t,L}$ and $C_{d,L}$ are constants defined in (30) and (33) respectively. Hence

Combining all these results, we finally obtain the bound for error term $B$ in Lemma C.3: for $h\leq\frac{1}{64C_{t,L}g(T-kh)^{2}}$,

C.4 Change of Measure

Define the Langevin diffusion w.r.t. $p(x)$:

Rearrange this inequality to obtain the desired result. ∎

In the setting of Lemma C.3, it holds that

Applying Lemma C.6 to the density $\psi_{t}q_{t}$ yields

Note that we cannot expect analogous results for a general $u(x)$ as in Lemma C.6. In the general case, we apply the Donsker-Varadhan variational principle, which states that for probability measures $p$ and $q$,

Towards this end, we first need to analyze $\operatorname{KL}(\psi_{t}q_{t}||p_{t})$.

Since $p_{t}$ satisfies LSI with constant $C_{t}$,

With this in hand, we are ready to bound the second moment of $\psi_{t}q_{t}$ as well as the variance of a Gaussian random vector with respect to this measure:

where $C_{t}$ is the LSI constant of $p_{t}$, which is bounded in Lemma E.6, and the second moment of $p_{t}$ is bounded in Lemma E.8.

Since $p_{t}$ has LSI constant $C_{t}$, by Donsker-Varadhan variational principle,

for any $s>0$. By Lemma E.1, for any $s\in[0,\frac{1}{C_{t}})$, we have

Now with the bound of $\operatorname{KL}(\psi_{t}q_{t}||p_{t})$ in Lemma C.8, we obtain

where $C_{t}$ is the LSI constant of $p_{t}$.

Note that $\int_{kh}^{t}g(T-s)dw_{s}$ is a Gaussian random vector with variance $\int_{kh}^{t}g(T-s)^{2}ds\cdot I_{d}$. Using the Donsker-Varadhan variational principle, for any random variable $X$,

where the last inequality is due to Lemma C.8. We have proved

C.5 Perturbation Error

where $p_{x,\sigma^{2}}$ denotes the probability density

where $y^{*}\in\operatorname{argmax}_{y}p_{x,\sigma^{2}}(y)$ is a mode of the distribution $p_{x,\sigma^{2}}$. We now bound each of these terms.

For the first term, note that $p_{x,\sigma^{2}}$ is $\left({\frac{1}{\sigma^{2}}-L}\right)$-strongly convex, so satisfies a Poincaré inequality with constant $\left({\frac{1}{\sigma^{2}}-L}\right)^{-1}$. Thus

For the second term, by Lemma E.3, noting that $V(y)+\frac{\left\|{x-y}\right\|^{2}}{2\sigma^{2}}$ is $\left({\frac{1}{\sigma^{2}}+L}\right)$-smooth,

where the last inequality uses $\sigma^{2}\leq\frac{1}{2L}$.

For the third term, we note that the mode satisfies

Putting these together and using $\left({\frac{1}{\sigma^{2}}-L}\right)^{-1}\leq 2$, we obtain

With the setting in Lemma C.11 and the notation $p_{\alpha}(x)=\alpha^{d}p(\alpha x)$ for $\alpha\geq 1$, we have that for $L\leq\frac{1}{2\alpha^{2}\sigma^{2}}$,

Without loss of generality, we can assume that $p(x)=e^{-V(x)}$; then $p_{\alpha}(x)=\alpha^{d}e^{-V(\alpha x)}$. Hence

Since $\alpha\nabla V(\alpha x)$ is $\alpha^{2}L$-Lipschitz, by Lemma C.11,

The result follows from combining the three inequalities above. ∎

In the setting of Lemma C.3, we have for $t\in[kh,(k+1)h]$,

In both SMLD and DDPM models, we have the following relationship for $t\in[kh,(k+1)h]$:

where $p_{\alpha}(x)=\alpha^{d}p(\alpha x)$. In SMLD, $\alpha=1$ and $\sigma^{2}=\int_{kh}^{t}g(T-s)^{2}\,ds$, while in DDPM, $\alpha=e^{\frac{1}{2}\int_{kh}^{t}g(T-s)^{2}\,ds}$ and $\sigma^{2}=1-e^{-\int_{kh}^{t}g(T-s)^{2}\,ds}$. Now for SMLD,

where in the last inequality we use the fact that $g$ is increasing, so that for $h\leq\frac{1}{4Lg(T-kh)^{2}}$,

Recall that to use Lemma C.11, it suffices that $L\leq\frac{1}{2\alpha^{2}\sigma^{2}}$, and so it suffices that $h\leq\frac{1}{4Lg(T-kh)^{2}}$ in SMLD.

For DDPM, observe that for $h\leq\frac{1}{4g(T-kh)^{2}}$,

By Lemma C.12, using the assumption that $L\geq 1$, we obtain

C.6 Auxiliary Lemmas

With the setting of Lemma C.3, we have the following bound of the second moment of estimated score function with respect to $\psi_{t}q_{t}$:

where $C_{t,L}$ and $C_{d,L}$ are constants defined in Lemma C.13.

Recall that we need to bound this second moment of estimated score function with respect to $\psi_{t}q_{T}$. For the first term, as $\left\|{s(x.T-kh)-\nabla\ln p_{kh}(x)}\right\|$ is $\varepsilon_{kh}$-bounded, we have trivial bound that

By Lemma C.13, the second term is bounded by

for constant $C_{t,L}$ and $C_{d,L}$ defined in (30) and (33) respectively. The last term is bounded in Corollary C.7 by

Combining these three inequalities, we obtain that for $h\leq\frac{1}{g(T-kh)^{2}}$,

where the next-to-last line is due to the fact that the estimated score function is $L_{s}$-Lipschitz. We also use the fact that $g(t)$ is an increasing function and hence $g(T-t)\leq g(T-kh)$ for any $t\in[kh,(k+1)h]$. Hence if $h\leq\frac{1}{3(L_{s}+1/2)g(T-kh)^{2}}$, then

Therefore, by the fact that $(a+b)^{2}\leq 2a^{2}+2b^{2}$ for any $a,b>0$,

With the results of Lemma C.14 and Lemma C.9, we have

Now plugging this and the result of Lemma C.10 into (34), we get that

C.7 Proof of Theorem 3.1

We state a more precise version of Theorem 3.1. The structure of the proof is similar to that of Theorem 2.1.

running (4) starting from $p_{\textup{prior}}$ for time $T=\Theta\left({\ln(C_{\textup{LS}}d)\vee C_{\textup{LS}}\ln\left({\frac{1}{\varepsilon_{\operatorname{TV}}}}\right)}\right)$ and step size $h=\Theta\left({\frac{\varepsilon_{\chi}^{2}}{C_{\textup{LS}}(C_{\textup{LS}}+d)(L\vee L_{s})^{2}}}\right)$ results in a distribution $q_{T}$ such that $q_{T}$ is $\varepsilon_{\operatorname{TV}}$-far in TV distance from a distribution $\overline{q}_{T}$, where $\overline{q}_{T}$ satisfies $\chi^{2}(\overline{q}_{T}||p_{\textup{data}})\leq\varepsilon_{\chi}^{2}$. In particular, taking $\varepsilon_{\chi}=\varepsilon_{\operatorname{TV}}$, we have $\operatorname{TV}(q_{T}||p_{\textup{data}})\leq 2\varepsilon_{\operatorname{TV}}$.

We first define the bad sets where the error in the score estimate is large,

Given $t\geq 0$, let $t_{-}=h\left\lfloor\frac{t}{h}\right\rfloor$. Given a bad set $B$, define the interpolated process by

In other words, simulate the reverse SDE using the score estimate as long as the point is in the good set (for the current $p_{t}$) at the previous discretization step, and otherwise use the actual gradient $\nabla\ln p_{t}$. Let $\overline{q}_{t}$ denote the distribution of $\overline{z}_{t}$ when $\overline{z}_{0}\sim q_{0}$; note that $q_{nh}$ is the distribution resulting from running LMC with estimate $b$ for $n$ steps and step size $h$. Note that this process is defined only for purposes of analysis, as we do not have access to $\nabla\ln p_{t}$.

We can couple this process with the predictor algorithm using $s$ so that as long as $x_{mh}\not\in B_{mh}$, the processes agree, thus satisfying condition 1 of Theorem 4.1.

Let $T=Nh$, and let $K_{\chi}=\chi^{2}(q_{0}||p_{0})$. Then by Theorem 4.3,

For this to be bounded by $\varepsilon_{\chi}^{2}$, it suffices for the terms to be bounded by $\frac{\varepsilon_{\chi}^{2}}{2},\frac{\varepsilon_{\chi}^{2}}{4},\frac{\varepsilon_{\chi}^{2}}{4}$; this is implied by

(We choose $h$ so that the condition in Theorem 4.3 is satisfied; note $\varepsilon_{\chi}\leq 1$.) By Theorem 4.1,

In order for this to be $\leq\varepsilon_{\operatorname{TV}}$, it suffices for

Supposing that we run for time $T=\Theta(T_{\min})$, we have that $n=\frac{T}{h}=O\left({\frac{C_{T}T_{\min}}{h}}\right)$. Thus it suffices for