The Ising model

The Ising model is a statistical physics model of magnetism, also known as the pairwise maxent model [13]. It consists of a set of spins {s_i} with 2 possible orientations (up and down), each responds to its own external magnetic field h_i and each pair is coupled to each other with pairwise coupling J_ij. The strength of the magnetic field determines the tendency of each of the spins to orient in a particular direction and the couplings determine whether the spins tend to point together (J_ij > 0) or against each other (J_ij < 0). Typically, neighbors are defined as spins that interact with one another given by some underlying network structure. Figure 1 shows a fully-connected example.

The energy of each configuration determines its probability via Eq (5),

such that lower energy states are more probable.

We can derive the Ising model from the perspective of maxent. Fixing the the means and pairwise correlations to those observed in the data

we go through the procedure of constructing the Langrangian from Eq 4

where the –1 in Eq 15 has been absorbed into the normalization factor

such that the probability distribution is normalized ∑_sp(s) = 1. Thus, the resulting model is exactly the Ising model mentioned earlier.

Despite the simplicity of the Ising model, the structure imposed by the discrete nature of the spins means that finding the parameters is challenging analytically and computationally. In the last few years, numerous techniques have been suggested for solving the inverse Ising problem exactly or approximately [14]. We have implemented some of them in ConIII and designed a package structure to make it easily extensible to include more methods. Here, we briefly describe the algorithms that are part of the first official version of the package. The goal is to give the user a sense of how they work without getting bogged down in heavy detail. For more detail, we suggest perusing the papers referenced in each section or the review [14]. For a complete beginner, it may be useful to first get familiar with a slower introduction like in the Appendices of Ref [6], Ref [15], or Ref [10].

Enumeration

The naїve approach that only works for small systems is to write out the equations from Eq 9 and solve them numerically. After writing out all K equations,

we can use any standard root-finding algorithm to find the parameters λ_k. This approach, however, involves enumerating all states of the system, whose number grows exponentially with system size.

For the Ising model, writing down the equations has a number of steps $\mathcal{O}\left(K^2{2}^N\right)$M21 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}\left( {{K^2}{2^N}} \right) \] \end{document} , where K is the number of constraints and N the number of spins. Each evaluation of the objective in the root-finding algorithm will be of the same order. For relatively small systems, around N ≤ 15, this approach is feasible on a typical desktop computer and is a good way to test the results of a more complicated algorithm. This approach is implemented by the Enumerate class.

Monte Carlo Histogram (MCH)

Perhaps the most straightforward and most expensive computational approach is Monte Carlo Markov Chain (MCMC) sampling. A series of states sampled from a proposed p(s) is produced by MCMC to approximate ⟨f_k⟩ and determine how close we are to matching ⟨f_k⟩_data. The parameters are then adjusted using a learning rule, and both sampling and learning are repeated until a stopping criterion is met. This can be combined with a variety of approximate gradient descent methods to reduce the number of sampling steps by predicting how the distribution will change if we modify the parameters slightly. The particular technique implemented in ConIII is the Monte Carlo Histogram (MCH) method [16].

Since the sampling step is expensive, the idea behind MCH is to reuse a sample for more than one gradient descent step [16]. Given that we have a sample with probability distribution p(s) generated with parameters λ_k, we would like to estimate the proposed distribution p′(s) from adjusting our parameters λ′_k = λ_k + Δλ_k. We can leverage our current sample to make this extrapolation.

To estimate the average,

To be explicit about the fact that we only have a sampled approximation to p, we replace p with the sample distribution.

Likewise, the ratio of the partition function can be estimated

At each step, we update the Lagrangian multipliers {λ_k} while being careful to stay within the bounds of a reasonable extrapolation. One suggestion is to update the parameters with some inertia [17]

This has a fixed point at the correct parameters.

In practice, MCH can be difficult to tune properly and one must check in on the progress of the algorithm often. One issue is choosing how to set the learning rule parameters η and ∈. One suggestion for η is to shrink it as the inverse of the number of iterations [17]. Another issue is that parameters cannot be changed by too much when using the MCH approximation step or the extrapolation to λ′_k will be inaccurate and the algorithm will fail to converge. In ConIII, this can be controlled by setting a bound on the maximum possible change in each parameter Δλ_max and restricting the norm of the vector of change in parameters ${\sum }_{\text{k}}\sqrt{\Delta {\lambda }_{\text{k}}^2}$M29 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \sum\nolimits_{\rm{k}} {\sqrt {\Delta \lambda _{\rm{k}}^2} } \] \end{document} . Another issue is setting the parameters of the MCMC sampling routine. Both the burn time (the number of iterations before starting to sample) and sampling iterations (number of iterations between samples) must be large enough that we are sampling from the equilibrium distribution. Typically, these are found by measuring how long the energy or individual parameter values remain correlated as MCMC progresses. The parameters may need to be updated during the course of MCH because the sampling parameters may need to change with the estimated parameters of the model. For some regimes of parameter space, samples are correlated over long times and alternative sampling methods like Wolff or Swendsen-Wang would vastly reduce time to reach the equilibrium distribution although these are not included in the current release of ConIII. We do not discuss these sampling details here, but see Refs [18, 19] for examples.

The main computational cost for MCH lies in the sampling step. For each iteration of MCH, the runtime is proportional to the number of samples n, number of MCMC iterations T, and the number of constraints for the Ising model N², $\mathcal{O}\left(\mathrm{Tn}{N}^2\right)$M30 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}\left( {Tn{N^2}} \right) \] \end{document} , whereas the MCH estimate is relatively quick $\mathcal{O}\left(\mathrm{tn}{N}^2\right)$M31 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}\left( {tn{N^2}} \right) \] \end{document} because the number of MCH approximation steps needed to converge is much smaller than the number of MCMC sampling iterations t << T.

MCH is implemented in the MCH class.

Pseudolikelihood

The pseudolikelihood approach is an analytic approximation to the likelihood that drastically reduces the computational complexity of the problem and is exact as N → ∞ [20]. We calculate the conditional probability of each spin s_i given the rest of the system {s_j≠i}

Taking the logarithm, we define the approximate log-likelihood by summing over data points indexed by r:

In the limit where the ensemble is well sampled, the average over the data can be replaced by an average over the ensemble:

To find the point of maximum likelihood for a single spin s_i, we calculate the analytical gradient and Hessian, ∂f/∂J_ij and ∂²f/∂J_ij∂J_i′j′ for a Newton conjugate-gradient descent method. After maximizing likelihood for all spins, the maximum likelihood parameters may not satisfy the symmetry J_ij = J_ji. We impose the symmetry by insisting that

Pseudolikelihood is extremely fast and often surprisingly accurate. Each calculation of the gradient is order $\mathcal{O}\left(R{N}^2\right)$M36 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}\left( {R{N^2}} \right) \] \end{document} and Hessian $\mathcal{O}\left(R{N}^3\right)$M37 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}\left( {R{N^3}} \right) \] \end{document} , which must be done for all N. With analytic forms for the gradient and Hessian, the conjugate-gradient descent method tends to converge quickly.

Pseudolikelihood for the Ising model is implemented in Pseudo.

Minimum Probability Flow (MPF)

Minimum probability flow involves analytically approximating how the probability distribution changes as we modify the configurations [21, 22]. In the methods so far mentioned, the approach has been to maximize the objective (the likelihood function) by immediately taking the derivative with respect to the parameters. With MPF, we first posit a set of dynamics that will lead the data distribution to equilibrate to that of the model. When these distributions are equivalent, then there is no “probability flow” between them. This technique is closely related to score matching, where we instead have a continuous state space and can directly take the derivative with respect to the states without specifying dynamics [23].

First note that Monte Carlo dynamics (satisfying ergodicity and detailed balance) would lead to equilibration to the stationary distribution. The dynamics are specified by a transition matrix, an example of which is given in Ref [22]:

with transition probabilities Γ_ss′ from state s′ to state s. The connectivity matrix ɡ_ss′ specifies whether there is edge between states s and s′ such that probability can flow between them. By choosing a sparse ɡ_ss′ while not breaking ergodicity, we can drastically reduce the computational cost of computing this matrix.

Imagine that we start with the distribution over the states as given by the data and run the Monte Carlo dynamics. When data and model distributions are different, probability will flow between them and indicate that the parameters must be changed. By minimizing a derivative of the Kullback-Leibler divergence, we measure how the difference between the model and the states in the data $\mathcal{D}$M40 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[\cal D\] \end{document} changes when the dynamics are run for an infinitesimal amount of time.

The idea is that this derivative is also minimized with optimal parameters: the MPF algorithm looks for a minimum of the objective function L.

For the Ising model, each evaluation of the objective function where Γ_ss′ connects each data state with G neighbors has runtime $\mathcal{O}(\mathrm{RG}{N}^2)$M42 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}(RG{N^2}) \] \end{document} . In a large fully connected system, G ∼ 2^N would be prohibitively large so a sparse choice is necessary.

MPF is implemented in the MPF class.

Regularized mean-field method

One attractively simple and efficient approach uses a regularized version of mean-field theory. In the inverse Ising problem, mean-field theory is equivalent to treating each binary individual as instead having a continuously varying state (corresponding to its mean value). The inverse problem then turns into simply inverting the correlation matrix C [24]:

where

and where p_i corresponds to the frequency of individual i being in the active (+1) state and p_ij is the frequency of the pair i and j being simultanously in the active state.

A simple regularization scheme in this case is to discourage large values in the interaction matrix J_ij. This corresponds to putting more weight on solutions that are closer to the case with no interactions (independent individuals). A particularly convenient form adds the following term, quadratic in J_ij, to the negative log-likelihood:

In this case, the regularized version of the mean-field solution in (33) can be solved analytically, with the slowest computational step coming from the inversion of the correlation matrix. For details, see Refs. [3, 25].

The idea is then to vary the regularization strength γ to move between the non-interacting case (γ → ∞) and the naively calculated mean-field solution (33) (γ → 0). While there is no guarantee that varying this one parameter will produce solutions that are good enough to “fit within error bars,” this approach has been successful in at least one case of fitting social interactions [3].

The inversion of the correlation matrix is relatively fast, bounded by $\mathcal{O}(N^3)$M46 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}({N^3}) \] \end{document} . Finding the optimal γ involves Monte Carlo sampling from the model distribution, which has computational cost similar to MCH. It is, however, much more efficient because we are only optimizing a single parameter.

This is implemented in RegularizedMeanField.

Cluster expansion

Adaptive cluster expansion [24, 25] iteratively calculates terms in the cluster expansion of the entropy S:

where the sum is over clusters Γ and in the exact case includes all 2^N – 1 possible nonempty subsets of individuals in the system. In the simplest version of the expansion, one expands around S₀ = 0. In some cases it can be more advantageous to expand around the independent individual solution or one of the mean-field solutions described in the previous section [25].

The inverse Ising problem is solved independently on each of the clusters, which can be done exactly when the clusters are small. These results are used to construct a full interaction matrix J_ij. The expansion starts with small clusters and expands to use larger clusters, neglecting any clusters whose contribution ΔS_Γ to the entropy falls below a threshold. To find the best solution that does not overfit, the threshold is initially set at a large value and then lowered, gradually including more clusters in the expansion, until samples from the resulting J_ij fit the desired statistics of the data sufficiently well.

The runtime will depend on the size of clusters included in the expansion. If the expansion is truncated at clusters of size n, the worst-case runtime would be $O\left(\left(\begin{array}{c}N\\ n\end{array}\right)2^n\right)$M48 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\cal O}\left( {\left( {\begin{array}{*{20}{c}} N\\ n \end{array}} \right){2^n}} \right) \] \end{document} . The point is that S can often be accurately estimated even when n ≪ N. The adaptive cluster expansion method is implemented in the ClusterExpansion class.