{'fraction_non_alphanumeric': 0.036348222328193715, 'fraction_numerical': 0.013935251417368727, 'mean_word_length': 4.8343122102009275, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerate Monte Carlo Simulations with Restricted Boltzmann Machines

Despite their exceptional flexibility and popularity, the Monte Carlo methods often suffer from slow mixing times for challenging statistical physics problems. We present a general strategy to overcome this difficulty by adopting ideas and techniques from the machine learning community. We fit the unnormalized probability of the physical model to a feedforward neural network and reinterpret the architecture as a restricted Boltzmann machine. Then, exploiting its feature detection ability, we utilize the restricted Boltzmann machine for efficient Monte Carlo updates and to speed up the simulation of the original physical system. We implement these ideas for the Falicov-Kimball model and demonstrate improved acceptance ratio and autocorrelation time near the phase transition point.

Despite their exceptional flexibility and popularity, the Monte Carlo methods often suffer from slow mixing times for challenging statistical physics problems. We present a general strategy to overcome this difficulty by adopting ideas and techniques from the machine learning community. We fit the unnormalized probability of the physical model to a feedforward neural network and reinterpret the architecture as a restricted Boltzmann machine. Then, exploiting its feature detection ability, we utilize the restricted Boltzmann machine for efficient Monte Carlo updates and to speed up the simulation of the original physical system. We implement these ideas for the Falicov-Kimball model and demonstrate improved acceptance ratio and autocorrelation time near the phase transition point.
Monte Carlo method is one of the most flexible and powerful methods for studying many-body systems [1,2]. Its application ranges from the physical sciences [3] including condensed matter physics [4,5], nuclear matter [6] and particle physics [7] all the way to the biological and social sciences [8][9][10]. Monte Carlo methods randomly sample configurations and obtain the answer as a statistical average. However, because of the configuration spaces are exceptionally large for many-body systems, it is typically impossible to perform direct sampling, therefore one resorts to the Markov chain random walk approach to explore the configuration space. In this case, one only needs to know the relative ratio between the probabilities of two configurations.
Designing efficient strategies to explore the configuration space efficiently is at the heart of Markov chain Monte Carlo algorithms. This is, however, a challenging endeavor. Not even mentioning the fundamentally difficult case of glassy energy landscapes, naive Monte Carlo samplings are usually painfully slow close to the phase transitions. These drawbacks motivated noteworthy algorithmic developments in the past decades [11][12][13][14][15][16]. In essence, those algorithms exploit various physical aspects of the problem for efficient Monte Carlo updates. It is however difficult to devise a general strategy to guide optimal Monte Carlo algorithm design.
We address these difficulties in a general setting with insights from the machine learning. Recently, there has been a rising interest in applying machine learning approaches to many-body physics problems. This includes classifying phases of matter [17][18][19][20][21], using the neural networks as variational wave functions [22,23], fitting the density functionals [24][25][26], and solving inverse problems in quantum many-body physics [27,28].
In this paper, we propose a general way to accelerate Monte Carlo simulations of statistical physics problems. We present two algorithmic innovations: a simple supervised learning approach to train the restricted Boltzmann machine (RBM) [29,30] as a proxy of the physical distributions, and an efficient Monte Carlo sampling strategy that exploits the latent structure of the RBM. The RBM is a building block for deep learning and plays an important role in its the recent renaissance [31]. The significance of using the RBM in the Monte Carlo simulations is that it automatically identifies relevant features (such as correlations and collective modes) in the physics model and proposes updates correspondingly with high acceptance rates and low autocorrelations. This approach makes better use of the sampled Monte Carlo data because in addition to estimate the physical observables, the RBM builds an adaptive model for the physical probability distribution and guides better explorations. We illustrate these general ideas using the Falicov-Kimball model [32] as an example. The model describes mobile fermions and localized fermions interact with onsite interactions. The Hamiltonian readŝ where x i ∈ {0, 1} is a classical binary variable representing the occupation number of the localized fermion at site i.ĉ i is the fermion annihilation operator and n i ≡ĉ † iĉ i is the occupation number operator of the mobile fermion. K is the kinetic energy matrix of the mobile fermions. In the following, we consider the model on a periodic square lattice with N sites. Thus, K ij = −t for nearest neighbors and is zero otherwise. The −1/2 offsets in Eq. (1) ensures that both the mobile and localized fermions are half-filled on average. Previous studies show that at U/t = 4 and temperature T /t ≈ 0.15 the system undergoes a phase transition to the checkerboard density wave (CDW) state [33][34][35].
Tracing out the mobile fermions, the occupation number of the localized fermions x ∈ {0, 1} N follow the probability distribution p FK (x) = e −FFK(x) /Z FK , where Z FK is a normalization factor. The negative "free energy" reads (omitting an unimportant constant βU N/4) where β = 1/T is the inverse temperature and for the second term of Eq. (2). The case of classical fields coupled to quadratic fermions represents a broad class of physics problems, including the double-exchange model [36], the mean-field model for phase fluctuated superconductors [37,38] and the Kitaev model after a transformation [39]. Moreover, if one allows imaginarytime dependence in Eq. (2), it covers an even broader class of condensed matter physics problems such as the Hubbard models [40] and spin-fermion models [41].
To compute the physical properties of the Falicov-Kimball model (1), one can perform the Monte Carlo sampling of the classical variables x. To this end, one designs an ergodic strategy to update the variables x and decide whether accept or reject each move. It is sufficient for the Markov chain to converge to the true distribution if the updates satisfy the detailed balance condition [42,43], where T (x → x ) is the proposal probability of an update and A(x → x ) is the acceptance probability of the update. A naive approach would randomly change the classical variables and recompute Eq. (2). To keep the acceptance rate high, one typically applies local updates such as randomly selects a site i and tries to flip the bit We offset the training data so that the minimum value is at zero.
the system size but also has long autocorrelation times. Several improved update schemes [44][45][46][47] have been developed by exploiting the specific features of the Falicov-Kimball model (2). We next present a general approach to propose efficient Monte Carlo updates T (x → x ) by training and simulating an RBM. The RBM is a classical statistical mechanics system defined by the following energy function where x ∈ {0, 1} N and h ∈ {0, 1} M are binary numbers referred as visible and hidden variables. Besides the local biases a i and b j , the visible and hidden units are coupled by the weights W ij . Crucially, there is no coupling within the visible or hidden units themselves, see Fig. 1(a). The joint probability distribution of the visible and hidden variables follows the Boltzmann distribution p(x, h) = e −E(x,h) /Z, where the partition function Z is a normalization factor.
Given sufficiently large number of hidden units, one can tune the parameters of an RBM to let the marginal distribution p(x) = h p(x, h) approximate any discrete distribution [48,49]. This training task is similar to solving the "inverse Ising problem" in statistical physics [50]. The machine learning community has developed practical approaches to train the RBM by minimizing the negative log-likehood − ln p(x) estimated on the samples drawn from the target distribution [51][52][53]. References [54,55] uses this approach to train the RBM for the thermal distribution of the classical Ising model.
However, training the RBM can be even simpler for the statistical physics problems. Two observations are crucial here. First, in addition to the dataset sampled from the target distribution, we do have access to the unnor-  Viewing the RBM as a functional mapping from x to F (x) in Eq. (5) [56], it is tempting to employ a supervised learning approach to train its parameters by fitting the free energy Eq. (5) to the one of the physical model Eq. (2). Here the second observation is important: since only the relative probability ratio matters in the Markov chain Monte Carlo sampling, one thus only cares about F (x) up to an additive constant. Thus, the intractable partition functions of the RBM and the physical model do not appear in the fitting. The overall constant offset of Eq. (2) and Eq. (5) can be chosen at our convenience in the supervised learning. We set up a feedforward neural net for Eq. (5) shown in Fig. 1(b). Its trainable biases and connection weights correspond to the parameters a i , b j and W ij of the RBM. The red hidden neurons activate via the softplus function f (z) = ln(1+e z ). The yellow output neuron sums up the outputs of the hidden neurons and the results coming directly from the input neurons for the final result. Physically, the mobile fermions induce effective interactions between the localized fermions. By training an RBM, we represent the fermionic environment by the classical binary fields represented by the hidden neurons. This supervised training approach is significantly simpler and more efficient compared to the conventional unsupervised learning approach [51][52][53]55].
To collect data and labels for the supervised learning task, we first run single bit-flip simulations of the model (1) to generate 50, 000 independent configurations together with the corresponding negative logprobabilities Eq. (2). In light of the similarity between Eq. (2) and Eq. (5), we set a i = βU/2 and focus on the fitting of b j and W ij . Since sum of the softplus functions Eq. (5) is always positive, we subtract the minimum value in the collected labels to make them nonnegative. This again uses the fact that the fitting is only up to an additive constant. We use 80% of the collected data for training and test on the remaining 20% data to check the generalizability of the fitted neural net. In addition, we apply L2 regulation to the connection weights. This not only prevents overfitting but also makes the sampling of the RBM easier [53]. Figure 2 shows the negative log-probability of the test samples in solid red line and the predictions of the neural net in blue dots. The fitting successfully captures the overall trend of the physical probability distribution. Moreover, the connection weights W ij shown in Fig. 3 also acquire appealing physical meaning by acting as feature detectors for the visible variables. One clearly sees many cross structures corresponding to the staggered density-wave pattern of the localized fermions. There are also several features extend throughout the lattice, indicating that the hidden unit is sensitive to the nonlocal features of the physical model. These learned weights change with the temperature [57], showing that the RBM can pick up the characteristic features of the physics model automatically. When simulating the RBM as a statistical mechanics system, an activated hidden neuron will stimulate the corresponding features in the visible layer. This is similar to the application of RBM to images of hand-written digits. There, the RBM can pick up features like pen strokes and draw new images containing digits [51,52].
After training of the RBM, we use it to generate efficient Monte Carlo updates for the physical variables. To this end, we simulate the RBM using the Monte Carlo method and ensure the updates satisfy the detailed balance condition T (x→x ) T (x →x) = p(x ) p(x) for the visible variables [57]. Therefore, the Metropolis-Hastisings [42,43] solution of the acceptance ratio in Eq. (3) reads Ideally, the acceptance ratio is one if the RBM fits the Falicov-Kimball model perfectly. In this case, one accepts all the proposals from the RBM as was attempted in the Ref. [55]. However, in practice the fitting of the RBM is never perfect given the limited number of hidden units. Equation (6) corrects this error by rejecting unlikely proposals and guarantees exact physical results even with an imperfectly trained RBM. A standard way to simulate the RBM is the block Gibbs sampling approach. Because of the RBM's bipartite architecture, the conditional probability of the hidden variables factorizes where σ(z) = 1/(1 + e −z ) is the sigmoid function. The block Gibbs sampler samples back-and-forth between the hidden and visible layers using Eqs. (7,8), shown in Fig. 4(a). When the simulation of the RBM is much cheaper than the original physical model, one can perform many of these block Gibbs sampling steps before evaluating Eq. (6). The RBM suggests nonlocal updates for the visible variables while still keeping the acceptance ratio high. Moreover, we argue that already a single Gibbs sampling step can be beneficial for the simulation of statistical physics model. Importantly, flipping a hidden variable with the local Gibbs sampling may have nonlocal effects to the physical variables. This is because the hidden neuron may control an extended region of visible variables as shown in Fig. 3. To further encourage this effect, one can perform additional sampling of the hidden variables in between the Gibbs sampling steps, shown in Fig. 4(b). We suggest to change the hidden variables h → h and accept the update with probability min [1, p(h ) p  Fig. 4(a,b) respectively. The critical temperature is at T /t ≈ 0.15 [33][34][35]. The estimated physical observables agree within errorbars for all sampling approaches.
Notice that the visible variables have been traced out in in Eq. (9), sampling according to p(h) captures the generic distribution of the hidden features and is unaffected by the current visible variables. This further improves the sampling of RBM by avoiding the visible and hidden variables to lock each other's feature.
We demonstrate the improvement of using RBM in the Monte Carlo simulation in Fig. 5. The acceptance ratio of the single bit-flip updates simulation of the Falicov-Kimball model (1) decreases monotonously with the temperature lowers, shown in Fig. 5(a). This is because when the system enters into the CDW phase it is harder to add or remove the fermions. In contrast, the acceptance ratio of RBM updates remain high in the whole temperature range across the phase transition. This is because the RBM correctly captures the distribution of the physical system. As a better measure of the improvement, Fig. 5(b) shows the autocorrelation time measured in the unit of Monte Carlo steps per lattice site [1]. The RBM updates reduce the autocorrelation time by at least a factor of two. The scheme of Fig. 4(b) with four additional bit-flip attempts of the hidden units further reduces the autocorrelation time. The overhead of performing the sampling using the RBM is O(M N ), which is negligible compared to the cost of computing Eq. (2) via diagonalizing the fermionic Hamiltonian.
The proposed approach is general. Besides the Falikov-Kimball model and its relatives mentioned after Eq. (2) it is straightforward to use the RBM in Monte Carlo simulations with binary degree of freedoms, such as the Ising and Z 2 gauge fields models, variational [58] and determinantal [40] Monte Carlo simulation of the Hubbard models, and the Fermi bag approach of lattice field theories [59]. For models with continuous variables, one can use the RBM with Gaussian variables [53]. The RBM sampling approach can also be used in combination with the other efficient sampling approaches developed for statistical mechanics problems [11][12][13][14][15][16][44][45][46][47].
To make the presentation cleaner, we divided the computational tasks into three phases: collecting the training data, fitting the RBM, and the actual Monte Carlo simulations. In future, one can use on-line learning to optimize the RBM progressively with newly collected configurations. After trained the RBM in the equilibration phase of the Monte Carlo simulation, one can use it to generate new samples with improved efficiency. One also needs to check the scalability of the proposed approach for larger and more complicated physical systems.
Another future extension is to explore the deep Boltzmann machines [60] and deep belief nets [61] for Monte Carlo simulations. Deeper hierarchical structures may allow even higher level abstraction of the physical degree of freedoms. There were observations indicating that deep structure improves the mixing time of the Monte Carlo sampling [62]. To further exploit the translational symmetry of the physical problems, one may consider to use the shift invariant RBM [22,63] or the convolutional RBM [64][65][66].
Last but not least, in this paper we view the RBM as a generative model and use it to accelerate the Monte Carlo simulation of physical systems. On the other hand, along the line of Ref. [18], it is also interesting to view the RBM as an unsupervised feature detector and explore the patterns in the weights and the latent variables for discovery of physical knowledge.
We end by noting an independent study [67]. It applies similar ideas to the Ising model.  The learned weights change drastically near the critical temperature. Figure 6 shows the weights learned by the RBM at lower temperature. Compared to Fig. 3 at T /t = 0.15, there are more hidden neurons controlling extended regions of the visible variables, indicating enlarged correlation length at lower temperature. The checkerboard pattern of the low temperature phase is also more visible.

Proof of the detailed balance conditions
We prove the simulation of the RBM shown in Fig. 4 of the main texts satisfies the detailed balance condition.
For the case of block Gibbs sampling shown in Fig. 4(a) For the case of Gibbs sampler with additional Metropolis steps for the hidden variables shown in Fig. 4 For the second equality we use that the Metropolis update of the hidden variable satisfies T (h → h )/T (h → h) = p(h )/p(h). This proof generalizes to compositions of several of such updates.