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

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.

heading: 1 INTRODUCTION
text: 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 ⊂ Rp, and a hidden or latent random vector, h ∈ H ⊂ Rr, with the density function
p(v,h;θ) = 1
Z(θ) exp{−E(v,h;θ)}, (1)
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, pv(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-Perpiñá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 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.
Under the framework of Glynn & Rhee (2014) and Jacob et al. (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.

heading: 2 A BRIEF REVIEW OF CONTRASTIVE DIVERGENCE
text: 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{pv(v;θ)} = log{ ∫ p(v,h;θ)dh}. Assume that E(v,h;θ) is continuously differentiable for θ, and then with n data points D = (v1, . . . ,vn), the derivative of the log-likelihood function `(θ;D) = ∑ `(θ;vi), also known as the score function, can be written as
∂`(θ;D) ∂θ
= −n [ E(v,h)∼p(D)p(h|v;θ) { ∂E(v,h;θ)
∂θ
} − E(v,h)∼p(v,h;θ) { ∂E(v,h;θ)
∂θ
}] , (2)
where p(D) stands for the empirical distribution ofD, 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 ED{f(x;θ)} and EM{f(x;θ)}, respectively, whereM := p(v,h;θ) is the complete model distribution. In many cases, for example the RBM model, ED{f(x;θ)} has a closed form, so the major computational difficulty comes from the EM{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 limt→∞ E{f(ξt;θ)} = EM{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:
∆(θ) := − [ED{f(x;θ)} − f(ξk;θ)] . (3) 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.
Despite its simplicity, various research articles have pointed out the weaknesses of the CD-k algorithm. For instance, Sutskever & Tieleman (2010) gave an example to show that E{∆(θ)} is not the gradient of any objective function, and Schulz et al. (2010); Fischer & Igel (2014) studied numerical experiments in which CD-k does not converge at all for small k values. Carreira-Perpiñán & Hinton (2005) considered the fixed points of ∆(θ), the θ values such that E{∆(θ)} = 0, and showed that they do not match the fixed points of ∂`(θ;D)/∂θ 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 an 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 are also some other training methods as extensions to CD, such as the multi-grid method (Gao et al., 2018) and the short-run MCMC (Nijkamp et al., 2019), but all these methods inherit the bias of CD.
To summarize, it is surprising that virtually none of the popular training methods for energy-based models, including CD and PCD, provide a 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 to design a training algorithm that completely removes the bias of CD.

heading: 3 THE UNBIASED CONTRASTIVE DIVERGENCE ALGORITHM

heading: 3.1 UNBIASED MCMC ESTIMATORS
text: Since CD highly relies on the MCMC method, the main ingredient of the proposed UCD algorithm is the theory of unbiased MCMC developed by Glynn & Rhee (2014) and Jacob et al. (2017). Consider the second term in (2), namely, EM{f(x;θ)}. In what follows we omit the dependence on θ for brevity if no confusion is caused. If a Markov chain {ξt} satisfies E{f(ξt)} → EM{f(x)} as t→∞, then under some regularity conditions, we can express the limit as a telescoping sum,
EM{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
EM{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 EM{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.

heading: 3.2 COUPLING OF MARKOV CHAINS
text: LetMt denote the marginal distribution of a Markov chain {ξt} at the t-th step. By construction, Mt converges toM as t→∞. To develop the unbiased estimator Hk(ξ, η), the second chain {ηt} must satisfy two conditions: (1) marginally ηt ∼Mt; (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 ∼ Mt and ηt−1 ∼ Mt−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}.
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 ∼ M2 and η1 ∼ M1, which can be a difficult task asMt 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 forM follow the natural division of blocks x = (v,h). That is, one can easily sample from the two transition distributions Tv(v|h) := p(v|h;θ) and Th(h|v) := p(h|v;θ). The more sophisticated cases, for example h consists of multiple layers h = (h1, . . . ,hL), can be dealt with similarly. In Algorithm 1, we describe the steps to sample two coupled chains {ξt = (vt,ht)} 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 Tv 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 (ṽ, h̃|v,h) = Tv(ṽ|h)Th(h̃|ṽ) is the transition density for a full update cycle. It is easy to show that T (ṽ, h̃|v′,h′)/T (ṽ, h̃|v,h) = Tv(ṽ|h′)/Tv(ṽ|h), so the Th part cancels in the ratio. (2) Once ξt and ηt−1 meet, they stay identical afterwards, because by then Tv(·|h′t−2) = Tv(·|ht−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 Tv(v|h) and Th(h|v), initial values ξ0 = (v0,h0) = η0 = (v′0,h ′ 0), Tmax Output: Coupled chains {ξt} and {ηt} 1: Sample v1 ∼ Tv(·|h0) and h1 ∼ Th(·|v1). Set ξ1 = (v1,h1) 2: for t = 2, 3, . . . do 3: Sample vt ∼ Tv(·|ht−1), ht ∼ Th(·|vt), and U ∼ Uniform(0, 1) 4: if U ≤ Tv(vt|h′t−2)/Tv(vt|ht−1) or t ≥ Tmax (maximum stopping time) then 5: Set ξt = (vt,ht), ηt−1 = ξt 6: else 7: Sample v′t−1 ∼ Tv(·|h ′ t−2), h ′ t−1 ∼ Th(·|v′t−1), and U ′ ∼ Uniform(0, 1)
until U ′ > Tv(v′t−1|ht−1)/Tv(v′t−1|h ′ t−2)
8: Set ξt = (vt,ht), ηt−1 = (v′t−1,h ′ t−1)
9: end if 10: end for

heading: 3.3 UNBIASED CONTRASTIVE DIVERGENCE
text: 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, ED{f(x;θ)}, can be computed exactly, and the second term, EM{f(x;θ)}, is approximated by an unbiased estimator g̃2(θ) := 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.
Algorithm 2 UCD Algorithm for estimating θ Input: T , {αi}, k, initial value θ0 Output: Parameter estimate for θ
1: for i = 0, 1, . . . , T − 1 do 2: Draw one data point v ∼ p(D), and sample h ∼ p(h|v;θi) 3: Set ξ0 = η0 = (v,h), and run Algorithm 1 with θ = θi until ξτi = ητi−1 4: g̃(θ)← −Eh∼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, L1, L2 > 0 such that γ1γ2 < 1 and
Ev∼p(v|h;θ)r(v) ≤ γ1l(h) + L1, Eh∼p(h|v;θ)l(h) ≤ γ2r(v) + L2, ∀v ∈ V, h ∈ H, θ ∈ Θ.
Also, there exist constants c > 0 andD > 0 such that |f(x;θ)|2+c ≤ l(h) and Eh∼ph(h;θ)l(h) ≤ D for all x = (v,h) ∈ X and θ ∈ Θ. Assumption 2. (Minorization condition) There exist constants d > 2(γ2L1 + L2)/(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 D1, D2 > 0 such that E{g̃(θ)} = ∂`(θ;v)/∂θ, E [ {g̃2(θ)}2 ] ≤ D1, and E(τi) ≤ D2 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).
The proof Corollary 1 is standard, see for example Bottou (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 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 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.

heading: 4 TRAINING RESTRICTED BOLTZMANN MACHINES
text: RBM is one of the most popular and widely-used energy models in machine learning, defined by the energy function E(v,h;θ) = −vTb − vTWh − hTc, where v ∈ {0, 1}m, h ∈ {0, 1}n, 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(σ(Wh+ b)) and h|{v = v} ∼ Bernoulli(σ(WTv + 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.
Let u,p ∈ Rr, and the notation y = 1{u ≤ p} stands for a binary vector such that yi = 1 if ui ≤ pi and yi = 0 otherwise. Also let Tv(v|h) = ∏m i=1 p vi i (1− pi)1−vi denote the transition density from h to v, where p = (p1, . . . , pm)T = σ(Wh+ b). Then the specialized coupling method for RBM is given in Algorithm 3.
Algorithm 3 Coupling method for RBM Input: Model parametersW , b, c, step-t states ξt = (vt,ht), ηt−1 = (v′t−1,h ′ t−1) Output: New states ξt+1 = (vt+1,ht+1), ηt = (v′t,h ′ t)
1: Sample U1 ∼ Uniform(0, 1), Z1 ∼ Uniform([0, 1]m), and set vt+1 = 1{Z1 ≤ σ(Wht + b)} 2: if U1 ≤ Tv(vt+1|h′t−1)/Tv(vt+1|ht) then 3: Set v′t = vt+1 4: else 5: repeat 6: Sample U2 ∼ Uniform(0, 1), U ′2 ∼ Uniform(0, 1), Z2 ∼ Uniform([0, 1]m) 7: if vt+1 has not been accepted then 8: Propose vt+1 = 1{Z2 ≤ σ(Wht + b)}, accept if U2 > Tv(vt+1|h′t−1)/Tv(vt+1|ht) 9: end if
10: if v′t has not been accepted then 11: Propose v′t = 1{Z2 ≤ σ(Wh ′ t−1 + b)}, accept if U ′2 > Tv(v′t|ht)/Tv(v′t|h ′ t−1) 12: end if 13: until vt+1 and v′t are both accepted 14: end if 15: Sample Z3 ∼ Uniform([0, 1]n) 16: Set ht+1 = 1{Z3 ≤ σ(WTvt+1 + c)}, h′t = 1{Z3 ≤ σ(W Tv′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 6= ηt}, ξt+1 and ηt are independent in Algorithm 1 but correlated in Algorithm 3, achieved by the use of common random variatesZ2 andZ3. 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;θ) = (σ(Wh+ b)hT, σ(Wh+ 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.

heading: 5 RELATED WORK
text: 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.

heading: 6 NUMERICAL EXPERIMENTS

heading: 6.1 BARS-AND-STRIPES DATA
text: 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
1We were informed of this method by a public comment on https://openreview.net/forum?id= r1eyceSYPr, after the current paper was accepted.
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 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.

heading: 6.2 SIMULATED RBM DATA
text: 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.

heading: 6.3 FASHION-MNIST DATA
text: 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.

heading: 7 DISCUSSION
text: 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.

heading: A DISTRIBUTION COUPLING METHODS
text: 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 (ξ, η),
P (ξ = η) ≤ 1− ‖p− q‖TV = ∫ min{p(x), q(x)}dx, (4)
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(‖vt+1− v′t‖2|vt+1 6= v′t) and E(‖ht+1 − h′t‖2|vt+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(‖ht+1 − h′t‖2|vt+1,v′t). Given the v variables, ht+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(‖ht+1 − h′t‖2|vt+1,v′t) is minimized by setting ht+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 {vt+1 6= v′t}, the maximal coupling algorithm generates vt+1 and v′t with marginal distributions proportional to p1(v)−min{p1(v), p2(v)} and p2(v)−min{p1(v), p2(v)}, respectively, where p1(·) = Tv(·|ht) and p2(·) = Tv(·|h′t−1). Algorithm 3 implements this by two rejection sampling steps (line 8 and line 11). If we use the common random variates Z2 ∼ Uniform([0, 1]n) for the proposals, then vt+1 and v′t are maximally coupled in the event that they are accepted in the same iteration. This property helps to reduce E(‖vt+1 − v′t‖2|vt+1 6= v′t). To demonstrate that the specialized Algorithm 3 improves the general Algorithm 1, we sample coupled Markov chains from an RBM model withm = 500, n = 100, and the elements of (W , b, c) are all generated from aN (0, 0.12) distribution. For both algorithms, the elements of the initial state v0 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.

heading: B ADDITIONAL RESULTS FOR NUMERICAL EXPERIMENTS
text: All experiments in this article were run on an Intelr Xeonr 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
3https://www.openblas.net/ 4https://www.openmp.org/
discarded samples in the rejection sampling step. Table 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.
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
V (f̂(θ)) = E {∥∥∥f̂(θ)− E(f̂(θ))∥∥∥2
F
} ,
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. 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, andW along the training process.
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.

heading: D PROOF OF THEOREMS
text: 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
Ex∼Tθ(x|x′)l(h) ≤ γl(h ′) + γ2L1 + L2 (5)
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(ṽ)p(h̃|ṽ;θ)dy
on x,x′ ∈ D × V, where y = (ṽ, h̃). Note that q(ṽ)p(h̃|ṽ;θ) is a joint density function, so∫ q(ṽ)p(h̃|ṽ;θ)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) ≤ κ1l(h0)ρt1. As a result, we obtain
E(τ) = ∞∑ t=0 P (τ > t) ≤ κ1l(h0) 1− ρ1 . (6)
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.1 of Jacob et al. (2017), need to be verified. Under the drift condition (5) and Assumption 2, Theorem 12 of Rosenthal (1995) shows that {ξt} is a geometrically ergodic Markov chain. Therefore, there exist constants κ2 > 0 and ρ2 ∈ (0, 1) such that
|Eφ(ξt)− EMφ| ≤ κ2l(h0)ρt2 (7) for all φ : |φ(v,h)| ≤ l(h). Therefore, there exists a constant M > 0 such that |Eφ(ξt)| ≤ EMl + κ2l(h0)ρt2 ≤ D + κ2l(h0) ≤M <∞. Since |f(x;θ)|2+c ≤ l(h) by Assumption 1, we get E{|f(ξt)|2+c} ≤ M for all t > 0. Moreover, (7) implies that E{f(ξt)} → EMf . These two results verify Assumption 2.1 of Jacob et al. (2017), and hence the unbiasedness of g̃(θ) is true by design.
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 g̃2(θ) = ∑∞ t=k ∆t. From the ergodicity property (7) we immediately get E(∆2k) ≤M . For t ≥ k + 1,
E(|∆t|2+c) = E(|f(ξt)− f(ηt−1)|2+c) ≤ 21+cE(|f(ξt)|2+c + |f(ηt−1)|2+c) ≤ 22+cM. Then by Hölder’s inequality,
E(∆2t ) = E(∆2t1{τ > t}) ≤ { E(|∆t|2+c) }2/(c+2) {P (τ > t)}c/(c+2) ≤ 4M2/(c+2) { κ1l(h0)ρ t 1 }c/(c+2) .
Since g̃2(θ) = ∑∞ t=k ∆t, we get
E [ {g̃2(θ)−∆k}2 ] = E  ( ∞∑ t=k+1 ∆t )2 ≤ ∞∑ t=k+1 ∞∑ s=k+1 E |∆t∆s1{τ > t}1{τ > s})| .
We also have E |∆t∆s1{τ > t}1{τ > s})| ≤ √
E(∆2t1{τ > t})E(∆2s1{τ > s}) ≤ C0ρs3ρt3 for some constants C0 > 0 and ρ3 ∈ (0, 1), so finally,
E [ {g̃2(θ)}2 ] ≤ 2E(∆2k) + 2E [ {g̃2(θ)−∆k}2 ] ≤ 2M + 2 ∞∑ t=k+1 ∞∑ s=k+1 C0ρ 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 L1 = L2 = 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.

Review summaries:
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.

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.