$10^{5}$ | 5.78$\times 10^{-3}$ secs | 1.59$\times 10^{-31}$ | 2.12$\times 10^{-1}$ secs | 6.05$\times 10^{-31}$
$10^{6}$ | 5.76$\times 10^{-3}$ secs | 2.23$\times 10^{-31}$ | 2.11$\times 10^{-1}$ secs | 4.69$\times 10^{-31}$
$10^{7}$ | 5.76$\times 10^{-3}$ secs | – | 2.11$\times 10^{-1}$ secs | –
$10^{8}$ | 5.82$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{9}$ | 5.72$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{12}$ | 5.72$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{15}$ | 5.70$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{18}$ | 5.70$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
Table 6:  The evaluation of the modal Green’s function in quadruple precision for a for varying $\kappa$ with large source-to-target distance ($\beta_{-}=1$). The error is evaluated by using adaptive Gaussian quadrature as the gold standard. Note for $\kappa>10^{6}$, the resource requirements of prior methods becomes excessive.  | $\beta_{-}=10^{-12}\,,m=10$ | $\beta_{-}=10^{-12}\,,m=1000$
---|---|---
$\kappa$ | Evaluation Time | Absolute Error | Evaluation Time | Absolute Error
$10^{-6}$ | 4.49$\times 10^{-5}$ secs | 3.08$\times 10^{-13}$ | 1.34$\times 10^{-3}$ secs | 2.90$\times 10^{-11}$
$10^{-3}$ | 4.45$\times 10^{-5}$ secs | 2.90$\times 10^{-13}$ | 1.37$\times 10^{-3}$ secs | 2.88$\times 10^{-11}$
$10^{0}$ | 4.77$\times 10^{-5}$ secs | 1.90$\times 10^{-13}$ | 1.40$\times 10^{-3}$ secs | 2.84$\times 10^{-11}$
$10^{1}$ | 4.79$\times 10^{-5}$ secs | 4.35$\times 10^{-14}$ | 1.41$\times 10^{-3}$ secs | 2.74$\times 10^{-11}$
$10^{2}$ | 4.61$\times 10^{-5}$ secs | 1.80$\times 10^{-14}$ | 1.43$\times 10^{-3}$ secs | 2.29$\times 10^{-11}$
$10^{3}$ | 3.57$\times 10^{-5}$ secs | 1.07$\times 10^{-14}$ | 1.44$\times 10^{-3}$ secs | 4.19$\times 10^{-12}$
$10^{4}$ | 3.49$\times 10^{-5}$ secs | 3.33$\times 10^{-14}$ | 1.43$\times 10^{-3}$ secs | 5.28$\times 10^{-13}$
$10^{5}$ | 3.46$\times 10^{-5}$ secs | 1.50$\times 10^{-13}$ | 1.44$\times 10^{-3}$ secs | 6.58$\times 10^{-13}$
$10^{6}$ | 3.41$\times 10^{-5}$ secs | 5.11$\times 10^{-13}$ | 1.45$\times 10^{-3}$ secs | 3.04$\times 10^{-13}$
$10^{7}$ | 3.48$\times 10^{-5}$ secs | – | 7.88$\times 10^{-4}$ secs | –
$10^{8}$ | 3.43$\times 10^{-5}$ secs | – | 7.87$\times 10^{-4}$ secs | –
$10^{9}$ | 3.38$\times 10^{-5}$ secs | – | 7.87$\times 10^{-4}$ secs | –
$10^{12}$ | 3.22$\times 10^{-5}$ secs | – | 7.85$\times 10^{-4}$ secs | –
$10^{15}$ | 3.38$\times 10^{-5}$ secs | – | 7.87$\times 10^{-4}$ secs | –
$10^{18}$ | 3.33$\times 10^{-5}$ secs | – | 7.85$\times 10^{-4}$ secs | –
Table 7:  The evaluation of the modal Green’s function in double precision for varying $\kappa$ with small source-to-target distance ($\beta_{-}=10^{-12}$). The error is evaluated by using adaptive Gaussian quadrature as the gold standard. Note for $\kappa>10^{6}$, the resource requirements of prior methods becomes excessive.  | $\beta_{-}=10^{-12}\,,m=10$ | $\beta_{-}=10^{-12}\,,m=1000$
---|---|---
$\kappa$ | Evaluation Time | Absolute Error | Evaluation Time | Absolute Error
$10^{-6}$ | 7.23$\times 10^{-3}$ secs | 7.06$\times 10^{-31}$ | 2.08$\times 10^{-1}$ secs | 5.75$\times 10^{-29}$
$10^{-3}$ | 7.53$\times 10^{-3}$ secs | 7.12$\times 10^{-31}$ | 2.08$\times 10^{-1}$ secs | 5.75$\times 10^{-29}$
$10^{0}$ | 7.16$\times 10^{-3}$ secs | 6.54$\times 10^{-31}$ | 2.10$\times 10^{-1}$ secs | 5.70$\times 10^{-29}$
$10^{1}$ | 7.53$\times 10^{-3}$ secs | 8.03$\times 10^{-32}$ | 2.11$\times 10^{-1}$ secs | 5.56$\times 10^{-29}$
$10^{2}$ | 7.09$\times 10^{-3}$ secs | 5.00$\times 10^{-32}$ | 2.11$\times 10^{-1}$ secs | 4.55$\times 10^{-29}$
$10^{3}$ | 6.25$\times 10^{-3}$ secs | 6.41$\times 10^{-32}$ | 2.12$\times 10^{-1}$ secs | 7.93$\times 10^{-30}$
$10^{4}$ | 6.23$\times 10^{-3}$ secs | 2.07$\times 10^{-32}$ | 2.11$\times 10^{-1}$ secs | 1.55$\times 10^{-31}$
$10^{5}$ | 6.21$\times 10^{-3}$ secs | 1.43$\times 10^{-31}$ | 2.12$\times 10^{-1}$ secs | 2.21$\times 10^{-31}$
$10^{6}$ | 6.21$\times 10^{-3}$ secs | 2.37$\times 10^{-31}$ | 2.12$\times 10^{-1}$ secs | 3.88$\times 10^{-31}$
$10^{7}$ | 6.25$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{8}$ | 6.22$\times 10^{-3}$ secs | – | 1.53$\times 10^{-1}$ secs | –
$10^{9}$ | 6.19$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{12}$ | 5.78$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{15}$ | 5.72$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
$10^{18}$ | 5.69$\times 10^{-3}$ secs | – | 1.52$\times 10^{-1}$ secs | –
Table 8:  The evaluation of the modal Green’s function in quadruple precision
for varying $\kappa$ with small source-to-target distance
($\beta_{-}=10^{-12}$). The error is evaluated by using adaptive Gaussian
quadrature as the gold standard. Note for $\kappa>10^{6}$, the resource
requirements of prior methods becomes excessive.

###  5.3 Performance of the Algorithm with Varying Fourier Mode ($m$)

We examined the performance of our algorithm over a wide range of Fourier
modes (represented by the parameter $m$). Because the number of points in the
quadrature scales linearly with $m$, as demonstrated by Table 9, evaluation
time scales linearly with the Fourier mode. Recall from the introduction of
this section that the evaluation was performed on a four-core processor using
two threads.

$m$ | Evaluation Time |
---|---|---
1 | 3.88$\times 10^{-5}$ secs |
10 | 5.56$\times 10^{-5}$ secs |
$10^{2}$ | 1.75$\times 10^{-4}$ secs |
$10^{3}$ | 1.46$\times 10^{-3}$ secs |
$10^{4}$ | 1.43$\times 10^{-2}$ secs |
$10^{5}$ | 1.37$\times 10^{-1}$ secs |
$10^{6}$ | 1.36$\times 10^{0}$ secs |
$10^{7}$ | 1.29$\times 10^{1}$ secs |
Table 9:  The evaluation time of the modal Green’s function in double
precision for varying $m$ ($\beta_{-}=10^{-12},\kappa=10,000$).

###  5.4 Parallelization of the Algorithm

The cost of our algorithm is $O(m)$ and does not depend on $\kappa$ or
$\beta_{-}$ (see Section 4.4). Because our algorithm is quadrature based, it
is embarrassingly parallelizable.

We measured the algorithm’s performance on a server with a 16-core Intel Xeon
2.9 GHz processor, where each core can run two threads for a total of
32-threads. We vary the number of threads from 1 to 32, and report the results
in Figure 11.

Figure 11:  Evaluation time of the modal Green’s function plotted against $m$
with varying numbers of threads ($\beta_{-}=10^{-7},\kappa=10,000$). The
calculation is performed in double precision. The evaluation times
corresponding to 32 threads are not plotted for small $m$.

##  6 Conclusions and Generalizations

We have developed an algorithm which evaluates the modal Green’s function for
the Helmholtz equation in $O(m)$ time, that is completely independent of both
the wavenumber and the source-to-target distance. Furthermore, our algorithm
is embarrassingly parallelizable. Our algorithm’s method can be readily
extended to several associated problems in computational electromagnetics,
described in Sections 6.1\- 6.4.

###  6.1 An $O(1)$ Evaluator for Small Wavenumber ($\kappa\ll m$)

Recall that our algorithm is independent of the wavenumber because we
integrate along Gustafsson’s contours, which are the steepest descent contours
with respect to the spherical wave component (see Section 3.1). When the
Fourier mode $m$ is larger than the scaled wavenumber $\kappa$, it is more
efficient to integrate along a different contour. If instead, we choose the
steepest descent contour on which $\exp{(im\phi)}$ does not oscillate, we
arrive at an alternative algorithm whose cost is $O(\kappa)$ and independent
of $m$ . When $\kappa$ is extremely small, this algorithm is essentially
$O(1)$. The case where $\beta_{-}$ is small (i.e., when the source and target
are close) is handled in an identical fashion to the method described in
Section 4.2. Thus, this alternative algorithm’s cost is completely independent
of both $m$ and $\beta_{-}$, and grows as $O(\kappa)$.

###  6.2 An $O(1)$ Evaluator of the Modal Green’s Functions for the Laplace
Equation

The same method described in Section 6.1 can be applied to the case where
$\kappa=0$ to yield an $O(1)$ evaluator of the modal Green’s function for the
Laplace equation, whose cost is independent of $\beta_{-}$ (i.e., the cost is
independent of the source-to-target distance).

###  6.3 Extension of the Algorithm to Complex $\kappa$

In this manuscript, we assumed $\kappa\in\mathbb{R}$ and $\kappa>0$, where
$\kappa$ is the scaled wavenumber. When the scaled wavenumber $\kappa$ is
complex (i.e., when the medium is attenuating), Gustafsson’s steepest descent
contours are rotated in the complex plane. The same algorithm described in
this manuscript applies in this case, with the only modification being a
change in the geometry of the steepest descent contours and the locations of
the intersection points of the contours with the Bernstein ellipse.

###  6.4 Extension to an $O(m)$ Evaluator for a Collection of Modal Green’s
Functions, with Amortized Cost $O(1)$

This manuscript presents an algorithm for the evaluation of a single modal
Green’s function for the Helmholtz Equation in $O(m)$ time, independent of
$\beta_{-}$ and $\kappa$, where $\beta_{-}$ is the scaled minimum source-to-
target distance and $\kappa$ is the scaled wavenumber. It is possible to use
this algorithm to compute all of the modal Green’s functions
$-M,-M+1,\dots,M-1,M$ in $O(M)$ time using the following method. In [17],
Matviyenko presents a five-term recurrence relation for the modal Green’s
functions for the Helmholtz equation. He observes that the recurrence relation
is stable upwards for one range of Fourier modes and stable downwards for
another range of modes. Furthermore, there exists a range of modes for which
the recurrence is bi-unstable. Thus, a classical Miller-type algorithm cannot
be applied. However, it was recently observed in [18] that if a recurrence
relation is represented as a banded matrix, then the inverse power method can
be used to find a solution, even when the stability behavior is mixed in the
sense just described. We thus apply the inverse power method, as described in
[18], to the resulting five-diagonal matrix corresponding to Matviyenko’s
recurrence relation. In this fashion, we obtain all the eigenvectors
corresponding to the zero eigenvalue; only one vector in this eigenspace
corresponds to the vector of modal Green’s functions. We thus use the $O(m)$
evaluator of this manuscript to select the vector corresponding to the modal
Green’s functions. The cost of performing the inverse-power method is $O(M)$,
and the cost of the evaluation of the $M$th modal Green’s function is $O(M)$,
meaning that all $M$ Fourier coefficients are obtained in $O(M)$ time.

## References

  * [1] Abdelmageed, Alaa K. “Efficient evaluation of modal Green’s function arising in EM scattering by bodies of revolution.” Pr. Electromag. Res. S. 27 (2000): 337–356.
  * [2] Abramowitz, Milton and Irene A. Stegun. Handbook of Mathematical Functions. National Bureau of Standards, 1964.
  * [3] Andreasen, M. “Scattering from bodies of revolution.” IEEE. T. Antenn. Propag. 13.2 (1965): 303–310.
  * [4] Bremer, James. “An algorithm for the numerical evaluation of the associated Legendre functions that runs in time independent of degree and order.” J. Comput. Phys. 360 (2018): 15-38.
  * [5] Bremer J., Z. Gimbutas, and V. Rokhlin. “A nonlinear optimization prcocedure for generalized Gaussian quadratures.” SIAM J. Sci. Comput. 32.4 (2010): 1761–1788.
  * [6] Cheng, Hongwei, W.Y. Crutchfield, Z. Gimbutas, L.F. Greengard, J.F. Ethridge, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao. “A wideband fast multipole method for the Helmholtz equation in three dimensions.” J. Comput. Phys. 216.1 (2006): 300–325.
  * [7] Cohl, H. and J. Tohline. “A Compact Cylindrical Green’s Function expansion for the Solution of Potential Problems.” Astrophys. J. 527.1 (1999): 86.
  * [8] Conway, J. and H.S. Cohl. “Exact Fourier expansion in cylindrical coordinates for the three-dimensional Helmholtz Green function.” Z. Angew. Math. Phys. 61.3 (2010): 425–443.
  * [9] Epstein, C., L. Greengard, and M. O’Neil. “A high-order wideband direct solver for electromagnetic scattering from bodies of revolution.” J. Comput. Phys. 387 (2019): 205–229.
  * [10] Gedney, S. and R. Mittra. “The use of the FFT for the efficient solution of the problem of electromagnetic scattering by a obdy of revolution.” IEEE. T. Antenn. Propag. (1988): 92–95.
  * [11] Gustafsson, Mats. “Accurate and efficient evaluation of modal Green’s functions.” J. of Electromagnet. Waves. 24.10 (2010): 1291–1301.
  * [12] Helsing, J. and A. Holst. “Variants of an explicit kernel-split panel based Nystrom discretization scheme for Helmholtz boundary value problems.” Adv. Comput. Math. 41.3 (2015): 691–708.
  * [13] Helsing, J. and A. Karlsson. “An explicit kernel-split panel-based Nystrom scheme for integral equations on axially symmetric surfaces.” J. Comput. Phys. 272 (2014): 686–703.
  * [14] Lai, J. and M. O’Neil. “An FFT-accelerated direct solver for electromagnetic scaterring from penetrable axisymmetric objects.” J. Comput. Phys. 390 (2019): 152–174.
  * [15] Liu, Y. and A. Barnett. “Efficient numerical solution of acoustic scattering from doubly-periodic arrays of axisymmetric objects.” J. Comput. Phys. 324 (2016): 226–245.
  * [16] Mason, J. Chebyshev polynomials. CRC Press, 2002
  * [17] Matviyenko, Gregory. “On the azimuthal Fourier components of the Green’s function for the Helmholtz equation in three dimensions.” J. Math. Phys. 36.9 (1995): 5159–5169.
  * [18] Osipov, Andrei. “Evaluation of small elements of the eigenvectors of certain symmetric tridiagonal matrices with high relative accuracy.” Appl. Comput. Harmon. A. 43.2 (2017): 173–211
  * [19] Trefethen, N. Approximation Theory and Practice. SIAM, 2019
  * [20] Trefethen, N. Spectral methods in MATLAB. SIAM, 2000
  * [21] Vaessen, Jean-Pierre A., and M. van Beurden. “Accurate and efficient computation of the modal Green’s function arising in the electric-field integral equations for a body of revolution.” IEEE T. Antenn. Propag. 60.7 (2012): 3294–3304.
  * [22] Wang, Peng and G. Xiao. “A Note on the Singularity Extraction Technique in Solving Scattering Problems for Bodies of Revolution.” Asia Pacif. Microwave. (2010): 2146–2148.
  * [23] Young, P., S. Hao, and P.G. Martinsson. “A high-order Nystrom discretization scheme for boundary integral equations defined on rotationally symmetric surfaces.” J. Comput. Phys. 40.1 (2014): 4142–4159.


# Training Energy-Based Models with Diffusion Contrastive Divergences

Weijian Luo1, Hao Jiang3, Tianyang Hu2, Jiacheng Sun2, Zhenguo Li2, Zhihua
Zhang1
1Peking University, 2Huawei Noah’s Ark Lab, 3Harbin Institute of Technology
(Shenzhen) This work was done when he was a research intern at Huawei Noah’s
Ark Lab. Email<EMAIL_ADDRESS>work was done when he was a
research intern at Huawei Noah’s Ark LabCorresponding to: Jiacheng Sun
<EMAIL_ADDRESS>
###### Abstract

Energy-Based Models (EBMs) have been widely used for generative modeling.
Contrastive Divergence (CD), a prevailing training objective for EBMs,
requires sampling from the EBM with Markov Chain Monte Carlo methods (MCMCs),
which leads to an irreconcilable trade-off between the computational burden
and the validity of the CD. Running MCMCs till convergence is computationally
intensive. On the other hand, short-run MCMC brings in an extra non-negligible
parameter gradient term that is difficult to handle. In this paper, we provide
a general interpretation of CD, viewing it as a special instance of our
proposed Diffusion Contrastive Divergence (DCD) family. By replacing the
Langevin dynamic used in CD with other EBM-parameter-free diffusion processes,
we propose a more efficient divergence. We show that the proposed DCDs are
both more computationally efficient than the CD and are not limited to a non-
negligible gradient term. We conduct intensive experiments, including both
synthesis data modeling and high-dimensional image denoising and generation,
to show the advantages of the proposed DCDs. On the synthetic data learning
and image denoising experiments, our proposed DCD outperforms CD by a large
margin. In image generation experiments, the proposed DCD is capable of
training an energy-based model for generating the Celab-A $32\times 32$
dataset, which is comparable to existing EBMs.

##  1 Introduction

Energy-Based Models (EBMs) are an important part of unsupervised learning [1,
2, 3]. Paired with the superb expressive power of deep neural networks, EBMs
draw great attention in the machine learning community and have broad
applications in many unsupervised learning tasks such as generative modeling
[4, 5, 6, 7, 8, 9], out-of-distribution detection [10, 11, 12], concept
learning [13, 14] and others [15, 16, 17, 18]. Despite the popularity, the
training of EBMs is challenging and remains an active field of research. One
dominant line of training methods of EBMs relies on sampling from the EBMs by
running MCMC chains [19, 20, 2, 8, 21, 5, 9], whose convergence can be
computationally expensive in practice. To improve efficiency, [20] proposed
the Contrastive Divergence (CD), which was calculated via _short-run_ MCMC
chains that are initialized from data samples. An overview of CD can be seen
in Figure 1(a), where the data distribution is transported with EBM-induced
MCMCs as the upper line of the figure illustrates. The CD was further
developed in many works [22, 23, 8, 6, 24, 8, 9] and has become a general
approach for training EBMs.

Nonetheless, the CD has its own drawbacks that are deeply rooted in the
employed MCMC mechanism. To be more specific, samples from MCMCs are induced
by EBMs, so these samples depend on EBMs’ parameters, leading to a non-
negligible gradient term that is difficult to handle as we introduced in
Section 2. Some works overlooked the parameter dependence for simplicity [20,
25]. As pointed out by Du et al. [21], such an omission leads to training
failures, e.g., non-convergence of training objectives. To address the
parameter-dependence issue, Du et al. [21] proposed to consider the non-
negligible gradient term through an additional non-parametric entropy
estimation component. However, the non-parametric entropy estimation is
neither efficient nor scalable for high-dimensional data.

(a) Contrastive Divergence

(b) Diffusion Contrastive Divergence

Figure 1: Illustration of DCD and CD. The yellow area represents the
corresponding divergence. The CD takes the EBM-induced Langevin dynamics to
transport data and EBM distribution to meet with the same EBM distribution.
The DCD considers a more general diffusion process to transport both data and
EBM distribution to meet with the same distribution.

In this work, we address the parameter-dependence issue of CD by extending the
Langevin diffusion, a commonly used MCMC for CD, to general diffusion
processes and propose a novel family of divergences — the _diffusion
contrastive divergence_ (DCD) family, as illustrated in Figure 1(b). Our
proposed DCD family is both theoretically sound and computationally efficient.
The contributions of our proposed DCD is three folded. First, the DCD
overcomes the non-negligible gradient issue of CD that influence the accuracy
of CD. Second, the DCD does not depend on EBM-induced MCMC so is efficient
when implemented. Third, the proposed DCD framework provides a unified view
that includes the CD as a special instance. The framework can potentially
benefit further understanding and developing algorithms for training EBMs. To
demonstrate the effectiveness and efficiency of the proposed DCD, we
instantiate the DCD with a special VE diffusion process and call it the DCD-VE
(or just DCD for short) algorithm. We conduct experiments with DCD-VE in three
experiments including synthetic data modeling, image denoising, and image
generation. On the synthetic data learning and high-dimensional image
denoising experiments, the proposed DCD-VE outperforms CD with a significant
margin. On the image generation experiment, we train a time-dependent energy-
based model on the CelebA dataset of a resolution of $32\times 32$. The
trained EBM is comparable to previous EBMs on generation. Besides, the
experiments demonstrate that the DCD is more efficient than CD, being 2-4
times faster in terms of the wall-clock time.

##  2 Background

#### Energy-based models.

Let $p_{d}$ represent the data distribution. An energy-based model specifies
the density with a neural-parametrized energy function $f_{\theta}(\bm{x})$
with the form

$\displaystyle
p_{\theta}(\bm{x})=\frac{\exp(f_{\theta}(\bm{x}))}{Z_{\theta}},$ (1)

where $f_{\theta}$ is usually a deep neural network and
$Z_{\theta}=\int\exp(f_{\theta}(\bm{u}))d\bm{u}$ is the unknown normalizing
constant. To make the derivation neat, we slightly abuse the conventions and
call $f_{\theta}(\bm{x})$ the energy function. In most cases, $Z_{\theta}$ is
so complicated that is intractable, making the likelihood intractable as well.
Previous works find out that the difficulty of estimating the normalizing
constant can be circumvented with consistent sampling from the EBM when
training with Maximum Likelihood Estimation (MLE). More precisely, the
derivative of EBM’s expected likelihood over data distribution has an
expression

$\displaystyle\frac{\partial}{\partial\theta}\mathbb{E}_{p_{d}}\log\frac{\exp(f_{\theta}(\bm{x}))}{Z_{\theta}}=\mathbb{E}_{p_{d}}\frac{\partial}{\partial\theta}f_{\theta}(\bm{x})-\mathbb{E}_{p_{\theta}}\frac{\partial}{\partial\theta}f_{\theta}(\bm{x}).$
(2)

This expression shows that the likelihood function’s parameter gradient can be
estimated with samples consistently drawn from data and the EBM. The Langevin
dynamics (LD), is a usual choice of MCMC for obtaining samples from EBMs. It
simulates the diffusion process

$\displaystyle\mathrm{d}\bm{x}_{t}=\frac{1}{2}\nabla_{\bm{x}_{t}}\log
p_{\theta}(\bm{x}_{t})\mathrm{d}t+\mathrm{d}\bm{w}_{t},$ (3)

in order to draw samples from the EBM. Under mild conditions [26], the
marginal distribution of equation 3 will converge to the target distribution
regardless of the initial distribution. Here $\bm{w}_{t}$ is an independent
Wiener process. Notice that the normalizing constant $Z_{\theta}$ in Equation
equation 1 is independent of $\bm{x}$, so we have

$\displaystyle\nabla_{\bm{x}_{t}}\log
p_{\theta}(\bm{x}_{t})\coloneqq\nabla_{\bm{x}_{t}}\big{[}f_{\theta}(\bm{x}_{t})+Z_{\theta}\big{]}=\nabla_{\bm{x}_{t}}\log
f_{\theta}(\bm{x}_{t}).$

This shows that the LD can take EBMs’ neural network without the influence of
the unknown normalizing constant.

#### Contrastive divergence and the non-negligible gradient term.

Training EBMs with MLE requires MCMC chains to run sufficiently long so as to
draw samples from the EBM. Some works studied the possibility of training EBMs
with un-converged MCMCs. Hinton [20] and Hinton et al. [2] observed that a few
MCMC steps which are initialized from data samples work well empirically so
they argued the MCMC chains do not need to fully converge when training EBMs.
They thus formally proposed the Contrastive Divergence as

$\displaystyle\mathcal{D}_{CD}(p_{d},p_{\theta})=\mathcal{D}_{KL}(p_{d},p_{\theta})-\mathcal{D}_{KL}(p_{d,\theta}^{(T)},p_{\theta}),$
(4)

where $p_{d,\theta}^{(T)}$ stands for the marginal distribution of a short-run
MCMC initialized from $p_{d}$ with transition time $T$ and the notation
$\mathcal{D}_{KL}$ denotes the Kullback–Leibler (KL) divergence. For such a
definition, the non-negativity $\mathcal{D}_{CD}(p,q)\geq 0$ holds and
$\mathcal{D}_{CD}(p_{d},p_{\theta})=0$ only when $p_{t}=q_{t}$ almost
everywhere. This makes $\mathcal{D}_{CD}$ a reasonable divergence, we put
detailed derivation on the non-negativity of CD in the Appendix. If we take
the parameter derivative, we have

$\displaystyle\frac{\partial}{\partial\theta}\mathcal{D}_{CD}(p_{d},p_{\theta})=\mathbb{E}_{p_{d,\theta}^{(T)}}\big{[}\frac{\partial}{\partial\theta}f_{\theta}(\bm{x})\big{]}-\mathbb{E}_{p_{d}}\big{[}\frac{\partial}{\partial\theta}f_{\theta}(\bm{x})\big{]}-\mathbb{E}_{p_{d,\theta}^{(T)}}\bigg{[}\log
p_{\theta}(\bm{x})\frac{\partial}{\partial\theta}\log
p_{d,\theta}^{(T)}(\bm{x})\bigg{]}.$ (5)

The third gradient term is difficult to handle because, for EBM, we do not
know the value of normalizing constant $Z_{\theta}$. So the gradient
$\frac{\partial}{\partial\theta}\log p_{d,\theta}^{(T)}(\bm{x})$ is also
unknown. Hinton [20] and Liu and Wang [25] proposed to omit the third gradient
term and simplify the CD equation 6 as

$\mathbb{E}_{\bm{x}_{T}\sim\operatorname{sg}[p_{d,\theta}^{(T)}]}f_{\theta}(\bm{x}_{T})-\mathbb{E}_{\bm{x}\sim
p_{d}}f_{\theta}(\bm{x}).$ (6)

Here the notation $\bm{x}_{T}\sim\operatorname{sg}[p_{d,\theta}^{(T)}]$
represents the sample $x_{T}$ is drawn from $p_{d,\theta}^{(T)}$ but omitting
the parameter dependence of $\theta$. In practice, there is always a non-
negligible term for the gradient of the contrastive divergence. Du et al. [21]
tried to address the non-negligible term by introducing an additional non-
parametric entropy estimation component together with the training of EBM,
viewing the non-negligible third term of equation 5 as a parameter derivative
of Shannon entropy that is estimated non-parametrically. Although technically
sound, the entropy estimation which Du et al. [21] brought in is
computationally intensive and not scalable in high dimensions.

#### Diffusion process.

A diffusion process is a stochastic process driven by a stochastic
differential equation (SDE) [27] with a drift vector $\bm{F}$ and a diffusion
matrix $\bm{G}$,

$\mathrm{d}\bm{x}_{t}=\bm{F}(\bm{x}_{t},t)\mathrm{d}t+\bm{G}(t)\mathrm{d}\bm{w}_{t},$
(7)

where $\bm{w}_{t}$ is a standard Wiener process. For simplicity, we assume
$\bm{G}$ to be a scalar function of time $t$ in the rest of the paper. If a
diffusion process is initialized with an initial distribution $p_{0}$, then
the evolution of marginal probability density is governed by the Fokker-Planck
equation [28]:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(\bm{x},t)$
$\displaystyle=-\langle\nabla_{\bm{x}},p(\bm{x},t)\bm{F}(\bm{x},t)\rangle+\frac{1}{2}\bm{G}^{2}(t)\Delta_{\bm{x}}p(\bm{x},t),\
p(\bm{x},0)=p_{0}(\bm{x}).$ (8)

The Langevin dynamics defined in equation 3 is an instance of diffusion
processes. The VE diffusion is a commonly used diffusion process in generative
modeling [29, 30, 31]. It writes

$\displaystyle\mathrm{d}\bm{x}_{t}=g(t)\mathrm{d}\bm{w}_{t}.$ (9)

The diffusion has explicit conditional distributions
$p_{t}(\bm{x}_{t}|\bm{x}_{0})$ and their marginal samples are cheap to obtain
as we put in the Appendix.

##  3 Diffusion contrastive divergences

Our goal is to propose novel training methods that overcome both the non-
negligible gradient term and the inefficiency issue caused by MCMC of CD, by
generalizing the definition of CD to other parameter-free diffusion processes,
named diffusion contrastive divergence (DCD). In this section, we first give
the formal definition of DCD. Then we establish the connections of DCD to
existing methods, namely the diffusion recovery likelihood and the KL-
contraction divergence. Later we proposed a practical algorithm, the DCD-VE
based on the VE diffusion equation 9 for training energy-based models.

###  3.1 CD with general diffusions

We follow the notations defined in Section 2 and take the LD as the MCMC which
defines the CD. Recall the definition of CD equation 6.

One of the most important reasons for taking LD to define the divergence is
that the KL divergence of the marginal distributions with LD is strictly
decreasing and converges to $0$ when $T\to\infty$ unless $p_{\theta}=p_{d}$
(We put in Appendix). This makes the CD a well-defined divergence. But the
definition of LD equation 3 incorporates the EBM and its parameters $\theta$,
giving rise to a hard-to-handle non-negligible gradient term as we pointed out
in 2. Besides, obtaining samples with LD also relies on the sequential
simulation of SDE which is computationally inefficient. So it would be ideal
if the Langevin dynamics that the CD uses are replaced with some parameter-
free alternatives.

Fortunately, other diffusion processes, such as the VE process equation 9 with
the properly defined function $g(t)$ also guarantee the strict decrease and
the convergence of the KL between marginal distributions as LD does. Besides,
the definition of such diffusion processes does not contain any EBM
parameters, and the marginal samples are efficient to obtain as we put in
discussions in Appendix.

Based on such an observation, we formally define the _Diffusion Contrastive
Divergence_ (DCD), as the KL difference between an initial distribution and
the transitional distribution under some pre-defined diffusion process.

######  Definition 1 (Diffusion Contrastive Divergence).

$\displaystyle\mathcal{D}_{DCD}^{(\bm{F},\bm{G},T)}(p_{d},p_{\theta}):=\mathcal{D}_{KL}(p_{d},p_{\theta})-\mathcal{D}_{KL}(p_{d}^{(T)},p_{\theta}^{(T)}).$
(10)

Here $p_{d}^{(T)}$ and $p_{\theta}^{(T)}$ stand for the marginal distributions
of the diffusion equation 7 that are initialized with $p_{d}$ and $p_{\theta}$
respectively.

To further study the properties of the proposed DCDs, we first give a theorem
to verify that the DCD is a well-defined probability divergence.

######  Theorem 2.

Let $\bm{F}(\bm{x},t)$ and $\bm{G}(t)$ be two pre-defined functions. For two
distributions $p$ and $q$, assume both $p,q$ evolve according to the same
diffusion process equation 7. Let $p^{(t)}$ and $q^{(t)}$ denote the time $t$
marginal distribution under SDE evolution. Then we have

$\displaystyle\mathcal{D}_{DCD}^{(\bm{F},\bm{G},T)}(p,q)=\frac{1}{2}\int_{0}^{T}\mathbb{E}_{\bm{x}_{t}\sim
p^{(\bm{F},\bm{G},t)}(x)}\bm{G}^{2}(t)\|\nabla_{\bm{x}_{t}}\log
p^{(\bm{F},\bm{G},t)}(\bm{x}_{t})-\nabla_{\bm{x}_{t}}\log
q^{(\bm{F},\bm{G},t)}(\bm{x}_{t})\|^{2}_{2}\mathrm{d}t.$

We give detailed proof in the Appendix. From Proposition 1, we see that DCD is
non-negative.

######  Proposition 1.

For any two distributions $p$ and $q$, any function $\bm{F},\bm{G}$ and any
diffusion time $T$, then

$\mathcal{D}_{DCD}^{(\bm{F},\bm{G},T)}(p,q)\geq 0.$

With a suitable choice of $\bm{F}$ and $\bm{G}$, the KL divergence between
marginal distributions is strictly decreasing, thus the defined
$\mathcal{D}^{(\bm{F},\bm{G},T)}$ does not degenerate, making
$\mathcal{D}^{(\bm{F},\bm{G},T)}(p,q)=0$ if and only if $p=q,a.e.$. As we show
in Appendix, the VE diffusion satisfies this property.

For a diffusion process that does not depends on EBM’s parameter, the
corresponding DCD avoids the parameter-dependence issues. Figure 1(b) gives
the concept of DCD. Both the data and EBM’s distribution evolve along the
diffusion process specified by $(\bm{F},\bm{G})$ as in equation 7. With
suitable choices, when $T\to\infty$, two involved distributions coincide with
the same stationary distribution. The yellow region accounts for what DCD
measures.

######  Remark 1.

Notice that $p_{\theta}$ itself is a stationary distribution of the above LD
as we put in Appendix. Hence, $p_{\theta}^{(t)}=p_{\theta}$ holds for any
$t\in[0,T]$. So if we choose a special
$\bm{F}(\bm{x},t)=\nabla_{\bm{x}}f_{\theta}(\bm{x})/2$ and
$\bm{G}(t)=\mathbf{I}$, the proposed $\mathcal{D}_{DCD}^{(\bm{F},\bm{G},T)}$
recovers CD (equation 6).

Table 1: Comparison of DCD and CD. Method | MCMC | Diffusion Process | One-step DCD Formula
---|---|---|---
CD | $\checkmark$ | $\mathrm{d}\bm{x}_{t}=-\nabla\frac{f_{\theta}(\bm{x}_{t})}{2}\mathrm{d}t+\mathrm{d}\bm{w}_{t}$ | Stationary
DCD-VE | ✗ | $\mathrm{d}\bm{x}_{t}=g(t)\mathrm{d}\bm{w}_{t}$ |  Eq.(12)

To be more concrete, we consider VE diffusion as a demonstration. Recall the
definition of VE diffusion 9. The conditional distribution of the VE diffusion
does not depend on EBM’s parameter $\theta$. The marginal samples can be drawn
with $\bm{x}_{0}\sim p_{d},\bm{x}_{t}\sim p_{t}(\bm{x}_{t}|\bm{x}_{0})$.

######  Theorem 3.

Minimizing the DCD is equivalent to minimizing the following divergence.

$\displaystyle\mathcal{L}_{DCD}(\theta)=$
$\displaystyle\mathbb{E}_{\bm{x}_{0}\sim p_{d},\bm{x}_{t}\sim
p(\bm{x}_{t}|\bm{x}_{0})}\big{[}f_{\theta}^{(\bm{F},\bm{G},T)}(\bm{x}_{t})\big{]}-\mathbb{E}_{\bm{x}_{0}\sim
p_{d}}\big{[}f_{\theta}(\bm{x}_{0})\big{]}.$ (11)

Here $f_{\theta}^{(\bm{F},\bm{G},T)}$ are time $T$ marginal energy under
diffusion process (equation 7).

Check the Appendix for detailed proof. The term $\log p_{d}(\bm{x})$ and $\log
p_{d}^{(\bm{F},\bm{G},T)}$ are independent of parameter $\theta$ since the
diffusion process is parameter-free. The equation equation 11 defines a
tractable objective that is equivalent to DCD.

The advantages of the DCD with VE diffusion over CD are two-fold. First,
recall that the CD is hindered by the parameter-dependence of both the
transitional distribution $\bm{x}_{T}\sim p_{d,\theta}^{(T)}$ and the $T$-time
evolved data distribution $\log p_{d,\theta}^{(T)}(\bm{x}_{T})$ in the MCMC
chains. These two terms are parameter-free if we choose a parameter-free
diffusion instead of Langevin dynamics.

Second, the sampling from VE diffusion gets significantly cheaper when taking
specially designed diffusions such as VE diffusion. However as a trade-off,
one needs to evaluate the time $T$ marginal energy of
$f_{\theta}^{(\bm{F},\bm{G},T)}(\bm{x}_{t})$, which can be easier to handle.
We provide further analysis of the energy evolution in Section 3.3. To
summarize, DCD is an MCMC-free method that overcomes the CD’s two difficulties
with one easier problem of estimating the energy evolution. Such MCMC-free
training methods for EBMs are a hot research area in EBM community[32]. We
give a brief summary of the differences between CD and the DCD that are
defined through the VE diffusion in Table 1.

###  3.2 Connections to existing methods

The DCD framework not only provides a new understanding of CD but also more
insights into existing works on training EBMs. For instance, DCD has inner
connections to two existing methods, the Diffusion Recovery Likelihood [5, 33]
and the KL-Contraction Divergence [34].

#### Connection to Diffusion Recovery Likelihood.

Let
$p^{(\sigma)}(\tilde{\bm{x}}|\bm{x})=\mathcal{N}(\tilde{\bm{x}};\bm{x},\sigma^{2}\mathbf{I})$
denotes a Gaussian perturbation on $x$. The recovery likelihood of a data
$\bm{x}$ is defined as the conditional probability to recover $\bm{x}$ from
noise perturbed observation $\tilde{\bm{x}}$, i.e.,
$p_{\theta}(\bm{x}|\tilde{\bm{x}})=p^{(\sigma)}(\tilde{\bm{x}}|\bm{x})p_{\theta}(\bm{x})/p_{\theta}(\tilde{\bm{x}})$,
which is proportional to
$\exp(f_{\theta}(\bm{x})-\frac{1}{2\sigma^{2}}\|\tilde{\bm{x}}-\bm{x}\|_{2}^{2})$.

Gao et al. [5] viewed recovery likelihood as a new EBM for $\bm{x}$ if
$\tilde{\bm{x}}$ is given as fixed and minimized the recovery likelihood
through a CD-like MCMC method for which negative samples are consistently
sampled from $p_{\theta}(\bm{x}|\tilde{\bm{x}})$-induced MCMC. Gao et al. [5]
also extended the recovery likelihood to multi-level Gaussian noise level
$\\{\sigma_{i}\\}$ to define a diffusion recovery likelihood. Surprisingly as
we show in this section, the recovery likelihood objective is a special case
of DCD when taking the diffusion process to be the VE diffusion. Revisit that
the definition of the recovery likelihood writes

$\mathbb{E}_{\bm{x}\sim p_{d},\tilde{\bm{x}}\sim
p_{d}^{(\sigma)}(\tilde{\bm{x}})}\log p_{\theta}(\bm{x}|\tilde{\bm{x}}),$

the $p^{(\sigma)}(\tilde{\bm{x}}|\bm{x})$ and $p_{d}(\bm{x})$ are independent
of parameter $\theta$, so maximizing the recovery likelihood is equivalent to
minimizing

$\displaystyle\mathcal{D}_{KL}(p_{d}(\bm{x}),p_{\theta}(\bm{x}))-\mathcal{D}_{KL}(p_{d}^{(\sigma)}(\tilde{\bm{x}}),p_{\theta}^{(\sigma)}(\tilde{\bm{x}})).$

We put the detailed derivation in the Appendix. Here
$p_{\theta}^{(\sigma)}(\tilde{\bm{x}})=\int
p_{\theta}(\bm{x})p(\tilde{\bm{x}}|\bm{x})d\bm{x}$ is the marginal density of
Gaussian perturbed distribution. The recovery likelihood and its diffusion
counterpart are special cases of DCD when taking the diffusion process to be
VE diffusion equation 9. When setting $\sigma_{i}^{2}=\int_{0}^{t_{i}}g(s)ds$,
the DCD-VE recovers the diffusion recovery likelihood. However, the
implementation of maximizing recovery likelihood in [5] is different. They
sample from $\log p_{\theta}(\bm{x}|\tilde{\bm{x}})$ through MCMC when
training, making the training procedure computationally expensive. In our
definition of the DCD, we do not require sampling from recovery likelihood. We
instead use contrastive mechanics between $p$ and $p^{(T)}$ to cancel out the
normalizing constant as we introduced in later sections. Besides, the DCD
framework can be generalized to other diffusion processes of which the
definition does not involve EBM’s parameters.

#### DCD as a KL-contraction divergence.

[34] proposed the so-called KL contraction divergence framework. They pointed
out that if an operator $\Phi(p)$ satisfies the KL contraction property,
meaning

$\mathcal{D}_{KL}(\Phi(p),\Phi(q))\leq\mathcal{D}_{KL}(p,q),$

a KL-contraction divergence can be defined as
$\mathcal{D}_{KL}(p,q)-\mathcal{D}_{KL}(\Phi(p),\Phi(q))$. As we mentioned in
the Theorem 2, the marginalization along any diffusion process is a KL
contraction operator, so the DCD can be viewed also as a KL-contraction
divergence. However, in our paper, we define the DCD through the motivation of
generalizing the CD. Besides, we propose a concrete divergence, the DCD-VE,
which is much different from the instances that have been studied in Lyu [34].

###  3.3 Evolution of the energy function

Since the definition of DCD equation 10 involves the computation of the
diffused density function $p_{d,\theta}^{(T)}(\bm{x})$ and corresponding
energy function $f_{\theta}^{(T)}(\bm{x})$, so in this section, we
characterize the evolution of the energy function $f_{\theta}^{(T)}(\bm{x})$
through a partial differential equation. Denote
$p_{\theta}^{(0)}(\bm{x})=e^{f_{\theta}(\bm{x})}/Z_{\theta}$ where
$Z_{\theta}$ is the normalizing constant. We show that the evolution of the
energy function under the diffusion process (equation 7) follows a PDE.

######  Proposition 2.

Assume $p_{\theta}^{(0)}(\bm{x})=e^{f_{\theta}(\bm{x})}/Z_{\theta}$ where
$Z_{\theta}$ is a parameter-dependent normalizing constant. Assume
$p_{\theta}^{(t)}$ denotes the evolved density along a diffusion process
equation 7, then for any fixed $\bm{x}$, the energy value
$p_{\theta}^{(t)}(\bm{x})$ evolves according to a PDE

$\mathrm{d}\log
p_{\theta}^{(t)}(\bm{x})/\mathrm{d}t=\mathcal{O}(\nabla_{\bm{x}}\log
p_{\theta}^{(t)}),$

where $\mathcal{O}(\nabla_{\bm{x}}\log p_{\theta}^{(t)})$ is the following
operator which is independent of the normalizing constant,

$\displaystyle\langle\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t),\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\rangle+\langle\nabla,\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t)\rangle.$

It is worth emphasizing that since the evolution operator $\mathcal{O}(.)$
does not depend on $Z_{\theta}$, the normalizing constant keeps unchanged in
the process and thus will be exactly canceled out when we substitute the
$T$-time KL and initial KL as in DCD expression. So the DCD is not bothered by
a parameter-dependent normalizing constant. We give a more detailed argument
in the Appendix.

In practice, we do not need many steps when training EBM. So we use a single
step as an approximation when implementing DCD. Our experiments show that the
single-step DCD works well in practice. Here we derive a one-step
approximation of DCD for practical implementations.

#### DCD-VE.

For VE diffusion equation 9, the time-change rate of energy can be
approximated with

$\displaystyle\mathcal{L}_{DCD}^{(VE)}(\theta)=\mathbb{E}_{p_{t}}\frac{1}{2}\bm{G}^{2}(0)\bigg{[}\|\nabla_{\bm{x}}f_{\theta}(\bm{x}_{t})\|^{2}+\Delta
f_{\theta}(\bm{x}_{t})\bigg{]}+\frac{1}{t}\bigg{[}\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})]\bigg{]}.$
(12)

The detailed derivations are put in Appendix. We formally define the DCD-VE
objective for training EBM as $\mathcal{L}_{DCD}^{(VE)}(\theta)$ in (12) with
a small perturbation level $t$. For one-step
$\mathcal{L}_{DCD}^{(VE)}(\theta)$, if data is low dimensional, the second
order derivative is computationally tractable. However, for high dimensional
data such as natural images, the second order derivative (the Laplacian term)
can be efficiently estimated by the widely-used Hutchinson’s trace estimation
techniques [35, 36, 37, 38, 39].

###  3.4 Train time-dependent EBM with DCD

Inspired by recent success on score-based diffusion models [5, 40, 29, 30],
learning a diffusion time-dependent EBM helps for better generative
performance. In this section, we modify our DCD-VE for training time-dependent
EBMs. Assuming $(\bm{F},\bm{G})$ denotes a pre-defined forward diffusion
process equation 7 (as we use when defining DCD). Let $p_{d}^{(0)}$ denotes
the data distribution, and $p_{d}^{(t)}$ denotes the $t$-time diffused data
distribution initialized with $p_{d}^{(0)}$. A time-dependent EBM is a
$f_{\theta}^{(t)}$ if a neural network that takes both $\bm{x}$ and time $t$
to output the energy function of a point $\bm{x}$ at diffusion time $t$. One
can train $f_{\theta}^{(t)}$ to model the diffused data energy $\log
p_{d}^{(t)}(\bm{x})$ at any time $t$. More precisely, at each training
iteration, we randomly pick a timestamp $t\sim Unif([0,T])$, and apply DCD
training at timestamp $t$ with a small diffusion perturbation $\delta$. In
practice, if we discretize the time interval of a diffusion process to
$\\{t_{i}\\}_{i=1,..,K}$, the perturbation $\delta$ can be chosen to be
$\delta_{i}=t_{i}-t_{i-1}$ for different time $t_{i}$. Such a setting combines
the DCD and diffusion process in a more natural way. We summarize the DCD
training for time-dependent EBM in an Algorithm in the Appendix.

##  4 Experiments

###  4.1 Energy modeling of 2D distributions

In this section, we validate our proposed DCD on 7 commonly used 2D synthetic
datasets. This experiment shows that DCD is capable of learning challenging
distributions such as the Checkerboard distribution whose distribution changes
rapidly (as shown in the left part of Figure 2).

#### Experiment Setting.

We use a 3-layer MLP with Gaussian Error Linear Unit (GELU) activations [41]
and 300 hidden units for implementation of the EBM. We compare the DCD-VE with
CD and Persistent Contrastive Divergence (PCD)[22], which is a well-known
variant of CD. Since the CD training requires many iterations of inference of
the EBM, we limit the times of score function evaluation to 10 times to make
an equal comparison. We set the training batch size to be 1000 and PCD’s
replay buffer size to be 10 times the batch size. All models share the same
architecture and the same training setting. We put detailed settings in
Appendix.

#### Evaluation metric.

We compute the score-matching loss over the training data as the evaluation
metric. The score matching loss is defined with

$\displaystyle\operatorname{L}(\theta)\coloneqq\mathbb{E}_{\bm{x}\sim
p_{d}}\bigg{[}\frac{1}{2}\|\nabla_{\bm{x}}f_{\theta}(\bm{x})\|_{2}^{2}+\Delta_{\bm{x}}f_{\theta}(\bm{x})\bigg{]}.$

So the smaller the SM loss is, the better the learning performance of the EBM.

Table 2: Estimated SM loss of learned EBM. Dataset | Swissroll | Circles | Rings | Moons | 8 Gaussians | 2 Spirals | Checkerboard
---|---|---|---|---|---|---|---
DCD-VE | -2398.81 | -131.37 | -758.33 | -200.67 | -120.09 | -470.92 | -178.43
CD | $+\infty$ | -130.03 | $+\infty$ | -195.72 | -117.29 | $+\infty$ | -67.22
PCD | $+\infty$ | -108.54 | $+\infty$ | -193.76 | -97.59 | $+\infty$ | -124.27

(a) Comparison of CD, PCD and DCD-VE

(b) Generated CelebA $32$ samples from EBM.

Figure 2: Left: 2D examples when CD and PCD fails to learn a correct EBM but
DCD-VE can learn successfully; Right: Generated CelebA $32$ samples from EBM
trained with DCD-VE.

#### Performance.

We estimate the Score Matching (SM) loss ([38, 42]) on training data to
evaluate the trained EBM. The smaller the SM loss, the better performance the
EBM behaves. Table 2 shows the resulting SM losses for EBMs that are trained
with DCD-VE, CD, and PCD. Since the SM loss is the training objective of SM-
related training methods, we do not include them in the comparison. As is
shown in Table 2, DCD-VE outperforms CD and PCD on all datasets by a
significant margin. Besides, CD and related methods do not converge on the
more challenging Swiss roll, Rings, and 2Spirals dataset, while the DCD-VE can
learn all data energy equally well. Figure 3 demonstrates the learned energies
on five datasets with DCD-VE.

Table 3: CelebA. Models |  FID $\downarrow$
---|---
ABP [43] | 51.50
ABP-SRI [44] | 36.84
VAE [45] | 38.76
Glow [46] | 23.32
DCGAN [47] | 12.50
EBM-FCE [48] | 12.21
GEBM [49] | 5.21
CoopFlow(T=30) [23] | 6.44
EBM-DCD | 13.85
Figure 3: Comparison of different training methods. Figure 4: The CD fails to
denoise large added noise, while the DCD (VE) can denoise successfully.

###  4.2 Image denoising with EBM

Image denoising is a common task to test explicit generative models [42].

In this section, we validate the proposed DCD for training EBM on high-
dimensional datasets and evaluate the image-denoising performance on four
datasets, the MNIST, FashionMNIST, CIFAR10, and the SVHN datasets.

#### Experiment Setting.

We train EBM with DCD-VE and compare it with CD. We added the Gaussian noise
with three strength levels on test images and evaluate the average root of the
mean of the squared error (RMSE) of non-noised and denoised images. For the
implementation of the EBM, we use the wide resnet[50] model with GELU
activations as our energy model. More details are put in Appendix.

Table 4 shows the denoising performance of Gaussian noise with different
scales of the noise (low for 0.3, middle for 0.6, and high for 0.9). The DCD-
VE performs consistently better than CD across different datasets and
different noise strengths. We also surprisingly find that one advantage of the
DCD is its impressive performance on large noise strength. As Figure 4 shows,
for a high noise scale of 0.9, the EBM trained with CD fails to denoise
successfully, while the EBM trained with DCD-VE still shows denoising ability.

Table 4: Average RMSE of clean and reconstructed input with Gaussian noise on datasets. (We set low, mid, and high-level noise as 0.3, 0.6, and 0.9.) Method | MNIST | FMNIST | CIFAFR10 | SVHN
---|---|---|---|---
low | mid | high | low | mid | high | low | mid | high | low | mid | high
DCD | 0.165 | 0.194 | 0.303 | 0.170 | 0.217 | 0.497 | 0.129 | 0.193 | 0.244 | 0.099 | 0.137 | 0.294
CD | $0.194$ | $0.390$ | $1.099$ | $0.171$ | $0.2792$ | $0.872$ | $0.154$ | $0.317$ | $8.572$ | $0.124$ | $0.293$ | $6.938$

###  4.3 Image generation with time-dependent EBM

#### Experiment Setting.

We use DCD to train EBMs for image generation on the CelebA dataset of a
resolution of $32\times 32$. We use a time-dependent neural network with
residual network architecture [51] as the implementation of the EBM. We use
the VE diffusion with the diffusion coefficient $g(t)=t$ as the forward
diffusion, which is the same as Karras et al. [52]. We train the time-
dependent energy-based model on the CelebA dataset which is downsampled to
have a resolution of 32x32. We evaluate the Frechet Inception Score (FID) [53]
as a metric of generation performance.

#### Performance.

Table 3 shows the performance of our trained EBMs for a generation. It shows
that the DCD (VE) is capable of handling complex image datasets. It
demonstrates that the proposed DCD is capable of training EBMs with comparable
performance to DCGAN [47] and other EBMs (i.e. EBM with FCE [48]), and
superior performance to normalizing flow models and VAE. However, the
performance is worse than EBM which requires more advanced tricks such as
cooperating with flow models (CoopFlow [54]) and cooperating with GAN models
(GEBM [49]). The right-hand side of Figure 2 shows some generated samples from
our trained EBM. In summary, the proposed DCD-VE is able to train time-
dependent EBM with comparable generative performance as existing training
methods.

##  5 Limitations and Future Works

In this paper, we propose a novel family of probability divergences, the
_diffusion contrastive divergence_ family. The DCD provides a special view
that unifies the contrastive divergence as a special instance of the DCD. It
also spurs new divergences for training EBM which overcomes two major
drawbacks of the contrastive divergence. We also establish the connection of
the proposed DCDs with existing recovery likelihood and the KL-contraction
divergences. We validate the efficiency and superior performance of our
proposed DCDs on several benchmark EBM tasks such as 2D energy modeling, image
denoising, and image generation.

However, the DCD also has its limitations. First, the calculation of DCD
requires the computation of a higher-order derivative of the energy function,
meaning that the energy-based model should be at least twice differentiable.
Second, the long-time DCD requires the calculation of the evolved energy
function. Such evolution is not easy to compute for the general diffusion
process. We plan to leave the research of the long-time energy evolution of
DCD in our further work.

## References

  * LeCun et al. [2006] Y. LeCun, S. Chopra, R. Hadsell, M. Ranzato, and F. Huang, “A tutorial on energy-based learning,” _Predicting structured data_ , vol. 1, no. 0, 2006\.
  * Hinton et al. [2006] G. E. Hinton, S. Osindero, and Y. W. Teh, “A fast learning algorithm for deep belief nets,” _Neural Computation_ , vol. 18, pp. 1527–1554, 2006.
  * Zhu et al. [2004] S.-C. Zhu, Y. N. Wu, and D. Mumford, “Filters, random fields and maximum entropy (frame): Towards a unified theory for texture modeling,” _International Journal of Computer Vision_ , vol. 27, pp. 107–126, 2004.
  * Xie et al. [2016] J. Xie, Y. Lu, S.-C. Zhu, and Y. Wu, “A theory of generative convnet,” in _International Conference on Machine Learning_. PMLR, 2016, pp. 2635–2644.
  * Gao et al. [2020a] R. Gao, Y. Song, B. Poole, Y. N. Wu, and D. P. Kingma, “Learning energy-based models by diffusion recovery likelihood,” _arXiv preprint arXiv:2012.08125_ , 2020.
  * Nijkamp et al. [2019] E. Nijkamp, M. Hill, S.-C. Zhu, and Y. N. Wu, “Learning non-convergent non-persistent short-run mcmc toward energy-based model,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * Zhao et al. [2021] Y. Zhao, J. Xie, and P. Li, “Learning energy-based generative models via coarse-to-fine expanding and sampling,” in _ICLR_ , 2021.
  * Du and Mordatch [2019] Y. Du and I. Mordatch, “Implicit generation and generalization in energy-based models,” _arXiv preprint arXiv:1903.08689_ , 2019.
  * Grathwohl et al. [2019a] W. Grathwohl, K.-C. Wang, J.-H. Jacobsen, D. Duvenaud, M. Norouzi, and K. Swersky, “Your classifier is secretly an energy based model and you should treat it like one,” _arXiv preprint arXiv:1912.03263_ , 2019.
  * Zhai et al. [2016] S. Zhai, Y. Cheng, W. Lu, and Z. Zhang, “Deep structured energy based models for anomaly detection,” in _International conference on machine learning_. PMLR, 2016, pp. 1100–1109.
  * Liu et al. [2020] W. Liu, X. Wang, J. Owens, and Y. Li, “Energy-based out-of-distribution detection,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 21 464–21 475, 2020.
  * Lee et al. [2020] K. Lee, H. Yang, and S.-Y. Oh, “Adversarial training on joint energy based model for robust classification and out-of-distribution detection,” _2020 20th International Conference on Control, Automation and Systems (ICCAS)_ , pp. 17–21, 2020.
  * Mordatch [2018] I. Mordatch, “Concept learning with energy-based models,” _arXiv preprint arXiv:1811.02486_ , 2018.
  * Du et al. [2020a] Y. Du, S. Li, and I. Mordatch, “Compositional visual generation with energy based models,” in _NeurIPS_ , 2020.
  * Haarnoja et al. [2017] T. Haarnoja, H. Tang, P. Abbeel, and S. Levine, “Reinforcement learning with deep energy-based policies,” in _ICML_ , 2017.
  * Xie et al. [2017] J. Xie, S.-C. Zhu, and Y. N. Wu, “Synthesizing dynamic patterns by spatial-temporal generative convnet,” _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 1061–1069, 2017.
  * Xie et al. [2018] J. Xie, Z. Zheng, R. Gao, W. Wang, S.-C. Zhu, and Y. N. Wu, “Learning descriptor networks for 3d shape synthesis and analysis,” _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 8629–8638, 2018.
  * Ingraham et al. [2019] J. Ingraham, A. J. Riesselman, C. Sander, and D. S. Marks, “Learning protein structure with a differentiable simulator,” in _ICLR_ , 2019.
  * Song and Kingma [2021] Y. Song and D. P. Kingma, “How to train your energy-based models,” _arXiv preprint arXiv:2101.03288_ , 2021.
  * Hinton [2002] G. E. Hinton, “Training products of experts by minimizing contrastive divergence,” _Neural computation_ , vol. 14, no. 8, pp. 1771–1800, 2002\.
  * Du et al. [2020b] Y. Du, S. Li, J. Tenenbaum, and I. Mordatch, “Improved contrastive divergence training of energy based models,” _arXiv preprint arXiv:2012.01316_ , 2020\.
  * Tieleman and Hinton [2009] T. Tieleman and G. Hinton, “Using fast weights to improve persistent contrastive divergence,” in _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 1033–1040.
  * Xie et al. [2022a] J. Xie, Y. Zhu, J. Li, and P. Li, “A tale of two flows: Cooperative learning of langevin flow and normalizing flow toward energy-based model,” _arXiv preprint arXiv:2205.06924_ , 2022.
  * Miyato et al. [2018] T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida, “Spectral normalization for generative adversarial networks,” _arXiv preprint arXiv:1802.05957_ , 2018\.
  * Liu and Wang [2017] Q. Liu and D. Wang, “Learning deep energy models: Contrastive divergence vs. amortized mle,” _arXiv preprint arXiv:1707.00797_ , 2017.
  * Pavliotis [2014] G. A. Pavliotis, _Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations_. Springer, 2014, vol. 60.
  * Särkkä and Solin [2019] S. Särkkä and A. Solin, “Applied stochastic differential equations,” 2019\.
  * Risken [1984] H. Risken, “Fokker-planck equation,” 1984.
  * Song et al. [2020] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” _arXiv preprint arXiv:2011.13456_ , 2020.
  * Song et al. [2021a] Y. Song, C. Durkan, I. Murray, and S. Ermon, “Maximum likelihood training of score-based diffusion models,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 1415–1428, 2021.
  * Karras et al. [2022a] T. Karras, M. Aittala, T. Aila, and S. Laine, “Elucidating the design space of diffusion-based generative models,” in _Proc. NeurIPS_ , 2022.
  * Grathwohl et al. [2021] W. S. Grathwohl, J. J. Kelly, M. Hashemi, M. Norouzi, K. Swersky, and D. Duvenaud, “No {mcmc} for me: Amortized sampling for fast and stable training of energy-based models,” in _International Conference on Learning Representations_ , 2021.
  * Bengio et al. [2013] Y. Bengio, L. Yao, G. Alain, and P. Vincent, “Generalized denoising auto-encoders as generative models,” _Advances in neural information processing systems_ , vol. 26, 2013.
  * Lyu [2011] S. Lyu, “Unifying non-maximum likelihood learning objectives with minimum kl contraction,” _Advances in Neural Information Processing Systems_ , vol. 24, 2011.
  * Hutchinson [1989] M. F. Hutchinson, “A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines,” _Communications in Statistics - Simulation and Computation_ , vol. 18, pp. 1059–1076, 1989.
  * Chen et al. [2019] R. T. Q. Chen, J. Behrmann, D. K. Duvenaud, and J.-H. Jacobsen, “Residual flows for invertible generative modeling,” _ArXiv_ , vol. abs/1906.02735, 2019.
  * Grathwohl et al. [2019b] W. Grathwohl, R. T. Q. Chen, J. Bettencourt, I. Sutskever, and D. K. Duvenaud, “Ffjord: Free-form continuous dynamics for scalable reversible generative models,” _ArXiv_ , vol. abs/1810.01367, 2019.
  * Song et al. [2019] Y. Song, S. Garg, J. Shi, and S. Ermon, “Sliced score matching: A scalable approach to density and score estimation,” in _UAI_ , 2019.
  * Song et al. [2021b] 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\. [Online]. Available: https://openreview.net/forum?id=PxTIG12RRHS
  * Song and Ermon [2020] Y. Song and S. Ermon, “Improved techniques for training score-based generative models,” _Advances in neural information processing systems_ , vol. 33, pp. 12 438–12 448, 2020.
  * Hendrycks and Gimpel [2016] D. Hendrycks and K. Gimpel, “Gaussian error linear units (gelus),” _arXiv: Learning_ , 2016.
  * Meng et al. [2020] C. Meng, L. Yu, Y. Song, J. Song, and S. Ermon, “Autoregressive score matching,” _ArXiv_ , vol. abs/2010.12810, 2020.
  * Han et al. [2017] T. Han, Y. Lu, S. Zhu, and Y. N. Wu, “Alternating back-propagation for generator network,” in _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence (AAAI)_ , San Francisco, CA, 2017, pp. 1976–1984.
  * Nijkamp et al. [2020] E. Nijkamp, B. Pang, T. Han, L. Zhou, S. Zhu, and Y. N. Wu, “Learning multi-layer latent variable model via variational optimization of short run MCMC for approximate inference,” in _Proceedings of the 16th European Conference on Computer Vision (ECCV, Part VI)_ , Glasgow, UK, 2020, pp. 361–378.
  * Kingma and Welling [2014] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” in _Proceedings of the 2nd International Conference on Learning Representations (ICLR)_ , Banff, Canada, 2014.
  * Kingma and Dhariwal [2018] D. P. Kingma and P. Dhariwal, “Glow: Generative flow with invertible 1x1 convolutions,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , Montréal, Canada, 2018, pp. 10 236–10 245.
  * Radford et al. [2016] A. Radford, L. Metz, and S. Chintala, “Unsupervised representation learning with deep convolutional generative adversarial networks,” in _Proceedings of the 4th International Conference on Learning Representations (ICLR)_ , San Juan, Puerto Rico, 2016.
  * Gao et al. [2020b] R. Gao, E. Nijkamp, D. P. Kingma, Z. Xu, A. M. Dai, and Y. N. Wu, “Flow contrastive estimation of energy-based models,” in _Proceedings of the 2020 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , Seattle, WA, 2020, pp. 7515–7525.
  * Arbel et al. [2021] M. Arbel, L. Zhou, and A. Gretton, “Generalized energy based models,” in _Proceedings of the 9th International Conference on Learning Representations (ICLR)_ , Virtual Event, 2021.
  * Zagoruyko and Komodakis [2016a] S. Zagoruyko and N. Komodakis, “Wide residual networks,” _ArXiv_ , vol. abs/1605.07146, 2016.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 770–778, 2016.
  * Karras et al. [2022b] T. Karras, M. Aittala, T. Aila, and S. Laine, “Elucidating the design space of diffusion-based generative models,” _ArXiv_ , vol. abs/2206.00364, 2022.
  * Heusel et al. [2017] 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,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 6626–6637.
  * Xie et al. [2022b] J. Xie, Y. Zhu, J. L. Li, and P. Li, “A tale of two flows: Cooperative learning of langevin flow and normalizing flow toward energy-based model,” _ArXiv_ , vol. abs/2205.06924, 2022.
  * Stein [1981] C. M. Stein, “Estimation of the mean of a multivariate normal distribution,” _The annals of Statistics_ , pp. 1135–1151, 1981.
  * Zagoruyko and Komodakis [2016b] S. Zagoruyko and N. Komodakis, “Wide residual networks,” _arXiv preprint arXiv:1605.07146_ , 2016.
  * Elfwing et al. [2017] S. Elfwing, E. Uchibe, and K. Doya, “Sigmoid-weighted linear units for neural network function approximation in reinforcement learning,” _Neural networks : the official journal of the International Neural Network Society_ , vol. 107, pp. 3–11, 2017.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _CoRR_ , vol. abs/1412.6980, 2014.

##  Appendix A Technical details

###  A.1 Proof of Theorem 2 (Section 3.1)

Before we prove Theorem 2, we give two lemmas to simplify the proof.

######  Lemma 4.

Assume function $p$ is positive and twice differentiable, then the following
identity holds

$\Delta p(\bm{x})=p(\bm{x})\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}+p(\bm{x})\Delta\log p(\bm{x}),$

where

$\displaystyle\nabla_{\bm{x}}\log p(\bm{x})=\sum_{i=1}^{D}\frac{\partial
p(\bm{x})}{\partial{x}_{i}},\ \Delta\log
p(\bm{x})=\sum_{i=1}^{D}\frac{\partial^{2}\log
p(\bm{x})}{\partial{x}_{i}^{2}}.$

Here $x_{i}$ represents the $i$-th covariate of vector $\bm{x}$ and $D$ is the
data dimension.

###### Proof.

With a slight abuse of notation, we write
$\langle\nabla_{\bm{x}},\bm{f}(\bm{x})\rangle\coloneqq\sum_{i=1}^{D}\partial\bm{f}_{i}(\bm{x})/\partial{x}_{i}$.
Since $p$ is twice differentiable, we have

$\displaystyle\Delta p(\bm{x})$
$\displaystyle=\langle\nabla_{\bm{x}},\nabla_{\bm{x}}p(\bm{x})\rangle$
$\displaystyle=\langle\nabla_{\bm{x}},p(\bm{x})\nabla_{\bm{x}}\log
p(\bm{x})\rangle=\langle\nabla_{\bm{x}}p(\bm{x}),\nabla_{\bm{x}}\log
p(\bm{x})\rangle+p(\bm{x})\langle\nabla_{\bm{x}},\nabla_{\bm{x}}\log
p(\bm{x})\rangle$ $\displaystyle=p(\bm{x})\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}+p(\bm{x})\Delta\log p(\bm{x})$

∎

######  Lemma 5.

Assume function $p(.)$ is a positive and twice differentiable probability
density. Then the following identity holds

$\displaystyle\mathbb{E}_{p(\bm{x})}\log p(\bm{x})\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}=-\mathbb{E}_{p(\bm{x})}\big{[}\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}+\log p(\bm{x})\Delta\log p(\bm{x})\big{]}.$

###### Proof.

Notice that

$\displaystyle\mathbb{E}_{p(\bm{x})}\log p(\bm{x})\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}$ $\displaystyle=\mathbb{E}_{p(\bm{x})}\langle\log
p(\bm{x})\nabla_{\bm{x}}\log p(\bm{x}),\nabla_{\bm{x}}\log p(\bm{x})\rangle$
$\displaystyle=\mathbb{E}_{p(\bm{x})}-\langle\nabla_{\bm{x}},\log
p(\bm{x})\nabla_{\bm{x}}\log p(\bm{x})\rangle.$

The above equality holds because of Stein’s identity [55], i.e.,

$\mathbb{E}_{p(\bm{x})}\langle\bm{f}(\bm{x}),\nabla_{\bm{x}}\log
p(\bm{x})\rangle=-\mathbb{E}_{p(\bm{x})}\langle\nabla_{\bm{x}},\bm{f}(\bm{x})\rangle$

for vector value function $\bm{f}$ which lies in Stein class of $p$ 111A
vector-value function $\bm{f}$ lies in Stein class of distribution $p$ means
three conditions hold:  • $\bm{f}$ is 2nd-order smooth; • both
$\|\bm{f}\|_{2}^{2}$ and $\|\nabla_{\bm{x}}\bm{f}^{T}\|_{F}^{2}$ is integrable
w.r.t. $p$. The notation $\|.\|_{F}$ represents the Frobenius norm. •
$p(\bm{x})\|\nabla_{\bm{x}}\bm{f}^{T}(\bm{x})\|_{F}\to 0$ when
$\|\bm{x}\|_{2}\to\partial support(p)$ . Thus the proof is finished with

$\displaystyle\mathbb{E}_{p(\bm{x})}\log p(\bm{x})\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}=-\mathbb{E}_{p(\bm{x})}\bigg{[}\|\nabla_{\bm{x}}\log
p(\bm{x})\|_{2}^{2}+\log p(\bm{x})\Delta\log p(\bm{x})\bigg{]}$

∎

We give the proof for Theorem 2 with the above two lemmas 4 and 5.

###### Proof.

Recall that the two distributions $p,q$ evolve along a general Ito’s diffusion
process

$d\bm{x}_{t}=\bm{F}(\bm{x}_{t},t)\mathrm{d}t+\bm{G}(t)\mathrm{d}\bm{w}_{t}.$

Here $\bm{F}(\bm{x},t)$ is a vector value function, and $\bm{G}(t)$ is a
scalar function of $t$. Note that $p_{0}=p,q_{0}=q$. We denote
$p^{(F,G,t)}_{t},q^{(F,G,t)}_{t}$ as $p_{t},q_{t}$ for short. The KL
divergence between $p_{t},q_{t}$ is defined as

$\mathcal{D}_{KL}(p_{t},q_{t})=\mathbb{E}_{p_{t}}\log\frac{p_{t}(\bm{x})}{q_{t}(\bm{x})}=\int
p_{t}\log\frac{p_{t}}{q_{t}}\mathrm{d}\bm{x}.$

We declare all integrals are w.r.t. $\bm{x}$ and omit the $\mathrm{d}\bm{x}$
in integral formulas for simplification. The change rate of KL divergence is

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{D}_{KL}(p_{t},q_{t})\frac{\mathrm{d}}{\mathrm{d}t}\int
p_{t}(\bm{x})\log\frac{p_{t}(\bm{x})}{q_{t}(\bm{x})}$
$\displaystyle=\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}\log
p_{t}-\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}\log
q_{t}+\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}-\int\frac{p_{t}}{q_{t}}\frac{\mathrm{d}q_{t}}{\mathrm{d}t}$
(13) $\displaystyle:=A+B+C+D.$

The third term

$C=\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\int
p_{t}=\frac{\mathrm{d}}{\mathrm{d}t}1=0.$

Hence the above equation remains 3 terms. By the Fokker-Planck equation
equation 8, the evolved density follows

$\displaystyle\frac{\mathrm{d}p_{t}}{\mathrm{d}t}$
$\displaystyle=-\langle\nabla_{\bm{x}},p_{t}\bm{F}\rangle+\frac{1}{2}\bm{G}^{2}(t)\Delta
p_{t}$ $\displaystyle=-p_{t}\langle\nabla_{\bm{x}}\log p_{t},\bm{F}\rangle-
p_{t}\langle\nabla_{\bm{x}},\bm{F}\rangle+\frac{1}{2}\bm{G}^{2}(t)p_{t}\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)p_{t}\Delta\log p_{t}.$

Substitute the above equation to equation (13), the first term $A$ becomes

$\displaystyle\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}\log p_{t}$ (14)
$\displaystyle=\int p_{t}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)\Delta\log
p_{t}-\bm{F}(x,t)^{T}\nabla_{\bm{x}}\log
p_{t}-\langle\nabla_{\bm{x}},\bm{F}(\bm{x},t)\rangle\bigg{]}\log p_{t}$
$\displaystyle=\mathbb{E}_{p_{t}}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\log
p_{t}\|\nabla_{\bm{x}}\log p_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)\log
p_{t}\Delta\log p_{t}$ $\displaystyle-\log
p_{t}\langle\bm{F}(\bm{x},t),\nabla\log p_{t}\rangle-\log
p_{t}\langle\nabla_{\bm{x}},\bm{F}(\bm{x},t)\rangle\bigg{]}$
$\displaystyle=\mathbb{E}_{p_{t}}\bigg{[}-\frac{1}{2}\bm{G}^{2}(t)\big{[}\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\log p_{t}\Delta\log p_{t}\big{]}+\frac{1}{2}\bm{G}^{2}(t)\log
p_{t}\Delta\log p_{t}$ $\displaystyle-\log p_{t}\langle\bm{F}(x,t),\nabla\log
p_{t}\rangle-\log
p_{t}\langle\nabla_{\bm{x}},\bm{F}(\bm{x},t)\rangle\bigg{]}.$

By Stein’s identity,

$\displaystyle\mathbb{E}_{p_{t}}\log p_{t}\langle\bm{G},\nabla_{\bm{x}}\log
p_{t}\rangle=\mathbb{E}_{p_{t}}\langle(\log p_{t})\bm{F},\nabla_{\bm{x}}\log
p_{t}\rangle$ $\displaystyle=-\mathbb{E}_{p_{t}}\langle\nabla_{\bm{x}},\log
p_{t}\bm{F}\rangle=-\mathbb{E}_{p_{t}}\bigg{[}\langle\nabla_{\bm{x}}\log
p_{t},\bm{F}\rangle+\log p_{t}\langle\nabla_{\bm{x}},\bm{F}\rangle\bigg{]}.$

This term (14) becomes

$\displaystyle\mathbb{E}_{p_{t}}\bigg{[}-\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}-\log p_{t}\langle\bm{F},\nabla_{\bm{x}}\log p_{t}\rangle-\log
p_{t}\langle\nabla_{\bm{x}},\bm{F}\rangle\bigg{]}$
$\displaystyle=\mathbb{E}_{p_{t}}\bigg{[}-\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\big{[}\langle\bm{F},\nabla_{\bm{x}}\log p_{t}\rangle+\log
p_{t}\langle\nabla_{\bm{x}},\bm{F}\rangle\big{]}-\log
p_{t}\langle\nabla_{\bm{x}},\bm{F}\rangle\bigg{]}$
$\displaystyle=\mathbb{E}_{p_{t}}\bigg{[}-\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\langle\bm{F},\nabla_{\bm{x}}\log p_{t}\rangle\bigg{]}.$

Next, we calculate the second term $B$ with a similar argument

$\displaystyle B=\int-\frac{\mathrm{d}p_{t}}{\mathrm{d}t}\log q_{t}$
$\displaystyle=-\mathbb{E}_{p_{t}}\log
q_{t}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)\Delta\log
p_{t}-\langle\bm{F},\nabla_{\bm{x}}\log
p_{t}\rangle-\langle\nabla_{\bm{x}},\bm{F}\rangle\bigg{]}$
$\displaystyle=-\mathbb{E}_{p_{t}}\bigg{[}-\frac{1}{2}\bm{G}^{2}(t)\langle\nabla_{\bm{x}}\log
p_{t},\nabla_{\bm{x}}\log q_{t}\rangle+\langle\bm{F},\nabla_{\bm{x}}\log
q_{t}\rangle\bigg{]}.$

The fourth term $D$ writes

$\displaystyle D$
$\displaystyle=\int-\frac{p_{t}}{q_{t}}\frac{\mathrm{d}q_{t}}{\mathrm{d}t}$
$\displaystyle=-\int\frac{p_{t}}{q_{t}}q_{t}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
q_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)\Delta\log
q_{t}-\langle\bm{F},\nabla_{\bm{x}}\log
q_{t}\rangle-\langle\nabla_{\bm{x}},\bm{F}\rangle\bigg{]}$
$\displaystyle=-\mathbb{E}_{p_{t}}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
q_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)\Delta\log
q_{t}-\langle\bm{F},\nabla_{\bm{x}}\log
q_{t}\rangle-\langle\nabla_{\bm{x}},\bm{F}\rangle\bigg{]}.$

With Stein’s identity,

$\mathbb{E}_{p_{t}}\langle\nabla_{\bm{x}}\log q_{t},\nabla_{\bm{x}}\log
p_{t}\rangle=-\mathbb{E}_{p_{t}}\langle\nabla_{\bm{x}},\nabla_{\bm{x}}\log
q_{t}\rangle=-\mathbb{E}_{p_{t}}\Delta\log q_{t}.$

Substitute the above equality to the fourth term, we have

$\displaystyle
D=-\mathbb{E}_{p_{t}}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
q_{t}\|_{2}^{2}-\frac{1}{2}\bm{G}^{2}(t)\langle\nabla_{\bm{x}}\log
q_{t},\nabla_{\bm{x}}\log p_{t}\rangle-\langle\bm{F},\nabla_{\bm{x}}\log
q_{t}\rangle+\langle\bm{F},\nabla_{\bm{x}}\log p_{t}\rangle\bigg{]}.$

Combine all three terms, we have

$\displaystyle\frac{d}{\mathrm{d}t}\mathcal{D}_{KL}(p_{t},q_{t})=\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}\log
p_{t}-\int\frac{\mathrm{d}p_{t}}{\mathrm{d}t}\log
q_{t}-\int\frac{p_{t}}{q_{t}}\frac{\mathrm{d}q_{t}}{\mathrm{d}t}$
$\displaystyle=-\mathbb{E}_{p_{t}}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}\|_{2}^{2}+\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
q_{t}\|_{2}^{2}-\bm{G}^{2}(t)\langle\nabla_{\bm{x}}\log
p_{t},\nabla_{\bm{x}}\log q_{t}\rangle\bigg{]}$
$\displaystyle=-\frac{1}{2}\mathbb{E}_{p_{t}}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}(x)-\nabla_{\bm{x}}\log q_{t}(x)\|_{2}^{2}$ (15)

So the integral representation writes

$\displaystyle\mathcal{D}_{KL}(p_{T},q_{T})-\mathcal{D}_{KL}(p_{0},q_{0})=\int_{0}^{T}\frac{d}{\mathrm{d}t}\mathcal{D}_{KL}(p_{t},q_{t})\mathrm{d}t$
$\displaystyle=-\int_{0}^{T}\frac{1}{2}\mathbb{E}_{p_{t}}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{t}(\bm{x})-\nabla_{\bm{x}}\log q_{t}(\bm{x})\|_{2}^{2}\mathrm{d}t.$

∎

###  A.2 Proof of Langevin dynamic’s stationary property

The stationary property states that $p_{\theta}$ is stationary under EBM-
induced Langevin dynamics.

###### Proof.

Notice that the evolution of a probability under EBM Langevin dynamics 3 is
governed by the Fokker-Planck equation equation 8

$\frac{\mathrm{d}}{\mathrm{d}t}p(\bm{x},t)=-\langle\nabla_{\bm{x}},\frac{1}{2}\nabla_{\bm{x}}\log
p_{\theta}(\bm{x})p(\bm{x},t)\rangle+\frac{1}{2}\Delta_{\bm{x}}p(\bm{x},t).$

Since $\Delta
p(\bm{x},t)=\langle\nabla_{\bm{x}},\nabla_{\bm{x}}p(\bm{x},t)\rangle$, we have

$\displaystyle\Delta_{\bm{x}}p(\bm{x},t)=\langle\nabla_{\bm{x}},\nabla_{\bm{x}}p(\bm{x},t)\rangle=\langle\nabla_{\bm{x}},p(\bm{x},t)\nabla_{\bm{x}}\log
p(\bm{x},t)\rangle$

Combining the above, we have the simplified Fokker-Planck equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p(\bm{x},t)=\frac{1}{2}\langle\nabla_{\bm{x}},\frac{1}{2}\big{[}\nabla_{\bm{x}}\log
p(\bm{x},t)-\nabla_{\bm{x}}\log p_{\theta}(\bm{x})\big{]}p(\bm{x},t)\rangle$

Substitute $p(\bm{x},t)=p_{\theta}(\bm{x})$, we have

$\frac{\mathrm{d}}{\mathrm{d}t}p(\bm{x},t)=0.$

So $p(\bm{x},t)=p_{\theta}(\bm{x})$ is stationary under $p_{\theta}$ induced
Langevin dynamics. ∎

###  A.3 Non-negativity of CD (Theorem 2)

Recall the definition of CD equation 4,

$\displaystyle\mathcal{D}_{CD}(p_{d},p_{\theta})=\mathcal{D}_{KL}(p_{d},p_{\theta})-\mathcal{D}_{KL}(p_{d,\theta}^{(T)},p_{\theta}),$

The non-negativity of CD in fact comes as a corollary of Theorem 2 as we have
proved in A.1.

###### Proof.

Recall the definition of CD,

$\mathcal{D}_{CD}(p_{d},p_{\theta})=\mathcal{D}_{KL}(p_{d},p_{\theta})-\mathcal{D}_{KL}(p_{d}^{(T)}(\theta),p_{\theta}).$

Here the $p_{d}^{(T)}(\theta)$ denote the $T$ time evolved EBM distribution
under EBM Langevin dynamcis

$d\bm{x}_{t}=\frac{1}{2}\nabla_{\bm{x}_{t}}\log
p_{\theta}(\bm{x}_{t})\mathrm{d}t+\mathrm{d}\bm{w}_{t}.$

Recall that $p_{\theta}(\bm{x})=p_{\theta}^{(T)}(\bm{x})$ as we show in A.2,
then in definition of CD and CD equals to

$\mathcal{D}_{CD}(p_{d},p_{\theta})=\mathcal{D}_{KL}(p_{d},p_{\theta})-\mathcal{D}_{KL}(p_{d}^{(T)}(\theta),p_{\theta}^{(T)}).$

By Theorem 2,

$\displaystyle\mathcal{D}_{CD}(p_{d},p_{\theta})=\frac{1}{2}\int_{0}^{T}\mathbb{E}_{\bm{x}_{t}\sim
p_{d}^{(t)}(\bm{x}_{t})}\|\nabla_{\bm{x}_{t}}\log
p_{\theta}^{(t)}(\bm{x}_{t})-\nabla_{\bm{x}_{t}}\log
q^{(t)}(\bm{x}_{t})\|^{2}_{2}\mathrm{d}t\geq 0$

∎

###  A.4 Non-negaligibility of the extra term of CD (Equation 5)

Recall the gradient formula of CD equation 5.

$\displaystyle\frac{\partial}{\partial\theta}\mathcal{D}_{CD}(p_{d},p_{\theta})=\mathbb{E}_{p_{d,\theta}^{(T)}}\big{[}\frac{\partial}{\partial\theta}f_{\theta}(\bm{x})\big{]}-\mathbb{E}_{p_{d}}\big{[}\frac{\partial}{\partial\theta}f_{\theta}(\bm{x})\big{]}-\mathbb{E}_{p_{d,\theta}^{(T)}}\bigg{[}\log
p_{\theta}(\bm{x})\frac{\partial}{\partial\theta}\log
p_{d,\theta}^{(T)}(\bm{x})\bigg{]}.$

The third term is

$\displaystyle(3)=-\mathbb{E}_{p_{d,\theta}^{(T)}}\bigg{[}\log
p_{\theta}(\bm{x})\frac{\partial}{\partial\theta}\log
p_{d,\theta}^{(T)}(\bm{x})\bigg{]}.$ (16)

For ease of expression, we may omit the notation $\mathrm{d}\bm{x}$ in the
integral. If $p_{d,\theta}^{T}(\bm{x})\to p_{\theta}(\bm{x})$ as we assumed,
the term 16 turns to

$\displaystyle(3)$ $\displaystyle=-\mathbb{E}_{p_{\theta}}\bigg{[}\log
p_{\theta}(\bm{x})\frac{\partial}{\partial\theta}\log
p_{\theta}(\bm{x})\bigg{]}$ $\displaystyle=-\int p_{\theta}(\bm{x})\log
p_{\theta}(\bm{x})\frac{1}{p_{\theta}(\bm{x})}\frac{\partial}{\partial\theta}p_{\theta}(\bm{x})\mathrm{d}\bm{x}$
$\displaystyle=-\int\log
p_{\theta}(\bm{x})\frac{\partial}{\partial\theta}p_{\theta}(\bm{x})\mathrm{d}\bm{x}$
$\displaystyle=-\frac{\partial}{\partial\theta}\int\log
p_{\theta}(\bm{x})p_{\theta}(\bm{x})\mathrm{d}\bm{x}+\int
p_{\theta}(x)\frac{\partial}{\partial\theta}\log p_{\theta}(\bm{x})$
$\displaystyle=-\frac{\partial}{\partial\theta}\mathbb{E}_{p_{\theta}}\log
p_{\theta}(\bm{x})+\int\frac{\partial}{\partial\theta}p_{\theta}(\bm{x})$
$\displaystyle=-\frac{\partial}{\partial\theta}\mathbb{E}_{p_{\theta}}\log
p_{\theta}(\bm{x})+\frac{\partial}{\partial\theta}\int p_{\theta}(\bm{x})$
(17)
$\displaystyle=-\frac{\partial}{\partial\theta}\mathbb{E}_{p_{\theta}}\log
p_{\theta}(\bm{x})+\frac{\partial}{\partial\theta}\mathbf{1}$
$\displaystyle=-\frac{\partial}{\partial\theta}\mathbb{E}_{p_{\theta}}\log
p_{\theta}(\bm{x})$

The equality 17 holds if $p_{\theta}(x)$ satisfies the conditions. Now if the
density function is satisfied the condition that (1). $p_{\theta}(x)$ is
Lebesgue integrable for $\bm{x}$ with each $\theta$; (2). For almost all
$\bm{x}\in\mathbf{R}^{D}$, the partial derivative $\partial
p_{\theta}(\bm{x})/\partial\theta$ exists for all $\theta\in\Theta$. (3) there
exists an integrable function $g(.):\mathbf{R}^{D}\to\mathbf{R}$, such that
$p_{\theta}(\bm{x})\leq g(\bm{x})$ for all $\bm{x}$ in its domain. Then the
derivative w.r.t $\theta$ can be exchanged with the integral over $\bm{x}$,
i.e.

$\displaystyle\int\frac{\partial}{\partial\theta}p_{\theta}(\bm{x})\mathrm{d}\bm{x}=\frac{\partial}{\partial\theta}\int
p_{\theta}(\bm{x})\mathrm{d}\bm{x}.$

###  A.5 Tractable form of DCD with parameter-free diffusion

Thus the DCD under the above SDE has the form

$\displaystyle\mathcal{D}_{DCD}^{(\bm{F},\bm{G},T)}(p_{d},p_{\theta})$
$\displaystyle=\mathbb{E}_{\bm{x}_{0}\sim p_{d}}\big{[}\log
p_{d}(\bm{x}_{0})-f_{\theta}(\bm{x}_{0})\big{]}-\mathbb{E}_{\bm{x}_{0}\sim
p_{d},\atop\bm{x}_{t}\sim p(\bm{x}_{t}|\bm{x}_{0})}\big{[}\log
p_{d}^{(\bm{F},\bm{G},T)}(\bm{x}_{t})-f_{\theta}^{(\bm{F},\bm{G},T)}(\bm{x}_{t})\big{]}.$

Here $f_{\theta}^{(\bm{F},\bm{G},T)}$ are time $T$ marginal energy under
diffusion process (equation 7). The term $\log p_{d}(\bm{x})$ and $\log
p_{d}^{(\bm{F},\bm{G},T)}$ are independent of parameter $\theta$ since the
diffusion process is parameter-free. As a result, we can drop them when using
gradient-based optimization algorithms. Thus, we have the final tractable
learning objective based on DCD as

$\displaystyle\mathcal{L}_{DCD}(\theta)=$
$\displaystyle\mathbb{E}_{\bm{x}_{0}\sim p_{d},\bm{x}_{t}\sim
p(\bm{x}_{t}|\bm{x}_{0})}\big{[}f_{\theta}^{(\bm{F},\bm{G},T)}(\bm{x}_{t})\big{]}-\mathbb{E}_{\bm{x}_{0}\sim
p_{d}}\big{[}f_{\theta}(\bm{x}_{0})\big{]}.$

###  A.6 More backgrounds on VE diffusion

#### The VE Diffusion.

Recall the VE diffusion equation 9,

$\displaystyle\mathrm{d}\bm{x}_{t}=g(t)\mathrm{d}\bm{w}_{t},$

is also favored for its easy-to-simulate property. The marginal transition of
VE writes

$\displaystyle
p(\bm{x}_{t}|\bm{x}_{0})=\mathcal{N}(\bm{x}_{0},\sigma(t)\mathbf{I}).$ (18)

$\sigma(t)=\int_{0}^{t}g(s)\mathrm{d}s$. Similar to the VP diffusion, marginal
samples of the VE diffusion are also cheap to obtain and parameter-free.

###  A.7 Proof of Proposition 1 (Section 3.1)

Recall the Proposition 3.3.

######  Proposition.

Assume $p_{\theta}^{(0)}(\bm{x})=e^{f_{\theta}(\bm{x})}/Z_{\theta}$ where
$Z_{\theta}$ is a parameter-dependent normalizing constant. Assume
$p_{\theta}^{(t)}$ denotes the evolved density along a diffusion process
equation 7, then for any fixed $\bm{x}$, the energy value
$p_{\theta}^{(t)}(\bm{x})$ evolves according to a PDE

$d\log p_{\theta}^{(t)}(\bm{x})/\mathrm{d}t=\mathcal{O}(\nabla_{\bm{x}}\log
p_{\theta}^{(t)}),$

where $\mathcal{O}(\nabla_{\bm{x}}\log p_{\theta}^{(t)})$ is the following
operator which is independent of the normalizing constant,

$\displaystyle\langle\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t),\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\rangle+\langle\nabla,\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t)\rangle.$

###### Proof.

Following the Fokker-Planck equation 8, the density
$p_{\theta}^{(t)}(\bm{x})=e^{f_{\theta}^{(t)}(\bm{x})}/Z_{\theta}$ evolves
with equation:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}p_{\theta}^{(t)}(\bm{x})=$
$\displaystyle=\langle\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t),\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\rangle+\langle\nabla,\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t)\rangle.$

Denote $f_{\theta}^{(t)}(\bm{x})=\log p_{\theta}^{(t)}(\bm{x})+\log
Z^{(t)}_{\theta}$. Then we have $\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})=\nabla_{\bm{x}}f_{\theta}^{(t)}(\bm{x})$. The
evolution of $f_{\theta}^{(t)}$ thus follows a partial differential equation
(Fokker-Planck equation), i.e.,

$\displaystyle\mathrm{d}f_{\theta}^{(t)}(\bm{x})/\mathrm{d}t=\langle\bm{G}^{2}(t)\nabla_{\bm{x}}f_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t),\nabla_{\bm{x}}f_{\theta}^{(t)}(\bm{x})\rangle+\langle\nabla_{\bm{x}},\bm{G}^{2}(t)\nabla_{\bm{x}}f_{\theta}^{(t)}(\bm{x})/2-\bm{F}(x,t)\rangle$
$\displaystyle=\mathcal{O}(f_{\theta}^{(t)},\bm{x}).$

In the above equation,

$\mathcal{O}(f)=\bm{G}^{2}\|\nabla_{\bm{x}}f\|^{2}-\langle\bm{F}(.),\nabla_{\bm{x}}\bm{F}\rangle+\bm{G}^{2}/2\Delta
f-\langle\nabla_{\bm{x}},\bm{F}(.)\rangle.$

Thus the $T$ time energy function equals

$f_{\theta}^{(T)}(\bm{x})=f_{\theta}(\bm{x})+\int_{t=0}^{T}\mathcal{O}(f_{\theta}^{(t)},\bm{x})\mathrm{d}t,$

where $f_{\theta}^{(t)}$ is the solution of the above energy diffusion ODE. We
can derive the change of normalizing constant with the following argument. By
writing

$f_{\theta}^{(t+\mathrm{d}t)}=\mathcal{O}(f_{\theta}^{(t)})\mathrm{d}t+f_{\theta}^{(t)},$

we have

$\exp(f_{\theta}^{(t+\mathrm{d}t)})=\exp(f_{\theta}^{(t)})\exp(\mathcal{O}(f_{\theta}^{(t)})\mathrm{d}t)=\exp(f_{\theta}^{(t)})(1+\mathcal{O}(f_{\theta}^{(t)})+o(\mathrm{d}t^{2})).$

Taking integral w.r.t $x$ on both sides, we have

$\displaystyle Z_{\theta}^{(t+\mathrm{d}t)}$
$\displaystyle=\int\exp(f_{\theta}^{(t+\mathrm{d}t)})\mathrm{d}\bm{x}$
$\displaystyle=\int\exp(f_{\theta}^{(t)})\bigg{[}1+\mathcal{O}(f_{\theta}^{(t)})+o(\mathrm{d}t^{2})\bigg{]}dx$
$\displaystyle=Z_{\theta}^{(t)}\bigg{[}1+\int\frac{\exp(f_{\theta}^{(t)})(\bm{x})}{Z_{\theta}^{(t)}}\mathcal{O}(f_{\theta}^{(t)})(\bm{x})\mathrm{d}\bm{x}\bigg{]}+o(\mathrm{d}t^{2})$
$\displaystyle=Z_{\theta}^{(t)}\bigg{[}1+\int\mathbb{E}_{p_{\theta}^{(t)}}\mathcal{O}(f_{\theta}^{(t)})(\bm{x})\mathrm{d}\bm{x}\bigg{]}+o(\mathrm{d}t^{2})$
$\displaystyle=Z_{\theta}^{(t)}\bigg{[}1+\mathbb{E}_{p_{\theta}^{(t)}}\mathcal{O}(f_{\theta}^{(t)})(\bm{x})\bigg{]}+o(\mathrm{d}t^{2}).$
(19)

Note that

$\displaystyle\mathbb{E}_{p_{\theta}^{(t)}}\mathcal{O}(f_{\theta}^{(t)})(\bm{x})=\mathbb{E}_{p_{\theta}^{(t)}}\mathcal{O}(\log
p_{\theta}^{(t)})(\bm{x})$
$\displaystyle=\mathbb{E}_{p_{\theta}^{(t)}}\bigg{[}\langle\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t),\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\rangle+\langle\nabla_{\bm{x}},\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})/2-\bm{F}(\bm{x},t)\rangle\bigg{]}$
$\displaystyle=\mathbb{E}_{p_{\theta}^{(t)}}\bigg{[}\frac{1}{2}\bm{G}^{2}(t)\|\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\|^{2}_{2}+\bm{G}^{2}(t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})+\bm{F}^{T}(\bm{x},t)\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})+\nabla_{\bm{x}}\bm{F}(\bm{x},t)\bigg{]}.$

Through Stein’s identity, we have

$\mathbb{E}_{p_{\theta}^{(t)}}\bigg{[}\|\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\|_{2}^{2}+\Delta_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})\bigg{]}=0,$
$\mathbb{E}_{p_{\theta}^{(t)}}\bigg{[}\bm{F}(\bm{x},t)^{T}\nabla_{\bm{x}}\log
p_{\theta}^{(t)}(\bm{x})+\nabla_{\bm{x}}\mathbf{F}(\bm{x},t)\bigg{]}=0.$

Thus we have

$\displaystyle\mathbb{E}_{p_{\theta}^{(t)}}\mathcal{O}(f_{\theta}^{(t)})(\bm{x})=0$
(20)

Substituting equation (20) into (19), we have

$\displaystyle
Z_{\theta}^{(t+\mathrm{d}t)}=Z_{\theta}^{(t)}+o(\mathrm{d}t^{2}).$

Thus

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}Z_{\theta}^{(t)}=0.$

The normalizing constant remains unchanged. Since the operator $\mathcal{O}$
only depends on the $\nabla_{\bm{x}}f$ term, the normalizing constant does not
influence the energy evolution. Then we have $\log
p_{\theta}^{(t)}(\bm{x})=f_{\theta}^{(t)}(\bm{x})-\log Z_{\theta}$. So the
normalizing constant $\log z_{\theta}$ can be abstracted in the loss function
as

$\displaystyle\mathcal{L}_{DCD}(\theta)=\mathbb{E}_{\bm{x}_{T}\sim
p_{T}(\bm{x}_{T})}\big{[}\log
p_{\theta}^{(T)}(\bm{x}_{T})\big{]}-\mathbb{E}_{\bm{x}_{0}\sim
p_{d}}\big{[}\log p_{\theta}^{(0)}(\bm{x}_{0})\big{]}$
$\displaystyle=\mathbb{E}_{\bm{x}_{T}\sim
p_{T}(\bm{x}_{T})}\big{[}f_{\theta}^{(T)}(\bm{x}_{T})-\log
Z_{\theta}\big{]}-\mathbb{E}_{\bm{x}_{0}\sim p_{d}}\big{[}\log
p_{\theta}^{(0)}(\bm{x}_{0})-\log Z_{\theta}\big{]}$
$\displaystyle=\mathbb{E}_{\bm{x}_{T}\sim
p_{T}(\bm{x}_{T})}\big{[}f_{\theta}^{(T)}(\bm{x}_{T})\big{]}-\mathbb{E}_{\bm{x}_{0}\sim
p_{d}}\big{[}\log p_{\theta}^{(0)}(\bm{x}_{0})\big{]}.$

∎

###  A.8 Detailed proof of connections to DRL (Section 3.2)

The expected recovery likelihood is

$\mathbb{E}_{\bm{x}\sim p_{d},\tilde{\bm{x}}\sim
p_{d}^{(\sigma)}(\tilde{\bm{x}})}\log p_{\theta}(\bm{x}|\tilde{\bm{x}}).$

Since $p^{(\sigma)}(\tilde{\bm{x}}|\bm{x})$ and $p_{d}(\bm{x})$ are
independent of parameter $\theta$, the objective is equivalent to minimizing

$\displaystyle-\mathbb{E}_{p_{d}(\bm{x})p^{(\sigma)}(\tilde{\bm{x}}|\bm{x})}\bigg{[}\log\frac{p_{\theta}(\bm{x}|\tilde{\bm{x}})p_{\theta}^{(\sigma)}(\tilde{\bm{x}})}{p_{d}(\bm{x}|\tilde{\bm{x}})p_{d}^{(\sigma)}(\tilde{\bm{x}})}-\log\frac{p_{\theta}^{(\sigma)}(\tilde{\bm{x}})}{p_{d}^{(\sigma)}(\tilde{\bm{x}})}\bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{p_{d}(\bm{x})p_{\sigma}(\tilde{\bm{x}}|\bm{x})}\bigg{[}\log\frac{p_{d}(\bm{x},\tilde{\bm{x}})}{p_{\theta}(\bm{x},\tilde{\bm{x}})}\bigg{]}-\mathcal{D}_{KL}(p_{d}^{(\sigma)}(\tilde{\bm{x}}),p_{\theta}^{(\sigma)}(\tilde{\bm{x}}))$
$\displaystyle=$
$\displaystyle\mathbb{E}_{p_{d}(\bm{x})p_{\sigma}(\tilde{\bm{x}}|\bm{x})}\bigg{[}\log\frac{p_{d}(\bm{x})p(\tilde{\bm{x}}|\bm{x})}{p_{\theta}(\bm{x})p(\tilde{\bm{x}}|\bm{x})}\bigg{]}-\mathcal{D}_{KL}(p_{d}^{(\sigma)}(\tilde{\bm{x}}),p_{\theta}^{(\sigma)}(\tilde{\bm{x}}))$
$\displaystyle=$
$\displaystyle\mathbb{E}_{p_{d}(\bm{x})p_{\sigma}(\tilde{\bm{x}}|\bm{x})}\bigg{[}\log\frac{p_{d}(\bm{x})}{p_{\theta}(\bm{x})}\bigg{]}-\mathcal{D}_{KL}(p_{d}^{(\sigma)}(\tilde{\bm{x}}),p_{\theta}^{(\sigma)}(\tilde{\bm{x}}))$
$\displaystyle=$
$\displaystyle\mathcal{D}_{KL}(p_{d}(\bm{x}),p_{\theta}(\bm{x}))-\mathcal{D}_{KL}(p_{d}^{(\sigma)}(\tilde{\bm{x}}),p_{\theta}^{(\sigma)}(\tilde{\bm{x}})).$

###  A.9 Proof of Proposition 2 (Section 3.3)

When $t$ is small and by the first-order Taylor approximation, we can write

$f_{\theta}^{(t)}(\bm{x})=f_{\theta}(\bm{x})+t\big{[}\frac{d}{\mathrm{d}t}f_{\theta}^{(t)}(\bm{x})\big{]}|_{t=0}+{o}(t),$

The corresponding DCD objective becomes

$\displaystyle\mathcal{L}_{DCD}^{(VE)}(\theta)=\mathbb{E}_{\bm{x}_{0}\sim
p_{d},\bm{x}_{t}\sim
p(\bm{x}_{t}|\bm{x}_{0})}[f_{\theta}^{(t)}(\bm{x}_{t})]-\mathbb{E}_{\bm{x}_{0}\sim
p_{0}}[f_{\theta}(\bm{x}_{0})]$
$\displaystyle=\mathbb{E}_{p_{t}(\bm{x}_{t})}[f_{\theta}^{(t)}(\bm{x}_{t})-f_{\theta}(\bm{x}_{t})]+\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})]$
$\displaystyle=\mathbb{E}_{p_{t}}t[\frac{d}{\mathrm{d}t}f_{\theta}^{(t)}(x_{t})]|_{t=0}+\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})]$
$\displaystyle=\mathbb{E}_{p_{t}}\frac{1}{2}\bm{G}^{2}(0)\big{[}\|\nabla_{\bm{x}}f_{\theta}(\bm{x})\|^{2}+\Delta_{\bm{x}}f_{\theta}(\bm{x})\big{]}+\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})].$

###  A.10 Derivation of DCD-VE (Equation 12)

$d\bm{x}_{t}=\bm{G}(t)d\bm{w}_{t}$, the energy evolution is

$df_{\theta}^{(t)}(\bm{x})/\mathrm{d}t=\frac{1}{2}\bm{G}^{2}(t)\big{[}\|\nabla_{\bm{x}}f_{\theta}^{(t)}(\bm{x})\|^{2}+\Delta_{\bm{x}}f_{\theta}^{(t)}(\bm{x})\big{]}.$
(21)

When $t$ is small and by the first-order Taylor approximation

$f_{\theta}^{(t)}(x)=f_{\theta}(x)+t\big{[}\frac{d}{\mathrm{d}t}f_{\theta}^{(t)}(x)\big{]}|_{t=0}+{o}(t),$

the corresponding DCD objective becomes

$\displaystyle t\mathcal{L}_{DCD}^{(VE)}(\theta)=$
$\displaystyle\mathbb{E}_{\bm{x}_{0}\sim p_{d},\bm{x}_{t}\sim
p(\bm{x}_{t}|\bm{x}_{0})}[f_{\theta}^{(t)}(\bm{x}_{t})]-\mathbb{E}_{\bm{x}_{0}\sim
p_{0}}[f_{\theta}(\bm{x}_{0})]$
$\displaystyle=\mathbb{E}_{p_{t}(\bm{x}_{t})}[f_{\theta}^{(t)}(\bm{x}_{t})-f_{\theta}(\bm{x}_{t})]+\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})]$
$\displaystyle=\mathbb{E}_{p_{t}}t[\frac{d}{\mathrm{d}t}f_{\theta}^{(t)}(\bm{x}_{t})]|_{t=0}+\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})]/t$
$\displaystyle=\mathbb{E}_{p_{t}}\frac{1}{2}\bm{G}^{2}(0)\big{[}\|\nabla_{x}f_{\theta}(x)\|^{2}+\Delta
f_{\theta}(x)\big{]}$
$\displaystyle+\mathbb{E}_{p_{t}}[f_{\theta}(\bm{x}_{t})]-\mathbb{E}_{p_{d}}[f_{\theta}(\bm{x}_{0})].$

###  A.11 Backgrounds on Skilling-Hutchison trick

Skilling-Hutchison’s (SH) [35] stochastic trace estimation trick is a commonly
used solution for efficient computation of trace of the Jacobian matrix for
high-dimensional problems. In our work, we adapt the SH trick to estimating
the trace of Jacobian which appears in equation 12. More precisely, we aim to
compute the trace of the Jacobian term

$\displaystyle\Delta_{\bm{x}}f_{\theta}(\bm{x})\coloneqq\nabla_{\bm{x}}\bm{s}_{\theta}(\bm{x}),$
(22)

where $\bm{s}_{\theta}\coloneqq\nabla_{\bm{x}}f_{\theta}(\bm{x})$ is the score
function of the EBM. The SH estimation uses a stochastic quadratic form to
estimate the trace term, i.e.

$\displaystyle\nabla_{\bm{x}}\bm{s}_{\theta}(\bm{x})=\mathbb{E}_{\epsilon\sim
p_{\epsilon}}\epsilon^{T}\nabla_{\bm{x}}\bm{s}_{\theta}(\bm{x})\epsilon=\mathbb{E}_{\epsilon\sim
p_{\epsilon}}(\epsilon^{T}\nabla_{\bm{x}}\bm{s}_{\theta}(\bm{x}))\epsilon.$
(23)

The distribution $p_{\epsilon}$ is assumed to be isotropic, i.e.
$\mathbb{E}_{\epsilon\sim p_{\epsilon}}{\epsilon\epsilon^{T}}=\mathbf{I}$. The
multivariate Gaussian distribution is a usual choice for $p_{\epsilon}$. The
vector-Jacobian-product term
$\epsilon^{T}\nabla_{\bm{x}}\bm{s}_{\theta}(\bm{x})$ is efficient to implement
with deep learning computation framework such as PyTorch with $\mathcal{O}(1)$
memory costs. More precisely, for a data $\bm{x}$, we first compute the score
function $\bm{s}_{\theta}(\bm{x})$ of the EBM by automatic gradient
computation functions of deep learning frameworks such as PyTorch. Then we
randomly sample a Gaussian vector and compute the Jacobian-vector product of
$\bm{v}^{T}\bm{s}_{\theta}(\bm{x})$. After that, we calculate the final
quadratic form
$\bm{v}^{T}\bm{s}_{\theta}(\bm{x})\bm{v}=(\bm{v}^{T}\bm{s}_{\theta}(\bm{x}))\bm{v}$.
However, Though the Skilling-Hutchison trace estimation trick can alleviate
the non-linear memory cost problem, frankly speaking, the DCD consumes more
GPU memory than MCMC-based methods. From this point of view, the DCD can be
understood as a method that trades memory costs for computational efficiency
when training EBMs.

###  A.12 Algorithm for training time-dependent EBM with DCD-VE

Input: dataset $\mathcal{D}=\left\\{x_{i}\right\\}_{i=1}^{n}$, time-dependent
EBM $f_{\theta}(x,t)$, diffusion process $(F,G)$, perturbation time $\delta$,
end timestamp $T$, mini-batch size B.

while  _not converge_ do

Sample time step $t\sim Unif[0,T]$,

Sample mini-batch uniformly
$\\{x_{i}^{(0)}\\}_{i=1}^{B}\sim\mathcal{D},i=1,..,B$,

Diffuse data sample with $x^{(t)}_{i}\sim p(x^{(t)}_{i}|x^{(0)}_{i})$,

Calculate DCD objective $\mathcal{L}_{DCD}(\theta)$(equation 12) with data
samples $\\{x^{(t)}_{i}\\}_{i=1}^{B}$,

Update $\theta$ with gradient decent according to minimize
$\mathcal{L}_{DCD}(\theta)$.

end while

return _$\theta$._

Algorithm 1 Training time-Dependent EBM with DCD

The available objective $\mathcal{L}_{DCD}(\theta)$ can be
$\mathcal{L}_{DCD}^{VE}(\theta)$ or $\mathcal{L}_{DCD}^{VP}(\theta)$ as
proposed in previous sections.

##  Appendix B More on experiments

###  B.1 Experiment Details on 2D Synthetic Modeling

#### Datasets.

We train EBMs on seven 2D datasets: Swissroll, Circles, Rings, Moons,
8Gaussians, 2Spirals and Checkerboard. The code to generate the dataset is
adapted from the open source codebase222https://github.com/wgrathwohl/LSD.

#### Model architecture

We use the multi-layer perceptron (MLP) with 4 layers and 300 hidden units in
each layer as the implementation of the energy-based model. We use the
Gaussian Error Linear Units (GELU) [41] as the activation function.

#### Hyper-parameters for DCD-VE.

We use the one-step DCD-VE (equation equation 12) for implementation. We use
$t=0.0005$ and $G(0)^{2}=1$. We train all models (with different methods) with
the same hyper-parameters: the optimizer is Adam optimizer with
$\beta=(0.9,0.99)$. The batch size is 1000, the learning rate is 0.001 and the
number of training iterations is 5000. For ablation training methods, i.e. CD
and PCD. For CD, we use 0.001 to be the step size of Langevin dynamics. The
number of iterations of the Langevin dynamics is set to be 10. For PCD, we use
a replay buffer with a size of 10000. The Langevin dynamic step size is set to
be 0.001 and the number of MCMC steps is 20. We set the update frequency of
the replay buffer to be 5%, which follows the setting of [8].

#### Evaluation metric.

We compute the score-matching loss over the training data as the evaluation
metric. The score matching loss is defined with

$\displaystyle\operatorname{L}(\theta)\coloneqq\mathbb{E}_{\bm{x}\sim
p_{d}}\bigg{[}\frac{1}{2}\|\nabla_{\bm{x}}f_{\theta}(\bm{x})\|_{2}^{2}+\Delta_{\bm{x}}f_{\theta}(\bm{x})\bigg{]}.$
(24)

So the smaller the SM loss is, the better the learning performance of the EBM.

###  B.2 Details on image denoising

In this experiment, we train EBM with CD and DCD-VE on four image datasets for
denoising: CIFAR10, SVHN, MNIST, and the FashionMNIST datasets.

#### Model architecture.

We use the Wide ResNet [56] with the Sigmoid-weighted Linear Units (SiLU) [57]
activations and no normalization as the implementation of the energy-based
model. For MNIST and the FashionMNIST model, we set the depth to 16 and the
widen factor to 8. For the CIFAR10 and SVHN datasets, we set the depth to 28
and the widen factor to 10.

#### Training details.

We first pre-process the data to scale the range of an image to $[-1,1]$. In
order to let the EBM learn the denoising ability of data samples, we pre-
process the training data by adding a Gaussian noise of amount $\sigma=0.3$.
We use the Adam optimizer [58] with $\beta_{0}=0.9$ and $\beta_{1}=0.99$ and
learning rate $0.0002$. For the DCD-VE training algorithm, we set the
diffusion strength $t=0.018$ and $G(0)^{2}=1$. To make a fair comparison, we
set the step size of the Langevin dynamics also to be $0.018$. For CD, we use
one-step of Langevin dynamics for implementing CD.

#### Evaluation metric.

To evaluate the denoising performance of trained EBM, we use the trained EBM
to denoise noisy images which are added Gaussian noise with three levels:
$\sigma=0.3$, $\sigma=0.6$ and $\sigma=0.9$.

###  B.3 Details on image generation

We train time-dependent EBM with the EDM [31] forward diffusion which is a
special instance of VE diffusion equation 9, for which the $g(t)=t$.

Samples of $\bm{x}_{t}$ are cheap to obtain by adding Gaussian noise to data
samples $\bm{x}_{0}\sim p_{d}$. We randomly choose a time
$t\sim\operatorname{LogNormal}(t;-1.2,1.2)$ following the same setting as the
EDM model and draw samples with

$\bm{x}_{t}=\bm{x}_{0}+\sigma(t)\epsilon,\epsilon\sim\mathcal{N}(\epsilon,\bm{0},\bm{I}).$

Here $\bm{x}_{0}\sim p_{d}$ denotes a data sample and $\epsilon$ is a standard
Gaussian vector of the same size as $\bm{x}_{0}$. Then we slightly diffuse
$\bm{x}_{t}$ to $\bm{x}_{t+\Delta t}$. This can be done by adding another
Gaussian noise of variance $\sqrt{\sigma(t+1)^{2}-\sigma(t)^{2}}$ and
calculate DCD with $\bm{x}_{t+\Delta t}$ and $\bm{x}_{t}$.

#### Network architecture.

We adopt a UNet encoder from the VP architecture of EDM model [31]. We add an
additional SiLU non-linearity to the layer before the last pooling layer.

#### Sampling Method.

We adapt the Heun sampling algorithm from Karras et al. [31] for sampling from
time-dependent EDM. We discretize the noise levels from 0.01 to 80.0 to 18
time-stamps with the same strategy of Karras et al. [31].