Transforming Generalized Ising Model into Boltzmann Machine

TL;DR

This paper presents an exact method to convert generalized Ising models with many-spin interactions into equivalent Boltzmann machines with only two-spin interactions, improving simulation efficiency.

Contribution

The authors develop a precise algebraic mapping using star-triangle and decoration-iteration transformations, enabling simpler models and enhanced Monte Carlo simulations.

Findings

Mapping reduces model complexity to two-spin interactions.

Application with Swendsen-Wang algorithm decreases critical slowing down.

Effective in models with multi-spin interactions on Kagomé lattice.

Abstract

We find an exact mapping from the generalized Ising models with many-spin interactions to equivalent Boltzmann machines, i.e., the models with only two-spin interactions between physical and auxiliary binary variables accompanied by local external fields. More precisely, the appropriate combination of the algebraic transformations, namely the star-triangle and decoration-iteration transformations, allows one to express the model in terms of fewer-spin interactions at the expense of the degrees of freedom. Furthermore, the benefit of the mapping in Monte Carlo simulations is discussed. In particular, we demonstrate that the application of the method in conjunction with the Swendsen-Wang algorithm drastically reduces the critical slowing down in a model with two- and three-spin interactions on the Kagom\'e lattice.

Equations

$$\pi (\mathbf{\sigma })$$

$$\tilde{\pi }(\mathbf{\sigma },\mathbf{h})$$

$$\pi (\mathbf{\sigma })$$

$$\tilde{\pi }(\mathbf{\sigma },\mathbf{h},\mathbf{d})$$

$$e^{J{\sigma }_1{\sigma }_2}=\Delta \sum _{h=\pm 1}\exp \left[W({\sigma }_1+\mathrm{sgn}(J){\sigma }_2)h\right],$$

$$W=\mathrm{arc}\cosh \left(e^{2|J|}\right)/2.$$

$$\exp \left[M{\sigma }_1{\sigma }_2{\sigma }_3+J_1{\sigma }_2{\sigma }_3+J_2{\sigma }_3{\sigma }_1+J_3{\sigma }_1{\sigma }_2\right]$$

$$=\Delta \sum _{h=\pm 1}\exp \left[\sum _{i=1}^3(W_{i}h+a_i){\sigma }_i+bh\right],$$

$$\exp (4{\chi }_i{a}_i)=\frac{\sinh (2(|J_i|+M))}{\sinh (2(|J_i|-M))},$$

$$\cosh (2W_i)=\frac{e^{2|J_i|}\cosh \left(2(|J_j|+|J_k|)\right)-e^{-2|J_i|}\cosh \left(2(|J_j|-|J_k|)\right)}{[2\cosh (4|J_i|)-2\cosh (4M){]}^{1/2}},$$

$$\sinh (2b)=\frac{-\sinh (2{\chi }_i{W}_i)\sinh (4M)}{[(\cosh (4J_j)-\cosh (4M))(\cosh (4J_k)-\cosh (4M)){]}^{1/2}},$$

$$|M|<|J_i|$$

$$\text{sgn}(J_1{J}_2{J}_3)$$

$${\pi }_{|\mathcal{C}|}({\mathbf{\sigma }}_{\mathcal{C}};M):=\exp \left(M\prod _{j\in \mathcal{C}}{\sigma }_j\right),$$

$$\pi (\mathbf{\sigma })=\prod _{\mathcal{C}}{\pi }_{|\mathcal{C}|}({\mathbf{\sigma }}_{\mathcal{C}};M_{\mathcal{C}})=\prod _{\mathcal{C}}\left(\sum _{{\mathbf{h}}_{\mathcal{C}},{\mathbf{d}}_{\mathcal{C}}}{\tilde{\pi }}_{\mathcal{C}}({\mathbf{\sigma }}_{\mathcal{C}},{\mathbf{h}}_{\mathcal{C}},{\mathbf{d}}_{\mathcal{C}})\right),$$

$${\pi }_3(\mathbf{\sigma };M)$$

$${\pi }_3({\mathbf{\sigma }}_{{\mathcal{C}}_1};M){\pi }_3({\mathbf{\sigma }}_{{\mathcal{C}}_2};M)$$

$${\pi }_4(\mathbf{\sigma };M^{(0)})$$

$$M^{(1)}=\mathrm{arc}\cosh \left(e^{2|M^{(0)}|}\right)/2.$$

$${\pi }_4(\mathbf{\sigma };M^{(0)})$$

$$\times {\tilde{\pi }}_1(\mathbf{\sigma },\mathbf{h},d),$$

$${\tilde{\pi }}_2^{\mathrm{STT}}(\mathbf{\sigma },\mathbf{h},d)$$

$${\tilde{\pi }}_2^{\mathrm{DIT}}(\mathbf{\sigma },\mathbf{h},d)$$

$$\displaystyle\ \ \ \ \ \ \ \ +W^{\prime}(\sigma_{1}-\sigma_{2})h_{5}+W^{\prime}(\sigma_{3}-\sigma_{4})h_{6}\Bigg{]},$$

$${\tilde{\pi }}_1(\mathbf{\sigma },\mathbf{h},d)$$

$$\exp (4a)$$

$$\cosh (2W)$$

$$\sinh (2b)$$

$$n_h^{\mathrm{STT}}=3(p-2),n_h^{\mathrm{DIT}}=p-2,n_d=p-3.$$

$$H=-\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j,$$

$$-\beta H(\mathbf{\sigma })=\beta \sum _{⟨i,i^{\prime }⟩\in E}{\sigma }_i{\sigma }_{i^{\prime }}+\sum _{△}{M}_{△}{\tau }_{△}+\sum _{▽}{M}_{▽}{\tau }_{▽},$$

$$\pi (\mathbf{\sigma })$$

$$=$$

$$\cosh (2W_{△(▽)})$$

$$\sinh (2b)$$

$$\tilde{\tilde{\pi }}(\mathbf{h})$$

$$W_h=\mathrm{arc}\cosh \left(e^{2W}\right)/2.$$

Full text

I Introduction

Inspired by the tremendous success in machine learning fields such as the image or speech recognition, climate prediction, and market analysis, the increasing number of studies have shown the validity of artificial neural network in physical problems. For instance, physicists have recognized that the phase classification can be done in a parallel fashion to the ordinary machine learning tasks by using the data of the classical, quantum, or auxiliary degrees of freedom Carrasquilla and Melko (2017); Wang (2016); Torlai and Melko (2016); Ohtsuki and Ohtsuki (2016, 2017); Yoshioka et al. (2018); Nieuwenburg et al. (2017); Ch’ng et al. (2017); Wetzel (2017); Huembeli et al. (2018); Araki et al. (2019); Rem et al. (2018). Another important direction focuses on the representability itself. The applications such as the variational wave functions describing ground states turn out to outperform the state-of-the-art methods Carleo and Troyer (2017); Nomura et al. (2017); Glasser et al. (2018); Carleo et al. (2018); Freitas et al. (2018), and are also valid in estimating the Boltzmann factor of a system to accelerate the classical Monte Carlo (MC) simulations Torlai and Melko (2016); Huang and Wang (2017). Furthermore, numerous studies reveal the ability to precisely express nonlinear functions such as the wave functions of stabilizer states and chiral topological states Deng et al. (2017a, b); Gao and Duan (2017); Chen et al. (2018); Jia et al. (2019); Wang (2017). This is not only insightful from the perspective of exploring soluble models but also fruitful in constructing new numerical techniques including optimization of variational functions, MC simulations, and so on.

Since the invention of the MC method, physicists have long made efforts to develop versatile and efficient simulation methods to investigate statistical models. In classical lattice systems, the single-spin flip (SSF) algorithm is undoubtedly one of the most widely-used techniques, as it is model independent. The locality of the variables involved in a single update procedure, however, inevitably leads to a severe slowing down near critical points or at low temperature for non-ordering systems. One of the solutions is to apply the global updates such as the cluster algorithms Swendsen and Wang (1987); Wolff (1989), worm algorithm Prokof’ev and Svistunov (2001), and loop algorithm Evertz (2003), but they are mostly restricted to two-body interacting systems. Developing a generic technique applicable to a wide variety of systems involving many-body interactions is highly challenging.

In this work, we establish a mapping from the generalized Ising model to the Boltzmann machine (BM). The former, which includes many-spin interactions, is used to describe magnetic and thermodynamic properties of solids, the effective model of alloys, spin glass models, and so on Biroli and Mézard (2001); Blum et al. (2005); Berthier and Biroli (2011). The BM, on the other hand, is an expression of probability distribution by two-spin interactions between physical (visible) and auxiliary (hidden) binary degrees of freedom. The mapping procedure goes as follows: one first decomposes a $p$-spin Ising interaction into a sum of three- and $(p-1)$-spin interactions by adding an auxiliary spin. This is repeated for the latter many-spin interaction until the total interaction is expressed in terms of two- and three-spin interactions. Finally, each three-spin interaction is transformed into two-spin interactions and single-spin terms. The idea of such consecutive decomposition can also be seen in the context of quantum annealing Biamonte (2008); Lechner et al. (2015); Leib et al. (2016); Choi (2011).

After constructing the rigorous mapping , we take advantage of it by presenting a novel global update scheme; the application of the Swendsen-Wang (SW) algorithm to the exactly surrogate BM. By comparing the autocorrelation times of the visible spin configurations, we demonstrate in a model with two- and three-spin interactions on the Kagomé lattice that the sampling efficiency of the cluster algorithm on the BM at the critical temperature is drastically improved compared to that of SSF performed on the original Hamiltonian. While the mapping introduced in this paper is applicable to the generalized Ising model in any number of dimensions, our results are presented in the above-mentioned model for clarity and simplicity.

The remainder of this paper is organized as follows. In Sec. II we first introduce the most primary transformation techniques, namely the decoration-iteration and star-triangle transformations. Embedding the hidden spins by combining the two, an arbitrary many-spin interaction is shown to be mappable to the BM. Furthermore, multiple interaction terms can be taken into account by simply considering the embedding procedures independently. In Sec. III, we see that the transformation is numerically beneficial since the existing cluster update can be applied to the BM. Finally, the summary for the current work and discussion concerning the future direction are given in Sec. IV. For the completeness of the paper, in Appendix A we derive the explicit expression of the BM mapped from the generalized Ising model. In Appendix B, we qualitatively see that the block Gibbs sampling for the RBM equivalent to the ferromagnetic Ising model is highly inefficient. The partition function of the pure three-spin interacting model on the Kagomé lattice, from which the absence of the phase transition follows, is calculated in Appendix C. Also, two schemes to perform the cluster update under magnetic field are given in Appendix D, and the observation of the physical quantity in the extended space is discussed in Appendix E.

II Algebraic Transformation of Boltzmann Factors

In this section, we find the equivalence of the generalized Ising models and the BMs with specific architectures, i.e., the restricted Boltzmann machine (RBM) and the deep Boltzmann machine (DBM). [See Fig. 1 for the graphic representation.] First, we state the definition of the RBM and DBM, which give some probability distributions of visible spins with auxiliary binary degrees of freedom, or the hidden spins. The Boltzmann factor of an RBM is given as Smolensky (1986); Hinton (2002)

[TABLE]

where $\pi({\bm{\sigma}})$ and $\widetilde{\pi}({\bm{\sigma}},{\bm{h}})$ are the Boltzmann factor of $N_{v}$ visible spins and the RBM with additionally $N_{h}$ hidden spins, respectively. Here, $\sigma_{i}$ and $h_{j}$ are the $i$-th visible and $j$-th hidden spins that take either $+1$ or $-1$, coupled via the interaction $W_{ij}$. The local external field is given as $a_{i}$ and $b_{j}$ for visible and hidden spins, respectively. We denote the visible and hidden spin configurations by ${\bm{\sigma}}:=(\sigma_{1},\sigma_{2},...,\sigma_{N_{v}})\in{\cal S}$ and ${\bm{h}}:=(h_{1},h_{2},...,h_{N_{h}})\in{\cal H}$, where ${\cal S}=\{-1,1\}^{N_{v}}$ and ${\cal H}=\{-1,1\}^{N_{h}}$ are the sets of all possible binary spin configurations for visible and hidden spins, respectively. Also, to discriminate between the spaces with and without the hidden spins, we call $\cal S$ as the “original space” and ${\cal S}\cup{\cal H}$ as the “extended space.” Note that the absence of the intralayer couplings is reflected in the bipartite structure as is depicted in Fig. 1(a). The BM without any restriction on the connectivity is not considered in the following.

A crucial notion in introducing the DBM is the “layer,” for which we provide a description based on the graphical understanding of Fig. 1. As aforementioned and also obvious from Fig. 1(a), no coupling is present between the visible or hidden spins in the RBM. Such a bipartite structure is understood as a BM with the “visible layer and one hidden layer.” Namely, the visible layer contains all the visible spins, and the hidden layer consists of hidden spins that are interacting with one or more visible spins. We may analogously construct another hidden layer on top as depicted in Fig. 1(b), which is a BM with the “visible layer and two hidden layers” in turn. In principle, one may add as many hidden layers as desired. Such an architecture is generally referred to as the DBM. The maximum number of layers obtained by mapping from the generalized Ising model is two as we discuss in Sec. II.2.3, and therefore we discriminate the first and second hidden layers as the “hidden layer” and “deep layer” in the following. The spins included in these layers are correspondingly referred to as the “hidden spins” and “deep spins.” To be concrete, the Boltzmann factor of a DBM is given as

[TABLE]

where $\widetilde{\pi}({\bm{\sigma}},{\bm{h}},{\bm{d}})$ is the Boltzmann factor for each spin configuration, $d_{k}$ is the $k$-th deep spin with the local field $b^{\prime}_{k}$, and $W^{\prime}_{jk}$ is the interaction between the $j$-th hidden spin and the $k$-th deep spin. As is the case for the visible and hidden spins, a configuration of the deep spins is denoted by ${\bm{d}}:=(d_{1},d_{2},...,d_{N_{d}})\in{\cal D}$, where ${\cal D}=\{-1,1\}^{N_{d}}$ is the set of all possible deep spin configurations. The union ${\cal S}\cup{\cal H}\cup{\cal D}$ is also referred to as the extended space in the following.

The probabilistic network given by quadratic terms as Eqs. (2) and (4) is known to be powerful to express non-linear functions and in fact applied widely in the field of machine learning, condensed matter physics, quantum physics, and so on Deng et al. (2017a, b); Gao and Duan (2017); Chen et al. (2018); Jia et al. (2019); Carleo and Troyer (2017); Nomura et al. (2017); Glasser et al. (2018). In the remainder of this section, we find the exact mapping from the generalized Ising model to the RBM or DBM.

II.1 Transformation techniques

We introduce two mapping techniques to embed hidden spins as is graphically described in Fig. 2: the decoration iteration transformation (DIT) and star-triangle transformation (STT). Note that the newly embedded spins are auxiliary, and the original interactions are realized by tracing out such degrees of freedom.

The DIT, depicted in Fig. 2(a), is a very simple transformation which embeds a hidden spin $h$ between two interacting visible spins as follows Fisher (1959); Wegner (1971),

[TABLE]

where $\Delta$ is the normalization factor. The new interaction $W$ is given as

[TABLE]

Since the DIT can be carried out for any $J$, an arbitrary Ising model with two-spin interactions including random spin-glass, frustrated system, and fully-connected models can be mapped into an equivalent RBM. Application of such transformation technique allows one to obtain the exact solution for a model on, for instance, two-spin interacting Ising model on a bond-decorated lattice that can be transformed into the soluble model on an undecorated lattice Syozi (1951, 1972); Fisher (1959); Antonosyan et al. (2009); Rojas et al. (2009); Lisnyi and Strečka (2014).

The other technique, known as the STT, embeds a hidden spin $h$ into three visible spins with both two- and three-spin interactions as is illustrated in Fig. 2(b) Wegner (1971); Wu (1990); Lin and Wu (1990); Lu and Wu (2005); Caravelli and Nisoli (2018). Expressed in the form of the Boltzmann weight, this can be written as Lin and Wu (1991)

[TABLE]

where $M$ and $J_{i}$ are the amplitudes of three- and two-spin interactions, respectively. The interaction between the visible and the hidden spins in the extended space is denoted by $W_{i}$, and the local fields are denoted by $a_{i}$ and $b$, respectively. It can be shown that Eq. (II.1) amounts to eight nonequivalent conditions which yields the solutions as

[TABLE]

where $\chi_{i}={\rm sgn}(J_{i})$ is the sign of the two-spin interaction. The subscripts in Eqs. (9) and (LABEL:STT_b), i.e., $i,j,$ and $k$, must be chosen such that none of them are identical to each other. Importantly, the STT is valid under the following conditions com

[TABLE]

The notion of the DIT and STT can be generalized to include many-spin interactions. A system in the original space with both four- and two-spin interactions, for instance, can be mapped to a model with three-spin interactions that involves a single hidden spin. This mapping, used to obtain the solution of the zero-field eight-vertex model, is referred to as the star-square transformation Baxter (1971, 1972); Wegner (1971); Strečka (2018). While the transformation from the extended space into the original space, known as “the star-polygon transformation,” is achieved by tracing out the hidden spins and is in general tractable Rojas et al. (2009); Strečka (2010), its inverse mapping exists in very limited cases. We note in passing that the DIT and STT have been extended to models with local quantum degrees of freedom such as Heisenberg spins and itinerant electrons Čanová et al. (2008); Strečka et al. (2008, 2009).

II.2 Generalized Ising model as Boltzmann Machine

Here, we show that the generalized Ising model, which consists of many-spin interactions, can be mapped to an equivalent model with two-spin interactions and local fields. First, let us consider some bare many-spin interaction within the model. The Boltzmann factor is given as

[TABLE]

where $M$ is the amplitude of the interaction. Here, the set of sites of the visible spins involved, dubbed as the “cell” in the following, is denoted as ${\cal C}$. Correspondingly, a visible spin configuration of the cell is denoted by ${\bm{\sigma}}_{\cal C}\in{\cal S_{\cal C}}$, where ${\cal S}_{\cal C}$ is the set of all possible configurations of the visible spins included in $\cal C$. The subscript of the Boltzmann factor, $|{\cal C}|$, is the number of the visible spins included in the cell. Note that arbitrary hidden and deep spins can be uniquely labeled due to the construction procedure we introduce in the following. It is straightforward to see that the Boltzmann weight of the original system can be written as

[TABLE]

where $\widetilde{\pi}_{\cal C}({\bm{\sigma}}_{\cal C},{\bm{h}}_{\cal C},{\bm{d}}_{\cal C})$ is the weight of the RBM obtained by transformation. Note that although a visible spin may be included in multiple cells, the auxiliary spins are uniquely allocated to their corresponding cells by construction. A single realization of hidden and deep spins labeled by ${\cal C}$ is denoted as ${\bm{h}}_{\cal C}\in{\cal H_{\cal C}}$ and ${\bm{d}}_{\cal C}\in{\cal D_{\cal C}}$, where ${\cal H}_{\cal C}$ and ${\cal D_{\cal C}}$ are the sets of all possible configurations of the hidden and deep spins included in $\cal C$, respectively.

As we see from Eq. (14), the Boltzmann factor of any model can be decomposed into products over cells and the independence of hidden spins holds. For simplicity, we focus on bare $p$-spin interactions in the following.

II.2.1 Three-spin interaction as RBM

First, we discuss the three-spin interaction, which is the simplest possible case. The STT cannot be applied straightforwardly to a bare three-spin interaction due to the conditions given in Eqs. (11) and (12). To avoid this problem, we introduce the “virtual two-spin interactions” that cancel each other out as follows,

[TABLE]

As is depicted in Fig. 3(a), this can be mapped into the RBM by applying Eq. (II.1) to the first and subsequently Eq. (5) to the second factors. Note that the amplitudes of the virtual two-spin interaction, $J_{i}$, can be taken arbitrarily as long as Eqs. (11) and (12) are satisfied.

Next, let us consider two sets of interacting spins as denoted in Fig. 3(b). Although the naive application of Eq. (15) yields eight hidden spins, two of them on the shared edges can be eliminated by modifying the signs of virtual two-spin interactions. For instance, by considering the virtual two-spin interactions with a homogeneous amplitude, we obtain

[TABLE]

from which six hidden spins emerge.

Furthermore, a system with pure three-spin interaction on the triangular lattice, known to be exactly soluble and dubbed as the Baxter-Wu model Baxter and Wu (1973), can also be mapped into the RBM merely without hidden spins generated by the DIT. In other words, we may choose the signs of the virtual two-spin interactions on the triangles carefully so that the sum at each edge would cancel out, resulting in a reduced number of the hidden spins. Still, there are exponential number of ways to represent this model by tuning the amplitudes and the signs of the virtual two-spin interactions. Shown in Fig. 3(c) is the mapping with 4-fold periodicity along the $x-$axis.

II.2.2 Four-spin interaction as DBM

Next, we show that four-spin interaction can be expressed by introducing the second hidden layer, or the “deep” layer. The illustration of the two-step transformation is shown in Fig. 4 (a). In the first step, we interpret the product of Ising variable as a single new binary variable and apply the DIT as

[TABLE]

where $\pi_{4}({\bm{\sigma}};M^{(0)})$ is the Boltzmann factor for four-spin interaction with the amplitude $M^{(0)}$. The interaction amplitudes $M^{(0)}$ and $M^{(1)}$ are related by Eq. (6) as

[TABLE]

Next, we use Eq. (15), the result for the three-spin interaction, under homogeneous virtual two-spin interactions $J>|M^{(1)}|$ for simplicity. Here, we obtain the expression as

[TABLE]

where

[TABLE]

Here, $\widetilde{\pi}_{2}^{\rm STT({\rm DIT})}({\bm{\sigma}},{\bm{h}},d)$ is the Boltzmann factor for the two-spin interaction in the extended space obtained by the STT (DIT), and $\widetilde{\pi}_{1}({\bm{\sigma}},{\bm{h}},d)$ is that for the local external fields. The interaction $W$ $(W^{\prime})$ is the interlayer coupling introduced by the STT (DIT). The magnetic field for the visible, deep, and hidden spins are denoted by $a,a_{d},$ and $b$, respectively. The equivalence of amplitudes of the three-spin interactions, $M^{(1)}$, leads to $a_{d}=2a$. Also, the application of DIT to the virtual two-spin interaction yields $W^{\prime}={\rm arc}\cosh(e^{2J})/2.$ Note that the other parameters are obtained by substituting $M=M^{(1)}$ and homogeneous $J_{i}=J$ into Eq. (II.1) as

[TABLE]

Another way to transform four-spin interaction is to apply the star-square transformation. Although it also requires a single deep spin, the architecture of the hidden spins would be symmetric and hence different from the aforementioned transformation. To keep the number of auxiliary spins minimum, we will not use the star-square transformation in the following.

Note that $N_{h}/N_{d}$, or the ratio of the number of the hidden spins to that of the visible spins, may be reduced in a larger system as well as in Sec. II.2.1; we may cancel out the virtual two-spin interactions by modifying their signs.

II.2.3 $p$-spin interaction as DBM

The discussion for four-spin interaction can be extended to $p$-spin interaction. Assume that Ising interactions up to $(p-1)$-spin terms can be mapped into a BM. The DIT splits the original model into two $(p/2+1)$-spin terms for even $p$ and $((p+1)/2+1)$- and $((p-1)/2+1)$-spin terms for odd $p$. The mapping exists for $p=3,4$ as we have shown in the previous sections, hence for arbitrary $p>4$ as well. For further discussion and explicit expression in the form of BM, see Appendix A.

The number of the hidden and deep spins can be computed as well. We denote the numbers of hidden spins in the first layer generated by the STT, DIT, and that of the deep spins in the second hidden layer as $n_{h}^{\rm STT},n_{h}^{\rm DIT},$ and $n_{d}$, which are shown in blue, red, and black filled circles in Fig. 4, respectively. It is straightforward to show that the following relations are satisfied,

[TABLE]

It is noteworthy from the perspective of the numerical cost that the number of hidden and deep spins increases only linearly with respect to $p$. The computational order remains to be the same for any update scheme.

III Monte Carlo sampling on Boltzmann machine

In the following, we utilize the BM obtained by the transformation to classical MC sampling. The block Gibbs sampling, which is frequently used in the machine learning community, turns out to lack efficiency in terms of autocorrelation although the numerical cost per a single MC step is low Goodfellow et al. (2016); Huang and Wang (2017). We alternatively apply the Swendsen-Wang algorithm to the extended space and demonstrate the speed up compared to the SSF on the original space. First, we take the square-lattice Ising model for a simplified description of our scheme, and then proceed to show the results in the model with ferromagnetic two-spin interactions and alternating three-spin interactions on the Kagomé lattice, which is one of the the most comprehensible models that includes the many-spin interaction and also suffers from the slowing down at the critical temperature.

III.1 Ising model on square lattice

We consider a Hamiltonian with two-spin interactions on a square lattice,

[TABLE]

which shows the ferro-paramagnetic transition at $T_{\rm c}=2/\ln(\sqrt{2}+1)$ as is widely known in statistical physics Onsager (1944). The study by the SSF algorithm, i.e., the Glauber dynamics in a broader sense, suffers from severe critical slowing down Glauber (1963); Ito (1993). The application of global updates, e.g., cluster updates Swendsen and Wang (1987); Wolff (1989), is one of the solutions under some circumstances.

Note that Eq. (25) itself can be regarded as the RBM since the square lattice is bipartite. Since the block Gibbs sampling method, which fix the configuration in a sublattice to enable independent sampling from another, is not beneficial as is qualitatively discussed in Appendix B, we alternatively apply the well-known Swendsen-Wang algorithm to take advantage of the expression. Although there is substantially no further gain for two-spin interacting model, in the following subsection, we demonstrate the speed up in terms of the autocorrelation time. As is shown in Sec. II, the number of the hidden and deep spins increases linearly with the system size, and therefore the computational order of a single MC step remains to be $O(N)$ in a model with short-range interactions, which is also the case in the current work.

III.2 Generalized Ising model on Kagomé lattice

Here, we consider a model with ferromagnetic two-spin interactions and also three-spin interactions on a Kagomé lattice. Let $E$ be the set of edges and $\triangle(\bigtriangledown)$ be an upward (downward) triangle on the lattice. The Hamiltonian is defined as

[TABLE]

where $\tau_{\triangle}=\prod_{i\in\triangle}\sigma_{i}$ and $\tau_{\bigtriangledown}=\prod_{i\in\bigtriangledown}\sigma_{i}$ is the product of the spin variables, $M_{\triangle}/\beta\ (M_{\bigtriangledown}/\beta)$ denotes the amplitude of the three-spin interactions for upward (downward) triangles, and $\langle i,i^{\prime}\rangle$ is the edge connecting sites $i$ and $i^{\prime}$. The symmetry that combines the spin inversion and the mirror inversion is present when $M_{\triangle}+M_{\bigtriangledown}=0$, and the model exhibit second-order transition at finite temperature for finite $M_{\triangle}$ [See Appendix C for discussion when only three-spin interactions are present.]. This can be understood from the correspondence between the present model and the antiferromagnetic Ising model with a uniform external field on the honeycomb lattice Wu et al. (1989).

In the following, we assume $M_{\triangle}=-M_{\bigtriangledown}=M>0$. The Boltzmann weight is transformed as

[TABLE]

where the parameters in the extended model are obtained by substituting $M_{\triangle(\bigtriangledown)}$ and $\beta$ into the STT, or Eq. (II.1). Note that the external fields on the visible spins is absent due to the cancellation caused by the alternating signs of the interactions. As is graphically described in Fig. 5, one hidden spin is embedded per triangle and denoted as $h_{\triangle}$ or $h_{\bigtriangledown}$. The explicit expressions for the parameters can be read off from

[TABLE]

where $W_{\triangle}=-W_{\bigtriangledown}=W>0$ and the signs of the parameters satisfy ${\rm sgn}(W_{\triangle(\bigtriangledown)})={\rm sgn}(M_{\triangle(\bigtriangledown)})$. Owing to the cancellation of the magnetic field on visible spins, we obtain a simple expression by tracing out the visible spins as

[TABLE]

where $\Delta^{\prime}$ is another renormalization factor, $\overline{E}$ and $\overline{V}$ are the sets of edges and vertices in the honeycomb lattice as is shown in the rightmost panel of Fig. 5. Note that the interaction is antiferromagnetic, reflecting the alternating signs of the three-spin interactions in the original space. The amplitude of the two-spin interaction is obtained by the DIT as

[TABLE]

Although the model defined by Eq. (31) is not soluble at $b\neq 0$, an approximate solution of the transition point can be obtained by imposing some assumption after mapping the original model to an eight-vertex model Wu et al. (1989). This turns out to be in a fairly good precision but not exact, and hence we determine the transition point from the finite size scaling of the Binder ratio Binder (1981); Tanaka et al. (2017).

As was introduced by Binder, the renormalization group theory leads us to assume the scaling law for the Binder ratio,

[TABLE]

where $m$ is the magnetization per site and $\langle\cdots\rangle$ denotes the thermal average. The scaling of this quantity in the vicinity of the critical temperature, $T_{\rm c}$, is given as follows,

[TABLE]

where $L$ is the linear system size, $\nu$ is the critical exponent for the correlation length, and $F$ is an appropriate polynomial function, which has been taken to be cubic in this work. Shown in Fig. 6 are the results for ferromagnetic two-spin and alternating-sign three-spin interacting model on the Kagomé lattice with $M/\beta=0.1$. From the data collapse in Fig. 6, we see that the scaling analysis is valid. The critical exponent $\nu$ is confirmed for numerous $M/\beta$ to be in good agreement with $\nu=1$, suggesting that the transition falls into the two-dimensional Ising universality class Fisher (1969). The simultaneously estimated quantity, i.e., the critical temperature $T_{\rm c}$, is summarized in Fig. 7(a) together with the approximate solution given by Ref. Wu et al. (1989).

At the critical temperature, we use the cluster algorithm in the extended space. The detailed description is given in Appendix D. Our main result is summarized in Fig. 7(b), in which we compare the autocorrelation time $\tau$ of the magnetization measured in units of Monte Carlo steps per site for the whole system. [See Appendix E for the calculation of physical quantities in the extended space. ] The magnetization at $t$-th Monte Carlo step is calculated from the spin configurations ${\bm{\sigma}}(t)=(\sigma_{1}(t),\sigma_{2}(t),...,\sigma_{N_{v}}(t))$ with $N_{v}$ being the number of visible spins as

[TABLE]

The estimation of $\tau$ is done by evaluating the decay of the equilibrium autocovariance Evertz (2003),

[TABLE]

where $\langle\cdots\rangle$ denotes the average over $t_{0}$, namely the MC steps. The critical slowing down is constantly observed in the SSF performed on the original space, while the application of the cluster update to the BM significantly improves the situation. We observe that the dynamical exponent $z$, which is the slope of data in Fig. 7(b), is also reduced, while the possibility that it gradually grows in the larger system sizes cannot be ruled out.

The increase in the autocorrelation time along the three-spin interaction $M/\beta_{\rm c}$ is understood in the following way. The virtual magnetic field induced in the hidden spins by the STT is amplified as $M/\beta_{\rm c}$ is increased, and thus the virtual magnetization per cluster increases, modifying the flipping probability of each cluster to be unbalanced. While the cluster is flipped randomly when the magnetic field is absent, finite-valued Zeeman energy results in unbalanced flipping probability due to the detailed balance condition. [See Appendix D for further discussion.] Such a situation prevents the system from exploring the spin configurations efficiently, and thus show a weaker speed up.

IV Conclusion and Discussion

In the current work, we found an algebraic transformation of the Ising model with many-spin interactions into the Boltzmann machine in which only two-spin interactions and virtual local fields are present. The decoration-iteration and star-triangle transformations were applied to embed hidden and deep spins, namely the auxiliary degrees of freedom to be traced out. At the expense of the dimension of the spin space, significant suppression of the critical slowing down is achieved by applying the cluster algorithm to the Boltzmann machine.

Our scheme is also capable of handling continuous classical spin systems with many-spin interactions. As in the case with two-spin interactions Wolff (1989), we may consider projecting each variable on some axis. Namely, we rewrite a continuous variable ${\bm{S}}_{i}$ on site $i$ by a new Ising variable $\sigma_{i}$ as ${\bm{S}}_{i}=\sigma_{i}|{\bm{S}}_{i}\cdot{\bm{n}}_{i}|{\bm{n}}_{i}+{\bm{S}}_{i}^{\perp}$, where ${\bm{n}}_{i}$ is the randomly chosen projection axis and ${\bm{S}}_{i}^{\perp}$ is orthogonal to ${\bm{n}}_{i}$. We can now apply our method by regarding the model for $\{\sigma_{i}\}$ obtained by projection as the generalized Ising model. The randomness of the projection axis at each Monte Carlo step would assure the ergodicity of the scheme.

Beyond our scope in the current work is the optimal transformation for simulation. The transformation is non-unique when the virtual interaction is required, and we may even consider infinitely strong coupling to express extreme situations such as decoupled or completely aligned pairs of spins. Non-uniqueness arises also when four-spin interaction is present. Although we excluded the application of the star-square transformation for clarity, comparison of the numerical efficiency between different transformations may be worth investigating. Switching into different Boltzmann machines for each step may allow us to explore the free energy landscape more efficiently.

In closing, we would like to note the applicability of the algebraic transformation to wider fields of research. One interesting direction is undoubtedly the pursuit of equilibrium statistical physics, which includes extending and exploring exactly soluble models and replacing the Swendsen-Wang algorithm with other global updates to tackle frustrated systems. Another problem lies in the field of computer science; the application to combinatorial optimization problems; our decompositions applicable also in the finite-temperature case may open a new way to introduce ancilla spins required in the experimental implementation of annealing process.

Acknowledgements

The authors acknowledge helpful discussions with Seiji Miyashita, Tsuyoshi Okubo, Shu Tanaka, Synge Todo, and Lei Wang. This work was supported by JSPS KAKENHI Grant Nos. JP17K14352. JP18H04220, 18K03445, 18H04478. N. Y. was supported by the JSPS through Program for Leading Graduate Schools (ALPS) and JSPS fellowship (JSPS KAKENHI Grant No. JP17J00743).

Appendix A Transformation of $p$-spin interaction

In the following, we discuss the transformation of a bare $p$-spin interaction into the BM. The procedure consists of two steps; embedding the deep spins, and then the hidden spins. A system with a bare $p$-spin interaction is mapped into an equivalent system with three-spin interactions in the former step and subsequently broken into systems with two-spin interactions and local fields in the latter step. Since the DIT is applied repeatedly, we define

[TABLE]

where $M$ determines the amplitude of the interaction as $\pi_{p}({\bm{\sigma}};M)$. The number of DITs applied, or $k$, is referred to the DIT transformation order.

Let us take four-spin interaction as the starting point. As we have shown in the main text, the Boltzmann factor can be expressed by embedding the deep spin as

[TABLE]

By replacing one of the Ising variables with a product of two, we obtain the expression for five-spin interaction as

[TABLE]

which is visually described in Fig. 8(a). Repeating the DIT such that the DIT transformation order is as homogeneous as possible, we can show a posteriori that the general expression is given as

[TABLE]

where $n$ is an integer satisfying $2^{n}\leq p-2<2^{n+1}$. Note that there are two factors consisting of two visible spins, whereas the others contain only one. The number of factors with $M^{(n+1)}$ is zero if $p-2=2^{n}$.

To consider the hidden spins, we substitute $M^{(n)}$ in the transformation introduced in Sec. II.2.1 as, for instance,

[TABLE]

which is described in the right-most panels of Fig. 8.

Appendix B Block Gibbs sampling from RBM equivalent to square Ising model

In the following, we shortly introduce the block Gibbs sampling and discuss its application to the RBM equivalent to the ferromagnetic Ising model on the square lattice. Let $A$ and $B$ be the sublattices of the square lattice. The spin configuration on a sublattice is denoted as $\sigma_{A(B)}=\left\{\{\sigma_{i}\}|i\in V_{A(B)}\right\}$ where $V_{A(B)}$ is the set of $A$($B$)-sublattice sites. The posterior distribution for spin configuration on the $A$ sublattice is written as

[TABLE]

where $\pi(\sigma_{A},\sigma_{B})$ denotes the Boltzmann factor for the total system and $\langle i,j\rangle$ denotes the edge connecting the sites $i$ and $j$. Since spins in the $A$ sublattice do not couple to each other, Eq. (43) can be factorized as

[TABLE]

where $\partial i$ denotes the set of sites adjacent to $i$. This allows us to sample each spin on $A$ sublattice independently without rejection. Update method that alternately switches the sublattice is called the block Gibbs sampling.

However, such an algorithm is not beneficial for the following reason. Consider a domain consisting of upward spins. In the bulk region of the domain, it is highly probable according to the Eq. (45) that the newly sampled spins remain upwards as well. In other words, the limited range of the connection results in the mutual locking structure in the bulk region, allowing only the peripheral region to flip. Such a problematic situation is exacerbated as the domain size grows, and turns out that the slowing down is much worse than the single-spin flip.

Appendix C Partition function of three-spin interacting model on Kagomé lattice

In this Appendix, we show that a model with only three-spin interactions on the Kagomé lattice is soluble. This model is a specific case of a broader class of models with crossing symmetry studied in Ref. Simon and Fendley (2013), which do not exhibit a phase transition at finite temperature. The partition function is written as

[TABLE]

where $M_{\triangle}$ denotes the three-spin interaction and $\tau_{\triangle}:=\prod_{j\in\triangle}\sigma_{j}$ the product of the Ising spins in an upward triangle in the lattice. Also $M_{\bigtriangledown}$ and $\tau_{\bigtriangledown}$ are defined similarly for a downward triangle. In order to compute Eq. (47), we introduce the identity for a binary variable $x=\pm 1$ as follows,

[TABLE]

Applying this identity to each triangle yields

[TABLE]

Next, let us consider taking the sum over $\sigma_{j}$ at some site $j$ in Eq. (47). Denoting the triangles touching the site $j$ as $\triangle(j)$ and $\bigtriangledown(j)$, the contribution from the spin at $j$ can be given as $\sum_{\sigma_{j}}\sigma_{j}^{n_{\triangle(j)}+n_{\bigtriangledown(j)}}$, which is nonzero only when $n_{\triangle(j)}=n_{\bigtriangledown(j)}$. This argument holds for arbitrary $j$, and therefore the requirement $n_{\triangle}=n_{\bigtriangledown}=0,1$ for all $\triangle$ and $\bigtriangledown$ imposed for nonzero contribution. Accordingly, we obtain the concise expression of the partition function as

[TABLE]

where $\displaystyle C=2^{N_{v}}\prod_{\triangle}\prod_{\bigtriangledown}\cosh(M_{\triangle})\cosh(M_{\bigtriangledown})$ with $N_{v}$ being the number of visible spins, i.e., the number of sites. The above expression clearly shows that the free energy in the thermodynamic limit is analytic, and hence the model does not show a phase transition at finite temperature. For $M_{\triangle}=M_{\bigtriangledown}$, Eq. (50) reproduces the partition function of the uniform model studied in Ref. Muttalib et al. (2017). We note that the present method is not limited to two-dimensional models. In fact, a similar model with four-spin interactions on a three-dimensional pyrochlore lattice is also soluble using the same technique Barry and Wu (1989).

Appendix D Cluster updates under magnetic field

In this Appendix, we introduce two flavors of cluster updates accompanied with magnetic fields. We assume a model defined on a graph $G=(V,E)$ as follows:

[TABLE]

where $J_{i,j}$ is the two-spin interaction between two binary degrees of freedom at sites $i$ and $j$, or $\sigma_{i}$ and $\sigma_{j}$, and $h_{i}$ is the external field on site $i$. One way to take the external field into account is to modify the probability of flipping the cluster that is formed using the information on the interactions, and the other is to extend the space to express the field terms by interactions with the auxiliary space.

To introduce the first approach, let us remind that the detailed balance condition is given as,

[TABLE]

where $\pi(A)$ is the Boltzmann weight corresponding to a state $A$, $P(A\rightarrow B)$ is the trial proposal probability from the state $A$ to $B$, and $\alpha(A\rightarrow B)$ is the corresponding acceptance probability. The Swendsen-Wang algorithm applied in the main text adopts the Metropolis-Hastings rule,

[TABLE]

which satisfies $\alpha=1$ under $h_{i}=0$. Since the external field modifies the Boltzmann weight as $\pi(A)\rightarrow e^{-\beta\sum_{i}h_{i}\sigma_{i}^{A}}\cdot\pi(A)$ at the inverse temperature $\beta$, the trial proposal must absorb such change to realize a rejection-free scheme under arbitrary external field. Concretely, the $k$-th cluster $C_{k}$ formed by the ordinary bonding process is flipped with probability $p_{k}=e^{-\beta m_{k}}/(e^{-\beta m_{k}}+e^{+\beta m_{k}})$ where $m_{k}=-\sum_{i\in C_{k}}h_{i}\sigma_{i}$ is the Zeeman energy by the external field. The additional computational effort per single MC step is ignorable.

In the second approach, known as the “ghost spin method,” one introduces an auxiliary spin that interacts with any spin exposed to the external (or virtual) field Coniglio et al. (1989); Kent-Dobias and Sethna (2018). Defining $\widetilde{G}=(\widetilde{V},\widetilde{E})$ with the ghost spin on the [math]-th site as

[TABLE]

we alternatively consider a Hamiltonian as follows,

[TABLE]

where

[TABLE]

Now that the new Hamiltonian consists solely of two-spin interactions, the ordinary cluster update can be applied.

Appendix E Observation of physical quantity in extended model

The transformation considered in the main text preserves the partition function, and moreover the Boltzmann factor for visible spin configurations. Therefore, to compute the expectation value of a physical observable $O({\bm{\sigma}})$ in the extended space, one may simply consider the identical mapping $\widetilde{O}(\sigma,h)=O(\sigma)$ to obtain

[TABLE]

In other words, one may simply ignore all the hidden spins and compute the quantities using the operator in the original space.