Score-based generative models have demonstrated significant practical success in data-generating tasks. The models establish a diffusion process that perturbs the ground truth data to Gaussian noise and then learn the reverse process to
transform noise into data. However, existing denoising methods such as Langevin dynamic and numerical stochastic differential equation solvers enjoy randomness but generate
data slowly with a large number of score function evaluations, and the ordinary differential equation solvers enjoy faster sampling speed but no randomness may influence the sample quality. To this end, motivated by the Stochastic Gradient Descent (SGD) optimization methods and the high connection between the model sampling process with the SGD,
we propose adaptive momentum sampling to accelerate the transforming process without introducing additional hyperparameters. Theoretically, we proved our method promises convergence under given conditions. In addition, we empirically show that our sampler can produce more faithful images/graphs in small sampling steps with 2 to 5 times speed up and obtain competitive scores compared to the baselines on image and graph generation tasks.

§ INTRODUCTION

Two successful classes of diffusion models involve sequentially corrupting training data with slowly increasing noise, and then learning to reverse this corruption to sample new data from a generative model. The Score Matching with Langevin Dynamics (SMLD) trained a network to predict the gradient of the data distribution, a technique known as score matching [9, 38]. Indeed, the network can be trained from data, and we can use the gradient (log probability density) obtained from data with Langevin Dynamic (LD) [29, 40] to generate new data samples from our learned data distribution.

The Denoising Diffusion Probabilistic Modeling (DDPM) [6, 30]) employs probability tools for training probabilistic models and reversing the noise corrupting and generating the images. The authors in [34] have provided a unified score-based Stochastic Differential Equation (SDE) framework for these two diffusion models. In this article, we refer to these two classes of generative models as Score-based Generative Models (SGMs). SGMs have shown superior performance in many generation tasks, e.g. image generation [6, 34], video generation [7], and graph generation [10, 24, 36].

In addition to their promising applications, SGMs mainly suffer from the expensive generation cost caused by the diffusion process. For high-quality sample synthesis, these models typically require over $1000$ sampling steps for high-quality images, resulting in a lower sampling speed compared to GANs [5]. Recently, most advancements in samplers have focused on accelerating the sampling process with the Ordinary Differential Equation (ODE) samplers [20, 31, 39] or necessitate additional model training targets [3]. However, the primary drawback of ODE samplers is the deterministic nature of the generation process. In other words, the generation process is deterministic as the initial point is determined. Stochastic sampling still presents a bottleneck for existing SGMs. It's worth noting that introducing randomness into the generation process can yield more faithful results and enable the computation of likelihood properties [3, 13]. Given that stochastic generation via Langevin dynamics bears a close relationship with Stochastic Gradient Descent (SGD)  [22, 34], this connection unveils new avenues for optimizing stochastic sampling techniques. Integrating techniques from SGD into the generation process of SGMs could potentially address this challenge.

SGD has emerged as an effective choice for addressing large-scale problems due to its low computational costs and high empirical performance [27], this method aims to minimize a given objective function through a gradual decline scheme. Polyak’s Heavy Ball (HB) method algorithm [25] extends the standard SGD method by incorporating a special term to capture the changing trend of variations, thereby guiding the next parameter to iterate toward the solution. This method is referred to as the Stochastic Heavy Ball (SHB) method, which introduces a new parameter $\beta$ to control the momentum term, necessitating additional fine-tuning. Previous research [35, 41] introduced an adaptive method for tuning $\beta$ based on historical update information, resulting in improved experimental performance and offering a means to reduce the overhead associated with additional hyperparameter tuning.

§.§ Contributions

This paper is built on the previous work by [33, 41]. After introducing all necessary backgrounds, we will propose theoretical improvements on sampling time and sampling quantification with our method. Our main contributions are as follows:

* We integrated the traditional SGD method into score-based generative models to develop a novel sampling method as an alternative to the Langevin sampling method. Opening a new frontier for gradient-based sampling optimization.

* We leveraged the mathematical tools in the Markov Chain to analyze the properties of the sampling process of SGMs and promise the convergence and speed of our proposed sampling method.

* We further evaluate our sampling method as an SDE solver, integrating this method into PC samplers [34] as a new corrector, performing various experiments showing that our proposed sampler can generate more faithful image/graph with faster sampling speed.

In addition, we empirically validate our proposed method's performance on some image/graph generalization tasks: CIFAR10 [16], CelabA [19], FFHQ [14], LSUN (church and bedroom) [44]. For the graph generation tasks, we choose Ego-small [24], Community-small [21], and a large graph datasets Grid. Showing that our approach can lead to significantly higher quality samples than all other methods tested. For reproducibility, we will make our code and the calculated samples publicly available. We present our experimental results in section <ref>.

§ BACKGROUND

§.§ Related Works

There are many accelerated samplers for existing SGMs. Paper [3] integrated velocity in Langevin dynamics and proposes the Critically-damped Langevin Dynamic (CLD), but tailored with a modified training target and additional hyperparameters. Additionally, most existing accelerated samplers focus on ODE optimization thus there are no additional noise terms in deterministic sampling, and a larger step size can be adopted [13, 20, 31, 39, 42]. The authors in [31] introduced Denoising Diffusion Implicit Models (DDIM), a non-Markov sampler derived from the Euler sampler. Paper [39] utilized adaptive momentum by root mean square propagation for DDIM, but the momentum parameter $\beta$ still requires fine-tuning. Due to the deterministic nature of the sampling process (non-Markovian), there is a lack of diversity in the generated images as no noise is introduced. Conversely, LD incorporates randomness and has a physical explanation. To the best of our knowledge, there are no training-free and hyperparameter-free accelerated Langevin samplers. In this paper, we mainly focus on stochastic samplers.

§.§ Score Matching with Langevin Dynamic

For a given continuously differentiable probability density $p({\mathbf{x}})$, we denote $\nabla_{\mathbf{x}} \text{log}p(\mathbf{x}) $ as its score function. Once the score function is known, we can utilize the Langevin dynamic to sample new data from an initialized sample by the following scheme

\begin{equation}\label{eq::langevin_sample}
\mathbf{x}^{t+1} = \mathbf{x}^{t} +\alpha \nabla \text{log} p(\mathbf{x}^{t}) + \sqrt{2\alpha} \varepsilon(\mathbf{x}^{t}),\ \ \ 1\le t\le T,
\end{equation}

where $\varepsilon(\mathbf{x}^{t})$ is a noise term dependent or independent on $\mathbf{x}^{t}$, $\alpha$ is a step size parameter chosen to be sufficiently small, $T$ is chosen to be sufficiently large. In the annealed Langevin dynamic [32], the noise term is a zero-centered Gaussian noise with standard deviation $\mathbf{I}$, which can be represented by $\varepsilon$. The convergence of $\mathbf{x}^{t}$ to $p(\mathbf{x})$ is guranteed under certain regularity conditions [29, 40].

Denoising Score Matching (DSM) [8] with Langevin dynamic involves estimating the score function for the objective data distribution $p(\mathbf{x})$. This can be achieved by training a Noise Conditional Score Network (NCSN) [33], denoted as $s_{\theta}(\mathbf{x},\sigma)$ by the following loss function [38]

\begin{equation}\label{eq::scorenetwork}
\frac{1}{2n}\sum_{i=1}^{n}\mathbb{E}_{p_{data}(\mathbf{x})}\mathbb{E}_{p_{\sigma_i}(\tilde{\mathbf{x}}|\mathbf{x})}\Big[  \big\| \sigma_i s_{\theta}(\tilde{\mathbf{x}},\sigma_i) +\frac{\tilde{\mathbf{x}}-\mathbf{x}}{\sigma_i}\big\|^2\Big].
\end{equation}

Here, $p_{data(\mathbf{x})}$ represents the underlying true data distribution of our target dataset, $ \{\sigma _i\}_{i=1}^{n}$ are pre-set noise scale satisfying $\sigma_n<\sigma_{n-1}<\cdots<\sigma_1$, and $p_{\sigma_i}(\tilde{\mathbf{x}}|\mathbf{x})$ represents the corrupted data distribution. In practice, we use the perturbation kernel $p_{\sigma}(\tilde{\mathbf{x}}|\mathbf{x}):=\mathcal{N}(\tilde{\mathbf{x}};\mathbf{x},\sigma^2\mathbf{I})$.

The authors in [34] established a connection between the forward process and the Stochastic Differential Equation (SDE), solving the reverse SDE using a reverse diffusion predictor. For the forward process in SMLD, we have
\begin{equation}\nonumber
\mathbf{x}^{t+1} = \mathbf{x}^{t} + \sqrt{\sigma_{t+1}^2-\sigma_{t}^2}\mathbf{n}^{t},
\end{equation}

where $\mathbf{n}^t \sim \mathcal{N}(0,\mathbf{I})$, this process can be formulated by an SDE with the form
$$d\mathbf{x} = f(\mathbf{x,t})dt+g(t)d\mathbf{w},$$

where $\mathbf{w}$ is the standard Wiener process, $f(\cdot, t): \mathbb R^d \rightarrow \mathbb R^d$ is a function dipict the diffusion process of $\mathbf{x}^{t}$, and $g(\cdot):\mathbb R\rightarrow R$ is a scalar function depend on the diffusion time $t$.

Besides sampling new data by Eq. (<ref>), paper [33] suggests integrating all the noise information to establish a step size sequence $\{\alpha_i \}_{i=1}^{n}$ with a schedule $\alpha_i = \alpha \sigma_{i}^2/\sigma_{n}^2$, which is called Annealed Langevin dynamics (ALS). We provide its pseudo-code in Algorithm <ref>.

§.§ Intergrate Momentum in SMLD

Annealed Langevin Sampling

hyperparameters $\epsilon>0$, score-net $s_{\theta}$, noise scale $(\sigma_i)_{i=1}^{n}$, $n_{\sigma}$

Initialize $\mathbf{x}$
$i \leftarrow 1$ to $n$
$\alpha_{i} = \epsilon \sigma_{i}^{2}/{\sigma_{n}^{2}}$
$k \leftarrow 1$ to $n_{\sigma}$
Draw Gaussian noise $\mathbf{\varepsilon} \sim \mathcal{N}(0,\mathbf{I})$
$\mathbf{x} \leftarrow \mathbf{x} +\alpha_i s_{\theta}(\mathbf{x},\sigma_i) +\sqrt{2\alpha_i}\varepsilon$

return $\mathbf{x}+\sigma_{n}^{2}s_{\theta}(\mathbf{x},\sigma_n)$
return $\mathbf{x}$

Adaptive Momentum Sampling (ours)

hyperparameters $\epsilon>0$, $\delta>0$, score-net $s_{\theta}$, noise scale $(\sigma_i)_{i=1}^{n}$, $n_{\sigma}$

Initialize $\mathbf{x}, \beta = 0$
$i \leftarrow 1$ to $n$
$\alpha_{i} = \epsilon \sigma_{i}^{2}/{\sigma_{n}^{2}}, \tilde{\alpha} = \alpha_{i} (1 + \beta)^2$
$k \leftarrow 1$ to $n_{\sigma}$
Draw Gaussian noise $\varepsilon \sim \mathcal{N}(0,\mathbf{I})$
Calculate $\beta$ and Update $\mathbf{x}$ by Eq. (<ref>), and (<ref>)
return $\mathbf{x}+\sigma_{n}^{2}s_{\theta}(\mathbf{x},\sigma_n)$
return $\mathbf{x}$

For SMLD, samples are generated by the reverse progress using Langevin dynamics in Algorithm <ref>. It is worth mentioning that this updating schedule can be regarded as a gradient ascent of the function $\nabla_{\mathbf{x}}\text{log}p(\mathbf{x})$.  This forms the primary connection between Langevin dynamics and Stochastic Gradient Descent (SGD). The key difference between Langevin dynamics and traditional Stochastic Gradient Descent (SGD) lies in whether the sampling process is stochastic or deterministic.

Consider the following Stochastic Gradient Descent (SGD)

\begin{equation}\label{eq::SGD}
\mathbf{x}^{t+1} = \mathbf{x}^{t} -\alpha\nabla f(\mathbf{x}^t),
\end{equation}

where $\alpha$ is the step size hyperparameter. We can regard the Langevin dynamic as a special form of SGD. Therefore, we can incorporate some improved techniques of SGD into the Langevin dynamic and develop a new sampling schedule for SGMs.

HB method introduced a momentum term for SGD Eq. (<ref>) and updates $\mathbf{x}$ with momentum and gradient

\begin{equation}\label{SHB}
\mathbf{x}^{t+1}= \mathbf{x}^{t} - \alpha\nabla f(\mathbf{x}^{t}) +\beta(\mathbf{x}^{t}-\mathbf{x}^{t-1}),
\end{equation}

where $\beta \in [0,1)$ is the momentum hyperameter and requires additional adjustment. If we add an additional normalization term to Eq. (<ref>), we derive the form of Normalized Stochastic Heavy Ball (NSHB) which integrates by the following

\begin{equation}\label{eq::NSHBupd}
{\mathbf{x}}^{t+1} ={\mathbf{x}}^{t}- \alpha (1-\beta)\nabla f(\mathbf{x}^{t}) +\beta ({\mathbf{x}}^{k} -{\mathbf{x}}^{k-1}).
\end{equation}

Moreover, we slightly modified the stepsize schedule in Algorithm <ref> with $\tilde{\alpha} = \alpha (1 + \beta)^2$. Now we use our trained score network from in Eq. (<ref>) and change our goal to maximize the objective function $\text{log} p_{\mathbf{x}}$ for stochastic gradient ascent, we can rewrite NSHB and formulate a new sampling method by replacing the gradient with the trained score network $s_{\theta} (\mathbf{x}^{t},\sigma)$ as follows

\begin{equation}\label{eq::NSHBsamp}
\mathbf{x}^{t+1} ={\mathbf{x}}^{t} + \tilde{\alpha} \mathbf{m}^{t} + \sqrt{2\alpha} \varepsilon,\
\mathbf{m}^{t} = \beta _t \mathbf{m}^{t-1} + (1-\beta_t) s_{\theta} (\mathbf{x}^{t},\sigma).
\end{equation}

§ TECHNICAL TOOLS

§.§ Notation

To start with our analysis, we introduce several notations to enhance the readability. Consider the finite-dimensional vector space $\mathbb R^d$ embedded with the canonical inner product  $<\cdot,\cdot>$. We use lowercase and lowercase boldface letters to separately denote scalars and vectors like  $\mathbf{x}=(x_1,x_2,\cdots,x_d)$, $\|\mathbf{x}\|$ denote its $\ell_2$ norm. The symbol $\otimes$ denotes the tensor product.  Uppercase boldface letter to denote a matrix, i.e., $\mathbf{A}$, and we use $\|\mathbf{A}\|$ to denote its $\ell_2$ norm, $\mathbf{I}$ denotes identity matrix. For $a =\mathcal{O}(b)$, we specify there exists a positive constant $0<C<+\infty$, satisfies that $a<Cb.\ \mathcal{B}(\mathbb{R}^d)$ stands for the Borel-$\sigma$ field of $\mathbb{R}^d,\ \delta_{\mathbf{x}}$ stands for the Diarc measure for $\mathbf{x}$. We use $\mathcal{P}(\mathbb{R}^d)$ to represent the space consisting all Borel probability measure $\theta$ in $\mathbb{R}^d$ which satisfy $\int_{{\mathbb{R}^d}} \|\mathbf{z}\|^2\theta(d\mathbf{z})<+\infty$. For the derivative operation and Hessian matrix for a function $f({\mathbf{x}}): \ \mathbb{R}^d \rightarrow \mathbb{R}$ by $ \nabla f({\mathbf{x}})$ and $\nabla^2 f({\mathbf{x}})$,   respectively.

For probably measures $\theta,\hat{\theta}$ in $\mathcal{P}(\mathbb{R}^d)$, define the Wasserstein distance of order 2
between $\theta$ and $\hat{\theta}$ as

\begin{equation}\nonumber
W_2(\theta,\hat{\theta}) :=\underset{\xi\in\prod(\theta,\hat{\theta})}{\inf} (\int\|\mathbf{x}-\mathbf{y}\|^2 \xi(d\mathbf{x},d\mathbf{y}))^{\frac{1}{2}},
\end{equation}

where $\prod(\theta,\hat{\theta})$ is a set of all joint probability distribution between $\theta$ and $\hat{\theta}$.

Now we collect some assumptions underlying our analysis, let $f : \mathbb R^{d} \rightarrow \mathbb R$ satisfies the following:

Assumption 1. The function $f$ is strong convex with a constant $\mu >0$, i.e., for all       $\mathbf{x},\mathbf{y} \in \mathbb R^{d}$, and $t\in[0,1) $ such that

\begin{equation}\nonumber
f(t\mathbf{x}+(1-t)\mathbf{y}) \le tf(\mathbf{x}) + (1-t)f(\mathbf{y}) - \frac{\mu}{2}t(1-t)\|\mathbf{x-y}\|^2.
\end{equation}

Assumption 2. The function $f$ is differentiable, and the gradient is $L$-Lipschiz, i.e., for all $\mathbf{x,y}\in\mathbb R^d$

\begin{equation}\nonumber
\|\nabla f(\mathbf{x})-\nabla f(\mathbf{y})\| \le L\|\mathbf{x-y}\|.
\end{equation}

§.§ Markov Chain Associated with Momentum

In this section, we formulate NSHB Eq. (<ref>) within the framework of Markov chains and introduce fundamental notions related to the theory [1, 23], we define $\mathbf{z}^{t+1} \in \mathbb{R}^{2d}= \begin{bmatrix}
\mathbf{x}^{t+1} \\
\mathbf{x}^{t}
\end{bmatrix},$ it is evident that $(\mathbf{z})^t$ constitutes a Markov chain,  from which we derive the subsequent Markov kernel $R_{\alpha}$ on $(\mathbb R^{2d},\mathcal{B}(\mathbb R^{2d}))$ associated with NSHB in Eq. (<ref>) iterates. We observe that $ R_{\alpha}(\mathbf{z}^{t},C) = P(\mathbf{z}^{t+1}\in C|\mathbf{z}^{t}).$

For all interger $t$, we define the Markov kernel $R_{alpha}^{t}$ recursively by $R_{\alpha}^{1} = R_{\alpha}$. For $t\ge 1$, $\hat{\mathbf{v}}\in \mathbb R^{2d}$, and $C\in \mathcal{B}(\mathbb R^{2d})$,

\begin{equation}\nonumber
R_{\alpha}^{t+1} (\hat{\mathbf{v}},C):= \underset{\mathbb R^{2d}}{\int}R_{\alpha}^{t}(\hat{\mathbf{v}},d\mathbf{v})R_{\alpha}(\mathbf{v},C).
\end{equation}

For any probability measure $\lambda$ on $(\mathbb R^{2d},\mathcal{B}(\mathbb{R}^{2d}))$, we define the composite measure by [37]

\begin{equation}\nonumber
\lambda R_{\alpha}^{t}(C) = \underset{\mathbb{R}^{2d}}{\int}\lambda (d\mathbf{v}) R_{\alpha}^{t}(\mathbf{v},C).
\end{equation}

There exist numerous technical tools for establishing a schedule for tuning the step size $\alpha_i$. Similar to [32], we are motivated by maintaining the magnitude of the "signal-to-noise" ratio. Introducing additional momentum can propagate noise from the previous step. Commencing with the initial step, our "signal-to-noise" ratio is given by$\alpha_i s_{\theta}(\mathbf{x}^i,\sigma_i)/(1+\beta)\sqrt{2\alpha_i}$. If we have $\alpha_i \propto \sigma_i^2(1+\beta)^2$, we can infer that $\mathbb E(\mathbf{x}^i,\sigma_i)/(1+\beta)\sqrt{2\alpha_i}\propto \frac{1}{2}.$

§ MATHEMATICAL ANALYSIS AND ALGORITHM

§.§ Quadratic Case

In this section, we investigate the inefficiency property when $f$ takes the form of a quadratic function:
$$f(\mathbf{x})=\frac{1}{2}\mathbf{xAx}^T +\mathbf{b}^T\mathbf{x}+c,\ \nabla f(\mathbf{x}) = \mathbf{Ax}+\mathbf{b}.$$
Furthermore, we denote the smallest and largest eigenvalues of $\mathbf{A}$ as $\mu$ and $L$ respectively. Hence, we can reformulate Eq. (<ref>) with an additional noise term as

\begin{equation}\nonumber
\begin{bmatrix}
\mathbf{x}^{t+1} \\
\mathbf{x}^{t}
\end{bmatrix} = \begin{bmatrix}
(1+\beta)\mathbf{I}-\alpha(1-\beta)\mathbf{A} &-\beta \mathbf{I} \\
\mathbf{I} & \mathbf{0}
\end{bmatrix}\begin{bmatrix}
\mathbf{x}^{t} \\
\mathbf{x}^{t-1}
\end{bmatrix}+\begin{bmatrix}
\end{bmatrix} +\varepsilon_{t+1}\Big(\begin{bmatrix}
\mathbf{x}^{t} \\
\mathbf{x}^{t-1}
\end{bmatrix}\Big).}
\end{equation}

For simplicity, we denote
\begin{equation}\nonumber
\mathbf{T}:=\begin{bmatrix}
(1+\beta)\mathbf{I}-\alpha(1-\beta)\mathbf{A} &-\beta \mathbf{I} \\
\mathbf{I} & \mathbf{0}
\end{bmatrix},
\end{equation}

we can have the following Markov chain

\begin{equation}\label{markovz}
\mathbf{z}^{t+1} = \mathbf{Tz}^{t} + \begin{bmatrix}
\end{bmatrix}+ \varepsilon_{t+1}(\mathbf{z}^{t+1}).
\end{equation}

For any $\alpha \in(0,\frac{2}{L}]$, the Markov chain $(\mathbf{z}^t)_{t\ge 0}$ defined by the recursion Eq. (<ref>) admits a unique stationary distribution $\pi_{\alpha}^*$ such that $\pi_{\alpha}^* R_{\alpha} = \pi_{\alpha}^*$. Additionally, for all $\mathbf{\hat{z}}\in\mathbb R^{2d},t\in \mathbb N^*$, we have

\begin{equation}\nonumber
W(R_{\alpha}^{t}(\mathbf{z},\cdot)) = (\rho(T)+\epsilon)^{t}\underset{\mathbb R^{2d}}{\int}(\|\mathbf{z-\hat{z}}\|^2\pi_{\alpha}^*(d\mathbf{z}))^{\frac{1}{2}}.
\end{equation}

We can regard $\hat{\mathbf{z}}$ in Proposition <ref> as the initialization of the sampling progress.  Thus it is straight forward to yield that $\pi_{\alpha}^{*} = \begin{bmatrix}
\pi_{\alpha}\\
\pi_{\alpha}
\end{bmatrix}.$

Next, we analyze the properties of the chain starting at $\mathbf{z}_{0}$, distributed according to $\pi_{\alpha}^*$. Consequently, we find that the mean of the distribution is given by $\overline{\mathbf{z}} = \int_{\mathbb R^{2d}} \mathbf{z} \pi_{\alpha}^*d\mathbf{z} = \mathbf{z}^{*} + \alpha\triangle +\mathcal{O}(\alpha^2) .$

Assume $\mathbf{A}$ is a positive defined matrix, $\alpha \in (0,\frac{L}{2})$. Then it holds that $\mathbf{I-T}\mathbf{T}^T$ is invertible and we have the following property of the expectation on the distribution $\pi_{\alpha}$

\begin{equation}\nonumber
\underset{\mathbb R^{2d}}{\int} (\mathbf{z}-\mathbf{z}^*)^{\otimes 2}\pi_{\alpha}^*d(\mathbf{z}) = (\mathbf{I}-\mathbf{T}\mathbf{T}^T)^{-1} \alpha\underset{\mathbb R^{2d}}{\int} \hat{\mathcal{C}}(\mathbf{z})\pi_{\alpha}^*d(\mathbf{z}).
\end{equation}

Where $\hat{\mathcal{C}}(\mathbf{z})$ is a diffierentiable function for variable $\mathbf{z}$. This proposition highlights the crucial fact that for a given quadratic function, the mean under the limit distribution is the optimal point.

§.§ Strong Convexity Case

In this section, we dive into analyzing the convergence of Eq. (<ref>) in the strong convex case. Assume assumptions 1 and 2 holds, we have the following proposition
Assume $\alpha \in (0,2/L)$, then $\mathbf{I} \otimes \nabla f(\mathbf{x}^{*})+\nabla f(\mathbf{x}^{*})\otimes\mathbf{I}-\alpha(1-\beta)^2\nabla f(\mathbf{x}^{*})^{\otimes 2}$ is invertible, we derive

\begin{equation}\nonumber
\begin{aligned}
&\underset{\mathbb R^d}{\int} (\mathbf{x}-\mathbf{x}^*)^{\otimes 2}\pi_{\alpha}d(\mathbf{x}) = \alpha(1-\beta) \mathbf{B}\mathcal{C}(\mathbf{x}^{*}) +\mathcal{O}(\alpha^2),\\
\overline{\mathbf{x}}&-\mathbf{x}^{*} =\alpha(1-\beta)\big(\nabla^2 f(\mathbf{x}^{*})\big)^{-1}\nabla^3 f(\mathbf{x}^{*})\mathbf{B}\mathcal{C}(\mathbf{x}^{*})+\mathcal{O}(\alpha^2),
\end{aligned}
\end{equation}

where $\mathcal{C}(\mathbf{x})$ is a differentiable function devoted to variable $\mathbf{x},$ and
$$\mathbf{B} = (\mathbf{I} \otimes \nabla f(\mathbf{x}^{*})+\nabla f(\mathbf{x}^{*})\otimes\mathbf{I}-\alpha(1-\beta)^2\nabla f(\mathbf{x}^{*})^{\otimes 2})^{-1}.$$

This proposition shows that in the case of strong convexity, i.e., when the function is non-quadratic and satisfies $\nabla^3 f(\mathbf{x})\ne 0$, there will be an additional error compared to the quadratic case.

§.§ Our Algorithm

It is straightforward that the hyperparameter $\beta$ can significantly influence the success of the sampling process and may necessitate additional parameter tuning. To address this, paper [41] established a scheme to approximate the optimal choice of $\beta$ and thereby reduce the parameter tuning cost by the following:

\begin{equation}\label{eq::The_best_beta}
\beta_{t+1} =\left\{
\begin{array}{lr}
{\rm Proj}_{[0,1-\delta]} \Bigg(\dfrac{1-\alpha\dfrac{\|{\mathbf{g}}^t-{\mathbf{g}}^{t-1}\|}{\|{\mathbf{x}}^t-{\mathbf{x}}^{k-1}\|}}{1+\alpha\dfrac{\|{\mathbf{g}}^t-{\mathbf{g}}^{t-1}\|}{\|{\mathbf{x}}^t-{\mathbf{x}}^{t-1}\|}}\Bigg)^2  &t\ge 2,\\
\end{array}
\right.
\end{equation}

Where ${\rm Proj}_{[0,1-\delta]}(\cdot):=\max(0, \min(\cdot, 1-\delta))$ operation with a threshold $\delta$. Putting it all together, we have developed a new sampling method suitable for SGMs. This method ensures convergence under specified assumptions and exhibits favorable stochastic properties. We present its pseudo-code in  Algorithm <ref>, and name it Adaptive Momentum Sampling (AMS).

[EM-LC sampler, 150 NFE]
[EM-MC sampler, 150 NFE]
[RD-MC sampler, 150 NFE]
Partial estimate of FID (2k samples) with different initial $\epsilon$. We can see that if an appropriate $\epsilon$ is chosen, MC can bring a considerable improvement in image quality.

§ NUMERICAL RESULTS

To evaluate the performance of our approach, we focus on image synthesis and implementation tasks on the NCSN architecture [32]. All experiments are conducted on a server with 8 NVIDIA 3090 GPUs and a host with 1 NVIDIA 4090 GPU. To aid reproducibility, we detail all experimental hyperparameters in Appendix <ref>.

§.§ Experimental Setup

Datasets and Models. For the image generation tasks, we evaluate our proposed algorithm (<ref>)  within the score-net architecture proposed by [34]. Our primary focus for image experiments is on the CIFAR-10 dataset [16]. Additionally, we assess our sampler on other image datasets; further details and additional samples can be found in the App. <ref>. For graph generation tasks, we utilize the score-net provided by [10], employing datasets such as Ego-small [24], Community-small [21] and the Grid dataset.

Baselines. When utilizing the score-net of NCSN2 and UNET, we choose the baselines ALS [32] and CLS [12] as proposed in their original work. We present the experimental results in the App. <ref>. For further comparison in NCSN++ cont./DDPM ++ cont., we include different SDE solvers such as the Reverse Diffusion (RD) predictor solver, Euler-Maruyama (EM) solver [15], each combined with our momentum sampler and the original Langevin sampler, DDIM [31], and Gotta Go Fast [11]. Moreover, we rewrite the Critically-Damped Langevin sampler (CDL) in [3] with a corrector form to enhance the extensiveness. The primary evaluation metric used is the Fréchet Inception Distance (FID) [5]. For graph quality comparison, we include various baselines, such as autoregression models like GraphRNN [43], DeepGME [17], GRAN [18], and GraphDF [21]. Diffusion-based methods, e.g. EDP-GNN [24], GDSS [10], and DiGress [36].

§.§ Image Generation Tasks

FIDs for different NFEs on CIFAR10 (VE). In those small NFEs, only DDIM outperforms our sampler.

Sampler/NFE     150      195      210     270     300

RD-MC (ours)     $6.75 $     $ \textbf{3.05}$    $ \textbf{2.86}$     $2.85$    $2.93$

RD-CDL    $ 38.6 $     $ 5.04$    $ 4.37$   $ 4.80$     $ 4.43$

150     270     500     1000     2000

RD-LC     $47.6$    $8.87$     $3.62$     $2.62$    $\textbf{2.38}$

RD     $40.3$     $7.51$    $3.72$    $2.98$    $2.89$

EM     $13.6$     $12.1$    $12.3$    $12.6$    $12.9$

Gotta Go Fast     $8.85$     $3.23$     $2.87$     -    -

EM (VP)     $30.3$     $13.1$    $ 4.42$     $ \textbf{2.46}$     $2.43$

DDIM (VP)     $\textbf{4.03}$     $ 3.53$     $3.26$    $3.09$     $3.01$

MC (ours)    $235$     $27.6$   $9.16$     $17.7$   $9.29$

CDL     $273$    $27.9$   $23.7$     $32.4$    -

LC    $286$     $202$    $121$    $57.8$     $25.9$

1 We present the corresponding IS and KID in the Appendix.

[ Bedroom, 500 NFE]

[ CelebA, 900 NFE]

In this section, we assess our momentum corrector combined with different predictors, such as the RD predictor and EM predictor, focusing on image synthesis using PC samplers in NCSN2 [32], UNET [12], NCSN++, and DDPM++ [34].

For the NCSN++ experiments, we ran the RD-MC sampler with a total NFE of $150$ for NCSN++ cont. (VE) [34], increase the hyperparameter $n_\sigma$ from $0$ (no corrector) to $3$. Our findings indicate that we can improve the FID score by incorporating correct steps for both the EM predictor and RD predictor. Specifically, choosing $n_\sigma = 2$ yields the best performance for both the EM and RD predictors. Conversely, adding additional LC correct steps yields contrasting results. For more details, please refer to Figure <ref>.

For our approach, we initially tested the best FID performance using $n_\sigma \in {1, 2, 3}$ with the RD predictor, and EM predictor. For these samplers without predictors, we adopted $n_\sigma = 5$ and tested the scores with different initial step size hyperparameters $\epsilon$. For predictors that include LC as a corrector, we used the suggested hyperparameter $n_\sigma = 1$ and chose the best $\epsilon$. Surprisingly, we observed that compared to predictors without a corrector, adding LC correction steps may worsen the FID in small NFEs.  However, [34] reported the best performance if NFE $= 2000$ with this PC sampler. Additionally, we tested different combinations of hyperparameters for this PC sampler,  details can be found in Appendix <ref>.

We observe that our proposed corrector MC can serve as an effective corrector for the RD predictor, as shown in Table <ref>, where it achieved a lower FID compared to the EM/RD sampler without a corrector. Moreover, the performance improvement for the MC corrector is notably higher than for the CDL corrector.  Specifically, for 150 NFEs, our MC substantially improves the FID performance for the RD predictor, from $40.3$ to $6.75$. In contrast, adding LC correction steps during the RD prediction times can have a negative effect on the FID in small NFEs, leading to an increase from $40.3$ to $47.6$. RD alone struggles to generate high-quality images if the NFE is $150$. However, when combined with our corrector, the RD predictor can produce images with a much lower FID even with a much smaller NFE.

For those small NFEs, e.g. $150$, $195$, our MC combined with RD predictor report an FID of $6.75, 2.86$, only DDIM outperforms RD-MC in $150$ NFE. We should figure out that, the DDIM sampler was designed for small NFE sampling, and can be seen as a method that dropped the Markov property to accelerate the denoising process derived from an ODE solver. Empirically, the EM predictor as an SDE solver performs better in the VP process than VE, the performance of this predictor is largely unchanged with different NFEs.

We found the optimal acceleration interval for RD-MC is less than $210$ with an FID of $2.86$, representing about $4$ to $5$ times speed-up compared to the PC sampler, and approximately $2$ times speed-up compared to the accelerated SDE solver Gotta Go Fast [11] (report an FID of $2.87$ with $490$ NFE). As we increase the NFE,  there is less performance gain and the FID score may worsen. We argue there may be an overfitting, which could occur in all gradient-based methods. Since our method adopts $n_\sigma=2$ throughout the whole experiment,  where every noise scale corresponds to a target function, only two momentum steps may lead to excessive acceleration if we use frequent predictor steps (large NFEs). Additionally, in the realm of non-predictor samplers, our method can still maintain a speed-up ratio of at least $4$ compared to LC.

Furthermore, we validated our MC sampler on the CIFAR10 dataset for Inception Score (IS) and Kernel Inception Distance (KID) using the same settings. We observed that using MC can also generate images with better performance in terms of IS and KID as an accelerator. The results are presented in Table <ref> and <ref>.

§.§ Graph Generation Tasks

Quantitative results on generic graph datasets. We can see that our method can generate higher-quality graphs or graphs of the same level as the testing methods under a small NFE.

2*Dataset         4cEgo-small     4cCommunity-small     4cGrid

4cReal     4cSynthetic     4cSynthetic

Methods     Deg. $\downarrow$     Clus. $\downarrow$     Orbit $\downarrow$    Avg.     Deg. $\downarrow$    Clus. $\downarrow$    Orbit $\downarrow$     Avg.     Deg. $\downarrow$    Clus. $\downarrow$    Orbit $\downarrow$     Avg.

4*Autoreg.     GraphRNN     $0.060$     $ 0.220$     $ 0.003$     $0.104$     $ 0.080$     $0.120$     $0.040$     $0.080$     $\textbf{0.064}$     $0.043$     $\textbf{0.021}$     $\textbf{0.043}$

DeepGMG     $0.040$     $ 0.100$     $ 0.020$     $0.053$     $ 0.220$     $0.950$     $0.400$     $0.523$    $ -$    $ -$    $ -$    $ -$

GRAN     $0.020$     $ 0.126$     $ \underline{0.010}$     $0.052$     $ 0.070$     $0.045$     $0.021$     $0.045$    $ -$    $ -$    $ -$    $ -$

GraphDF     $0.040$     $ 0.130$     $\underline{0.010}$     $0.060$     $ 0.060$     $0.120$     $0.030$     $0.070$    $ -$    $ -$    $ -$    $ -$

6*Diffusion     EDP-GNN (1000 NFE)    $0.052$     $ 0.093$     $ \textbf{0.007}$     $0.051$     $ 0.053$     $0.144$     $0.026$     $0.074$    $0.455$     $0.238$     $0.328$    $0.340$

GDSS (1000+ NFE)     $\textbf{0.021}$     $ 0.024$     $ \underline{0.007}$     $\textbf{0.017}$     $ 0.045$     $0.086$     $0.007$     $0.046$    $0.378$    $\underline{0.010}$     $0.442$    $0.276$

GDSS-seq (1000+ NFE)     $0.032$     $ 0.027$     $ 0.011$     $0.023$     $ 0.090$     $0.123$     $0.007$     $0.073$     $\underline{0.111}$    $\textbf{0.005}$     $\underline {0.070}$    $\underline{0.062}$

DiGress (500 NFE)    $ 0.120$     $0.180$     $0.034$    $ 0.111$     $0.089$     $\textbf{0.019}$     $0.748$     $0.285 $    $0.957$    $0.026$     $1.030$    $0.671$

MC (ours, 300 NFE)     $0.024$    $\textbf{0.009}$    $0.013$   $\textbf{0.015}$     $0.067$   $0.109$    $0.031$     $0.069$    $ -$    $ -$    $ -$    $ -$

MC (ours, 750 NFE)     $0.024$    $\underline{0.012}$    $0.014$     $0.020$     $\textbf{0.030}$     $\underline{0.042}$     $\textbf{0.003}$     $\textbf{0.025}$     $ 0.127 $     $ 0.016 $    $ 0.203 $     $ 0.115 $

1 The best and second results for each metric are separately bolded and underlined.

In this section, we extend our Algorithm <ref> to the filed or graph generation. Following the previous diffusion-based graph generation works [10, 36], we evaluate the quality of the generated graphs using maximum mean discrepancy (MMD) metrics. Specifically, we compare the distributions of degree, clustering coefficient, and the number of occurrences of orbits with 4 nodes across different samplers. For the hyperparameters, we search the sampling $\epsilon$ of our method in $\{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0\}$ and choose the result with the lowest average metric. For the $n_{\sigma}$, we follow the previous analysis and choose $n_{\sigma} = 2$ during all the experiments.

We present the experimental results in Table <ref>. It's important to note that all baseline diffusion methods, except DiGress, use an NFE of more than 1000, typically employing an NFE of 2000. We evaluated our method with NFEs of 300 and 750. In the Ego-small dataset, an NFE of 300 is sufficient for our momentum-corrector to achieve state-of-the-art performance, resulting in a speed-up of more than 3 times compared to the baseline diffusion samplers. Given that the graph scale of the Community-small dataset is larger than that of the Ego-small dataset, an NFE of $750$ is adequate for our algorithm to generate high-quality samples, resulting in a speed-up of approximately 2 times.

Quantitative results on Community-small. See that our MC does not incur additional computational time overhead but can improve quality.

NFE     Methods     Deg.     Clus.     Orbit      Time/s

4*300     EM     $0.074$    $0.111$    $\underline{0.007}$    $2.99$

EM-LD    $ \underline{0.061}$    $ 0.122$    $ 0.031$     $ 2.86$

EM-CLD     $0.065$    $0.122$   $0.031$    $2.94$

EM-MC (ours)    $0.067$   $\underline{0.109}$    $0.031$     $2.83$

4*450     EM     $0.049$    $0.121$    $0.013$    $4.28$

EM-LD     $0.042$    $0.071$    $0.005$    $4.09$

EM-CLD     $0.043$    $0.068$   $0.005$    $4.27$

EM-MC (ours)     $\underline{0.041}$    $\underline{0.053}$    $\underline{0.004}$    $4.10$

4*600     EM     $0.066$    $0.173$    $0.011$    $5.19$

EM-LD     $0.025$    $0.124$   $\underline{0.004}$     $5.21$

EM-CLD     $0.046$    $0.092$   $0.006$    $5.40$

EM-MC (ours)    $\underline{0.037}$    $\underline{0.048}$    $\underline{0.004}$    $5.26$

4*750     EM     $0.451$    $0.123$    $0.006$    $6.63$

EM-LD     $0.050$    $0.065$   $0.004$    $6.41$

EM-CLD    $0.047$    $0.072$   $\underline{\textbf{0.003}}$    $7.03$

EM-MC (ours)   $\underline{\textbf{0.030}}$     $\underline{\textbf{0.042}}$     $\underline{\textbf{0.003}}$     $6.88$

In addition to the datasets with small $|V|$, we also applied our method on the dataset Grid with $100\le|V|\le 400$. Upon examining the experimental results, we see that only GDSS-seq [10] and GraphRNN [43] outperform our method in the quantitative results, while the other baselines perform worse (average MMD larger than 0.2). Only our method with $750$ NFE achieves comparable results with GDSS-seq, with a speedup of about $2$ times. It's important to note that GDSS-seq generates node features and the weighted adjacency matrix sequentially, potentially doubling the generation consumption. By combining all the results from the above datasets, we can conclude that our method is a more general and efficient choice for graph generation tasks.

Besides directly comparing our algorithm to existing methods, we conducted an additional experiment to compare our method with different SDE correctors, such as LD, CLD, and no-corrector, across different NFEs. We used the same metrics and measured the sampling time for generating 128 graphs using a single NVIDIA 4090 GPU. The experimental results are shown in Table <ref>. By comparing the sampling times, we note that the additional $\beta$ computing step does not significantly impact the computation speed. As the NFE increases, our MC shows superior performance compared to all the tested Markov samplers in all three metrics while maintaining a sampling time close to that of the EM-LD sampler.

§ CONCLUSION

Motivated by traditional SGD optimization, we developed a novel gradient-based sampling method with the Markov technique for existing SGMs. This sampler requires no information beyond the score function with a free application.  Theoretically, we proved the convergence properties and speed of our approach in given conditions with the Markov Chain, thereby establishing a connection between SGD and gradient-based sampling and opening a new frontier for acceleration.

In addition, we test our sampler in NCSN (VE), empirically shows that our proposed approach performs better than the tested baseline samplers in small NFEs without incurring additional computational costs. For image generation tasks, our main results were obtained from NCSN architecture with the CIFAR10 dataset [16], while our sampler can be extended to all gradient-based models and corresponding datasets. We present the corresponding samples in Appendix <ref>. Moreover, we extended our method to graph generation tasks, demonstrating that our sampler can produce more faithful graphs with a lower NFE.

[1]
E. Çinlar and E. ðCınlar.
Probability and stochastics, volume 261.
Springer, 2011.

[2]
A. Dieuleveut, A. Durmus, and F. R. Bach.
Bridging the gap between constant step size stochastic gradient descent and markov chains.
The Annals of Statistics, 2017.

[3]
T. Dockhorn, A. Vahdat, and K. Kreis.
Score-based generative modeling with critically-damped langevin diffusion.
In International Conference on Learning Representations (ICLR), 2022.

[4]
C. Godsil and G. Royle.
Matrix Theory, pages 163–192.
Springer New York, New York, NY, 2001.

[5]
M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler, and S. Hochreiter.
Gans trained by a two time-scale update rule converge to a local nash equilibrium.
Advances in neural information processing systems, 30, 2017.

[6]
J. Ho, A. Jain, and P. Abbeel.
Denoising diffusion probabilistic models.
Advances in neural information processing systems, 33:6840–6851, 2020.

[7]
J. Ho, T. Salimans, A. Gritsenko, W. Chan, M. Norouzi, and D. J. Fleet.
Video diffusion models.
ArXiv, abs/2204.03458, 2022.

[8]
A. Hyvärinen.
Estimation of non-normalized statistical models by score matching.
J. Mach. Learn. Res., 6:695–709, dec 2005.

[9]
A. Hyvärinen, J. Hurri, and P. O. Hoyer.
Estimation of Non-normalized Statistical Models, pages 419–426.
Springer London, London, 2009.

[10]
J. Jo, S. Lee, and S. J. Hwang.
Score-based generative modeling of graphs via the system of stochastic differential equations.
In International Conference on Machine Learning, 2022.

[11]
A. Jolicoeur-Martineau, K. Li, R. Piché-Taillefer, T. Kachman, and I. Mitliagkas.
Gotta go fast when generating data with score-based models.
arXiv preprint arXiv:2105.14080, 2021.

[12]
A. Jolicoeur-Martineau, R. Piché-Taillefer, I. Mitliagkas, and R. T. des Combes.
Adversarial score matching and improved sampling for image generation.
In International Conference on Learning Representations, 2021.

[13]
T. Karras, M. Aittala, T. Aila, and S. Laine.
Elucidating the design space of diffusion-based generative models.
ArXiv, abs/2206.00364, 2022.

[14]
T. Karras, S. Laine, and T. Aila.
A style-based generator architecture for generative adversarial networks.
In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4401–4410, 2019.

[15]
P. E. Kloeden, E. Platen, P. E. Kloeden, and E. Platen.
Stochastic differential equations.
Springer, 1992.

[16]
A. Krizhevsky, G. Hinton, et al.
Learning multiple layers of features from tiny images.

[17]
Y. Li, O. Vinyals, C. Dyer, R. Pascanu, and P. W. Battaglia.
Learning deep generative models of graphs.
ArXiv, abs/1803.03324, 2018.

[18]
R. Liao, Y. Li, Y. Song, S. Wang, C. Nash, W. L. Hamilton, D. K. Duvenaud, R. Urtasun, and R. S. Zemel.
Efficient graph generation with graph recurrent attention networks.
In Neural Information Processing Systems, 2019.

[19]
Z. Liu, P. Luo, X. Wang, and X. Tang.
Deep learning face attributes in the wild.
In Proceedings of the IEEE international conference on computer vision, pages 3730–3738, 2015.

[20]
C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, and J. Zhu.
Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps.
ArXiv, abs/2206.00927, 2022.

[21]
Y. Luo, K. Yan, and S. Ji.
Graphdf: A discrete flow model for molecular graph generation.
In International Conference on Machine Learning, 2021.

[22]
S. Mandt, M. D. Hoffman, and D. M. Blei.
Stochastic gradient descent as approximate bayesian inference.
ArXiv, abs/1704.04289, 2017.

[23]
S. P. Meyn and R. L. Tweedie.
Markov chains and stochastic stability.
Springer Science & Business Media, 2012.

[24]
C. Niu, Y. Song, J. Song, S. Zhao, A. Grover, and S. Ermon.
Permutation invariant graph generation via score-based generative modeling.
In International Conference on Artificial Intelligence and Statistics, 2020.

[25]
B. Polyak.
Some methods of speeding up the convergence of iteration methods.
USSR Computational Mathematics and Mathematical Physics, 4:1–17, 1964.

[26]
Y. Qian, Q. Cai, Y. Pan, Y. Li, T. Yao, Q. Sun, and T. Mei.
Boosting diffusion models with moving average sampling in frequency domain.
ArXiv, abs/2403.17870, 2024.

[27]
H. Robbins and S. Monro.
A stochastic approximation method.
Annals of Mathematical Statistics, 22:400–407, 1951.

[28]
H. E. Robbins.
An empirical bayes approach to statistics.
In Breakthroughs in Statistics: Foundations and basic theory, pages 388–394. Springer, 1992.

[29]
G. O. Roberts and R. L. Tweedie.
Exponential convergence of langevin distributions and their discrete approximations.
Bernoulli, pages 341–363, 1996.

[30]
J. Sohl-Dickstein, E. Weiss, N. Maheswaranathan, and S. Ganguli.
Deep unsupervised learning using nonequilibrium thermodynamics.
In International conference on machine learning, pages 2256–2265. PMLR, 2015.

[31]
J. Song, C. Meng, and S. Ermon.
Denoising diffusion implicit models.
arXiv preprint arXiv:2010.02502, 2020.

[32]
Y. Song and S. Ermon.
Generative modeling by estimating gradients of the data distribution.
In Neural Information Processing Systems, 2019.

[33]
Y. Song and S. Ermon.
Improved techniques for training score-based generative models.
Advances in neural information processing systems, 33:12438–12448, 2020.

[34]
Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole.
Score-based generative modeling through stochastic differential equations.
In International Conference on Learning Representations, 2021.

[35]
T. Sun, H. Ling, Z. Shi, D. Li, and B. Wang.
Training deep neural networks with adaptive momentum inspired by the quadratic optimization.
ArXiv, abs/2110.09057, 2021.

[36]
C. Vignac, I. Krawczuk, A. Siraudin, B. Wang, V. Cevher, and P. Frossard.
Digress: Discrete denoising diffusion for graph generation.
In The Eleventh International Conference on Learning Representations, 2023.

[37]
C. Villani et al.
Optimal transport: old and new, volume 338.
Springer, 2009.

[38]
P. Vincent.
A connection between score matching and denoising autoencoders.
Neural computation, 23(7):1661–1674, 2011.

[39]
X. Wang, A.-D. Dinh, D. Liu, and C. Xu.
Boosting diffusion models with an adaptive momentum sampler.
ArXiv, abs/2308.11941, 2023.

[40]
M. Welling and Y. W. Teh.
Bayesian learning via stochastic gradient langevin dynamics.
In Proceedings of the 28th international conference on machine learning (ICML-11), pages 681–688, 2011.

[41]
Z. Wen, X. Deng, T. Sun, and D. Li.
Normalized stochastic heavy ball with adaptive momentum.
In European Conference on Artificial Intelligence, volume 372, pages 2615 – 2622, 2023.

[42]
S. Wizadwongsa, W. Chinchuthakun, P. Khungurn, A. Raj, and S. Suwajanakorn.
Diffusion sampling with momentum for mitigating divergence artifacts.
ArXiv, abs/2307.11118, 2023.

[43]
J. You, R. Ying, X. Ren, W. L. Hamilton, and J. Leskovec.
Graphrnn: Generating realistic graphs with deep auto-regressive models.
In International Conference on Machine Learning, 2018.

[44]
F. Yu, Y. Zhang, S. Song, A. Seff, and J. Xiao.
Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop.
arXiv preprint arXiv:1506.03365, 2015.

§ APPENDIX / SUPPLEMENTAL MATERIAL

§ LEMMAS AND PROOFS

We have that for any $\alpha \in (0,2/L)$, the integral $\int_{\mathbb R^d}\nabla f(\mathbf{x})\pi_{\alpha}d(\mathbf{x}) = 0$ [2].

Let $f(\mathbf{x})$ be a strongly convex function and twice-differentiable, and let ${\{\mathbf{x}^k}\}_{k\ge 0}$ generated by Eq. (<ref>). Assume that
$$0<\mu\le\min_\mathbf{x}\{\lambda_{\min}(\nabla^2 f(\mathbf{x}))\}\le \max_\mathbf{x}\{\lambda_{\max}(\nabla^2 f(\mathbf{x}))\}\le L$$
and Assumptions 1, 2, 3 hold. Then we have

\mathbb{E}\|{\mathbf{x}^K}-{\mathbf{x}^*}\|^2= (1-2\alpha\delta\mu)^K\mathbb{E} \|{\mathbf{x}^1-\mathbf{x}^*}\|^2 + \mathcal{O}(\alpha).

This Lemma indicates that $\int_{\mathbb R^d}\|\mathbf{x}-\mathbf{x}^{*}\|^2\pi_{\alpha}^*d(\mathbf{x}) = \mathcal{O}({\alpha})$ if we have $\mathbf{x}^{0}$ is distributed according to $\pi_{\alpha}^*$ and the Markov Kernel is stationary. Also, we have
$\int_{\mathbb R^d}\|\mathbf{x}-\mathbf{x}^{*}\|^3\pi_{\alpha}^*d(\mathbf{x}) = \mathcal{\mathcal{O}(\alpha}^{3/2})$, respectively.

Proof of proposition <ref>
For two distributions $\theta,\hat{\theta}\in\mathcal{P}_2(\mathbf{R}^{2d})$, we have the fact that there exist a couple of random variables $\mathbf{z}^{0}, \hat{\mathbf{z}}^{0} \ \text{satisfy}\ \mathcal{W}_{2}^{2}(\theta,\hat{\theta}) = \mathbf{E}[\|\mathbf{z}^{0}-\hat{\mathbf{z}}^0\|^2]$ [Theorem 4.1, <cit.>].

Asssume $\mathbf{z}^0, \hat{\mathbf{z}}^0$ are initializations and (${\mathbf{z}^{t})_{\ge 1}}, (\hat{\mathbf{z}}^t)_{\ge 1}$ are generated by Eq. (<ref>) for all $t\ge 1$, the distribution of $(\mathbf{z}^{t},\hat{\mathbf{z}}^t)$ belongs to $\prod (\theta R_{\alpha}^{t},\hat{\theta}R_{\alpha}^{t})$.

Note we assume $\mathbf{z}^0, \hat{\mathbf{z}}^0$ are independent of $\varepsilon_{1}$, we have the fact

\mathbb{E}[<\hat{\mathbf{z}}^0, \varepsilon(\mathbf{z}_0)>] =0.

Leverage the above we derive

\begin{equation}\nonumber
\begin{aligned}
& \mathcal{W}^{2}_{2}(\theta R_{\alpha}^{t},\hat{\theta}R_{\alpha}^{t}) \leq \mathbf{E}\|\mathbf{z}^t- \hat{\mathbf{z}}^t\|^2 \\
&= \mathbb{E} \|\mathbf{T}\mathbf{z}^{t-1}+\begin{bmatrix}
-\alpha(1-\beta)\mathbf{b} \\
\end{bmatrix} -\mathbf{T}\hat{z}^{t-1}-\begin{bmatrix}
\end{bmatrix}\|^2   \\
&=\mathbb{E}\|\mathbf{T}(\mathbf{z}^{t-1}-\hat{\mathbf{z}}^{t-1})\|^2 =\dots =\mathbb{E}\|\mathbf{T}^t(\mathbf{z}^0-\hat{\mathbf{z}}^0)\|^2 \\ & \leq\|\mathbf{T}^t\|^2\mathbb{E}\|\mathbf{z}^0,\hat{\mathbf{z}}^0\|^2 =\|\mathbf{T}^t\|\mathcal{W}^2(\theta,\hat{\theta}).
\end{aligned}
\end{equation}

Using Gelfland's theory [4], given any small $\epsilon >0$, there exists $K$ which is dependent of $T$ and $\epsilon$ that

\|\mathbf{T}^{t}\|^2 \leq (\lambda_{max}(\mathbf{T})+\epsilon)^{2t},

if we have $t\ge K$.
Hence, we are then led to

\begin{equation}\nonumber
\sum_{t=K}^{+\infty}\mathcal{W}(\theta R_{\alpha}^{t},\hat{\theta}R_{\alpha}^t)\leq (\sum_{t=K}^{+\infty}(\lambda_{max}(\mathbf{T})+\varepsilon)^t)\mathcal{W}(\theta,\hat{\theta})<+\infty.
\end{equation}

Which indicates

\begin{equation}\label{different_mu_equal_0}
\underset{t\rightarrow +\infty}{\lim}(\theta R_{\alpha}^t, \hat{\theta}R_{\alpha}^t) = 0,
\end{equation}
for any $\theta, \hat{\theta}$.
By setting $\hat{\theta} = \theta R_{\alpha}$, we derive $\underset{t\rightarrow +\infty}{\lim} \mathcal{W}(\theta R_{\alpha}^{t},\theta R_{\alpha}^{t+1})=0$ for any $\theta$. Note that $\mathcal{W}$ is a polish space [Theorem 6.16, [37]], we can see that sequence $(\theta R_{\alpha}^{t})_{t \ge 1}$ converges to some distribution $\pi_{\alpha}^{\theta}$, i.e.,

\begin{equation}\label{R_convergeto_pi}
\underset{t}{\lim}\mathcal{W}(\pi_{\alpha}^{\theta}, \theta R_{\alpha}) = 0.
\end{equation}

In the following, we prove $\pi_{\alpha}^{\theta}$ is not related to $\theta $ and unique. Assume that there exists a distribution $\pi_{\alpha}^{\hat{\theta}}$ satisfies that $\underset{t\rightarrow +\infty}{\lim} \mathcal{W}(\hat{\theta}R_{\alpha}^{t},\pi_{\alpha}^{\hat{\theta}}) = 0$.
Using the triangle inequality, we are led to

\begin{equation}\nonumber
\mathcal{W}(\pi_{\alpha}^{\theta},\pi_{\alpha}^{\hat{\theta}})\leq \mathcal{W}(\pi_{\alpha}^{\theta},\theta R_{\alpha}^{t})+\mathcal{W}(\theta R_{\alpha}^{t},\hat{\theta}R_{\alpha}^{t})+\mathcal{W}(\hat{\theta}R_{\alpha}^{t},\pi_{\alpha}^{\hat{\theta}}).
\end{equation}

Taking the limits as $t\rightarrow +\infty$ and use Eq. (<ref>), we conclude $ \mathcal{W}(\pi_{\alpha}^{\theta},\pi_{\alpha}^{\hat{\theta}}) = 0$. We see that $\pi_{\alpha}^{\theta} = \pi_{\alpha}^{\hat{\theta}}$. Therefore, we can denote $\pi_{\alpha}^{\theta}$  as $\pi_{\alpha}^{*}$. It is straightforward $\pi_{\alpha}^{*}$ is unique for given $\alpha$.

Moreover, we see that
$$\mathcal{W}(\pi_{\alpha}^{*},\pi_{\alpha}^{*}R_{\alpha})\leq\mathcal{W}(\pi_{\alpha}^{*},\pi R_{\alpha}^t)+\mathcal{W}(\pi_{\alpha}^{*}R_{\alpha}^t, \pi_{\alpha}^{*}R_{\alpha})$$
and leverage Eq. (<ref>), limit $t\rightarrow +\infty$. This indicates $\pi_\alpha =\pi_{\alpha}^{*} R_{\alpha}$.

To this end, we begin to investigate the spectral radius of matrix $\mathbf{T}$, denote the eigenvalues of $\mathbf{T}$  as $\lambda$:

\begin{equation}\nonumber
\begin{aligned}
& \text{det} \begin{bmatrix}
(1+\beta)\mathbf{I}-\alpha(1-\beta)\mathbf{A}-\lambda\mathbf{I} & -\beta\mathbf{I} \\
\mathbf{I} & -\lambda\mathbf{I}
\end{bmatrix} = 0 \\ &\Rightarrow \text{det}\begin{bmatrix}
(1+\beta-\lambda-\frac{\beta}{\lambda})\mathbf{I}-\alpha(1-\beta)\mathbf{A} & \mathbf{0} \\
\mathbf{I} & -\lambda\mathbf{I}
\end{bmatrix}=0 \\
& \Rightarrow \text{det}
\end{aligned}
\end{equation}

This indicates $\lambda +\frac{\beta}{\lambda}$ is an eigenvalue of $(1+\beta)\mathbf{I}-\alpha(1-\beta)\mathbf{A}$, we denote $\lambda^*$ being any eigenvalue of $(1+\beta)\mathbf{I}-\alpha(1-\beta)\mathbf{A}$.
By [41], we have

\begin{equation}\label{lambda_max_T}
\lambda_{\text{max}}\leq \frac{1-\alpha\mu}{1+\alpha\mu}
\end{equation}

Proof of Proposition <ref>
We denote the optimize choice of $\mathbf{x}$ as $\mathbf{x}^*$, $\hat{\mathbf{z}}^{k+1}$ as $\begin{bmatrix}
\mathbf{x}^{k+1}-\mathbf{x}^* \\
\mathbf{x}^{k} -\mathbf{x}^*
\end{bmatrix}$, and $\mathbf{z}^*$ as $\begin{bmatrix}
\mathbf{x}^*\\
\mathbf{x}^*
\end{bmatrix}$.

In the following, we start to prove the mean of the stationary distribution $\overline{\mathbf{z}} =\underset{\mathbb{R}^{2d}}{\int}\mathbf{z}\pi_{\alpha}^{*}(d\mathbf{z}) = \mathbf{z}^*+\alpha\triangle +\mathcal{O}(\alpha^2)$.

\begin{equation}\nonumber
\hat{\mathbf{z}}^{k+1} =\begin{bmatrix}
\mathbf{x}^{k+1}-\mathbf{x}^{*}\\
\mathbf{x}^{k}-\mathbf{x}^{*}
\end{bmatrix} = \mathbf{T}\hat{\mathbf{z}}^{k}+\begin{bmatrix}
\mathbf{0}
\end{bmatrix}+\varepsilon_{k}(\mathbf{z}^{k}).
\end{equation}

\begin{equation}\nonumber
(\mathbf{z}^{k+1})^{\otimes 2} = \Big( \mathbf{T}\hat{\mathbf{z}}^{k}+
\begin{bmatrix}
\alpha(1-\beta)\mathbf{b}\\
\mathbf{0}
\end{bmatrix}+\varepsilon_{k}(\hat{\mathbf{z}}^{k})\Big )^{\otimes 2}
\end{equation}

Taking expectation on both sides, and $\mathbf{z}^{k}$ is distributed according to $\pi_{\alpha}^{*}$, we derive:

\begin{equation}\nonumber
\begin{aligned}
\underset{\mathbb{R}^{2d}}{\int}(\mathbf{z}-\mathbf{z}^{*})^{\otimes 2}\pi_{\alpha}^{*}(d\mathbf{z}) = \mathbf{T}\underset{\mathbb{R}^{2d}}{\int}(\mathbf{z}-\mathbf{z}^{*})^{\otimes 2}\pi_{\alpha}^{*}(d\mathbf{z})\mathbf{T}^T +\alpha\int\mathcal{C}(\mathbf{z})\pi_{\alpha}^{*}d(\mathbf{z})
\end{aligned}
\end{equation}

\begin{equation}\nonumber
(\mathbf{I}-\mathbf{T}\mathbf{T}^T) \Big[\underset{\mathbb{R}^{2d}}{\int}(\mathbf{z}-\mathbf{z}^{*})^{\otimes 2}\pi_{\alpha}^{*}(d\mathbf{z})\Big] =
\alpha\int\mathcal{C}(\mathbf{z})\pi_{\alpha}^{*}d(\mathbf{z}).
\end{equation}

It remains to show that $\mathbf{I}-\mathbf{T}\mathbf{T}^T$ is invertible.

Which further yields

\begin{equation}\nonumber
\overline{\mathbf{z}} =\underset{\mathbb{R}^{2d}}{\int}\mathbf{z}\pi_{\alpha}^{*}(d\mathbf{z}) = \mathbf{z}^*+\alpha\triangle +\mathcal{O}(\alpha^2).
\end{equation}

Proof of Proposition <ref>.
Let $\alpha \in (0,L/2)$, and $(\mathbf{x}^k)_{k\ge 0}$ generated by Eq. (<ref>), and $\mathbf{x}_{0}$ distributed according to $\pi_{\alpha}$. For given function $f$, we do a Taylor expansion arround $x^*$, then we derive

\begin{equation}\nonumber
\nabla f(\mathbf{x}) = \nabla^2 f(\mathbf{x}^{*})(\mathbf{x}-\mathbf{x}^{*})+\frac{1}{2}\nabla^3 f(\mathbf{x}^{*})(\mathbf{x}-\mathbf{x}^*)^{\otimes 2}+\mathcal{R}_3(\mathbf{x}),
\end{equation}

where $\mathcal{R}_3: \mathbb R^d\rightarrow \mathbb R^d$ satisfies

\begin{equation}\label{taylor_3rd_expand}
\underset{\mathbf{x}\in \mathbb R^d}{\text{sup}}\{ \|\mathcal{R}_3(\mathbf{x})\|/\|\mathbf{x}-\mathbf{x}^{*}\|\}<+\infty
\end{equation}

Leverage Lemma <ref>, take the integral on both sides of Eq. (<ref>), with respect to $\pi_{\alpha}$, we are led to

\begin{equation}\nonumber
0 = \underset{\mathbb R^d}{\int}\Big(\nabla^2 f(\mathbf{x}^*)(\mathbf{x}-\mathbf{x}^*) +\frac{1}{2}\nabla^3(\mathbf{x}-\mathbf{x}^{*})^{\otimes 2} +\mathcal{R}_3(\mathbf{x})\Big)\pi_{\alpha}d(\mathbf{x}) .
\end{equation}

Then using Lemma <ref> and Holder's inequality, we are led to

\begin{equation}\label{expoalpha^3/2}
\nabla^2 f(\mathbf{x})(\overline{\mathbf{x}}-\mathbf{x}^{*})+\frac{1}{2}\nabla^2 f(\mathbf{x})\Big[\underset{\mathbb R^d}{\int}(\mathbf{x}-\mathbf{x}^{*})^{\otimes 2}\pi_{\alpha}d(\mathbf{x})\Big]=\mathcal{O}(\alpha^{\frac{3}{2}}).
\end{equation}

Besides, it is straightforward

\begin{equation}\nonumber
\mathbf{x}^{2} = \mathbf{x}^{1} - \alpha(1-\beta)\nabla f(\mathbf{x})+\beta(\mathbf{x}^{1}-\mathbf{x}^{0})+\varepsilon_2(\mathbf{x}^{1}).
\end{equation}

By Taylor's second moment expansion, we have

\begin{equation}\nonumber
\begin{aligned}
\mathbf{x}^{2} =\mathbf{x}^{1} -\alpha(1-\beta)\Big[ \nabla^2 f(\mathbf{x}^{1})(\mathbf{x}^{1}-\mathbf{x}^{*})+\mathcal{R}_2(\mathbf{x}^{1})\Big]+\beta(\mathbf{x}^{1}-\mathbf{x}^{0})+\varepsilon_2(\mathbf{x}).
\end{aligned}
\end{equation}

Taking the second moment of the above equation, $\mathbf{x}_{1}$ is independent of $\varepsilon_2$, and Holder's inequality, we derive

\begin{equation}\nonumber
\begin{aligned}
&\underset{\mathbb R^d}{\int}(\mathbf{x}-\mathbf{x}^{*})^{\otimes 2}\pi_{\alpha}d(\mathbf{x}) \\&= \Big(\mathbf{I}-\alpha(1-\beta)\nabla^2 f(\mathbf{x}^{*})\Big) \times \underset{\mathbb R^d}{\int}(\mathbf{x}-\mathbf{x}^{*})^{\otimes 2}\pi_{\alpha}d(\mathbf{x})\Big(\mathbf{I}-\alpha(1-\beta)\nabla^2 f(\mathbf{x}^{*})\Big) \\ &+ \alpha^2\underset{\mathbb R^d}{\int}\mathcal{C}(\mathbf{x})\pi_{\alpha}d(\mathbf{x})+\mathcal{O}(\alpha^{\frac{5}{2}}).
\end{aligned}
\end{equation}

This indicates

\begin{equation}\label{oalpha^3/2}
\begin{aligned}
& \underset{\mathbb R^d}{\int}(\mathbf{x}-\mathbf{x}^{*})^{\otimes 2}\pi_{\alpha}d(\mathbf{x}) \\&= \Big(\mathbf{I} \otimes \nabla f(\mathbf{x}^{*})+\nabla f(\mathbf{x}^{*})\otimes\mathbf{I}-\alpha(1-\beta)^2\nabla f(\mathbf{x}^{*})^{\otimes 2}\Big)^{-1} \times \underset{\mathbb R^d}{\int}\mathcal{C}(\mathbf{x})\pi_{\alpha}d(\mathbf{x})+\mathcal{O}(\alpha^{\frac{3}{2}}).
\end{aligned}
\end{equation}

Now we show that the matrix
$$\mathbf{I} \otimes \nabla f(\mathbf{x}^{*})+\nabla f(\mathbf{x}^{*})\otimes\mathbf{I}-\alpha(1-\beta)^2\nabla f(\mathbf{x}^{*})^{\otimes 2}$$
is invertible, by the Assumption <ref>, and $\mu \in (0,2/L)$ we have that $\mathbf{I}-(\alpha/2)(1-\beta)^2 \nabla f(\mathbf{x})\otimes$ is symmetric positive definite.

\begin{equation}\nonumber
\begin{aligned}
&\mathbf{I} \otimes \nabla f(\mathbf{x})+\nabla f(\mathbf{x})\otimes\mathbf{I}-\alpha(1-\beta)^2\nabla f(\mathbf{x})^{\otimes 2} \\&= \nabla f(\mathbf{x})\otimes\Big(\mathbf{I}-\frac{\alpha}{2}(1-\beta)^2\nabla f(\mathbf{x})\Big)+\Big(\mathbf{I}-\frac{\alpha}{2}(1-\beta)^2\nabla f(\mathbf{x})\Big)\otimes\nabla f(\mathbf{x}).
\end{aligned}
\end{equation}

We can see that the eigenvalues of $\mathbf{I} \otimes \nabla f(\mathbf{x}^{*})+\nabla f(\mathbf{x}^{*})\otimes\mathbf{I}-\alpha(1-\beta)^2\nabla f(\mathbf{x}^{*})^{\otimes 2}$ is positive, then this matrix is invertible. Combine Eq. (<ref>) and (<ref>), we have Proposition <ref>.

§ IMPLEMENTION DETAILS

§.§ Additional Denoising Step

Non-denoised/denoised FIDs and the initial step size with different NFEs. Results obtained from score net NCSN2.

Steps      50     100    150    200    250     500     1000

2*ALS    $72.75/42.74$     $49.78/18.35$   $43.95/3.87$     $40.50/12.04$   $38.14/10.98$    $33.55/\mathbf{9.33}$     $30.17/\mathbf{8.74}$

$1.07\times 10^{-4}$    $4.5\times 10^{-5}$    $3\times 10^{-5}$    $2.2\times 10^{-5}$    $1.8\times 10^{-5}$    $1\times 10^{-5}$    $6\times 10^{-6}$

2*AMS    $64.11/\mathbf{30.52}$     $41.74/\mathbf{13.46}$    $ 34.18 /\mathbf{11.35}$    $ 29.37/\mathbf{10.95}$    $27.14/\mathbf{10.84}$    $18.94/10.67$     $14.51/9.82$

$7.5 \times 10^{-5}$    $2.6\times 10^{-5}$    $1.2\times 10^{-5}$    $1.2\times 10^{-5}$    $5\times 10^{-6}$    $2.2\times 10^{-6}$    $1\times 10^{-6}$

Time     18.64     34.21    48.87     64.00     79.03     154.28     306.78

* We present the whole running time (second) for 500 image generation. We set $n_\sigma = 5$.

Following the suggestion by [12, 32], and added a denoise step at the end of Langevin sampling. Assume a noise-conditional score net $s_{\theta}(\mathbf{x},\sigma)$, and a noisy sample data $\mathbf{x}$. We can compute the Expected Denoised Sample (EDS) $\tilde{\mathbf{x}}$ by Empirical Bayes mean [28]

\begin{equation}
s_{\theta}(\mathbf{x},\sigma) = \frac{\tilde{\mathbf{x}} - \mathbf{x}}{\sigma^2} \rightarrow \tilde{\mathbf{x}} = \sigma^2 s_{\theta}(\mathbf{x},\sigma) + \mathbf{x},
\end{equation}

if the score net is well-trained. For the noise-conditional score net, the denoising step is straightforward to apply.

Many experiments [12, 32] have shown adding a denoise step at the end of the Langevin dynamic can improve the score of generated samples in terms of FID (perhaps there will be a negative impact on inception score). This trick also works in AMS, we show the synthetic 2D comparative experiments in Figure <ref>. However, our proposed AMS can lower the gap between those FIDs in non-denoised/denoised figures, we conclude there is no need for the additional denoising step if we adopt large NFEs.

[Swiss Roll dataset]

[25 Gaussians dataset]
AMS generated samples after denoising on synthetic 2D experiments. For each subplot, the left is the data point without denoise, the right is the expected denoised samples. The orange point is real data, and the green is generated points.

§.§ Parameters

Used hyperparameters for Graph sampling/evaluation.

Ego-small ($300$ NFE)     $\epsilon = 0.11,\ \varepsilon= 0.5$

Ego-small ($750$ NFE)     $\epsilon = 0.02,\ \varepsilon = 1.0$

Community-small ($300$ NFE)     $\epsilon = 0.02,\ \varepsilon = 0.8$

Community-small ($750$ NFE)    $\epsilon = 0.08,\ \varepsilon = 1.0$

Grid     $\epsilon = 0.15,\ \varepsilon = 0.7$

Used hyperparameters for NCSN++ cont. (VE) sampling/evaluation.

Sampler     NFE $=150$      NFE $=195$      NFE $=210$     NFE $=270$     NFE $=300$

RD-MC     $\epsilon=0.20,n_\sigma=2$    $\epsilon=0.19,n_\sigma=2$    $\epsilon=0.21,n_\sigma=2$      $\epsilon=0.29,n_\sigma=2$    $\epsilon=0.265,n_\sigma=2$

RD-CDL     $\epsilon=0.07,n_\sigma=2$    $\epsilon=0.06,n_\sigma=2$    $\epsilon=0.05,n_\sigma=2$      $\epsilon=0.04,n_\sigma=2$    $\epsilon=0.03,n_\sigma=2$

NFE $=150$     NFE $=270$     NFE $=500$     NFE $=1000$     NFE $=2000$

MC    $\epsilon=0.30,n_\sigma=5$     $\epsilon=0.20,n_\sigma=5$    $\epsilon=0.085,n_\sigma=5$     $\epsilon=0.05,n_\sigma=5$    $\epsilon=0.03152,n_\sigma=5$

LC    $\epsilon=0.44,n_\sigma=5$     $\epsilon=0.38,n_\sigma=5$    $\epsilon=0.32,n_\sigma=5$    $\epsilon=0.26,n_\sigma=5$   $\epsilon=0.22,n_\sigma=5$

CDL    $\epsilon=0.12,n_\sigma=5$     $\epsilon=0.08,n_\sigma=5$    $\epsilon=0.0561,n_\sigma=5$     $\epsilon=0.03975,n_\sigma=5$    -

RD-LC     $\epsilon=0.18,n_\sigma=1$    $\epsilon=0.20,n_\sigma=1$    $\epsilon=0.18,n_\sigma=1$     $\epsilon=0.18,n_\sigma=1$    $\epsilon=0.16,n_\sigma=1$

EM-LC     $\epsilon=0.24,n_\sigma=1$    $\epsilon=0.18,n_\sigma=1$    $\epsilon=0.18,n_\sigma=1$     $\epsilon=0.18,n_\sigma=1$    $\epsilon=0.16,n_\sigma=1$

In this section, we present the necessary hyperparameter setting to enhance the reproducibility. The initial step size is a key parameter for NCSNv2 sampling, to find out the best initial step size for FID comparison, we swept over the sampling step size with approximately 2 significant numbers of precision, the parameter used in Table <ref> is as follows. For all the FID tests, we used a sampling batch size of 1k a total FID sample number of 50k, and partial FID for 2k samples in total. The evaluation for our AMS sampler adopted $n_\sigma = 5$ during the whole experiment process, thus we found using $n_\sigma = 5$ instead of $n_\sigma = 1$ can lead to better FID performance for both ALS and AMS sampler (we should divide the total sampling step by five if we use $n_\sigma = 5$ instead of $n_\sigma = 1$ ), see Table <ref>. For all the sampling/evaluation processes, we added a denoising step using Tweedie's formula [28].

For the SDE sampler, we use the same score net provided by [34], e.g. NCSN++ cont., and DDPM++ cont. The hyperparameters in the unified SGM framework are mainly two: $\epsilon, n_\sigma$, $\epsilon$ can be an initial step size hyperparameter, like the tactics we used in NCSN2 evaluation, we first ran partial FID estimation to find out the best $\epsilon$ hyperparameter setting and adopt this best $\epsilon$ for the full FID evaluation. Similarly, the $n_\sigma$ can leverage this method too, and $n_\sigma=2$ may be the best choice for the small NFEs sampling. For the MC/LC samplers, we directly use $n_\sigma=5$ during the whole experiments, since we have tested the corrector-only sampler in NCSN2 experiments. We found our momentum corrector requires a much smaller $\epsilon$ than the Langevin corrector, this is mainly because of our S-N-R schedule, we multiplied the step size by $(1+\beta)^2$ during our momentum sampling/evaluation experiments. The Langevin corrector only improves the performance of the predictor in large NFEs, the more Langevin correction steps we use, the worse FID we will gain in those small NFEs, e.g. 150, 270. Note that for those samplers only with predictors, e.g. EM or RD, there is no impact for those hyperparameters.

§.§ Choosing the  and

In NCSN2 experiments, it is clear to see that an increase $n_\sigma$ will lead to better performance if we adopt ALS or AMS as a sampler, see Table <ref>. We test the $\epsilon$ and $n_\sigma$ by running partial FID estimation (1k or 2k samples) and find out the best choice for these two hyperparameters. For the Langevin corrector, we found that adding additional Langivin correct steps does not benefit the FID scores in NCSN++ cont. for any $n_\sigma$ settings in those small NFEs, while the Momentum corrector derives the contrast, improved the FID compared to those only predictor samplers. We present the details in Figure <ref>. We note the $\epsilon$-FID learning curve of RD-MC sampler with NFE $=270$ waves between an interval, our experiment adopts the second trough since this $n_\sigma$ can lead to a better inception score and similar FID.

§.§ Extend to SDE

RD-MC sampler for VE-SDE

: hyperparameters $(\epsilon)_{i=0}^{n}>0$, $\delta>0$, score-net $s_{\theta}$, noise scale $(\sigma_i)_{i=1}^{n}$, $n_{\sigma}$.

Initialize $\mathbf{x}$.
$i \leftarrow 1$ to $n$

Draw Gaussian noise $\varepsilon \sim \mathcal{N}(0,\mathbf{I}).$
$\mathbf{x} \leftarrow \mathbf{x} + (\sigma_{i+1}^2-\sigma_{i}^2)s_{\theta}(\mathbf{x},\sigma_{i+1}).$
$\mathbf{x} \leftarrow \mathbf{x} + \sqrt{\sigma_{t+1}^2-\sigma_{t}^2}\mathbf{\varepsilon}.$
$k \leftarrow 1$ to $n_{\sigma}$
Draw Gaussian noise $\mathbf{z} \sim \mathcal{N}(0,\mathbf{I}).$
$\alpha_i = 2 \big(\epsilon_i (1 + \beta) ^ 2\| \mathbf{m}^{i}\|_2 / \| \mathbf{z} \|_2\big )^{2}.$
Calculate $\beta$ and Update $\mathbf{x}$ by Eq. (<ref>), and (<ref>).
return $\mathbf{x}+\sigma_{n}^{2}s_{\theta}(\mathbf{x},\sigma_n).$
return $\mathbf{x}.$

RD-MC sampler for VP-SDE

: hyperparameters $(\epsilon)_{i=0}^{n}>0$, $\delta>0$, score-net $s_{\theta}$, $n_{\sigma}$, noise scale $(\sigma_i)_{i=1}^{n}$, and $(\beta_i)_{i=1}^{n}$.
Initialize $\mathbf{x}$.
$i \leftarrow 1$ to $n$
Draw Gaussian noise $\varepsilon \sim \mathcal{N}(0,\mathbf{I}).$
$\mathbf{x} \leftarrow (2-\sqrt{1-\beta_{i+1}})\mathbf{x}+\beta_{i+1}s_{\theta}(\mathbf{x},i+1).$
$\mathbf{x} \leftarrow \mathbf{x}+\sqrt{\beta_{i+1}}\varepsilon.$
$k \leftarrow 1$ to $n_{\sigma}$
Draw Gaussian noise $\mathbf{z} \sim \mathcal{N}(0,\mathbf{I}).$
$\alpha_i = 2 \epsilon_i \big (\epsilon_0 (1 + \beta) ^2 \| \mathbf{m}^{i}\|_2 / \| \mathbf{z} \|_2 \big )^{2}.$
Calculate $\beta$ and Update $\mathbf{x}$ by Eq. (<ref>), and (<ref>).
return $\mathbf{x}+\sigma_{n}^{2}s_{\theta}(\mathbf{x},\sigma_n).$
return $\mathbf{x}.$

[34] connected the forward process to SDE, and solved the reverse SDE with reverse diffusion predictor. Moreover [34] modified the traditional annealed Langevin dynamic in the step size schedule and provided the new annealed Langevin dynamic corrector. For the forward process in SMLD, we have
\begin{equation}\nonumber
\mathbf{x}^{t+1} = \mathbf{x}^{t} + \sqrt{\sigma_{t+1}^2-\sigma_{t}^2}\mathbf{n}^{t},
\end{equation}

where $\mathbf{n}^t \sim \mathcal{N}(0,\mathbf{I})$, this process can be formulated by an SDE with the form

\begin{equation}\nonumber
d\mathbf{x} = f(\mathbf{x,t})dt+g(t)d\mathbf{w},
\end{equation}

where $\mathbf{w}$ is the standard Wiener process, $f(\cdot, t): \mathbb R^d \rightarrow \mathbb R^d$ is a function dipict the diffusion process of $\mathbf{x}^{t}$, and $g(\cdot):\mathbb R\rightarrow R$ is a scalar function depend on the diffusion time $t$. [34] provide a unified framework to solve the reverse SDE, we present the pseudo-code when using the predictor of reverse diffusion predictor and corrector of MC for NCSN (VE) and DDPM (VP) in Algorithm <ref>, and Algorithm <ref>. We present the perturbation kernel of SMLD (VESDE) and DDPM (VPSDE) in [34], and the corresponding SDE as follows:

\begin{equation}\nonumber
\begin{aligned}
\mathbf{x}_{i} = \mathbf{x}_{i-1} + \sqrt{\sigma_i^2 - \sigma_{i-1}^2}\mathbf{\varepsilon_{i-1}} &\Rightarrow d\mathbf{x} = \sqrt{\frac{d[\sigma^2(t)]}{dt}}d\mathbf{w}, \\
\mathbf{x}_{i} = \sqrt{1-\beta_i}\mathbf{x}_{i-1} + \sqrt{\beta_1} \mathbf{\varepsilon_{i-1}}&\Rightarrow d\mathbf{x} = -\frac{1}{2}\beta(t)\mathbf{x}dt + \sqrt{\beta(t)}d\mathbf{w}.
\end{aligned}
\end{equation}

Besides, [34] also proposed a new type of SDE that performs particularly well on the likelihood named Sub-VP SDE

\begin{equation}\nonumber
\mathrm{d}\mathbf{x}=-\frac12\beta(t)\mathbf{x}\mathrm{~d}t+\sqrt{\beta(t)(1-e^{-2\int_0^t\beta(s)\mathrm{d}s})}\mathrm{d}\mathbf{w}.
\end{equation}

§ LIMITATIONS

Our paper proposed AMS, extending this sampler into an SDE corrector without introducing any hyperparameters compared to the standard Langevin dynamic. However, the theoretical development of SDE is not as extensive as that of ODE, which makes the generation speed of stochastic samplers significantly slower compared to deterministic samplers [20, 39]. When using extremely small NFEs e.g., 30, and 50, the performance gap between stochastic samplers and deterministic samplers becomes significant. Additionally, we observe that the predictor plays a more crucial role than the corrector in the generation process, although our MC substantially enhances generation quality in those small NFEs. It is important to note that this improvement necessitates a finely tuned hyperparameter $\epsilon$. Similar to standard Langevin dynamics, the performance of our MC is highly sensitive to hyperparameters.

We identify several potential avenues to address the aforementioned challenges. Beyond Langevin dynamics, we observe that the predictor also shares common iterations with SGD. This suggests that integrating momentum into the predictor, rather than the corrector, could be beneficial. Furthermore, the hyperparameter $\epsilon$ is typically kept constant throughout the sampling process. Developing a specialized step-size schedule for existing Langevin dynamics might yield improved results.  Additionally, momentum has shown promise in enhancing ODE-based methods as well  [26, 39, 42]. Exploring these directions could provide valuable insights and advancements in sampling efficiency and quality.

§ SUPPLEMENT EXPERIMENTAL RESULTS

§.§ NCSN2 Experiments

We test our AMS (Eq. (<ref>)) sampling in score-net architecture in [33] which is named NCSN2, and similar score-net architecture with an unconditional version (UNET) in [6, 12]. Since the score-net architecture can obtain lower FID than NCSN2 but may take more time to train and sample. Using the ALS (Eq. (<ref>)) as a baseline method, we ran our experiments in CIFAR-10 [16]. We compared those factors like sampling step size, $n_{\sigma}$, initial step size, and denoising. Details on how those experiments were conducted can be found in Appendix <ref>.

For analyzing the performance of our sampler in different Numbers of score Function Evaluations (NFEs) [11], we ran AMS and ALS in NFEs e.g. $50,\ 100,\ 150,\ 200,\ 250,\ 500,\ \text{and} \ 1000$. FID results are shown in Table <ref>. We see that AMS is more suitable in smaller NFEs e.g. $50,\ 100$. AMS outperforms ALS heavily in these step sizes, e.g. $42.74,\ \text{v.s.}\ 30.52$. Besides, in all the tested NFEs, non-denoised FID with AMS can be far ahead of those FID with ALS. Since the FID score in large NFE sizes is small enough, all those two samplers reached the upper limit of what the NCSN2 can reach. As the NFE is larger than $250$, both AMS and ALS have limited gains in denoised FID, the gap between non-denoised and denoised FID with AMS shrinking as the number of steps increases, improving the consistency of those FIDs. AMS obtains higher FIDs in large NFEs, this result coincident with those in [41] that momentum can accelerate the training process. However, there are no visible differences between the figures sampled by AMS and ALS. Thus, AMS can act as an accelerator sampler that works better in small NFEs and saves sampling time.

Following [32], we focus on the widely used CIFAR10 unconditional SGM benchmarks, e.g. NCSN [32] with ALS and NCSN2 [33] with ALS and Consistent Annealed Sampling (CAS) [12], and those GANs from the same period with NCSN2. We present our result in Table <ref>. Our AMS reaches an FID of $14.51/9.82$ for without/with denoising in NCSN2. Note that NCSN2 with CAS and NCSN with ALS used the NFE of $1160\ (232 \times 5)$, and our results were obtained by $1000\ (200 \times 5)$ steps. AMS outperforms ALS in consistency and CAS in the final denoised FID, more like ALS meets CAS. We have to figure out that AMS is designed for small NFEs and can lead to better performance with those NFEs, e.g. $50,\ 100,\ 150,\ 200$, see Table <ref> and Figure <ref> for details.

For choosing noise scale $n_\sigma$, we test $n_\sigma = 1$, and $n_\sigma = 5$ for NFEs of $50,\ 150,\ 250$ with ALS, if we attribute $n_\sigma = 5$, we will have $n=10$ for total of $50$ NFEs. As shown in Table <ref>, we can derive that if we adopt $n_\sigma = 5$, we can achieve a lower FID within the same step level. This conclusion is available for the AMS sampler.

FIDs obtained from NCSN2 with different NFEs. AMS can lead to better FIDs with those small NFEs, e.g. 50, 100, 150, 200, 250. Besides, AMS is only slightly behind ALS in large NFEs. For experimental details, please refer to Appendix <ref>.

[50 NFE]
[100 NFE]
[150 NFE]

[200 NFE]
[500 NFE]
[1000 NFE]
CIFAR10 samples for different NFEs, result obtained with our sampler from NCSN2.

Non-denoised/Denoised FIDs with $n_\sigma$ in $\{1,5\}$ for the NCSN2.

50     150     250

$n_\sigma = 1$     $79.39/57.71$    $48.10/17.16$    $41.19/12.29$

$n_\sigma = 5$     $72.75/\textbf{42.74}$    $43.95/\textbf{13.87}$    $38.14/\textbf{10.98}$

Unconditional CIFAR-10 generative performance.

Model     Sampler     FID $\downarrow$

2*NCSN2      AMS (ours), w/o denoising      $14.51$

AMS (ours), w denoising       $9.82$

2*UNET      AMS (ours), w/o denoising     $6.89$

AMS (ours), w denoising      $4.84$

4*NCSN2      CAS, w/o denoising       $12.7$

CAS, w denoising       $11.2$

ALS, w/o denoising       $30.17$

ALS, w denoising       $8.74$

2*UNET      ALS, w/o denoising       $16.03$

ALS, w denoising      $4.24$

Quantitative results on Community-small with different solver/NFEs.

NFE     Methods     Deg.     Clus.     Orbit      Time/s

4*300     EM     $0.074$    $0.111$    $\underline{0.007}$    $2.99$

EM-LD    $ \underline{0.061}$    $ 0.122$    $ 0.031$     $ 2.86$

EM-CLD     $0.065$    $0.122$   $0.031$    $2.94$

EM-MC (ours)    $0.067$   $\underline{0.109}$    $0.031$     $2.83$

4*450     EM     $0.049$    $0.121$    $0.013$    $4.28$

EM-LD     $0.042$    $0.071$    $0.005$    $4.09$

EM-CLD     $0.043$    $0.068$   $0.005$    $4.27$

EM-MC (ours)     $\underline{0.041}$    $\underline{0.053}$    $\underline{0.004}$    $4.10$

4*600     EM     $0.066$    $0.173$    $0.011$    $5.19$

EM-LD     $0.025$    $0.124$   $\underline{0.004}$     $5.21$

EM-CLD     $0.046$    $0.092$   $0.006$    $5.40$

EM-MC (ours)    $\underline{0.037}$    $\underline{0.048}$    $\underline{0.004}$    $5.26$

4*750     EM     $0.451$    $0.123$    $0.006$    $6.63$

EM-LD     $0.050$    $0.065$   $0.004$    $6.41$

EM-CLD    $0.047$    $0.072$   $\underline{\textbf{0.003}}$    $7.03$

EM-MC (ours)   $\underline{\textbf{0.030}}$     $\underline{\textbf{0.042}}$     $\underline{\textbf{0.003}}$     $6.88$

IS, KID $\times 10^{-3}$ for different NFEs on CIFAR10, results obtained with different samplers in model NSCN++ cont. (VE).

Sampler     NFE = $150$     NFE = $270$    NFE $=500$    NFE $=1000$    NFE $=2000$

MC (ours)     $2.77,\ 246$    $8.65,\ 19.7$    $9.97,\ 7.16$     $10.1,\ 13.2$   $9.74,\ 5.71$

EM-LC     $9.12,\ 12.4$    $9.89,\ 9.51$    $10.3,\ 10.4$    $10.4,\ 10.8$   $11.6,\ 11.1$

RD-LC     $7.46,\ 34.9$    $9.11,\ 4.68$    $9.56,\ 1.25$    $9.73,\ 0.60$   $9.87,\ 0.48$

LC    $1.85,\ 308$     $3.31,\ 198$   $5.37,\ 99.6$     $7.55,\ 41.1$   $9.20,\ 16.9$

EM     $9.54,\ 11.1$     $9.89,\ 10.6$    $10.0,\ 11.0$   $10.2,\ 11.4$    $10.1,\ 1.70$

RD     $7.71,\ 29.8$    $9.17,\ 4.11$    $9.39,\ 1.55$    $9.40,\ 1.22$    $9.46,\ 1.25$

Gotta Go Fast     $9.87,\ -$    $9.64,\ -$     $9.57,\ -$    $-,\ -$    $-,\ -$

CLD ($\beta = 0.9$)     $2.03,\ 298$    $9.45,\ 24.0$    $9.63,\ 17.4$    $9.08,\ 18.4$    -

FID, IS, KID $\times 10^{-3}$ scores obtained from sampler RD-MC with different NFEs in model NSCN++ cont. (VE).

Scores     Sampler     NFE $=150$      NFE $=195$      NFE $=210$     NFE $=270$     NFE $=300$

2*FID      RD-MC    $6.75 $     $ 3.05$    $ 2.86$     $2.85$    $2.93$

RD-CDL     $ 38.6 $     $ 5.04$    $ 4.37$   $ 4.80$     $ 4.43$

2*IS     RD-MC     $9.07 $     $ 9.51$    $ 9.55$     $9.69$    $9.70$

RD-CDL     $7.69 $     $ 9.25$    $ 9.15$     $9.32$    $9.25$

2*KID    RD-MC     $3.13 $     $ 1.06$    $ 1.04$     $0.91$    $0.86$

RD-CDL     $25.3 $     $ 2.99$    $ 2.27$     $3.22$    $2.13$

KID (lower is better) $\times 10^{-3}$ for different NFEs on CIFAR10, results obtained with different samplers in model NSCN++ cont. (VE).

Sampler     NFE = $150$     NFE = $270$    NFE $=500$    NFE $=1000$    NFE $=2000$

MC (ours)    $246$     $19.7$    $7.16$    $13.2$    $5.71$

EM-LC     $12.4$    $9.51$    $10.4$    $10.8$   $11.1$

RD-LC     $34.9$    $4.68$    $1.25$   $0.60$   $0.49$

LC     $308$    $198$    $99.6$    $41.1$   $16.9$

EM     $11.1$     $10.6$    $11.0$    $11.4$    $1.70$

RD     $29.8$    $4.11$     $1.55$    $1.22$    $1.25$

CLD ($\beta = 0.9$)     $298$   $24.0$    $17.4$     $18.4$    -

[RD-LC sampler, 270 NFE]
[RD-MC sampler, 270 NFE]
[EM-MC sampler, 270 NFE]
Partial estimate of FID (2k samples) with different initial snrs
, NFE $=270$.

FIDs (50k samples) for different NFEs on CIFAR10, results obtained with different samplers in model NSCN++ cont. (VE).

Sampler     NFE = 150     NFE = 270     NFE $= 500$     NFE $=1000$     NFE $=2000$

MC (ours)    $235$     $27.6$   $9.16$    $17.7$    $9.29$

EM-LC     $17.2$     $11.3$    $11.5$    $11.6$    $11.6$

RD-LC     $47.6$    $8.87$     $3.62$     $2.62$    $2.38$

LC    $286$     $202$    $121$    $57.8$     $25.9$

EM     $13.6$     $12.1$    $12.3$    $12.6$    $12.9$

RD     $40.3$     $7.51$    $3.72$   $2.98$     $2.89$

CLD ($\beta = 0.9$)     $273$    $27.9$   $23.7$     $32.4$    -

Gotta Go Fast     $8.85$     $3.23$     $2.87$     -    -

EM (VP)     $30.3$     $13.1$    $ 4.42$     $ 2.46$     $2.43$

DDIM (VP)     $4.03$     $ 3.53$     $3.26$    $3.09$     $3.01$

§ UNCURATED SAMPLES

Additional samples on CIFAR10 from NCSN++ cont. (VE) using RD-MC sampler with NFE $=210$, we report an FID of $2.86$ with this setting.

Additional samples on CIFAR10 from NCSN++ cont. (VE) using MC sampler with NFE $=500$, we report an FID of $9.16$ with this setting.

Additional samples on CIFAR10 from NCSNv2 using AMS sampler with NFE $=1000$. We report an FID of $9.82$ with this setting.

Additional samples on CIFAR10 from DDPM++ cont. (VP) using RD-MC sampler with NFE $=2000$.

Additional samples on $256 \times 256$ CelebA-HQ from NCSN++ cont. (VE) using RD-MC sampler with NFE $=450$.

Additional samples on $256 \times 256$ CelebA-HQ from NCSN++ cont. (VE) using RD-MC sampler with NFE $=900$.

Additional samples on $256 \times 256$ bedroom from NCSN++ cont. (VE) using RD-MC sampler with NFE $=500$.

Additional samples on $256 \times 256$ church from NCSN++ cont. (VE) using RD-MC sampler with NFE $=300$.