Adaptive Cluster Expansion for Ising spin models

TL;DR

This paper introduces an adaptive cluster expansion algorithm to efficiently approximate solutions of the direct Ising problem, enabling computation of free energy and observables for complex spin systems.

Contribution

It adapts the cluster expansion method, originally for inverse problems, to directly solve the Ising model, improving computational efficiency and applicability.

Findings

Effective in one and two-dimensional systems

Handles random and fully connected graphs

Works with homogeneous and heterogeneous parameters

Abstract

We propose an algorithm to obtain numerically approximate solutions of the direct Ising problem, that is, to compute the free energy and the equilibrium observables of spin systems with arbitrary two-spin interactions. To this purpose we use the Adaptive Cluster Expansion method, originally developed to solve the inverse Ising problem, that is, to infer the interactions from the equilibrium correlations. The method consists in iteratively constructing and selecting clusters of spins, computing their contributions to the free energy and discarding clusters whose contribution is lower than a fixed threshold. The properties of the cluster expansion and its performance are studied in detail on one dimensional, two dimensional, random and fully connected graphs with homogeneous or heterogeneous fields and couplings. We discuss the differences between different representations (Boolean and Ising) of the spin variables.

Figures

Equations

$$H_{\mathcal{J}}(\mathbf{\sigma })=-\sum _i{h}_i{\sigma }_i-\sum _{i<j}{J}_{ij}{\sigma }_i{\sigma }_j.$$

$$F_{\mathcal{J}}=-{\beta }^{-1}\log Z(\mathcal{J}),$$

$$\mathcal{C}=\{\{⟨{\sigma }_i⟩\},\{⟨{\sigma }_i{\sigma }_j⟩\}\}.$$

$$Z(\mathcal{J})=\sum _{\mathbf{\sigma }}{\text{rme}}^{-\beta H_{\mathcal{J}}(\mathbf{\sigma })},$$

$$\begin{array}{rl}& \Delta W_i({\mathcal{J}}_i)=W_i({\mathcal{J}}_i),\\ & \Delta W_{\Gamma }({\mathcal{J}}_{\Gamma })=W_{\Gamma }({\mathcal{J}}_{\Gamma })-\sum _{{\Gamma }^{\prime }\subset \Gamma }\Delta W_{{\Gamma }^{\prime }}({\mathcal{J}}_{{\Gamma }^{\prime }}),\end{array}$$

$$W_{\Gamma }({\mathcal{J}}_{\Gamma })=\sum _{{\Gamma }^{\prime }\subseteq \Gamma }\Delta W_{{\Gamma }^{\prime }}({\mathcal{J}}_{{\Gamma }^{\prime }}).$$

$$W_{\Omega }(\mathcal{J})=\sum _{\Gamma \subseteq \Omega }\Delta W_{\Gamma }({\mathcal{J}}_{\Gamma }).$$

$$\left(\genfrac{}{}{0ex}{}{N}{1}\right)+\left(\genfrac{}{}{0ex}{}{N}{2}\right)+\dots +\left(\genfrac{}{}{0ex}{}{N}{N}\right)=\sum _{k=0}^N\left(\genfrac{}{}{0ex}{}{N}{k}\right)-\left(\genfrac{}{}{0ex}{}{N}{0}\right)=2^N-1.$$

$$W_{\Gamma }({\mathcal{J}}_{\Gamma })=\log Z_{\Gamma }({\mathcal{J}}_{\Gamma })=\log \sum _{{\sigma }_i,i\in \Gamma }{e}^{{\sum }_{i\in \Gamma }\beta h_i{\sigma }_i+{\sum }_{\{i,j\}\subseteq \Gamma }\beta J_{ij}{\sigma }_i{\sigma }_j},$$

$$\begin{array}{rl}⟨{\sigma }_i⟩& =\frac{1}{\beta }\frac{d{W}_{\Omega }(\mathcal{J})}{d{h}_i}=\frac{1}{\beta }\sum _{\{i\}\subseteq \Gamma \subseteq \Omega }\frac{d\Delta W_{\Gamma }({\mathcal{J}}_{\Gamma })}{d{h}_i},\\ ⟨{\sigma }_i{\sigma }_j⟩& =\frac{1}{\beta }\frac{d{W}_{\Omega }(\mathcal{J})}{d{J}_{ij}}=\frac{1}{\beta }\sum _{\{i,j\}\subseteq \Gamma \subseteq \Omega }\frac{d\Delta W_{\Gamma }({\mathcal{J}}_{\Gamma })}{d{J}_{ij}}.\end{array}$$

$$W_{\Omega }(\mathcal{J})\sim \sum _{\Gamma \in L}\Delta W_{\Gamma }({\mathcal{J}}_{\Gamma }).$$

$$H=-J\sum _{i=1}^N{s}_i{s}_{i+1},s_{N+1}=s_1,s_i\in \{-1,1\}.$$

$$H=-4J\sum _{i=1}^N{\sigma }_i{\sigma }_{i+1}+4J\sum _{i=1}^N{\sigma }_i-NJ,{\sigma }_{N+1}={\sigma }_1,{\sigma }_i=\{0,1\}.$$

$$\Delta W_n=\log Z_n-2\log Z_{n-1}+\log Z_{n-2}=0,\text{mbox}forn>2,$$

$$e^{\beta J{s}_i{s}_j}=\cosh (\beta J)+s_i{s}_j\sinh (\beta J)=\cosh (\beta J)(1+w{s}_i{s}_j),w=\tanh (\beta J),$$

$$Z=\sum _{\mathbf{s}}{e}^{{\sum }_{(ij)}\beta J{s}_i{s}_j}=\cosh (\beta J{)}^{\ell }\sum _{\mathbf{s}}\prod _{(ij)}(1+w{s}_i{s}_j),$$

$$\sum _{\mathbf{s}}{s}_{i_1}{s}_{i_2}\dots s_{i_n}=\{\begin{array}{ll}{2}^N& \text{mbox}ifallspinsappearanevennumberoftimes,\\ 0& \text{mbox}otherwise.\end{array}$$

$$\begin{array}{rl}Z& =2^N\cosh (\beta J{)}^{\ell }\left[1+\sum _{p>0}{w}^p{g}_N(p)\right],\\ W_{\Omega }& =\log Z=N\log 2+\ell \log \cosh (\beta J)+\log \left[1+\sum _{p>0}{w}^p{g}_N(p)\right],\end{array}$$

$$W_{\Omega }=N\log 2+N\log \cosh (\beta J)+\log (1+w^N).$$

$$H=-J\sum _{i=1}^N{s}_i{s}_{i+1}-h\sum _{i=1}^N{s}_i$$

$$H=-4J\sum _{i=1}^N{\sigma }_i{\sigma }_{i+1}-2(h-2J)\sum _{i=1}^N{\sigma }_i+hN-NJ.$$

$$H=-J\sum _{(ij)}{s}_i{s}_j,s_i=\{-1,+1\}.$$

$$\begin{array}{rl}& k\doteq \frac{1}{\sinh {\left(2\beta J\right)}^2}\\ & W=-N\beta f=N\left(\frac{\log (2)}{2}+\frac{1}{2\pi }{\int }_0^{\pi }\log \left[\cosh {\left(2\beta J\right)}^2+\frac{1}{k}\sqrt{1+k^2-2k\cos (2\theta )}\right]d\theta \right)\\ & m=⟨s_i⟩=[1-{\sinh }^{-4}(2\beta J){]}^{1/8},⟨{\sigma }_i⟩=(1+m)/2.\end{array}$$

$$H=-4J\sum _{(ij)}{\sigma }_i{\sigma }_j+8J\sum _i{\sigma }_i-2NJ,{\sigma }_i=\{0,1\}.$$

$$H=-J\sum _{i=1}^N{s}_i{s}_{i+1}-\sum _i^N{h}_i{s}_i,h_i\sim \mathcal{N}(0,1),$$

$$H=\sum _{i=1}^N(4J-2h_i){\sigma }_i-4J\sum _{i=1}^N{\sigma }_i{\sigma }_{i+1}+NJ-\sum _{i=1}^N{h}_i.$$

$$H=-\sum _{(ij)}{J}_{ij}{s}_i{s}_j-\sum _{i=1}^N{h}_i{s}_i,s_i=\{-1,1\},$$

$$H=\sum _i\left(2\sum _{j\in \partial i}{J}_{ij}-2h_i\right){\sigma }_i-4\sum _{(ij)}{J}_{ij}{\sigma }_i{\sigma }_j-\sum _{(ij)}{J}_{ij}+\sum _i{h}_i,{\sigma }_i=\{0,1\}.$$

$$\frac{\partial W}{\partial \beta }=-⟨H{⟩}_{\beta }\text{mbox}soW(\beta )=W(0)-{\int }_0^{\beta }⟨H{⟩}_{{\beta }^{\prime }}d{\beta }^{\prime },$$

$$H=-\sum _{i<j}\frac{J_{ij}}{\sqrt{N}}{s}_i{s}_j-\sum _i{h}_i{s}_i,s_i=\pm 1.$$

Full text

Keywords Cluster expansion, Ising model

1. 1 Introduction

2. 2 Cluster expansion

3. 2.1 Adaptive Cluster Expansion: general construction and selection rules

4. 2.2 Significance threshold

5. 2.3 Observables

6. 3 Applications to ordered systems

7. 3.1 One-dimensional Ising model without external field

8. 3.2 One-dimensional Ising model with homogeneous external field

9. 3.3 Two-dimensional Ising model

10. 3.4 Summary

11. 4 Applications to disordered systems

12. 4.1 One-dimensional model with random fields

13. 4.2 Two-dimensional Edwards-Anderson spin glass model

14. 4.3 Sherrington-Kirkpatrick model

15. 4.4 Erdös-Rényi random graph model

16. 4.5 Application on biological data

17. 5 Conclusions

1 Introduction

We consider a finite number $N$ of binary variables, either Boolean variables $\sigma_{i}\in\{0,1\}$ or Ising spin variables $s_{i}\in\{-1,1\}$, interacting via a set of $N$ local fields $h_{i}$ and of $\frac{1}{2}N(N-1)$ pairwise couplings $J_{ij}$. We define $\Omega\doteq\{1,\ldots,N\}$ and $\mathcal{J}=\{\{h_{i}\},\{J_{ij}\}\}$. The energy of a configuration $\boldsymbol{\sigma}=\{\sigma_{1},\ldots,\sigma_{N}\}\in\{0,1\}^{N}$ is given by the Hamiltonian:

[TABLE]

In statistical mechanics, one is usually interested in computing the expectation values of certain physical quantities, like the free energy of the system in thermal equilibrium:

[TABLE]

and the set $\mathcal{C}$ of the equilibrium correlations, meaning the average values of the spins and the spin-spin correlations over the Gibbs measure corresponding to the Hamiltonian (1):

[TABLE]

Both calculations require the computation of the partition function

[TABLE]

the required computational time being $\mathcal{O}(2^{N})$. Practically speaking, for a large system the calculation then becomes impossible, and approximation methods are required.

This is the problem we are interested in this paper, and we refer to it as the direct Ising problem. The inverse Ising problem consists, on the contrary, in finding the interaction set $\mathcal{J}$ from the measure of the correlations $\mathcal{C}$ [1, 2, 3, 4]. The Adaptive Cluster Expansion (ACE) method [1, 5] has been originally developed to solve the inverse problem. The aim of this paper is to adapt the ACE method to solve the direct Ising problem: in particular, we want to compute the free energy and the correlations $\mathcal{C}$ of the model from the set $\mathcal{J}$.

Let us stress from the very beginning that, of course, several very efficient methods to solve the direct Ising problem are available in the literature. In the following, and in the bulk of the paper, we review them shortly and discuss what kind of advantages (and disadvantages) the ACE method could offer with respect to those methods. As a general comment, we believe that the ACE method has the advantage of being very general and simple to implement; but on the other hand, it is certainly not competitive for large-scale simulations of ordered systems. Its efficiency for disordered systems should be discussed on a case-to-case basis. If one is interested in obtained a good (but not necessarily excellent) approximation of the free energy and correlations for a large ensemble of rather small systems with generic random couplings and fields, then we believe that the ACE method offers a good alternative to existing methods and could outperform Monte Carlo schemes if disorder is strong enough. Note that the situation described above is quite common in applications to biological problems, see e.g. [6].

Several cluster methods have been used before to solve the direct Ising problem. Indeed, for ordered models defined on a regular lattice, the set of clusters that should be kept in the expansion can be chosen ad hoc by following the symmetries of the problem, as in the Ising model [7]. In fully connected model or large-dimensional models, both ordered and disordered, only a finite number of clusters contribute in the thermodynamic limit [8]; examples are the Sherrington-Kirkpatrick model [9, 10] or even more general spin systems such as the Heisenberg model [11]. However, these methods cannot be easily extended to general disordered systems where underlying symmetries are missing and all clusters are potentially relevant. The ACE algorithm, on the contrary, is designed to automatically find the most relevant clusters. Other expansion schemes such as belief propagation can provide an estimation of the free energy, at least for sparse graphs with few loops [12]. Generalisations of belief propagation [12, 13] and cluster expansions [14, 8, 15, 16, 17] have been used to construct several more sophisticated expansion methods. Compared to these classical cluster expansions, the ACE method is simpler to implement, because it does not require any consistency equation between the marginals calculated from different clusters, or between a cluster and the whole system.

A major problem in cluster expansion methods is how to choose clusters and how to sum them up in an appropriate way to obtain the partition function. Here, we base our cluster expansion on a construction rule which selectively builds up larger clusters from smallest ones, and on a selection rule to sum up clusters according to the amplitudes of their free energy contributions [1, 5]. It is important to point out that the method is not exact: for instance, because we do not impose any consistency equation, the choice of the representations of the variables (e.g. Ising or Boolean) is very important for the rapid convergency of the algorithm, as we will discuss in the following. Nevertheless, we find that the method has good performance in interesting cases.

The Monte Carlo (MC) technique is also very efficient to compute the spin correlations, but the calculation of the free energy via Monte Carlo is less efficient, because it requires a thermodynamic integration [18, 19, 20]. Several biased Monte Carlo sampling schemes have been developed in this context to better evaluate the free energy of the model. An exhaustive review of these methods goes beyond the scope of the paper; we mention here as examples importance sampling [19], the Wang-Landau approach [21], the component-distribution approach [22], cluster methods [19, 23], umbrella sampling [19, 20] and tethered methods [24]. In the following, we will limit ourselves to compare our ACE implementation with the simplest MC scheme, namely a standard MC using thermodynamic integration.

In the following, we first present the ACE method (section 2), and then apply it to some classical examples of Ising models (sections  3 and 4), both on ordered and disordered systems and at different temperatures. We conclude by some comments on the comparison with Monte Carlo performances, and perspectives for future applications of the method (section 5).

2 Cluster expansion

Let us call a cluster $\Gamma$ a subset of $\Omega$ and consider a function $W_{\Gamma}(\mathcal{J}_{\Gamma})$ that depends only on the interaction parameters $\mathcal{J}_{\Gamma}=\{h_{i},J_{ij}\}_{i,j\in{\Gamma}}$ that involve variables in the cluster $\Gamma$. We can derive the contribution of the cluster $\Gamma$ to such function, $\Delta W_{\Gamma}$, from the following recursive formula:

[TABLE]

and therefore:

[TABLE]

For the whole system we thus get:

[TABLE]

Eq. (7) is an exact expansion, called Möbius transform, and the total number of clusters $\Gamma$ that appear in the sum is:

[TABLE]

The expansion holds for an arbitrary function, but here we are interested in

[TABLE]

in such a way that $W_{\Omega}(\mathcal{J})$ is the log partition function. Here, $\beta$ is the inverse temperature and fixes the overall scale of the couplings.

The free energy is simply related to the log partition function by $F_{\Omega}(\mathcal{J})=-\beta^{-1}W_{\Omega}(\mathcal{J})$. From the log partition function, differentiating with respect to $h_{i}$ and $J_{ij}$ respectively, we can compute the correlations in Eq. (3). Because the functions $W_{\Gamma}$ and $\Delta W_{\Gamma}$ depend only on the couplings in the cluster $\Gamma$, we have:

[TABLE]

The expansion in Eq.(7) is useless if one has to compute all the clusters: the total number of clusters is $2^{N}-1$ and computing large clusters is equivalent to a full computation of the partition function. We therefore need an accurate truncation scheme, that allows one to keep only a certain number of small enough clusters and yet provides a good approximation for the free energy and the correlations. First, we note that only connected clusters contribute to the expansion for the log partition function. A cluster $\Gamma$ is disconnected if it can be separated in two disjoint subsets, $\Gamma=\Gamma_{1}\cup\Gamma_{2}$ with $\Gamma_{1}\cap\Gamma_{2}=\emptyset$, in such a way that there is no interaction $J_{ij}$ connecting spins in the two subsets, i.e. $J_{ij}=0$ if $i\in\Gamma_{1}$ and $j\in\Gamma_{2}$. The log partition function satisfies an additive property: if $\Gamma$ is disconnected, then $W_{\Gamma}(\mathcal{J}_{\Gamma})=W_{\Gamma_{1}}(\mathcal{J}_{\Gamma_{1}})+W_{\Gamma_{2}}(\mathcal{J}_{\Gamma_{2}})$. Thanks to this property, it is easy to show that $\Delta W_{\Gamma}=0$ for disconnected clusters.

Second, we note that $\Delta W_{\Gamma}$ becomes small when the size of the cluster grows. In fact, when the size of the cluster becomes larger than the correlation length of the model, the far away spins have very weak interactions and the cluster behaves almost like a disconnected cluster. In this situation, it is reasonable to expect that $\Delta W_{\Gamma}$ vanishes quickly. Hence, we can truncate the series by selecting a list $L$ of small enough connected clusters, thus obtaining an approximation of the function $W_{\Omega}$:

[TABLE]

The quality of the approximation of Eq.(11) depends on the value of the neglected clusters, i.e. on the procedure used to construct the list $L$, which we discuss in the next section.

2.1 Adaptive Cluster Expansion: general construction and selection rules

The main problem of a cluster expansion is how to choose the list $L$ of clusters. One possibility is to sum them by size, i.e. to consider all clusters whose length is lower than a given $K_{max}$. Because the calculation of the partition function grows exponentially with the size of the clusters, $K_{max}$ is fixed by computational limitations. Here, using a normal personal computer, we can reach a $K_{max}=18$ . However, clusters having the same length can have very different $\Delta W_{\Gamma}$ and clusters having different length can have comparable and sometimes opposite in sign contributions, so this is not a well-converging procedure.

A more efficient one is to discard the clusters whose $|\Delta W_{\Gamma}|$ is smaller than a fixed threshold $\Theta$ (we will call them non significant clusters) [1, 5]. Clearly, even for rather small systems, it is not possible to compute all of the $2^{N}-1$ clusters $\Delta W_{\Gamma}$ to find the significant ones. Thus we use a recursive construction rule for the clusters (see the pseudocode of Algorithm 1, readapted from Ref. [5]). We begin with the computation of $\Delta W_{(i)}=W_{(i)}=\log(\sum_{\sigma_{i}=0,1}e^{\beta h_{i}\sigma_{i}})$ for all $N$ clusters of size $k=1$. Then we compute $\Delta W_{(i,j)}$ for all ${N\choose 2}$ clusters of size $k=2$. From Eq. (5), that gives $\Delta W_{(i,j)}=W_{(i,j)}-W_{(i)}-W_{(j)}$. Each subsequent step follows the same pattern. Firstly, clusters with $\mid\Delta W_{\Gamma}\mid<\Theta$ are discarded and we add them to the list of non significant clusters. We include in the next step of the expansion all clusters which are unions of two significant clusters of size $k$, $\Gamma=\Gamma_{1}\cup\Gamma_{2}$, such that the new cluster $\Gamma$ contains $k+1$ spins. Note that this rule prevents any combinatorial explosion of the computational time, as it only selects a relatively small subset of the total number of clusters. At the end of the procedure, we have a list $L$ of significant clusters and we use them in the expansion of Eq. (11) to compute an approximation of the true free energy111Note that we considered also stricter versions of this rule, where a higher number of significant subclusters is required for a cluster to be significant, but we found that such stricter choices prevent the inclusion of relevant clusters, so in the following we stick to this rule..

Algorithm 1, which contains the main iteration loop that selects the clusters, calls as a subroutine Algorithm 2, which calculates the contribution $\Delta W_{\Gamma}$ of the cluster $\Gamma$ and, in turn, calls a subroutine, Algorithm 3, to calculate the free energy of the clusters, $W_{\Gamma}$. Note that the derivatives of the cluster free energies are required to compute the correlations according to Eq. (10) and Algorithm 1. Explicit expressions of these derivatives can be obtained analytically, and their numerical evaluation requires only a minor modification of Algorithm 3. The time to calculate $W_{\Gamma}$, for a cluster of size $k$, is $O(2^{k})$, but the time to select clusters and compute $\Delta W_{\Gamma}$ is approximately linear in the number of clusters.

2.2 Significance threshold

Choosing a good value of $\Theta$ is of crucial importance: if we choose a too small $\Theta$ we will likely consider lots of long clusters and the computational time will increase. On the other hand if we choose a too high $\Theta$ we will consider too few cluster and the estimated value of $W_{\Omega}$ will be far from the true value. Because we do not know a priori the optimal value of $\Theta$, we run the algorithm several times decreasing progressively the threshold:

• •

Start with a large value of $\Theta$.

• •

Once all possible new clusters are constructed calculate $M=\max_{\Gamma\in\overline{L}}|\Delta W_{\Gamma}|$, where $\overline{L}$ is the list of all the non-significant clusters that have been constructed and tested in the previous run.

• •

Set a new value $\Theta=0.99M$ of the significance threshold.

In this way, at each step we are sure that new clusters will be considered significant, and new clusters will be constructed. In this paper, in order to investigate the convergence properties of the algorithm, we stop the program when $M=0$, which indicates that all the connected components of the system have been considered and no more significant clusters can be constructed; or, we set a minimal value of $\Theta_{min}$ and we stop the program when $\Theta<\Theta_{min}$. Obviously, in practical applications, one can stop the code as soon as convergence of the relevant observables is reached.

For ordered systems, as we will see in more detail in the following, some clusters give a zero free energy contribution because of symmetries. As an example, it can be shown (see section 3) that in ordered lattices with no external field and in the Ising spin representation, non zero contributions for clusters containing more than 2 spins come only from closed loops. The construction rule must then be adapted, for such ordered cases, in order to avoid that the expansion stops because the list of clusters of a given size is empty.

2.3 Observables

In the following we report the results of our algorithm when applied to Ising models with different interaction structures. Interesting quantities are:

• •

the approximate value of $W_{\Omega}=\log Z_{\Omega}(\mathcal{J})$ as a function of $\Theta$, given by $W_{\Omega}(\Theta)=\sum_{\Gamma\in L(\Theta)}\Delta W_{\Gamma}$, where the list $L(\Theta)$ includes all the significant clusters,

• •

$N_{con}(k)$, the number of clusters of size $k$ constructed by the algorithm, and the total $N_{con}=\sum_{k}N_{con}(k)$,

• •

$N_{sel}$, the number of significant clusters, i.e. the clusters $\Gamma\in L(\Theta)$,

• •

$K_{max}$, the size of the largest significant clusters $\Gamma\in L(\Theta)$,

• •

$N_{op}=\sum_{k}N_{con}(k)2^{k}$ which is an estimate of the number of operations performed by the algorithm, because constructing a cluster of size $k$ requires $2^{k}$ operations.

We will examine the behavior of these quantities in several examples of Ising models. Along the way, we illustrate several important properties of the cluster expansion. Because most of the examples we consider are diluted, it is convenient to introduce the following notation. The interaction graph ${\cal G}$ is such that a link $(ij)\in{\cal G}$ if and only if $J_{ij}\neq 0$. The size $\ell=|{\cal G}|$ is the number of non-zero couplings. Also, we call $\partial i$ the set of neighbors of $i$, i.e. $\partial i=\{j:J_{ij}\neq 0\}$.

3 Applications to ordered systems

3.1 One-dimensional Ising model without external field

We first consider the one-dimensional Ising model with Ising spins and nearest-neighbour interactions, $J_{i,i+1}=J$, and zero external field, $h_{i}=0$. The spins are on a ring with $N$ sites (periodic boundary conditions). The Hamiltonian is:

[TABLE]

Using Boolean variables $\sigma_{i}=(s_{i}+1)/2$, the Hamiltonian becomes:

[TABLE]

In Figure 1 we report results for the case $N=100$ and $J=1$ for different $\beta$.

The first observation is that, although the final result for $W_{\Omega}$ must be independent of the representation of the binary variables (Ising or Boolean), the intermediate steps of the cluster expansion strongly depend on the representation. The reason is that interactions are classified differently in the two representations. In the Ising spin representation, the one-spin clusters have zero field and coupling, hence $W_{(i)}=\log 2$. Also, in the Ising spin representation, all clusters keep the symmetry under inversion $s_{i}\to-s_{i}$ of the original Hamiltonian, so the contribution to the average magnetization $\langle s_{i}\rangle$ is zero for all clusters. On the contrary, in the Boolean representation the one-spin clusters have a field $h_{i}=-4J$, hence $W_{(i)}=\log\sum_{\sigma=0,1}e^{-4\beta J}=\log(1+e^{-4\beta J})$. Also, the one-spin clusters contribute to the average a quantity $\langle\sigma_{i}\rangle=e^{-4\beta J}/(1+e^{-4\beta J})<0.5$. Despite the fact that the Hamiltonian is globally invariant under $\sigma_{i}\to 1-\sigma_{i}$, the symmetry results from a cancellation between fields $h_{i}$ and couplings $J_{ij}$, and is broken by the cluster expansion. While the cluster expansion in $s_{i}$ is an expansion around independent paramagnetic spins, for large enough $\beta J$, the cluster expansion in $\sigma_{i}$ is an expansion around independent strongly magnetized Boolean spins with $\langle\sigma_{i}\rangle\sim 0$. For this reasons, in a system with inversion symmetry (which is not spontaneously broken), the Ising spin representation might seem more appropriate.

There is, however, another important consequence of the symmetry in the Ising-spin variable. In the one-dimensional case, clusters made by non-consecutive spin have $\Delta W=0$, while due to translational invariance the $n$ spins clusters are all identical. The $\Delta W_{n}$ of a cluster made by $n$ consecutive spins is $\Delta W_{n}=W_{n}-2\Delta W_{n-1}-3\Delta W_{n-2}-\ldots$, which implies $\Delta W_{n}=W_{n}-2W_{n-1}+W_{n-2}$ for $n>2$. We have $Z_{1}=2$, $Z_{2}=2(e^{\beta J}+e^{-\beta J})$, $Z_{3}=2(e^{\beta J}+e^{-\beta J})^{2}=2^{2-1}Z_{2}^{2}$ and in general $Z_{n}=2^{2-n}Z_{2}^{n-1}$ for $n>2$. Then

[TABLE]

which shows that all clusters with $n\geq 3$ have vanishing $\Delta W$, except for the cluster made by all $N$ spins which has an additional contribution, exponentially small in $N$, coming from the $N$ consecutive spins. Summing only one- and two-spin clusters is therefore accurate with an exponentially small error in $N$. Moreover, as we will see in another example below, the fact that all three-spin clusters have $\Delta W=0$ can become a problem, if there are higher order clusters that have non-zero $\Delta W$. Therefore, we conclude that all clusters with $n>2$ spins have $\Delta W=0$, and the program will stop after constructing the two-spin clusters. Finally, it is important to stress that in order to obtain $k$-spin connected correlations we need clusters of size $k$ even if the free energy has already reached convergence; in fact, these clusters give zero contribution to the free energy, but a finite contribution to spin correlations. Our selection rule on the free energy is then not appropriate to compute correlations in this case.

As shown in Figure 2, one possibility to solve this problem is to add a small magnetic field to the Hamiltonian. In this way all clusters with $n>2$ acquire a $\Delta W$ which, while being exponentially small, is not zero anymore. Therefore the program will not stop after constructing the two-spin clusters, but can reach clusters of arbitrary dimension $k$ and obtain the $k$-spin connected correlations.

Note that the above result is a simple instance of a more general result for binary variables on ordered interaction graphs. Using the identity

[TABLE]

the partition function can be rewritten in the following way (which is essentially a high-temperature expansion):

[TABLE]

where $\ell$ is the total number of links $(ij)$ carrying a non-zero coupling $J$. Furthermore, the sum of a product of spins is

[TABLE]

Developing the product $\prod_{(ij)}(1+ws_{i}s_{j})$ one obtains products of spins over links. The only way to have each spin appearing twice is to form a closed polygon. Hence

[TABLE]

where $g_{N}(p)$ the number of closed polygons in $\cal G$ having perimeter $p$. Specializing this to a one-dimensional ring of $N$ sites, we have $N$ links and the only closed polygon is the ring itself, with perimeter $p=N$, so that

[TABLE]

The $N$ one-spin clusters contribute $W_{1}=\log 2$, which is the first term. The $N$ two-spin clusters made of neighboring spins contribute $\Delta W_{2}=\log\sum_{s_{1}s_{2}}e^{\beta Js_{1}s_{2}}-2W_{1}=\log[2(e^{\beta J}+e^{-\beta J})]-2\log 2=\log\cosh(\beta J)$, which is the second term. The ring contributes the last term, $W_{N}=\log(1+w^{N})$, which vanishes exponentially for $N\to\infty$.

In summary, as shown in Figure 1, for a one dimensional system with inversion symmetry:

• •

The Ising spin representation converges extremely fast: only one-spin and two-spins clusters are computed, after which the expansion stops. The error in the free energy is due to the $N$-spin clusters, which gives a contribution exponentially small in $N$. The magnetization is zero at all orders, but only nearest-neighbor correlations can be computed.

• •

The Boolean spin representation is slower, because it breaks the inversion symmetry, and it is an expansion around a magnetized state. One has to include several clusters for the free energy and the magnetizations to converge to their asymptotics value. However, in this case one can compute all the $n$-point correlations, provided a sufficient number of clusters is included.

3.2 One-dimensional Ising model with homogeneous external field

The inversion symmetry (and the consequent cancellation of cluster contributions) of the previous model is broken if we add an external field to the system. We then study the case of a homogeneous external field,

[TABLE]

which in the Boolean representation becomes,

[TABLE]

Results for this system are reported in Figure 3 for $\beta=1$, $N=100$ and two values of $h$. The cluster contribution $\Delta W_{\Gamma}$ decreases exponentially with the length of the clusters and the correlation length of the system: $\Delta W_{\Gamma}\sim\exp[-L_{\Gamma}/\xi]$ where $\xi\approx 0.49$ for $h=1$ and $\xi\approx 0.96$ for $h=0.5$.

3.3 Two-dimensional Ising model

We now study the ferromagnetic 2d Ising model on a square lattice with open boundary conditions. The couplings are $J_{ij}=J$ if the distance between the spin $i$ and the spin $j$ is one and zero otherwise, and the fields are absent, $h_{i}=0$, so the model has inversion symmetry. The Hamiltonian is:

[TABLE]

This model has been solved by Onsager in the limit $N\rightarrow\infty$ [25]. The free energy and magnetization are given by:

[TABLE]

In Figure 4 we plot the cluster contributions $\Delta W_{\Gamma}$ as a function of the cluster size. Because the model has inversion symmetry, the polygon expansion in Eq. (18) applies and indeed we observe that only clusters formed by closed polygons have a non zero $\Delta W$. The figure illustrates that our algorithm for Ising spins does not work well in this case. Indeed, new clusters are created by adding a new spin to significant clusters (i.e. with $|\Delta W|>0$), but because all clusters with 3 spins have $\Delta W=0$ (so they are never significant) we cannot create larger clusters using our construction rule. Consequently the expansion is interrupted after clusters of size one and two are constructed, which is not enough to obtain a good value for the free energy. Clearly, this is a pathology due to the particular structure of this problem. The solution is to take into account the clusters having the largest contributions $|\Delta W|$, see Figure 4, with their multiplicity, depending on the system size and of the symmetries. In this way, one can obtain very accurate estimations of the free energy even in the thermodynamic limit, as illustrated in Figure 5.

As in the one-dimensional case, another way to solve this problem can be to study the same system but using Boolean spins, for which the Hamiltonian in Eq. (22) becomes

[TABLE]

The inversion symmetry is then broken by the cluster expansion, and there are significant clusters of any size. In Figure 6 we show the results for a system with $N=10\times 10$ spins on a square lattice, with periodic boundary conditions, $J=1$ and $\beta=0.2,0.5,0.7$ (the Onsager critical point is at $\beta_{c}\sim 0.44$), and we compare the free energy and the magnetization with the Onsager solution, Eq. (23). Note that the expansion is around the state with all $\sigma\sim 0$, because for the Hamiltonian in Eq. (24), the one-spin clusters have a strong field $h=-8J$ that favors the $\sigma=0$ configuration. Therefore, the expansion works better in the ferromagnetic phase, where the free energy and the magnetization converge immediately to one of the two possible values, because the system is very strongly magnetized. On the contrary, in the paramagnetic state where one has to include more clusters to reach zero magnetization. As shown in Figure 7, the algorithm also properly converges when approaching the phase transition temperature $\beta_{c}\approx 0.44$.

Note that Osanger’s solution is only correct in the thermodynamic limit ($N\rightarrow\infty$), where the inversion symmetry can be spontaneously broken, while at finite $N$ one must always find $\langle\sigma\rangle=0.5$ for any $\beta$. However, zero magnetization is recovered only by considering very large cluster (of size comparable with the system size). To show this, in Figure 8 we plot the magnetization $\langle\sigma\rangle$ as a function of $\Theta$ for a small system made by $4\times 4$ spins; it is observed that the magnetization converges to $\langle\sigma\rangle<0.5$, but when the cluster with 16 spins is included the symmetry is restored and it jumps to $\langle\sigma\rangle=0.5$.

3.4 Summary

The main message we learn from the application of the ACE method to ordered systems is that, when the correlation length is finite, the cluster contribution ${\Delta W_{\Gamma}\sim\exp[-L_{\Gamma}/\xi]}$ decays exponentially with the cluster size. Hence, the free energy converges exponentially fast (in the size of the largest clusters) to the thermodynamic limit. However, some symmetries can impose the vanishing of some $\Delta W_{\Gamma}$, which can spoil the efficiency of the construction rule; symmetries should then be carefully taken into account. Finally, both the speed of convergence and the role of symmetries depend on the spin representation (Ising or Boolean).

4 Applications to disordered systems

We now discuss the performances of the ACE when applied to disordered systems. Note that while in many spin glass studies, external fields are absent and the inversion symmetry is thus preserved, we prefer to consider the most general case in which external fields are present. In this way, the inversion symmetry is explicitly broken and the problems discussed in section 3.3 are not present. We also note that in applications to biological systems, external fields are always present.

4.1 One-dimensional model with random fields

We now test the algorithm on a system with the same geometry, $N=100$ spins, ferromagnetic nearest-neighbour interactions, but random external fields. In the Ising spin representation, the Hamiltonian is:

[TABLE]

where ${\cal N}(0,1)$ denotes a Gaussian random variable with zero mean and unit variance. With Boolean spins, Eq. (25) becomes:

[TABLE]

The exact solution of this model can be obtained by the transfer matrix method.

In this case, both with Ising and Boolean spins, the algorithm constructs clusters with arbitrarily large number of spins upon lowering $\Theta$ (Figure 9). As in the case without random fields, we observe that the free energy expansion converges faster with Ising spins than with Boolean spins. However, when looking at the magnetizations and correlations, the convergence of the Ising expansion is only slightly better than the Boolean expansion. Both expansions achieve very good precision in reproducing the magnetizations and correlations when $\Theta\sim 10^{-5}$.

4.2 Two-dimensional Edwards-Anderson spin glass model

We also applied our algorithm to a disordered Ising model on a two-dimensional square lattice, with periodic boundary conditions, the so-called Edwards-Anderson spin-glass model with random external fields [26, 27]. The Hamiltonian is

[TABLE]

where the couplings $J_{ij}$ and the fields $h_{i}$ are i.i.d. random variables chosen from a normal distribution $\mathcal{N}(0,1)$. We also considered the same model with a Boolean spin representation:

[TABLE]

Note that in this case no exact solution is available. Thus we compared our expansion with the result of Monte Carlo simulations, which provide an alternative way to access to the magnetizations and correlations of the system. With Monte Carlo, one can also compute the free energy $W$ by thermodynamic integration, using the relation:

[TABLE]

with $W(0)=N\log(2)$ (where $N$ is the number of spins of the system). We computed the average energy $\langle H\rangle$ with a Monte Carlo-Metropolis algorithm [19] and performed numerically the integration using the Simpson’s rule. The computational time is of order $\mathcal{O}(MC\cdot INT)$ with $MC$ being the number of Monte Carlo steps (including the thermalization phase) and $INT$ the number of intervals needed to perform the numerical integration.

In Figure 10 we plot the results for a system with $10\times 10$ spins. We observe that for the free energy both the Ising and Boolean expansions converge reasonably well, but the Ising expansion is faster. On the contrary we note that in the case $\beta=0.7$, even if the free energy seems to converge, this is not the case for the connected correlations.

4.3 Sherrington-Kirkpatrick model

In this section, we report the results of our algorithm when applied to disordered Ising models on a complete graph. We consider the case with $N=100$ spins, i.i.d. random external fields $h_{i}\sim{\cal N}(0,1)$ and scaled couplings $J_{ij}\sim{\cal N}(0,1)$ (Sherrington-Kirkpatrick or SK model) [27]:

[TABLE]

An exact solution of this model is known in the thermodynamic limit $N\to\infty$ [27]. At the “replica symmetric” level, the solution is

[TABLE]

Here, ${\cal D}z=dz\frac{e^{-z^{2}/2}}{\sqrt{2\pi}}$ is a Gaussian integration measure.

In this case, the algorithm seems to converge towards the correct result, although with large oscillations, see Figure 11. Large oscilllations are due to the fact that the model is fully connected so at each cluster size there is a large number of relevant clusters that need to be evaluated. Therefore, although $K_{max}$ remains relatively small ($K_{max}\leq 6$), the number of operations becomes quickly very large because of the large number of clusters, preventing an efficient evaluation of the free energy. We conclude that the algorithm is not very efficient for densely connected interaction graphs. Note that we did not explore the low-temperature spin glass phase below the de Almeida-Thouless line [27].

4.4 Erdös-Rényi random graph model

In this section we report the results of our algorithm when applied to disordered Ising models on diluted random graphs with finite connectivity. The random networks are generated from the Erdös-Rényi ensemble, where $M=\frac{c}{2}N$ edges are drawn, uniformly and at random, between $N$ points with $c$ the average degree of a vertex. On the selected bonds $(ij)$ the couplings $J_{ij}$ are chosen from $\mathcal{N}(0,1)$. All other couplings and the fields $h_{i}$ are set to zero (Viana-Bray model) [28, 29]:

[TABLE]

In the paramagnetic phase, the solution can be obtained via the cavity method [29] and reads

[TABLE]

This solution holds if $c\,\overline{\tanh^{2}(\beta J_{ij})}\leq 1$. For example, for $c=5$, the paramagnetic phase is stable for $\beta<\beta_{c}=0.5531\cdots$. In the spin glass phase, no exact solution is available but $W$ can be estimated via Monte Carlo for small enough $N$.

When applying the ACE algorithm to this model we noted that the cluster expansion is truncated at small $K_{max}$. The reason is that on this particular graph, there are few closed polygons which are the only ones with $\Delta W_{\Gamma}\neq 0$ according to Eq. (18). All the others clusters have $\Delta W_{\Gamma}=0$ and are thus discarded, similarly to the one-dimensional model. To avoid this problem we set a a minimal cluster length $K_{min}=4$, meaning that all clusters with size less than or equal to 4 are considered significant and used in the ACE algorithm. This procedure allows to consider longer cluster thus improving the convergence of the expansion (see Figure 12). Note that setting $K_{min}=4$ we can construct clusters of length 5, thus the minimal number of operation is $\mathcal{O}{{50}\choose{5}}\sim 10^{7}$ and the size of the largest significant clusters $K_{max}$ is at least 5 for each $\Theta$.

As in the previous cases, we also considered the case where random fields $h_{i}\sim\mathcal{N}(0,1)$ are added to the Hamiltonian. In Figure 13 we report results for this model with $N=100$, temperature $T=2$, and several connectivities. In this case the system is paramagnetic, but due to the presence of external fields the cavity solution is more complicated [29]. We thus compare the ACE with a Monte Carlo simulation, finding good agreement.

4.5 Application on biological data

We also applied our algorithm on data coming from the prefrontex cortex. A population of $N=37$ neurons is recorded and their interactions are represented by a fully connected Boolean spin Hamiltonian, Eq. (28), inferred using the inverse ACE algorithm [30]. Note that the inferred fields are of the order of $h\sim-4$, while the inferred couplings are of the order of $J\sim 0.1$. Here we apply our algorithm on the inferred model. We note that in less than $10^{8}$ operations, the free energy converges and the one- and two-points (non-connected) frequencies $m_{i}=\langle\sigma_{i}\rangle$, $f_{ij}=\langle\sigma_{i}\sigma_{j}\rangle$ are well correlated with those coming from the data. On the contrary, a Monte Carlo sampling with the same number of operations is less correlated with the statistics from the data (see Figure 14). Hence, in this case, the ACE algorithm appears to be faster and more reliable than the Monte Carlo.

5 Conclusions

In this paper we have presented an Adaptive Cluster Expansion (ACE) algorithm for the direct Ising problem, and we discussed its application to some benchmark systems, both ordered and disordered, and with different topologies of the interaction network. Although each case has its own specificities, we can draw some general conclusions. In ordered cases, the simple recursive construction rule of the ACE algorithm can get stuck because of the vanishing of cluster contributions due to symmetries, but it is easy to adapt the construction rule to take symmetries into account. Indeed, one can classify the relevant clusters and find numerically their contributions, which can be summed with the correct multiplicity to obtain an accurate estimation of the free energy. An example has been given in section 3, for the two-dimensional ferromagnetic Ising model, in which only one-spin, two-spin and closed-loop clusters contribute. Another example is given by the Erdös-Rényi model without external field, in which all three-spin cluster contributions vanish and we had to set a $K_{min}=4$ in the construction rule. The simple, recursive, construction rule is particularly adapted to treat cases with disorder where there are no special symmetries that make some clusters vanishing (such as the inversion symmetry). Also, diluted systems are more efficiently treated by the algorithm. On the contrary, the ACE scheme runs into difficulties for densely connected systems.

Because of the locality of the cluster expansion, the choice of variable representation, e.g. Boolean versus Ising, turns out to be important for the convergence of the expansion. As we have shown in the one-dimensional Ising model, when the representation is changed from Ising to Boolean variables, even in absence of external fields in the original Ising model, non zero fields arise from the change to Boolean variables. As a consequence, while the expansion in Ising variables is around independent paramagnetic spins, the expansion for Boolean variables is around strongly magnetized spins. While Ising variables seem more appropriate for the free energy expansion in the particular case of zero field and total inversion symmetry, this is not the case for the correlation functions which are set to zero for spins at distance larger than two in the Ising representation, where clusters larger than two have zero free energy contribution. The quality of the approximation for the free energy at small cluster size depends on the representation, while the cluster contributions decay is controlled in both cases by the correlation length of the model. To calculate correlations between two sites we need clusters which include these sites, which implies that depending on where the expansion stops, some correlations are set to zero, and an error is made which depends on the correlation length of the system. In these cases, one can add a small magnetic field and take a numerical derivative to obtain the correlations. In summary, the Boolean representation is more adapted for systems with strong magnetic fields, such as the neuronal system studied in Figure 14, while the Ising representation is more adapted for systems with weak magnetic fields. Note that the Ising and Boolean representations are only two special cases in a continuum of possible representations where the spin variables take two generic values $a$ and $b$. Given a model, one could, at least in theory, identify an “optimal” representation, meaning two values $a$ and $b$ such that the ACE algorithm converges most quickly. The performance of the algorithm depends in a non-trivial way on $a$ and $b$, because to determine the optimal representation one should a-priori compute the Hamiltonian and cluster symmetries. Here we did not explore the problem further, but see Ref. [31] for a discussion of this point in the context of the inverse problem.

In the case of disordered systems, for which no exact solution is available, and cluster contributions are heterogeneous, we compared our ACE algorithm to the standard Monte Carlo method. While Monte Carlo simulation is a well-established and powerful tool in statistical mechanical calculations, it is well known that the calculation of free energies is a difficult problem for this algorithm, as it requires a thermodynamic integration, e.g. as in Eq. (29). This method requires the computation of the mean energy of a system at different temperatures and the computational time is of order $\mathcal{O}(MC\cdot INT)$ where $MC$ is the number of Monte Carlo steps (including the thermalization phase) and $INT$ the number of intervals needed to perform the numerical integration correctly. These two parameters, which strongly influence the quality of the simulation, depend on the system size and on the system Hamiltonian, and optimizing their value is not straightforward because of the difficulty of checking the equilibration of the Monte Carlo. On the contrary the ACE algorithm depends on a single parameter $\Theta$, the cluster expansion threshold, which can be lowered until convergence, if the number of operations does not exceed the available computational time. Convergence of the observables can thus be monitored directly. We found that its performances are comparable, or in some cases better, than Monte Carlo. We therefore believe that the ACE could be a useful tool for the free energy calculation (possibly adapting the cluster summation rule to the symmetries of the system), in cases where the network is not too dense, no symmetries are present, and disorder is strong enough to slow down the MC equilibration considerably. Note that in this paper we did not test more advanced MC schemes, that might be more efficient than the ACE in some cases, especially on large systems. We believe that no general conclusion can be drawn, and one should test in each case which is the best performing algorithm.

Note that the ACE algorithm can be easily adapted for spin variables taking more than two values. Indeed, the expansion in Eq. (7) is independent of the Hamiltonian of the system and is valid for any arbitrary function $W$. However, to this end, one should modify the routines Algorithm 2 and 3, and this extension can be computationally expensive. Indeed, the ACE algorithm performs an exact computation of $W$ for each cluster; it requires a computational time $\mathcal{O}(2^{k})$ (with $k$ the size of the cluster) which becomes $\mathcal{O}(q^{k})$ when considering spin variables taking $q$ values. In our numerical experiments (with $q=2$) we checked that we can construct maximum clusters containing 20 spins. Therefore, for a $q$-state spin variable, the maximum size of the cluster is $K_{max}={20\log 2}/{\log q}$. As a reference, we can consider the case of $q=21$ which is particularly relevant since it is the model used for the prediction of structural contacts from protein sequences (the inverse ACE algorithm was originally developed for this purpose [1]). In this case $K_{max}\sim 5$, and, in some situations, the direct ACE algorithm could fail to reach convergence with such small clusters. Yet, “compression” tricks can be implemented to reduce the effective number of states [31]. In principle, the ACE algorithm could also be adapted to more general spin systems such as the Heisenberg model [11, 25], classical XY model [25] or Blume-Capel model [32]. Nevertheless, the computation of the cluster contribution $\Delta W$ with continuous spin variables can become more problematic and we did not explore further the possibility.

The ACE algorithm could be improved in several ways. One could for example perform the expansion in Eq.(7) in $W_{\Omega}(\mathcal{J})-W_{0}$, where $W_{0}$ is a reference value obtained, for example, via mean field [5]. This procedure might be helpful for fully connected systems, for which we have seen that the adaptive cluster expansion does not converge well.

It is important to notice that the cluster expansion can be applied to any function, not only the log partition function of an Ising model with pairwise couplings on different interaction graph, as discussed in this paper. An example is the cross-entropy, as originally implemented for the inverse problem [1]. For instance, the cluster expansion can also be used to study the stationary states of a dynamical systems, by considering spatio-temporal evolutions of configurations. Dynamical cluster expansions have been recently built to study for example the spreading of epidemic diseases [33, 34, 35]. The adaptive cluster expansion can also be adapted to include temporal correlations [36].

Acknowledgements

The code we implemented for the direct Ising problem is available online at https://github.com/GiancarloCroce/directACE. It is based on code used for the inverse problem in Ref. [5], https://github.com/johnbarton/ACE. We thank John Barton and Remi Monasson for useful discussions.

References