Can Boltzmann Machines Discover Cluster Updates ?

TL;DR

This paper demonstrates that Boltzmann machines can not only learn physical distributions but also discover new cluster Monte Carlo algorithms by leveraging their latent representations, inspiring innovative computational methods.

Contribution

It reveals that Boltzmann machines can autonomously discover cluster algorithms for Monte Carlo simulations, a novel application beyond traditional generative modeling.

Findings

Boltzmann machines can identify clusters in physical systems.

They can discover new Monte Carlo algorithms.

Effective in models with complex interactions.

Abstract

Boltzmann machines are physics informed generative models with wide applications in machine learning. They can learn the probability distribution from an input dataset and generate new samples accordingly. Applying them back to physics, the Boltzmann machines are ideal recommender systems to accelerate Monte Carlo simulation of physical systems due to their flexibility and effectiveness. More intriguingly, we show that the generative sampling of the Boltzmann Machines can even discover unknown cluster Monte Carlo algorithms. The creative power comes from the latent representation of the Boltzmann machines, which learn to mediate complex interactions and identify clusters of the physical system. We demonstrate these findings with concrete examples of the classical Ising model with and without four spin plaquette interactions. Our results endorse a fresh research paradigm where intelligent machines are designed to create or inspire human discovery of innovative algorithms.

Equations

$$A(\mathbf{s}\to {\mathbf{s}}^{\prime })=\min \left[1,\frac{p(\mathbf{s})}{p({\mathbf{s}}^{\prime })}\cdot \frac{\pi ({\mathbf{s}}^{\prime })}{\pi (\mathbf{s})}\right].$$

$$\pi (\mathbf{s})=\exp \left(\beta J\sum _{\ell }\prod _{i\in \ell }{s}_i\right),$$

$$E(\mathbf{s},\mathbf{h})=-\sum _{\ell }\left(W\prod _{i\in \ell }{s}_i+b\right)h_{\ell }.$$

$$p(h_{\ell }=1|\mathbf{s})=\sigma \left(W\prod _{i\in \ell }{s}_i+b\right),$$

$$\frac{1+e^{b+W}}{1+e^{b-W}}=e^{2\beta J}$$

$$\pi (\mathbf{s})=\exp \left(\beta J\sum _{\ell }\prod _{i\in \ell }{s}_i+\beta K\sum _{\wp }\prod _{i\in \wp }{s}_i\right),$$

$$\exp \left(\beta K\prod _{i\in \wp }{s}_i\right)=\frac{e^{-\beta K}}{2}\sum _{h_{\wp }\in \{0,1\}}\exp \left[W\left(h_{\wp }-\frac{1}{2}\right){\mathcal{F}}_{\wp }(\mathbf{s})\right],$$

$$E(\mathbf{s},\mathbf{h})=-\sum _{\ell }\left[\beta J+W\sum _{\wp }\left(h_{\wp }-\frac{1}{2}\right)\left({\delta }_{\ell {\ell }_{\wp }}+{\delta }_{\ell {\bar{\ell }}_{\wp }}\right)\right]\prod _{i\in \ell }{s}_i$$

$$p(h_{\wp }=1|\mathbf{s})=\sigma \left(W{\mathcal{F}}_{\wp }(\mathbf{s})\right).$$

$$E(\mathbf{s},\mathbf{h})=-E(\mathbf{s})-\sum _{\alpha }\left[W_{\alpha }{\mathcal{F}}_{\alpha }(\mathbf{s})+b_{\alpha }\right]h_{\alpha },$$

$$A(\mathbf{s}\to {\mathbf{s}}^{\prime })=\min \left[1,\frac{T({\mathbf{s}}^{\prime }\to \mathbf{s})}{T(\mathbf{s}\to {\mathbf{s}}^{\prime })}\cdot \frac{\pi ({\mathbf{s}}^{\prime })}{\pi (\mathbf{s})}\right].$$

$$F_{\mathbf{h}}(\mathbf{s}\to {\mathbf{s}}^{\prime })p(\mathbf{s},\mathbf{h})=F_{\mathbf{h}}({\mathbf{s}}^{\prime }\to \mathbf{s})p({\mathbf{s}}^{\prime },\mathbf{h}).$$

$$\frac{T(\mathbf{s}\to {\mathbf{s}}^{\prime })}{T({\mathbf{s}}^{\prime }\to \mathbf{s})}$$

Full text

Introduction–

It is intriguing to wonder whether artificial intelligence can make scientific discoveries Langley et al. (1987); Alai (2004); Nilsson (2009); Gil et al. (2014) ever since the dawn of AI. With recent rapid progress in machine learning, AI is reaching human-level performance in many tasks and it is becoming even more optimistic that AI can indeed achieve the lofty goal of making scientific discoveries. Focussing on the computational study of many-body physical problems with relatively well-defined rules and targets, e.g. distinguishing phases of matter or finding the lowest energy state, the above question is addressed recently in Refs. Carrasquilla and Melko (2017); van Nieuwenburg et al. (2017); Wang (2016); Carleo and Troyer (2017).

An equally interesting question is whether AI can invent, or at least, inspire human discovery of new problem-solving strategies, i.e. new algorithms. In this respect, the two examples from DeepMind, where computers discover optimal strategy in the video game Breakout Mnih et al. (2015) and master the board game Go Silver et al. (2016), are particularly inspiring. These are vivid examples of an intelligence agent which discovers nontrivial problem-solving strategies even unknown to its programmer.

Along this line, it is highly desirable to devise new efficient Monte Carlo algorithms to sample configuration spaces of physical problems more rapidly. The existing ones such as hybrid Monte Carlo Duane et al. (1987), cluster algorithms Swendsen and Wang (1987); Wolff (1989), loop algorithm Evertz et al. (1993), and worm algorithm Prokof’ev and Svistunov (2001) are all landmark achievements in computational physics and find wide applications in physical, statistical and biological problems. There are some recent efforts to improve the Monte Carlo sampling Huang and Wang (2017); Liu et al. (2017, 2016); Huang et al. (2016); Xu et al. (2016) using ideas and techniques from machine learning. Similar approaches were also discussed in statistics literature Neal (1996); Liu (2008); Rasmussen et al. (2003) where one uses surrogate functions to guide and accelerate hybrid Monte Carlo calculations Duane et al. (1987). As suggested in Huang and Wang (2017), using intelligent Boltzmann Machines (BM) Ackley et al. (1985); Hinton and Sejnowski (1986) opens possibilities of algorithmic innovations because they can discovery unknown algorithmic strategies instead of merely acting as cheaper surrogate functions.

BM is an energy based model consists of stochastic visible ($\mathbf{s}$) and hidden ($\mathbf{h}$) variables illustrated in Fig 1. The BM architecture is specified by an energy function $E(\mathbf{s},\mathbf{h})$ which depends on the connectivity, connection weights, and biases of the units. The joint probability distribution of the units follows the Boltzmann distribution $p(\mathbf{s},\mathbf{h})=e^{-E(\mathbf{s},\mathbf{h})}$. One can train the BM by tuning its energy function such that the marginal distribution of the visible variables $p(\mathbf{s})=\sum_{\mathbf{h}}p(\mathbf{s},\mathbf{h})$ approximates the target probability distribution $\pi(\mathbf{s})$ of a dataset. The hidden units of a BM mediate interactions between the visible units and serve as internal representations of the data.

A successfully trained BM can capture salient features of the input data. For example, the BM learns about pen strokes from an image dataset of handwritten digits 111http://deeplearning.net/tutorial/rbm.html. Once trained, the BM can generate new samples from the learned distribution. To keep the physical simulation unbiased, the BM recommended update of the visible units $\mathbf{s}\rightarrow\mathbf{s}^{\prime}$ is accepted with the probability according to Metropolis-Hastings rule Metropolis et al. (1953); Hastings (1970); SM ,

[TABLE]

Equation (1) shows that the BM guides the simulation and increases the acceptance rate by exploiting the knowledge learned from data. In particular, one can even achieve a rejection free Monte Carlo simulation scheme if the BM perfectly describes the target probability distribution. This is in principle possible because BM is a universal approximator of discrete probability distributions Freund and Haussler (1994); Le Roux and Bengio (2008); Montufar and Ay (2011). The expressive powers of BM were studied recently from physics perspectives Chen et al. (2017); Deng et al. (2017); Gao and Duan (2017); Huang and Moore (2017). See Refs. Carleo and Troyer (2017); Torlai and Melko (2016a, b) for other recent applications of BM to quantum and statistical physics problems.

Reference Huang and Wang (2017) argues that the efficient simulation of the BM, and more importantly, its ability to capture high level features make BM ideal recommender systems to accelerate Monte Carlo simulation of challenging physical problems. In particular, Ref. Huang and Wang (2017) employs a restricted architecture of BM where the connections are limited to be between the visible and hidden units. Such a restricted BM can be sampled efficiently by blocked Gibbs sampling alternating between the hidden and visible units. The generative sampling of the restricted BM already gives rise to nonlocal updates because the BM learns about collective density correlations from the Monte Carlo data. This suggests that besides being used as a general purpose recommender engine for accelerating Monte Carlo simulations, the BM may discover conceptually new efficient updates.

In this paper, we demonstrate the BM’s creative power by exact constructions and then present a general framework exploit the power. The crucial insight is that the hidden units of the BM can learn to mediate complex interactions between the visible units and decouple the visible units into disconnected clusters. The generative sampling of the BM then automatically propose efficient cluster updates. To encourage these discoveries, it is crucial to design the BM in a suitable architecture and allow its parameters adapt to the physical distribution via learning.

Example: Ising model–

To make the discussions concrete, we start with the classical Ising model and show that the generative sampling of the BM encompasses a wide range of celebrated cluster algorithms Swendsen and Wang (1987); Niedermayer (1988); Kandel and Domany (1991); Kawashima and Gubernatis (1995); Higdon (1998). The Boltzmann weight of the Ising model reads

[TABLE]

where $\beta=1/T$ is the inverse temperature and $J$ is the coupling constant. We consider ferromagnetic coupling $J>0$ in the following for clarity. The considerations are nevertheless general and valid for the antiferromagnetic case as well. Equation (2) consists of a summation over links $\ell$ of a lattice and a product over Ising spins $s_{i}\in\{-1,1\}$ reside on the vertices connected by the link.

To devise a BM inspired cluster update of the Ising model, we consider the architecture illustrated in Fig. 2(a). We view the Ising spins as visible variables and introduce binary hidden variables $h_{\ell}\in\{0,1\}$ on the links of the lattice. These units are coupled according to the following energy function

[TABLE]

Equation (3) is a high-order BM Sejnowski (1986) because the interaction consists of three-spin interactions (one hidden unit and two visible units). Similar architectures were discussed in the machine learning literature under the name three-way Boltzmann Machines Memisevic and Hinton (2007, 2010); Krizhevsky and Hinton (2010). In light of the translational invariance of the Ising model (2) we use the same connection weight $W$ and bias $b$ for all the links. Therefore the BM energy function Eq. (3) only contains two free parameters.

To perform generative sampling of the BM (3) we proceed in two steps by exploiting its particular architecture shown in Fig. 2. First, given a set of visible Ising spins, we can readily perform direct sampling of the hidden units. This is because the conditional probability factorizes into products over each link $p(\mathbf{h}|\mathbf{s})=p(\mathbf{s},\mathbf{h})/p(\mathbf{s})=\prod_{\ell}p(h_{\ell}|\mathbf{s})$, where

[TABLE]

and $\sigma(z)=1/(1+e^{-z})$ is the sigmoid activation function. As shown in Fig. 2(a) the inactive hidden units (white squares) divide the lattice into disconnected components since $h_{\ell}=0$ in Eq. (3) decouples the visible Ising spins reside on the link $\ell$. One can therefore identify connected components using union-find algorithm Sedgewick and Wayne (2011); Gubernatis et al. (2016) and flip all the visible Ising spins within each component collectively at random. This cluster move respect the statistical weight of the BM Eq. (3) due to the $Z_{2}$ symmetry of visible Ising spins in the energy function.

Combining the two steps in Fig. 2(b) forms an update of the visible units of the BM. Recommending the update to the Ising model Monte Carlo simulation, it is accepted with the probability Eq. (1) SM . Matching the unnormalized marginal distribution of the BM (3) $p(\mathbf{s})=\prod_{\ell}\left(1+e^{W\prod_{i\in\ell}s_{i}+b}\right)$ and the Ising model Boltzmann weight Eq. (2) gives rise to a rejection free Monte Carlo scheme. The resulting condition

[TABLE]

can always be satisfied with appropriate chosen $W$ and $b$. It is instructive to examine the BM recommended updates in two limiting cases.

In the limit of $b\rightarrow-\infty$, the solution of Eq. (5) reads $W+b=\ln(e^{2\beta J}-1)$. Thus, the conditional sampling of the hidden units Eq. (4) will set $h_{\ell}=1$ with probability $\sigma(W+b)=1-e^{-2\beta J}$ if the link connects to two parallel spins $\prod_{i\in\ell}s_{i}=1$. While it will always set the hidden unit to inactive $h_{\ell}=0$ if the link connects to anti-parallel spins. Combined with the random cluster flip of visible units, this BM recommended update shown in Fig. 2(b) exactly reproduces the Swendsen-Wang cluster algorithm Swendsen and Wang (1987) of the Ising model.

While in the opposite limit $b\rightarrow\infty$, the solution of Eq. (5) approaches to $W=\beta J$. In this limit, all the hidden units are frozen to $h_{\ell}=1$ because the activation function in Eq. (4) saturates no matter whether the visible Ising spins are aligned or not. The BM Eq. (3) then trivially reproduces the Ising model statistics by copying its coupling constant $\beta J$ to the connection weight $W$. In this limit the BM recommended update shown in Fig. 2(b) is a trivial global flip of the visible Ising spins.

In between the above two limiting cases, the BM still recommends valid rejection free Monte Carlo updates for the Ising model. These updates correspond to the Niedermayer’s cluster algorithm Niedermayer (1988) where the sites are randomly connected into clusters according to Eq. (4) and the clusters may contain misalligned visible spins. The bias parameter $b$ in Eq. (4) controls the activation threshold of the hidden units and thus affects the average cluster size. In essence, the BM (3) encompasses several general cluster Monte Carlo frameworks including the Kandel-Domany Kandel and Domany (1991) and the dual Monte Carlo Kawashima and Gubernatis (1995); Higdon (1998) algorithms. In these algorithms, the Monte Carlo sampling alternates between the physical degrees of freedom and auxiliary graphical variables corresponding to the hidden units of the BM.

Example: Ising model with plaquette interactions–

The potential of BM goes beyond reproducing existing algorithmic frameworks Swendsen and Wang (1987); Niedermayer (1988); Kandel and Domany (1991); Kawashima and Gubernatis (1995); Higdon (1998). By further exploiting its power from latent representations one can make nontrivial algorithmic discoveries. We illustrate this using the plaquette Ising model Liu et al. (2017) as an example. The Boltzmann weight reads

[TABLE]

where the second term contains four-spin interactions on each square plaquette denoted by $\wp$. We consider $K>0$ for concreteness. Since no simple and efficient cluster algorithm is known, Ref. Liu et al. (2017) fits the Boltzmann weight Eq. (6) to an ordinary Ising model Eq. (2) and propose Monte Carlo updates by simulating the latter model with cluster algorithms Wolff (1989); Barkema and Marko (1993). However, the acceptance rates decrease for large systems due to imperfect fittings and the approach ends up to show similar scaling behavior as the single spin flip update.

Here we construct a BM which suggests an efficient, unbiased, and rejection free cluster Monte Carlo algorithm for Eq. (6). First, we decompose the four-spin plaquette interaction using the Hubbard-Stratonovich (HS) transformation Stratonovich (1957); Hubbard (1959)

[TABLE]

where $W=\mathrm{acosh}(e^{2\beta K})$ is the coupling strength between the binary HS field $h_{\wp}$ and the sum of two-spin products $\mathcal{F}_{\wp}(\mathbf{s})=\prod_{i\in\ell_{\wp}}s_{i}+\prod_{i\in\bar{\ell}_{\wp}}s_{i}$ defined for the plaquette. The two parallel links $\ell_{\wp}$ and $\bar{\ell}_{\wp}$ constitute the plaquette $\wp$, see Fig. 3. Equation (7) is equivalent to the discrete HS transformation widely adopted for the Hubbard models Hirsch (1983). Regarding the HS field $h_{\wp}$ as a hidden unit, the following BM

[TABLE]

exactly reproduces Eq. (6) after marginalization. Since Eq. (7) holds for arbitrary partition of the plaquette into two links $\ell_{\wp}\cup\bar{\ell}_{\wp}=\wp$ and $\ell_{\wp}\cap\bar{\ell}_{\wp}=\varnothing$, we choose vertical or horizontal break up at random for each plaquette.

Simulation of the BM Eq. (8) suggests an efficient cluster update for the original plaquette Ising model (6). First of all, sampling the hidden variables given the visible Ising spins is straightforward since the conditional probability factorizes over plaquettes $p(\mathbf{h}|\mathbf{s})=\prod_{\wp}p(h_{\wp}|\mathbf{s})$, where

[TABLE]

Therefore, the hidden unit of each plaquette activates independently given the local feature $\mathcal{F}_{\wp}(\mathbf{s})$. Next, once the hidden variables are given, the BM Eq. (8) corresponds to an Ising model with two-spin interactions only, shown in Fig. 3. One can sample it efficiently using the cluster updates Swendsen and Wang (1987) by taking into account of the randomly modified coupling strengths. As discussed in the above, this amounts to introduce another set of hidden units which plays the role of auxiliary graphical variables. Finally, according to Eq. (1) the updates of the visible Ising spins are always accepted because the BM Eq. (8) exactly reproduces the statistics of the plaquette Ising model Eq. (6).

To demonstrate the efficiency of the discovered cluster update, we simulate the plaquette Ising model (6) in the vicinity of the critical point and compare performance to the simple local update algorithm. Figure 4(a) shows the Binder ratio $\braket{\left(\sum_{i}s_{i}\right)^{4}}/\braket{\left(\sum_{i}s_{i}\right)^{2}}^{2}$ for various system sizes at $K/J=0.2$, which indicates a critical temperature $T/J=2.4955(5)$. The black dashed line indicates the universal critical value of the Binder ratio $1.1679$ corresponding to the two-dimensional Ising universality class Salas and Sokal (2000). Figure 4(b) shows the energy autocorrelation times Ambegaokar and Troyer (2010) of the local updates and the cluster updates at the critical point, both measured in the unit of Monte Carlo sweeps of the visible spins Newman and Barkema (1999). The local updates exhibit the same scaling for the Ising model ($K=0$) and the plaquette Ising model ($K=0.2$). While the cluster updates are orders of magnitude more efficient than the local updates. The dynamic exponent of the cluster algorithm is also significantly reduced compared to the local update.

General framework–

To sum up, we outline a general framework of discovering cluster updates using the following BM

[TABLE]

where $\mathcal{F}_{\alpha}(\mathbf{s})$ is a feature of the visible units and $h_{\alpha}\in\{0,1\}$ is the corresponding hidden variable. $W_{\alpha}$ and $b_{\alpha}$ are connection weight and bias, and $\alpha$ is the index for various features. For example, Eq. (3) and Eq. (8) used the features defined on links ($\alpha=\ell$) and on plaquettes ($\alpha=\wp$) respectively. In general, one is free to design features consist of long-range interactions or even multispins interactions 222Introducing two set of hidden units coupled respectively to the feature $\mathcal{F}_{\ell}(\mathbf{s})=\prod_{i\in\ell}s_{i}$ and $\mathcal{F}_{\wp}(\mathbf{s})=\prod_{i\in\wp}s_{i}$ gives an alternative cluster algorithms for the plaquette Ising model Eq. (6). The resulting algorithm is a simple generalization of the Swendsen-Wang algorithm Swendsen and Wang (1987) to the case of plaquette unit. The satisfied plaquette $\mathcal{F}_{\wp}(\mathbf{s})=1$ is connected with probability $1-e^{-2\beta K}$ while the unsatisfied plaquette is always disconnected. However, this algorithm has larger averaged clusters size compared to the one presented in the texts. This is because the visible units are more likely to be connected with the plaquette decomposition. . There are several crucial points in the design of Eq. (10). First, one can easily sample the hidden units conditioned on these features since there is no interaction between the hidden variables, i.e. Eq. (10) is a semi-restricted BM. The activation probability of each hidden unit is $p(h_{\alpha}=1|\mathbf{s})=\sigma\left(W_{\alpha}\mathcal{F}_{\alpha}(\mathbf{s})+b_{\alpha}\right)$. Second, once the hidden units are given, Eq. (10) reduces to an effective model for the visible spins, which should be easier to sample compared to the original problem. For example, one can randomly flip each disconnected component separated by the inactive hidden units if $E(\mathbf{s})=0$ and $\mathcal{F}_{\alpha}(\mathbf{s})=\mathcal{F}_{\alpha}(-\mathbf{s})$. Alternatively, one can build another BM to simplify the sampling of Eq. (10) given the hidden units. One can even apply this idea iteratively and build a hierarchy of BMs. Overall, the key design principle is to choose appropriate features in Eq. (10) such that the BM correctly reproduce the distribution of the physical problem and is easy to simulate. A good design is likely to exploit the knowledge of the original physical problem 333When all the hidden units are frozen to the active state by infinitely strong bias, the marginal probability of the BM (10) approaches to $p(\mathbf{s})=e^{E(\mathbf{s})+\sum_{\alpha}W_{\alpha}\mathcal{F}_{\alpha}(\mathbf{s})}$. In this limit the BM reduces to an effective Hamiltonian of the visible spins with designed features. One thus loses the flexibility empowered by the latent structures. Moreover, lacking the nonlinearity in the marginal distribution induced by the dynamical hidden units also reduces the fitting capacity of the BM..

The BM Eq. (8) provides a general paradigm to discover efficient cluster updates automatically from data because its parameters are learnable, see Fig. 1. Although we adopted a constructive approach in this paper, in general, the hidden units of the BMs can learn to be auxiliary graphical variables or HS fields. The BM learning can be done in multiple ways, either through unsupervised learning of the configuration data Hinton (2002); Torlai and Melko (2016a), or through supervised learning of the unnormalized target distribution $\pi(\mathbf{s})$ Huang and Wang (2017), or even through reinforcement learning Sutton and Barto (1998) by optimizing the autocorrelation time of the Monte Carlo samples.

In closing, we note many cluster quantum Monte Carlo algorithms Evertz (2003); Kawashima and Harada (2004) share the framework of Kandel and Domany (1991); Kawashima and Gubernatis (1995); Higdon (1998). Generalization of the hidden units to higher integers or even continuous variables is likely to increase the capacity of the BM. One can include higher order self-interactions of the hidden variables in Eq. (10) in this case. To this end, these BMs provide concrete parametrization of valid Monte Carlo update policies which can be optimized through learning. This approach opens a promise of discovering practically useful Monte Carlo algorithms for a broad range of problems, such as frustrated magnets or correlated fermions where known efficient cluster updates are rare. Exploring more general and powerful BM architectures in these settings may lead to even more exciting algorithmic discoveries.

Acknowledgements.

L.W. thanks Li Huang and Yi-feng Yang for collaborations on Huang and Wang (2017); Huang et al. (2016) and acknowledges Jing Chen, Junwei Liu, Yang Qi, Zi-Yang Meng, and Tao Xiang for useful discussions. L.W. thanks Li Huang for comments on the manuscript. L.W. is supported by the Ministry of Science and Technology of China under the Grant No.2016YFA0302400 and the start-up grant of IOP-CAS. The simulation is performed at Tianhe-2 supercomputer in National Supercomputer Center in Guangzhou. We use the ALPS library Bauer et al. (2011) for the Monte Carlo data analysis.

Appendix A Detailed balance condition of Eq. (1)

The acceptance probability of the recommender update from the restricted Boltzmann Machine is derived in Huang and Wang (2017). We repeat the derivation for the BM considered in the main texts for the convenience of the readers.

First of all, the Metropolis-Hastings Metropolis et al. (1953); Hastings (1970) acceptance rate of the physical model satisfies

[TABLE]

The transition probability is recommended from the simulation of the BM, where we sample alternatingly between the hidden and visible units, i.e, $T(\mathbf{s}\rightarrow\mathbf{s}^{\prime})=F_{\mathbf{h}}(\mathbf{s}\rightarrow\mathbf{s}^{\prime})p(\mathbf{h}|\mathbf{s})$. Here $F_{\mathbf{h}}$ denotes update of the visible units given the hidden variables $\mathbf{h}$. In Ref. Huang and Wang (2017) we have used $F_{\mathbf{h}}(\mathbf{s}\rightarrow\mathbf{s}^{\prime})=p(\mathbf{s}^{\prime}|\mathbf{h})$ as the conditional probability given the hidden units of the restricted BM is simply tractable. For the BMs consider in this paper, we adopted more sophisticated update such as cluster flip transition of the visible units. The requirement on $F_{\mathbf{h}}$ is that it respects the joint probability distribution of the BM given the hidden units

[TABLE]

The ratio of the transition probability thus satisfy

[TABLE]

Substitute Eq. (13) into the Metropolis-Hastings acceptance probability Eq. (11), we obtain Eq. (1) in the main texts.