Entropy, Free Energy, and Work of Restricted Boltzmann Machines

Let us start with a brief introduction of the RBM Smolensky (1986); Hinton (2012); Fischer and Igel (2014). 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, $\mathbf{v}=(v_{1},\dots,v_{N})\in\{0,1\}^{N}$ and $\mathbf{h}=(h_{1},\dots,h_{M})\in\{0,1\}^{M}$, respectively. The interaction between the visible and hidden layers is given by the so-called weight matrix $w\in\mathbb{R}^{N}\times\mathbb{R}^{M}$ where the weight $w_{ij}$ is the connection strength between a visible unit $v_{i}$ and a hidden unit $h_{j}$. The biases $b_{i}\in\mathbb{R}.$ and $c_{j}\in\mathbb{R}$ are applied to visible unit $i$ and hidden unit $j$, respectively. Given random vectors $\mathbf{v}$ and $\mathbf{h}$, the energy function of the RBM is written as an Ising-type Hamiltonian

where the set of model parameters is denoted by $\theta\equiv\{w_{ij},b_{i},c_{j}\}$. The joint probability of finding $\mathbf{v}$ and $\mathbf{h}$ of the RBM is given by the Boltzmann distribution

where the partition function, $Z(\theta)\equiv\sum_{\mathbf{v},\mathbf{h}}e^{-E(\mathbf{v},\mathbf{h};\theta)}$, is the sum over all possible configurations. The marginal probabilities $p(\mathbf{v};\theta)$ and $p(\mathbf{h};\theta)$ 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 $\theta$ such that the marginal probability of the visible layer $p(\mathbf{v};\theta)$ becomes as close as possible to the unknown probability $p_{\rm data}(\mathbf{v})$ that generate the training data. Given identically and independently sampled training data ${\cal D}\in\{\mathbf{v}^{(1)},\dots,\mathbf{v}^{(D)}\}$, the optimal model parameters $\theta$ can be obtained by maximizing the likelihood function of the parameters, ${\cal L}(\theta|{\cal D})=\prod_{i=1}^{D}p(\mathbf{v}^{(i)};\theta)$, or equivalently by maximizing the log-likelihood function $\ln{\cal L}(\theta|{\cal D})=\sum_{i=1}^{D}\ln p(\mathbf{v}^{(i)};\theta)$. Maximizing the likelihood function is equivalent to minimizing the Kullback-Leibler divergence or the relative entropy of $p(\mathbf{v};\theta)$ from $q(\mathbf{v})$ Kullback and Leibler (1951); Cover and Thomas (2006)

where $q(\mathbf{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 $\mathbf{v}^{(i)}$ and a reconstructed visible vector $\bar{\mathbf{v}}^{(i)}$ of the RBM,

The stochastic gradient ascent method for the log-likelihood function is used to train the RBM. Estimating the log-likelihood 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. Hinton (2012); Fischer and Igel (2014); Melchior et al. (2016). Here we employ the CD-$k$ method.

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(\mathbf{v},\mathbf{h};\theta)$

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(\mathbf{v},\mathbf{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\times 2$ Bar-and-Stripe data are considered  Hinton and Sejnowski (1986); MacKay (2002). Fig. 2 shows the 6 possible $2\times 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_{ij}$, 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_{ij}$ 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(\mathbf{v})$ of the visible layer and $p(\mathbf{h})$ of the hidden layer before and after training. Note that the marginal probability $p(\mathbf{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 monitored 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_{\rm 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 $\sim 10,000$, the cost function $C$ continues approaching zero while the entropy $S$ and the Kullback-Leibler divergence $D_{\rm 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(\mathbf{v};\theta)$ to an unknown probability $q(\mathbf{v})$ which produces the data. Eq. (4) is expanded as

The second term $-\sum_{\mathbf{v}}q(\mathbf{v})\ln p(\mathbf{v};\theta)$ depends on the parameter $\theta$. As the training proceeds, $p(\mathbf{v};\theta)$ becomes close to $q(\mathbf{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(\mathbf{v})=p(\mathbf{v};\theta)$ that leads to $D_{\rm 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 $\mathbf{v}$ and the hidden random vector $\mathbf{h}$ are independent, i.e., $p(\mathbf{v},\mathbf{h};\theta)=p(\mathbf{v};\theta)p(\mathbf{h};\theta)$, 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}$, Reif (1965),

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 Araki and Lieb (1970); Nielsen and Chuang (2000). This property can be proved using the log inequality, $-\ln x\geq-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},\dots,a_{n}$ and $b_{1},\dots,b_{n}$,

In other words, Eq. (12) can be regarded as the negative of the relative entropy or Kullback-Leibler divergence of the joint probability $p(\mathbf{v},\mathbf{h})$ to the product probability $p(\mathbf{v})\cdot p(\mathbf{h})$,

For the 2$\times 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_{\rm KL}(q\,||\,p)$ as a function of epoch. Fig. 6 (a) shows that the Kullback-Leibler divergence, $D_{\rm 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.

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 $\theta$ slowly. If the parameters $\theta$ 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 after-training, $W_{\infty}=F_{1}-F_{0}\,.$ For switching $\theta$ 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 Jarzynski (1997, 2011) proved that for any switching rate the free energy difference $\Delta 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 $\theta$ through a sequence $\{\theta_{0},\theta_{1},\dots,\theta_{\tau}\}$, 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 $(\mathbf{v},\mathbf{h})$ of the RBM in configuration space. From the initial configuration, $(\mathbf{v}_{0},\mathbf{h}_{0})$ which is sampled from the initial Boltzmann distribution, Eq. (2), the trajectory $(\mathbf{v}_{0},\mathbf{h}_{0})\to(\mathbf{v}_{1},\mathbf{h}_{1})\to\cdots\to(\mathbf{v}_{\tau},\mathbf{h}_{\tau})$ is obtained using the Metropolis-Hastings algorithm of the Markov chain Monte-Carlo method Metropolis et al. (1953); Hastings (1970). Assuming the evolution is Markovian, the probability of taking a specific trajectory is the product of the transition probabilities at each step,

The transition $(\mathbf{v},\mathbf{h})\to(\mathbf{v}^{\prime},\mathbf{h}^{\prime})$ can be implemented by the Metropolis-Hastings algorithm based on the detailed balance condition for the fixed parameter $\theta$,

The work done on the RBM at epoch $i$ may be given by

The total work $W=\sum\delta W_{i}$ performed on the system is written as Crooks (1998)

Given the sequence of the model parameter $\{\theta_{0},\theta_{1},\dots,\theta_{\tau}\}$, the Markov evolution of the visible and hidden vectors $(\mathbf{v},\mathbf{h})\in\{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(\mathbf{v},\mathbf{h};\theta)$ of the RBM for the $2\times 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 $\mathbf{v}$ as a function of epoch. Before training, the visible vector $\mathbf{v}$ is distributed over all possible configurations, represented by the number $(0,\cdots,15)$. As the training progresses, the visible vector $\mathbf{v}$ becomes trapped into one of the six possible outcomes $(0,3,5,10,12,15)$.

In order to examine the relation between work done on the RBM during the training and the free energy difference, the Monte-Carlo simulation is performed to calculate the average of the work over paths generated by the Metropolis-Hastings algorithm of the Markov chain Monte-Carlo method. Each path starts from an initial state sampled from the uniform distribution over configuration space, as shown in Fig. 4 (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 $\langle{W}\rangle\approx-5.481$, and its standard deviation is $\sigma_{W}\approx 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. Jarzynski (2011) 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, $\langle{e^{-W}}\rangle_{\rm path}$ to check the Jarzynski equality, Eq. (16). The free energy difference can be estimated as

where $N_{\rm 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 $\Delta 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 Jarzynski (2006); Zuckerman and Woolf (2002); Lechner et al. (2006); Lechner and Dellago (2007); Halpern and Jarzynski (2016). As shown in Fig. 9, the distribution of work is given by the Gaussian distribution $\rho(W)$ with the mean $\langle W\rangle$ and the standard deviation $\sigma_{W}$. If the standard deviation $\sigma_{W}$ becomes larger, the peak position of $\rho(W)e^{-W}$ moves to the long tail of the Gaussian distribution. So the main contrition of the integration of $\langle e^{-W}\rangle$ comes from the rare realizations. Fig. 10 shows that the standard deviation $\sigma_{W}$ grows with epoch, so the error of the Monte-Carlo estimation of the exponential average of the work grows quickly.

If $\sigma_{W}^{2}\ll 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., $\sigma_{W}^{2}\gg k_{B}T$ $(=1)$, so the central limit theorem does not work and the above equation can not be applied Hendrix and Jarzynski (2001). Fig. 10 shows how the average of work, $\langle W\rangle_{\rm 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}=\Delta F$. The difference between the actual work and the reversible work is called the dissipative work, $W_{d}=W-W_{r}$ Crooks (1998). As depicted in Fig. 10, the magnitude of the dissipative work grows with training epoch.