Moment Matching Denoising Gibbs Sampling
Authors: Mingtian Zhang
Alex Hawkins-Hooker
Brooks Paige
David Barber
Venue: NeurIPS 2023 poster (2023)
Keywords: denoising score-matching
gibbs sampling
diffusion model

Abstract: Energy-Based Models (EBMs) offer a versatile framework for modelling complex data distributions. However, training and sampling from EBMs continue to pose significant challenges. The widely-used Denoising Score Matching (DSM) method for scalable EBM training suffers from inconsistency issues, causing the energy model to learn a noisy data distribution. In this work, we propose an efficient sampling framework: (pseudo)-Gibbs sampling with moment matching, which enables effective sampling from the underlying clean model when given a noisy model that has been well-trained via DSM.  We explore the benefits of our approach compared to related methods and demonstrate how to scale the method to high-dimensional datasets.

Moment Matching Denoising Gibbs Sampling
Mingtian Zhang' Centre for Artificial Intelligence University College London m. m.zhang ucl.c....uk
Alex Hawkins-Hooker Centre for Artificial Intelligence University College London a. a.hawkins-hookerG cs. as..uk.ac.uk
Brooks Paige Centre for Artificial Intelligence University College London pppige@uul..cc
David Barber Centre for Artificial Intelligence
University College London
david.barber @ @ucc.uucl.uu. acl.aac.uk
Abstract
Energy-Based Models (EBMs) offer versatile framework for modeling com- plex data distributions. However, training and sampling from EBMs continue pose significant challenges. The widely-used Denoising Score Matching (DSM) method [41] for scalable EBM training suffers from inconsistency issues, causing the energy model to learn 'noisy" data distribution. In this work, we propose efficient sampling framework, (pseudo)-Gibbs sampling with moment matching, which enables effective sampling from the underlying clean model when given a "noisy' model that has been well-trained via DSM. We explore the benefits our approach compared to related methods and demonstrate how to scale the method to high-dimensional datasets.
1 Energy-Based Models
Energy-Based Models (EBMs) have attracted a lot of attention in the generative model literature 25
46. (5) EBMs are type of non-oooall probabilistic model that determines the probability density function without a known normalizing constant. For continuous data x, the density function of an EBM is specified as (x(x) exp( - (f())/// where the fe(z) is nonlinear function with parameter and ((0 = exp((fo(z) dæ is the normalization constant that is independent of r. The energy parameterization allows for greater flexibility in model parameterization and the ability to model a wider range of probability distributions. However, the lack of a known normalizing constant makes training these models challenging. We start by giving brief introduction of how estimate in EBMs and refer the reader to [37 for a detailed overview of different training techniques for continuous EBMs.
Likelihood based training: A classic method to learn is to minimize the KI. divergence between the data distribution pa and the model density qe. which is defined
KL(pallga) f pa(x) log qe (z) dr
/ pa(x) d
log Z(0),
(I)
where we use = to denote the equivalence up to constant that is independent of The integration of pa(ar) be approximated by Monte Carlo with the training dataset Xiaai fx1,.... .NN) pa(x); in this case, it is equivalent to the maximum likelihood estimate (MLF) [5] However, for EBMs, minimizing the KI. divergence requires the estimation of Z()), which intractable for nonlinear (x) defined by a neural network. Various methods have been proposed to alleviate the
"This work was partially done during an internship Huawei Noah's Ark Lab.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

intractability by introducing techniques like Markov chain Monte Carlo (MCMC) 13 26 8
or
adversarial training [88 47 6
based training: Alternatively, which is defined
16 proposes to minimize the Fisher divergence to learn 0,
FD(pallg))
(x)1l
(2)
where use (z( to denote the score function of distribution pi sp(xx) = Ve p(r). Under certain
regularity conditions, the Fisher divergence is equivalent to the score-matching (SM) objective [6
1 FD(pulllee> 2 pa(x)) (11s (z)1 + Tr(V S00 (x))) dr,
which does not require estimation of the intractable Z(0). However, this objective needs to caleulate the Hessian trace V Vz (z) - fo(z) 2)) in every gradient step during training, which is computa- tionally expensive and does not scale to high dimensional data requires approximation 38 In this paper, we will focus on another training method, denoising score matching [41), which overcomes the tractability and scalability issues mentioned above, and is introduced in the next section.
1.1 Denoising Score Matching
For the target data density Pa(x), a noise distribution p(o)2) N(x, o²I) is introduced to construct a noised data distribution (²) J (a(x)) (((())((()))) dz. Denoising score matching (DSM) [4] minimizes the Fisher divergence between the noised data distribution (ā) and an energy-based model (ū) = exp(-fo(( (0), with
/ ((t))lsp.(t) sne (f)l
FD(pallão) 2
/ p((l=)pa(=)) V logp(()>)
2
Sqo (i) di dx
2 // p(F(x))a(())
1
i - + s5g (ū)
di dæ
(4)
d²
where the last equation is due to Vz V=lo(((ū(r) being tractable for the Gaussian distribution p(7(xr).
Compared to the KL SM objectives, the DSM objective is scalable and well-defined when the data
distribution is singulat [49 and can alleviate the blindness problem of score matching [35] 45 51
On the other hand, there is a notable disadvantage associated with the DSM objective: for fixed a >> O, the DSM objective is not a consistent objective for learning the underlying data distribution Pa since FD( 0
Pa = 90 * Pa- A common solution is anneal training. However, Equation 4 is not defined when a = O since the division in Equation
0 during
4 will make (ã - =x//o² unbounded, which results in an inconsistent objective. Annealing a increases the variance of the training gradients (37 42). which makes the optimization challenging in practice.
To overcome the challenges, propose alternative data generation scheme: we use DSM with fixed σ > 0 to train a "noisy" energy model and then construct a sampler which targets the underlying 'clean' model. Specifically, contributions are summarized follows:
We demonstrate that for an EBM that learns noisy data distribution, there exists unique underlying clean model which recovers the true data distribution.
We introduce Beeudo-Gibbs sampling scheme incorporating analytical moment- matching approximation of the denoising distribution. This allows us to sample from the underlying clean model without requiring additional training.
We illustrate how to scale our method for high-dimensional data and demonstrate the generation of high-quality images using only single level of fixed noise. Furthermore, showcase the application of our proposed method in multi-level noise scenarios, closely resembling a diffusion model.
"The singular distribution is absolutely continuous (a.c.) with respect the Lehesgue measure, thus doesn't 't allow a density function 40 p. 172). A typical example a data distribution supported on a lower- dimensional manifold. In this case, the KL divergence is ill-defined and cannot be used to train the models. When using DSM, the distribution after Gaussian convolution would always be a.c. thus can be a valid training objective, see Zhang et al. 49 or Arjovsky al. [[] for detailed introduction.

---

2 Clean Model Identification
For fixed >> O, DSM only learn "noisy' data distribution (#) in the ideal case where the Fisher divergence is exactly minimized, since FP ( = 0 - Pa do o Pa- In this case, the following theorem shows that there exists 'clean' model that implicitly defined that learns the true data distribution.
Theorem 2.1 (Existence of the underlying clean model for optimal @o (3)). When the Fisher diver- gence goes to 0, FD (6 0 + Ps qo. there exists an unique underlying clean model g(x) such that qo (>) q(( p .)) dr and ((x) pa(z).
See Appendix A.1 for proof. This theorem shows that despite training an EBM on noisy data, there is an implicit model within that can recover the true data distribution. Therefore, instead of annealing the noise σ - 0 to recover the true data distribution, we will demonstrate how to directly sample from the implicityy-defined clean model given the noisy energy-based model in the next section.
We want to highlight that the 'perfect fit' assumption, i.e. achieving FD [āo) - O, may not hold for complex data distribution Pa or underpowered EBMs. Therefore, we provide general sufficient conditions for the existence of the clean model for an imperfect EBM in Appendix A.2
2.1 Gibbs Sampling with Gaussian Moment Matching
Given a well-trained noisy energy-based model (+) pa(T), the clean model has the form
g(x) - p (T)dū
(5)
where the denoising distribution can be written as p(z)ž) & p())=))(Z)) We notice that, since the noise distribution p(o|x) Mxx I2 is known, Gibbs sampling scheme can be constructed to sample from the underlying clean model if we know the denoising distribution p(z|ž), with
- p(ā|= = xx-1),,
p(r|f #x-1).
(6)
where the initial sample Zo (x)) can be drawn from a standard Gaussian pa(a)) = (N 7). However, as the denoising distribution p(x|z) usually intractable for complex Pa. we propose an analytical Gaussian moment matching approximation of p(z(ā).
Denote the mean and covariance of p(x(ā)
((ū) = (T)
= (x))=
(a)p(()))
The classic Gaussian moment matching method (2 specifies a Gaussian approximation p(a) (t) R N (H) Σ())). which matches the first and second moment of p(z[()) ( When d = Pd. the first
mean of the denoised distribution has a well-known analytical form
29
((ū) i
= ā + o²s sqy
(8)
we include the derivation in Appendix
A.3 Using this identity, we can rewrite Equation as
//
FD go)
20 p(a)x)pa(z)) ||2 m()) dā
where we can see that the Fisher divergence only depends on μ()). Since FD( = Pa (Theorem 2.1 the function (x) fully characterizes the distribution q. Therefore, m(() and p( 2:) (xx o² I) can provide sufficient information to determine p(x (ā) & q(z)p(ā(z). As consequence, the following theorem shows that the covariance function also be analytically derived.
Theorem 2.2 (Analytical Covariance Identity). Given a clean model q(x) such that pa(ā) with p(((z) N(z,0²1), the m(z) and E(t) of the p(z(z) x (()))(;(x) has the following relations
[g(z)p())))) dz (z())
[ o'Vzu(z) V (e) +
(10)
See Appendix A.3 for proof. This analytical covariance identity can be seen as a high-dimensional generalization of the 2nd-order Tweedie's 's Formula 9 29 Therefore, the analytical full-covariance moment matching approximation can be written
p(z|+)= ~ A(# + o log qo ( V log (+) + 021).
(11)
3

---

We want to highlight that since the Gaussian moment matching is only an approximation of p(z|z) the sampling scheme in Equation 6 is a 'pseudo' Gibbs sampler unless the true p(z[ā) is also a Gaussian distribution which is not true for general non-Gaussian pa- However, since uu(ā) and p(0[z) Max T² are already sufficient to specify p(æ](), it should be possible to derive expressions for higher-order moments which themselves involve only uuã) and we leave this to future work. To our knowledge, the z-conditioned full covariance Gaussian moment matching approximation to p(x[ā) has not been derived previously. In the next section, we briefly discuss the connections between our method and other related approaches.
2.2 Connection Covariance Learning Approaches
Bengio et al. 44] proposes to approximate the true posterior p(x[z) oc pa(z)p(ūjz)) with variational distribution qo ( ()) The parameter 9 is then learned by minimizing the joint KL divergence
1 /f pa(z)p(((z)) log go(x) dē dz,
= --
(12)
where (x( = f pa(z))(()))) dr. The joint KL divergence in Equation 12 encourages qo (z) ( to match the moments of the true posterior pt (i), and defines an upper bound of the marginal KI.
(L (()((()))((()()))))))))))))()))) Igo (< (x(=)) (z))) 2 (L(pa(())llox))),
(13)
where the model is implicitly defined as the marginal of the joint qa(x) = f (()))) When qa [0] is a consistent estimator of p(z[ this asymptotic distribution of the Gibbs sampling will For continuous data, the variational distribution (T) Ex)),, where the mean (e(-)) and the covariance E-(-) are parameterized by neural networks. We note that the only difference between the KL and DSM objective (Equation P is that the KL objective additionally learns the covariance. We thus show that the optimal covariance under KI. minimization is the proposed analytical covariance.
converge the true data distribution (x)
is chosen as a Gaussian distribution go (z)(²)
Theorem 2.3 (Optimal Gaussian Approximation). Let p(ā/zz = NN o² I) and assume Gaussian distribution qo (x(z) = (11, (i)( then the optimal qº such that
q-' = arg min (c)p((( 9
(14)
has the mean and covariance with the form
M T = p(a))))
[=)
= o² VHC
(15)
see Appendix A.4 for proof. Therefore, when the optimal mean function is learned (=( = (z).
the optimal Σ ( can be analytically derived, making the learning of Eo( redundant. In addition to the training inefficiency caused by more parameters, the amortized covariance network may suffer from poor generalization [50). Moreover, the KL objective is also not well-defined for learning data distributions which lie on a fow-dimensional manifold, e.g. MNIST, see Section 4 for a detailed discussion. In this case, the learned (T) may be a degenerate matrix, making the Gaussian density function (xx())) ill-defined 49 which impedes the training, see Figure 5 for an example.
Paper 23 proposes higher-order score- matching loss to simultancously learn both the first order score V log go (ō) and the second order score V log qo (z). However, our findings indicate that the mean function (=)) (or the first order score V= log p(:)) already contains all the moment information of the underlying true distribution Pa, and the optimal moment can be derived using the mean function. Therefore, learning the second-order score is redundant and may lead to sub-optimal inference.
2.3 Connection to Analytic DDPM
The recent paper Bao et al. 2 considers a constrained variational family = N(uo(T), in the context of diffusion model and derive the optimal og as
-
1 (Tr (Cov ((CVv((11))1t)) (41) pa(a))
0 = arg min (L(((E
=
(16)
'For example, when (a()) = NN(ua, ou), the true posterior p(x[ā) x p() will be a Gaussian with mean u²²) = (o² ua + 000)/(0² + oz)) and variance ( = 000² /(aa + +z)). we can verify that (() o² (z)

---

which can also be rewritten using the score function
-² =²= (d (i)112) pa(2)
(+)]]
(17)
In Appendix B we provide detailed derivation show how this approximation can be linked our method using the Fisher information identity [10]. This approximation has two potential limitations: first, compared full covariance moment matching. the assumed isotropic covariance structure may be insufficiently flexible to capture the true posterior; second, the covariance is independent of ž.
The second assumption only holds when μ(()) a linear function of a (e.g. when Pa( is Gaussian) and does not hold for other non- Gaussian Pz(z). Therefore, our ï-dependent full-covariance approximation offers a more ver- satile approximation family, which ultimately results in a more precise estimation. However, in certain applications such as accelerating the sampling procedure of diffusion model [2], is advantageous to use -7-independent isotropic covariance due to inexpensive estimation. On the other hand, our ë-dependent covari- ance necessitates the computation of the Hes- sian for each I. making it inefficient for high- dimensional data. Section we will explore approaches to mitigate this limitation.
S pa(x)p(ūlz) da
Pa(
Pa(
(a) Gaussian mixtures visualizations
100000
(b) Density of p(ar ( (blue) and four ø' points (green)
2.4 Posterior Approximation Comparison
(c) Analytical full covariance moment matching
We now consider a toy example to compare the three denoising posterior approximations discussed above. Let pa(x) be Mixture
of Gaussians (MoG) (a( Σ k= (x) whose components 91:4 are 2D Gaussians [-1,
(d) Analytical isotropic covariance moment matching
with means
[1] and isotropic covariance 7², with 0g = 0.2. The noise distribution is p(][2) = W(z,a21) with 0.2, so (()) = p(=)p(() da is an MoG with the same component means and diagonal covariance (a) T²) see Figure for a visu- alization. In this case the true posterior p(r[ū)) does not allow tractable form. Fortunately, given a noisy sample j' and an evaluation point ar, we can evaluate the true density p(a) Z') us- ing Bayes rule: p(a'p |
(c) Learned diagonal covariance
Figure I: Figure (a) shows the clean data distribu- tion pa(#)) and the corresponding noisy distribution ()))). Figure (b) shows 4 conditioned samples in the noisy space. Figures (c. d, e) visualize the true posterior p(x(ū) (green) and three posterior approximations (orange). We find that only the proposed ē-dependent analytical full-covariance moment matching can capture the variance of the true posterior, whereas the other two methods un- derestimate the variance.
K') p(f' x
Figure I shows the true pos- teriors given four different x' where we use grid data in z-space to visualize the density.
(pa((')/Pa(3)
To train the model, we sample 10,000 data points from Pa our training data. For the Kraiineeed
Gibbs sampler described Section 2.2 we use network with hidden layers with 400 hidden units, Swish activation 28 and output size to generate both mean and log standard deviation of the Gaussian approximation.
the moment-matcl matching Gibbs sampler (including both full and isotropic covariance), we use the same network architecture but with output size 1 to get the scalar energy and DSM as the training objective. Both networks are trained with batch size 100 and Adam [19 optimizer with learning rate 1x 10 for 100 epochs. For the ř-independent isotropic covariance, use the Monte Carlo approximation to estimate the variance with 10000 samples from pa()).
"Since when ((z) is a linear function of Z. using Theorem 2.2 we have ()) = o VVVap(Z) will not depend F. see also footnote for an example.
5

---

Figure I visualizes the approximations the denoising posterior p(z|ă) es- timated by each of the three meth- ods described in the previous sections. We surprisingly find that although the KL objective in Equation [2 encour- ages (z] to match the moments of p(z[ū), the learned covariance in Figure lc still underestimates the variance of the posterior: This shows the redundancy of covariance learning degrade the variational approximation performance. Additionally, the ž-independent covariance fails to account for the relative positions of 2' and lacks the ability to predict the posterior's 's elliptical shape due its isotropic nature. In contrast, t-dependent full covariance approximation overcomes these limitations, enabling more accurate predictions that capture the intricate geometry of the posterior distribution.
Table 1: MMD evaluations of single chain
Data
Learn diag.
Analytic iso.
Analytic full
MoGs
0.929 + 0.343
0.724 + 0.361
0.305 = 0.141
Rings
0.364 + 0.044
0.006 = 0.002
0.005 = 0.001
Roll
0.053 + 0.011
0.030 = 0.001
0.016 = 0.002
We then use the estimated posterior to conduct (pseudo) Gibbs sampling to generate sample
Specifically, we initialize the first sample 20 (0, 0.1) and run one Markov Chain with 10,000 time steps generate 10,000 samples. mixture of Gaussian datasets, also train and generate samples from the 2D Swiss roll and two-ring datasets. For numerical evaluation, we calcu- late the Maximum Mean Discrepancy (MMD) [2] between 10k samples gen- erated by single-chain Gibbs sam- pler and 10k samples from the train- ing dataset respectively. The kernel insides MMD is a sum over Gaus- sian kernels with bandwidth ranging over [2-2 2-1 20 2 2²). The MMD) results (including both mean and std) are calculated using 5 random seeds. find that Gibbs sampling with the proposed analytical full covariance achieves the best results; numerical results are in Table I with visual comparison in Figure
True data
(b) Leam. diag.
Analytic iso.
Analytic full
addition to the
(e) True data
(l) Learn. diag.
(g) Analytic iso. (h) Analytic full
True data
(j) Learn. diag. Analytic iso. Analytic full
Figure 2: Samples from single chain Gibbs sampling
3 Scalable Implementations for Image Data
Scalable Diagonal Hessian Approximation As we discussed in Section 2.3 the proposed full covariance Gaussian approximation in Equation [0reqqurrs calculating D DDDDHesiiii V log (()) for each with size D. which brings both memory and computation difficulties for high-dimensional data. A naive diagonal Hessian method (only using the diagonal entries in the Hessian) will address the memory bottleneck but still needs times backward passes for the exact computation of the diagonal term [22]. this paper, the following diagonal Hessian approximation 3.
Diag(H) ~ 1/S Σ ==1 U. O HU..
(18)
where ~ p(u) (u) is Rademacher random variable with entries +1 and O denotes the element-wise produc This estimator will converge to the exact Hessian diagonals when S + oo [ The computation for each v, be computed by two that our ž-dependent diagonal moment matching approach provides a comparable level of flexibility to the variational method proposed in 4 while eliminating the need for additional training of the
passes. It is worth emphasizing
*The code of the experiments can be found in httph:tt///// / htt///gitttt httpgttttuub..../.gg wtt.//tttoorroo///ot..oooo This estimation should be distinguished from the Hutchinson's Trace estimation 15 Tr(
Neur
S of Hus where dot dot-product is used between v, and Hu..

---

0344
4 9 2
4 2 9 7
$ 28 8
6 ~ 2 0
6 780
2 2 523
2 - 0
69
69387 3
L2 9
0 5
02095
a026 8
4 2
40
6 2
(a) Diagonal ž-dependent covariance moment matching (Ours)
6
3
) A 856
9
0
43
5 4
6
5
30
2
5
0 2
1 10520
2
3
2
3 29
3 90
2
06
K 2
(b) Isotropic ž-independent covariance moment matching 2
64
5
D
से
22
3 2 R 2
8
14
2 2 M
2 a
2
4
6 y 0
O 2
& 2 $
4
6 8 2
θ d
9
(c) Diagonal covariance learned by KL minimization
Figure 3: Figures (a,b,c) show the MNIST experiment comparisons, where we compare samples generated by pseudo-Gibbs sampling with three different g(x(ā). We plot samples from 25 indepen- dent Markov Chains with E {0, 1, 10, 20) time steps. We can find the samples generated by the proposed analytical covariance moment matching with diagonal approximation achieved the best sample quality.
diagonal covariance. Furthermore, method remains more flexible than the isotropic ă ē-independent moment matching method proposed by 2].
Energy or Score Parameterization For the full-covariance moment matching in Equation 10 we require V2 log go (ā) to symmetric to obtain valid Gaussian approximation. However, we learn the score function V log p()) = so (i) using a network so(+) (-) R - TR.,. its Jacobian not guaranteed to be symmetric. In this case, we follow [34] and directly parameterize the density function (o(-) with a neural network fo(-) R - IR and let the score function V log ())
(G)).
This can be obtained by AutoDiff packages like PyTorch (27]. and this parameterization guarantees (Z) to be symmetric. We also notice that when using the diagonal Hessian approximation (Equation 18) we only need entries in Diag H)) to be positive in order to obtain valid Gaussian approximation. this case, the score parameterization remains applicable and offers more efficient training compared to the energy parameterization. Therefore, the combination of full/diagonal covariance and energy/score parameterization provides a tradcoff between flexibility and inference speed, allowing for flexible approach while maintaining computational efficiency during training.
4 Image Generation with Single Noise Level
We then apply the proposed method to model the grey-scale MNIST 21 dataset. We use the standard U-Net architecture [35|] 31 with single fixed noise level 0.5; the effect of varying is explored in Appendix I For the KL training objective, the output channel size is 2 to generate both mean
and log-std at the same time. For DSM training, we take the sum of the U-Net output obtain the scalar energy evaluation which also relates to the product-of-experts model described in (34. We train both networks for 300 epochs with learning rate xx10 x and batch-size 100.
7

---

As discussed in Section 1.1 the KL divergence not well-defined for manifold data distributions. This limitation becomes evident when working with MNIST, where the presence of constant black pixels in the boundary areas leads to rapid decrease towards in the variance of qo (2) (z) during training. Consequently, the likelihood value log (=(=()) tends to approach infin- ity, resulting in unstable training. In contrast, the DSM objective is well- defined for manifold data, providing stable training process even in the presence of such boundary effects. [reuvides visual comparison of the two training procedures, demonstrating the improved stability and effectiveness of the DSM objective in handling manifold data distributions.
(a) ä-dep diag (Ours) (b) ž-dep diag. (KL)
(c) j-ind iso.
Figure 4: Figures (a,b,c) visualize the covariance approxi-
oe), e NO, o² I) on 25 ø samples. We use sigmoid function to map the real value noise into grayscale pixels for the visualization.
mations q
For the sample generation process, calculating the full-covariance Gaus- sian posterior becomes challenging. We therefore apply the scalable diago- nal Hessian approximation described in Section 3 to approximate the diag- onal Gaussian covariance of p(=: We find that the estimated diagonal Hessian occasionally contains small negative values due approximation error: we, therefore, use the max(", e) function with e itivity of the diagonal covariance. The ž-independent isotropic covariance and the proposed i- dependent diagonal covariance share the same mean function.
2040
2000
280 1M0 1999 2000
(a) KL loss
(b) DSM loss
Figure 5: Training loss comparison of two objectives. We
plot the training loss every iteration during total 300 epochs.
ensure the pos-
We first visualize the covariance esti- mated by three different methods in Figure [ We use 100 Rademacher samples in estimating the diagonal Hessian (Equation T8p and 50,000 samples in estimating the isotropic variance (Equation] I7). We find that both ä-dependent diagonal covariance approximations can capture the pos- terior structure whereas the isotropic f- f-independent covariance is just Gaussian noise since the variance is shared between different digit and pixel locations. In Figure we plot the sample comparison for three methods.
(a) Data x
(b) sample
(c) 10 samples (d) 100 samples
Figure 6: Visualizations of the diagonal covariance q(
+ σε) with different number of Rademacher samples.
Since the isotropic covariance has the same variance in each dimension, the generated samples in Figure 3b contain white noise in the black background, whereas the proposed full-covariance sampler
Figure 7: Samples from three Markov chains. We plot the samples every 10 Gibbs steps.

---

Algorithm 1 Sampling with Langevin Dynamics
Algorithm 2 Sampling with the proposed pseudo Gibbs Sampling
Require: {or) 6, K
Initialize 2/1 fort + T do
Require: {or)E-o- K Initialize z1/
0t s 1-007/00 for k 1 to K do Draw z ~ ((,
for + T to do
for 4- 1 to K do
Draw xit ~ p(0++1 x1)
x k- +a k-1 000 22,2
Draw x ~ ( = =x(t1)
end for
end for 2/--1 end for return 2 + o0s0(z0, σo)
J/-11 c nr end for return = + o3s0(=B. σn)
t xx
can generate a clean black background Figure 3a On the other hand, the samples generated by the KL-trained Gibbs sampler (Figure 3cb have worse sample quality due to the unstable training.
We then apply the same method model the more complicated CIFAR 10 [20] dataset. We use the same U-Net structure as used in [36] and directly parameterize the score function rather than the energy function to speed up the training. The noise level is fixed at 0.3. We train the model using Adam optimizer with learning rate 1 1x10 and batch size 100 for 1000 epochs. We visualize the denoising posterior diagonal covariance Figure when using different numbers of Rademacher samples (Equation 18}. We observe that better covariance estimation can be obtained by increasing the number of samples. To balance efficiency and accuracy, we use a sample number of 10 in the subsequent Gibbs sampling stage. Figure 7 shows three independent Markov chains with the samples plotted every 10 Gibbs steps, which demonstrates that sharp images can be generated with even one fixed level of noise.
Limitation: In the CIFAR experiment, we observe mode collapse phenomenon when running multiple independent Markov chains for a longer time. This phenomenon is likely due to the small noise level a = 0.3, which prevents the sampler from exploring the full space, as commonly found with MCMC methods (30). This effect is visually represented in Figure S where we assess the Fréchet In- ception Distance (FID) values for 50,000 images sampled with varying numbers of Gibbs steps. Notably, the FID increases beyond 40 Gibbs steps, and visual evidence of mode collapse is observed (Figure 12). In the ensuing sec- tion, we will demonstrate the application of method to settings with multiple noise levels, an approach that may help mitigate the issue of mode collapse.
420
250
150
20
60
Figure FID evaluation with Increased Gibbs Steps. We can find the FID in- creases after 40 Gibbs steps.
5 Image Generation Using Multiple Noise Levels
The success of diffusion models and lessons from prior work on score-based generative models point to the importance of using multiple noise levels 14 355 multi-modal distributions. Intuitively, by learning denoise data at a range of noise levels, single network can learn both the fine and global structure of the distribution, which in turn allows for more effective sampling algorithms capable of efficiently exploring diverse modes (35]. We therefore propose to adapt the denoising Gibbs sampling procedure to sample from distributions corrupted with multiple noise levels. For this purpose, we use a coise-conditioned score network trained by Song and Ermon [36], who generated high-quality samples using a procedure inspired by annealed Langevin dynamics 43 This procedure involves generating samples from a sequence of distributions pT (#T), Po (foo). corrupted by progressively decreasing levels of Gaussian noise (parameterized via standard deviations oe, with or >,... ,> oo). At a given step the sequence, Langevin dynamics is used to sample from the corresponding noised distribution pe (2xe, using the
when modelling data with complex

---

score network (T σt) to approximate the gradient of the noised distribution. The outputs of this Langevin dynamies run are then used to initialize the same procedure at the next noise level, leading the sampling procedure to converge gradually towards the data distribution as the noise level tends to zero (i.e. (o)) R Pa(x). do R 0). The algorithm of the annealed Langevin dynamics with multi-level noise used in 35
36
is summarized in Algorithm
We show that the proposed Gibbs sampling scheme can be directly applied to a pre-trained score- based generative model as a drop-in replacement for Langevin dynamics MCMC in the generation stage, At each noise level, we use samples 2+41 from the previous noise level to initialise a Gibbs sampling chain targeting the marginal distribution at the current noise level pa (1z), which now plays the role of the 'clean' distribution in Equation [5] Therefore, the noisy distribution time step t is a Gaussian p( [ = Max,, o 1 - of). The optimal denoising distribution g(2, (T++1) is thus
a function of the level of noise at step & relative to the level of noise the previous step / + 1. For
generation efficiency, we employ 3 Gibbs steps at each noise level, using 3 Rademacher samples to approximate the diagonal Hessian. The sampling procedure is summarized in Algorithm 2 In Figure
9 and 10 we visualize the samples from models that are trained on CIFARIO and CelebA separately.
Further experimental details can be found in Appendix
C.2
For direct comparison with the results of [36] on CIFARIO, we retain the same schedule of
Table 2: CIFARIO Inception and FID Scores
noise levels used to generate samples with
Model
Inception
FID
Langevin dynamics. We generate 50000 sam-
NCSNv2 (Langevin [36
8.40 + 0.07
10.87
ples using this approach and report FID and
NCSNv2 (Gibbs, Ours)
8.28 + 0.07
10.75
Inception scores Table [ Our multi-level Gibbs sampling scheme produces samples of equivalent quality to the multi-level Langevin dynamics of [36). confirming its applicability to complex natural image data. The FID is also notably superior to that of the single-noise level Gibbs sampling, and the samples exhibit significant visual diversity (Figure 9 This underlines the importance of employing multi-level noise in our approach. Recent advances in sampling strategies for score-based models leveraging the framework of stochastic differential equations [39] have led to significant further improvements in generation quality as shown in Table 2 we leave the exploration of possible applications of our method to this framework to future work.
DDPM [4]
9.46 + 0.11
3.17
NCSN++ 39
9.89
2.20
Figure 9: CIFAR Samples
Figure 10: CelebA Samples
6
Conclusion
This paper focuses on addressing the inconsistency problem in training energy- based models (EBMs) using denoising score matching. Specifically, we identify the presence of an underlying clean model within a 'noisy" EBM and propose an efficient sampling scheme for the clean model. We demonstrate how this method can be effectively applied to high-dimensional data and showcase image generation results in both single and multi-level noise settings. More broadly, we hope our more accurate denoising posterior opens new avenues for future work on score-based methods in machine learning.
10

---

References
[1] M. Arjovsky, S. Chintala, and L. Bottou. Wasserstein generative adversarial networks. In International conference machine learning, pages 214-223. PMLR, 2017.
[2] Bao, C. Li, Zhu, and B. Zhang. Analytie-dpm: an analytic estimate of the optimal reverse variance in diffusion probabilistti models. arXiv preprint arXiv:2201.06500 2022.
[3] C. Bekas, E. Kokiopoulou, and Y. Saad. An estimator for the diagonal of matrix. Applied numerical mathematics, 57(11-12)12114-222) 2007.
[4] Y. Bengio, L. Yao, G. Alain, and P. Vincent. Generalized denoising auto-encoders generative models. Advances in neural information processing systems, 26, 2013.
[5] C. M. Bishop and N. M. Nasrabadi. Pattern recognition and machine learning, volume 4. Springer, 2006.
[6] A. J. Bose, H. Ling, and Y. Cao. Adversarial contrastive estimation. arXiv preprint arXiv:1805.03642, 2018.
[7] Y. Du and I. Mordatch. Implicit generation and modeling with energy based models. Advances in Neural Information Processing Systems, 32, 2019.
[8] Y. Du, S. Li, Tenenbaum, and Mordatch. Improved contrastive divergence training of energy based models. arXiv preprint arXiv:2012.01316, 2020.
[9] B. Efron. Tweedie's ' formula and selection bias. Journal of the American Statistical Association, 106(496);1602-1614, 2011.
[10] R. A. Fisher. Theory of statistical estimation. In Mathematical proceedings of the Cambridge philosophical society, volume 22, pages 700-725. Cambridge University Press, 1925.
[11] R. Gao, Y. Lu, J. Zhou, S.-C. Zhu, and Y. N. Wu. Learning generative convnets via multi-grid modeling and sampling. In Proceedings of the IEEE Conference Computer Vision and Puttern Recognition, pages 9155-9164, 2018
[12] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola. A kernel two-sample test. The Journal of Machine Learning Research, 13(1):723-773, 2012.
[13] G. E. Hinton, Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):177--1800 2002.
[14] J. Ho, A. Jain, and Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-660 6851,, 2020.
[15] M. Hutchinson. stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statistics-Simulation and Computation, 19(2):433 450, 1990.
[16] A. Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
[17] A. sarial score matching and improved sampling for image generation. In International Con- ference on Learning Representations, 2021. URL https://opanreviee wtt.eenrreviee.netforum.... nwt./frmufii eLfqM13z31q
R. Piché-Taillefer. I. Mitliagkas, and R. T. des Combes. Adver-
[18] Kim and Y. Bengio. Deep directed generative models with energy- based probability estima- tion. arXiv preprint 160i 03439, 2016.
[19] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 2014.
[20] A. Krizhevsky, G. Hinton, et al. Learning multiple layers of features from tiny images. 2009.

---

[21] Y. LeCun. The mnist database of handwritten digits. http://yann. lecun, com/exdbhmmist/, 1998.
[22] J. Martens, I. Sutskever, and K. Swersky. Estimating the hessian by back-propagating curvature. arXiv preprint arXiv: 220.66644 2012.
[23] C. Meng, Y. Song, W. Li, and Ermon. Estimating high order gradients of the data distribution by denoising. Advances in Neural Information Processing Systems, 34:25359 253699,99936,, 2021.
[24] T. P. Minka. Expectation propagation for approximate bayesian inference. arXiv preprint arXiv:1301.2294, 2013.
[25] J. Ngiam, Z. Chen, P. W. Koh, and A. Y. Ng. Learning deep energy models. In Proceedings of the 28th international conference machine learning (ICML-II), pages 1105-1112, 2011.
[26] E. Nijkamp, M. Hill, S.-C. Zhu, and N. Wu. Learning non-convergent non-persistent short-run meme toward energy-based based model. Advances in Neural Information Processing Systems, 32, 2019.
[27] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, 7. Vito, 7. Lin, A. Desmaison, L.. Antiga, and A. Lerer. Automatic differentiation in pytorch. 2017.
[28] P. Ramachandran, B. Zoph, and Q. V. Le. Searching for activation functions. arXiv preprint arXiv:1710.05941, 2017.
[29] H. E. Robbins. An empirical bayes approach statistics. In Breakthroughs in Statistics: Foundations and basic theory, pages 388-394. Springer, 1992.
[30] C. P. Robert, G. Casella, and G. Casella. Monte Carlo statistical methods, volume 2. Springer, 1999.
[31] O. Ronneberger, P. Fischer, and T. Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer-asiisted intervention, pages 234-241. Springer, 2015.
[32] T. Salimans and J. Ho. Should ebms model the energy or the score? In Energy Based Models Workshop-ICLR 2021, 2021.
[33] S. Saremi. On approximating V with neural networks. arXiv preprint 1rXiv:1910.10 12744. 2019.
[34] S. Saremi, A. Mchrjou, B. Schölkopf. and A. Hyvärinen. Deep energy estimator networks. arXiv preprint arXiv: 2018.
[35] Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019.
[36] Y. Song and S. Ermon. Improved techniques for training score-based generative models. In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, ed- itors, Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurlPS 2020, December 6-12, 2020, virtual, 2020. URL ftttps://proceedingg hetgis.nurris... ceeurips.cc/pppeeeeee /2000/hash/ 100099661111555799997776666
cobstract.html.
[37] Y. Song and D. P. Kingma. How train your energy-based models. arXiv preprint arXiv:2101.03288, 2021.
[38] Y. Song, S. Garg, J. Shi, and S. Ermon. Sliced score matching: A scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence, pages 574 PMLR, 2020.
[39] Y. Song, Sohl-Dickstein, D. Kingma, A. Kumar, S. Ermon, and B. Poole. Score-based gener- ative modeling through stochastic differential equations. In International Conference on Learn- ing Representations, 2021. //tttst:///ppnreeeeeesete
[40] T. Tao. An introduction to measure theory, volume 126. American Mathematical Society Providence, 2011.
12

---

[41] P. Vincent. A connection between score matching and denoising autoencoders, Neural compu- tation, 23(7):1611-1674, 2011.
[42] Z. Wang, S. Cheng, L. Yueru, Zhu, and B. Zhang. wasserstein minimum velocity approach to learning unnormalized models. In International Conference on Artificial Intelligence and Statistics, pages 3728-372 3738. PMLR, 2020.
[43] M. Welling and Y. W. Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th international conference machine learning (ICML-II), pages
68 1688, 2011.
[44] H. Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.
[45] L. Wenliang and H. Kanagawa. Blindness of score- based methods to isolated components and mixing proportions. arXiv preprint arXiv:2008. /0087, 2020.
[46] Xie, Lu, S.-C. Zhu, and Y. Wu. A theory of generative convnet. International Conference on Machine Learning, pages
PMLR, 2016.
[47] S. Zhai, Y. Cheng, R. Feris, and Z. Zhang. Generative adversarial networks as variational training of energy based models. arXiv preprint arXiv: 2016.
[48] M. Zhang, T. Bird, R. Habib, T. Xu, and D. Barber. Variational f-divergence minimization. arXiv preprint arXiv: 1907. 1891, 2019.
[49] M. Zhang. P. Hayes, Bird, R. Habib, and D. Barber. Spread divergence. In International Conference Machine Learning, pages 11106-11116. PMLR, 2020.
[50] M. Zhang, P. Hayes, and D. Barber. Generalization gap in amortized inference. Advances in neural information prucessing systems, 2022.
[51] M. Zhang, O. Key, P. Hayes, D. Barber, B. Paige, and F.-X. Briol. Towards healing the blindness of score matching. arXiv preprint arXiv:2209.07336, 2022.
13

---

A Proof and Derivations
A.1 Proof of Theorem 2.1
The existence is straightforward, since FD >= 0 - pu = go-, we can simply let g(x) = Pa(æ), which makes f (())p(((()) dx = (=))p)(( 2:) dz Pa- To show the uniqueness, denote density k(c) A o I),, go (z) and pz())) can be written convolutions
go(#) = q k pa(T) = Pa * k.
(19)
we then have
= pd 2o + q * = pa k + =
(20)
where F denotes the Fourier transform. Since the Fourier transform of Gaussian is also Gaussian, so F(k) > 0 everywhere, we have
Pa ςo- + - F(g) F F(pa) + q - Pa-
(21)
Therefore, q Pa is the unique distribution that makes pa qo- This technique has also been used to construct spread KL divergence (we denote as KL) [49], which defined as (L(pullgg = KL(ps * kllgo k) where k(e) = N (0, I2), to train implicit model qa. Different from the DSM situation, when (Lallll9)) O, the underlying model go = Pa is directly available, whereas the EBM trained by DSM learns to be the noisy distribution f0 Pa * k.
A.2 General Conditions Characterising the Existence of the Clean Model
In the previous section, we assume for flexible neural network parameterized fo, the energy energy-based based model (ā) exp(-f(z)) (())-(())//2()) trained by Equation cann recover the target noisy data distribution Çe- = Pa so there exists an underlying model q such that qo- q k and q = pa- This assumption is commonly used in the literature on score-based methods. For example, in the score-based diffusion for any data x RDD the score function V log qo (ă) is usually parameterized by a neural network NNo(.) IR + i. However, this parameterization cannot guarantee NN (ō) is a conservative vector field, or in other words, there doesn' 't exist a distribution (ō) such that V (ā) score estimation V log Pa( 2) V log @o (i) is implicitly assumed to allow an EBM interpretation.
models literature (35 14 2.
V = log go (ă) and V log q(ī) is symmetric [32]33] Therefore, perfect
However, the underlying clean model doesn't always exist for imperfect model go * Pa- We here provide the sufficient and necessary conditions which guarantee the existence of the underlying clean model.
Theorem A.1 (Necessary and Sufficient conditions for the existence of the underlying clean model.). For a model ğo with the convolutional noise distribution k(c) NO0 o²!), there exists an underlying model q such that q k q if and only if F(āe)/ is positive semi-definite Additionally, the underlying distribution q can be written as
q == (F(o)
(22)
where F-1 is the inverse Fourier transform. This theorem is a straightforward corollary of Bochner's Theorem [ However, for the energy model (Ø) exp (-fo it's difficult to design func- tioning family of f that satisfies the positive semi-definite condition and have the tractable score function at the same time | We thus leave the design of better energy function parameterizations as promising future direction.
"A continuous function f Rd - C is positive semi-definite if for all E N, all sets of pairwise distinct centers X = (#1...., ax) E Bd and all E C.V -E'", Σn= (ra - æ;)) 2 see 44 Definition 6]
Bochner's Theorem [ Theorem 6.6]: A continuous function f RRd + C positive semi-definite if and only if is the Fourier transform a finite non-negative Borel measure on Rd
Forr example, one can define noisy enargy-based model ğe = (2)) with fa(
log I (=)g(()) 1/o²|[0 - ==12) dr, which always allows underlying clean energy- based model (e(2) exp (z))()))Z(()) such that gz) = qe(z) & with k(e) = (100 I1)). However, the score function V: ( fofo()) intractable in this case.
14

---

A.3 Proof of Theorem 2.2
Derivation of the Mean Identity
We let (0( = f k(0)z))0(()) dř, where k(0/2) = NN o2I), we have
V=)o(T)
S Vakk(ā (x) do
(o
te (i)
(x())
1
/
(ā (ā- (( dz
qo
k
x (00l)))a(() (x)
V² log (2) +
dr (m)
9o (z)
where we define the model denoising posterior using Bayes rule go (x| second equality due to the following Gaussian distribution property
g( The
V 2
-G -
202
E
262
I I
k(ālz).
(23)
V
0²
2TT
x²
Derivations of the Analytical Full Covariance Identity
We have derived the mean identity
(ā) = (1) = V²V2 log qo (ū) + a.
(24)
Taking the gradient over in both side and scaling with 0², we have
o² )²V = ~ log (σ)) + I
(25)
We can also expand the hessian of the log go (E):
1
/ V
k(ālz ((
V log (o(
(E x)
dx
σ²
(8)
1(3))))0(x)
0² /
1
Vzk(Tlan go go
aago (z)
dr +
(# - x)
dx
σ²
(T)
(
V Vak(T 44 V log qo (Z) dr
log 2o (T) + +1 =
(i x))
()(T)
= =
xz (7- x) x))(f()))))) (1) + _
(x) qo ( (k)
(i- -x)
dx
(x)
To To ( + 002 = / (- (- + x² - (z)))(()))) go (x(T) dr
(x (0)))) (x)² (=|=)) = Eg
Therefore, we obtain the analytical full covariance identity.
(3) = 22
(26)
A.4 Proof of Theorem
Lemma A.2 (KL to Gaussian (2)). Let p(z) be distribution with mean pp and covariance Σ.p and ((x) = (Hg En), denote the differential entropy as H(p) =-f log p(x) da, we have
KL(pllg) KL Ep)lla) + (HA Ep)) H(p)
(27)
The proof can be found in 2 Lemma 2.
We then prove Theorem 2.3
Since p(E(=)pu(x) p(z)#))u()), where Pa(Ø)
J (a()))(ā)x) da, we have
( ((a)
(28)
15

---

Assume Gaussian distribution ((x]ā) (z), (a))and denote the mean and covariance of the true posterior are (z) and (z), then the optimal q"
9º arg min KLlotlllchpaae
(29)
arg min (KL60())))lo((
(30)
((3)
arg min (KL 22) (uN E,)) H(p(z(z)) (())
(31)
arg min (KL 21)1(((())))) p())
const..
(32)
Therefore, the optimal u(x]ā)= = N (M, (x) under the joint KL has the mean and covariance N6 ( = (p(T), (G) (T)
B Connection to Analytical DDPM
Paper 2
considers the constrained variational family qo (0(2)) (i and derive the
optimal o as
σ arg min KLll(((loopp(()))aaaaal))) - d (Tr (Cov (Covg(el=)())) )
(33)
which can also be rewritten using the score function
og = 0² d (1sspo (i)]1 pa(2)
(34)
To make deep connection, we can also plug our analytical full covariance (Equation 10) into Equation 16
2 +
V2 log go( Pu(a))
d
Pa( = =0 0a (1)m3))))) Ba(a)
Tr (s (z) (i)
(35)
d
d
which recovers Equation 17 where the first equality is due to the well-known Fisher information identity [10).
C Experiments
All the experiments conducted in this paper run single NVDIA GTX 3090.
C.1 Effeet of the Single Noise Choice MNIST
Figure II shows the samples generated by our method with the EBM trained with difference E {0.3, 0.5, 0.8} in the noise distribution p(ž|z), we can find the image quality also heavily depends on the choice of the noise scale and 0.5 achieves the best visual quality, we then use this hyper-parameter in the subsequent comparisons.
C.2 Multi-level Noise Details
For full details on the architecture and noise schedule used in the multi-level noise experiments in Section we refer to Appendix B of [36]. For our multi-level Gibbs sampling procedure, we used Gibbs steps each noise level and Rademacher samples for each diagonal Hessian computation. Following (36), we used total of 232 noise levels, distributed according to their proposed geometric schedule, and applied a final denoising step which the mean of the clean distribution conditioned on the final output of the sampling procedure is returned (the final output of the sampling procedure a sample from the noised distribution from the noise distribution at the smallest noise level). This denoising step previously found to improve FID scores [17] significantly.
16

---

C7529
77 2 0 3
8 5 808
7 4 22
5 0
5
3 2
3 2 306
6 1553
5 57818
3 88058
0 0
8 89068 8
9 6 3 8 2
(a) 0.3
(b) 0.5
(c) a 0.8
Figure 11: Sample comparisons with different a value.
Figure 12: Mode Collapse visualization of 25 Markov chains, we plot the samples every 20 Gibbs steps, we can find less modes are covered if we run the Gibbs sampling for a longer time.
17


Reviews:
[{"review_id": "8Hdkuu0Poj", "reviewer": "Reviewer_ddyX", "summary": "This paper proposes a sampling framework, i.e., (pseudo)-Gibbs sampling with moment matching, which aims at sampling from the implicitly-defined clean model given the noisy energy-based model well-trained by the denoising score matching method. It first provides a theoretical justification of the existence of the underlying clean model for a well-trained energy-based model. It then introduces the gibbs sampling framework based on the proposed analytical Gaussian moment matching approximation method for the denoising distribution. Experiments on several toy datasets and MNIST and CIFAR10 datasets are conducted to demonstrate the effectiveness of the proposed approximation and sampling method.", "questions": "Please see the weaknesses section.", "limitations": "The authors have mentioned the limitations in Sec 3. Scalable Implementations for Image Data and Sec 4. Image Generation with a Single Noise Level.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The paper provides theoretical justification of the existence of the underlying clean model of a well-trained noisy EBM, and discusses the existence of the clean model for an imperfect EBM. \n- It derives an analytical Gaussian moment matching approximation method for the denoising distribution and validates the design on several toy datasets, and discusses some scalable implementations of the proposed method for image data.\n- The presentation in general is clear.", "weaknesses": "- The full covariance moment matching requires calculating the Hessian matrix, which could be computationally expensive even with diagonal approximation methods.\n- The experiments are conducted on toy datasets and MNIST and CIFAR-10 datasets. More results on other image dataset could provide better intuitive understanding of the proposed method.\n- The performance of the proposed sampling method is close to the existing methods. The significance of the proposed sampling method is not very clear.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and valuable feedback that improves the quality of our work. We are encouraged by the positive comments regarding our theoretical contribution and presentation. \n\n### Hessian computation\nWe want to highlight that our method can directly obtain the optimal covariance by **only learning the mean function** during training. In contrast, the classic KL training Gibbs Sampling approach necessitates the learning of **both the mean and covariance functions** in the training time. This not only demands extra parameters and increases computational costs but also poses a risk of yielding a sub-optimal solution. Therefore, our method is both **theoretically optimal** and **more computationally efficient in the training stage** (our method doesn’t require any Hessian computations in the training).\n\nHowever, in the generation stage, we do agree with the reviewer that calculating the Hessian or even diagonal Hessian is computationally expensive for high-dimensional data, but it completes the following spectrum of **flexibility and efficiency trade-offs**:\n\n- Data-independent isotropic Gaussian approximation proposed [1]: This is theoretically the least flexible approach but is the most computationally efficient approach.\n- Proposed data-dependent Diagonal Gaussian approximation: More flexible than the data-independent isotropic method, but less computationally efficient.\n- Proposed data-dependent full covariance Gaussian approximation: This method represents the theoretical optimal point on the flexibility spectrum, but is the least computationally efficient of the three options.\n\nWe posit that these trade-off methods will expand the repertoire of ML toolboxes, catering to various applications, each with its distinct demands for accuracy and efficiency.\n\n\n### Additional Experiments\nWe appreciate the reviewer's suggestion to test our method on an additional dataset.  In response to this, we have extended our methodology to include the CelebA dataset, notable for its heightened dimensionality (64x64x3), see the attached PDF in the global rebuttal for the sample results. The generated samples further highlight our method's capability to produce high-quality human face images, other experimental details and numerical evaluations will be added to the revised version.\n\n### Contribution and Effectiveness\n We would like to clarify that the primary contribution of this paper is to unravel the theoretical optimal full-covariance Gaussian approximation of the denoising distribution. This is an open problem in the realm of the denoising score-matching method. The recent paper [1] preliminary derived an analytical form in a constrained data-independent isotropic  Gaussian approximation setting. Our investigation has led to a surprising discovery that this optimal approximation can be elegantly expressed in an analytical form for the full-covariance data-dependent moment matching, providing a solution that is both more flexible and accurate.\n\nWe've demonstrated that, in the single-noise level setting, our proposed method **outperforms** the classic denoising Gibbs sampling, and it does so with fewer parameters and lower training costs. For multi-level noise, our primary goal is to provide a proof of concept that the proposed Gibbs sampler can serve as a **competitive drop-in replacement** for the widely used Langevin dynamics MCMC method in the score-based diffusion model.  Our results have successfully validated this goal. We posit that delving into the relationship between denoising Gibbs sampling and Langevin dynamics MCMC will deepen our comprehension and refine sampling within the diffusion model. This exploration stands as an enticing avenue for future research.\n\n\n### Reference\n[1] Bao F, Li C, Zhu J, et al. Analytic-dpm: an analytic estimate of the optimal reverse variance in diffusion probabilistic models. ICLR, 2022.\n\n\n"}, {"review_id": "bYV6pj03Zv", "reviewer": "Reviewer_oYaT", "summary": "This work proposes a sampling method from implicitly defined clean model given a noisy model without retraining and demonstrates its applicability to high-dimensional data / images. ", "questions": "Please address the comments in the Weaknesses section. \n\n- What if log scale is used in Fig 3(a)? Then, it may show similar curves like Fig 3(b) unlike the claims in this work.\n- It is unclear why the comparable results of NCSNv2 with Langevin / the proposed method is ok. Are they using the same multi-level noise? Then, what is the advantage of the proposed method over the conventional Langevin dynamics?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- this work presents interesting theoretical results such as the uniqueness of underlying clean model given noisy data distribution, analytical covariance identity, and optimal Gaussian approx. \n- this work presents interesting corrections to prior works such as analytical DDPM, KL training Gibbs sampling.", "weaknesses": "- It is unclear if the proposed method (e.g., Eq. (12)) is effective in general. Some of the toy examples look promising, but some image examples may not look good (see Table 2 for large performance gaps with DDPM and NCSN++ / possibly more with recent methods). While this work emphasized that it will work well with a single level of fixed noise, eventually the proposed method with multi-level noise was used for more complex image cases without the issue of mode collapse. This inconsistency is my major concern.\n- Weak experimental results: some of the experiments may not effectively demonstrate the performance of the proposed scheme well. For example, it is unclear what this work wanted to show in Figures 5, 6 and 8. Comparison experiments may be needed.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We express our gratitude to the reviewer for recognizing the depth and value of our research. It's gratifying to see the acknowledgement of our novel theoretical findings. Here are our responses to the questions posed.\n\n### Effectiveness of the Proposed Method\nWe wish to emphasize that the central thrust of this paper is the elucidation of the optimal full-covariance Gaussian approximation within the denoising distribution. As outlined in Theorems 2.2 and 2.3, we've shown the theoretical superiority of our proposed approximation. Practically, our approach offers enhanced outcomes in single-noise level Gibbs sampling, achieving this with reduced parameter demands and training overhead, eliminating the need to learn the covariance.\n\nFor multi-level noise, our goal is to show that the proposed Gibbs sampler can be a competitive drop-in replacement to the widely used Langevin dynamics method in the score-based diffusion model. We also conduct an experiment on CelebA ($64{\\times}64{\\times}3$) dataset, which further shows that our method can generate high-quality images, see the PDF in the global rebuttal for the sample visualizations. Overall, these findings serve as a successful validation of our proposed proof of concept.\n\nIn reference to Table 2, our results are derived from a pre-trained NCSN v2, which exhibits performance on par with Langevin dynamics sampling. This comparison is fair (we use the same pre-trained network and noise schedule) and substantiates our experimental objective. Discrepancies in performance between other diffusion models can be attributed to variations in architectures and training procedures. The adaptation of these models in conjunction with our method represents a promising avenue for future research.\n\n\n### Inconsistency between applications\nThe paper's primary contribution lies in deriving the theoretical optimal Gaussian moment matching for Gaussian denoising distribution. We demonstrate that the proposed approximation is applicable in both single-noise and multi-level noise scenarios, treating them as separate applications. This ensures the proposed moment-matching technique is consistent across varied contexts. We appreciate the reviewer's feedback on this matter and will clarify further on this aspect in the revised manuscript.\n\n\n\n### Interpretation of the figures: \nWe are grateful to the reviewer for highlighting the need for a clearer objective for the figures presented. We will add the following clarification in our revised paper:\n\n- Figure 5 offers a visualization of the covariance estimation quality with different Rademacher samples, illustrating a discernible trade-off between estimation quality and computational cost.\n- Figure 6 portrays samples from three independent Gibbs Chains trained with single-level noise, adhering to standard plot conventions in MCMC research, and proving that high-quality images can be generated with only fixed single-level noise.\n- Figure 8 presents samples from the multi-level noise scenario, providing further insight into this specific case.\nIn the revised version of the paper, we will include more detailed clarifications for these figures, and we extend our gratitude to the reviewer for drawing attention to this need. \n\n### Figure comparison 4a vs 4b:\nFigure 4a depicts the diagonal covariance estimated using our method, while Figure 4b represents the diagonal covariance obtained through KL minimization. As shown in Theorem 2.3, the optimal diagonal covariance learned by KL (Figure 4b) should resemble that in Figure 4a, this is why they should look similar. However, the crucial distinction is that our method (Figure 4a). directly estimates the optimal covariance without the necessity for training. In contrast, the KL learning approach (Figure 4b) necessitates additional parameters and training, potentially leading to sub-optimal covariance estimation. This is further illustrated in Figures 4d and 4f, demonstrating that our method yields superior sample generation compared to the covariance learning approach. We'll provide a more detailed clarification in our revised paper.\n"}, {"review_id": "1yGf0aoD76", "reviewer": "Reviewer_DnqP", "summary": "This paper addressed the inconsistency problem in training energy-based models (EBMs) using denoising score matching (DSM). An efficient sampling framework, called (pseudo)-Gibbs sampling with moment matching, was proposed to enable effective sampling from the underlying clean model when given a well-trained 'noisy' model obtained through DSM. The presence of an underlying clean model within the noisy EBM was identified, and a pseudo-Gibbs sampling scheme incorporating analytical moment matching approximation was introduced for sampling from the clean model without requiring additional training. The effectiveness of the method was demonstrated in generating high-quality images for both single and multi-level noise settings. Additionally, the scalability of the method to high-dimensional datasets was illustrated. Comparisons with related methods highlighted the advantages of the proposed approach.", "questions": "In experiments, the authors used a noise-conditioned score network trained in [32], does the proposed model also excel in other settings?\n", "limitations": "None", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The theoretical derivation and presentation are clear and well-defined.", "weaknesses": "The limited number of experiments conducted to validate the proposed method. While the paper provides some demonstration of the effectiveness of the approach, further experiments across a wider range of scenarios and datasets would enhance the robustness and generalizability of the findings.\n\nThe authors should provide more details for those experiments in Supplementary Material.\n\nThe authors failed to demonstrate the wide applications of the proposed method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We're grateful to the reviewer for their valuable feedback, which has undeniably enhanced the quality of our work. The commendation on our theoretical contribution and presentation is heartening and much appreciated. Below, we respond to the queries presented.\n\n### Wider scenarios and applications \nWe want to thank the author for the feedback, emphasizing the need for broader experimentation to solidify the robustness and generalizability of our findings. In response to this, we have extended our methodology to include the CelebA dataset, notable for its heightened dimensionality (64x64x3), see the attached PDF in the global rebuttal for the sample results. The generated samples further highlight our method's capability to produce high-quality human face images. \n\nIn addition to the image generation setup, our method is generally applicable for any continuous Energy-based models or denoising score matching training, see [1] for a reference of different energy-based model applications. We will also include a discussion of wider applications in our revised paper.\n\n###  Experimental details\n The suggestion to provide more detailed information about the experiments in the Supplementary Material is well-received. We will include expanded details and the code in the revised paper. \n\n\n###  Applicability in other diffusion models\n Our proposed method boasts a versatile application spectrum, suitable for any diffusion model employing a forward Gaussian convolution process. This includes but is not limited to, the hierarchy latent variable model and the score-based model perspective. The profound interrelation between these perspectives is elaborated upon in [2]. We appreciate the reviewer's suggestion and will incorporate a detailed discussion on this topic in the updated version of our paper.\n\n### Reference\n[1] Luo C. Understanding diffusion models: A unified perspective. arXiv preprint arXiv:2208.11970, 2022."}, {"review_id": "JqBB5qVOOe", "reviewer": "Reviewer_afYY", "summary": "This paper proposed a sampling scheme called (pseudo)-Gibbs sampling with moment machine. It can be applicable to the sampling problem from  from the underlying clean model when a noisy model hat has been well-trained via DSM. The proposed sampling technique is used to sample images for  image generation and showed good efficiency on toy dataset.", "questions": "NIL", "limitations": "NIL", "rating": 5, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is well written\n2. The problem is a basic problem with many applications\n3.  The idea is new and seems working well.", "weaknesses": "The improvement in image generation seems minor.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer for acknowledging the strengths of our paper.\n\n### Contribution and Effectiveness\nWe wish to emphasize that this paper's central contribution is to elucidate the theoretical optimal full-covariance Gaussian approximation for the denoising distribution. As delineated in Theorems 2.2 and 2.3, our proposed approximation stands as theoretically optimal. In practical applications, our approach offers enhanced performance in single-noise level Gibbs sampling, all while reducing both parameter needs and training expenses—particularly, our method obviates the necessity of learning the covariance.\n\nIn the context of multi-level noise, we aim to position our Gibbs sampler as a viable alternative to the prevalent Langevin dynamics method within the score-based diffusion model. Our experiment on the CelebA dataset (64x64x3) bolsters this claim, illustrating our method's proficiency in generating high-quality images. For a glimpse of these visuals, kindly refer to the accompanying PDF in the comprehensive review. Such findings effectively echo our foundational aspirations. We are inclined to believe that a deeper dive into the interplay between denoising Gibbs sampling and Langevin dynamics MCMC could shed additional light, further refining the diffusion model's sampling process. This presents a promising avenue for future exploration."}, {"review_id": "Pj0t9SoDFX", "reviewer": "Reviewer_1SZs", "summary": "This paper proposes a new denoising scoring matching (DSM) technique to train energy-based models. They address the well-known inconsistency problem of DSM objectives by training a \"noisy\" model and then sampling \"clean\" samples from the \"noisy\" model by learning a denoising distribution.\n\nThe contributions of the paper are as follows:\n\n1) Theoretically, they show that a \"clean\" model is uniquely identified by an EBM learning a noisy data distribution.\n2) A sampling strategy inspired from Gibbs sampling to approximate the denoising distribution, which allows to sample from the underlying clean model without the need to retrain the noisy model.\n3) Practically, they show how to scale their approach to high-dimensional data by overcoming the quadratic cost of evaluating the Hessian of the DSM for every $\\tilde{x}$. ", "questions": "As stated above in the weaknesses section, my main doubt revolves around how does the Hessian approximation affect the training stability and overall output samples quality.", "limitations": "The authors have partly discussed the limitations in section 3, but an experiment comparing the quality of the approximation of the generated samples with the full Hessian and the approximated Hessian is welcome to strengthen the demonstration.\n\nNo potential negative societal impact at this stage.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Well-written paper, easy to follow and most claims are justified either by citing previous work or in appendix / math proofs\n- Prior work is duly cited\n- Experiments on synthetic and real-world (MNIST/CIFAR) datasets are appreciated to assess the quality of the technique proposed", "weaknesses": "I'm not an expert in this field and I have a hard time evaluating the overall contribution of this paper to the field of \"Generative Models\". In particular, it feels like this denoising techniques introduce a scalability problem (although partly discussed and solved by contribution 3) in order to solve the inconsistency problem of score-based matching training. A discussion on the quality of the approximation and comparing the full Hessian with its approximation could help answer my doubts.\n\nBesides, I think that describing the full sampling algorithm could help the reader get a better picture of the whole training and sampling procedure. This would clarify what has been made in the experiments section.\n\nI list below a few comments / typos I have found while reviewing the paper:\n\n1) l.33: there is a missing $\\nabla_x$ in the definition of the score matching function. Otherwise, I don't see how $\\nabla_x s_{p_{\\theta}}(x) = -\\nabla^2_x f_{\\theta}(x)$.\n2) l.26: \"it is equivalent to MLE\". A reference is welcome. \n3) l.34: \"under certain regularity conditions, the Fisher divergence is equivalent to SM\". Reference?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We're grateful for the constructive feedback.  We're glad the clarity, references, and experiments in our paper resonated with you. Your feedback greatly enhances the paper's quality.\n\n### Typos and Reference\nWe want to thank the reviewer for pointing out the typo $s(x)\\equiv \\nabla_x\\log p(x)$. The reference for the equivalent of MLE and KL training can be found in [1] and the connection between the Fisher divergence and the Score matching can be found in [2], which we will include in the revised version.\n\n### Scalability\nTo demonstrate our method's scalability, we further conduct an experiment using the CelebA dataset, characterized by its dimensions of $64{\\times}64{\\times}3$. For visual evidence of the generated samples, kindly refer to the PDF in the global rebuttal. The results highlight our method's capacity to produce high-resolution, quality human faces.\n\n### Stability and Approximation Quality\nFurthermore, we would like to provide a clear explanation regarding our methodology: the estimation of the Hessian is required only at the time of generation, not during training. The training process adheres to the traditional framework of denoising score-matching. As visualized in Figure 5, our depiction of covariance estimation with different Rademacher samples will be supplemented in the revised version with a study on how varying samples might influence the quality of the results. We deeply appreciate this insightful suggestion, and it will certainly enhance the comprehensiveness of our research.\n\n### Reference\n[1]: Bishop C M, Nasrabadi N M. Pattern recognition and machine learning. New York: springer, 2006.\n\n[2]: Hyvärinen A, Dayan P. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 2005.\n"}]

Author rebuttal:
We want to thank all the reviewers again for the time spent and their valuable feedback that improves the quality of our work. Their recognition of our theoretical contribution, which forms the cornerstone of our paper, is deeply encouraging. Here's a succinct overview of our primary contributions.

### Contributions
1. We introduced an analytical data-dependent full covariance moment-matching objective for the posterior within the Denoising Score Matching framework and validated its **theoretical optimality**.
2. We pinpointed the **existence of the clean model** in a noisy EBM trained via denoising score matching and developed a moment-matching Gibbs sampler for the energy-based models. This approach drastically cuts down on parameters, training expenses, and boosts the performance of the conventional denoising Gibbs sampling [1,2] with amortized covariance learning [3].
3. We further scale moment-matching estimation for high-dimensional data, which provides a **trade-off between flexibility and efficiency**.  We also demonstrate the scalable method can:
     - generate high-quality natural images with even **one fixed-level noise**.
     - achieve **competitive performance** compared to the widely used Langevin Dynamics sampling in multiple-noise settings.

### Effectiveness
A shared concern among the reviewers centres on the effectiveness of the proposed method. We want to highlight that, as mentioned in contribution point 2, we've demonstrated that our method not only yields enhanced performance in single-noise level settings compared to traditional denoising Gibbs sampling but also achieves this with reduced training costs, notably eliminating the need to learn the covariance.

For multi-level noise, our goal is to provide a proof of concept that the proposed Gibbs sampler can serve as a competitive drop-in replacement for the widely used Langevin dynamics method in the score-based diffusion model. Our preliminary results have successfully validated this goal for CIFAR10. Since submission, we have run further experiments on the CelebA dataset, demonstrating scaling to image sizes $64{\times}64{\times}3$, see the attached PDF for representative samples. We can see our method is also able to generate high-quality images in a higher dimension space. Experimental details and numerical evaluations will be added in the camera-ready version.

### Reference
[1] Bengio Y, Yao L, Alain G, et al. Generalized denoising auto-encoders as generative models. Advances in neural information processing systems, 2013, 26.

[2] Alain G, Bengio Y. What regularized auto-encoders learn from the data-generating distribution. The Journal of Machine Learning Research, 2014

[3] Kingma D P, Welling M. Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114, 2013.


Meta-review: This paper proposes a new sampling framework, i.e., (pseudo)-Gibbs sampling with moment matching, for the sake of drawing samples from the implicitly-defined clean model given the noisy energy-based model, which is well-trained via denoising score matching method. The paper provides both theoretical and empirical justification of the proposed method.  In response to the concerns initially raised by the reviewers in their initial assessments, the rebuttal systematically addressed these issues. Following extensive discussions, a consensus was achieved among all five reviewers, culminating in the decision to accept the paper. This agreement is attributed to the paper's profound theoretical justification, and sufficient empirical validation that it offers for this particular study. The AC agrees with the reviewers and recommends accepting the paper. To further improve the paper quality, AC suggests the authors to revise their paper by taking into account all the suggestions provided by the reviewers, including integrating extra experimental results conducted during the rebuttal phase.
Decision: Accept (poster)