RAPID COMMUNICATIONSPHYSICAL REVIEW E 91, 050101(R) (2015)

Advanced mean-field theory of the restricted Boltzmann machine
Haiping Huang* and Taro Toyoizumi
RIKEN Brain Science Institute, Wako-shi, Saitama 351-0198, Japan(Received 4 February 2015; published 18 May 2015)
Learning in restricted Boltzmann machine is typically hard due to the computation of gradients of log-likelihood
function. To describe the network state statistics of the restricted Boltzmann machine, we develop an advanced
mean-field theory based on the Bethe approximation. Our theory provides an efficient message-passing-based
method that evaluates not only the partition function (free energy) but also its gradients without requiring statisticalsampling. The results are compared with those obtained by the computationally expensive sampling-based
method.
DOI: 10.1103/PhysRevE.91.050101

PACS number(s): 02.50.Tt, 87.19.L−, 75.10.Nr

Restricted Boltzmann machine (RBM) forms building
blocks of a deep belief network [1,2], which is able to learncomplex internal representations of structured objects (suchas nature image, speech, or handwriting). RBM also has wide
applications in computational biology problems, for example,
modeling high-dimensional neural activity data from cortical
microcolumns [3].
However, learning in RBM is computationally difficult,since gradients of the log-likelihood function need to becomputed at every iteration step to update the model parameters. This computation is usually accomplished by the
Gibbs-sampling-based method or its variants [4,5], for which
the tradeoff between accuracy and convergence speed requires
careful considerations. Furthermore, an efficient way to evaluate the partition function (e.g., log-likelihood function forcross-validation analysis) remains unknown.
Here, we develop a mean-field theory for the RBM basedon the cavity method (Bethe approximation) [6], whichyields an efficient and fully distributed algorithm to evaluate
the free energy (partition function) of a RBM of interest.
The remarkable efficiency is confirmed by comparing thecomputation results of gradients of log-likelihood function by
Gibbs sampling and the proposed mean-field theory.
A RBM [7,8] consists of one hidden layer and one visible
layer without lateral connections between nodes in each layer.
We assume the hidden layer has M nodes, while the visible
layer has N nodes. Hidden node a with external field ha isconnected to visible node j with field φj by a symmetriccoupling w
function forRBM is thus defined
aj . The energy 
by E = − i,a σi wai sa − i σi φi − a ha sa , where σi andsa are used to specify the state of visible node i and hidden
node a, respectively. Due to the conditional independence ofhidden nodes’ state given σ , the state of the hidden nodes canbe marginalized. This leads to the following probability of avisible state:1 P (σ ) =[2 cosh(wa σ + ha )]eφi σi ,(1)
Z a
iwhere wa denotes the ath row vector of the coupling matrix w.
Z is a normalization constant (also called the partition
function) of the model. As a model study, we assume theelement of the matrix w is independently and identically*

physhuang@gmail.com; sites.google.com/site/physhuang

1539-3755/2015/91(5)/050101(5)distributed with a normal distribution with mean zero andvariance g/N . We assume that the external field for both layers
follows a normal distribution with mean zero and variance v.
We denote the ratio between the number of hidden nodes andthat of visible nodes by α = M/N, where M and N canbe arbitrarily large. A schematic representation of a RBM(M = 3, N = 5) is shown in Fig. 1.
An exact computation of Z requires an exponential
computational complexity (2N ), which becomes impossible
for a relatively large N . However, advanced mean-fieldapproximation can be used to compute approximate values
under certain condition, and its prediction should be comparedwith numerical simulations. Here, we propose the Bethe
approximation [6] to tackle this problem. First, we transformthe original model (left panel of Fig. 1) into a factor graph(right panel of Fig. 1) [9], where each square node indicatesa Boltzmann factor 2 cosh(wa σ + ha ) in Eq. (1). Then, we
introduce the cavity probability Pi→a (σi ) that the visible
node i takes state σi in the absence of the contribution
from factor node a [10], and Pi→a (σi ) satisfies the following
self-consistent equations:μb→i (σi ),(2a)Pi→a (σi ) ∝ eφi σib∈∂i\aμb→i (σi ) =

{σj |j ∈∂b\i}

2 cosh (w b · σ + hb )

Pj →b (σj ),j ∈∂b\i(2b)where the symbol ∝ indicates a normalization constant, ∂i\adenotes the neighbors of node i except factor node a, ∂b\i
denotes the neighbors of factor node b except visible node i,and the auxiliary quantity μb→i (σi ) represents the contribution
from factor node b to visible node i given the value of σi [10].
With these definitions, the products in Eq. (2) are reasonable
under the weak correlations assumption, whereas, the validityof this Bethe approximation should be checked by a stability
analysis.
Note that the computation in Eq. (2b) is still intractable
due to the summation
 over all possible σ except σi . However,
because Ub→i ≡ j ∈∂b\i wbj σj is a sum of a large number ofnearly independent random variables, the central-limit theorem implies that the distribution of Ub→i is well
 characterized
by its mean and variance [11], i.e., Gb→i = j ∈∂b\i wbj mj →b2and 2b→i  j ∈∂b\i wbj(1 − m2j →b ), respectively, where

050101-1©2015 American Physical Society

RAPID COMMUNICATIONSHAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 91, 050101(R) (2015)

FIG. 1. (Color online) A RBM is composed of one hidden layer
and one visible layer. No lateral connections exist within both hidden
and visible layers. Each hidden node is connected to all visible nodes
with symmetric coupling weights and is responsible for capturinghigh-order dependence. The original RBM is shown in the left panel as
an example of three hidden nodes (solid circles) and five visible nodes(empty circles). The right panel shows a transformed factor graphafter marginalization of hidden states for theoretical analysis. Eachfactor node (square node) represents the probabilistic normalizationof a hidden node given the state of all visible nodes (see the maintext).

random values, and then iterate Eq. (3) until it converges withina prescribed accuracy. Note that the overall time complexity
is of the order O(N 2 ); furthermore, the algorithm is fully
distributed and thus amenable to large-scale applications.
In the remaining part, we demonstrate the computationof the free energy on large single instances by applying themessage passing algorithm and confirm the accuracy of theresults by comparing gradients of the log-likelihood with thoseobtained by the Gibbs sampling method. A stability analysisof the message-passing algorithm is also presented.
We run the message-passing equations on single instancesof RBMs as size N , hidden-node density α, and couplingstrength g are varied. As displayed in Fig. 2(a), the free-energy
density decreases as α increases. Note that the density doesmj →b ≡ σj σj Pj →b (σj ) denotes the cavity magnetization(the average of σj in the absence of factor node b).
Because we consider the binary spin variable σi = ±1,
Pi→a (σi ) and μb→i (σi ) can be parametrized by mi→a and cavb→i (σi =1)
ity bias ub→i , respectively. ub→i is defined as 12 ln μμb→i.
(σi =−1)
The practical recursive equations, the so-called message
passing equations, are thus derived as⎛⎞mi→a = tanh ⎝φi +ub→i ⎠ ,(3a)b∈∂i\aub→i =

1 cosh(hb + Gb→i + wbi )ln,
2 cosh(hb + Gb→i − wbi )(3b)where the b→i dependency in Eq. (3b) drops because of thesymmetry of cosh. The cavity magnetization can be understoodas the message passing from the visible node to the factor node,
while the cavity bias is interpreted as the message passing
from the factor node to the visible node. This message passing
based computation is much more accurate than naive meanfield approximation [12], which assumes a fully factorized
distribution for Eq. (1). In contrast, Eq. (3) captures nearest
neighbors’ correlations.
Once the iteration of Eq. (3) converges, the free energy ofthe model can be computed from the fixed-point solution.
Under the Bethe approximation, the Bethe free energy isexpressed as [10,11]:F =−ln Zi + (N − 1)ln Za ,
(4)
i

a
−φi

where Zi = eb∈∂i μb→i (+1) + eb∈∂i μb→i (−1),
2b→i /2
cosh(hb + Gb→i + wbi σi ).
in which μb→i (σi ) = 2e22Za = 2ea /2 cosh(ha + G
a are given bya ). Ga2 and 
2
w
m
and
w
(1
−
m
j ∈∂a aj j →a
j ∈∂a aj
j →a ), respectively.
Each stable solution of the message passing algorithm ofEq. (3) corresponds to a local minimum of the free-energyfunction in Eq. (4) [13].
RBM defined in Fig. 1 is basically a densely connected
graphical model. Our mean-field theory provides a practical
way to estimate the free energy of single instances (typicalexamples of the model). More precisely, we initialize the cavity
magnetization and bias on each link of the factor graph by
φiFIG. 2. (Color online) Free-energy density (f = F /N ) of single
instances of RBMs. Iterations of the message passing equations
are always converged to produce the data points. The error bars
give statistical errors across ten random instances. (a) Free-energy
density as a function of α (density of hidden nodes); g = 1.0 andv = 0.05. The inset shows the absolute difference, |f |, of thefree-energy density estimated by the exact enumeration and theBethe approximation (BA) for N = 20 (see a comparison for ten
instances with α = 0.6 in the inset, where the line indicates equality).
(b) Free-energy density as a function of weight strength g and fieldstrength v (shown in the inset). α = 0.5.

050101-2

RAPID COMMUNICATIONSADVANCED MEAN-FIELD THEORY OF THE RESTRICTED . . .

PHYSICAL REVIEW E 91, 050101(R) (2015)

not change significantly at two large sizes (N = 100 and
N = 1000). Furthermore, the inset of Fig. 2(a) shows that thetheoretical result even matches well with the exact enumeration
result for small size N = 20. As the variance parameter gof weights increases, the free-energy density also decreases[Fig. 2(b)]. The same property also holds when the variancev of external field increases [the inset of Fig. 2(b)]. In the
explored range of g (or v) and α, Eq. (3) converges in a
few steps to a single fixed point on which the free energy
is calculated. Therefore, the Bethe approximation provides
an accurate estimation of free energy much faster than othersampling based procedures, which are typically slow to reachan equilibrium state.
The stability of Eq. (3) can also be studied on single
instances. Apart from the cavity magnetization, we introduce
its variance as an extra message denoted by Vi→a [14,15]. The
evolution of Vi→a follows:(tmax = 500). λ grows with α and g, and the fluctuation acrossinstances becomes strong near to the critical point (λ = 1).
Note that increasing α has an equivalent effect of increasingg. In the right inset, two typical examples are shown. Near
to the critical point, some instances have decaying variancestrength (λ < 1), while some have growing strength (λ > 1).
Equation (5) thus tells us how stable (unstable) the iteration of
Eq. (3) is for a particular RBM. The algorithm converges in a
few iteration steps to a solution unless the recursive process is
close to the instability boundary.
RBM can be used to model the real data, and the parameters
are fitted to maximize the probability of observing thetraining data [2]. This would lead to computation of thefollowing quantities, mi ≡ σi , m̂a ≡ tanh(wa σ + ha ), andCaj ≡ tanh(wa σ + ha )σj , where the average, ·, is takenover the distribution defined in Eq. (1), which is intractablewithout approximations. Here, mi is the average of visible
state σi , m̂a is the average of hidden state sa , and Caj is
the correlation between sa and σi . Although an accurate
evaluation of the above quantities requires a sufficiently
long Gibbs sampling (so-called the most difficult negative
phase in the machine-learning community [2]), we cancompute them by the message passing equations and thencompare the results with those obtained by Gibbs sampling tocheck the consistency of the theory.
Following the same spirit, our theory gives the theoreticalevaluation of the above quantities as

Vi→a =

1 − m2i→a
4

2



Pb→i

b∈∂i\a× [tanh(b→i ) − tanh(b→i − 2wbi )]2 ,(5)

where b→i ≡ hb + Gb→i + wbi and Pb→i ≡2. The stability is measured by the total
j ∈∂b\i wbj Vj →bvariance S(t) = (i,a) Vi→a (t) summed over all connectedpairs (i,a), where t is the iteration step. The explosion ofS(t) indicates the instability of the message passing equations,which is related to the divergence of the (nonlinear) spin-glasssusceptibility [16], and thus the Bethe approximation becomesinconsistent. We study this effect on single instances of
RBM as shown in Fig. 3. A relative strength is denoted as
λ = S(t + 1)/S(t), where t denotes the step at which theiteration converges or exceeds a prefixed maximal number

FIG. 3. (Color online) Stability parameter λ as a function ofmodel parameters. N = 1000, α = 0.5, and v = 0.05. The error bars
give statistical errors across ten random instances. The left inset gives
stability versus α with g = 1.0 and v = 0.05. The right inset gives
two examples (black and red) taken from the main figure at g = 2.1.
S(0) is the initial total variance.

mi = tanh φi +

ub→i ,(6a)˜ a x + G̃a ),Dx tanh((6b)Caj  m̂a mj + waj 1 − m2j Aa ,

(6c)


m̂a =

b∈∂i

√
2
where
Dx ≡ e−x /2 / 2π dx is aGaussian measure. Aa ≡ 1 −
2 ˜
˜ 2a 
(a x + G̃a ), G̃a = k∈∂a wak mk + ha , and 
Dx tanh22k∈∂a wak (1 − mk ) [17]. Equation (6) is computed based on
the fixed point of the iterative algorithm [Eq. (3)].
We used alternating Gibbs sampling [5] to evaluate theequilibrium average of the gradients. More precisely, thehidden nodes are updated in parallel according to P (sa =1|σ ) = ewa σ +ha /[2 cosh(wa σ + ha )], while the visible nodes
are then all updated in parallel according to P (σi = 1|s) =Tewi s+φi /[2 cosh(wTi s + φi )], where wi is the ith column of theweight matrix. Note that the visible nodes are conditionallyindependent given the hidden states and vise versa [5]. These
two steps of updates form one full step of the alternating Gibbs
sampling. If this Markov chain is run for a sufficiently longtime, the stationary (equilibrium) distribution is expected tobe reached, from which the averages can be estimated. We
test our theory in a system with N = 100 visible nodes andrun the Markov chain with 106 steps for thermal equilibration and the other 4 × 106 steps to collect a total numberof 105 samples to calculate the average. We measure theperformance by the root-mean-square (RMS) error between
the Gibbs sampling (GS) result and the Bethe approximation(BA) result, which is shown in Fig. 4. The RMS error is050101-3

RAPID COMMUNICATIONSHAIPING HUANG AND TARO TOYOIZUMIPHYSICAL REVIEW E 91, 050101(R) (2015)

FIG. 4. (Color online) Evaluation performance of mean-field
theory in comparison with the Gibbs sampling. Iterations of therecursive equations always converged to produce the data points.
The error bars give statistical errors across ten random instances.
(a) RMS error as a function of g; N = 100 and v = 0.02.
(b) Scatter plot for a typical example of N = 100, α = 0.5, g = 0.1,
and v = 0.02. The inset shows an example of g = 0.55 (otherparameters do not change). The line indicates equality. (c) RMS
error δm reached by CD-k as a function of k in comparison to theBethe approximation; N = 100, α = 0.2, g = 0.55, and v = 0.02.
The result is averaged over five random instances.

 
|Y |
defined as δY ≡ |Y1 | i=1 (YiGS − YiBA )2 , where Y takes m,m̂, or C, and |Y | indicates the number of these parameters.
Small RMS error indicates that the Bethe approximation isaccurate.
As shown in Fig. 4(a), all evaluation errors grow with theweight strength g, which is reasonable since our mean-field
theory will break down when the network enters a strongly
correlated state, as already shown by the stability analysis.
In a similar manner, the error grows with the hidden nodedensity α, because each hidden node puts a constraint to
the network and all constraints compete with each other togive an equilibrium state, resulting in strong correlations
with high α. However, the magnitude of all errors is small,
implying that one can acquire accurate estimation of gradientsof log-likelihood function by passing messages on a factor
graph as well. We show this point more clearly with ascatter plot in Fig. 4(b) (Bethe approximation result versus
Gibbs sampling result). This accuracy is obtained by requiringmuch fewer computational costs compared with the Gibbs
sampling. The comparison further confirms the efficiency ofthe proposed mean-field method across a wide range of model
parameters.
Note that the Gibbs sampling result serves as the groundtruth here, since we run the Markov chain for a long time.
More practically, one can estimate the statistics by k-steps
contrastive divergence (CD-k) algorithm [5] that requires
a time complexity of O(kT MN ), where T denotes thenumber of sample particles. However, to reach a similar
accuracy as the Bethe approximation, it typically requiresk  10 and T ∼ 105 under the current setup [Fig. 4(c)]. In
contrast, the Bethe approximation yields a time complexityof O(nMN ) with n < 100, where n is the number of
iterations and one iteration involves the update of MN cavity
messages.
In conclusion, we propose a mean-field theory for the RBM,
a widely used model in the machine-learning communityand biological data analysis. The theory captures nearest
neighbors’ correlations by operating on the cavity factor graph(by removing factor nodes), leading to an approximate estimation of the free-energy function (log-likelihood function) for
single instances of large-size networks, for which the standard
Gibbs sampling procedure becomes prohibitively slow to
get a reliable result (e.g., for evaluating the log-likelihood
function for cross-validation analysis). Moreover, we replace
the normal Gibbs sampling with a mean-field computation
based on message passing algorithm, to estimate the gradientsof log-likelihood function and show its efficiency by extensive
numerical simulations on single instances. The nature of thisfast inference lies in the fact that the information is exchanged
locally between factor nodes and visible nodes, to reacha coherent fixed point, which may provide a computation
paradigm for probabilistic inference in neural networks. We
expect the mean-field theory inspired calculation will be usefulin practical applications and bring more insights to understandthe RBM and its role in deep learning [18].
This work was supported by RIKEN Brain Science Instituteand the Brain Mapping by Integrated Neurotechnologies for
Disease Studies (Brain/MINDS) by the Ministry of Education,
Culture, Sports, Science and Technology of Japan (MEXT).

050101-4

RAPID COMMUNICATIONSADVANCED MEAN-FIELD THEORY OF THE RESTRICTED . . .

PHYSICAL REVIEW E 91, 050101(R) (2015)[1] G. E. Hinton and R. R. Salakhutdinov, Science 313, 504 (2006).
[2] Y. Bengio, A. Courville, and P. Vincent, IEEE Trans. PatternAnal. Machine Intell. 35, 1798 (2013).
[3] U. Köster, J. Sohl-Dickstein, C. M. Gray, and B. A. Olshausen,
PLoS Comput. Biol. 10, e1003684 (2014).
[4] G. Hinton, Neural Comput. 14, 1771 (2002).
[5] G. Hinton, S. Osindero, and Y. Teh, Neural Comput. 18, 1527
(2006).
[6] M. Mézard and G. Parisi, Eur. Phys. J. B 20, 217 (2001).
[7] P. Smolensky, Information Processing in Dynamical Systems:
Foundations of Harmony Theory (MIT Press, Cambridge, MA,
1986), pp. 194–281.
[8] Y. Freund and D. Haussler, Technical Report (Santa Cruz, CA,
1994).
[9] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, IEEE Trans.
Inf. Theory 47, 498 (2001).

[10] M. Mézard and A. Montanari, Information, Physics,
and Computation (Oxford University Press, Oxford,2009).
[11] H. Huang, K. Y. M. Wong, and Y. Kabashima, J. Phys. A: Math.
Theor. 46, 375002 (2013).
[12] H. J. Kappen and F. B. Rodriguez, Neural Comput. 10, 1137
(1998).
[13] T. Heskes, Neural Comput. 16, 2379 (2004).
[14] Y. Kabashima, J. Phys. A 36, 11111 (2003).
[15] A. Montanari and F. Ricci-Tersenghi, Phys. Rev. B 70, 134406(2004).
[16] O. Rivoire, G. Biroli, O. Martin, and M. Mézard, Eur. Phys. J.
B 37, 55 (2004).
[17] H. Huang and Y. Kabashima, J. Stat. Mech.: Theory Exp. (2014)
P05020.
[18] P. Mehta and D. J. Schwab, arXiv:1410.3831.

050101-5