Entropy, Free Energy, and Work of Restricted Boltzmann Machines

A restricted Boltzmann machine is a generative probabilistic graphic network. A probability of finding the network in a certain configuration is given by the Boltzmann distribution. Given training data, its learning is done by optimizing the parameters of the energy function of the network. In this paper, we analyze the training process of the restricted Boltzmann machine in the context of statistical physics. As an illustration, for small size bar-and-stripe patterns, we calculate thermodynamic quantities such as entropy, free energy, and internal energy as a function of the training epoch. We demonstrate the growth of the correlation between the visible and hidden layers via the subadditivity of entropies as the training proceeds. Using the Monte-Carlo simulation of trajectories of the visible and hidden vectors in the configuration space, we also calculate the distribution of the work done on the restricted Boltzmann machine by switching the parameters of the energy function. We discuss the Jarzynski equality which connects the path average of the exponential function of the work and the difference in free energies before and after training.


I. INTRODUCTION
A restricted Boltzmann machine (RBM) [1] is a generative probabilistic neural network.RBMs and general Boltzmann machines are described by a probability distribution with parameters, i.e., the Boltzmann distribution.An RBM is an undirected Markov random fields, and is considered a basic building block of deep neural network.RBMs have been applied widely, for example, dimensional reduction, classification, feature learning, pattern recognition, topic modeling, and so on [2][3][4].
As its name implies, a RBM is closely connected to physics, and they shares some important concepts such as entropy, free energy, etc. [5].Recently, RBMs have renewed much attention in physics since Carleo and Troyer [6] showed that a quantum many-body state could be efficiently represented by the RBM.Gabré et al. and Tramel et al. [7] employed the Thouless-Anderson-Palmer mean-field approximation, used for a spin glass problem, to replace the Gibbs sampling of contrast-divergence training.Amin et al. [8] proposed a quantum Boltzmann machine based on quantum Boltzmann distribution of a quantum Hamiltonian.More interestingly, there is a deep connection between Boltzmann machine and tensor networks of quantum many-body systems [9][10][11][12][13].Xia and Kais combined the restricted Boltzmann machine and quantum algorithms to calculate the electronic energy of small molecules [14].
While the working principles of RBMs have been well established, it may be still needed to understand the RBM better for further applications.In this paper, we investigate the RBM from the perspective of statistical physics.As an illustration, for bar-and-stripe pattern data, the thermodynamic quantities such as the entropy, the internal energy, the free energy, and the work, are calculated as a function of epoch.Since the RBM is a bipartite system composed of visible and hidden layers, it may be interesting, and informative, to see how the correlation between the two layers grows as the training goes on.We show that the total entropy of the RBM is always less than the sum of the entropies of visible and hidden layers, except at the initial time when the training begins.This is the so-called subadditivity of entropy, indicating that the visible layer becomes correlated with the hidden layer as the training proceeds.The training of the RBM is to adjust the parameters of the energy function, which can be considered as the work done on the RBM, from a thermodynamic point of view.Using the Monte-Carlo simulation of the trajectories of the visible and hidden vectors in configuration space, we calculate the work of a single trajectory and the statistics of the work over the ensemble of trajectories.We also examine the Jarzynski equality that connects the ensemble of the work done on the RBM and the difference in free energies before and after training of the RBM.
The paper is organized as follows.In Section II, the detailed analysis of the RBM from the statistical physics point of view is described.In Section III, we presents the summary of the result together with discussions.
Let us start with a brief introduction of the RBM [1][2][3].As shown in Fig. (1), the RBM is composed of two layers; the visible layer and the hidden layer.Possible configurations of the visible and hidden layers are represented by the random binary vectors, v = (v 1 , . . ., v N ) ∈ {0, 1} N and h = (h 1 , . . ., h M ) ∈ {0, 1} M , respectively.The interaction between the visible and hidden layers is given by the so-called weight matrix w ∈ R N × R M where the weight w i j is the connection strength between a visible unit v i and a hidden unit arXiv:2004.04900v1[cond-mat.dis-nn]10 Apr 2020

Visible layer
FIG. 1: Graph structure of a restricted Boltzmann machine with the visible layer and the hidden layer.
h j .The biases b i ∈ R. and c j ∈ R are applied to visible unit i and hidden unit j, respectively.Given random vectors v and h, the energy function of the RBM is written as an Ising-type Hamiltonian where the set of model parameters is denoted by θ ≡ {w i j , b i , c j }.The joint probability of finding v and h of the RBM is given by the Boltzmann distribution where the partition function, Z(θ) ≡ v,h e −E(v,h;θ) , is the sum over all possible configurations.The marginal probabilities p(v; θ) and p(h; θ) for visible and hidden layers are obtained by summing up the hidden or visible variables, respectively, The training of the RBM is to adjust the model parameter θ such that the marginal probability of the visible layer p(v; θ) becomes as close as possible to the unknown probability p data (v) that generate the training data.Given identically and independently sampled training data D ∈ {v (1) , . . ., v (D) }, the optimal model parameters θ can be obtained by maximizing the likelihood function of the parameters, L(θ|D) = D i=1 p(v (i) ; θ), or equivalently by maximizing the log-likelihood function ln L(θ|D) = D i=1 ln p(v (i) ; θ).Maximizing the likelihood function is equivalent to minimizing the Kullback-Leibler divergence or the relative entropy of p(v; θ) from q(v) [15,16] where q(v) is an unknown probability that generates the training data.Another method of monitoring the progress of training is the cross-entropy cost between the input visible vector v (i) and a reconstructed visible vector v(i) of the RBM, The stochastic gradient ascent method for the log-likelihood function is used to train the RBM.Estimating the loglikelihood function requires the Monte-Carlo sampling for the model probability distribution.Well-known sampling methods are the contrast-divergence, denoted by CD-k, and the persistent contrast divergence PCD-k.For the detail of the RBM algorithm, please see Refs.[2][3][4].Here we employ the CD-k method.
B. Free energy, entropy, and internal energy From physics point of view, the RBM is a finite classical system composed of two subsystems, similar to an Ising spin system.The training of the RBM is considered the driving of the system from an initial equilibrium state to the target equilibrium state by switching the model parameters.It may be interesting to see how thermodynamic quantities such as free energy, entropy, internal energy, and work, change as the training progresses.
It is straightforward to write down various thermodynamic quantities for the total system.The free energy F is given by the logarithm of the partition function Z, The internal energy U is given by the expectation value of the energy function E(v, h; θ) The entropy S of the total system comprising the hidden and visible layers is given by Here, the convention of 0 ln 0 = 1 is employed if p(v, h) = 0.The free energy ( 6) is related to the difference between the internal energy (8) and the entropy ( 9) where T is set to 1. Generally, it is very challenging to calculate the thermodynamic quantities, even numerically.The number of possible configurations of N visible units and M hidden units grow exponentially as 2 N+M .Here, for a feasible benchmark test, the 2×2 Bar-and-Stripe data are considered [17,18].Fig. 2 shows the 6 possible 2 × 2 Bar-and-Stripe patterns out of 16 possible configurations, which will be used as the training data in this work.We take the sizes of the visible and the hidden layers as N = 4 and M = 6, respectively.In order to understand better how the RBM is trained, the thermodynamic quantities are calculated numerically for this small benchmark system.Fig. 3 shows how the weight w i j , the bias b i on the visible unit i and the bias c j on the hidden unit j change as the training goes on.The weight w i j are clustered into 3 classes.The evolution of the bias b i on the visible layer is somewhat different from that of the bias c j on the hidden layer.The change in c i are larger than that in b i .Fig. 4 shows the change in the marginal probabilities p(v) of the visible layer and p(h) of the hidden layer before and after training.Note that the marginal probability p(v) after training is not distributed exclusively over 6 possible outcomes according to the training data set in Fig. 2.
Typically, the progress of learning of the RBM is moni- tored by the loss function.Here the Kullback-Leibler divergence, Eq. ( 4) and the reconstructed cross entropy, Eq. ( 5) are are used.Fig. 5 plots the reconstructed cross entropy C, the Kullback-Leibler divergence D KL , the entropy S , the free energy F, and the internal energy U as a function of epoch.As shown in Fig. 5 (a), it is interesting to see that even after a large number of epochs ∼ 10, 000, the cost function C continues approaching zero while the entropy S and the Kullback-Leibler divergence D KL become steady.On the other hand, the free energy F continues decreasing together with the internal energy U, as depicted in Fig. 5 (b).The Kullback-Leibler divergence is a well-known indicator to the performance of RBMs.Then, our result implies that the entropy may be another good indicator to monitor the progress of RBM while other thermodynamic quantities may be not.In addition to the thermodynamic quantities of the total system of the RBM, Eqs. ( 6), (8), and (7), it is interesting to see how the two subsystems of the RBM evolve.Since the RBM has no intra-layer connection, the correlation between the visible layer and the hidden layer may increase as the training proceeds.The correlation between the visible layer and the hidden layer can be measured by the difference between the total entropy and the sum of the entropies of the two subsystems.The entropies of the visible and hidden layers are given by The entropy S V of the visible layer is closely related to the Kullback-Leibler divergence of p(v; θ) to an unknown probability q(v) which produces the data.Eq. ( 4) is expanded as The second term − v q(v) ln p(v; θ) depends on the parameter θ.As the training proceeds, p(v; θ) becomes close to q(v) so the behavior of the second term is very similar to that of the entropy S V of the visible layer.If the training is perfect, we have q(v) = p(v; θ) that leads to D KL (q || p) = 0 while S V remains nonzero.The difference between the total entropy and the sum of the entropies of subsystems is written as Eq. ( 12) tells that if the visible random vector v and the hidden random vector h are independent, i.e., p(v, h; θ) = p(v; θ)p(h; θ), then the entropy S of the total system is the sum of the entropies of subsystems.In general, the entropy S of the total system is always less than or equal to the sum of the entropy of the visible layer, S V , and the entropy of the hidden layer, S H , [19], This is called the subadditivity of entropy, one of the basic properties of the Shannon entropy, which is also valid for the von Neumann entropy [20,21].This property can be proved using the log inequality, − ln x ≥ −x+1.In other way, Eq. ( 13) may be proved by using the log-sum inequality, which states that for the two sets of nonnegative numbers, a 1 , . . ., a n and b In other words, Eq. ( 12) can be regarded as the negative of the relative entropy or Kullback-Leibler divergence of the joint probability p(v, h) to the product probability p(v) • p(h), For the 2×2 Bar-and-Stripe pattern, the entropies of visible and hidden layers, S V , S H are calculated numerically.Fig. 6 plots the entropies, S V , S H , S , and the Kullback-Leibler divergence D KL (q || p) as a function of epoch.Fig. 6 (a) shows that the Kullback-Leibler divergence, D KL (q || p) becomes saturated, though above zero, as the training proceeds.Similarly, the entropy S V of the visible layer is saturated.This implies that the entropy of the visible layer, as well as the total entropy shown in Fig. 5, can be an indicator to learning better than the reconstructed cross entropy C, Eq. ( 5).The same can also be said about the entropy of the hidden layer, S H .
The difference between the total entropy and the sum of the entropies of the two subsystems, S − (S V + S H ), becomes less than 0, as shown in Fig. 6 (b).Thus it demonstrates the subadditivity of entropy, i.e., the correlation between the visible and the hidden layer as the training proceeds.As it is saturated just as the total entropy and the entropies of the visible and hidden layers after a large number of epoch, the correlation between the visible layer and the hidden layer can also be a good quantifier of the RBM progress.

C. Work, free energy, and Jarzynski equality
The training of the RBM may be viewed as driving a finite classical spin system from an initial equilibrium state to a final equilibrium state by changing the system parameters θ slowly.If the parameters θ are switched infinitely slowly, the classical system remains in quasi-static equilibrium.In this case, the total work done on the systems is equal to the Helmholtz free energy difference between the before-training and the aftertraining, W ∞ = F 1 − F 0 .For switching θ at a finite rate, the system may not evolve immediately to an equilibrium state, the work done on the system depends on a specific path of the system in the configuration space.Jarzynski [22,23] proved that for any switching rate the free energy difference ∆F is related to the average of the exponential function of the amount of work W over the paths The RBM is trained by changing the parameters θ through a sequence {θ 0 , θ 1 , . . ., θ τ }, as shown in Fig. 3. To calculate the work done during the training, we perform the Monte Carlo simulation of the trajectory of a state (v, h) of the RBM in configuration space.From the initial configuration, (v 0 , h 0 ) which is sampled from the initial Boltzmann distribution, Eq. ( 2), the trajectory (v 0 , h 0 ) → (v 1 , h 1 ) → • • • → (v τ , h τ ) is obtained using the Metropolis-Hastings algorithm of the Markov chain Monte-Carlo method [24,25].Assuming the evolution is Markovian, the probability of taking a specific trajectory is the product of the transition probabilities at each step, The transition (v, h) → (v , h ) can be implemented by the Metropolis-Hastings algorithm based on the detailed balance condition for the fixed parameter θ, The work done on the RBM at epoch i may be given by The total work W = δW i performed on the system is written as [26]  The vertical and the horizontal axes represent each configuration of the visible and the hidden layers, respectively.The black tiles represent the lowest energy configurations among all configurations, thus the probability of finding that configuration is high.
Given the sequence of the model parameter {θ 0 , θ 1 , . . ., θ τ }, the Markov evolution of the visible and hidden vectors (v, h) ∈ {0, 1} N+M may be considered the discrete random walk.Random walkers move to the points with low energy in configuration space.Fig. 7 shows the heat map of energy function E(v, h; θ) of the RBM for the 2 × 2 Bar-and-Stripe pattern after training.One can see the energy function has deep levels at the visible vectors corresponding to the Bar-and-Stripe patterns of the training data set in Fig. 2, representing probability of generating the trained patterns is high.Also note that the energy function has many local minima.Fig. 8 plots a few Monte-Carlo trajectories of the visible vector v as a function of epoch.Before training, the visible vector v is distributed over all possible configurations, represented by the  (a).Since the work done on the system depends on the path, the distribution of the work is calculated by generating many trajectories.Fig. 9 shows the distribution of the work over 50000 paths at 5000 training epoch.The Monte-Carlo average of the work is W ≈ −5.481, and its standard deviation is σ W ≈ 3.358.The distribution of the work generated by the Monte-Carlo simulation is well fitted to the Gaussian distribution, as depicted by the red curve in Fig. 9.This agrees with the statement of in Ref. [23] that for the slow switching of the model parameters the probability distribution of work is approximated to the Gaussian.We perform the Monte-Carlo calculation of the exponential average of work, e −W path to check the Jarzynski equality, Eq. ( 16).The free energy difference can be estimated as where N mc is the number of the Monte-Carlo samplings.At small epoch number, the Monte-Carlo estimated value of free energy difference is close to ∆F calculated from the partition function.However, this Monte-Carlo calculation gives rise to the poor estimation of the free energy difference if the epoch is greater than 5000.This numerical errors can be explained by the fact that the exponential average of the work is dominated by rare realization [27][28][29][30][31].As shown in Fig. 9, the distribution of work is given by the Gaussian distribution ρ(W) with the mean W and the standard deviation σ W .If the standard deviation σ W becomes larger, the peak position of ρ(W)e −W moves to the long tail of the Gaussian distribution.So the main contrition of the integration of e −W comes from the rare realizations.Fig. 10 shows that the standard deviation σ W grows with epoch, so the error of the Monte-Carlo estimation of the exponential average of the work grows quickly.If σ 2 W k B T , the free energy is related to the average of work and its variance as Here, the case is opposite, the spread of the value of work is large, i.e., σ 2 W k B T (= 1), so the central limit theorem does not work and the above equation can not be applied [32].

III. SUMMARY
In summary, we analyzed the training process of the RBM in the context of statistical physics.In addition to the typical loss function, i.e., the reconstructed cross entropy, the thermodynamic quantities such as free energy F, internal energy U, and entropy S were calculated as a function of epoch.While the free energy and the internal energy decrease rather indefinitely with epochs, the total entropy and the entropies of the visible and the hidden layers become saturated together with the Kullback-Leibler divergence after a sufficient number of epochs.This result suggests that the entropy of the system may be a good indicator to the RBM progress along with Kullback-Leibler divergence.It seems worth investigating the entropy for other larger data sets, for example, MNIST handwritten digits [33], in future works.
We have further demonstrated the subadditivity of the entropy, i.e., the entropy of the total system is less than the sum of the entropies of the two layers.This manifested the correlation between the visible and hidden layers growing with the training progress.Just as the entropies are well saturated together with Kullback-Leibler divergence, so is the correlation that is determined by the total and the local entropies.In this sense, the correlation between the visible and the hidden layer may become another good indicator to the RBM performance.
We also investigated the work done on the RBM by switching the parameters of the energy function.The trajectories of the visible and hidden vectors in configuration space were generated using Markov chain Monte-Carlo simulation.The distribution of the work follows the Gaussian distribution and its standard deviation grows with training epochs.We discussed the Jarzynski equality, which connects the free energy difference and the average of the exponential function of the work over the trajectories.
A more detailed analysis from a full thermodynamics or statistical physics point of view can bring us useful insights into the performance of RBM.This course of study may enable us to come up with possible methods for a better performance of RBM for many different applications in the long run.Therefore, it may be worthwhile to further pursue our study, e.g. a rigorous assessment of scaling behavior of thermodynamic quantities with respect to epochs as the sizes of the visible and hidden layers increase.We also expect that a similar analysis on a quantum Boltzmann machine can be valuable as well.

FIG. 3 :
FIG. 3: (a) Bias b i on the visible unit i and bias c j on the hidden unit j are plotted as a function of epoch.(b) Weight w i j connecting the visible unit i and the hidden unit j are plotted as a function of epoch.

FIG. 4 :
FIG. 4: Marginal probabilities p(v) of visible layer and p(h) of hidden layer are plotted (a) before training and (b) after training.The binary vector v or h in x-axis is represented by the decimal number as noted in the caption of Fig. 2. The visible and the hidden layers have total number of configurations given by 2 4 = 16 and 2 6 = 64, respectively.The learning rate is 0.15, the training epoch 20000, and CD-k 100.

FIG. 5 :
FIG. 5: For 2 × 2 bar-and-stripe data, (a) cost function C, entropy S , and the Kullback-Leibler divergence D KL (q || p) are plotted as a function of epoch.(b) Free energy F, entropy S , and internal energy U of the RBM are calculated as a function of epoch.
FIG. 6: (a) Kullback-Leibler divergence D KL (q || p), entropy S V , and their difference are plotted as a function of epoch.(b) Entropy S of the total system, entropy S V of the visible layer, entropy S H of the hidden layer, and the difference S − S H − S V are plotted as a function of epoch.

10 FIG. 7 :
FIG. 7: Heat map of energy function E(v, h; θ), representing the energy level of each configuration, after training of 2 × 2 Bar-and-Stripe patterns for 50000 epochs.The sizes of visible and hidden layers are N = 4 and M = 6, respectively.The learning rate is r = 0.15 and the value of CD k is k = 100.The vertical and the horizontal axes represent each configuration of the visible and the hidden layers, respectively.The black tiles represent the lowest energy configurations among all configurations, thus the probability of finding that configuration is high.

FIG. 8 :FIG. 9 :
FIG. 8: Markov chain Monte-Carlo trajectories of the visible vector v i are plotted as a function of epoch.The visible vector jumps frequently in the early state of training and becomes trapped into one of target states as the training proceeds.

FIG. 10 :
FIG. 10: Average of work done with standard deviation and free energy difference ∆F = F(epoch) − F(epoch = 0) as a function of epoch.The error bar of the work represent the standard deviation of the Gaussian distribution.

Fig. 10
shows how the average of work, W path , over the Markov chain Monte-Carlo paths changes as a function of epoch.The standard deviation of the Gaussian distribution of the work also grows as a function of training epoch.The free energy difference between before-training and after-training is called the reversible work W r = ∆F.The difference between the actual work and the reversible work is called the dissipative work, W d = W − W r [26].As depicted in Fig. 10, the magnitude of the dissipative work grows with training epoch.