Multimeasurement Generative Models
Authors: Saeed Saremi
Rupesh Kumar Srivastava
Venue: ICLR 2022 Poster (2022)
Keywords: energy based models
Langevin MCMC
score matching
denoising autoencoders
empirical Bayes

Abstract: We formally map the problem of sampling from an unknown distribution with a density in $\mathbb{R}^d$ to the problem of learning and sampling a smoother density in $\mathbb{R}^{Md}$ obtained by convolution with a fixed factorial kernel: the new density is referred to as M-density and the kernel as multimeasurement noise model (MNM). The M-density in $\mathbb{R}^{Md}$ is smoother than the original density in $\mathbb{R}^d$, easier to learn and sample from, yet for large $M$ the two problems are mathematically equivalent since clean data can be estimated exactly given a multimeasurement noisy observation using the Bayes estimator. To formulate the problem, we derive the Bayes estimator for Poisson and Gaussian MNMs in closed form in terms of the unnormalized M-density. This leads to a simple least-squares objective for learning parametric energy and score functions. We present various parametrization schemes of interest including one in which studying Gaussian M-densities directly leads to multidenoising autoencoders—this is the first theoretical connection made between denoising autoencoders and empirical Bayes in the literature.  Samples in $\mathbb{R}^d$ are obtained by walk-jump sampling (Saremi & Hyvarinen, 2019) via underdamped Langevin MCMC (walk) to sample from M-density and the multimeasurement Bayes estimation (jump). We study permutation invariant Gaussian M-densities on MNIST, CIFAR-10, and FFHQ-256 datasets, and demonstrate the effectiveness of this framework for realizing fast-mixing stable Markov chains in high dimensions.

Published as a conference paper at ICLR 2022
MULTIMEASUREMEN GENERATIVE MODELS
Saeed Saremil.2& Rupesh Kumar Srivastava! 'NNAISENSE Inc.
Redwood Center, UC Berkeley
(saeed, rupesh]@nnat wwwaeisseseecc.mmm com
ABSTRACT
We formally map the problem of sampling from an unknown distribution with density in IR to the problem of learning and sampling a smoother density in IR Md obtained by convolution with a fixed factorial kernel: the new density is referred to M-density and the kernel multimeasurrment noise model (MNM). The M-density in Idd is smoother than the original density in R easier to learn and sample from, yet for large M the two problems mathematically equivalent since clean data can be estimated exactly given a multimeasurement noisy observation using the Bayes estimator. To formulate the problem, we derive the Bayes estimator for Poisson and Gaussian MNMs in closed form in terms of the unnormalized M-density. This leads a simple least-squares objective for learning parametric energy and score functions. We present various parametrization schemes interest including one in which studying Gaussian M-densities directly leads to multidenoising autoencoders-tthie is the first theoretical connection made between denoising autoencoders and empirical Bayes in the literature. Samples in Rd are obtained by julk--jump sampling (Saremi & Hyvärinen, 2019) via underdamped Langevin MCMC (walk) to sample from M-density and the multimeasurement Bayes estimation (jump). We study permutation invariant Gaussian M-densities MNIST, CIFAR-10, 10, and FFHQ-256 datasets, and demonstrate the effectiveness of this framework for realizing fast-mixing stable Markov chains in high dimensions.
INTRODUCTION
Consider a collection of i.i.d. samples {a:) assumed have been drawn from an unknown distribution with density px Rd. An important problem in probabilistic modeling is the task of drawing independent samples from px. which has numerous potential applications. This problem is typically approached in two phases: approximating Px, and drawing samples from the approximated density. In unnormalized models the first phase is approached by learning the energy function fx associated with the Gibbs distribution Px exp(-fx fx ). and for the second phase one must resort to Markov chain Monte Carlo methods, such as Langevin MCMC, which are typically very slow to mix in high dimensions. MCMC sampling is considered an "art" and we do not have black box samplers that converge fast and are stable for complex (natural) distributions. The source of the problem is mainly attributed to the fact that the energy functions of interest are typically highly nonconvex.
A broad sketch of our solution to this problem is to model smoother density in an M-fold expanded space. The new density denoted by p(y), called M-density, is defined in R Ma where the boldfaced y is shorthand for (V1...., . yMs) M-density is smoother in the sense that its marginals Pm (ym) are obtained by the convolution Pmu (ym) = -P (ym (xr (p(()) (r) dd with smoothing kernel pm (ym (x) which for most of the paper we take to be the isotropic Gaussian:
Y,. = X + N(0,, o Ta).
Although we bypass learning p(z), the new formalism allows for generating samples from p(x) since X can be estimated exactly given Y y (for large M). To give a physical picture, approach here is based on "taking apart" the complex manifold where the random variable X is concentrated in Rdd and mapping
smoother manifold in IR Md where Y (Yi.... Ym) is now concentrated.
Smoothing a density with a kernel is a technique in nonparametric density estimation that goes back to Parzen (1962). In kernel density estimation, the estimator of p(x) is obtained by convolving the

---

Published as a conference paper at ICLR 2022
empirical measure with kernel. In that methodology, the kernel bandwidth (o, for Gaussian kemels)
is adjusted to estimate p(æ) in IR given a collection of independent samples {2x }=.) This estimator, like most nonparametric estimators, suffers from a severe curse of dimensionality (Wainwright, 2019).
But what if the kernel bandwidth is fixed: how much easier is the problem of estimating p(y)? This question is answered in (Goldfeld et al., 2020), where they obtained the rate of convergence eo 1/2 (measured using various distances) in remarkable contrast to the well-known n rate for estimating p(æ). This nonparametric estimation result is not directly relevant here, but formalizes the intuition that learning p(y) J p(y|=)p(()))) is lot simpler than learning p(xr).
With this motivation, we start with an introduction to the problem of learning unnormalized p(y), based on independent samples from p(a)). This problem was formulated by Vincent (201 1) using score matching (Hyvärinen, 2005). It was approached recently with the more fundamental methodology of empirical Bayes (Saremi & Hyvärinen, 2019). The idea to use the Bayes estimator of X given Y the study of which is at the root of the empirical Bayes approach to statistics (Robbins, 1956), in least-squares objective. This machinery builds on the fact that the estimator F(y) = E| = yl can be expressed in closed form in terms of unnormalized p(y) (Sec. 3.1). For Gaussian kernels, the learning objective arrived at in (Saremi & Hyvärinen, 2019) identical (up to a multiplicative constant) to the denoising score matching formulation (Vincent, 2011) 1), but with new insights rooted in empirical Bayes which is the statistical framework for denoising.
The main problem with the empirical Bayes methodology is that p(xly) remains unknown and cannot be sampled from. The estimator ((y) = E| = y) can be computed, but the concentration of the posterior p(rly) around the mean is not in our control. Our solution to this problem starts with an observation that is very intuitive from a Bayesian perspective: one can sharpen the posterior by simply taking more independent noisy measurements. This scheme is formalized by replacing p(y p(y|æ) with the factorial kernel p(y
M
plylx)
II
Pm (ym (z), where y = yar):
m=1
which we name multimeasurement noise model (MNM). Now, the object of interest is different density which we call M-density obtained by convolving p(z) with the factorial kernel:
p(y) = p(yl=)p(z))))
(2)
This formally maps the original problem of drawing samples from p(æ: to drawing samples from p(y) for any fixed noise level since the estimator of X given Y = y is asymptotically exact. We quantify this for Gaussian MNMs using the plug-in estimator (the empirical mean of measurements).
Smooth & Symmetric! Consider Gaussian MNMs with equal noise level a in the regime of large o. large M such that a /d/M is "small". 1 In that regime, the complex manifold associated with the data distribution is mapped to a very smooth symmetric manifold in a much higher dimensional space. The original manifold can be reconstructed via a single step by computing (y)) Due to equal noise levels, the manifold associated with M-density is symmetric under the permutation group:
p(y1..... Mar) p(Na(1), M*(a)));
where is a permutation of indices (Fig. 1). Although we develop general methodology for studying M-densities, in the later part of the paper we focus on permutation invariant Gaussian M-densities.
The paper is organized as follows. In Sec. 2, we derive Bayes estimators for Poisson and Gaussian MNMs. In Sec. 3, we present the least-squares objective for learning Gaussian M-densities. We also give weaker formulation of the learning objective based score matching. Sec. 4 is devoted the important topic of parametrization, where we introduce multidenoising autoencoder (MDAE) in which we formally connect M-densities to the DAE literature. DAEs have never been studied for factorial kernels and the emergence of MDAE as generative model should be of wide interest. In addition, we introduce metaencoder formulated in an unnormalized latent variable model, which is mainly left as a side contribution. In Sec. 5, we present the sampling algorithm used in the paper. In Sec. 6, we present our experiments on MNIST, CIFAR-10 10, and FFHQ-256 datasets which were focused on permutation invariant M-densities. The experiments demonstrating the effectiveness of this method in generating fast mixing and very long Markov chains in high dimensions. Related works are discussed in Sec. 7, and we finish with concluding remarks.
mainly of qualitative nature
The regime a dd/M & is obtained in our analysis of the highly suhoptimal plug-in estimator (Sec. 2.3).
2

---

Published as a conference paper at ICLR 2022
0.25
0.20
0.10
0.05
-6
-4
-2
0
2
(a)
p(2)
(b) p(y)
(c) logp(y)
(d) V logp(y)
Figure 1: (M-density) (a) A mixture of Gaussian in 1d. (b,c,d) The M-density (M 2. 02), the corresponding log-density and score function are visualized (based on calculations in Appendix A).
Notation. The subscripts dropped from densities and energy functions when it is clear from their arguments: p(y) py (y), ply x)) - PY)X== f(y) (y), , etc. Bold fonts are reserved for multimeasurement random variables:
= (Yi....... Yar). The following are shorthand notations: [M] = (I...., .............. M) and V = V ym. Throughout, V is the gradient with respect to inputs (in RMd), not parameters. The following convention is used regarding parametric functions: fo(.) = f(;0). Different parametrization schemes come with difffernnt set of parameters, the collection of which we denote by For all the datasets used the paper, X takes values the hypercube [0, 1]d.
FORMALISM: MULTIMEASUREMENT BAYES ESTIMATORS
This work is based generalizing the empirical Bayes methodology MNMs. It is well known that the least-squares estimator of given Y (for any noise model) is the Bayes estimator:
f ())(()
(y)) = f p(y(x)p(z)dz
(4)
Next we study this estimator la Robbins (1956) for Poisson (the Poisson kernel was the first example studied in 1956) and Gaussian MNMs. In both cases the estimator F(y) is derived to be functional of the joint density p.y). In addition, 7(y) is invariant to scaling p(y) ) by a constant, therefore one ignore the partitioon function in this estimation problem. This is the main appeal of this formalism. Analytical results for Poisson MNMs are included to demonstrate the generality the new formalism, but will not pursue challenges regarding sampling discrete distributions high dimensions (see Remark 7).
generative model in our experiments for technical reasons due the
2.1
POISSON MNM
Let X be random variable taking values in R.. The Poisson MNM is defined by:
2: yi! 31
p(ylx) c Mæ
E.
The numerator in r.h.s. of Eq. is computed next. The measurement index below is an arbitrary index [M] used for absorbing x such that xp(y|z) has the same functional form as p(y[æ))
f xppyys))(x)de f
- MM m + +1) (ym + 1)!
x ((0m+1) I)
II
2 p(x)dd = (ym + 1)p(y + 1m),
yn!
where 1m is defined as vector whose component is 8mm. Using Eq. 4, it immediately follows
Fly) = (ym + ply + 1mm
E [M].
(5)
p(y)
We emphasize that the dependency F(y) on the noise channel (measurement index) m that appears on the right hand side of the expression above is only artifact of the calculation (we observe this again for Gaussian MNMs). The result above holds for any measurement index m E [M], therefore
Vu 1)p(y + 1,m) (Ym + 1)p(y 1 mm) for all m' [M].
3

---

Published as a conference paper at ICLR 2022
Example.
We can derive the estimator (() analytically for p(x) We first derive p(y):
(2) vt)! II wn! (M + -1-1
p(y)
where the sums/products are over the measurement indices E [M]. Using Eq. 5, follows
F(y) (ym + 1) ply + 1m)
(ym 1) E, U1 +
Σ, D2 +1
p(y)
ym M +1
M +1
As expected one arrives at the same result by computing Eq. for any measurement index m.
2.2 GAUSSIAN MNM
Let X be random variable in Rd. The Gaussian MNM is defined by:
M
p(y|))=
II Pm (ym [x), where Pm (ym (x) = N(ym; 2,, o Ta).
(6)
mmm
follows (as in the Poisson MNM, m the equation below is an arbitrary index in [M]):
om V mp(y(z:) (x -mm)pyy (z)).
We multiply both sides of the above expression by p(z) and integrate over æ. The derivative Vm (with respect to ym) and the integration over commute, and using Eq. 4 it follows
o mpy) 2(y ((y) - ymp(y);
which we simplify by dividing both sides by p(y):
Fly) = Dm + o m. log ply). m E [M].
(7)
This expression is the generalization of the known result due to Miyasawa (1961), As in the Poisson MNM, the result above holds for any m E [M], therefore:
Vm 2 m. log p(y) Ym + V. m logp(y) for all m, m' [M].
(8)
Example. We also studied the M-density for Gaussian MNMs analytically. The calculations are insightful and give more intuitions on the new formalism. We refer Appendix A for the results.
Remark (#g (m)) (y) notation). For parametric M- M-densities, Eq. 8 in general an approximation. We use the superscript m to emphasize that there are M ways to compute the parametric estimator:
a/m² (y) Ym + o V log poly)
Remark (a M notation). Gaussian MNMs with the constant noise level a are denoted by a & M. For a & M models, the M-density ply) (resp. the score function V log p(y)) is invariant (resp. equivariant) with respect to the permutation group π. [M] + [M]. See Fig. 1 for an illustration.
2.3 CONCENTRATION OF THE PLUG-IN ESTIMATOR
The plug-in estimator of X is the empirical mean of the multiple noisy measurements we denote by Tmean (y) = M-1 Σn Vm. This estimator highly suboptimal but analysis high dimensions for Gaussian MNMs is useful. Due to the concentration of measure phenomenon (Tao, 2012) we have
|| meaan (v) || Ve,
(9)
where den ("eff" for effective) is given by
1/2
1 Ceff M
Σ
(10)
The calculation is straightforward since Vm x (y) = -M Em Em which has the same law as N(0, 2. This calculation shows that the estimator of X in Rd concentrates at the true value at worst-case rate O( O(vd/M) (consider replacing the sum in Eq. 10 by Mo This analysis is very conservative, as it ignores the correlations between the components y in Md (within and across noise channels), and one expects (much) tighter concentration for the optimal estimator. In the next section, we present an algorithm to learn the multimeasurement Bayes estimator based on independent samples from p(x).
+ Ems where Em N La). It follows:
Thhe derivation is straightforward using So o o da =a- -1-A for 0,8 B
4

---

Published as a conference paper at ICLR 2022
3 LEARNING GAUSSIAN M-DENSITIES
3.1 NEURAL EMPIRICAL. BAYES
In this section, we focus on learning Gaussian M-densities using the empirical Bayes formalism. We closely follow the approach taken by Saremi & Hyvärinen (2019) and extend it to M-densities. The power of empirical Bayes lies in the fact that it is formulated in the absence of any clean data. This reflected by the fact that F(y) is expressed in terms of p(y) which can in principle be estimated without observing samples from p(æ) (Robbins, 1956); we generalized that factorial kernels in Sec. 2. What if we start with independent samples [ from p(z) and our goal is to learn p(y)?
A key insight in neural empirical Bayes was that the empirical Bayes machinery can be turned on its head in the form of a Gedankenexperiment (Saremi & Hyvärinen, 2019): we can draw samples (indexed by j) from the factorial kernel Yiz ~ p(y|0+)) and feed the noisy data to the empirical Bayes "experimenter" (the word used in 1956). The experimenter's task (our task!) is to estimate X, but since we observe X It, the squared Eg norm ||=; (y()))) ||2 serves as a signal to learn p(y). and also (y) for unseen noisy data (as a reminder the Bayes estimator is the least-squares estimator).
Next, present the least-squares objective to learn p(y) x e for Gaussian M-densities. It is important to note that F(y) is expressed in terms of unnormalized p(y), without which we must estimate the partition function (or its gradient) during learning. Here, we can choose any expressive family of functions to parametrize F(y) which is key to the success of this framework. Using our formalism (Sec. 2.2), the Bayes estimator takes the following parametric form (see Remark 1):
2((²) (y) = 3m - o foly,, m E [M].
(11)
There are therefore M least-squares learning objectives
((m) (0) E (r y ((m) (z. 0), where (zm) y: 0) ||+ 3/(m) (y) 12
(z.
(12)
that in principle need to be minimized simultaneously, since as a corollary of Eq. 8 we have:
L (m) 8*) = L ((mm)( y; 8*) for all m, m' [M], z Rd, y Rddd
(13)
The balance between the M losses can be enforced during learning by using a softmax-weighted sum of them in the learning objective, effectively weighing the higher losses more in each update. However, in our parametrization schemes, coming up, simply taking the mean of the M losses as the learning objective proved to be sufficient for a balanced learning across all the noise channels:
M
(())
Σ
(14)
M
m=1
The above learning objective the one we use in the remainder of the paper.
3.2 DENOISING SCORE MATCHING
One can also study M-densities using score matching with the following objective (Hyvärinen, 2005):
a(0) Ey fyyeyy Vy p(y))
(15)
In Appendix B. we show that the score matching learning objective is equal (up to an additive constant independent of 0) to the following multimeasurement denoising score matching (MDSM) objective:
M ((0) Σ (0), where (m) (z,y) V Vmff(y))+ + ym -
12
(16)
This is simple extension of the result by Vincent (2011) 1) to Gaussian MNMs. The MDSM objective and the empirical Bayes' (Eq. 14) are identical (up to multiplicative constant) for @ M models.
Remark 3. Compared to neural empirical Bayes (NEB), denoising score matching (DSM) takes very different approach regarding learning M-densities. DSM starts with score matching (Eq. 15). NEB starts with deriving the Bayes estimator of X given Y y for a known kernel p(y (z) (Eq. 7). NEB is a stronger formulation in the sense that two goals are achieved at once: learning M-density and learning Zyy). What remains unknown DSM is the latter, and knowing the estimator is key here. Without it, we cannot draw formal equivalence between px and pv (see Sec. 1). We return to this discussion at the end of the paper from a different angle with regards to denoising autoencoders.
5

---

Published as a conference paper at ICLR 2022
4 PARAMETRIZATION SCHEMES
We present three parametrization schemes for modeling Gaussian M-densities. Due to our interest in a & M models we switch to lighter notation. Here, the learning objective (Eq. 14) takes the form
1 M
[())
( (e @ M y + o VoVSs
(17)
where x @ M denotes (z....,,2) .)) repeated M times. Parametrization schemes fall under two general groups: multimeasurement energy model (MEM) and mulrimeasurement score model (MSM). In what follows, MDAE is an instance of MSM, MEM & MUVB are (closely related) instances of MEM.
4.1 MDAE (MULTIDENOISING AUTOENCODER)
The rigorous approach to learning Gaussian M-densities is to parametrize the energy function In that parametrization. V is computed with automatic differentiation and used in the objective (Eq. 17). That is computational burden, but it comes with a major advantage as the learned score function is guaranteed to be gradient field. The direct parametrization of V f is problematic due to this requirement analyzed by (Saremi, 2019); see also (Salimans & Ho, 2021) for a recent discussion. Putting that debate aside, in MDAE we parametrize the score function explicitly go : Radd - Riddd Then, go replaces- Lo in Eq. 17. In particular, we consider the following reparametrization
go(y) = (ve(y) - -y)/o²
(18)
simply motivated by the fact would like cancel y in Eq. 17, otherwise the loss starts at very high values at the initialization. This is especially so in the regime of large o, M, and d. It follows:
( M ((r.))-p(yle))(e)) x @ M - ue (1)
(19)
This is very intriguing and it is worth taking moment to examine the result: modeling M-densities is now formally mapped to denoising multimeasurement noisy data explicitly. It is a DAE loss with a multimeasurement twist. Crucially, due our empirical Bayes formulation, the MDAE output (y) becomes a parametrization of the Bayes estimator(s) (combine Eq. 18 and Eq. 11):
(1m) (y: 0) = I'm (y; 0).
(20)
This makes strong theoretical connection between our generalization of empirical Bayes to factorial kernels and new learning paradigm under multidenoising autoencoders, valid for any noise level
4.2 2EM (MEM with QUADRATIC form with an optional METAENCODER)
Can we write down an energy function associated with the score function in Eq. 18? The answer is no, since go is not gradient field in general, but we can try the following (1. ()):
sy = 20 |yy - v, h (,
(21)
Ignoring he for now, by calculating -o we do get both terms in Eq. 18, plus other terms. The function hy which we call metaencoder is optional here. Intuitively, it captures the "higher order interactions" between and v, beyond the quadratic term (see below for another motivation). The metaencoder is implicitly parametrized by (via Myy), while having its own set of parameters Ç.
4.3 MUVB (MULTIMEASUREMENT UNNORMALIZED VARIATIONAL BAYES)
The expression above (Eq. 21) is simplification of an unnormalized latent variable model that one can set up, where we take the variational free energy to be the energy function. This builds on recent studies towards bringing together empirical Bayes and variational Bayes (Saremi, 2020a;b). We outline the big picture here and refer to Appendix C for details. Latent variable models operate on the principle of maximum likelihood, where one is obliged to have normalized models. Since our model is unnormalized we consider setting up latent variable model with unnormalized conditional density. Essentially both terms in Eq. 21 arise by considering (z the vector of latent variables)
(ylcc)(((2) xx exp
ly - un he he(y,
(22)
20
named metalikelihood which further underscores the fact that it is unnormalized. The full expression for the energy function also involves the posterior (2ly)) As remark, note that what is shared between all three parametrization schemes is v, although vastly different in how F(y) is computed.

---

Published as a conference paper at ICLR 2022
5 SAMPLING ALGORITHM
Our sampling algorithm is an adaptation of walk-jump sampling (WIS) (Saremi & Hyvärinen, 2019). We run MCMC algorithm sample M-density by generating a Markov chain of multimeasurement noisy samples. This is the walk phase of WJS schematized by the dashed arrows in Fig. 2. At arbitrary discrete time clean samples are generated by simply computing F(yx) ) (# is dropped for a clean notation). This is the jump phase of WJS schematized by the solid arrow in Fig. 2. What is appealing about multimeasurement generative models the fact that for large M. p(æly=) is highly concentrated around its mean Ty))), therefore this scheme is exact sampling scheme-thi was in fact our original motivation to study M-densities. For sampling the M-density (the walk phase), we consider Langevin MCMC algorithms that are based on discretizing the underdamped Langevin diffusion:
dv, -V, 2 V fyy, dtt + (V2yu)dB.
(23)
dy: = v,dt.
Here y is the friction, u. the inverse mass, and B, the standard Brownian motion in Raid, Discretizing the Langevin diffusion and their analysis are challenging problems due to the non-smooth nature of the Brownian motion (Mörters & Peres, 2010). There has been a significant progress being made however in devising and analyzing Langevin MCMC algorithms,, e.g. Cheng et al. (2018) introduced an algorithm with a mixing time that scales as O(vd) in notable contrast to the best known (d) for MCMC algorithms based on overdamped Langevin diffusion. The dimension dependence of the mixing time is of greal interest here since expand the dimension M-fold. To give an idea regarding the dimension, for the 4 & model on FFHQ-256, Md 105. implemented (Cheng et al., 2018, Algorithm 1) in this paper which to our knowledge is its first use for generative modeling. In addition, we used Langevin MCMC algorithm due to Sachs et al. (2017), also used by Arbel et al. (2020). Note, in addition to y and Langevin MCMC requires ö, the step size used for discretizing Eq. 23.
Figure 2: (WJS schematic) In this schematic, the Langevin walk denoted by the dashed arrow. The jump is denoted by the solid arrow which is deterministic. The jumps in WJS are asynchronous (Remark 5). In presenting long chains in the paper we show jumps at various frequencies denoted by Ak (Remark 6). We use the same MCMC parameters for all noise channels due to permutation symmetry in M models. See Appendix D for more details.
Yk-1
Y&
Yk41
a(yc)
(a) & 4, MDAE
(b) 4 8, MDAE
Figure 3: (WJS chains on FFHQ-256) The chains are shown skipping Ak = 10 steps (no warmup).
We used Algorithm 1 with 8 = 2,y 1/2, 1. Transitions are best viewed electronically.
7

---

Published as a conference paper at ICLR 2022
Figure 4: (4 & gallery) Samples from our FFHQ-256, MDAE, & 8 model.
EXPERIMENTAL RESULTS
For all models with the MDAE and MEM parametrizations, we used U2-Net (Qin et al., 2020) network architecture, a recent variant of UNet (Ronneberger et al., 2015), with a few simple changes (see Appendix E). For the MUVB parametrization for MNIST (evaluated in Sec. G.3), we used Residual networks for the encoder, decoder and metaencoder. Despite the inherent complexity of the task of learning high-dimensional energy models, it is notable that our framework results a single least-squares objective function, which we optimized using Adam (Kingma & Ba, 2014) without additional learning rate schedules. Additional details of the experimental setup are Appendix E.
Arguably, the ultimate test for an energy model is whether one can generate realistic samples from it with single fast mixing MCMC chain that explores all modes of the distribution indefinitely, starting from random noise. We informally refer to them as lifelong Markov chains. To give a physical picture, a gas of molecules that has come in contact with a thermal reservoir does not stop mid air after thermalizing- arriving at its "first sample" distribution as long as the physical system exists. To meet this challenge, the energy landscape must not have any pathologies that cause the sampling chain to break or get stuck in certain modes. In addition, we need fast mixing (Langevin) MCMC algorithms, very active area of research by itself.
continues generating samples from the Boltzmann
To put this goal in context, in recent energy models for high-dimensional data (Xie et al., 2016; 2018; Nijkamp et al., 2019; Du & Mordatch, 2019; Zhao et al., 2020; Du et al., 2020; Xie et al., 2021), sampling using MCMC quickly breaks or collapses to a mode and chains longer than a few hundred steps were not reported. Thus, evaluation in prior work relies on samples from independent MCMC chains, in addition by using heuristics like "replay buffer" (Du & Mordatch, 2019). In this work, we report FID scores obtained by single MCMC chains, the first result of its kind, which we consider as a benchmark for future works on long run MCMC chains (see Table 4 for numerical comparisons).
For MNIST, we obtain fast-mixing chains for over 1M steps using MDAE and MUVB. On CIFAR- 10 and FFHQ-256, we obtain stable chains up to 1 M and 10K steps, respectively, using MDAE The results are remarkable since energy models-that learn single energy/score function-have never been scaled to 256 x 256 resolution images. The closest results are by Nijkamp et al. (2019) and Du et al. (2020) to (Nijkamp et al., 2019, Figs. 1, 2). For CIFAR-10, we report the FID score of 43.95 for 1 & 8 model, which we note is achieved by a single MCMC chain (Fig. 13); the closest result FID score obtained for long run MCMC is 78. /2 by Nijkamp et al. (2022) which is not on single chains, but on several parallel chains (in that paper the "long run" chains are K steps vs. M steps in this work).
CelebA (128x 128) in particular, results in Figs. 3, 4 can be compared
Our detailed experimental results are organized in appendices as follows. Appendix F is devoted to an introduction to Langevin MCMC with demonstrations using our FFHQ-256 4 & 8 model. Appendix G is devoted to MNIST experiments. At first we compare 1/4 @ 1 and 1 @ 16 (Sec. G.1). These models are statistically identical regarding the plug-in estimators (Eq. 10) and they arrive at similar training losses, but very different as generative models. This sheds light on the question why we considered such high noise levels in designing experiments (Fig. 5). We then present the effect of increasing M for fixed a, the issue of time complexity, and mixing time vs. image quality trade-off (Sec. G.2). The discussions in Sec. G.I and Sec. G.2 are closely related. We then compare MDAE and MUVB for @ 4 with the message that MDAE generates sharper samples while MUVB has better mixing properties (Sec. G.3). We then present one example
lifelong Markov chain (Sec. G.4). Our MDAE model struggles the 10 challenge due to optimization problems. The results and further discussion are presented in Appendix H. Finally, in Appendix we show chains obtained for FFHQ-256, including the ones that some of the images in Fig. 4 were taken from.
8

---

Published as a conference paper at ICLR 2022
Figure 5: (single-step multidenoising) A clean data (æ), a sample from MNM (y) for the @ 4 model (Fig. 3a), the outputs of MDAE (V), and meann (y) are visualized from left to right. The outputs (y)) correspond to the Bayes estimator(s) which as expected from theory (Eq. 8) are indistinguishable.
RELATED WORK
Despite being fundamental framework for denoising, empirical Bayes has been surprisingly absent in the DAE/DSM literature. Raphan & Simoncelli (201 I) first mentioned the connection between empirical Bayes and score matching, but this work is practically unknown in the literature; for much of its development, DSM was all tied to DAEs (Vincent, 201 I; Alain & Bengio, 2014); see (Bengio et al., 2013) for an extensive survey. Abstracting DSM away from DAE is due to Saremi et al. (2018) where they directly parametrized the energy function. In this work, we closed the circle, from empirical Bayes to DAEs. One can indeed use the MDAE objective and learn a generative model: highly generalized DAE naturally arise from empirical Bayes. The important subtlety here is, as we remarked in Sec. 3.2, we can only arrive at this connection starting from our generalization of empirical Bayes, not the other way around. Of course, this core multimeasurement aspect of our approach, without which we do not have a generative model, does not exist in the DAE literature."
To address the problem of choosing noise level in DSM (Saremi al., 2018), Song & Ermon (2019) studied it with multiple noise levels by summing up the losses using a weighing scheme. See (Li et al., 2019; Chen et al., 2020; Song & Ermon, 2020; Kadkhodaie & Simoncelli, 2020; Salicoer-Marrtneen et al., 2020) in that direction. The learning objectives in these models are based on heuristics in how different noise levels are weighted. Denoising diffusion models (Ho et al.. 2020; Song al., 2020; Gao et al., 2020) follow the same philosophy while being theoretically sound. Sampling in these models are based on annealing or reversing a diffusion process. Our philosophhy here is fundamentally different. Even "noise levels" have very different meaning here, associated with the M noise channels of the factorial kernel. All noise levels (which we took to be equal in later parts, a meaningless choice in other methods) encapsulated in single energy/score function. Using Langevin MCMC we only sample highly noisy data y. Clean data in Rd is generated via a single step by computing Byy).
8 CONCLUSION
We approached the problem of generative modeling by mapping a complex density in Rd to a smoother "dual" density, named M-density, in R MM. Permutation invariance is design choice which we found particularly appealing. Using factorial kernels for density estimation and generative modeling have never been explored before and this work should open up many new avenues. We believe the topic of parametrization will play an especially important role in future developments. There is a unity in the parametrization schemes we proposed and studied in the paper, but much remains to be explored in understanding and extending the relationship between MEMs and MSMs.
Our MDAE parametrization scheme (an instance of MSM) is especially appealing for its simplicity and connections to the past literature. DAEs have a very rich history in the field of representation learning. But research on them as generative models stopped around 2015. Our hypothesis that one must use very large noise to make them work as generative models for complex datasets in high dimensions. Our permutation invariant multimeasurement approach an elegant solution since allows for that choice at the cost of computation (large M), with only one parameter left to tune: af Crucially, the probabilistic interpretation ingrained here, given by the algebraic relations between the MDAE output p the Bayes estimator (y)), and the score function V log p(y).
3 'There are similarities between Eq. in (Alain & Bengio, 2014) and Eq. 18 here. However, one is on barely noisy data (a + 0) in ddd the other on multimeaaurementss in H Md for any (high) noise level a. There they arrived at Eq. 5 (with some effort) starting from a DAE objective. Here, we start with Fq. 18 (it defines the score function) and arrive at MDAE objective (Eq. 19), in one line of algebra.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENT
We would like to thank Aapo Hyvärinen, Francis Bach, and Faustino Gomez for their valuable comments on the manuscript, and Vojtech Micka for help in running experiments.
ETHICAL CONSIDERATIONS
By their very nature, generative models assign a high (implicit or explicit) likelihood to points in their training set. When samples are generated from a trained generative model, quite possible that some samples are very similar or identical to certain data points in the training set. This is important to keep in mind when the data used to train a generative model is confidential contains personally identifiable information such pictures of faces in the FFHQ-256 dataset (Karras et al., 2019). Care should be taken to abide by license terms and ethical principles when using samples from generative models. For FFHQ 256 in particular, please terms of use ht ttsps: //githtn thtb.uu..... com Nom/NVlan aalabb/// / ffq-ddatsse /entt//luu /wllb//bb/////////oo////oo hoo/mmtterr /REAEEEE .md. Additionally, we recommend considering the in-depth discussion this subject by Prabhu & Birhane (2020) before deploying applications based on generative models.
REPRODUCIBILITY STATEMENT
All important details of the experimental setup are provided Appendix E. Our code is publicly available at bttps://gittiie wwwmnnaisense//eemnmaaaaaaaeeeee mmmms.
REFERENCES
Guillaume Alain and Yoshua Bengio. What regularized auto-encoders learn from the data-generating distribution. The Journal of Machine Learning Research, 15( 15(1):3563-3593, 2014.
Philip Anderson. More is different. Science, 177(4047):393-33 1972.
Christophe Andrieu, Nando De Freitas, Arnaud Doucet, and Michael I Jordan. An introduction to MCMC for machine learning. Machine learning, 50(1-2):5 2003.
Michael Arbel, Liang Zhou, and Arthur Gretton. Generalized energy based models. arXiv preprint arXiv:2003.05033, 2020.
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new perspectives. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(8):1798-18) 1828, 2013.
Nanxin Chen, Yu Zhang, Heiga Zen, Ron Weiss, Mohammad Norouzi, and William Chan. Wavegrad: Estimating gradients for waveform generation. arXiv preprint arXiv:2009.00713 2020.
Xiang Cheng, Niladri S Chatterji, Peter L Bartlett, and Michael I Jordan. Underdamped Langevin MCMC: A non-asymptotic analysis. Conference on Learning Theory, pp. 300-323. PMLR. 2018.
Rewon Child. Very deep VAEs generalize autoregressive models and can outperform them on images. arXiv preprint arXiv:2011.10633 2020.
Yilun Du and Igor Mordatch. Implicit generation and generalization in energy- models. arXiv 2019.
preprint 1903v
Yilun Du, Shuang Li. Joshua Tenenbaum, and Igor Mordatch. Improved contrastive divergence training of energy based models. arXiv preprint arXiv:2012.01316. 2020.
Stefan Elfwing. Eiji Uchibe, and Kenji Doya. Siegoid-weighted linear units for neural network function approximation in reinforcement learning. arXiv preprint 1r02v:1702.03/11 2017.
William Falcon et al. Pytorch Lightning. URL httt:s///gghhbbb..m///tttttttttt htg.
ingg.
Ruiqi Gao, Yang Song, Ben Poole, Ying Nian Wu, and Diederik P Kingma. Learning energy-based models by diffusion recovery likelihood. arXiv preprint arXiv:20/2.08/22 2020.
10

---

Published as a conference paper at ICLR 2022
Ziv Goldfeld, Kristjan Greenewald, Jonathan Niles- Weed, and Yury Polyanskiy. Convergence of smoothed empirical measures with applications to entropy estimation. IEEE Transactions on Information Theory, 66(7):436--438 2020.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition, In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770-778, 2016.
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. GANs trained by two time-scale update rule converge local nash equilibrium. Advances in neural information processing systems, 30, 2017.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. arXiv preprint arXiv:2006. 1239, 2020.
Aapo Hyvärinen.
Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(Apr):695-709, 2005.
Alexia Jolicoeur-Martineuu Rémi
Rémi Tachet des Combes, and Ioannis Mitliagkas. Adversarial score matching and improved sampling for image generation. arXiv preprint arXiv:2009.05475. 2020.
Michael I Jordan, Zoubin Ghahramani, Tommi Jaakkola, and Lawrence Saul. An introduction to variational methods for graphical models. Machine learning, 37(2):183-25) 1999.
Zahra Kadkhodaie and Eero Simoncelli. Solving linear inverse problems using the prior implicit in denoiser. arXiv preprint arXiv:2007.136400 2020.
Tero Karras, Samuli Laine, and Timo Aila. A style-based generator architecture for generative adversarial networks. Proceedings of the IÉEE/CVF Conference on Computer Vision and Puttern Recognition, pp. 4401-4410, 2019.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv:1412.6980 2014.
Diederik P Kingma and Max Welling, Auto-encoding variational Bayes, arXiv preprint arXiv:1312.6114 2013.
Alex Krizhevsky. Learning multiple layers of features from tiny images. Master's thesis, University of Toronto, 2009.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86 86(11):228--2324 1998.
Zengyi Li, Yubei Chen, and Friedrich T Sommer. Learning energy-based models in digh-dimessional spaces with multi-scale denoising score matching. arXiv preprint arXiv:1910.07762. 2019.,
Min Lin, Qiang Chen, and Shuicheng Yan. Network in network. arXiv preprint arXiv:1312.4400 2013.
Koichi Miyasawa. An empirical Bayes estimator of the mean of normal population. Bulletin of the International Statistical Institute, 3814)1181-188 1961.
Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral normalization for generative adversarial networks. arXiv preprint arXiv:1802.05955 2018.
Peter Mörters and Yuval Peres. Brownian motion, volume 30. Cambridge University Press, 2010.
Erik Nijkamp. Mitch Hill, Song-Chun Zhu, and Ying Nian Wu. Learning non-convergent non- persistent short-run MCMC toward energy-based model. arXiv preprint arXiv: 1904,0970, 2019.
Erik Nijkamp. Ruiqi Gao, Pavel Sountsov, Srinivas Vasudevan, Bo Pang, Song-Chun Zhu, and Ying Nian Wu. MCMC should mix: Learning energy-based model with flow-based backbone. In International Conference Learning Representations, 2022.
11

---

Published as a conference paper at ICLR 2022
Emanuel Parzen. On estimation of probability density function and mode. The annals of mathemat- ical statistics, 33(3):0065-1076, 1962.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch. 2017.
Vinay Uday Prabhu and Abeba Birhane. Large image datasets: A pyrrhic win for computer vision? arXiv preprint arXiv:2006. /69223,36.66,,, 2020.
Xuebin Qin, Zichen Zhang, Chenyang Huang. Masood Dehghan, Osmar Zaiane, and Martin Jagersand. U2-Net: Going deeper with nested U-structure for salient object detection. volume 106, pp. 107404, 2020.
Prajit Ramachandran, Barret Zoph, and Quoe Le. Swish: self-gated activation function. arXiv preprint arXhv:1710.05944 7, 2017.
Martin Raphan and Eero P Simoncelli. Least squares estimation without priors supervision. Neural Computation, 23(2):374-420, 2011.
Herbert Robbins. An empirical Bayes approach to statistics. In Proc. Third Berkeley Symp., volume 1, pp. 157-163, 1956.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer- assisted intervention, pp. Springer, 2015.
Matthias Sachs, Benedict Leimkuhler, and Vincent Danos. Langevin dynamics with variable coeffi- cients and nonconservative forces: from stationary states to numerical methods. Entropy, 19(12): 647,2017.
Tim Salimans and Jonathan Ho. Should EBMs model the energy or the score? In Energy Based Models Workshop-ICI.8 2021, 2021.
Saeed Saremi. On approximating V with neural networks. arXiv preprint arXiv:1910. 127444, 2019.
Saced Saremi. Learning and inference in imaginary noise models. arXiv preprint arXiv:2005.0900 2020a.
Saced Saremi. Unnormalized variational Bayes. arXiv preprint arXiv:2007.151333 2020b.
Saeed Saremi and Aapo Hyvärinen. Neural empirical Bayes. Journal of Machine Learning Research, 20(181):1-23, 2019.
Saeed Saremi, Arash Mehrjou, Bernhard Schölkopf, and Aapo Hyvärinen. Deep energy estimator networks. arXiv preprint arXiv: 2018.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. arXiv preprint arXiv: 1907:v/0000000000 2019.
Yang Song and Stefano Ermon. Improved techniques for training score-based generative models. arXiv preprini arXiv:2006.09011, 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456. 2020.
Rupesh Kumar Srivastava, Klaus Greff, and Jürgen Schmidhuber. Training very deep networks. In Proceedings of the 28th International Conference on Neural Information Processing Systems- Volume 2, pp. 2377-2385, 2015.
Terence Tao. Topics in random matrix theory. American Mathematical Society, 2012.
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Computa- fion, 23(7):1661-1674 201 I.
12

---

Published as a conference paper at ICLR 2022
Martin J Wainwright. High-dimensional statistics: non-asymptotic viewpoint, volume 48, Cam- bridge University Press, 2019.
Norbert Wiener. Extrapolation, interpolation, and smoothing of stationary time series: with engi- neering applications, volume 8. MIT press Cambridge, MA, 1964.
Jianwen Xie, Yang Lu, Song-Chun Zhu, and Yingnian Wu. theory of generative convnet. In International Conference Machine Learning, pp. (635-2644, PMLR, 2016.
Jianwen Xie, Yang Lu, Ruiqi Gao, Song-Chun Zhu, and Ying Nian Wu. Cooperative training of descriptor and generator networks. IEEE transactions on pattern analysis and machine intelligence, 42(1):27-45, 2018.
Jianwen Xie, Zilong Zheng, and Ping Li. Learning energy-based model with variational auto-encoder as amortized sampler. In The Thirty-Fifth AAAI Conference on Artificial Intelligence (AAAI), volume 2, 2021.
Yang Zhao, Jianwen Xie, and Ping Li. Learning energy-based generative models via coarse-to-fine expanding and sampling. In International Conference Leurning Representations, 2020.
13

---

Published as a conference paper at ICLR 2022
A GAUSSIAN M-DENSITIES
We study the Bayes estimator for Gaussian MNMs analytically for X = of, where μe IR As in the example for Poisson MNM in the paper, the first step is to derive an expression for the joint density
ply) / II
Pm (Jm (x))d.)
(24)
We start with M 2. Next we perform the integral above:
1111/ + 03) U2 100 + oî) + |1|12(00 + +
log ((v.()))))
2 + 0200 + oboo)
(25)
2(th 12) 3/,,32)0f + 2/,,1100
log Z(0 02,
2(0703 3 + 0220 o
where Z is the partition function:
Z(oo, σ1. 02) (100)000))))) +o 0 4/2
The multimcasuuemment Bayes estimator of X is computed next via Eq.
f y1 y20got o, + ob0z + t
(yn, U2)
(26)
As sanity check, F(y) is the same whether Eq. 7 is computed using m = 1 or 2. Another sanity check the limit σi - 0, which we recover ( y) = I. The linear relation (linear in y for μ 0) observed in Eq. 26 is generalization of the well known results for M known as Wiener,fiter (Wiener, 1964).
visualization of this result for a mixture of Gaussian is given in Fig. 1. Next we state the results for general M, and after that we derive the expression for the Bayes estimator 2(y).
is convenient define the following notations:
ΣM
II m==
oh ZM, where ZM
II o m=0
Σ 1000
M
Σ Ma where ym == Um
II [e[M]Am
of and yo == μ.
m==
With these notations, the expression for the M-density is given by:
logply' Σ
dlog (2T) M VVVVNNN) %MM)
(27)
2x
2
We now derive 2(y) using Eq. 7. Given the expression above, the algebra is straightforward:
a(y) ym + a Vm logp(y)
Um + o2 (2E > 1V
2
Vm
M
(28)
(2 II of ZM) - -220
II
of ū
t
LE[M]]m
g
Note that the final expression is the same by carrying the algebra for any measurement index
14

---

Published as a conference paper at ICLR 2022
B M-DENSITY SCORE MATCHING = MULTIMEASUREMENT DSM
Here we provide the score matching formalism (Hyvärinen, 2005) on M-densities and we arrive at multmmassuremett denoising score matching (MDSM) learning objective. The derivation closely follows (Vincent, 2011) using our notation. In what follows equalities are modulo additive constants (not functions of parameters 0) which we denote by the color gray:
a(0) Ey Vy.oly) Vy log p(y))
Ey Voy 2 + VV Vy log logp(y))+ V log p(y))
Ey Vy 112 +2 / p(y) V(y)(Vy,o()), Vy logply))dy
Fy VyJooy>) 11 +2 / (Vyb(y), vyp(y))
Ey Vyfo(y) +2 2 / Vyff(l)) Vyp(yl=))p(()dzz)
Ey Vyfo(yy 11 + 2 / (Vyfoly), ((y)) (otla///ooopolll))yoaoaayaaaaaa/
Ey fy)o 112 + E VV v logp(y/x))
F( ( Vy 2 R (x (V (x)):: Vy + Vy logp(y/(z))
E(ayl vy Vy logg((y/<>)11
(29)
The M-density score matching learning objective is given by the first identity above and the last equality is the MDSM objective:
8(0) B V fely) Vy b
(30)
For Gaussian MNMs, the MDSM learning objective simplifies to:
M
a(0) Σ 3m (0), where o. = E(xx)) mofoly) + 7/m x
(31)
m==
C MUVB: AN UNNORMALIZED LATENT VARIABLE MODEL
In this section we give formulation of the energy function parametrization scheme which is named mullimeasurement unnormalized variational Bayes (MUVB). Conceptually, MUVB motivated by the goal of connecting empirical Bayes' denoising framework with variational Bayes' inference machinery. There seem to be fundamental challenges here since one model is based on least- squares estimation, the other based on the principle of maximum likelihood. In what follows, this "fundamental conffict" takes the shape of latent variable model that is unnormalized. There have been recent sudies along these lines (Saremi, 2020a:b) that expand on; in particular introduce the novel metaencoder. At a high level, the idea here is to set up latent variable model for M-densities and the variational free energy as the energy function. We refer to (Jordan et al., 1999) for an introduction to variational methods. In what follows, we use the lighter notation for a @ M models (Remark 2) that we also used in Sec.
In its simplest form, the latent variable model for M-densities is set up by the following choices:
The prior over latent variables which we take be Gaussian
p(2) = N(2; 0, Ta,)
The approximate posterior over latent variables (=(=)y) which is taken to be the factorized Gaussian, a standard choice in the literature.
For the conditional density P., (yl (z) we can consider the following
py (y|2) N(y; (2), (T Tasa):
15

---

Published as a conference paper at ICLR 2022
which may seem especially "natural" choice since samples from the M-density obtained by adding multimeasurement noise to clean samples.
In the variational autoencoder parlance, & are the encoder's parameters, l the decoder's parameters, the decoder's 's outputs. With this setup, the variational free energy is easily derived which we take to be the energy function as follows (0 = 7)):
and V , PM)
fo(y) 2/2 (
1
Mr (=)
2x2
where the KL divergence derived in closed form, and the expectation is computed by sampling (=(y))) via the reparametrization trick (Kingma & Welling, 2013).
In this machinery, in the inference step we are forced to have normalized (and tractable) posteriors since must take samples from (= What about the conditional density p(y[z)? If we were to formalize a VAE for modeling M-densities we had to have normalized conditional density p(y|z) since that framework is based on the principle of maximum likelihood. What is intriguing about our setup is that since our goal is to learn the unnormalized p(y), we can consider conditional densities p(y 2) that are unnormalized. An intuitive choice is the following that we name metalikelihood:
Pin.c) (y|z) x exp (-2Ø2 lly (n(2)112 (c(())o(z)))
Now, the energy function takes the following form (A = (6. C)):
soly) Exgl))
2o 20l -
(L(gg(zly))ppa))
(32)
This the final form of MUVB energy function. We refer to hę by metaencoder. Its input is (y, n()) (note that = itself is function of y via inference), and can simply use any encoder architecture to parametrize it. In our implementation, we used the standard VAE pipeline and we reused the encoder architecture parametrize the metaencoder.
D TWO WALK-JUMP SAMPLING ALGORITHMS
In this section we provide the walk-jump sampling (WJS) algorithm based on different discretizations of underdamped Langevin diffusion (Eq. 23). The first due to Sachs et al. (2017) which has also been used by Arbel et al. (2020) for generative modeling, The second is based on a recent algorithm analyzed by Cheng et al. (2018) which we implemented and give the detailed algorithm here. In addition, we extended their calculation to general friction (they set y 2 in the paper) and provide an extension of the (Cheng et al., 2018, Lemma 11) for completeness (the proof closely follows the calculation in the reference). In both cases we provide the algorithm for M models.
D.1 WALK-JUMP SAMPLING ALGORITHM
Walk-jump sampling (Saremi Hyvärinen, 2019) using the discretization of Langevin diffusion by Sachs et al. (2017). The for loop corresponds to the dashed arrows (walk) Fig. and line ((-) is computed over the measurement indices m)) to the solid arrow (jump).
Algorithm 1:
1: Input 5 (step size), (inverse mass), y (friction), K (steps taken)
2: Input Learned energy function f(y) or score function g(y) = V dogp(y)
3: Ouput Xx
4: Yo Unif( 10 1]
5: Vo +-
6: for = [0, [0,...,K) do
7: YET1 f- -V++0V/22
8: V++1 -VVy/(Y or V k+1 + g(Y++1)
9:
V t V++ ud k+1/2
10: BR+1 ~ N(0,IMa)
11:
VeLI + exp(-70) V 4++1 1
Vu exp(-270) B
Yo41 Y VV++122
13:
end for
Xx
(YK.m- ,2 fm or Ax + (Ykm m + o² (Y
16

---

Published as a conference paper at ICLR 2022
D.2 WALK-JUMP SAMPLING ALGORITHM II
Before stating the algorithm, we extend the calculation in Cheng et al. (2018, Lemma 11) to general friction, which used in implementation of their Langevin MCMC algorithm. The algorithm is stated following the proof of the lemma which closely follows Cheng et al. (2018).
The starting point is solving the diffusion Eq. 23. With the initial condition (yo, vo). the solution (ye, ve) to the discrete underdamped Langevin diffusion is (t is considered small here)
Vr = Voe voe (Le ( + /2Tu L. e (
(33)
Yr yo + ( Veds,
These equations above are then used in the proof of the following lemma. Lemma Conditioned on Vy), the solution of underdamped Langevin diffusion (Eq. 23) integrated up time is a Gaussian with conditional means
E[v1] = Voe - u(1 - e Vf(()
(34)
y yo + + - - "'voo - (1- Vf(yo),
and with conditional covariances
E((,, (xr - E[y.)) 1 7 - (21 47 47-44-1000 -
-1- - e + 7 E -2))) IMds
E[[vv (vv.))(v, B[v.)) =u(1 - -21) ) Insa,
(35)
[ [E(y)))(v B[v.])]] - [[1
2e rt + 21t). Imdd-
Proof. Computation of conditional means is straightforward as the Brownian motion has zero means:
E[v:] = Voe > - - - - - - Vf(())
Ely, = yo + + - - ">o - 7
7 (1 (1-e f((o),
All the conditional covariances only involve the Brownian motion terms. We start with the simplest:
E (v; - - - = 2 E (f e--llt-s) dB.) ( e dB.
2 % e-27(t-o) Ins
(36)
= =u(1- - -20) Iasd
From Eq. 33, the Brownian motion term for yr is given by
v 27 (1% (--((--)) dB, dr = 22 (- eT (f" e- dB.
2 1'' (1- dB..
The conditional covariance for y, follows similar to Eq. 36:
B( Ely,))(y+ E[y.)) 2x -/u (J,² (1 e -((--s) 2 ds -nsa
(37)
-- (2t 47 -1 E =7 + - (11 - e e22))) -ard
Finally the conditional covariance between yr and v. is given by
E E Ely.)) B[v.)) 2u [ (1 (1----10000) dB.) (1" ---lt--) ddB
2u (f/² (1 (1- (1- e -11l-0) TMd
= 1 - -4 + t LMd
17

---

Published as a conference paper at ICLR 2022
Algorithm 2: Walk-Jump Sampling (WJS) using the discretization of Langevin diffusion from Lemma 1. Below, we provide the Cholesky decomposition of the conditional covariance matrix.
1: Input 5 (step size), (inverse mass), Y (friction), K (steps taken)
2: Input Learned energy function f(y) or score function g(y) R V Vlgg((y)
3: Ouput Xx
4: Yo Unif( (0, 1jaad)
5: Vo f0
6: for k ......,. K) do
7: VA +- - VyJ(YR) v t g(Y))
8: VA+1 t VAe +1 u(1 e e-7 y
9: YRUT & Y (1- VV + + u (6 71 - (1--0))) V
-
10: BR+1 ~ ((0, IZMd)
YR+1
Ys+1
+ +LB++1
11 see Eq. 38
Vk+1
Va+1
12:
end for
Xx f (YKm - (mvm((K)) or Xx + (Y + om (Ym))
13:
The Cholesky decomposition E = LL for the conditional covariance in Lemma 1 is given by:
ENy 1/2
1/2 .Mdd
L
(38)
Egy
Eny: IMdd (Evv - -Eye E -Iad)
Eyy
where
Eyy = --'u (28 - 47 (1 - 110) -E -2100)
(39)
Eiv = u(1 - 1---6-
(40)
Eyv = 7 1-1 2e 21 + 2266
(41)
Algorithm vs. Algorithm The two algorithms presented in this section have very different properties in our high dimensional experiments. In general found Algorithm 1 (A1) to behave similarly to overdamped version (see Appendix F for some illustrations) albeit with much faster mixing. We found friction 7 to play very different role in comparing the two algorithms, also apparent in comparing the mathematical expressions in both algorithms: in A1, the friction only appears in the form ex)(->)) it more complex A2, Similar discrepancy found regarding u. Remark 4 (Initialization scheme), The initialization Yo ~ Unif[[0, (0, 1jāad) used in the algorithms is the conventional choice, but for very large we we found effective to initialize the chain by further adding the Gaussian noise in RR Md.
Yo Eo + Eo: where Eo 1] O'0,
Otherwise, with large step sizes 8 = O(a) the chains starting with the uniform distribution could break early on. We found this scheme reliable, and well motivated since M-density p(y) is obtained by convolving p(x) with the Gaussian MNM. With this initialization one starts "relatively close" " to the M-density manifold and this is more pronounced for larger a. However, this initialization scheme is based heuristics (as is the conventional one).
Remark 5 (The jump in WJS is asynchronous). The jump in WJS is asynchronous (Fig. 2) and it can also be taken inside the for loop in the algorithms provided without affecting the Langevin chain.
Remark (The Ak notation). In running long chains we set to be large, and we report the jumps (taken inside the for loops in the algorithms provided) with certain fixed frequency denoted by Åk. Remark (WJS for Poisson M-densities). WJS is a general sampling algorithm, schematized in Fig. 2. However, for Poisson MNMs it will involve sampling a discrete distribution in Nadd which we do not know how to do efficiently in high dimensions (note that Langevin MCMC cannot be used since the score function is not well defined for discrete distributions). We can use Metropolis- Hastings algorithm and its popular variant Gibbs sampling but these algorithms are very slow in high dimensions since at each iteration of the sampling algorithm in principle only I dimension out of Md shall be updated, so in the general case the mixing time is at best of (Md); one can use "blocking" techniques but they are problem dependent and typically complex (Andrieu et al., 2003).
18

---

Published as a conference paper at ICLR 2022
E EXPERIMENTAL SETUP
Datasets Our experiments were conducted on the MNIST (LeCun et al., 1998). CIFAR-10 (Krizhevsky, 2009) and FFHQ-256 (Karras et al., 2019; pre-processed version by Child, 2020) datasets. The CIFAR-10 and FFHQ-256 training sets were augmented with random horizontal flips during training.
Network Architectures For all experiments with the MDAE and MEM parametrization, the Ne network architecture (Qin et al., 2020) was chosen since a recent and effective variant of the UNet (Ronneberger et al., 2015) that has been successfully applied to various pixel-level prediction problems. However, we have not experimented with alternatives. Minor adjustments were made the Net from the original paper: removed all normalization layers, and switched the activation function to и + I sigmoid(z) a:) (Elfwing et al., 2017; Ramachandran et al., 2017). To approximately control the variance of inputs the network in each dimension, even if the noise level divided by scaling factor v0.225² a². We also adjusted
is large, the noisy input values the number of "stages" in the network, stage heights and layer widths (see Table 1 for details), in order fit reasonable batch sizes on available GPU hardware and keep experiment durations limited to few days (see Table 2).
For experiments with the MUVB parametrization (the MNIST model in Sec. G.3), we used Residual convolutional networks utilizing bottleneck block design with skip connections and average pooling layers (Lin al., 2013; Srivastava al., 2015; He et al., 2016). The architecture of the encoder and metaencoder were kept the same, while the decoder used the encoder architecture in reverse, with the pooling layers substituted with nearest neighbour upsampling. See Table 3 for a detailed description of the encoder architecture.
Remark 8. Strictly speaking, for a @ M models the energy function is permutation invariant and should be treated {y1...., yM }. a set consisting of M Euclidean vectors. This however difficult to enforce while keeping the model expressive. We note that in our model the permutation symmetry "fragile symmetry" which one can easily break by making om just slightly different from each other. Therefore, designed the architecture to be used for the general setling of am different from each other, and in experiments we just set them to be equal.
Training All models were trained using the Adam (Kingma & Ba, 2014) optimizer with the default setting 3 = 0.9, 33 = 0.999 and e = 10-8, Table lists the main hyperparameters used and hardware requirements for training MDAE models for each dataset. The resulting training curves shown in Fig. 6, showing the stability of the optimization with a relatively simple setup.
Software All models were implemented in the CPython (v3.8) library PyTorch v1.9.0 (Paszke et al., 2017) using the PyTorch Lightning framework v1.4.6 (Falcon et al.).
Table 1: Net architecture used for MDAE and 2EM parametrizations for all datasets. All layer widths (number of channels) were expanded by a width factor for CIFAR-10 and FFHQ-256.
Stage
Height
In channels
Mid channels
Out channels
Side
Encoder
3 or
64
64
Encoder 2
5
64
128
128
Encoder
128
128
256
Encoder 4
3
256
256
512
Encoder
512
256
512
512
Decoder 4
1024
128
256
256
Decoder 3
512
128
128
128
Decoder 2
256
128
64
64
Decoder
128
64
64
64
19

---

Published as a conference paper at ICLR 2022
20
50
1.500
Loss
1400
20
1000
10
200
400
800
800
1000
200
400
600
500
1000
50
100
1.00
200
250
Epoch
Epoch
Epocth
(a) MNIST & 16
(h) CIFAR-10 10
(c) FFHQ-256 60
Figure 6: MDAE training plots demonstrate that optimization is stable for all datasets. Loss values on the y-axis in these plots represent the value of our objective function (Eq. 14).
Table 2: Main hyperparameters and computational resources used for training MDAE models.
MNIST & 16
CIFAR-10 @
FFHQ-256 4 @ 8
Width factor
2
Learning rate
0.00002
0.00002
0.00002
Total batch size
128
256
GPUs
GTX Titan X
4xGTX Titan
V10000
Run Time
R=2.3 days (1000 epochs)
R22 days (1000 epochs) Rv2.75 days (250 epochs)
Table 3: MUVB Encoder architecture for MNIST. Each row indicates sequence of transforma- tions in the first column, and the spatial resolution of the resulting output in the second column. Block(N] x (D) denotes a sequence of blocks, each consisting of four convolutional layers with skip connection using "bottleneck" design: the layers have widths [N, 0.25 N, 0.25 N, NJ and kernel sizes 11,333, 1], except when the input resolution becomes smaller than in which case all the kernel sizes are 1.
Module
Output resolution
Input
x28
Block128x2
28x28
AveragePool
14x 14
Block256x2
14x
AveragePool
7x7
Block512x2
7x7
AveragePool
4x4
Block 1024
4x4
AveragePool
2x2
2x2
AveragePool
1x1
Block1024xl
1x1
20

---

Published as a conference paper at ICLR 2022
F LANGEVIN MCMC LABORATORY
In this section, we give light tutorial on Langevin MCMC used in the walk phase WJS. Perhaps the most interesting result enabled by multimeasurement generative models is that Langevin MCMC in high dimensions becomes stable enough for us to examine and dissect the performance of various Langevin MCMC variants in practice.
The reader might ask why we used underdamped Langevin MCMC versus its well-known overdamped version. The reason is the recent theoretical developments discussed in Sec. 5 that show much faster dimension dependence for the convergence (mixing time), O O(vi) Langevin MCMC. To give an idea regarding the dimension, for the & 8 model on FFHQ-256, Md = 10 Our trained model for 4 @ 8 is used throughout this section. In addition, the MCMC is initialized with the fixed random seed throughout this section.
O(d), for underdamped
F.1 OVERDAMPED LANGEVIN MCMC
We start with overdamped Langevin MCMC which is simpler and give intuitions the behavior of Langevin MCMC in general. The algorithm is based on the following stochastic iterations (setting inverse mass u 1):
Yk+1 t Yk 2
vf(yx)
+ 1. Ek. where EN ~ N(, Iara)
(42)
B
Comparing above with Eq. 7, it is intuitive to consider setting 8 Oo)). We start with ó in this case Aş pulls noisy data to clean manifold (note the important extra factor of 1/2 however) and B& term corresponds to adding multimeasurement noise sampled from N(0, a 2 TMa). The overall effect is that one starts "walking" on the manifold of the M-density. And if that is the case, the jumps computed by T(yx) should give samples close the manifold of the (clean) data distribution samples from p(r). We start with that experiment by setting 5 = a. Note that in walk-jump sampling we sample the M-density in Ma generating highly noisy data for
we only show clean (jump) samples in Rd. This is shown next by skipping every 30 steps, i.e. Ak 30:
Our intuition partly correct in that the sampler manages to "capture modes" but the step size is clearly too high. Now consider 8 = 0/2 and Ak = 100:
Continuing the chain, now shown for Ak 1000, we arrive at
Theoretically, converge the distribution (associated with the M-density) the step size should be small, where "small" is dictated by how "close" one wishes to be to the true distribution, e.g. (Cheng al., 2018, Theorem 1). As we see above, even for after 1000 steps, but a fixed 6 = a a/2 appears to be too high for this model. By continuing with the chain with the step size a/2 one starts gradually move away from the manifold shown next:
the sampler does not converge
21

---

Published as a conference paper at ICLR 2022
F.2 UNDERDAMPED LANGEVIN MCMC
In this section we focus on underdamped Langevin MCMC. The algorithms are based on discretizing underdamped Langevin diffusion (Eq. 23). In this paper, considered two such discretization schemes, by Sachs et al. (2017) that we used in Algorithm I, and the one due to Cheng et al. (2018) that we implemented and used Algorithm 2. The first algorithm is not analyzed but very easy to show that in the limit 7 - 3o it will be reduced to the overdamped Langevin MCMC we just discussed. Second algorithm is more complex in how friction ? and inverse mass behave. For comparison, we present results obtained by both algorithms, showing Algorithm / first. As in the previous section, we start with 0; we set 7 = 1 (u unless stated otherwise) with Ak = 30:
As in the previous section ò = a is too high for FFHQ-256, 4 & 8 model. (However, note the relative stability of Algorithm for this random seed.) Next, we consider 0/22, 7 1, shown with Ak = 100:
Note the remarkable contrast between the results above and the corresponding figure (for Ak 100) in the previous section. We continue the chain with the same Ak = 100:
Both algorithms start to diverge for this choice of parameters (5 = 0/2,7= 1). add more stability we can consider increasing the "mass" (lowering u). In (Cheng et al., 2018), u is replaced with 1/L. where L is the Lipschitz constant the score function. By lowering u. we are telling the algorithm that the score function is less smooth and MCMC uses more conservative updates. This can also be seen in the equations in both algorithms, where the step size 5 multiplied by u in several places. However, note this effect is more subtle than just simply scaling 8, especially for Algorithm 2 as can be seen by inspection. Next, we consider u = 1/2. This will affect the mixing time and the chains now shown with Ak 1000 (instead of Ak 100 earlier):
This experiment concludes this section. We come back to more FFHQ-256 chains in Appendix I.
22

---

Published as a conference paper at ICLR 2022
G MNIST
G.1 MORE IS DIFFERENT
There is notion of congruence which arise in our analysis in Sec. 2.3 which we denote by "=": two different models, a @ M and o' @ M', "agree with each other" in terms of the plug-in estimator
x M Mf
0
VM
VM'
However these models are vastly different, one is in I the other in R And if o >> o' we expect that a @ M model to have better properties as a generative model. This clear in the regime (a' < 1, M' = 1). In that regime, almost no learning occurs: at test time MCMC either gets stuck in a mode for a long time (the best case scenario) or simply breaks down. We refer to (Saremi & Hyvärinen, 2019) for a geometric intuition why large noise is essential in these models regarding the walk-jump sampling. We should point out that this is a general phenomenon in all denoising objectives (Alain & Bengio, 2014; Song & Ermon, 2020). In nutshell, when noise is small the Langevin sampler finds it difficult to navigate in high dimensions and the forces of the Brownian motion eventually takes first experiment below). If we are "lucky" we can find the manifold initially (the second experiment below) and MCMC may remain stable for while but model with o' will have better mixing properties. We should highlight that in the analysis of Cheng et al. (2018) the mixing time scales as R2 where i is the condition number (see Sec. 1.4.1) proportional to the Lipschitz constant: this quantifies the intuition that MCMC on smoother densities (smaller Lipschitz constant) mixes faster.
part of space unknown our model and MCMC breaks down (the
To validate these intuitions we consider comparing 1 @ 16 and 1/4 @ In addition to statistically identical plug-in estimators, these models also arrive at similar training losses with the same training schedules (1.03 for & 16 and 1.19 for 1/4 @ 1). We start with 1/4 @ 1. Below we show WJS using 1/2,
1) (the setting ā a 0/2,1 y 1/2, 1 is simple choice we used in testing all our a @ M models). Here, MCMC breaks down shortly after the last step shown here (this chain is shown with Ak 30, in the remainder Ak- 500):
Algorithm 2 with the setting (6 1/s,
2
Algorithm is more stable here and we arrive at the following chain (total of 10 steps):
9 7
O
0 O
0
6
8 4
8
0
Now, for 1 @ 16 model, using Algorithm 2 with parameters (ō 1/2,7 1/2, R
1),
we arrive at:
3 4 9 0 4
4
3
8
3
3
3
3
K
6
3
9
9
7
9
5
9
9
9
S
4
0
0
0
9
9
9
4
4
4
4
5
9
9
9
4
9
5
9
8
9
0
0
8
9
9
9
8
\
4
9
4
4
4
9
1
9
3 8
8
3
2
8
3
6
8
8
3
3
3
3
3
S
3
8
| 7
5
L
l
9
9
9
6
9 9
6
9
9
4
4
4
8
6
8
$
8
8
8
4
4
4
4
4
9
q
And using Algorithm with the same parameters (and initial seed) we arrive at:
3
4
4 1
\ \ 1 8 3
3
3
5
<
(
6 2 l
4
Y
1 7
4 3
3
3
2
2
2
\
5 5 9
S 9
3 4
/
1
\ 3
3
3
4
Y
q
9
9
9
9 9
9
9
5
5 o
O O
O -
O
0 4
5
(Note the stark differences in terms of mixing between Algorithms and 2 for this model.)
In summary, more is different: quantitative differences become qualitative ones (Anderson, 1972) and this is highlighted here for "congruent" a @ M models that differ in terms of the noise level. If computation is not an issue, one should always consider models with larger but this remains conjecture here in its limiting case (as a and M grow unbounded, while a 9/VM remains "small").
23

---

Published as a conference paper at ICLR 2022
G.2 INCREASING M: TIME COMPLEXITY, SAMPLE QUALITY, AND MIXING TIME TRADE-OFFE
Here we extend the ablation study in the previous section by studying @ M models with a fixed where we increase M. In the previous section studied two models where T has the same statistical properties as measured by the loss (and the plug-in estimator), and we validated the hypothesis that MCMC sampling has better mixing properties for the smoother M-density (with larger o). Now, the question is what if we keep a fixed and increase A/?
As we motivated in Sec. 1, by increasing M the posterior p(rly) concentrates on the mean (yy,, thus increasing the sample quality in WJS. However there are trade-offss: (i) The first is the issue of time complexity which is an open problem here. In our architecture, the . measurements are passed through convolutional layer and after that models with different M have approximately the same time complexity. The problem is this may not be an ideal architecture. Presumably, one should increase the network capacity when increasing M (for the same σ) but by how much we do not know (this clearly depends on how complex px is). Here we ignore this issue and we assume that our architecture has enough capacity for the largest M = 16 considered. (ii) The second trade-off is as we increase M, the price we pay for higher sample quality is that we will have longer mixing times, both in terms of generating the first sample, and more importantly in the mixing between modes. This trade-off is intuitive since .y.,...., JM) becomes more complex for larger M (for fixed o). Below ..... we compare WJS chains (as before, without warmup) for 1 and M 1, 2, 4, 16:
5
5
3
3 3 3
3
3
3
3
3 6
2
2
2 6 6 8 4 2 22
2 2.2 2 9 9 9 9 8
9
9
1
1
1
1
1
1
7
7
7
7
7
7
7
6
6
6
6 6 5 5 5 5 5 3
3
3
3
3
3
3
3
3
3
3
8
9
9
9
9
9
9
9 888888888
3
8 6
8
8 8
6 8 4 2
8888 0 8
8
3
8
8
8
5
5
S
S
S
5
S
0
8 88888684449111111
8
7
7
7
9
9
9 9
9 9
4
4
4
0944406688884
2
2
2
2
2
2
2 2
2 2
2
8 8
2
2
2
2 2
2
2 2 2
6 6
0
6 6
66
55555
555
333
3 3333
3
0555555
5688
5 F 366555552
55
8
1
3
3
1
11
3
10
01
01
01
8
6
0
Figure 7: (mixing time vs. image quality trade-off) WJS chains for <9 M models are presented in
"real time" (Ak = 1) starting from noise (320 steps in total), for M = 1, 2, 4, 16 in order from top pannel to the bottom one. We used Algorithm 2 (8 1/2, 7 = 1/2, u = 1) with the same initial seed. Note that, at the cost of sample quality, all classes are visited in @ 1 in the short chain presented.
24

---

Published as a conference paper at ICLR 2022
G.3 MUVB vs. MDAE
Below we make a one to one comparison between MUVB and MDAE for the @ setting we trained on MNIST. In both the Langevin MCMC algorithm by Sachs et al. (2017) with parameters (ā = 7 = 1/4, 1) (see Appendix D). This an "aggressive" choice of parameters for Langevin MCMC, designed for fast mixing, yet the chains do not break up to 1M+ steps (we stopped them due to computational reasons).
ran single chain for 1 million+ steps. For sampling M-density we used
Below show the chains at discrete time resolution of 5 steps (the first sample, top left corner, is obtained after only 5 steps). This comparison visually demonstrates that MUVB has better mixing and MDAE better sample quality. All digit classes are visited in MUVB steps shown here. MDAE chain is cleaner but the classes {0, 1,5, 6) are visited scarcely.
variety of styles in 4000
8
8
8
8
8
8
8
8
8
8
8
8
8
8
5
3
3
3
3
3
8
8
8
8
8
6
8
8
8
8
9
9
6
5
3 3
3
3
3
3
3
3
3
3
3
3
3
5
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3 3 3
3
3
3
5
5
5
S
5
5
6
5
4 
4
4 +11
/
444444444444444
4G444444444444
2 22222
2 a 2
2
3 03333333333333888335333
5 9
4
4
4
4
5
5 0
O
0 6
6
33333
Figure 8: MUVB, 1 (9 model on MNIST, A = 5 steps.
6 6 005988888
8 8 8 8 0218888655833388815
3
323335333333355598888
3
3
3
3
9 "
7
7
7
7
7
7
9 9
9
9
9
9
9
9
4
4 4 44
4
4
4
4
Y 4
1
9
8
8
8
7777777777777777
9
9
9
3
8 8 1
3
7
9
8
8
8
3
9
3
3
3
3333333333333
3
3
3
3
3
8
& 8 8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
9
8
8
8
8
8
8
8 8 8 8
8
8
8
8
8
8
Figure 9: MDAE, 4 model MNIST, Ak 5 steps.
25

---

Published as a conference paper at ICLR 2022
G.4 LIFELONG MARKOV CHAIN
We demonstrate the WJS chain with I million+ steps obtained in MDAE (1 @ 4) in its entirety (viewed left to right, top bottom). We informally refer these long chains that never break as lifelong Markov chains. The chain is reported here in its "entirety" to demonstrate the fact that it does not get stuck in a mode and remains stable through the end (we stopped due to computational reasons).
2 f 8 7
| 8 22220(22233332223333333333333333
8 
0 "oc""..e66<6333)33333333......."
8 01210111222210011099900000000000
6 2.021222200077022222333777700000
4 ".2",31001222266233331111.1....3
3 04222112222222355552223011111111
2 03220011100001111000443333111110
2 "92.2212288111111133331111111119
O 332451186602122222333333331111"3
8 312112/2222111777776622333111113
2 30051111021133333331000000110010
0 9 0F..00/60.00TTT11111120000100111
3 4 01100118111111110000000000011010
8 0 (2221331222333333333302001111111
2 4 3 115-1111122222223331700009991113
6 0A511000221133887188622000333969
2 y 02215515932000000000000000333313
0 9 (SS122A1113331111)223331111131
4 (.5"1111666600111111113336999999
7 "FF"11122222211111111122111""1""
6 01223115552222220000000111111113
2 04421665501111233102222333333311
4 3r211111111111111333336111111111
2 1772-18824422111199999999999
5 21200011161113333333333333333363
3 388 3 0
4 4 +1
Figure 10: MDAE, 1 8 model on MNIST, A = 500 steps, million+ steps in total.
26

---

Published as a conference paper at ICLR 2022
H CIFAR-10: SINGLE-CHAIN FID EVALUATION & SOME FAILURE MODES
The main "fear" in parametrizing the score function directly in MSM, which MDAE an instance of, we not guaranteed learn conservative score function (see Sec. 4.1). More importantly, we do not have control over how the learned model is failing in that regard. These fears were realized in our CIFAR experiments on the 1@8 model presented below. The main results use Algorithm 2. Algorithm 1 simply fails for the Langevin MCMC parameters we experimented with below (for completeness we present short chain the bottom panel). In 2D toy experiments the non-conservative score functions were quantified in (Saremi et al., 2018, Fig. 1h), but this quantification is difficult in high dimensions. However, what is intriguing about the experiments here is that the "non-conservative nature" of the score function is manifested- -tiis is merely conjecture- in the cyclic mode visits apparent the Markov chain.
Figure 11: Top panels use Algorithm 2 at two different checkpoints. The bottom panel uses Algorithm
1. MCMC parameters are (6 0.1, y 2, u 10) for the top and bottom panel, and и = 20 for the middle panel. In all cases Ak 500 (=> 10 steps). The FID scores are respectively 91 & 79.
27

---

Published as a conference paper at ICLR 2022
We encountered the same problem using spectral normalization (Miyato et al., 2018): the chain
below was obtained using Algorithm 2 with the setting (6 = 1/2, Y = 1/2, u = 1) where we got the FID score of 99, We ran the chain for M steps; the first 600 K steps visualized below:
Figure 12: WIS chain for MDAE 8 8, trained with spectral normalization, AA = 500.
28

---

Published as a conference paper at ICLR 2022
We also experimented with training MDAE using SGD optimizer. Using Algorithm with the setting (6 = 0.7. = 1, u = 1) the MCMC chain was stable for M steps where we stopped the chain. We
picked 50 K images at equal interval of Ak 20 from the chain resulting in an FID of 43.95:
Figure 13: The first 13 K steps of a WJS chain with 1 M steps for MDAE 1 & 8 on CIFAR-10 with
A 10. We obtained the FID score 43.95 by selecting 50 K samples skipping Ak = 20 steps.
29

---

Published as a conference paper at ICLR 2022
We finish this section with discussion numerical comparisons the sample quality between MDAE @ 8 model and other recent MCMC-baser methods, summarized Table
The table excludes denoising diffusion models since the MCMC chain in such models are "conditional" in the sense that sampling is via reversing diffusion process based sequence of conditional distributions which is learned during training (see the discussion in Sec. 7). By comparison, all the papers in the table below there is only one energy/score function that is being learned.
The use of the term "long chain" below is based on the current status of the field, chains of order 1,000 steps: the MCMC chains used for generation and reporting FID scores in most prior works has been "short", of order 50 to 100 steps by comparison which we have indicated in a column in the table. The FID scores in those papers were obtained by 50 K short-run parallel chains.
Our work is unique in this literature in that for the first time we report competitive FID scores obtained from single MCMC chain. As emphasized in Sec. 6 the quest for generating very long "life long" chains with good sample quality and mixing properties, as measured by the FID score (Heusel et al., 2017), is a scientific challenge and we believe it is an ultimate test to demonstrate that one has found good approximation to the true energy/score function. We should point out that the only competitive paper here obtaining FID scores with long MCMC chains is by Nijkamp et al. (2022) where the "long chains" have been limited to 2,000 steps and the FID obtained is much worse (78. 77)12) than what we obtained (43.95) with our @ 8 model from single MCMC chain of 1,000,000 steps (see Fig. 13).
Table 4: FID results for unconditional MCMC-based sample generation on CIFAR-10
FID
long chain
single chain
MCMC steps
Xie et al. (2018)
35.25
X
X
N/A
Nijkamp et al. (2019)
23.02
X
X
100
Du & Mordatch (2019)
40.58
X
60
Zhao (2020)
16.71
60
Xic et al. (2021)
36.20
X
50
Nijkamp ct al. (2022)
78.12
X
2,000
MDAE, @ (our work)
43.95
1,000,000
I FFHQ-256
this section, we provide several WJS chains for our MDAE 4 @ 8 model on FFHQ-256 dataset. We refer Algorithm by A1 and Algorithm by A2, and MCMC parameters are listed as (8, 7, u).
A1; (2, 1): AK 300
A2: (2, 1): Ak 200
30

---

Published as a conference paper at ICLR 2022
A2; (1/10, 1, 10); Ak: 500
A2; (4, 1/2, 1/): Ak 150
A1: (2, 1): AK 100
A1; (2, 1); A = 20
31


Reviews:
[{"review_id": "kB4LQivevPk", "reviewer": "Reviewer_dnHy", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "1. The author claim that the smoothed M-density $p_\\mathbf{Y}$ is easier to learn from than the original distribution $p_\\mathbf{X}$ in high-dimensional settings. However, there are no simple examples or existing literature included in the manuscript to support this statement. The difficulty of sampling from a high-dimensional distribution is a core problem in designing generative models, so the authors should consider making this motivation clear and better supported.\n2. In the introduction, the authors stated that the MCMC is considered an art that cannot converges fast for complex distribution, but later in Section 5, the authors adopted Langevin MCMC to sample the M-density. It is encouraged to compare the MCMC with the Langevin MCMC for M-density since the sampling efficiency is one of the selling point of this paper. For example, is the convergence rate of Langevin MCMC faster? The authors should at least show it in the empirical study.\n3. The authors did not provide a good motivation for why considering the Poisson and Gaussian MNMs in Section 2. What about other MNMs? Does the choice of the MNM kernel related to the properties of $p_\\mathbf{X}$? Is there an optimal MNM kernel to pick? In Section 2 the authors focused on the derivation of the closed-form expressions of the MNMs, but the reasoning on why they chose them is more important in my opinion.\n4. The authors proposed several parameterization schemes in Section 4, including MDAE, MEM${}^2$, and MUVB. What are the pros and cons of them over each other? A discussion is needed or at least added in Appendix.\n5. The permutation invariance is one the appealing points of using MNMs, and in the Introduction, the authors stated ``In Sec. 6, we present our experiments on MNIST, CIFAR-10, and FFHQ-256 datasets which were focused on permutation invariant M-densities.''. However, in Section~6, the discussion on the permutation invariant M-densities is missing.\n6. An empirical comparison between MDAE and other generative models is lacking. How to compare the quality of the generated samples? In the illustrations, the authors focused on facial images, and in some figures, the hair or background in the images are not generated well, why? It is unclear about the advantages of the proposed MDAE over other generative methods, and a discussion is needed.\n7. The abbreviation DAE for denoising auto-encoder is used without definition in the abstract. Although DAE is well-known in the field, it is suggested that the authors still defined all the abbreviations before usage in this manuscript. Another example is MCMC for Monte-Carlo Markov chain.", "summary_of_the_paper": "This paper introduced an alternative sampling method, with an application on generative models, by convolving an unknown distribution $p_x$ with a factorial kernel called multi-measurement noise model (MNM). The resulting M-density $p_y$ is smoother (easier to sample from), and is permutation invariant. Two factorial kernels, Poisson and Gaussian MNMs, are introduced for the convolution, and can be connected to Bayes estimator, and the learning of parametric energy and score functions. Two parameterization schemes are proposed for modeling the energy and score functions of the Gaussian M-density respectively. Empirical results on FFHQ-256 dataset are very impressive. The main contribution is the usage of factorial kernels, which is a direct extension of smoothing a density in non-parametric density estimation, and is not used for generative models before.", "summary_of_the_review": "The proposed MNMs and MDAE are very appealing, and novel. I believe this will be a great publication if the authors could please address the issues listed above. In particular, I think the authors should give (1) a better motivation on why this methodology outperforms other generative methods, (2) how to compare different parameterization schemes, (3) a more thorough discussion on the experiments regarding the permutation invariance, efficiency over MCMC or other methods, and (4) companion on the quality of the generated samples with existing methods. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "Qg5XI2Aohh", "reviewer": "Reviewer_mPcN", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper combined DAE and empirical Bayes to learn and sample from some unknown distribution. This method is a novel method and is very different from previous methods that use convolution to smooth the distributions. While other methods are based on weighted noise or annealing, different noise levels are equal in this method.\n\nThe main body of the paper is in general well-written. I have the following questions on the paper.\n1. Sec 2.1 derived the formula for poisson MNM, but it seems the other parts of the paper are mainly based on Gaussian NMN. Are there any further results on Poisson MNM? How does it compare with the Gaussian MNM?\n2. Sec 2.3 gives a concentration rate of the estimator, but it is not straightforward to see how sigma_eff is derived. \n3. Sec 5 mentioned that underdamped langevin is better than langevin because UL scales as sqrt d theoretically. It is also mentioned in the appendix that in your experiments, UL mixes faster. Are there any experimental results on how much faster it mixes? \n", "summary_of_the_paper": "This paper proposed a multimeasurement noise model to learn and sample from some unknown distribution and demonstrated its effectiveness on datasets like MNIST, CIAR-10 and FFHQ-256. It studies the convolution of the distribution with different levels of noise and uses the Bayes estimator in each noise level to recover the original sample. ", "summary_of_the_review": "This paper is in general a good paper. The new generative method it proposed can be of interest to the ICLR community.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "79QmnoGMOTi", "reviewer": "Reviewer_6P94", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "${\\bf \\text{Strengths}}$\n\n1-The problem addressed in this paper (generating new independent samples of a given selection of samples drawn from an unknown distribution) is of great interest to the machine learning community.\n\n2-The approach used (going through a variety where the density estimation is easier) is very original. The intuitions and connection with literature are interesting and allow a new point of view of the algorithms already used (DAE as an example). Several experiments have been designed to justify the relevance of the approach.\n\n3-The article is quite well written and well structured.\n\n4-Although the algorithm uses already known sampling methods and energy function optimization methods, the general idea of the Bayesian estimator is an interesting contribution to the machine learning community.\n\n${\\bf \\text{Weaknesses}}$\n\n1-The motivation of presenting the Poisson MNM is unclear specifically since it is never used afterward and more importantly no motivation to discard it and prefer the Gaussian MNM. More specifically, isn't it to present a simpler version of MNM before tackling the more involved Gaussian MNM? But more importantly what are the disadvantages of Poisson MNM to discard it?\n\n2-The discussion of the higher noise level is also not convincing. Moreover, the advantage/disadvantage of increasing M is not clearly discussed (only the advantages are pointed out). I am quite convinced that there is a great connection between $M$ and $\\sigma$. More concretely, there is no discussion about the time complexity of the algorithm which may be related to increasing $M$. It seems clear that with a low level of noise we do not need to have too many channels measurement (independent samples). Therefore for me, low values of $\\sigma$ and low value of $M$ are comparable to high values of $\\sigma$ and high value of $M$. So in the comparison for MNIST (between two different values of $\\sigma$), it is not fair to change the value of $M$. In an ablation study, you normally need to keep $M$ fixed and change $\\sigma$ since the two are connected.  Having big values of $\\sigma$ is great to smooth well the manifold but in compensation, we need high values of $M$, which may be problematic for time complexity (which is not discussed).\n\n3-In the abstract, I have the impression that the sentence \"Samples from $p(x)$ are obtained by walk jump\nsampling, using underdamped Langevin MCMC to sample from $p(y)$\" is incomplete since we should close the loop and mention that the empirical Bayes estimator is used after this underdamped Langevin.\n\n4-I the abstract, it is mentioned fast mixing. But the term fast is never concretely shown (no running time,...)\n\n5-In the concentration rate of the empirical estimate, it is mentioned that the optimal estimator is expected to give a better rate than the rate of Tao. This is not proved, right? More concretely what is the rate of the optimal estimator? What is the intuition to say that the rate will be better than just a constant?\n\n6-In the definition of the loss $\\sum_{m=1}^{M}\\mathcal{L}^{m}(\\theta)$, why a normalization with $\\frac 1M\\sum_{m=1}^{M}\\mathcal{L}^{m}(\\theta)$ is not used since the large $M$ case is included in the framework? More concretely is the limit $\\lim_{M\\rightarrow\\infty}\\sum_{m=1}^{M}\\mathcal{L}^{m}(\\theta)$ correctly defined?\n\n7-In section Denoising score matching it is written 'empircal' instead of empirical\n", "summary_of_the_paper": "Given $n$ independent samples $x_i$ in a space of dimension $d$ drawn from an unknown distribution $p(x)$, this paper is interested in drawing new samples independent of $x_i$ but coming from the same unknown distribution p(x). \n\nThe classical approach consists of learning $p(x)$ from $x_i$ (i.e., approximated by $\\tilde{p}(x)$) and sampling new points according to the approximated version $\\tilde{p}(x)$, which remains a difficult task given sometimes the highly non-convex character of $p(x)$. \n\nFor these reasons, the authors propose to noise the data $x$ using $M$ different channels whose noise levels are chosen to obtain data denoted $y$. Since the level noise of each channel is known, this allows using a Bayesian estimator $\\hat x(y)$ to find the samples in the original space. The advantage of this approach is that the Bayesian estimator depends on the density $p(y)$, which allows, thanks to a chosen parameterization of $p(y)$ (which therefore leads to a parametrization of $\\hat{x}(y)$), to find the optimal parameters through a least-square objective by minimizing $\\|x-\\hat{x}(y)\\|^2$ on the training data. The Bayesian estimator thus allows generating new samples using the optimized density $p^{\\star}(y)$ with the corresponding optimal parameters. \n\nThe authors make a connection between the proposed method and some methods of the literature constituting new intuitions, different views of algorithms proposed in the literature. The experiments conducted suggest the efficiency of the algorithm to generate the most diverse samples", "summary_of_the_review": "The paper proposes an original work by offering a way out of the problem of approximating very complex densities by passing through a smoother variety. The method can be inspiring since it makes some connection and gives another point of view of the Denoising AutoEncoder proposed in the literature. The experiment shows that the method allows for generating stable and mixed samples and does not get stuck to some mode. Even if some discussions are not well-motivated, I recommend acceptance for the paper.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "1: The main claims of the paper are incorrect or not at all supported by theory or empirical results.", "technical_novelty_and_significance": "1: The contributions are neither significant nor novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies the problem of generative modeling by convolving an unknown complex density with a factorial kernel called multi-measurement noise model (MNM) to obtain a smoother density that is easier to sample from. Poisson and Gaussian MNMs are proposed for the convolution. Experiment regarding image synthesis are conducted to demonstrate the effectiveness of the proposed framework. The paper studies a problem that is of great interest to the machine learning community, and the results are impressive and promising. However, the paper in the current form lack a comparative discussion and a quantitative comparison with some related works. After the rebuttal, all three reviewers tend to accept the paper. After several rounds of internal discussion among AC, reviewers and authors, the AC agrees with the reviewers, and recommends accepting the paper, given the changes the authors promised to make.  

In summary, the AC recommends an acceptance and urges the authors to further revise their paper by adding a comparative discussion with those closely related works regarding generative models using MCMC sampling.
Decision: Accept (Poster)