Using Restricted Boltzmann Machines to Model Molecular Geometries

The RBM consists of two layers of neurons, a visible layer and a hidden layer (shown in Figure 2). Every visible node $\displaystyle v_{i}$ is connected to every hidden node $\displaystyle h_{j}$ by a set of weights, and each node has its own offset, or bias. The state or value of a given node is dependent on the state of the nodes it is connected to, as well as its bias. For ease of computation, the values of the weights are stored in a weight matrix $\displaystyle W$, where $\displaystyle W_{ij}$ represents the connection from $\displaystyle v_{i}$ to $\displaystyle h_{j}$. Meanwhile, the values of the visible and hidden biases are stored in bias vectors $\displaystyle a$ and $\displaystyle b$, respectively. The RBM is “restricted” in the sense that there are no connections between nodes in the same layer, which simplifies learning [16].

The visible layer serves as an input to the machine, where the number of visible nodes corresponds to the number of variables that make up the data. Values for the hidden nodes are then calculated by multiplying the visible nodes by the weights, adding the biases, and then applying some sort of activation function. This inference step can be viewed as a transformation from the space of observable parameters to the space represented by the hidden nodes. The method used for calculating between layers is formalized by a set of conditional probabilities that appear in the section below. In this paper we refer to a $\displaystyle n$-$\displaystyle k$-RBM as an RBM with $\displaystyle n$ visible nodes and $\displaystyle k$ hidden nodes.

A notable feature of the RBM, in contrast with other machine learning models, is that it does not have an “output” in the normal sense. Depending on the situation, the output of an RBM can be the visible layer, hidden filters, or energy. RBMs are trained in an unsupervised fashion (without labelling the data), so the hidden nodes identify their own labels during the course of training [5]. In a famous application of RBMs by Netflix to provide movie recommendations (see the Netflix prize [15]), the RBM was trained on a large set of movie ratings obtained from individual users, where the value of each visible node corresponded to the rating of a given movie. Based on the simple pattern that users who like movies from a certain director or genre are likely to enjoy other movies from that same category, the RBM came to associate hidden nodes with movie genres or directors. In this way, each hidden node elucidates a connection or correlation between visible nodes, which is true for continuous data as well. Note that understanding what each hidden node represents requires additional post hoc analysis.

As an energy based model, the RBM associates a scalar measure of energy to each state of the machine. Usually the energy equation is a combination of each layer times its respective biases and a term which relates the two layers. The goal of training is to minimize the overall energy of the RBM with respect to the training data. Once trained, the RBM associates lower energies with inputs that fall within the training distribution and higher energies with those that fall outside the training distribution. This becomes useful when generating new configurations because the RBM estimates the energy of a proposed configuration based on where it fits with the training distribution.

The simplest and most widespread version of the RBM is the binary-binary RBM (BBRBM), which has binary visible and binary hidden units. In the binary-binary RBM, an active unit is represented as 1 while inactive units are represented with 0. For a BBRBM, the energy is given by:

$$\displaystyle E(v,h)=-a^{T}v-b^{T}h-v^{T}Wh$$

where $\displaystyle v$ is a vector of visible states, $\displaystyle h$ is the hidden states, $\displaystyle a$ is the visible biases, $\displaystyle b$ is the hidden biases, and $\displaystyle W$ is the weight matrix representing the connections between visible and hidden nodes. The visible states $\displaystyle v$ serve as an input to the RBM whereas $\displaystyle a$, $\displaystyle b$, and $\displaystyle W$ are all parameters learned during training.

The joint probability of states $\displaystyle v$ and $\displaystyle h$ is taken from the Boltzmann distribution:

$$\displaystyle p(v,h)=\frac{1}{Z}e^{-E(v,h)}$$

where $\displaystyle Z$ is the partition function, defined as:

$$\displaystyle Z=\sum_{v}\sum_{h}e^{-E(v,h)}$$

which is a sum of over all possible states of the machine. The partition function $\displaystyle Z$ serves to normalize the joint probability distribution so that probabilities sum to 1. From this, we can deduce the marginal probability of a visible configuration $\displaystyle v$, given by the sum of probabilities over all possible hidden configurations:

	$\displaystyle\displaystyle p(v)$	$\displaystyle\displaystyle=\sum_{h}p(v,h)$		(1)
	$\displaystyle\displaystyle p(v)$	$\displaystyle\displaystyle=\frac{1}{Z}\sum_{h}e^{-E(v,h)}$		(2)

Because the state of the hidden layer depends on the state of the visible layer, we must use a conditional probability $\displaystyle p(h|v)$ to calculate the hidden states. This conditional probability is derived directly from the energy equation [20], giving us the activation of hidden states:

$$\displaystyle p(h=1|v)=\sigma(b+W^{T}v)$$

where $\displaystyle\sigma(x)$ represents the sigmoid activation function. Performing this operation only gives a set of probabilities that each hidden neuron is active; one must then sample a Bernoulli distribution with the given probabilities to reach the actual binary states of the hidden layer. Then, the visible states are computed again from the hidden states using:

$$\displaystyle p(v=1|h)=\sigma(a+Wh)$$

which gives the probabilities of the visible neurons adopting a value of one. Figure 3 shows how conditional probabilities are used to sample back and forth between layers.

A major drawback of the binary RBM is that it is only able to represent binary states or bit strings. As most empirical data is real-valued, the Gaussian RBM was developed to model continuous variables.

The Gaussian RBM (GRBM), also called the Gaussian-binary RBM, is an effective tool for modeling real-valued data. While the hidden layer remains binary, the visible layer can now adopt real values, with an additional parameter $\displaystyle\sigma$ representing the standard deviation for each visible node. The equation for the joint energy becomes:

$\displaystyle\displaystyle E_{GRBM}(v,h)=\frac{\left\lVert v-a\right\rVert^{2}}{2\sigma^{2}}-b^{T}h-\frac{v^{T}}{\sigma^{2}}Wh$

(3)

where $\displaystyle\left\lVert\cdot\right\rVert$ is the Euclidean norm. This equation is similar to the binary energy except that the visible states are divided by $\displaystyle\sigma^{2}$ as a form of normalization, and the first term is replaced by $\displaystyle\frac{\left\lVert v-a\right\rVert^{2}}{2\sigma^{2}}$ which serves as a parabolic containment of the visible states. This means that the overall energy increases the further the visible states $\displaystyle v$ are from the visible biases $\displaystyle a$. Whereas the states in a binary RBM are bounded by 0 and 1, it is important that the visible states of a Gaussian RBM are restrained by this parabolic term in order to prevent them from trailing off.

Similar to the binary RBM, the conditional probability of the hidden states given the visible states is:

$$\displaystyle p(h=1|v)=\sigma(b+W^{T}\frac{v}{\sigma^{2}})$$

whereas the conditional probability of the visible states given the hidden states is:

$$\displaystyle p(v\ |\ h)=\mathcal{N}(a+Wh,\sigma^{2})$$

where $\displaystyle\mathcal{N}(\mu,\sigma^{2})$ is a Gaussian function with mean $\displaystyle\mu$ and variance $\displaystyle\sigma^{2}$. In practice this amounts to adding Gaussian noise after calculating visible states. A detailed derivation of the conditional distributions can be found in [20].

While initially there appears to be nothing inherently Gaussian about this energy equation, the probability starts to take on the form of a Gaussian function when substituted into the Boltzmann distribution:

$\displaystyle\displaystyle p(v,h)$

$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{\left\lVert v-a\right\rVert^{2}}{2\sigma^{2}}+b^{T}h+\frac{v^{T}}{\sigma}Wh}$

If we ignore the second and third terms of the energy equation, we see that the distribution of visible states follows a Gaussian distribution centered at the visible bias $\displaystyle a$ with variance $\displaystyle\sigma^{2}$. Moreover, it has previously been shown using these equations that the RBM is equivalent to a mixture of Gaussian (MoG) model, with the locations of each Gaussian represented by the column vectors of the weight matrix (for a detailed derivation see [11]).

The energy and sampling equations are usually simplified by normalizing the data to zero mean and unit variance, and then setting $\displaystyle\sigma=1$. However, one can choose to calculate the $\displaystyle\sigma$ of the data in advance, or alternatively one can use $\displaystyle\sigma$ as a separate parameter which is optimized during training. While training data need not follow a Gaussian distribution, the GRBM works best when modeling this type of distribution. Fortunately, most distributions found in nature are Gaussian due to the central limit theorem.

The goal of training the RBM is to maximize the likelihood of the data under the model. By tweaking the weights and biases, the RBM can associate higher probabilities with configurations found in the training data. For practical purposes, we work with the logarithm of the likelihood which allows us to write products as the sum of logarithms. Since $\displaystyle\log(x)$ is a monotonically increasing function, maximizing the log-likelihood is the same as maximizing the likelihood. The log-likelihood for a training sample $\displaystyle x$ is given by:

	$\displaystyle\displaystyle\mathcal{L}(x)$	$\displaystyle\displaystyle=\log p(x)$
		$\displaystyle\displaystyle=\log\sum_{h}p(x,h)$
		$\displaystyle\displaystyle=\log\sum_{h}\frac{e^{-E(x,h)}}{Z}$
		$\displaystyle\displaystyle=\log\sum_{h}e^{-E(x,h)}-\log Z$

using (1) and (2). By optimizing the log likelihood of the training data, we maximize the likelihood of the training samples over all possible samples in our visible space.

In the case where there is an entire set of training samples, we take the average log-likelihood $\displaystyle\hat{\ell}$ by computing the expectation of the log-likelihood over all the samples:

$$\displaystyle\hat{\ell}=\bigg{\langle}\log\sum_{h}e^{-E(x,h)}\bigg{\rangle}-\log Z$$

Average log-likelihood is a rigorous way of monitoring the training and convergence of an RBM and demonstrating its modeling capacity in statistical terms. The most difficult part about estimating this likelihood is calculating the partition function $\displaystyle Z$, which is intractable in most cases. In this study we use importance sampling to estimate $\displaystyle Z$ whenever calculating log-likelihood (see [8] for more details).

The gradient of $\displaystyle\hat{\ell}$, which is given by the derivative of $\displaystyle\hat{\ell}$ with respect to the model parameters $\displaystyle\theta$, comes out to be the difference between the data-based distribution and the distribution given by the entire model:

$$\displaystyle\frac{d\hat{\ell}}{d\theta}\propto\bigg{\langle}\frac{dE(x,h)}{d\theta}\bigg{\rangle}_{x}-\bigg{\langle}\frac{dE(v,h)}{d\theta}\bigg{\rangle}_{v}$$

where $\displaystyle x$ are explicit training samples and $\displaystyle v$ are samples from the model distribution.

Computing the partial derivatives of the energy function with respect to each parameter provides us with the gradient approximations for the weights and biases:

$$\displaystyle\frac{\partial\hat{\ell}}{\partial a}\propto\langle{x}-a\rangle_{{x}}-\langle v-a\rangle_{v}$$

$$\displaystyle\frac{\partial\hat{\ell}}{\partial b}\propto\langle p(h=1|{x})\rangle_{{x}}-\langle p(h=1|v)\rangle_{v}$$

$$\displaystyle\frac{\partial\hat{\ell}}{\partial W}\propto\langle{x}\ p(h=1|{x})^{T}\rangle_{{x}}-\langle v\ p(h=1|{v})^{T}\rangle_{v}$$

$$\frac{\partial\hat{\ell}}{\partial\sigma}\propto\bigg{\langle}\frac{\left\lVert x-a\right\rVert-2x^{T}Wp(h=1|x)}{\sigma^{3}}\bigg{\rangle}_{{x}}\\ -\bigg{\langle}\frac{\left\lVert v-a\right\rVert-2v^{T}Wp(h=1|v)}{\sigma^{3}}\bigg{\rangle}_{v}$$

These gradients are used as the update rules for the RBM. Typically training is done in batches, and the gradients are computed after each batch. The gradients are then added to the existing RBM parameters which gives rise to a new set of weights and biases. While the first term for each gradient can be calculated directly from the training data, the second term is almost always intractable, as it requires independent samples from an unknown model distribution.

Several algorithms are currently available for estimating the second term. The most widely used algorithm is Contrastive Divergence. In Contrastive Divergence, sampling between visible and hidden layers is used to create a Markov chain (as shown in Figure 3). The Markov chain is initialized at a training point, and the conditional probabilities are used to get the visible and hidden states after $\displaystyle k$ iterations, in a process known as Gibbs sampling. The value of these layers after sampling serves as a sufficient estimation of the model’s expectation. Typically $\displaystyle k$ is set to 1, as it provides a good estimation of the gradient and also minimizes computation time.

Persistent Contrastive Divergence (PCD) is another algorithm which has been proposed for estimating the model distribution. Instead of restarting the Markov chain for each data point, one persistent Markov chain is retained in memory and extended after each batch. While its success has mostly been reported in binary RBMs [18], some have had success in using PCD to train Gaussian RBMs [3]. A reported improvement on PCD is the Parallel Tempering algorithm (PT), which uses multiple Markov chains at different temperatures with a certain probability that states from different chains will swap [2]. The reasoning behind PT is that having different temperatures will ensure that the Markov chain is fully exploring high energy states. Though these algorithms have all been successfully implemented in binary RBMs, there is little evidence that they work in Gaussian RBMs. In practice, we found that both PCD and PT led to divergence unless a very small learning rate was used ($\displaystyle\alpha=0.0001$), which led to long training times with marginal improvements in likelihood. Overall, CD-1 was simplest to use and almost never led to divergence. Our results match those in [11] which conclude that CD is the best algorithm currently suited for training GRBMs.

Because the RBM provides us with a measure of energy, we can sample new configurations from a trained RBM using the Metropolis Monte Carlo algorithm. This is done by initializing a random configuration and evaluating its energy. Then for each iteration, a new sample $\displaystyle x^{\prime}$ is generated by adding a random displacement vector to the previous sample $\displaystyle x$ and evaluating the energy of the new sample. If the energy is lower than the previous, the new sample is accepted as part of the simulation. If not, the acceptance ratio is calculated using the Boltzmann distribution:

$$\displaystyle\frac{p(x^{\prime})}{p(x)}=\exp{\left(\frac{E(v,h)-E(v^{\prime},h^{\prime})}{k_{B}T}\right)}$$

Finally we generate a random number between 0 and 1, and if our number falls below this acceptance ratio, we keep the configuration as part of our ensemble. The Boltzmann constant $\displaystyle k_{B}$ and temperature parameter $\displaystyle T$ are usually set to one, though they may be useful in simulating higher temperatures or exploring higher energy states.

The benefit of this Metropolis Monte Carlo method is that the normalization for the RBM need not be known. Whereas the true probability density of the joint states requires knowledge of the partition function $\displaystyle Z$, which is intractable in most cases, the term cancels when calculating the acceptance ratio at each step in the simulation. The result is a Monte Carlo simulation which follows the probability density given by the model.

As we can see, the energy calculations for an RBM are much simpler than the energy calculations for a quantum system. A typical Hartree-Fock calculation has computational complexity $\displaystyle\mathcal{O}(n^{4})$, where $\displaystyle n$ is the total number of basis functions, since the number of two-electron integrals necessary to build the Fock matrix is $\displaystyle n^{4}$ [4]. In practice this often becomes $\displaystyle\mathcal{O}(n^{3})$ as most programs will ignore integrals that are close to zero. Depending on the choice of basis set and the size of the atoms, one can have anywhere from 3 to 100 basis functions for one given atom. Therefore the complexity becomes $\displaystyle\mathcal{O}(b^{3}a^{3})$ where $\displaystyle b$ is the average number of basis functions per atom and $\displaystyle a$ is the number of atoms. Because the complexity scales at least cubically with system size, these calculations become quite expensive as one increases either the atoms or basis functions included.

Meanwhile, RBM energy computations have complexity $\displaystyle\mathcal{O}(NH)$ where $\displaystyle N$ is the number of visible nodes and $\displaystyle H$ is the number of hidden nodes. In the proposed method, only three nodes are needed per additional atom, as compared to the large number of basis functions needed in Hartree-Fock calculations. Thus the complexity is $\displaystyle\mathcal{O}(3aH)$, expressed in terms of the number of atoms $\displaystyle a$. If we take the number of hidden nodes to be the same as the number of visible nodes, then the computational complexity grows quadratically with system size, and increasing atom size has no impact on the complexity of calculations. Moreover these calculations are performed using matrix multiplication which is simpler than integration.

To assess the validity of the methodology and accuracy of training over a wide range of distributions, we trained the RBM on a set of toy datasets. In the first toy dataset, random samples were taken from a Laplacian distribution and fitted to the line $\displaystyle y=x$, with Gaussian noise of unit variance added in both dimensions. The second dataset was generated by sampling from two independent Laplacian distributions and using a mixing matrix to merge the sources. An observable $\displaystyle x$ is generated from independent Laplacian sources $\displaystyle s=(s_{1},s_{2})$ using a mixing matrix $\displaystyle A$:

$$\displaystyle x=As$$

where $\displaystyle A$ is a linear transformation from the source space to the observable space. Given a set of observables, the RBM is capable of estimating the mixing matrix and recovering source data [11]. In the context of this task, independent component analysis is the standard computational method for separating observed data into source components. For this reason, the independent components of the original data were compared with the independent components of the data generated by RBM. Independent component analysis was performed using the FastICA algorithm.

Moving to more complex data, different RBMs were also evaluated on natural image patches taken from [11]. These image patches come from randomly sampling the greyscale images found in the van Hateran natural image database [19]. Each patch consists of 14 x 14 pixels, and for our purposes we used only 10,000 of these patches during training.

Meanwhile, molecular training data for H₂O and ethane was generated by performing molecular Monte Carlo simulations and evaluating the energy of each proposed sample using the restricted Hartree-Fock method. The Hartree-Fock calculations made use of the basis set cc-pVDZ because it had the best agreement with literature values for the bond angle and bond length of H₂O. These calculations were run using the PySCF module for quantum chemistry.

There are many types of RBMs apart from those already discussed. By altering the equation for the overall energy, a variety of hidden and visible units can be implemented. Table 1 summarizes the energy equations and conditional probabilities for a few different types of RBMs. These RBMs are mostly defined by the activation functions used to calculate the hidden layer, including sigmoid, Gaussian, ReLU.

Noticing the absence of an RBM with Tanh activation function, we derived a new energy function and conditional probabilities for an RBM with Tanh hidden units using the method outlined in [14] (see appendix for derivation). After developing the Tanh RBM, we implemented the full range of RBMs defined in Table 1, which use common activation functions from machine learning. This implementation can be found at https://github.com/peter1255/RBM_chem.

Depending on the correlations found within the data, different activation functions may be better for modeling different data. It has been shown that non-linear hidden units expand the capabilities of the RBM by allowing the model to represent non-linear correlations between visible units [14]. The type of non-linearity chosen may have an effect on the model’s representation of the data, though a more detailed comparison of RBMs is needed.

The RBMs found in Table 1 were implemented and trained on both natural image data and H₂O geometries. During 20 epochs of training, the performances of RBMs were compared through thorough calculations of the log-likelihood. These results are shown in Figure 4.

First proposed by [12], the ReLU RBM trains well but lacks an appropriate joint energy equation [17], making it impossible to calculate log-likelihood or generate simulations. Though we tried using the Gaussian joint energy equation (3) in test simulations not shown here, values quickly diverged as large hidden terms outweighed the parabolic containment term. Making a slight modification of ReLU into Leaky ReLU gives a viable energy equation which can be successfully trained and used for both sampling and calculating log-likelihood [10]. While this version of the RBM has not yet reached widespread use, we show in 4 that it performs on a level comparable with the previously established GRBM.

Out of all the RBMs tested, the Tanh RBM achieved the best maximum log-likelihood on both natural image data and H₂O molecular geometries. However, when trained on data sets with few parameters, like H₂O data, the Tanh RBM has trouble maintaining convergence and tends to diverge after reaching its maximum. While the Tanh RBM shows great promise for representing a variety of data, additional research must be done to learn how to prevent divergence in general cases.

The model that achieves second-best performance in terms of log-likelihood is the Gaussian RBM, which is quite close in the case of molecular data (Figure 4). For this reason, we selected this version of the model to test the proposed simulation method on molecular data.

Figure 5 compares toy datasets with the datasets generated by a trained GRBM using the Monte Carlo method outlined above. The left column displays the original datasets, while the middle column displays the data generated by RBM. The two datasets appear quite similar to one another, though some of the pointedness in the original data is lost through the RBM’s Gaussian probability distribution. After performing independent component analysis on both training and simulated data, we can see that the independent components (shown using arrows) from the simulated data (middle column) match those from the original data (left column).

Then, evaluating the energy of each point in the visible space using equation (3), we graphed the energy contour represented by the GRBM (right column). By moving along this energy contour we are able to generate the states represented in the middle column. Because the Monte Carlo method favors lower energy states, most of the samples lie within the middle of the contour. Furthermore, the arrows represented by each hidden node of the RBM match the independent components of the reconstructed data, showing that the Monte Carlo sampling method preserves the original components of the data.

Furthermore, very few hidden nodes were required for reproducing the data. In the case of the linear distribution which has only one independent component, only two hidden nodes were needed. In the case of the cross shaped distribution, only four hidden units were needed. In general, two hidden nodes are needed for each independent component in the data: because the hidden nodes are binary, they can only represent their positive vector span. Therefore an additional hidden node is needed to represent the opposite direction. This also provides a useful heuristic for determining the number of hidden nodes to include in an RBM. In this case, using any more arrows would be redundant, as the arrows begin to overlap.

After learning toy distributions, the RBM was trained on molecular geometries of H₂O generated by Hartree-Fock. Since the geometry of H₂O is defined by three parameters (two bonds and one angle), three visible nodes were included in the GRBM. The GRBM was trained on a set of 10000 geometries. Figure 3 displays the overall bond and angle energies evaluated using Hartree-Fock (top row) and GRBM (bottom row).

Not only is the energy contour of the RBM similar to the energy determined by quantum methods, the RBM was able to pick up on subtleties including the rotated parabola which defines the bond and angle energies (first and second columns). Looking at the bond and angle energies given by the RBM, the left arm of the parabola shows the energy associated with atom-atom repulsion, while the right side of the graph shows the energy associated with atom-atom attraction. Both methods match the literature values for the bond length and angle of H₂O.

Similar to the previous figure, energy contours (third and fourth columns) were generating by fixing one parameter to its given minimum and evaluating the energy with respect to the other two parameters. As we can see, the energy contours inferred by GRBM are quite similar to the original contours given by Hartree-Fock. Though the stratification of the contour layers is slightly different (a fundamental limitation of the RBM model structure) the overall shape of the contours is quite similar. The GRBM accurately represents the relationship between different data parameters.

Finally, we use the RBM to model a larger molecule: ethane. Though it only consists of 6 atoms, ethane requires 18 internal coordinates to represent all its various bonds, angles, and dihedrals. In general, a molecule with $\displaystyle n$ atoms requires $\displaystyle 3n-6$ unique coordinates, which can be stored in the form of a z-matrix. Therefore we train an RBM with 18 visible nodes on the internal ethane coordinates, and then sample from it to generate new geometries.

Figure 7 shows the difference between the original dataset (black) and the GRBM generated data (red). The reconstructed data is able to approximate molecular geometries by fitting Gaussian curves to the original data. For simple unimodal data, the GRBM does a good job of reproducing the existing shape of the data. However, for more complex distributions (i.e. the bimodal dihedral angle found in coordinate 15 of 7) the GRBM is unable to provide a fit for this complex pattern and instead regresses to the mean in between the two modes. Nevertheless, the rest of the marginal distributions are well represented, and the GRBM captures the essence of the joint distributions. Through our Monte Carlo sampling algorithm, the GRBM provides an accurate picture of the dynamics of an ethane molecule.

Using the RBM it is easy to evaluate the energies of any configurations at a fast speed. By scaling the energy units of the RBM to proper Hartree units, we could use the RBM to represent bond and angle energies without having to perform any quantum calculations. Moreover, because of the low computational cost in generating states of a water molecule, the technique here could also be extended to model an entire body of liquid water in an aqueous solution.

The software package used for this study along with a tutorial for how to use it to model molecular data is available on Github at https://github.com/peter1255/RBM_chem.

First we define the activation function $\displaystyle f(\eta_{j})=\tanh(\eta_{j})$ where $\displaystyle\eta_{j}=b_{j}+\sum_{i=1}^{I}W_{ij}v_{j}$ assuming data is normalized to unit variance. The inverse function is $\displaystyle f^{-1}(h_{j})=\tanh^{-1}(h_{j})$. The corresponding antiderivatives are:

	$\displaystyle\displaystyle F(\eta_{j})$	$\displaystyle\displaystyle=\int f(\eta_{j})d\eta$
		$\displaystyle\displaystyle=\ln{\cosh{\eta_{j}}}+C$
	$\displaystyle\displaystyle F^{*}(h_{j})$	$\displaystyle\displaystyle=\int f^{-1}(h_{j})dh$
		$\displaystyle\displaystyle=h_{j}\tanh^{-1}(h_{j})+\frac{1}{2}\ln(1-h_{j}^{2})+C$

Meanwhile, the activation for the visible units is linear so that $\displaystyle\overline{f}(\nu_{i})=\nu_{i}$ where $\displaystyle\nu_{i}=a_{i}+\sum^{J}_{j=1}W_{ij}h_{j}$. Following the same step we get that $\displaystyle\overline{F}(\nu_{i})=\frac{1}{2}\nu_{i}^{2}$ and $\displaystyle\overline{F}^{*}(v_{i})=\frac{1}{2}v_{i}^{2}$.

From Ravanbakhsh et al. [14], the conditional distribution is defined as:

$$\displaystyle p(h_{j}|\eta_{j})=\exp(h_{j}\eta_{j}-F(\eta_{j})-F^{*}(h_{j})+g(h_{j}))$$

Replacing the values of $\displaystyle F(\eta_{j})$ and $\displaystyle F^{*}(h_{j})$ we get the conditional probability:

	$\displaystyle\displaystyle p(h_{j}|\eta_{j})$	$\displaystyle\displaystyle=\exp(h_{j}\eta_{j}-\ln{\cosh{\eta_{j}}}$
	$\displaystyle\displaystyle-$	$\displaystyle\displaystyle h_{j}\tanh^{-1}(h_{j})-\frac{1}{2}\ln(1-h_{j}^{2})+g(h_{j}))$
		$\displaystyle\displaystyle=\exp(\tanh(\eta_{j})\eta_{j}-\ln{\cosh{\eta_{j}}}$
	$\displaystyle\displaystyle-$	$\displaystyle\displaystyle\tanh(\eta_{j})\tanh^{-1}(\tanh(\eta_{j})+g(h_{j}))-\frac{1}{2}\ln(1-h_{j}^{2})+g(h_{j}))$
		$\displaystyle\displaystyle=\exp(\tanh(\eta_{j})\eta_{j}-\ln{\cosh{\eta_{j}}}$
	$\displaystyle\displaystyle-$	$\displaystyle\displaystyle\tanh(\eta_{j})\tanh^{-1}(\tanh(\eta_{j}))-\frac{1}{2}\ln(1-h_{j}^{2})+g(h_{j}))$
		$\displaystyle\displaystyle=\exp(-\ln{\cosh{\eta_{j}}}-\frac{1}{2}\ln(1-h_{j}^{2})+g(h_{j}))$
		$\displaystyle\displaystyle=\frac{1}{\cosh{\eta_{j}}}\frac{1}{\sqrt{1-\tanh^{2}{\eta_{j}}}}\exp(g(h_{j}))$
		$\displaystyle\displaystyle=\frac{1}{\cosh{\eta_{j}}\operatorname{sech}{\eta_{j}}}\exp(g(h_{j}))$
		$\displaystyle\displaystyle=\exp(g(h_{j}))$

For simplicity we define $\displaystyle g(h_{j})=-\frac{1}{2}h_{j}^{2}-\log\sqrt{2\pi}$ so that the conditional probability distribution becomes a Gaussian with mean $\displaystyle h$ and unit variance, meaning that Gaussian noise is added to the hidden units after applying the Tanh activation. Similarly we get that $\displaystyle p(v|\nu)=\mathcal{N}(\nu,1)$, as expected for the Gaussian visible layer.

The joint energy equation was generalized in Yang et al.[22] through the following:

$$\displaystyle E(v,h)=-v^{T}Wh+\overline{F}^{*}(v)+F^{*}(h)$$

Replacing $\displaystyle\overline{F}^{*}$ and $\displaystyle F^{*}$ we get the energy for a joint configuration of the tanh RBM:

$$\displaystyle E(v,h)=-v^{T}Wh+\frac{1}{2}||v||^{2}+\sum^{J}_{j=1}(h_{j}\tanh^{-1}(h_{j})+\frac{1}{2}\ln(1-h_{j}^{2}))$$

From the joint energy we can derive the marginal probability of a given visible configuration $\displaystyle v$:

	$\displaystyle\displaystyle p(v)$	$\displaystyle\displaystyle=\frac{1}{Z}\int_{h}e^{-E(v,h)}dh$
		$\displaystyle\displaystyle=\frac{1}{Z}\int_{h}e^{v^{T}Wh-\frac{1}{2}||v||^{2}-\sum^{J}_{j=1}(h_{j}\tanh^{-1}(h_{j})-\frac{1}{2}\ln(1-h_{j}^{2}))}dh$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\int_{h}e^{v^{T}Wh-\sum^{J}_{j=1}(h_{j}\tanh^{-1}(h_{j})+\frac{1}{2}\ln(1-h_{j}^{2}))}dh$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\int_{h}e^{\sum^{J}_{j=1}\eta_{j}h_{j}-h_{j}\tanh^{-1}(h_{j})-\frac{1}{2}\ln(1-h_{j}^{2})}dh$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\int_{h}\prod^{J}_{j}e^{\eta_{j}h_{j}-h_{j}\tanh^{-1}(h_{j})-\frac{1}{2}\ln(1-h_{j}^{2})}dh$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\prod^{J}_{j}\int_{h_{j}}e^{\eta_{j}h_{j}-h_{j}\tanh^{-1}(h_{j})-\frac{1}{2}\ln(1-h_{j}^{2})}dh_{j}$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\prod^{J}_{j}e^{\eta_{j}\tanh(\eta_{j})-\tanh(\eta_{j})\tanh^{-1}(\tanh(\eta_{j}))-\frac{1}{2}\ln(1-\tanh^{2}(\eta_{j}))}$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\prod^{J}_{j}e^{-\frac{1}{2}\ln(1-\tanh^{2}(\eta_{j}))}$
		$\displaystyle\displaystyle=\frac{1}{Z}e^{-\frac{1}{2}||v||^{2}}\prod^{J}_{j}e^{-\frac{1}{2}\ln(1-\tanh^{2}(b_{j}+\sum_{i=1}^{I}W_{ij}v_{j}))}$

where Z is the partition function.

Since the marginal energy is given by

$$\displaystyle E(v)=-\log p(v)=F(v)-\log(Z)$$

where $\displaystyle F(v)$ is the free energy of the given configuration, we have the following expression for the free energy

$$\displaystyle F(v)=\frac{||v||^{2}}{2}+\sum^{J}_{j}\frac{1}{2}\log(1-\tanh^{2}(b_{j}+\sum_{i=1}^{I}W_{ij}v_{j}))$$