Title: UNBIASED CONTRASTIVE DIVERGENCE ALGORITHM FOR TRAINING ENERGY-BASED LATENT VARIABLE MODELS

Abstract: The contrastive divergence algorithm is a popular approach to training energybased 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.

1. INTRODUCTION:

Energy-based latent variable models cover a broad class of generative models that are frequently used to characterize sophisticated distributions of high-dimensional data. Popular examples of this class include the restricted Boltzmann machines (RBM, Smolensky, 1986; Hinton, 2012) , deep belief nets (Hinton et al., 2006) , and exponential family harmoniums (Welling et al., 2005) , among many others. Energy-based models are complementary to directed generative models such as the variational autoencoders (Kingma & Welling, 2014), and have gained great success in synthesizing realistic data samples (Xie et al., 2016; 2018b) . 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 a visible random vector, v ∈ V ⊂ R p , and a hidden or latent random vector, h ∈ H ⊂ R r , with the density function

EQUATION

where θ ∈ Θ is the unknown parameter vector, E(v, h; θ) is the energy function, and Z(θ) is a normalizing constant to ensure that p(v, h; θ) is a legitimate probability density or mass function.

The model distribution, p v (v; θ), is defined to be the marginal distribution of p(v, h; θ).

Similar to many other machine learning models, the standard approach to estimating the parameter vector θ is the maximum likelihood method. It can be shown that the derivative of the log-likelihood function can be expressed as 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 as 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 are typically consistent in the limiting case, but biased on finite steps, so one needs to run MCMC for a long time to obtain an accurate gradient, which would take tremendous amount of computing time.

To reduce the computational complexity, Hinton (2002) proposed a 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 θ, where k is a fixed integer as small as one. Such an approach is usually referred to as the CD-k algorithm. The simplicity and computational 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 is applied to directed generative models (Ruiz & Titsias, 2019), where CD is used to define a 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-Perpiñán & Hinton, 2005) , and diverges even in some simple models (Schulz et al., 2010; Fischer & Igel, 2010) . The main issue of CD is that the truncation of MCMC produces a biased stochastic gradient for the log-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.

In 2017), the proposed UCD algorithm is a Gibbs-sampler-based training method for energy-based latent variable models. We prove that the stochastic gradient generated by UCD is unbiased with a finite variance, which implies the convergence 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 approach with a solid convergence guarantee. The implementation of the UCD algorithm is available at https://github.com/yixuan/cdtau. The highlights of this article are as follows:

• We develop a new training algorithm for general energy-based latent variable models that include many popular models (e.g. RBM) as 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 theoretical justifications are developed to guarantee its convergence.

• We have tailored a specialized algorithm for RBM, which is shown to significantly reduce the computational cost.

2. A BRIEF REVIEW OF CONTRASTIVE DIVERGENCE:

In this section we 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

(θ; v) = log{p v (v; θ)} = log{ p(v, h; θ)dh}.

Assume that E(v, h; θ) is continuously differentiable for θ, and then with n data points D = (v 1 , . . . , v n ), the derivative of the log-likelihood function (θ; D) = (θ; v i ), also known as the score function, can be written as

EQUATION

where p(D) stands for the empirical distribution of D, and p(h|v; θ) is the conditional distribution of the latent variable h given v = v. A simple derivation of (2) can be found in Fischer & Igel (2014) . Throughout this article we denote x = (v, h) ∈ X := V × H and f (x; θ) = ∂E(v, h; θ)/∂θ. Then the two expectations in (2) can be abbreviated as E D {f (x; θ)} and E M {f (x; θ)}, respectively, where M := p(v, h; θ) is the complete model distribution.

In many cases, for example the RBM model, E D {f (x; θ)} has a closed form, so the major computational difficulty comes from the E M {f (x; θ)} term. A common scheme to approximate this expectation is to run a Markov chain ξ 0 → ξ 1 → • • • with M as the invariant distribution, and then under mild conditions we have lim t→∞ E{f (ξ t ; θ)} = E M {f (x; θ)}. Of course, such a limit cannot be reached in finite steps, so the CD-k algorithm truncates the Markov chain at the k-th step, resulting in the following approximation:

EQUATION

It is easy to see that ∆(θ) ≈ n -1 ∂ (θ; D)/∂θ is a stochastic approximation to the true gradient, so one can use SG to update θ via the iteration

θ i+1 = θ i + α i ∆(θ i )

, where θ i is the parameter estimate in the i-th iteration, and α i is the step size. 

E M {f (x)} = E{f (ξ k )} + ∞ t=k+1 [E{f (ξ t )} -E{f (ξ t-1 )}]

for any fixed k ≥ 0. Now assume that there exists another Markov chain {η t } such that ξ t and η t have the same marginal distributions for all t ≥ 0, and ξ t = η t-1 for all t ≥ τ , where τ is some random time. If we allow the exchange of expectation and summation, then we would get

E M {f (x)} = E f (ξ k ) + ∞ t=k+1 {f (ξ t ) -f (η t-1 )} = E f (ξ k ) + τ -1 t=k+1 {f (ξ t ) -f (η t-1 )} ,

where the first identity holds since E{f (ξ t )} = E{f (η t )} for all t ≥ 0, and the second one is due to the fact that ξ t = η t-1 for t ≥ τ . As a consequence, the quantity f

(ξ k )+ τ -1 t=k+1 {f (ξ t )-f (η t-1

)} is an unbiased estimator for E M {f (x)}. Such an idea seems rather simple, but the construction of the chain {η t }, which we describe in the next section, is a highly non-trivial task.

3.2. COUPLING OF MARKOV CHAINS:

Let M t denote the marginal distribution of a Markov chain {ξ t } at the t-th step. By construction, M t converges to M as t → ∞. To develop the unbiased estimator H k (ξ, η), the second chain {η t } must satisfy two conditions: (1) marginally η t ∼ M t ; (2) {ξ t } and the lag-one sequence {η t-1 } will meet and stay identical after some random time τ . Condition (1) can be trivially met if {ξ t } and {η t } are sampled independently. However, in this way the probability that ξ t = η t-1 may be extremely small, or even be zero for continuous random variables. Therefore, a special joint distribution for (ξ t , η t-1 ) needs to be assigned subject to ξ t ∼ M t and η t-1 ∼ M t-1 . Such a pair of random variables under the marginal distribution constraints is called a coupling, and for our purpose we attempt to seek a coupling scheme such that P (ξ t = η t-1 ) > 0. Figure 1 illustrates the coupling process of two Markov chains {ξ t } and {η t-1 }.

Figure 1 : An illustration of the coupling process. {ξ t } and {η t } start from the same value, and have the same marginal distribution M t at each step. The two chains are correlated in such a way that the event ξ t = η t-1 occurs with a positive probability for each t. After a random time τ (τ = 5 in the illustration), {ξ t } meets {η t-1 } and they stay identical afterwards.

To implement such a coupling, first let {ξ t } and {η t } start from the same initial value ξ 0 = η 0 , and additionally draw ξ 1 ∼ T (•|ξ 0 ), where T (y|x) stands for the transition density function from state x to state y. Next, we need to draw (ξ 2 , η 1 ) such that marginally ξ 2 ∼ M 2 and η 1 ∼ M 1 , which can be a difficult task as M t may not have closed forms. Fortunately, it is much simplified for Markov chains: due to the Markov property, ξ 2 and η 1 will have the requested marginal distributions if we sample ξ 2 |ξ 1 ∼ T (•|ξ 1 ) and η 1 |η 0 ∼ T (•|η 0 ) conditional on ξ 1 and η 0 . 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(•) = T (•|ξ 1 ) and q(•) = T (•|η 0 ), and then the problem reduces to drawing a coupling (ξ, η) such that ξ ∼ p(•), η ∼ q(•), and P (ξ = η) > 0, which can be accomplished via the maximal coupling technique (Appendix A.1). Specific to our problem (2), we need to sample x = (v, h) from p(v, h; θ). In energy-based latent variable models, the most widely-used MCMC method is 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 = h} and h|{v = v} follow multivariate Bernoulli distributions with independent components, which are very easy to sample from. The coupling for Gibbs samplers was briefly mentioned in Jacob et al. (2017) as a special case of the Metropolis-Hastings scheme (Metropolis et al., 1953; Hastings, 1970) , but next we show that some specific structure of Gibbs samplers can be utilized to simplify the process.

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 T v (v|h) := p(v|h; θ) and T h (h|v) := p(h|v; θ). The more sophisticated cases, for example h consists of multiple layers h = (h 1 , . . . , h L ), can be dealt with similarly. In Algorithm 1, we describe the steps to sample two coupled chains {ξ t = (v t , h t )} and {η t = (v t , h t )} based on the Gibbs sampler.

Three remarks are made for Algorithm 1: (1) The meeting event (line 4) only depends on the T v transition density. To verify this, note that at the t-th step, we need to draw

ξ t |ξ t-1 ∼ T (•|ξ t-1 ) and η t-1 |η t-2 ∼ T (•|η t-2 )

, where T (ṽ, h|v, h) = T v (ṽ|h)T h ( h|ṽ) is the transition density for a full update cycle. It is easy to show that T (ṽ, h|v , h )/T (ṽ, h|v, h) = T v (ṽ|h )/T v (ṽ|h), so the T h part cancels in the ratio. (2) Once ξ t and η t-1 meet, they stay identical afterwards, because by then

T v (•|h t-2 ) = T v (•|h t-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 are rejected, so its cost is tiny.

Algorithm 1 Coupling method for the Gibbs sampler:

Input: Densities T v (v|h) and T h (h|v), initial values ξ 0 = (v 0 , h 0 ) = η 0 = (v 0 , h 0 ), T max Output: Coupled chains {ξ t } and {η t } 1: Sample v 1 ∼ T v (•|h 0 ) and h 1 ∼ T h (•|v 1 ). Set ξ 1 = (v 1 , h 1 ) 2: for t = 2, 3, . . . do 3: Sample v t ∼ T v (•|h t-1 ), h t ∼ T h (•|v t )

, and U ∼ Uniform(0, 1)

4: if U ≤ T v (v t |h t-2 )/T v (v t |h t-1 ) or t ≥ T max (maximum stopping time) then 5: Set ξ t = (v t , h t ), η t-1 = ξ t 6: else 7: Sample v t-1 ∼ T v (•|h t-2 ), h t-1 ∼ T h (•|v t-1 ), and U ∼ Uniform(0, 1) until U > T v (v t-1 |h t-1 )/T v (v t-1 |h t-2 ) 8: Set ξ t = (v t , h t ), η t-1 = (v t-1 , h t-1 ) 9:

end if 10: end for

3.3. UNBIASED CONTRASTIVE DIVERGENCE:

The technical tools introduced in Sections 3.1 and 3.2 enable us to develop a new algorithm to train model (1) . Recall that the true gradient of the log-likelihood function is given by (2) . The first term, E D {f (x; θ)}, can be computed exactly, and the second term, E M {f (x; θ)}, is approximated by an unbiased estimator g2 (θ

) := f (ξ k ) + τ -1 t=k+1 {f (ξ t ) -f (η t-1

)}, where the coupled Markov chains {ξ t } and {η t } are generated by Algorithm 1. Assume that the parameter vector θ lies in a closed convex set Θ, and let P Θ (•) denote the projection onto Θ. Putting the pieces together, Algorithm 2 illustrates the UCD algorithm for training energy-based latent variable models. The initial chain length k can be any fixed number, and in this article we take k = 1 for all the numerical experiments. Set ξ 0 = η 0 = (v, h), and run Algorithm 1 with θ = θ i until ξ τi = η τi-1 4:

g(θ) ← -E h∼p(h|v;θ) {f (v, h; θ)} + f (ξ k ) + τi-1 t=k+1 {f (ξ t ) -f (η t-1 )} 5: θ i+1 ← P Θ (θ i + α i • g(θ i )) 6: end for 7: return θ = T -1 T i=1 θ i

Next, we analyze the theoretical property of Algorithm 2 and state the conditions for it to converge. As a standard setting, we assume that the Markov chains generated by the Gibbs sampler are ϕ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 of functions r :

V → [1, +∞), l : H → [1, +∞) and constants γ 1 , γ 2 , L 1 , L 2 > 0 such that γ 1 γ 2 < 1 and E v∼p(v|h;θ) r(v) ≤ γ 1 l(h) + L 1 , E h∼p(h|v;θ) l(h) ≤ γ 2 r(v) + L 2 , ∀v ∈ V, h ∈ H, θ ∈ Θ.

Also, there exist constants c > 0 and

D > 0 such that |f (x; θ)| 2+c ≤ l(h) and E h∼p h (h;θ) l(h) ≤ D for all x = (v, h) ∈ X and θ ∈ Θ. Assumption 2. (Minorization condition) There exist constants d > 2(γ 2 L 1 + L 2 )/(1 -γ 1 γ 2 ),

ε > 0, and a density function q(•) such that p(v|h; θ) ≥ εq(v) for all h ∈ D, v ∈ V, and θ ∈ Θ, where D = {h ∈ H : l(h) ≤ d}.

In the following theorem we show three important facts about the proposed stochastic gradient g(θ):

(1) g(θ) is unbiased for the true score function; (2) it has a bounded second moment uniformly in θ;

(3) in expectation it can be computed in finite time.

Theorem 1. Under Assumptions 1 and 2, there exist constants

D 1 , D 2 > 0 such that E{g(θ)} = ∂ (θ; v)/∂θ, E {g 2 (θ)} 2 ≤ D 1 , and E(τ i ) ≤ D 2 for all θ ∈ Θ and i = 1, 2, . . . , T -1.

Theorem 1 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 a typical setting, in the following corollary we consider a convex log-likelihood function.

Corollary 1. Assume that (θ; v) is convex and L-Lipschitz continuous in θ ∈ Θ, and Θ is a closed and bounded convex set. Then by choosing

α i = α 0 / √ i for some constant α 0 > 0, we have * -( θ; v) ≤ O(1/ √ T )

, where * is the maximum value of (θ; v). 

4. TRAINING RESTRICTED BOLTZMANN MACHINES:

RBM is one of the most popular and widely-used energy models in machine learning, defined by the energy function and h|{v = v} ∼ Bernoulli(σ(W T v + c)). The coupling method in Algorithm 1 directly works for RBM, but here we show an improved version that is tailored for RBM and is more efficient.

E(v, h; θ) = -v T b -v T W h -h T c, where v ∈ {0, 1} m , h ∈ {0, 1} n

Let u, p ∈ R r , and the notation y = 1{u ≤ p} stands for a binary vector such that y i = 1 if u i ≤ p i and y i = 0 otherwise. Also let T v (v|h) = m i=1 p vi i (1 -p i ) 1-vi denote the transition density from h to v, where p = (p 1 , . . . , p m ) T = σ(W h + b). Then the specialized coupling method for RBM is given in Algorithm 3.

Algorithm 3 Coupling method for RBM:

Input: Model parameters W , b, c, step-t states ξ t = (v t , h t ), η t-1 = (v t-1 , h t-1 ) Output: New states ξ t+1 = (v t+1 , h t+1 ), η t = (v t , h t ) 1: Sample U 1 ∼ Uniform(0, 1), Z 1 ∼ Uniform([0, 1] m ), and set v t+1 = 1{Z 1 ≤ σ(W h t + b)} 2: if U 1 ≤ T v (v t+1 |h t-1 )/T v (v t+1 |h t ) then 3: Set v t = v t+1 4: else 5: repeat 6: Sample U 2 ∼ Uniform(0, 1), U 2 ∼ Uniform(0, 1), Z 2 ∼ Uniform([0, 1] m ) 7:

if v t+1 has not been accepted then

8: Propose v t+1 = 1{Z 2 ≤ σ(W h t + b)}, accept if U 2 > T v (v t+1 |h t-1 )/T v (v t+1 |h t ) 9: end if 10: if v t has not been accepted then 11: Propose v t = 1{Z 2 ≤ σ(W h t-1 + b)}, accept if U 2 > T v (v t |h t )/T v (v t |h t-1 ) 12: end if 13:

until v t+1 and v t are both accepted 14: end if

15: Sample Z 3 ∼ Uniform([0, 1] n ) 16: Set h t+1 = 1{Z 3 ≤ σ(W T v t+1 + c)}, h t = 1{Z 3 ≤ σ(W T v t + c)}

The intuition behind Algorithm 3 is the following: line 2 indicates that it is also a maximal coupling method, so the probability P (ξ t+1 = η t ) is the same as Algorithm 1. However, in the event {ξ t+1 = η t }, ξ t+1 and η t are independent in Algorithm 1 but correlated in Algorithm 3, achieved by the use of common random variates Z 2 and Z 3 . The correlation between ξ t+1 and η t helps to make P (ξ t+2 = η t+1 ) larger, thus accelerating the meeting of {ξ t } and {η t-1 }. A more rigorous justification of this algorithm is given in Appendix A.2.

Finally, it is known that in the gradient expression (2), f (x; θ) = (σ(W h + b)h T , σ(W h + b), h) for RBM, corresponding to the parameters θ = (W , b, c). With the coupled chains {ξ t } and {η t-1 }, RBM can then be trained using UCD given by Algorithm 2.

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 a 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 to estimate the normalizing constant, which is consistent with a large sample. However, it still induces a 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 statistics and machine learning. Some background knowledge in this article, for example Section 3.1, is taken from Jacob et al. ( 2017), which established a general framework for unbiased MCMC. Our new contributions are in the following aspects. First, we have developed Algorithm 1 and Theorem 1 exclusively for the Gibbs sampler, taking into account the special structure of Gibbs MCMC. Second, our theoretical results, including Theorem 1 and Corollary 2, have more practical assumptions than the ones in Jacob et al. (2017) . For example, one of their key assumptions is that E{|f (ξ t )| 2+c } is uniformly bounded for every finite step, which is quite abstract and hard to verify in practice compared with our Assumption 1. Third, Jacob et 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 4 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 method1 (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 is 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-STRIPES 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-likelihood value can be evaluated exactly, and it demonstrates the divergence of CD-based training algorithms. In our study, k is set to 1 for CD (more experiments with larger k are given in Appendix B.1), and each algorithm is run for 100 times, accounting for the randomness in the training process. A common learning Published as a conference paper at ICLR 2020 rate α = 0.01 is set, and 1000 parallel Markov chains are used to approximate the gradient in each iteration. The results are shown in Figure 2 .

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 2 shows the following findings. First, we 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 an adaptive choice of the stopping time in the Markov chain, compared to the fixed k in CD. In the BAS data, the stopping time has a steep increase around the 1200th iteration. Interestingly, this is exactly where CD begins to fail. An interpretation of this phenomenon is that UCD automatically uses a 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 is 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 B.1 for more discussions). Third, the cost of the rejection sampling step in UCD (line 7 of Algorithm 1) is tiny, as the number of rejected samples rarely goes above two. Finally, UCD does not see a 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 we show that the findings for the BAS data can be observed in other model settings. We simulate a data set from an RBM model with 200 visible units and 20 hidden units, where the entries of weight and bias parameters are all generated from a N (0, 1) distribution. The sample size of the simulated data set is 1000, and we fit an RBM model using 100 hidden units, which is larger than the true model since we intend to mimic the common practice of overparameterization 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 are given in Appendix B.2.

Figure 3 : Approximate log-likelihood values for each algorithm on the simulated RBM data set.

6.3. FASHION-MNIST DATA:

Next we consider the Fashion-MNIST data set2 , a 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 a 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 4 demonstrates the training trajectories of the three algorithms. The training patterns for CD and PCD are surprising: the log-likelihood values of CD are decreasing, and PCD seems to bounce back and forth among three different paths. In fact, in Appendix B.3 we show that CD-k does not converge even if k is as large as 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 analysis and the numerical experiments at different scales demonstrate that the impact of bias elimination is huge.

It is expected that UCD may have slightly larger variance compared with CD and PCD, but we emphasize that the value of UCD is not a simple question of bias and variance trade-off. This is because 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 a similar log-likelihood level. Therefore, UCD is not meant to completely replace CD or PCD, but rather to serve as an important addition to the existing training algorithms. In practice, a very useful technique is to first run the fast CD-1 or PCD to the near-optimum, and then proceed with UCD for guaranteed convergence.

A DISTRIBUTION COUPLING METHODS:

A.1 THE MAXIMAL COUPLING ALGORITHM Let (ξ, η) be a pair of random variables defined on the same probability space, and p and q be two distributions. (ξ, η) is called a coupling of p and q if marginally ξ ∼ p and η ∼ q. It is well known (see e.g. Proposition 4.7 of Levin & Peres, 2017) that among all possible joint distributions of (ξ, η),

EQUATION

where p -q TV is the total variation distance between p and q. The maximal coupling algorithm, given in Algorithm 4, generates a coupling (ξ, η) that achieves the bound in (4). That is, it maximizes the probability that two random variables are equal subject to their marginal distributions.

Algorithm 4 Maximal coupling of two distributions p(•) and q(•), from Jacob et al. (2017).

Input: Density functions p(•) and q(•)

Output: A coupling (ξ, η) with ξ ∼ p(•) and η ∼ q(•) 1: Sample ξ ∼ p(•) and U ∼ Uniform(0, 1) independently 2: if U ≤ q(ξ)/p(ξ) then 3:

return (ξ, ξ) 4: else 5:

Sample η ∼ q(•) and U ∼ Uniform(0, 1) independently until U > p(η)/q(η)

6: return (ξ, η) 7: end if A.2 COUPLING METHOD FOR RBM

In the general maximal coupling method (Algorithm 4), let p(•) = T (•|ξ t ) and q(•) = T (•|η t-1 ) be the transition densities of RBM, and then it generates the new states ξ t+1 and η t with the probability P (ξ t+1 = η t |ξ t , η t-1 ) maximized. If the event {ξ t+1 = η t } does not happen, then ξ t+1 and η t are sampled independently.

However, if ξ t+1 and η t are close to each other, then in the next iteration the probability P (ξ t+2 = η t+1 |ξ t+1 , η t ) would be large, which helps to shorten the stopping time. Therefore, we are motivated to minimize some type of distance between ξ t+1 and η t . For RBM, we characterize it by E( v t+1v t 2 |v t+1 = v t ) and E( h t+1 -h t 2 |v t+1 , v t ). 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.

First consider E( h t+1 -h t 2 |v t+1 , v t ). Given the v variables, h t+1 and h t follow Bernoulli distributions elementwisely, with individual mean vectors denoted as µ 1 and µ 2 , respectively. It can be shown that the maximal coupling of two Bernoulli variables is achieved by using the same random variate, so

E( h t+1 -h t 2 |v t+1 , v t ) is minimized by setting h t+1 = 1{Z ≤ µ 1 } and h t = 1{Z ≤ µ 2 }, where Z ∼ Uniform([0, 1] n ).

This leads to lines 15-16 of Algorithm 3.

Next, conditional on {v t+1 = v t }, the maximal coupling algorithm generates v t+1 and v t with marginal distributions proportional to p

1 (v) -min{p 1 (v), p 2 (v)} and p 2 (v) -min{p 1 (v), p 2 (v)}, respectively, where p 1 (•) = T v (•|h t ) and p 2 (•) = T v (•|h t-1

). Algorithm 3 implements this by two rejection sampling steps (line 8 and line 11). If we use the common random variates Z 2 ∼ Uniform([0, 1] n ) for the proposals, then v t+1 and v t are maximally coupled in the event that they are accepted in the same iteration. This property helps to reduce

E( v t+1 -v t 2 |v t+1 = v t ).

To demonstrate that the specialized Algorithm 3 improves the general Algorithm 1, 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 N (0, 0.1 2 ) distribution. For both algorithms, the elements of the initial state v 0 are sampled from independent Bernoulli(0.5) distributions, and we generate coupled Markov chains until the stopping time τ is reached or the chain length exceeds 1000. Figure 5 shows the distributions of τ for both algorithms based on 1000 replications. It is very clear that the general algorithm leads to a surprisingly long stopping time, and even has a long tail beyond 1000. In contrast, 65.4% of the stopping times in the specialized algorithm are smaller than or equal to 10. In this sense the improvement brought about by Algorithm 3 is huge.

B ADDITIONAL RESULTS FOR NUMERICAL EXPERIMENTS:

All experiments in this article were run on an Intel Xeon Gold 6126 processor with 12 cores and 24 threads. CD and PCD algorithms used the OpenBLAS library3 for parallel matrix computations, and UCD used OpenMP4 to generate Markov chains in parallel.

B.1 BAS DATA:

For the 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 6 ). 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 a fixed-k CD algorithm with a similar performance. The running time for each training algorithm, which is given in Table 1 , also supports this claim. 

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 8 gives the average stopping time for UCD in each iteration, and the number of 2 illustrates the computational time for each algorithm. We can find that the computational cost of UCD is slightly larger than but very close to CD-1 and PCD. One useful technique to combine the computational efficiency of 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-MNIST data, Figure 10 gives the training trajectories of CD-10, CD-20, and CD-30 algorithms, but unfortunately none of them have a convergent pattern. Figure 11 shows the average stopping time and the number of discarded samples in the UCD algorithm. Again, the cost for rejection sampling can be ignored, and the value of τ gets stable around 30 after 2000 iterations. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 9: Combine PCD and UCD for efficient and convergent training.

Finally in Table 3 we demonstrate the running time for each training algorithm. Although UCD takes more time than other methods, the additional cost seems to be the necessary price to achieve a convergent result, as all other methods compared have a divergent log-likelihood. 

C VARIANCE OF GRADIENT ESTIMATES:

To compare the variances of stochastic gradients generated by different algorithms, in each optimization iteration we estimate the following gradient variance defined by where f (θ) is an estimator for the second term of (2), and • F is the Frobenius norm. Among the training methods considered in this article, f (θ) is biased in CD and PCD, and is unbiased in UCD.

V ( f (θ)) = E f (θ) -E( f (θ)) 2 F ,

We consider the example in Section 6.2, and estimate V ( f (θ)) using the 1000 parallel Markov chains in each iteration. Figure 12 shows the estimated gradient variances for b, c, and W along the training process.

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 a 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 1:

Omitting the iteration index i for brevity, we first prove that the stopping time τ has a finite expectation. Let T θ (x|x ) := T θ (v, h|v , h ) = p(v|h ; θ)p(h|v; θ) denote the transition density for a full update cycle. Under Assumptions 1 and 2, Lemma 2 of Johnson & Burbank (2015) shows that there exists a constant γ ∈ [γ 1 γ 2 , 1) such that

EQUATION

for all x ∈ X. Therefore, the drift condition holds for the transition density T θ (x|x ).

Next, since Algorithm 1 is a maximal coupling algorithm, we have (see for example Jacob et al., 2017)

P (ξ t+1 = η t |ξ t = x, η t-1 = x ) = min{T θ (y|x), T θ (y|x )}dy.

By Assumption 2, p(v|h; θ) ≥ εq(v) for all h ∈ D, v ∈ V, and θ ∈ Θ, so P (ξ t+1 = η t |ξ t = x, η t-1 = x ) ≥ εq(ṽ)p( h|ṽ; θ)dy on x, x ∈ D × V, where y = (ṽ, h). Note that q(ṽ)p( h|ṽ; θ) is a joint density function, so q(ṽ)p( h|ṽ; θ)dy = 1, and hence P (ξ t+1 = η t |ξ t = x, η t-1 = x ) ≥ ε for all h ∈ D and θ ∈ Θ.

Then by Proposition 3.4 of Jacob et al. (2017) , there exist constants κ 1 > 0 and ρ 1 ∈ (0, 1) such that for all t > 0, P (τ > t) ≤ κ 1 l(h 0 )ρ t 1 . As a result, we obtain

EQUATION

The right hand side of (6) does not depend on the value of θ, so E(τ ) is uniform in the iteration index i.

To make sure that g(θ) is well defined, a few regularity conditions, formulated as Assumption 2. Finally, we need to show that the second moment of the stochastic gradient is bounded uniformly in θ. Let ∆ k = f (ξ k ) and ∆ t = f (ξ t ) -f (η t-1 ) for t ≥ k + 1, and then g2 (θ) = ∞ t=k ∆ t . From the ergodicity property (7) we immediately get E(∆ 2 k ) ≤ M . For t ≥ k + 1,

E(|∆ t | 2+c ) = E(|f (ξ t ) -f (η t-1 )| 2+c ) ≤ 2 1+c E(|f (ξ t )| 2+c + |f (η t-1 )| 2+c ) ≤ 2 2+c M.

Then by Hölder's inequality,

E(∆ 2 t ) = E(∆ 2 t 1{τ > t}) ≤ E(|∆ t | 2+c ) 2/(c+2) {P (τ > t)} c/(c+2)

≤ 4M 2/(c+2) κ 1 l(h 0 )ρ t 1 c/(c+2) .

Since g2 (θ) = ∞ t=k ∆ t , we get

E {g 2 (θ) -∆ k } 2 = E    ∞ t=k+1 ∆ t 2    ≤ ∞ t=k+1 ∞ s=k+1 E |∆ t ∆ s 1{τ > t}1{τ > s})| .

We also have

E |∆ t ∆ s 1{τ > t}1{τ > s})| ≤ E(∆ 2 t 1{τ > t})E(∆ 2 s 1{τ > s}) ≤ C 0 ρ s 3 ρ t 3

for some constants C 0 > 0 and ρ 3 ∈ (0, 1), so finally,

E {g 2 (θ)} 2 ≤ 2E(∆ 2 k ) + 2E {g 2 (θ) -∆ k } 2 ≤ 2M + 2 ∞ t=k+1 ∞ s=k+1 C 0 ρ s 3 ρ t 3 < ∞,

and the bound does not depend on θ.

D.2 COROLLARY 2:

Since E(v, h; θ) is continuous in θ, every conditional distribution of p(v, h; θ) is also continuous in θ. This implies that T θ (x|x ) is continuous in θ as well. We have assumed that T θ is irreducible and aperiodic, so for each θ, T θ (x|x ) > 0 for every x, x ∈ X. By the compactness of Θ, there exist a constant ε > 0 such that T θ (x|x ) ≥ ε for all x, x ∈ X and θ ∈ Θ.

Similarly, f (x; θ) = ∂E(v, h; θ)/∂θ is continuous in θ, so for any c > 0, there is a constant M > 1 such that |f (x; θ)| 2+c ≤ M for all x ∈ X and θ ∈ Θ. Then by choosing constant functions l(h) = M , r(v) = M and constants γ 1 = γ 2 = 1/2 and L 1 = L 2 = M/2, we make Assumption 1 hold.

For Assumption 2, the density function q(•) can be chosen as a uniform distribution over the finite space X. Then the proof is complete.

:

We were informed of this method by a public comment on https://openreview.net/forum?id= r1eyceSYPr, after the current paper was accepted.

https://research.zalando.com/welcome/mission/research-projects/ fashion-mnist/

https://www.openblas.net/

https://www.openmp.org/

Figures and Tables:

Algorithm for estimating θ Input: T , {α i }, k, initial value θ 0 Output: Parameter estimate for θ 1: for i = 0, 1, . . . , T -1 do Draw one data point v ∼ p(D), and sample h ∼ p(h|v; θ i ) 3:

, and θ = (W , b, c) are model parameters. The Gibbs sampler for RBM has a nice structure: let σ(x) = 1/(1 + exp(-x)) be the sigmoid function, and then v|{h = h} ∼ Bernoulli(σ(W h + b))

Figure 4: Approximate log-likelihood values for each algorithm on the Fashion-MNIST data set.

Figure 5: Distribution of the stopping time τ for the general (Algorithm 1) and specialized (Algorithm 3) coupling methods.

Figure 6: CD-k algorithms for training the BAS data with different k.

Figure 7: Approximate log-likelihood values for each algorithm on the simulated RBM data set.

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.

Figure 10: Approximate log-likelihood values for CD algorithms on the Fashion-MNIST data set.

Figure 11: Left: average stopping time τ for UCD in training the model on Fashion-MNIST data. Right: average number of rejected samples in the coupling algorithm for UCD.

this article, we propose a new unbiased contrastive divergence (UCD) algorithm based on recent advances in unbiased MCMC theory, which offers new possibilities for solving the model training problem. In the seminal work Glynn & Rhee (2014), the authors developed an unbiased estimator for the expectation with respect to the invariant distribution of a Markov chain. More recently, this estimator was further extended to the MCMC setting by Jacob et al. (2017), using a 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.

Table content: <table/>

In what follows we omit the dependence on θ for brevity if no confusion is caused. If a Markov chain {ξ t } satisfies E{f (ξ t )} → E M {f (x)} as t → ∞, then under some regularity conditions, we can express the limit as a telescoping sum,

Table content: <table><tr><td>Despite its simplicity, various research articles have pointed out the weaknesses of the CD-k algo-</td></tr><tr><td>rithm. For instance, Sutskever &amp; Tieleman (2010) gave an example to show that E{∆(θ)} is not the</td></tr><tr><td>gradient of any objective function, and Schulz et al. (2010); Fischer &amp; Igel (2014) studied numerical</td></tr><tr><td>experiments in which CD-k does not converge at all for small k values. Carreira-Perpiñán &amp; Hin-</td></tr><tr><td>ton (2005) considered the fixed points of ∆(θ), the θ values such that E{∆(θ)} = 0, and showed</td></tr><tr><td>that they do not match the fixed points of ∂ (θ; D)/∂θ in general. This implies that even if CD-k</td></tr><tr><td>converges, the resulting parameter estimate may not be a local minimum of the likelihood function.</td></tr><tr><td>Another variant of CD is the persistent contrastive divergence (PCD, Tieleman, 2008; Tieleman &amp;</td></tr><tr><td>Hinton, 2009), which has been reported to improve CD in many numerical experiments. However,</td></tr><tr><td>it is still an approximation method, and its convergence properties are more difficult to analyze, as</td></tr><tr><td>the stochastic gradients generated by PCD become correlated across iterations. In fact, Schulz et al.</td></tr><tr><td>(2010); Fischer &amp; Igel (2010) also gave examples in which PCD failed to converge. There are also</td></tr><tr><td>some other training methods as extensions to CD, such as the multi-grid method (Gao et al., 2018)</td></tr><tr><td>and the short-run MCMC (Nijkamp et al., 2019), but all these methods inherit the bias of CD.</td></tr><tr><td>To summarize, it is surprising that virtually none of the popular training methods for energy-based</td></tr><tr><td>models, including CD and PCD, provide a solid convergence guarantee. The major defects of CD</td></tr><tr><td>stem from the fact that ∆(θ) is a biased estimator for the true log-likelihood gradient, and SG may</td></tr><tr><td>fail with uncontrolled bias accumulation. To this end, the ultimate solution is to design a training</td></tr><tr><td>algorithm that completely removes the bias of CD.</td></tr><tr><td>3 THE UNBIASED CONTRASTIVE DIVERGENCE ALGORITHM</td></tr><tr><td>3.1 UNBIASED MCMC ESTIMATORS</td></tr><tr><td>Since CD highly relies on the MCMC method, the main ingredient of the proposed UCD algorithm is</td></tr><tr><td>the theory of unbiased MCMC developed by Glynn &amp; Rhee (2014) and Jacob et al. (2017). Consider</td></tr><tr><td>the second term in (2), namely, E M {f (x; θ)}.</td></tr></table>

The proof Corollary 1 is standard, see for exampleBottou (2010);Bottou et al. (2018). When the log-likelihood function is nonconvex as in the RBM model, there are also other versions of the convergence result for SG, for example Theorem 4.10 ofBottou et al. (2018). Such directions can be studied separately and are omitted here.Finally, we shall point out an important special case of Theorem 1, 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 l(•) functions. We summarize this useful fact in the following corollary.Corollary 2. If X is a finite state space and Θ is compact, then Assumptions 1 and 2 hold, and Theorem 1 applies.

Table content: <table/>

Running time (in seconds) for each training algorithm on the BAS experiment. CD-1 w/ log-likelihood values CD-1 CD-2 CD-5 CD-10 CD-20 PCD UCD

Table content: <table><tr><td>65.01</td><td>5.38</td><td>8.44 17.89 33.17</td><td>63.51</td><td>5.18 15.99</td></tr><tr><td colspan="3">discarded samples in the rejection sampling step. Table</td><td/><td/></tr></table>

Running time (in seconds) for each training algorithm on the experiment in Section 6.2.

Table content: <table><tr><td colspan="2">CD-1 w/ log-likelihood values CD-1 CD-2 CD-5 CD-10 PCD UCD</td></tr><tr><td>216.46</td><td>16.52 24.44 48.54 89.09 16.18 28.47</td></tr></table>

Running time (in minutes) for each training algorithm on the Fashion-MNIST experiment.

Table content: <table><tr><td colspan="2">CD-1 w/ log-likelihood values CD-1 CD-10 CD-20 CD-30 PCD</td><td>UCD</td></tr><tr><td>65.14</td><td colspan="2">25.63 114.85 212.22 311.89 25.53 494.01</td></tr></table>

1 ofJacob et al. (2017), need to be verified. Under the drift condition (5) and Assumption 2, Theorem 12 ofRosenthal (1995) shows that {ξ t } is a geometrically ergodic Markov chain. Therefore, there exist constants κ 2 > 0 and ρ 2 ∈ (0, 1) such that|Eφ(ξ t ) -E M φ| ≤ κ 2 l(h 0 )ρ t : |φ(v, h)| ≤ l(h). Therefore, there exists a constant M > 0 such that |Eφ(ξ t )| ≤ E M l + κ 2 l(h 0 )ρ t 2 ≤ D + κ 2 l(h 0 ) ≤ M < ∞. Since |f (x; θ)| 2+c ≤ l(h) by Assumption 1, we get E{|f (ξ t )| 2+c } ≤ Mfor all t > 0. Moreover, (7) implies that E{f (ξ t )} → E M f . These two results verify Assumption 2.1 of Jacob et al. (2017), and hence the unbiasedness of g(θ) is true by design.

Table content: <table><tr><td>2</td><td>(7)</td></tr><tr><td>for all φ</td><td/></tr></table>

Peer review (iclr2020):
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. I 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. Pros 1. 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. 2. 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. 3. 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). Cons 1. 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. 2. 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. 3. The experimental section may be significantly improved. * It is unclear what value of k (number of initial Gibbs steps) from Algorithm 2 is used. * 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. * 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. * I would also like to see the wall-clock time comparison of the methods. Minor comments * Page 1. Of this kind -> of this class. The data distribution p_v (v; theta) -> The model distribution * Page 2. Property -> properties. CD-	au -- 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. * Page 3. Provides -> provide. Likelihood gradient -> log-likelihood gradient * 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. * 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. * 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?