Variational Pseudolikelihood for Regularized Ising Inference

I propose a variational approach to maximum pseudolikelihood inference of the Ising model. The variational algorithm is more computationally efficient, and does a better job predicting out-of-sample correlations than $L_2$ regularized maximum pseudolikelihood inference as well as mean field and isolated spin pair approximations with pseudocount regularization. The key to the approach is a variational energy that regularizes the inference problem by shrinking the couplings towards zero, while still allowing some large couplings to explain strong correlations. The utility of the variational pseudolikelihood approach is illustrated by training an Ising model to represent the letters A-J using samples of letters from different computer fonts.

Statistical mechanical models constructed from experimental observations provide a valuable tool for studying complex systems. The utility of the statistical mechanics approach to data-driven modeling is especially apparent in biology, providing insights into the behavior of flocking birds [1], the organization of neural networks in the brain [2][3][4], the structure and evolution of proteins [5][6][7][8][9], and many other topics [10,11]. Generally, the 'inverse' statistical mechanics approach refers to the construction of statistical models using the principle of maximum entropy [12,13]. In this approach, one constructs the probability distribution with the maximum entropy subject to constraints on its moments, which are derived from observations. Although maximum entropy modeling has been applied successfully to many different problems, it is still not clear how to estimate the parameters of the resulting probability distribution in an optimal way.
In this work, I consider the problem of estimating the parameters of an Ising model (see [14] for a review). The Ising model describes the statistics of a vector σ of N spin variables σ i ∈ {−1, +1}, and can be derived from the principle of maximum entropy with constraints on the moments σ i and σ i σ j . The probability distribution for σ is given by P ( σ) = Z −1 exp(−U ( σ)) where the energy is U ( σ) = − i h i σ i − i<j J ij σ i σ j and Z is a normalization constant. Here, h i is a local field that biases σ i towards the +1 or −1 configuration, and J ij is a coupling that describes the strength of the interaction between spins i and j. The goal of the inverse Ising problem is to infer the parameters, i.e. h i and J ij , from a set of n observed configurations σ (l) for l = 1 . . . n.
The Ising model presents two obstacles that make the inverse problem quite difficult. First, computing the partition function (or its derivatives) is computationally intractable for large spin systems. As a result, the inverse Ising problem has to be solved approximately. A number of approximate methods for inferring the parameters of the Ising model have been introduced including naive mean field theory [15], the Thouless-Anderson-Palmer (TAP) approximation [16], the isolated spin pair approximation [14], the Sessak-Monasson expansion [17], and others [5,[18][19][20][21][22][23][24][25][26][27]. The second obstacle -overfitting -is more fundamental and generally affects all high dimensional problems in statistical inference.
Overfitting is a simple concept to illustrate for the Ising model. All of the moments, i.e. σ i and σ i σ j , are noisy because they are computed from a finite sample of n observed configurations. Thus, fitting all of the moments exactly necessarily incorporates noise from the finite sample size. As a consequence, parameters obtained by fitting the moments of one dataset may not provide a good description of a new dataset derived from the same system. In general, overfitting becomes a serious problem when the number of parameters (i.e. N 2 ) exceeds the number of independent observations (i.e. n).
A common approach to mitigating the effects of overfitting is to penalize parameters with large values [28]. For example, the Ising model can be 'regularized' by adding an L 2 penalty λ i<j J 2 ij (as in [28]) or an L 1 penalty λ i<j |J ij | (as in [18,21]) to the objective function that describes the fit to the data. Alternatively, the empirical moments can be modified using a 'pseudocount' according to the rules σ i → (1 − α) σ i and σ i σ j → δ ij + (1 − δ ij )(1 − α) σ i σ j [28]. The intuition for why regularization works comes from the law of total variance, which implies that the error in a statistical estimator is error 2 = variance + bias 2 . Thus, a regularization method that decreases the variance in an estimator more than it increases the squared bias will have a smaller error and, therefore, better out-of-sample predictive ability.
The free parameter (λ ≥ 0 or 0 ≤ α ≤ 1) that penalizes large couplings must be carefully chosen to ensure that the regularization method actually improves out-of-sample performance. Typically, the choice of the regularization parameter and subsequent testing of predictive ability are performed using cross-validation. To perform cross-validation, the observed data are randomly partitioned into three mutually exclusive datasets. The first dataset, usually called the 'training' sample, is used for fitting the model. Here, the training sample will be denoted 'F' for 'fitting'. The resulting parameters are used to predict the data contained in the second dataset, called the 'validation' sample ('V'). The optimal value FIG. 1.
for the regularization parameter is chosen to maximize the predictive performance of the model on the validation sample. Finally, the out-of-sample performance of the resulting model is tested by measuring the agreement between the predictions of the model and the data contained in the third dataset, called the 'test' sample ('T').
Focusing on out-of-sample predictive ability provides a different metric than most previous papers on the inverse Ising model, which have studied the ability to reconstruct a coupling matrix from sampled configurations (e.g. [14]). In general, the performance of various algorithms for the inverse Ising model have been studied by defining a set of couplings, simulating configurations from the resulting Ising model, and then asking how well the couplings inferred by various methods reproduce the true J ij 's. By contrast, this work will compare the performance of different approaches to the inverse Ising model by their outof-sample predictive performance using a dataset of images developed for testing machine learning algorithms.
In order to compare the predictive performance of various inference methods for the Ising model, it is necessary to have a metric that quantifies the agreement between the model and the data in the test sample that can be easily computed. Here, I will use the negative log-pseudolikelihood, L data p (h, J) = − log i P (σ i |σ j =i ) data , where the angular brackets denote an average over the spin configurations in dataset 'F', 'V', or 'T'. Recently, Aurell and Ekeberg [18] demonstrated that the couplings of the Ising model can be inferred by minimizing the negative log-pseudolikelihood, which is given by (ignoring constant terms): is the effective field acting on spin i due to the configurations of the other n − 1 spins. The derivatives of L p (h, J) are easy to calculate: Note that these expressions assume J ii = 0 and J ij = J ji . Thus, the Ising model can be fit to the training sample by minimizing L F p (h, J) using gradient descent. While the computing the pseudolikelihood is certainly tractable it is, neverthless, computationally challenging to compute the averages over the training sample at every step of gradient descent if the sample size is large. Therefore, I propose a variational approach to maximum pseudolikelihood inference of the Ising model that is more computationally efficient and, as I will show, performs better out of sample. The variational pseudolikelihood approach minimizes an energy that provides an upper bound on L F p (h, J), given by: The upper bound follows from log(cosh(φ i )) ≤ log(cosh(µ i )) + log(cosh(ν i )) if φ i is Gaussian distributed with mean µ i and standard deviation ν i ; this inequality is illustrated in Fig. 1. Since φ i is a sum of many random variables, I will assume that it is approximately Gaussian distributed with mean After plugging in the estimate for h i , it is only necessary to minimize a function of J given by: The derivatives are given by: where D ij = tanh(ν i )δ ij /ν i and the diagonal elements of J are constrained to J ii = 0. The couplings can be estimated by minimizing Eq. 3 using the gradient descent update J ij (t + 1) = J ij (t) − (G ij (t) + ρG ij (t − 1)), where is a small step size and ρ is a 'momentum'. The variational pseudolikelihood approach was compared to ( [14]): direct pseudolikelihood maximization (with and without an L 2 penalty on the couplings), naive mean field (NMF) inversion (with and without pseudocount regularization), and the isolated spin pair (ISP) approximation (with and without pseudocount regularization). For the regularized methods, the penalty parameters were chosen by cross-validation to minimize L V p (h, J) on the validation sample. The performances of these methods was assessed using training, validation, and test samples con structured from the 'notMNIST' dataset, which consists of 28 × 28 pixel images of the letters A-J compiled from a large number of different fonts (see Appendix). The 'notMNIST' dataset contains roughly 50000 images for each letter, and the samples of each letter are very diverse. Examples of binarized images for letters A-J are shown in Fig. 2. Mutually exclusive training, validation, and test samples were constructed, each containing 500 randomly chosen images for each letter. Table I shows a comparision of these methods according to L F p (h, J) for the training set images, and L T p (h, J) for the test set images. As one would expect, ranking by the pseudolikelihood ensures that direct pseudolikelihood maximization is the best performing method on the training data. However, the variational pseudolikelihood approach presented here outperforms all other methods, including direct pseudolikelihood maximization, when compared by out-of-sample performance on the test set. In addition, a simple implementation of variational pseudolikelihood inference was approximately an order of magnitude faster than direct pseudolikelihood maximization running on a desktop computer, even ignoring the time spent choosing the penalty parameter for L 2 regularized pseudolikelihood.
Comparing the out-of-sample results for the unregularized and regularized versions of maximum pseudolikelihood, NMF, and ISP estimates in Table I demonstrates the importance of regularization for the inverse Ising model. Histograms of the couplings inferred with direct pseudolikelihood maximization, pseudlikelihood maximization with an L 2 penalty, and the variational pseudolikelihood method (Fig. 2) illustrate that the methods which performed better out-of-sample had smaller estimates for the couplings. However, adding the L 2 penalty to the pseudolikelihood resulted in an inferred Ising model with many weak couplings, and no strong couplings. The variational approach, by contrast, produced both strong and weak couplings, which resulted in improved out-of-sample prediction performance relative to other methods. The analysis presented above demonstrates that variational pseudolikelihood inference can be used to train an Ising model so that it makes good out-of-sample predictions. Here, I provide an example to illustrate how the variational pseudolikelihood method can be used in practice. I examined a simple question: what features does the Ising model learn when trained on images of the letters A-J? Using the variational pseudolikelihood approach, it is simple to infer coupling matrices with a specific structure, such as that of a Hopfield neural network [5,21,26,29]. The couplings in a Hopfield network are described by a Hebbian rule such that where ξ (s) for s = 1 . . . p are 'patterns' stored in the memory. These patterns describe features in the correlation structure of the data. The derivatives of the variational energy with respect to the patterns can be computed using the chain rule: Thus, the patterns can be inferred quite easily by minimizing the variational energy using gradient descent. Note, however, that the patterns are unconstrained with respect to orthogonal transformations [30], and the symmetry is broken by the initial conditions of the optimization.
The letters in the notMNIST dataset are quiet diverse, and hundreds of patterns are required to approach the performance of the unconstrained variational pseudolikelihood method (Fig. 4a). Nevertheless, the p = 10 Hopfield model provides a useful visual example. Starting from randomly initialized patterns (shown in the top row of Fig. 4b), the patterns converge after less than 1000 steps of gradient descent to images with recognizable features of letters (shown in the second row of Fig. 4b).
The notMNIST images present a challenging dataset for the Ising model. State-of-the-art machine learning methods for modeling similar types of images include a variety of pre-processing and fine tuning steps [31], and the best performing algorithms are based on deep neural networks that are believed to excel at capturing higher order correlation structure in the data that cannot be described by the Ising model [32]. Nevertheless, Fig.  4 clearly demonstrates an Ising model based on a few Hopfield patterns is able to capture recognizable features from the data.
In summary, I presented a variational approach to maximum pseudolikelihood inference of the Ising model. The variational pseudolikelihood method outperformed all other tested methods at describing the out-of-sample correlation structure in a dataset consisting of images of the letters A-J from various computer fonts. This method has many characteristics that make it attractive in practice: it is fast, accurate, can be used with large spin systems, and can easily be adapted to learn coupling matrices with specific structures such as Hopfield models. Moreover, it should be easy to extend the approach presented here to construct variational pseudolikelihood algorithms for inferring Potts models [6], Restricted Boltzmann Machines [33], and other inverse models for spin systems.
I would like to thank Alex H. Lang, Thomas C. Lang, Pankaj Mehta, Ilya Nemenman, and Javad Noorbakhsh for helpful discussions. Funding for this research was provided by a Simons Investigator Award to Pankaj Mehta.
Appendix: The 'notMNIST' dataset consists of 28×28 pixel images of the letters A-J compiled from a large number of different fonts. The images are available at http://yaroslavvb.blogspot.com/2011/09/notmnistdataset.html. Each image in the original dataset is grayscale, with pixels taking on values between 0 and 255. The images were binarized in the simplest way by defining pixels with intensity less than or equal to 127 as spin −1, and those with intensity greater than 127 as spin +1. Finally, each image was flattened into a vector of N = 784 spins. * charleskennethfisher@gmail.com