Unbiased Contrastive Divergence Algorithm for Training Energy-Based Latent Variable Models
Authors: Yixuan Qiu
Lingsong Zhang
Xiao Wang
Venue:  (2020)
Keywords: energy model
restricted Boltzmann machine
contrastive divergence
unbiased Markov chain Monte Carlo
distribution coupling

Abstract: The contrastive divergence algorithm is a popular approach to training energy-based latent variable models, which has been widely used in many machine learning models such as the restricted Boltzmann machines and deep belief nets. Despite its empirical success, the contrastive divergence algorithm is also known to have biases that severely affect its convergence. In this article we propose an unbiased version of the contrastive divergence algorithm that completely removes its bias in stochastic gradient methods, based on recent advances on unbiased Markov chain Monte Carlo methods. Rigorous theoretical analysis is developed to justify the proposed algorithm, and numerical experiments show that it significantly improves the existing method. Our findings suggest that the unbiased contrastive divergence algorithm is a promising approach to training general energy-based latent variable models.

Published as a conference paper at ICLR 2020
UNBIASED CONTRASTIVE DIVERGENCE ALGORITHM FOR TRAINING ENERGY-BASED LATENT VARIABLE MODELS
Yixuan Qiu
Lingsong Zhang Xiao Wang
Department of Statisties and Data Science Cameeii Mellon University Pitsburgh, PA 15213, USA y xuang@andrew cmu.cmc.ec......o eduu
Department of Statistics Purdue University West Lafayette, IN 47907, USA {1 lingsong, wangxiao) (mmuxiie))purduu..ee
ABSTRACT
The contrastive divergence algorithm is popular approach to training energy- based latent variable models, which has been widely used in many machine learn- ing models such the restricted Boltzmann machines and deep belief nets. De- spite its empirical success, the contrastive divergence algorithm also known to have biases that severely affect its convergence. In this article we propose biased version the contrastive divergence algorithm that completely removes its bias in stochastic gradient methods, based recent advances on unbiased Markov chain Monte Carlo methods. Rigorous theoretical analysis is developed to justify the proposed algorithm, and numerical experiments show that significantly im- proves the existing method. Our findings suggest that the unbiased contrastive divergence algorithm is a promising approach to training general energy- based latent variable models.
INTRODUCTION
Energy-based latent variable models cover broad class of generative models that are frequently used to characterize sophisticated distributions of high-ddmensional data. Popular examples this class include the restricted Boltzmann machines (RBM, (Smolensky 1986 Hinton 2012). deep be- lief nets (Hinton et al. 2006), and exponential family harmoniums (Welling et al. (2055). among many others. Energy-basee based models are complementary directed generative models such as the variational autoencoders (Kingma & Welling 2014), and have gained great success in synthesizing realistic data samples (Xie al., 20188) They can also be combined with directed models to build more sophisticated structures (Xie et al. 2018a). In this article we focus on the energy-based latent variable model, whose general form can be expressed in terms of the joint distribution of visible random vector, v E V C RP, and hidden or latent random vector, h E H C RR, with the density function
1
p(u, h: 0)
Z(0) exp( (xp(-F(1,.) h; 0)).
(1)
where 0 θ is the unknown parameter vector, Elv, h;; 0) is the energy function, and Z()) is normalizing constant to ensure that pvv, h; 0) is a legitimate probability density or mass function. The model distribution, pv 0), is defined to be the marginal distribution of p(v. h: 0).
Similar to many other machine learning models, the standard approach to estimating the parameter vector 0 is the maximum likelihood method. It can shown that the derivative of the log- logllikeldoood function can be expressed the difference of two expectations, and hence Monte Carlo methods, especially the Markov chain Monte Carlo (MCMC, Gilks et al., 1995), can be used to approximate the gradient. Various optimization techniques, such the stochastic gradient method (SG, Robbins & Monro, 1951 Bottou| 2010), can then proceed to iteratively update the parameter estimate. This strategy, though elegant in theory, is not without limitations. In particular, MCMC estimators typically consistent the limiting case, but biased on finite steps, so one needs to run MCMC for a long time to obtain accurate gradient, which would take tremendous amount of computing time.

---

Published as a conference paper at ICLR 2020
To reduce the computational complexity, HHinton 2002) proposed simple and fast algorithm, called the contrastive divergence (CD) algorithm. The basic idea of CD is to truncate MCMC at the k-th step, and use the resulting approximate gradient to update 9. where k is fixed integer small as one. Such an approach is usually referred to as the CD-k algorithm. The simplicity and computa- tional efficiency of CD makes it widely used in many popular energy-based models, and there was also numerous empirical evidence to illustrate the effectiveness of CD. More recently, the CD idea applied to directed generative models (Ruiz & Titsias 2019), where CD is used to define loss function that combines variational inference and MCMC.
However, the success of CD also raised a lot of questions regarding its convergence properties. Both theoretical and empirical results show that CD in general does not converge to a local minimum of the likelihood function (Carreira- Perpiinn models (Schulz et al. [ Fischer & Igel 2010). The main issue of CD is that the truncation of MCMC produces biased stochastic gradient for the tog-likelihood function in every iteration, and such uncontrolled biases may be accumulated to distort the true ascent direction. Due to this reason, the training of energy-based models has been a longstanding challenge in machine learning research.
Hinton 2005), and diverges even in some simple
In this article, we propose new unbiased contrastive divergence (UCD) algorithm based on recent advances in unbiased MCMC theory. which offers new possibilities for solving the model training unbiased estimator for the expectation with respect to the invariant distribution of a Markov chain. More recently, this estimator further extended to the MCMC setting by Jacob et al. (2017), using technique called coupling. At a high level, by carefully designing the MCMC algorithm, one is able to get an unbiased MCMC estimator with only finite number of Markov transitions.
problem. In the seminal work Glynn & Rhee (2014). the authors developed
Under the framework of Glynn & Rhee (2014) and andJJacob et al. (0117). the proposed UCD algorithm a saabbs-smlllr--bbseeeee training method for energy-basse latent variable models. We prove that the stochastic gradient generated by UCD is unbiased with finite variance, which implies the conver- gence of SG based on it. Similar to CD-k, UCD generates Markov chains to compute the gradient, but the chain in UCD stops at a random time, instead of a fixed one as in CD-k. The theoretical analysis indicates that the stopping time has a finite expectation, so on average the computation can be completed in finite time. Besides theoretical justifications, our numerical experiments show that UCD significantly improves existing training algorithms, suggesting that it is a promising ap- proach with solid convergence guarantee. The implementation of the UCD algorithm is available at pitttps:/// /t/ptittuuubb.t com/ yixuann///// cditau The highlights of this article are as follows:
We develop a new training algorithm for general mvergy-based latent variable models that include many popular models (e.g. RBM) special cases. To our best knowledge, this is the first algorithm that has a solid convergence guarantee for such models.
The proposed algorithm resolves a longstanding problem of the CD algorithm, the bias in approximating the gradient. In particular, our method is completely unbiased, and theoret- ical justifications are developed to guarantee its convergence.
We have tailored specialized algorithm for RBM, which is shown to significantly reduce the computational cost.
2 A BRIEF REVIEW OF CONTRASTIVE DIVERGENCE
In this section briefly review the CD algorithm, and point out some of its weaknesses that have been studied in the existing literature. For a single observation v, the marginal data log-likelihood function is f(0;0) v) log [p (v; 0)} log[_ p(v, h; d) Assume that E(v, h; O) is contin- uously differentiable for 0, and then with n data points D = (Vi,..., V.)). ) the derivative of the log- likellhoooo function ((0: = Σ E(O: v:). also known as the score function, can be written as
() -n (Dv.)i-p(()))(Aeoo @dE(v, h; 0) 00
(vv.h))p(v....)) {
@o h: 0)
(2)
00
d0
where p(D) stands for the empirical distribution of D, and 2 0) is the conditional distribution of the latent variable h given v = v. A simple derivation of (2) can be found in Fischer & Igel (0014).
Throughout this article denote (v, h) X V x H and f(;; 0) OF.(v, h; 0)/30. Then
2

---

Published as a conference paper at ICLR 2020
the two expectations in (2) can be abbreviated as Ex(f(xx)) 0)} and EMif(x: 0)). respectively, where M = p(v, V, h:; 0) is the complete model distribution.
In many cases, for example the RBM model, ED 0)} has closed form, so the major com- putational difficulty comes from the M (x: 0)} term. A common scheme to approximate this expectation to then under mild conditions we have o ((((();)) = [x 9)). Of course, such a limit cannot be reached in finite steps, the CD-A algorithm truncates the Markov chain the k-th step, resulting in the following approximation:
Markov chain Eu {1
with M as the invariant distribution, and
(( = [ 0)} - (Ex: 0)
is easy to see that ( = n n- (e /e0 a stochastic approximation to the true gradient, so one can SG to update f via the iteration 0,+1 0, ar (0:), where Og is the parameter estimate in the i-th iteration, and as is the step size.
Despite its simplicity, various research articles have pointed out the weaknesses of the CD-k algo- rithm. For instance. /Sutskever& & Tieleman (2010) gave an example to show that (3((()) [0)} is not the gradient of any objective function, and Schulz et al. (2010): Fischer Igel (2014) studied numerical experiments in which CD-& does not converge at all for small A values. Carreira-Perpihann & Hin- ton (2005) considered the fixed points of A A()), the @ values such that E [00} 0, and showed that they do not match the fixed points of ( ; D) /a0 in general. This implies that even if CD-k converges, the resulting parameter estimate may not be a local minimum of the likelihood function.
Another variant of CD is the persistent contrastive divergence (PCD, Tieleman 2008 Tieleman & Hinton 2009), which has been reported to improve CD in many numerical experiments. However, it is still approximation method, and its convergence properties are more difficult to analyze, as the stochastic gradients generated by PCD become correlated across iterations. In fact, Schulz et al. (2010): Fischer & Igel (2010) also gave examples in which PCD failed to converge. There some other training methods as extensions to CD, such as the multi-grie method (Gao et al. 2018) and the short-run MCMC (Nijkamp et al. 2019), but all these methods inherit the bias of CD.
also
To summarize, it is surprising that virtually none of the popular training methods for energy-based models, including CD and PCD, provide solid convergence guarantee. The major defects of CD) stem from the fact that (()) is a biased estimator for the true log-likelihood gradient, and SG may fail with uncontrolled bias accumulation. To this end, the ultimate solution is design training algorithm that completely removes the bias of CD.
3 THE UNBIASED CONTRASTIVE DIVERGENCE ALGORITHM
3.1 UNBIASED MCMC ESTIMATORS
Since CD highly relies the MCMC method, the main ingredient of the proposed UCD algorithm is the theory of unbiased MCMC developed by/Glnnn& & Rhee (0144) Jacob et al. (2017). Consider the second term in 2 namely, [ 0)}. In what follows we omit the dependence on 0 for brevity if no confusion is caused. If a Markov chain { satisfies [E((())) + EM [x() as I + oo, then under some regularity conditions, we express the limit as telescoping sum,
Em(((x) B((
Σ [B((())) -E((-))))
for any fixed & 2 O. Now assume that there exists another Markov chain {ni such that & and ne have the same marginal distributions for all t 2 0, and & 7パ-1 for all 2 where is some random time. If we allow the exchange of expectation and summation, then we would get
Σ (f(E)) - f(n. - = =E f(EN) + Σ (f(E) - f(n-1)1
Emif(x)) = E f(Ex) +
(=#+1
1=
where the first identity holds since [[( E(s for all 1.111200 0, and the second one due to the fact that & = 7/-1 fort 2 r. As consequeee, the quantity f(Ex)+ ET= (f(()) (f(E))-((==1)) an unbiased estimator for [x)] Such an idea seems rather simple, but the construction of the chain (n which we describe in the next section, is a highly non-trivial task.
3

---

Published as a conference paper at ICLR 2020
3.2 COUPLING OF MARKOV CHAINS
Let Mr denote the marginal distribution of Markov chain {Ec} at the 1-th step. By construction, Mr converges to M as t + To develop the unbiased estimator H2 (C, the second chain {ma must satisfy two conditions: (1) marginally f ~ Me; (2) { and the lag-one sequence {7)-1} will meet and stay identical after some random time Condition (1) can be trivially met if (5x) and (n: are sampled independently. However, in this way the probability that & = 7/--1 may be extremely small, or even be zero for continuous random variables. Therefore, a special joint distribution for (Ee, 17-1) needs to be assigned subject to & MM, and T-- ~ M Such a pair of random variables under the marginal distribution constraints is called coupling. and for our purpose we attempt to seek a coupling scheme such that P(s = 77-1) > O. Figure | illustrates the coupling process of two Markov chains {Er} and (7t-1))
t 170
Figure 1: An illustration of the coupling process. {&r} and {me} start from the same value, and have the same marginal distribution MM, at each step. The two chains are correlated in such way that the event & = 1/t-1 occurs with a posi- tive probability for each t. After ran- dom time r (r = 5 in the illustration), { meets }m-13 and they stay identi- cal afterwards.
77s E
176 ε
174
E
M1
M2
M3
M4
177 Es
To implement such a coupling, first let }E+} and {nr. } start from the same initial value Eo = 70- and additionally draw &1 ~ T-l,,)) where T(y (z)) stands for the transition density function from state z to state y. Next, we need to draw (E2, 77.) such that marginally &n ~ -M2 and 7/2 ~ MM. which can be difficult task as Mr may not have closed forms. Fortunately, it is much simplified for Markov chains: due to the Markov property. &2 and 7/ will have the requested marginal distributions if we sample E2 T(1(s) and 7/1 M% -((110) conditional on & and 7/-- That is, the coupling of Markov chains can be achieved by the coupling of one-step transitions, which is a much simpler task. Define two density functions p(-) = =(((k) and g(-) (-lmo)). and then the problem reduces to drawing coupling (5.1)) such that & ~ p(-),
g(). and P(g n) > 0, which can be accomplished via the maximal coupling technique (Appendix(A.I)
Specific to our problem ( we need to sample x (v. h) from p(v, h; 0). In energy- based latent variable models, the most widely-used MCMC method the Gibbs sampler (Geman & Geman 1984), which sequentially updates one block of x based on the conditional distribution of this block given the rest. As an example, in RBM models v = h} and h|[v = 4} follow multivariate Bernoulli distributions with independent components, which coupling for Gibbs samplers was briefly mentioned in [Jacob et al. (2017) as a special case of the Metropolis- Hearopolis-Hastings scheme (Metropolis et al. 1953, Hastings| 1970), but next we show that some specific structure of Gibbs samplers can be utilized to simplity the process.
very easy sample from. The
For simplicity and clarity, we assume that the Gibbs sampler for M follow the natural division of blocks x = (v, h). That is, one can easily sample from the two transition distributions (.( == p(vl h: 0) and (h)b) v) = p(hles lu; 0). The more sophisticated cases, for example h consists of multiple layers h = (h)..... hy). can be dealt with similarly. In Algorithm I we describe the steps to sample two coupled chains {E. (v,, ht)) and (n (o) h() based on the Gibbs sampler.
Three remarks are made for Algorithm D (1) The meeting event (line 4) only depends on the T, transition density. To verify this, note that at the t-th step, we need draw εelE--1 ~ T(-IE(-1) and
--1/1--2 ~ T--1
where T( h v. h) = T, (6) h Tr (h] () is the transition density for a full update cycle. It is easy to show that T(ö, ,/o'''')/TTT" /T' hlj h) =
so the T
part cancels in the ratio. (2) Once & and 7/-1 meet, they stay identical afterwards, because by then T. T. = To 1).. and the event in line 4 always happens. (3) Line 7 is a rejection sampling step. In our numerical experiments we find that very few samples rejected, its cost is tiny.
4

---

Published as a conference paper at ICLR 2020
Algorithm 1 Coupling method for the Gibbs sampler
Input: Densities T.(v ft) and Ti (hlo),, initial values Eo (vo, ho) - 110 (v/s hlg), Tmax
Output: Coupled chains {Er} and {ne
1: Sample U1 ~ To(-1ho) and h1 ~ (-(-))1). Set {1 = (v1.h1)
2: for t= 2.3.... do
3: Sample 02 ~ TTo[-hh 1), h ~ (ee) and U ~ Uniform(O,) 1)
4: if U To (v,( (To h1-1)) or T max (maximum stopping time) then
5: Sc = (Ur, ht), 11-11
6: else
7:
Sample 1 To(- h²2 TT(()//11). and U Uniffor(O,))
until U (o ( To//TT((// (1/p(-/1)///))
8: Set & = (Ur: hr), M1- (1/-1; hi-1)
9: end if'
end for
3.3 UNBIASED CONTRASTIVE DIVERGENCE
The technical tools introduced in Sections 3.1 and 3.2 enable us to develop new algorithm to train model ( Recall that the true gradient of the bog-likelihood function is given by (2) The first term, Ef })}, can be computed exactly, and the second term, E (x; 0) is approximated by an unbiased estimator (0) == (Exx + Ei=k+
(n--1),, where the coupled Markov chains {&+} and {me) are generated by Algorithm Assume that the parameter vector 0 lies in closed convex set θ, and let (e(.) denote the projection onto e. Putting the pieces together, Algorithm 2 illustrates the UCD algorithm for training energy-based latent variable models. The initial chain length & can be any fixed number, and in this article we take k = 1 for all the numerical experiments.
Algorithm UCD Algorithm for estimating
Input: T, {.x). k, initial value Og Output: Parameter estimate for 0
1: for = 0, ,,...,,T' 1 do
2: Draw one data point v ~ p(D). and sample h ~ o(h/u;;0 0,)
3: Set &o 7o (v, h), and run Algorithm with 0 = 0, until ET. 17--1
4: g(0)
[h) f.f(v. h; 0)} f(E)) + Ei=4 f(( 5: 0++1 + Pe (0, + a; ((.))
6: end for
θ= T Σ T 0:
7: return
Next, we analyze the theoretical property of Algorithmmaan state the conditions for it to converge. As a standard setting, we assume that the Markov chains generated by the Gibbs sampler are g-- irreducible and aperiodic (Meyn & Tweedie 2012). This is a very mild assumption that every practical Gibbs sampler should satisfy. Then we make the following two assumptions that guarantee the convergence of Gibbs samplers.
Assumption 1. (Drift condition) There exist a pair functions T [1, +00), / H (00) and constants 71, 72, L1, L2 > O such that 71112 < and
ev=p(vlh:o)) S (h) + L1: Enmp(h)e:0) S (2)(v)) + L2, E E H. 0 E
Also, there exist constants e > and D > such that If (x; 0) < I(h) and En (h. ((h) S D for all x = (v, h) E X and E e.
Assumption 2. (Minorization condition) There exist constants d > 2(12L1 + L2)((1 712), E > 0. and density function q(-) such that p(u 0) 2 =q(v) for all h E D. v E V, and 0 E θ, where DD (h H. 1(h) < d)
In the following theorem show three important facts about the proposed stochastic gradient ğ(0):
(1) g(0) is unbiased for the true score function; (2) it has a bounded second moment uniformly in 0; (3) in expectation it can be computed in finite time.
5

---

Published as a conference paper at ICLR 2020
Theorem 1. Under Assumptions ! and there exist constants D1, D2 0 such that [{ğ(0)}
88(0; v)/00, E [()((0))) < D1, and E (i) I D2 for all O E 2 and i 1,2, 2....... 1
Theorem ! provides the building blocks for the convergence analysis of Algorithm 2 With the unbiased gradient estimator and the bounded second moment, we establish a solid convergence guarantee for the proposed algorithm. As typical setting, in the following corollary we consider convex log-lik log-likelihood function.
Corollary 1. Assume that 8(0;0) is convex and L-Lipschitz continuous in E θ and θ closed and bounded convex set. Then by choosing a = 00/ vinv/iffo fir some constant au > 0, we have e e(0: v) < O(1/VT), where is the maximum value of v).
The proof Corollary [ is standard, see for example Bottou (010);; Bottou et al. 20018) When the log-likelihood function is nonconvex in the RBM model, there also other versions of the convergence result for SG, for example Theorem 4.10 of Bottou et al. (2018). Such directions can be studied separately and are omitted here.
Finally, we shall point out an important special case of Theorem [ i.e., if the Markov chain {&t has finite states, then the two assumptions are automatically satisfied. This shows that many widely-used models, for example RBM, can directly use the UCD algorithm without the need to find such r(-) and I(.) functions. We summarize this useful fact in the following corollary.
Corollary 2. If X is a finite state space and θ compact, then Assumptions and 2 hold, and Theorem appliees.
TRAINING RESTRICTED BOLTZMANN MACHINES
RBM is one of the most popular and widely-used energy models in machine learning, defined by the energy function E(v, h; 0) = Tb Wh h hTc where e E {0,1)" E {)" and 0 (W, b, are model parameters. The Gibbs sampler for RBM has nice structure: let (xx) = 1/(1 + (x)) be the sigmoid function, and then v = h} ~ Bernoulli (Wh b)) and h| {v = /} ~ Bernoulli(a (W"T + c)). The coupling method in Algorithm | directly works for RBM, but here we show improved version that is tailored for RBM and more efficient.
Let u, p E RT, and the notation and y otherwise. Also let (o h) I p. p.)' denote the transition density from v. where p (pI..... Pm. T a h b). Then the specialized coupling method for RBM is given in Algorithm
1(u S p} stands for binary vector such that y; if M S Ps
Algorithm 3 Coupling method for RBM
Input: Model parameters W, b, c, step-/ states & = (vr, Ie). 1/-11 (b/-11)(/-1)
Output: New states E = (V++1, 1441), (v/, h/))
Sample U1- Uniform 1). Z1 Uniform 1]m), and set 1141 1(2 s (WW, b)}
2: if U1 - (v (). (h) then
3: Set v/ U++1
4: else
5: repeat
6: Sample U2 Uniform(0, 1). U3- Uniform(0. 1), Z₂ Uniform([0, 1]
7: Vt+1 has not been accepted then
8: Propose U++1 1[2, S oWh,, b)), accept if U2 Tulluu:llbi
9: end if
if vf has not been accepted then
11: Propose 0y 1 < (Wh, h + b)), accept if TT(oo//o))//TTT00/...
end if
13: until vi+1 and v/ are both accepted
end if
15:
Sample Z3 ~ Uniform((0, 1|*)
16: Set h4+1 1[2 o WW
U1+i + c)). h1 1[2s s so(W"ot

---

Published as a conference paper at ICLR 2020
The intuition behind Algorithm is is the following: line 2 indicates that is also maximal coupling method, so the probability P(E++1 = m1) is the same as Algorithm[[] However, in the event {Eo+1 o Te) and 1/. are independent in A common random variates 22 and Z3- The correlation between Ê1+1 and 1/: helps to make P(E++2 = n--1) larger, thus accelerating the meeting of {Ee} and (n--1) A more rigorous justification of this algorithm is given in Appendix A..
correlated in Algorithm achieved by the use of
bo)"", aaWW + b), h) for RBM, corresponding to the parameters 0 = (W, b, c). With the coupled chains (fr and {n-1)-
Finally, it is known that in the gradient expression (2) [x; 0) =
RBM can then be trained using UCD given by Algorithm
5 RELATED WORK
In this section we highlight the novelty of our article and clarify its overlap with prior art. In literature there were several attempts to prove the convergence of CD in special cases, or to reduce the bias of CD using other sampling techniques, all with undesirable results. For example, Yuille (2005) gave conditions for CD to converge, which unfortunately can hardly be satisfied in any realistic models. Jiang et al. (2018) showed convergence result of CD for the exponential families, but consequently the model is restrictive and does not include the latent variable model. Krause et al. (2018) used importance sampling estimate the normalizing constant, which is consistent with large sample. However, it still induces bias in the finite case, and the bias heavily depends on the choice of the importance weights. In contrast, the UCD method proposed in this article directly fixes the bias of CD, and hence bypasses the challenges in algorithm convergence.
Unbiased MCMC is a relatively new topic in statisties and machine learning. Some background knowledge in this article, for example Section 3.1 is taken from Jacob et al. (2017), which estab- lished general framework for unbiased MCMC. Our new contributions are in the following aspects. First, we have developed Algorithm and [hecorem[Jxxclusvvel) for the Gibbs sampler, taking into account the special structure of Gibbs MCMC. Second, our theoretical results, including Theorem 1 aand Corollary 2 have more practical assumptions than the ones in Jacob al. (2017). For exam- ple, one of their key assumptions is that E ((+)) is uniformly bounded for every finite step, which is quite abstract and hard to verify in practice compared with our Assumption [ Third, [aiob eet al. (2017) studied MCMC with a fixed target distribution, whereas we need to control the variance of estimators that evolve with parameter updates. Finally, in Section we develop a specialized coupling algorithm for RBM, which is shown to be more efficient than the general one.
Another related work is the Markov chain Las Vegas method" (MCLV, Savarese et al. 2018 which also constructs an unbiased estimator for the score function of RBM. The main differences bewteen our method and MCLV are as follows: (1) MCLV based on the regeneration theory of Markov chains, whereas UCD is built upon the coupling technique. (2) MCLV is exclusively designed for RBM, but UCD applies to a broader range of models. (3) In its current state, MCLV cannot handle continuous random variables, but UCD can. More comprehensive comparisons between MCLV and UCD are left for future exploration.
6 NUMERICAL EXPERIMENTS
6.1 BARS-AND-STRIPER DATA
We compare CD-k, PCD, and the proposed UCD algorithm for training RBM models on different data sets. In the first experiment we reproduce the results for the bars-and-stripes (BAS) data that have been studied by [Schulz et al. (2010); Fischer & Igel 2010 2014). It is a small data set with 36 data points and 16 binary variables, and is fit by a small model with 16 hidden units. However, it is one of the most important benchmark data sets for RBM since its log- log-likelihood value can be evaluated exactly, and demonstrates the divergence of CD-based training algorithms. In our study, k is set to 1 for CD (more experiments with larger & are given in Appendix BB..). and each algorithm run for 100 times, accounting for the randomness in the training process. common learning
We were informed of this method by public comment on httpss:/ ew www.net/forum? rlyyceSYPr after the current paper was accepted.
7

---

Published as a conference paper at ICLR 2020
rate 0.01 is and 1000 parallel Markov chains are used to approximate the gradient in each iteration. The results are shown in Figure
Log-lkellheedd Function Value
2.0
-150-
Averaigo Stopping Time
Rejected Samples per
2.5-
200
-250
Methed
-300
FC3
-350
0.0
5000
7500
10000
10000
Iteration
Iteration
Iteration
Figure 2: Left: exact log-likelihood values in each iteration. The shaded bands stand for the 2.5% and 97.5% quantiles across 100 runs, and the three trajectories in darker colors are sample learning curves in one run. Middle: average stopping time for UCD in each iteration. Right: average number of rejected samples in the coupling algorithm for UCD.
Figure shows the following findings. First, reproduce the results in Fischer & Igel (2014) that CD and PCD fail to converge to the true maximum likelihood value. In contrast, UCD does converge. Second, UCD has adaptive choice of the stopping time in the Markov chain, compared to the fixed in CD. In the BAS data, the stopping time has a steep increase around the 1200th iteration. Interestingly, this is exactly where CD begins fail. An interpretation this phenomenon is that UCD automatically uses large MCMC sample for parameter values that result in a "hard" distribution. An even more surprising fact is that the average stopping time for UCD 2.40, making it computationally more efficient than the CD-20 algorithm, where 20 is the smallest k such that CD-k training is comparable to UCD (see Appendix [..]] for more discussions). Third, the cost of the rejection sampling step in UCD (line 7 of Algorithm [ is tiny, as the number of rejected samples rarely goes above two. Finally, UCD does not see massive increase in the variance. In fact, at the end of training the quantile band for UCD is much narrower than those of CD and PCD. All these findings further highlight the advantages of UCD.
6.2 SIMULATED RBM DATA
In the second example show that the findings for the BAS data be observed in other model settings. We simulate data set from an RBM model with 200 visible units and 20 hidden units, where the entries of weight and bias parameters sample size of the simulated data set is 1000, and we fit an RBM model using 100 hidden units, which is larger than the model since we intend to mimic the common practice of overparame- terization in RBM training. We use a common learning rate = 0.2 and 1000 Markov chains in each iteration for all three algorithms. The log-likelihood values are approximated by Monte Carlo averages. The result is given in Figure 3 which shows similar patterns to the BAS data: CD and PCD eventually diverge, whereas UCD follows the typical behavior of SG. Additional experiments and discussions given in Appendix B.2
all generated from A 1) distribution. The
CD-1
PCD
Unblased CD
Log Log-liken nigliiklihood Value Value
-60
-65
-75
-80
500
1000
1500
2000
0
500
1000
1500
2000
500
1000
1500
2000
Iteration
Figure 3: Approximate log-likelihood values for each algorithm on the simulated RBM data set.
8

---

Published as a conference paper at ICLR 2020
6.3 AASHION-MNISS DATA
Next we consider the Fashion-MNIST data se 2 replacement for the well-known but overused MNIST data set of handwritten digits (LeCun et al. 1990). Each data point in Fashion-MNIST contains 784 values within [0, 1]. representing a 28 28 greyscale image. The whole data set contains 60000 images, and we binarize the data by treating original values as probabilities and sampling from a Bernoulli distribution for each element. On the binarized data, we fit an RBM with 1000 hidden units, and train the model with different algorithms using mini-batch size of 1000 and a learning rate 0.1. For each training algorithm, 1000 parallel Markov chains are used to compute the gradient. Figure demonstrates the training trajectories of the three algorithms.
CD-1
PCD
Unbiased CD
Log-lik log-likelhoood Value
300
-400
-500
-600
1000
2000
3000 0
1000
2000
3000 0
1000
2000
3000
Iteration
Figure 4: Approximate log-likelihood values for each algorithm on the MMSIST data set.
The training patterns for CD and PCD are surprising: the lig-liiel values of CD are decreasing, and PCD seems to bounce back and forth among three different paths. fact, Appendix [B.3] we show that CD-k does not converge even if & is as large 30, whereas the average stopping time for UCD is about 27, These results further demonstrate the superior performance of UCD.
7 DISCUSSION
In this article we use the unbiased MCMC technique to estimate the score function of energy-based latent variable models, which effectively fixes the bias of CD algorithms. Both the theoretical analy- sis and the numerical experiments at different scales demonstrate that the impact of bias elimination huge.
is expected that UCD may have slightly larger variance compared with CD and PCD, but we emphasize that the value of UCD is not simple question of bias and variance trade-off. This is be- cause for MCMC-based methods, the variance can always be reduced by running more independent Markov chains and taking the average, while removing the bias is highly non-trivial. Moreover, in the context of SG methods, bias is typically more harmful than variance, as the former may lead to divergence of the algorithm.
In terms of computational cost, it is true that UCD is in general slower than CD-1 or PCD, but it is comparable to or even faster than CD-k that achieves similar log- nog-likelhhod level. Therefore, UCD is not meant to completely replace CD or PCD, but rather to serve important addition to the existing training algorithms. In practice, a very useful technique to first run the fast CD-1 or PCD to the near-optimum, and then proceed with UCD for guaranteed convergence.
REFERENCES
Léon Bottou. Large-scale machine learning with stochastic gradient descent. In Proceedings of COMPSTAT'2010, pp. 177-186. Springer, 2010. 3.3
Léon Bottou, Frank E Curtis, and Jorge Nocedal. Optimization methods for large-scale machine learning. SIAM Review, 60 2018. 3.3
httpss://essarran za www.o.cooo/......eee come/m Fashion-mme

---

Published as a conference paper at ICLR 2020
Miguel Á. ceedings of the Tenth International Workshop on Artificial Intelligence and Statistics, 2005.
and Geoffrey E. Hinton. On contrastive divergence learning. In Pro-
Asja Fischer and Christian Igel. Empirical analysis of the divergence of gibbs sampling based learning algorithms for restricted boltzmann machines, In International Conference Artificial Neural Networks, pp. 208-217. Springer, 2010.
6.1
Asja Fischer and Christian Igel. Training restricted boltzmann machines: An introduction. Pattern Recognition, 47(1):25 39, 2014.
2 6.1 6.1
Ruiqi Gao, Yang Lu, Junpei Zhou, Song-Chun Zhu, and Ying Nian Wu. Learning generative con- vnets via multi multi-grid modeling and sampling. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 9155-9164, 2018. 2
Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6 (6):721-741, 1984. 3.2
Walter R Gilks, Sylvia Richardson, and David Spiegelhalter. Markov Chain Monte Carlo in Prac- fice. Chapman and Hall/CRC, 1995. !
Peter Glynn and Chang-han Rhee. Exact estimation for markov chain equilibrium expectations. Journal of Applied Probability, 51(A):377-389, 2014. 3.1
W. Keith Hastings. Monte carlo sampling methods using markov chains and their applications. Biometrika, 57(1):97-109, 1970. 3.2
Geoffrey E Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):1771-1800, 2002.
Geoffrey E Hinton. A practical guide training restricted boltzmann machines. In Neural networks: Tricks of the trade, pp. 599-619. Springer, 2012. [
Geoffrey E Hinton, Simon Osindero, and Yee Whye Teh. A fast learning algorithm for deep belief' nets. Neural computation, 187751157- 2006.
Pierre E Jacob, John O'Leary, and Yves F Atchadé. Unbiased markov chain monte carlo with
couplings. arXiv preprint 1708.036222 2017. [ 3.1 3.2 D.I D.I D.I
Bai Jiang, Tung- Yu Wu, Yifan Jin, Wing H Wong, et al. Convergence of contrastive divergence algorithm in exponential family. The Annals of Statistics, 46(6A):306733098 2018.
Alicia A Johnson and Owen Burbank. Geometric ergodicity and scanning strategies for two- component gibbs samplers. Communications in StatisticsThheryy and Methods, 44(15):3125- 3145, 2015. D.I
Diederik P Kingma and Max Welling. Stochastic gradient vb and the variational auto-encoded In Proceedings of the 2nd International Conference on Learning Representations, 2014. []]
Oswin Krause, Asja Fischer, and Christian Igel. Populaaion-conrratiieddddvengne Does consis- tency help with rbm training? Pattern Recognition Letters, 102:1-7, 2018.
Yann LeCun, Bernhard E Boser, John S Denker, Donnie Henderson, Richard E Howard, Wayne E Hubbard, and Lawrence Jackel. Handwritten digit recognition with a pack-propagation net- work. In Advances in Neural Information Processing Systems 2, pp. 396-404, 1990. 6.3
David A Levin and Yuval Peres. Markov Chains and Mixing Times, volume 107. American Mathe- matical Soc., 2017. A.I
Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 1086:100 1953. 3.2
10

---

Published as a conference paper at ICLR 2020
Sean P Meyn and Richard L Tweedie. Markov Chains and Stochastic Stability, Springer Science & Business Media, 2012. 3.3
Erik Nijkamp, Mitch Hill, Song-Chun Zhu, and Ying Nian Wu. Learning non-convergent non- persistent short-run meme toward enrgyy-based model. In Advances in Neural Information Pro- cessing Systems 32, pp. 0233-5243, 2019.
Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathemat- ical Statistics, (2(3):400-407 1951.
Jeffrey S Rosenthal. Minorization conditions and convergence rates for markov chain monte carlo. Journal of the American Statistical Association, 90(430):558-566 1995. D.1
Francisco JR Ruiz and Michalis K Titsias. contrastive divergence for combining variational infer- ence and meme. arXiv preprint arXhv:1905.04062. 2019.
Pedro HP Savarese, Mayank Kakodkar, and Bruno Ribeiro. From monte carlo to las vegas: Im- proving restricted boltzmann machine training through stopping sets. In Thirty-Second AAAI Conference on Artificial Intelligence, 2018.
Hannes Schulz, Andreas Müller, and Sven Behnke. Investigating convergence restricted boltz- mann machine learning. In NIPS 2010 Workshop on Deep Learning and Unsupervised Feature Learning, 2010.
2.6.6 6.1
Paul Smolensky. Information processing in dynamical systems: Foundations of harmony theory. Technical report, Colorado Univ at Boulder Dept of Computer Science, 1986.
Ilya Sutskever and Tijmen Tieleman. On the convergence properties of contrastive divergence. In Proceedings of the Thirteenth International Conference Artificial Intelligence and Statistics, pp. 789-795, 2010.
Tijmen Tieleman. Training restricted boltzmann machines using approximations to the likelihood gradient. In Proceedings of the 25th International Conference on Machine Learning, pp. 1064- 1071. ACM, 2008. 2
Tijmen Tieleman and Geoffrey Hinton. Using fast weights to improve persistent contrastive diver- gence. In Proceedings the ACM, 2009.
the 26th International Conference on Machine Learning, pp. 1033-1040.
Max Welling, Michal Rosen-Zvi, and Geoffrey E Hinton. Exponential family harmoniums with an application to information retrieval. In Advances in Neural Information Processing Systems 17, pp. 1481-1488, 2005.
Jianwen Xie, Yang Lu, Song-Chur Zhu, and Yingnian Wu. A theory of generative convnet. In Proceedings of the 33rd International Conference Machine Learning, pp. 2635-266-42 2016.
Jianwen Xie, Yang Lu, Ruiqi Gao, and Ying Nian Wu. Cooperative learning of energy-based model and latent variable model via meme teaching. In Thirty-Second AAAI Conference on Artificial Intelligence, 2018 8a.
Jianwen Xie. Zilong Zheng, Ruiqi Gao, Wenguan Wang, Song-Chun Zhu, and Ying Nian Wu. Learning descriptor networks for 3d shape synthesis and analysis. In Proceedings of the IEEE Conference on Computer Vision and Puttern Recognition, pp. 8629-8638, 10.
Alan L Yuille. The convergence of contrastive divergences. Advances Neural Information Processing Systems 17, pp. 1593-1600, 2005.

---

Published as a conference paper at ICLR 2020
A DISTRIBUTION COUPLING METHODS
A.1 THE MAXIMAL COUPLING ALGORITHM
Let (E) be pair of random variables defined on the same probability space, and p and q be two distributions. (E. 7) is called a coupling of p and q if marginally 5 ~ p and 1j ~ q. is well known (see e.g. Proposition 4.7 of Levin & Peres, 2017) that among all possible joint distributions of (E,
f min(p(a)) ((x))dr,
P(g = n) S 11 - a||Tv parv=/
where ||p g1 is the total variation distance between and 9. The maximal coupling algorithm, given A generates a coupling (5. 17) that achieves the bound (4). That is, it maximizes the probability that two random variables are equal subject to their marginal distributions.
Algorithm Maximal coupling of two distributions p(-) and q()), from Jacob al. (0117).
Input: Density functions p(-) and q(-)
Output: A coupling (g, 7)) with & ~ p(-) and v ~ g())
1: Sample ç ~~ p(-) and U ~ Uniform(0, 1) independently
2: if U S ((E)/(()) then
3: return (E.E)
4: else
5: Sample 17 ~ q(-) and U Uniform 1) independently until U p(n)/q())
6: return (E.7)
7: end if
A.2 COUPLING METHOD FOR RBM
In the general maximal coupling method (Algorithm 4h. let p(.) (-1ss) and g(:) (- be the transition densities of RBM, and then it generates the new states ε and nhe with the probability P(E++1 = T) &c, 7h-1) maximized, If the event {E++1 = mh does not happen, then Er+1 and 1/s are sampled independently.
However, if E++1 and 1/, are close to each other, then the next iteration the probability P(E+22 = M--|EE+1, 77.) would be large, which helps to shorten the stopping time. Therefore, are motivated to minimize some type of distance between Ei+1 and Mt. For RBM, we characterize by [[[[vv-11 v2 Ve+1 # vz) and h(4 h 11+1,, vy)). Since all these variables are binary vectors, the norm of difference is basically the number of unequal components, and the problem reduces to the maximal coupling of Bernoulli variables.
[[ll)/ee1, o(). Given the v variables, h and hy follow Bernoulli distributions elementwisely, with individual mean vectors denoted as μa and μ12, respectively. It can be shown that the maximal coupling of two Bernoulli variables is achieved by using the same random variate, so h[ hill2 U241, vt) is minimized by setting hu+1 1(Z < p13 and h; 1(2 S μ2), where Z ~ Uniform((0, 1)²). This leads to lines 15-16 of Algorithm
First consider E(
Next, conditional on {ve+1 A v(). the maximal coupling algorithm generates V(+1 and v/ with marginal distributions proportional to (v) min(pr p2(v)} and P2(v) min P1 (v), (2(v)), respectively, where (1(-) = (- (4)) and P2(-) To(-Ih, (1/-1h Algorithm 3 implements this by two rejection sampling steps (line 8 and line 11). If we use the common random variates Z2 ~ Uniform ([0, 11²) for the proposals, then V and v/ are maximally coupled in the event that they accepted in the same iteration. This property helps to reduce E[[ v'll2/ve41 v().
To demonstrate that the specialized Algorithm 3 improves the general Algorithm I we sample coupled Markov chains from an RBM model with m 500, n 100, and the elements of (W, b, c) are all generated from a (0, 0. distribution. For both algorithms, the elements of the initial state vo are sampled from independent (ernoulli(0.5) distributions, and generate coupled Markov chains until the stopping time 7 is reached or the chain length exceeds 1000. Figure shows the distributions of for both algorithms based on 1000 replications.
12

---

Published as a conference paper at ICLR 2020
General (Algorithm 1)
Specialized (Algorithm 3)
200
Frequency
150
100
50
0
500
750
20
50
Stopping Time t
Figure 5: Distribution of the stopping time 7 for the general (Algorithm and specialized (Algo- rithm B coupling methods.
very clear that the general algorithm leads to a surprisingly long stopping time, and even has long tail beyond 1000. In contrast, 65.4% of the stopping times in the specialized algorithm are smaller than equal 10. In this sense the improvement brought about by Algorithm 3 Bis huge.
B ADDITIONAL RESULTS FOR NUMERICAL EXPERIMENTS
All experiments in this article were run on an Intel Xcon @ Gold 6126 processor with 12 cores and 24 threads. CD and PCD algorithms used the OpenBLAS library? for parallel matrix computations, and UCD used OpenMP to generate Markov chains in parallel.
B. BAS DATA
For BAS data, we gradually increase the value of k in the CD-k algorithm, and plot their training trajectories for the log-likelihood values (Figure (). It can be seen that the smallest k to make the result comparable to UCD is about 20. However, the stopping time for UCD has an average of 2.40 across all the iterations, which means that it is more efficient than fixed-l CD algorithm with similar performance. The running time for each training algorithm, which given in Table also supports this claim.
Log likellhood Function Value
1.0-
Method
CD-2
5
CD-10
CD-20
-350
2500
5000
7500
10000
Iteration
Figure 6: CD-k algorithms for training the BAS data with different k.
B.2 SIMULATED RBM DATA
For the experiment in Section 6.2 Figure 7 shows the training processes of CD-k algorithms with larger k. Figure gives the average stopping time for UCD each iteration, and the number of
https://wwer httss://w..gnnblas...
https:
13

---

Published as a conference paper at ICLR 2020
Table 1: Running time (in seconds) for each training algorithm the BAS experiment.
CD-1 w/ log-likelihood values
CD-1
CD-2
CD-5
CD-10
CD-20
PCD
UCD
65.01
5.38
8.44
17.89
33.17
63.51
5.18
15.99
discarded samples in the rejection sampling step. Table illustrates the computational time for each algorithm. We find that the computational cost of CCD is slightly larger than but very close to CD-1 and PCD.
CD 2
CD
CD 10
Log-liken Log-likelihodd Value
-60
-65
-70
-80
500
1000
1500
2000
500
1000
1500
2000
500
1000
1500
2000
Iteration
Figure 7: Approximate log-likelihood values for each algorithm on the simulated RBM data set.
2.2
Rejected Sampples per Clainn
20-
Average Stapaing
25
13-
1.5
1.4
500
1000
1500
2000
500
1000
1500
2000
Iteration
Iteration
Figure 8: Left: average stopping time for UCD in training the model on simulated RBM data. Right: average number of rejected samples in the coupling algorithm for UCD.
Table Running time (in seconds) for each training algorithm on the experiment in Section 6.2
CD-1 w/ log- likelihod values
CD-1
CD-2
CD-5
CD-10
PCD
UCD
216.46
16.52
24.44
48.54
89.09
16.18
28.47
One useful technique to combine the computational efficiency CD/PCD and the convergence of UCD is to use CD/PCD in the early stage of the optimization, and then take the resulting parameter values as initial starts for the UCD algorithm. Figure 9 demonstrates this idea: in the first 500 iterations the model is trained using PCD, and then the parameter values are fine-tuned by UCD with a small number of iterations and a reduced learning rate.
B.3 FASHION-MNIST DATA
For the Fashion-MNIS" data, Figure 10 gives the training trajectories CD-10, CD-20, and CD-30 algorithms, but unfortunately none of them have convergent pattern. Figure shows stopping time and the number of discarded samples in the UCD algorithm. Again, the cost for rejection sampling be ignored, and the value of gets stable around 30 after 2000 iterations.
shows the average
14

---

Published as a conference paper at ICLR 2020
-55
-60
Log- lilellhhodd Value
Stage
-70
PCD
UCD
-75
-80
200
400
Iteration
Figure 9: Combine PCD and UCD for efficient and convergent training.
Finally in Table 3 we demonstrate the rumning time for each training algorithm. Although UCD takes more time than other methods, the additional cost seems to be the necessary price to achieve convergent result, as all other methods compared have a divergent log-likelihood.
CD-10
CD-20
CD-30
Log-likelihode lihood Value
-300
-400
-500
-600
1000
3000
1000
3000
1000
3000
Iteration
Figure 10: Approximate log- eegliiklihhod values for CD algorithms the Fashion-MNIST data set.
Rejected Samp es per Chein
1.75
Average Stopping
60-
1.50
10
100
Iteration
Iteration
Figure 1I: Left: average stopping time 7 for UCD in training the model on Fashion-MNIS data. Right: average number of rejected samples in the coupling algorithm for UCD.
C VARIANCE OF GRADIENT ESTIMATES
To compare the variances of stochastic gradients generated by different algorithms, in each opti- mization iteration estimate the following gradient variance defined by
v(((0) = [|
15

---

Published as a conference paper at ICLR 2020
Table 3: Running time (in minutes) for each training algorithm the Fashion-MNIST experiment.
CD-1 w/ log likelikeiiood values
CD-1
CD-10
CD-20
CD-30
PCD
UCD
65. 14
25.63
114.85
212.22
311.89
25.53
494.01
where f(0 is an estimator for the second term of (2) and - || is the Frobenius norm. Among the training methods considered in this article, f(0) is biased in CD and PCD, and unbiased in UCD. We consider the example Section 6.2 and estimate V(f(0) using the 1000 parallel Markov chains in each iteration. Figure 12 shows the estimated gradient variances for b, and W along the training process.
Unbrased
Estimated Gradient Variance Variance
30
20
4000
1000
1000
20000
1000
20000
1000
20000
1000
20000
1000
20000
1000
2000
teration
Figure 12: Estimated gradient variance in each optimization iteration for different methods.
It can be observed from Figure 12 that the gradient variance of UCD is generally comparable to CD and PCD, with only few extreme cases. This finding suggests that the eliminated bias does not significantly increase the variance, which is important for the convergence speed of SG.
D PROOF OF THEOREMS
D.1 THEOREM
Omitting the iteration index for brevity, we first prove that the stopping time 7 has a finite expecta-
tion. Let (=(a=))) = (v, h/w', hh) (o)(hh))))(heee 0) denote the transition density for full
update cycle. Under Assumptions [and] 2 Lemma 2 ofJnonnon Burbank (2015) shows that there
exists constant 7 E [1172, 1) such that
Ex- S >((h') + 7/21 L2
(5)
for all x E X. Therefore, the drift condition holds for the transition density (a)la))).
Next, since Algorithm a maximal coupling algorithm, we have (see for example Jacob et al, 2017
P(E.+1 = melE: x, mh- x') = min((o( To(y)='))dd)
By Assumption 2 p(ulh; 0) 2 =q(v) for all h IB, v E V, and 0 O, so
P(E+1 1e x, M1- x') >
sq pp 6: O)dy
on x. E D V. where y = (ö, h). Note that g h 0 is joint density function, so
(h(o) O)dy 1, and hence P(E++1 m x, 71- x²) for all h T and O E O.
16

---

Published as a conference paper at ICLR 2020
Then by Proposition 3.4 of Jacob et al. 2017), , there exist constants 11 > 0 and Pr E (0, such that for all t > 0, P(T > t) I kil(ho)pp pi. As a result, we obtain
Ki k1((h))
B(7) - Σ
P(T -
(6)
1 P
t=
The right hand side of (6) does not depend on the value of 0, B(r) is uniform in the iteration index é.
To make sure that g()) is well defined, few regularity conditions, formulated Assumption 2.1 oaloa et al. (017)) need to be verified. Under the drift condition (5) and Assumption 2 Theorem 12 of [Rosentaal [195) shows that {&t} is geometrically ergodic Markov chain. Therefore, there exist constants K2 > 0 and P2 E 10 such that
[Bo(E)) EmØl S kgl(ha)p2
for all Φ : 10(u, h)l < ((h). Therefore, there exists constant M 0 such that [E((E)] I
Emil + h2) < D + k2) x M < 80.
Since [=; 0) 2+c < I(h) by Assumption [ we get (((S(E)))))) S M for all & > 0. Moreover, 7 implies that {(f(Ee) - E... These two results verify Assumption 2.1 of Jacob et al. 2017))
and hence the unbiasedness of g(0) is true by design.
Finally, we need to show that the second moment of the stochastic gradient is bounded uniformly in
O. Let A =((x) and A. [(.) - f()-1) for k 1, and then g2((0) Li A.. From
the ergodicity property 7 we immediately get R( M. For &
E(|A, (++)= (()) - (n--1) s 2 1 B(f((( 2+c + If [f(n=11) 22 = 22+
Then by Hölder's 's inequality,
E [AA)) B( > +3) (2(l) 2/(ot2) (P(T > 1)3c/let2)
S 4M2/(c+22 {mil(bo)pi} (o(02)
Since 52 (0) t= Do we get
(0) - As 1 =E
Σ 1=k+1
E/ > >/)((()))((>>>> > s]]]--
1= k=k+1
We also have
E |/ t)1(1 > s})1 S [TA ()) s}) S Copips
for some constants Co 0 and Ps (0, 1). finally,
E (0) < 2 2E 1oz 0) Ax] < 2M +2
Σ Σ
Copyp's < 00,
#+++1=kk11
and the bound does not depend on 0.
D.2 COROLLARY
Since E(v, h;; 0) is continuous in O, every conditional distribution of p(v, h; 0) is also continuous in O. This implies that To( is continuous in 0 as well. We have assumed that To is irreducible and aperiodic, so for each 0. To(x x') > 0 for every x, x' E X. By the compactness of θ there exist a constant
0 such that To(x xx)
for all x, x' E X and 0 E e.
Similarly, fx; 0) = dE(v, f;; 0)/00 is continuous in 0, so for any c > O. there is a constant M > 1 such that If (z:; 2 < M for all x X and 0 E O. Then by choosing constant functions ((h) = M, r(v) = M and constants 71 7/2 = 1/2 and L1 = L2 = M/2, we make Assumption hold.
For Assumption 2 the density function g(.) can be chosen as a uniform distribution over the finite space X. Then the proof is complete.
17


Reviews:
[{"review_id": "rkgmSvqotB", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper introduces an efficient, unbiased contrastive divergence-like algorithm for training energy-based generative models on the example of Restricted Boltzmann Machine.\nThe proposed algorithm is built upon a very interesting work on unbiased finite-step MCMC approximations by Jacob et al, 2017.\nDespite the actual theory being published some time ago, the submitted paper popularises these ideas in the machine learning community and contains optimised variants of the existing algorithms for training of RBMs.\n\nThe paper is mostly written well and does a good job of introducing unbiased MCMC estimators.\nAuthors evaluate their method on rather toyish datasets (by modern standards), however, their empirical analysis is thorough. The improvement upon the standard CD and persistent CD is clear.\nIt also appears that the algorithm actually does not require too many steps and generally does not introduce a lot of computational overhead. \nThe only question I have is why CD has only been tried with k=1 steps? \nI would be interested in its performance for different number of steps including the dynamically chosen number provided by the empirical \\tau in UCD for a given iteration.\nEven though I do not expect a significant improvement to be obtained, this would separate the effect of the number of steps chosen “right” from unbiasedness of the gradient estimator.\nOther baselines, including those mentioned in the related work, could also make the comparison more complete.\n\nI would also suggest including https://arxiv.org/abs/1905.04062, as it seems to be relevant in the spirit.", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}, {"review_id": "r1xivrLjYS", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes an algorithmic improvement that significantly simplifies training of energy-based models, such as the Restricted Boltzmann Machine. The key issue in training such models is computing the gradient of the log partition function, which can be framed as computing the expected value of f(x) = dE(x; theta) / d theta over the model distribution p(x). The canonical algorithm for this problem is Contrastive Divergence which approximates x ~ p(x) with k steps of Gibbs sampling, resulting in biased gradients. In this paper, the authors apply the recently introduced unbiased MCMC framework of Jacob et al. to completely remove the bias. The key idea is to (1) rewrite the expectation as a limit of a telescopic sum: E f(x_0) + \\sum_t E f(x_t) - E f(x_{t-1}); (2) run two coupled MCMC chains, one for the “positive” part of the telescopic sum and one for the “negative” part until they converge. After convergence, all remaining terms of the sum are zero and we can stop iterating. However, the number of time steps until convergence is now random.\n\nOther contributions of the paper are:\n1. Proof that Bernoulli RBMs and other models satisfying certain conditions have finite expected number of steps and finite variance of the unbiased gradient estimator.\n2. A shared random variables method for the coupled Gibbs chains that should result in faster convergence of the chains.\n3. Verification of the proposed method on two synthetic datasets and a subset of MNIST, demonstrating more stable training compared to contrastive divergence and persistent contrastive divergence.\n\nI am very excited about this paper and strongly support its acceptance, since the proposed method should revitalize research in energy-based models. While I find the experiments to be somewhat lacking, this is sufficiently offset by the theoretical contributions of the paper.\n\nPros\n1. The paper reads well and introduces all the necessary preliminaries to understand the method. This is important, since I expect many readers to be unfamiliar with the technique.\n2. The proposed method solves an important problem which, as far as I understand, has been the roadblock in large-scale training of RBMs and related models. It is also elegant and fairly straightforward to implement.\n3. The proof of finite computation time and variance is very nice to have. This is because in some cases removing the bias leads to infinite variance, e.g. a parallel submission on SUMO (https://openreview.net/forum?id=SylkYeHtwr).\n\nCons\n1. I don’t think Corollary 1 (convergence of gradient descent to the global optimum) is true for RBMs, as stated on Page 6. This is because the log-likelihood of RBM, or indeed any latent-variable model with permutation-invariant latents, is non-convex. I would suggest removing this corollary and simplifying Algorithm 2 to be regular SGD, as used in the experiments.\n2. There is no experimental comparison of Algorithm 1 (the general version) and Algorithm 3 (the specialized RBM version). It seems intuitive that the specialized version should have lower computation time, but this must be confirmed.\n3. The experimental section may be significantly improved.\n* It is unclear what value of k (number of initial Gibbs steps) from Algorithm 2 is used.\n* The experiments on just the “0” digits of MNIST seem a bit simplistic for the year 2019. It is also not clear what binarization protocol is used.\n* It would be very helpful to provide estimates of the gradient (not log-likelihood) variance of each method to better understand the trade-off between the bias and the variance.\n* I would also like to see the wall-clock time comparison of the methods.\n\nMinor comments\n* Page 1. Of this kind -> of this class. The data distribution p_v (v; theta) -> The model distribution\n* Page 2. Property -> properties. CD-\\tau -- I don’t think you can correctly refer to your method in this way, since it has at least double the computation time of CD for the same number of iterations.\n* Page 3. Provides -> provide. Likelihood gradient -> log-likelihood gradient\n* Algorithm 1 is an infinite loop with no break clause. It would be good to add a break statement after line 5. This would also simplify the discussion of the method.\n* Page 7. I wouldn’t call the fact that CD doesn’t converge on the BAS dataset remarkable, given that it’s been reported by Fischer & Igel 2014.\n* Page 9. The last paragraph stating that the proposed method is not a replacement for CD is confusing. Can you add a short experiment to demonstrate that this combination makes sense?", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly."}, "extra_text": {}}, {"review_id": "BkxfqjaaFH", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Based on recent progress in unbiased MCMC sampling the paper proposes an unbiased contrastive divergence (UCD) algorithm for training energy based models. Specifically they developed an unbiased version of the gibbs sampling contrastive divergence algorithm for training restricted Boltzman machines. The authors demonstrate their method on a toy dataset, simulated data, as well as a reduced version (only the zero digits) of the MNIST dataset and compare the results with the standard Contrastive divergence and Persistent Contrastive Divergence methods.\n\nScore:\nI find the line of work on unbiased estimators important and the (although i’m not an expert) the theory in the paper seems sound. Further the paper is well written and relatively easy to follow. However I do not find the experimental section completely comprehensive and some of the results seem to achieve worse performance than what is reported in the litterature for both the proposed method and baselines (see detailed questions below). Overall I currently score the paper as a weak reject although I can be convinced to bump the score depending on the author feedback.\n\nDetailed Questions:\nExperimental Results:\nQ1) In [Tieleman2008] log-likelihood values for the full MNIST dataset using a) a small model (25 hidden units) where the likelihood is computed exactly and b) an bigger model (500 hidden units) where the likelihood is approximated. On the full MNIST dataset they train using PCD, CD-1, CD-10 and report approximately Log-Likelihoods of -130 and -85 for the small and large models respectively. My questions are:\nQ1.1) In figure 4 you report approximate log-likelihood values on MNIST (only digits zero) of -150 for the different samplers using an RBM with100 hidden units. That seems to be lower performance than the models in [Tieleman2008] while training on a presumably easier dataset?\n\nQ1.2) In figure 4. Can you comment a bit on the variance of your method which seems to be higher, Is there a Bias/Variance trade-off between UCD and e.g PCD?\n\nQ1.3) [Tieleman2008] Reports training times of 1 to 9 Hours for training in on the full MNIST dataset in 2008 and [Hinton 2006] trained large RBMs in 2006. Why is that setting then computationally time-consuming today in your setup - Is there some difference in the setup that I'm missing? \n\nQ1.4) I highly value enlightening small scale experiments and do understand that computational resources are not available everywhere however I think it would benefit the paper greatly if the proposed method is demonstrated on some reasonably sized dataset (at the very least one of full MNIST, Fashion MNIST, FreyFaces).\n\nQ1.5) In Figure 2 you show some interesting figures for the average stopping time and number of rejected samples on the BAS toy dataset. How does these results look on a real dataset like the MNIST zero digit data?\n\n[Tieleman 2008], Training Restricted Boltzmann Machines using Approximations to the Likelihood Gradient,\n[Hinton 2006] Reducing the Dimensionality of Data with Neural Networks\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I did not assess the derivations or theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: Main content:

Blind review #1 summarizes it well:

The paper proposes an algorithmic improvement that significantly simplifies training of energy-based models, such as the Restricted Boltzmann Machine. The key issue in training such models is computing the gradient of the log partition function, which can be framed as computing the expected value of f(x) = dE(x; theta) / d theta over the model distribution p(x). The canonical algorithm for this problem is Contrastive Divergence which approximates x ~ p(x) with k steps of Gibbs sampling, resulting in biased gradients. In this paper, the authors apply the recently introduced unbiased MCMC framework of Jacob et al. to completely remove the bias. The key idea is to (1) rewrite the expectation as a limit of a telescopic sum: E f(x_0) + \sum_t E f(x_t) - E f(x_{t-1}); (2) run two coupled MCMC chains, one for the “positive” part of the telescopic sum and one for the “negative” part until they converge. After convergence, all remaining terms of the sum are zero and we can stop iterating. However, the number of time steps until convergence is now random.

Other contributions of the paper are:
1. Proof that Bernoulli RBMs and other models satisfying certain conditions have finite expected number of steps and finite variance of the unbiased gradient estimator.
2. A shared random variables method for the coupled Gibbs chains that should result in faster convergence of the chains.
3. Verification of the proposed method on two synthetic datasets and a subset of MNIST, demonstrating more stable training compared to contrastive divergence and persistent contrastive divergence.

--

Discussion:

The main objection in reviews was to have meaningful empirical validation of the strong theoretical aspect of the paper, which the authors did during the rebuttal period to the satisfaction of reviewers.

--

Recommendation and justification:

As review #1 said, "I am very excited about this paper and strongly support its acceptance, since the proposed method should revitalize research in energy-based models."
Decision: Accept (Spotlight)