Inverse statistical problems: from the inverse Ising problem to data science

TL;DR

This paper reviews inverse statistical physics problems, especially the inverse Ising problem, highlighting methods for inferring model parameters from data across various scientific fields.

Contribution

It provides a comprehensive overview of methods and applications for the inverse Ising problem, including equilibrium and non-equilibrium cases, with emphasis on recent advances and practical algorithms.

Findings

Pseudolikelihood is a powerful method for inverse Ising inference.

Applications include neural network reconstruction and protein structure prediction.

Various approximate solutions are effective for different data regimes.

Abstract

Inverse problems in statistical physics are motivated by the challenges of `big data' in different fields, in particular high-throughput experiments in biology. In inverse problems, the usual procedure of statistical physics needs to be reversed: Instead of calculating observables on the basis of model parameters, we seek to infer parameters of a model based on observations. In this review, we focus on the inverse Ising problem and closely related problems, namely how to infer the coupling strengths between spins given observed spin correlations, magnetisations, or other data. We review applications of the inverse Ising problem, including the reconstruction of neural connections, protein structure determination, and the inference of gene regulatory networks. For the inverse Ising problem in equilibrium, a number of controlled and uncontrolled approximate solutions have been developed in the statistical mechanics community. A particularly strong method, pseudolikelihood, stems from statistics. We also review the inverse Ising problem in the non-equilibrium case, where the model parameters must be reconstructed based on non-equilibrium statistics.

Tables

Table 1. Table 1: Classification of reconstruction methods based on the thermodynamic potentials used and the approximations employed to evaluate them. Pseudolikelihood maximisation falls outside this classification scheme, as it is not based on an approximation of the Boltzmann measure.

associated potential	exact	variational approximation	perturbative expansion
$F(𝑱,𝒉)$	convex nonlinear optimisation
$G(𝑱,𝒎)$		mean-field, Bethe–Peierls	Plefka, mean-field, TAP, Bethe–Peierls
$S(𝝌,𝒎)$		independent-pair approximation, Cocco–Monasson	Sessak–Monasson

Equations

$$P(\{s_{i}\})=\frac{1}{Z}\exp\left[\sum_{i}h_{i}s_{i}+\sum_{i<j}J_{ij}s_{i}s_{j}\right]\$$

$$H(\{s_i\})=-\sum _i{h}_i{s}_i-\sum _{i<j}{J}_{ij}{s}_i{s}_j$$

$$S[P]=-\sum _{\{s_i\}}P(\{s_i\})\ln P(\{s_i\})$$

$$C\frac{d{V}_i}{dt}=\sum _{j\ne i}{J}_{ij}\sum _{l}K(t-t_{jl})+I_i-g{V}_i+{\xi }_i(t),$$

$$P(\{x_i\})=\exp \left[\sum _{i\le j}{J}_{ij}{x}_i{x}_j+\sum _i{h}_i{x}_i\right]/Z,$$

$$\text{sign}(x_i^{t+1})=\{\begin{array}{ll}1& \text{if }{\sum }_j{J}_{ij}{x}_j^t>\kappa \\ -1& \text{if }{\sum }_j{J}_{ij}{x}_j^t<\kappa \end{array},$$

$$H=-\sum _{i<j}{J}_{ij}(s_i,s_j)-\sum _i{h}_i(s_i).$$

$$D{I}_{ij}=\sum _{s_i,s_j}{p}_{ij}(s_i,s_j)\ln \left(\frac{p_{ij}(s_i,s_j)}{p_i(s_i)p_j(s_j)}\right).$$

$$P(\{{\vec{s}}_i\})=\exp \left[\sum _{i,j}{J}_{ij}{\vec{s}}_i\cdot {\vec{s}}_j\right]/Z,$$

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

$$p(\mathbf{s})=\frac{1}{Z}{e}^{-H_{\mathbf{J},\mathbf{h}}(\mathbf{s})},$$

$$Z(\mathbf{J},\mathbf{h})=\sum _{\mathbf{s}}{e}^{-H_{\mathbf{J},\mathbf{h}}(\mathbf{s})}.$$

$$⟨Q(\mathbf{\sigma })⟩\equiv \sum _{\mathbf{s}}p(\mathbf{s})Q(\mathbf{s}),$$

$${\theta }^{\mathrm{ML}}={\mathrm{argmax}}_{\theta }p(x_1,x_2,\dots ,x_M|\theta )$$

$$p(\theta |x_1,x_2,\dots ,x_M)$$

$$L_{\mathrm{D}}(\mathbf{J},\mathbf{h})$$

$${⟨Q⟩}^{\mathrm{D}}=\frac{1}{M}\sum _{\mu }Q({\mathbf{s}}^{\mu }).$$

$$\mathrm{KL}(p|q)=\sum _{\mathbf{s}}p(\mathbf{s})\ln \frac{p(\mathbf{s})}{q(\mathbf{s})},$$

$$\mathrm{KL}(p^{\mathrm{D}}|p)$$

$$\ln p(\mathbf{J},\mathbf{h}|\mathrm{D})=M{L}_{\mathrm{D}}(\mathbf{J},\mathbf{h})-\gamma \sum _{i<j}|J_{ij}|$$

$$\{{\mathbf{J}}^{\mathrm{ML}},{\mathbf{h}}^{\mathrm{ML}}\}=\mathrm{argmax}{L}_{\mathrm{D}}(\mathbf{J},\mathbf{h})$$

$$\frac{\partial L_{\mathrm{D}}}{\partial h_i}(\mathbf{J},\mathbf{h})$$

$$\frac{\partial L_{\mathrm{D}}}{\partial J_{ij}}(\mathbf{J},\mathbf{h})$$

$$⟨{\sigma }_i⟩$$

$$⟨{\sigma }_i{\sigma }_j⟩$$

$$h_i^{n+1}$$

$$J_{ij}^{n+1}$$

$$p(\mathbf{s})=\frac{1}{Z(\mathbf{\lambda })}{e}^{{\sum }_k{\lambda }_k{Q}_k(\mathbf{s})}.$$

$$-\frac{{\partial }^2{L}_{\mathrm{D}}}{\partial {\lambda }_i\partial {\lambda }_j}(\mathbf{\lambda })=⟨Q_i{Q}_j⟩-⟨Q_i⟩⟨Q_j⟩.$$

$$W=\frac{M!}{{\prod }_{r=1}^R{n}_r!}$$

$$\ln W/M\approx -\sum _{r=1}^R{q}_r\ln q_r.$$

$$\frac{M!}{{\prod }_{r=1}^R{n}_r!}\prod _{r=1}^R(1/R{)}^{n_r}.$$

Full text

1. I Introduction and applications

2. I.1 Modelling neural firing patterns and the reconstruction of neural connections

3. I.2 Reconstruction of gene regulatory networks

4. I.3 Protein structure determination

5. I.4 Fitness landscape inference

6. I.5 Combinatorial antibiotic treatment

7. I.6 Interactions between species and between individuals

8. I.7 Financial markets

9. II Equilibrium reconstruction

10. II.0.1 Definition of the problem

11. II.1 Maximum likelihood

12. II.1.1 Exact maximization of the likelihood

13. II.1.2 Uniqueness of the solution

14. II.1.3 Maximum entropy modelling

15. II.1.4 Information theoretic bounds on graphical model reconstruction

16. II.1.5 Thermodynamics of the inverse Ising problem

17. II.1.6 Variational principles

18. II.1.7 Mean-field theory

19. II.1.8 The Onsager term and TAP reconstruction

20. II.1.9 Couplings without a loop: mapping to the minimum spanning tree problem

21. II.1.10 The Bethe–Peierls ansatz

22. II.1.11 Belief propagation and susceptibility propagation

23. II.1.12 The independent-pair approximation and the Cocco–Monasson adaptive-cluster expansion

24. II.1.13 The Plefka expansion

25. II.1.14 The Sessak–Monasson small-correlation expansion

26. II.2 Logistic regression and pseudolikelihood

27. II.3 Comparison of the different approaches

28. II.4 Reconstruction and non-ergodicity

29. II.5 Parameter reconstruction and criticality

30. III Non-equilibrium reconstruction

31. III.1 Dynamics of the Ising model

32. III.1.1 Sequential Glauber dynamics

33. III.1.2 Parallel Glauber dynamics

34. III.2 Reconstruction from time series data

35. III.2.1 Maximisation of the likelihood

36. III.2.2 Mean-field theory of the non-equilibrium steady state

37. III.2.3 The Gaussian approximation

38. III.2.4 Method comparison

39. III.3 Outlook: Reconstruction from the steady state

40. IV Conclusions

I Introduction and applications

The primary goal of statistical physics is to derive observable quantities from microscopic laws governing the constituents of a system. In the example of the Ising model, the starting point is a model describing interactions between elementary magnets (spins), the goal is to derive observables such as spin magnetisations and correlations.

In an inverse problem, the starting point is observations of some system whose microscopic parameters are unknown and to be discovered. In the inverse Ising problem, the interactions between spins are not known to us, but we want to learn them from measurements of magnetisations, correlations, or other observables. In general, the goal is to infer the parameters describing a system (for instance, its Hamiltonian) from extant data. To this end, the relationship between microscopic laws and observables needs to be inverted.

In the last two decades, inverse statistical problems have arisen in different contexts, sparking interest in the statistical physics community in taking the path from model parameters to observables in reverse. The areas where inverse statistical problems have arisen are characterized by (i) microscopic scales becoming experimentally accessible and (ii) sufficient data storage capabilities being available. In particular, the biological sciences have generated several inverse statistical problems, including the reconstruction of neural and gene regulatory networks and the determination of the three-dimensional structure of proteins. Technological progress is likely to to open up further fields of research to inverse statistical analysis, a development that is currently described by the label ‘big data’.

In physics, inverse statistical problems also arise when we need to design a many-body system with particular desired properties. Examples are finding the potentials that result in a particular single-particle distribution kunkin1969inverse ; chayes1984a , interaction parameters in a binary alloy that yield the observed correlations maysenholder1987determination , the potentials between atoms that lead to specific crystal lattices zhang2013a , or the parameters of a Hamiltonian that lead to a particular density matrix changlani2015 . In the context of soft matter, a question is how to design a many-body system that will self-assemble into a particular spatial configuration or has particular bulk properties rechtsman2006a ; torquato2009a . In biophysics, we may want to design a protein that folds into a specified three-dimensional shape kuhlman2003a . For RNA, even molecules with more than one stable target structure are possible flamm2001a . As a model of such design problems, distasio2013a ; marcotte2013a study how to find the parameters of Ising Hamiltonian with a prescribed ground state.

In all these examples, ‘spin’ variables describe microscopic degrees of freedom particular to a given system, for instance, the states of neurons in a neural network. The simplest description of these degrees of freedom in terms of random binary variables then leads to Ising-type spins. In the simplest non-trivial scenario, correlations between the ’spins’ are generated by pairwise couplings between the spins, leading to an Ising model with unknown parameters (couplings between the spins and magnetic fields acting on the spins). In many cases of interest, the couplings between spins will not all be positive, as is the case in a model of a ferromagnet. Nor will they couplings conform to a regular lattice embedded in some finite-dimensional space.

For a concrete example, we look at a system of $N$ binary variables (Ising spins) $\{s_{i}\},i=1,\,\ldots,N$ with $s_{i}=\pm 1$. These spins are coupled by pairwise couplings $J_{ij}$ and are subject to external magnetic fields $h_{i}$.

[TABLE]

is the Boltzmann equilibrium distribution $P(\{s_{i}\})=e^{-H(\{s_{i}\})}/Z$, where we have subsumed temperature into the couplings and fields. (We will discuss this choice in section II.0.1). The Hamiltonian

[TABLE]

specifies the energy of the spin system as a function of the microscopic spin variables, local fields, and pairwise couplings. The inverse Ising problem is the determination of the couplings $J_{ij}$ and local fields $h_{i}$, given a set of $M$ observed spin configurations. Depending on the particular nature of the system at hand, the restriction to binary variables or pairwise interactions may need to be lifted, or the functional form of the Hamiltonian may be altogether different from the Ising Hamiltonian with pairwise interactions (2). For non-equilibrium systems, the steady state is not even described by a Boltzmann distribution with a known Hamiltonian. However, the basic idea remains the same across different types of inverse statistical problems: even when the frequencies of spin configurations may be under-sampled, the data may be sufficient to infer at least some parameters of a model.

The distribution (1) is well known not only as the equilibrium distribution of the Ising model. It is also the form of the distribution which maximizes the (Gibbs) entropy

[TABLE]

under the constraint that $P(\{s_{i}\})$ is normalized and has particular first and second moments, that is, magnetisations and correlations. We will discuss in section II.1.3 how this distribution emerges as the ‘least biased distribution’ of a binary random variable with prescribed first and second moments jaynes1957a . The practical problem is then again an inverse one: to find the couplings $J_{ij}$ and local fields $h_{i}$ such that the first and second moments observed under the Boltzmann distribution (1) match the mean values of $s_{i}$ and $s_{i}s_{j}$ in data. In settings where third moments differ significantly from the prediction of (1) based on the first two moments, one may need to construct the distribution of maximum entropy given the first three moments, leading to three-spin interactions in the exponent of (1).

Determining the parameters of a distribution such as (1) is always a many-body problem: changing a single coupling $J_{ij}$ generally affects correlations between many spin variables, and conversely a change in the correlation between two variables can change the values of many inferred couplings. The interplay between model parameters and observables is captured by a statistical mechanics of inverse problems, where the phase space consists of quantities normally considered as fixed model parameters (couplings, fields). The observables, such as spin correlations and magnetisations on the other hand, are taken to be fixed, as they are specified by empirical observations. Such a perspective is not new to statistical physics; the analysis of neural networks in the seventies and eighties of the last century led to a statistical mechanics of learning watkin1993a ; hertz199a ; engel2001a , where the phase space is defined by the set of possible rules linking the input into a machine with an output. The set of all rules compatible with a given set of input/output relations then defines a statistical ensemble. In the inverse statistical problems, however, there are generally no explicit rules linking the input and output, but data with different types of correlations or other observations, which are to be accounted for in a statistical model.

Inverse statistical problems fall into the broader realm of statistical inference mackay2003 ; bishop2006 , which seeks to determine the properties of a probability distribution underlying some observed data. The problem of inferring the parameters of a distribution such as (1) is known under different names in different communities; also emphasis and language differ subtly across communities.

• •

In statistics, an inverse problem is the inference of model parameters from data. In our case, the problem is the inference of the parameters of the Ising model from observed spin configurations. A particular subproblem is the inference of the graph formed by the non-zero couplings of the Ising model, termed graphical model selection or reconstruction. In the specific context of statistical models on graphs (graphical models), the term inference describes the calculation of the marginal distribution of one or several variables. (A marginal distribution describes the statistics of one or several particular variables in a many-variable distribution, for example, $P(x_{1})=\sum_{x_{2},x_{3}}P(x_{1},x_{2},x_{3})$.)

• •

In machine learning, a frequent task is to train an artificial neural network with symmetric couplings such that magnetisations and correlations of the artificial neurons match the corresponding values in the data. This is a special case of what is called Boltzmann machine learning; the general case also considers so-called hidden units, whose values are unobserved Ackley1985a .

• •

In statistical physics, much effort has been directed towards estimating the parameters of the Ising model given observed values of the magnetisation and two-point correlations. As we will see in section II.1, this is a hard problem from an algorithmic point of view. Recently, threshold phenomena arising in inference problems have attracted much interest from the statistical physics community, due to the link between phase transitions and the boundaries separating different regimes of inference problems, for instance solvable and unsolvable problems, or easy and hard ones mezard2009information ; ZdeborovaKrzakala2015 .

Common theme across different applications and communities is the inference of model parameters given observed data or desired properties. In this review we will focus on a prototype inverse statistical problem: the inverse Ising problem and its close relatives. Many of the approaches developed for this problem are also readily extended to more general scenarios. We will start with a discussion of applications of the inverse Ising problem and related approaches in biology, specifically the reconstruction of neural and genetic networks, the determination of three-dimensional protein structures, the inference of fitness landscapes, the bacterial responses to combinations of antibiotics, and flocking dynamics. We will find that these applications define two distinct settings of the inverse Ising problem; equilibrium and non-equilibrium. Part II of this review treats the inverse Ising problem in an equilibrium setting, where the couplings between spins are symmetric. Detailed balance holds and results from equilibrium statistical physics can be used. This setting arises naturally within the context of maximum entropy models, which seek to describe the observed statistics of configurations with a simplified effective model capturing, for instance, collective effects. We introduce the basics of statistical inference and maximum entropy modelling, discuss the thermodynamics of the inverse Ising problem, and review different approaches to solve the inverse Ising problem, pointing out their connections and comparing the resulting parameter reconstructions. Part III of this review considers asymmetric coupling matrices, where in the absence of detailed balance couplings can be reconstructed from time series, from data on perturbations of the system, or from detailed knowledge of the non-equilibrium steady state.

We now turn to applications of the inverse Ising problem, which mostly lie in the analysis of high-throughput data from biology. One aim of inverse statistical modelling is to find the parameters of a microscopic model to describe this data. A more ambitious aim is achieved when the parameters of the model are informative about the processes which produced the data, this is, when some of the mechanisms underlying the data can be inferred. The data is large-scale measurements of the degrees of freedom of some system. In the language of statistical physics these describe the micro-states of a system: states of neurons, particular sequences of DNA or proteins, or the concentration levels of RNA. We briefly introduce some of the experimental background of these measurements, so their potential and the limitations can be appreciated. The models are simple models of the microscopic degrees of freedom. In the spirit of statistical physics, these models are simple enough so the parameters can be computed given the data, yet sufficiently complex to reproduce some of the statistical interdependences of the observed microscopic degrees of freedom. The simplest case, consisting of binary degrees of freedom with unknown pairwise couplings between them, leads to the inverse Ising problem, although we will also discuss several extensions.

I.1 Modelling neural firing patterns and the reconstruction of

neural connections

Neurons can exchange information by generating discrete electrical pulses, termed spikes, that travel down nerve fibres. Neurons can emit these spikes at different rates, a neuron emitting spikes at a high rate is said to be ‘active’ or ‘firing’, a neuron emitting spikes at a low rate or not at all is said to be ‘inactive’ or ‘silent’. The measurement of the activity of single neurons has a long history starting in 1953 with the development of micro-electrodes for recording dowben1953a . Multi-electrodes were developed, allowing to record multiple simultaneous neural signals independently over long time periods spira2013a ; obien2014a . Such data presents the intriguing possibility to see elementary brain function emerge from the interplay of a large number of neurons.

However, even when vast quantities of data are available, the different configurations of a system are still under-sampled in most cases. For instance, consider $N$ neurons, each of which can be either active (firing) or inactive (silent). Given that the firing patterns of thousands of neurons can be recorded simultaneously schwarz2014a , the number of observations $M$ will generally be far less than the total number of possible neural configurations, $M\ll 2^{N}$. For this reason, a straightforward statistical description that seeks to determine directly the frequency with which each configuration appears will likely fail.

On the other hand, a feasible starting point is a simple distribution, whose parameters can to be determined from the data. For a set of $N$ binary variables, this might be a distribution with pairwise interactions between the variables. In Schneidman2006a , Bialek and collaborators applied such a statistical model to neural recordings. Dividing time into small intervals of duration $\Delta\tau=20\text{~{}}\mathrm{ms}$ induces a binary representation of neural data, where each neuron $i$ either spikes during a given interval ($s_{i}=1$) or it does not ($s_{i}=0$). The joint statistics observed in $40$ neurons in the retina of a salamander was modelled by an Ising model (1) with magnetic fields and pairwise symmetric couplings. Rather than describing the dynamics of neural spikes, this model describes the correlated firing of different neurons over the course of the experiment. The symmetric couplings $J_{ij}$ in (1) describe statistical dependencies, not physical connections. The synaptic connections between neurons, on the other hand, are generally not symmetric.

In this context, the distribution (1) can be viewed as the form of the maximum entropy distribution over neural states, given the observed one- and two-point correlations Schneidman2006a . In Shlens2006a ; Cocco2009a , a good match was found between the statistics of three neurons predicted by (1) and the firing patterns of the same neurons in the data. This means that the model with pairwise couplings provides a statistical description of the empirical data, one that can even be used to make predictions. Similar results were obtained also from cortical cells in cell cultures tang2008a .

The mapping from the neural data to a spin model rests on dividing time into discrete bins of duration $\Delta\tau$. A different choice of this interval would lead to different spin configurations; in particular changing $\Delta\tau$ affects the magnetisation of all spins by altering the number of intervals in which a neuron fires. In roudinirenberglatham2009a , Roudi, Nirenberg and Latham show that the pairwise model (1) provides a good description of the underlying spin statistics (generated by neural spike trains), provided $N\,\Delta\tau\,\nu\ll 1$, where $\nu$ is the average firing rate of neurons. Increasing the bin size beyond this regime leads to an increase in bins where multiple neurons fire, as a result couplings beyond the pairwise couplings in (1) can become important.

As a minimal model of neural correlations, the statistics (1) has been extended in several ways. Tkačik et al. tkavcik2010 and Granot-Atedgi et al. granot2013a consider stimulus-dependent magnetic fields, that is, fields which depend on the stimulus presented to the experimental subject at a particular time of the experiment. Ohiorhenuan et al. looks at stimulus-dependent couplings ohiorhenuan2010a . When the number of neurons increases to $\sim 100$, limitations of the pairwise model (1) become apparent, which has be addressed by adding additional terms coupling triplets, etc. of spins in the exponent of the Boltzmann measure (1) ganmor2011a .

The statistics (1) serves as a description of the empirical data: the couplings between spins in the Hamiltonian (2) do not describe physical connections between the neurons. The determination of the network of neural connections from the observed neural activities is thus a different question. Simoncelli and collaborators paninski2004a ; pillow2008a and Cocco, Leibler, and Monasson Cocco2009a use an integrate-and-fire model burkitt2006review to infer how the neurons are interconnected on the basis of time series of spikes in all neurons. In such a model, the membrane potential of neuron $i$ obeys the dynamics

[TABLE]

where the first term on the right-hand side encodes the synaptic connections $J_{ij}$ and a memory kernel $K$; $t_{jl}$ specifies the time at which neuron $j$ emitted its $l^{\textrm{th}}$ spike. The remaining terms describe a background current, voltage leakage, and white noise. Finding the synaptic connections $J_{ij}$ that best describe a large set of neural spike trains is a computational challenge; paninski2006a ; Cocco2009a develop an approximation based on maximum likelihood, see section II.1. A related approach based on point processes and generalized linear models (GLM) is presented in truccolo2005point . We will discuss this problem of inferring the network of connections the context of the non-equilibrium models in section III.

Neural recordings give the firing patterns of several neurons over time. These neurons may have connections between them, but they also receive signals from neural cells whose activity is not recorded macke2011a . In tyrcha2013 , the effect of connections between neurons is disentangled from correlations arising from shared non-stationary input. This raises the possibility that the correlations described by the pairwise model (1) in Schneidman2006a and related works originate from a confounding factor (connections to a neuron other than those whose signal is measured), rather than from connections between recorded neurons kulkarni2007common .

I.2 Reconstruction of gene regulatory networks

The central dogma of molecular biology is this: Proteins are macromolecules consisting of long chains of amino acids. The particular sequence of a protein is encoded in DNA, a double-stranded helix of complementary nucleotides. Specific parts of DNA, the genes, are transcribed by polymerases, producing a single-stranded copy called m(essenger)RNAs, which are translated by ribosomes, usually multiple times, to produce proteins.

The process of producing protein molecules from the DNA template by transcription and translation is called gene expression. The expression of a gene is tightly controlled to ensure that the right amounts of proteins are produced at the right time. One important control mechanism is transcription factors, proteins which affect the expression of a gene (or several) by binding to DNA near the transcription start site of that gene. This part of DNA is called the regulatory region of a gene. A target gene of a transcription factor may in turn encode another transcription factor, leading to a cascade of regulatory events. To add further complications, the binding of multiple transcription factors in the regulatory region of a gene leads to combinatorial control exerted by several transcription factors on the expression of a gene buchlergerlandhwa2003a . Can the regulatory connections between genes be inferred from data on gene expression, that is, can we learn the identity of transcription factors and their targets?

Over the last decades, the simultaneous measurement of expression levels of all genes have become routine. At the centre of this development are two distinct technological advances to measure mRNA levels. The first is microarrays, consisting of thousands of short DNA sequences, called probes, grafted to the surface of a small chip. After converting the mRNA in a sample to DNA by reverse transcription, cleaving that DNA into short segments, and fluorescently labelling the resulting DNA segments, fluorescent DNA can bind to its complementary sequence on the chip. (Reverse transcription converts mRNA to DNA, a process which requires a so-called reverse transcriptase as an enzyme.) The amount of fluorescent DNA bound to a particular probe depends on the amount of mRNA originally present in the sample. The relative amount of mRNA from a particular gene can then be inferred from the fluorescence signal at the corresponding probes hoheisel2006a . A limitation of microarrays is the large amount of mRNA required: The mRNA sample is taken from a population of cells. As a result, cell-to-cell fluctuations of mRNA concentrations are averaged over. To obtain time series, populations of cells synchronized to approximately the same stage in the cell cycle are used gasch2000a .

The second way to measure gene expression levels is also based on reverse transcription of mRNA, followed by high-throughput sequencing of the resulting DNA segments. Then the relative mRNA levels follow directly from counts of sequence reads wang2009a . Recently, expression levels even in single cells have been measured in this way wills2013a . In combination with barcoding (adding short DNA markers to identify individual cells), $10^{4}$ cells can have their expression profiled individually in a single sequencing run macosko2015a . Such data may allow, for instance, the analysis of the response of target genes to fluctuations in the concentration of transcription factors. However, due to the destructive nature of single-cell sequencing, there may never be single-cell data that give time series of genome-wide expression levels.

Unfortunately, cellular concentrations of proteins are much harder to measure than mRNA levels. As a result, much of the literature focuses on mRNA levels, neglecting the regulation of translation. Advances in protein mass-spectrometry picotti2013a may lead to data on both mRNA and protein concentrations. This data would pose the additional challenge of inferring two separate levels of gene regulation: gene transcription from DNA to mRNA and translation from mRNA to proteins.

As in the case of neural data discussed in the preceding section, gene expression data presents two distinct challenges: (i) finding a statistical description of the data in terms of suitable observables and (ii) inferring the underlying regulatory connections. Both these problems have been addressed extensively in the machine learning and quantitative biology communities. Clustering of gene expression data to detect sets of genes with correlated expression levels has been used to detect regulatory relationships. A model-based approach to the reconstruction of regulatory connections is Boolean networks. Boolean networks assign binary states to each gene (gene expression on/off), and the state of a gene at a given time depends on the state of all genes at a previous time through a set of logical functions assigned to each gene. See Dhaeseleer2000a for a review of clustering and Boolean network inference and hickman2009a for a review of Boolean network inference.

A statistical description that has also yielded insight into regulatory connections is Bayesian networks. A Bayesian network is a probabilistic model describing a set of random variables (expression levels) through conditional dependencies described by a directed acyclic graph. Learning both the structure of the graph and the statistical dependencies is a hard computational problem, but can capture strong signals in the data that are often associated with a regulatory connection. In principle, causal relationships (like the regulatory connections) can be inferred, in particular if the regulatory network contains no cycles. For reviews, see friedman2004a ; koller2009a . Both Boolean or Bayesian networks have been applied to measurements of the response of expression levels to external perturbations of the regulatory network or of expression levels, see  ideker2000a ; peer2001a . A full review of these methods is beyond the scope of this article, instead we focus on approaches related to the inverse Ising problem.

For a statistical description of gene expression levels, lezon2006 applied a model with pairwise couplings

[TABLE]

fitted to gene expression levels. The standard definition of expression levels $\{x_{i}\}$ is $\log_{2}$-values of fluorescence signals with the mean value for each gene subtracted. Since (5) is a multi-variate Gaussian distribution, the matrix of couplings $J_{ij}$ must be negative definite. These couplings can be inferred simply by inverting the matrix of variances and covariances of expression levels. In lezon2006 , the resulting couplings $J_{ij}$ were then used to identify hub genes which regulate many targets. The same approach was used in locasale2009 to analyse the cellular signalling networks mediated by the phosphorylation of specific sites on different proteins. Again, the distribution (5) can be viewed as a maximum entropy distribution for continuous variables with prescribed first and second moments. This approach is also linked to the concept of partial correlations in statistics baba2004a ; krumsiek2011a .

Again the maximum-entropy distribution (5) has symmetric couplings between expression levels, whereas the network of regulatory interactions is intrinsically asymmetric. One way to infer the regulatory connections is time series sima2009a . Bailly2010b uses expression levels measured at different times to infer the regulatory connections, based on a minimal model of expression dynamics with asymmetric regulatory connections between pairs of genes. In this model, expression levels $x_{i}^{t}$ at successive time intervals $t$ obey

[TABLE]

where $\kappa$ is a threshold. The regulatory connections $J_{ij}$ are taken to be discrete, with the values $-1,1,0$ denoting repression, activation and no regulation of gene $i$ by the product of gene $j$. The matrix of connections is then inferred based on Bayes theorem (see section II.1) and an iterative algorithm for estimating marginal probabilities (message passing, see section II.1.11).

A second line of approach that can provide information on regulatory connections is perturbations tegner2003a . An example is gene knockdowns, where the expression of one or more genes is reduced, by introducing small interfering RNA (siRNA) molecules into the cell dorsett2004a or by other techniques. siRNA molecules can be introduced into cells from the outside; after various processing steps they lead to the cleavage of mRNA with a complementary sequence, which is then no longer available for translation. If that mRNA translates to a transcription factor, all targets of that transcription factor will be upregulated or downregulated (depending on whether the transcription factor acted as a repressor or an activator, respectively). Knowing the responses of gene expression levels to a sufficient number of such perturbations allows the inference of regulatory connections. molinelli2013a considers a model of gene expression dynamics based on continuous variables $x_{i}$ evolving deterministically as $\partial_{t}x_{i}=a_{i}\tanh(\sum_{j}J_{ij}x_{j})-c_{i}x_{i}$. The first term describes how the expression level of gene $j$ affects the rate of gene expression of gene $i$ via the regulatory connection $J_{ij}$, the second term describes mRNA degradation. The stationary points of this model shift in response to perturbations of expression levels of particular genes (for instance through knockdowns), and these changes depend on regulatory connections. In molinelli2013a , the regulatory connections are inferred from perturbation data, again using belief propagation.

I.3 Protein structure determination

Tremendous efforts have been made to determine the three-dimensional structure of proteins. A linear amino acid chain folds into a convoluted shape, the folded protein, thus bringing amino acids into close physical proximity that are separated by a long distance along the linear sequence.

Due to the number of proteins (several thousand per organism) and the length of individual proteins (hundreds of amino acid residues), protein structure determination is a vast undertaking. However, the rewards are also substantial. The three-dimensional structure of a protein determines its physical and chemical properties, and how it interacts with other cellular components: broadly, the shape of a protein determines many aspects of its function. Protein structure determination relies on crystallizing proteins and analysing the X-ray diffraction pattern of the resulting solid. Given the experimental effort required, the determination of a protein’s structure from its sequence alone has been a key challenge to computational biology for several decades dill2012a ; deJuan2013a . The computational approach models the forces between amino acids in order to find the low-energy structure a protein in solution will fold into. Depending on the level of detail, this approach requires extensive computational resources.

An attractive alternative enlists evolutionary information: Suppose that we have at our disposal amino acid sequences of a protein as it appears in different related species (so-called orthologs). While the sequences are not identical across species, they preserve to some degree the three-dimensional shape of the protein. Suppose a specific pair of amino acids that interact strongly with each other and bring together parts of the protein that are distal on the linear sequence. Replacing this pair with another, equally strongly interacting pair of amino acids would change the sequence, but leave the structure unchanged. For this reason, we expect sequence differences across species to reflect the structure of the protein. Specifically, we expect correlations of amino acids in positions that are proximal to each other in the three-dimensional structure. In turn, the correlations observed between amino acids at different positions might allow to infer which pairs of amino acids are proximal to each other in three dimensions (the so-called contact map). The use of such genomic information has recently lead to predictions of the three-dimensional structure of many protein families inaccessible to other methods ovchinnikov2017 , for a review see coccoetal2017 .

Early work looked at the correlations as a measure of proximity gobel1994 ; lockless1999a ; socolich2005 ; halabi2009a . However correlations are transitive; if amino acids at sequence sites $i$ and $j$ are correlated due to proximity in the folded state, and $j$ and $k$ are correlated for some reason, $i$ and $k$ will also exhibit correlations, which need not stem from proximity. This problem is addressed by an inverse approach aimed at finding the set of pairwise couplings that lead to the observed correlations or sequences Weigt2009a ; burger2010 ; Morcos2011a ; marks2011a ; sulkowska2012 ; dago2012 ; hopf2012a ; cocco2013a ; ekeberg2013 . Since each sequence position can be taken up by one of $20$ amino acids or a gap in the sequence alignment, there are $21^{2}$ correlations at each pair of sequence positions. In Weigt2009a ; Morcos2011a a statistical model with pairwise interactions is formulated, based on the Hamiltonian

[TABLE]

This Hamiltonian depends on spin variables $s_{i}$, one for each sequence position $i=1,\ldots,N$. Each spin variable can take on one of $21$ values, describing the $20$ possible amino acids at that sequence position as well as the possibility of a gap (corresponding to an extra amino acid inserted in a particular position in the sequences of other organisms). Each pair of amino acids $s_{i},s_{j}$ in sequence position $i,j$ contributes $J_{ij}(s_{i},s_{j})$ to the energy. The inverse problem is to find the couplings $J_{ij}(A,B)$ for each pair of sequence positions $i,j$ and pair of amino acids $A,B$, as well as field $h_{i}(A)$, such that the amino acid frequencies and correlations observed across species are reproduced. The sequence positions with strong pairwise couplings are then predicted to be proximal in the protein structure. A simple measure of the coupling between sequence positions is the matrix norm (Frobenius norm) $\sum_{s_{i},s_{j}}(J_{ij}(s_{i},s_{j}))^{2}$. The so-called direct information Weigt2009a is an alternative measure based on information theory. A two-site model is defined with $p_{ij}(s_{i},s_{j})=\exp\{J_{ij}(s_{i},s_{j})+h_{i}(s_{i})+h_{j}(s_{j})\}/Z_{ij}$. Direct information is the mutual information between the two-site model and a model without correlations between the amino acids

[TABLE]

The Boltzmann distribution resulting from (7) can be viewed as the maximum entropy distribution with one- and two-point correlations between amino acids in different sequence positions determined by the data. There is no reason to exclude higher order terms in the Hamiltonian (7) describing interactions between triplets of sequence positions, although the introduction of such terms may lead to overfitting. Also, fitting the Boltzmann distribution (7) to sequence data uses no prior information on protein structures; for this reason it is called an unsupervised method. Recently, neural network models trained on sequence data and protein structures (supervised learning) have been very successful in predicting new structures jones2015metapsicov ; wang2017accurate .

The maximum entropy approach to structure analysis is not limited to evolutionary data. In zhangwolynes2015a Zhang and Wolynes analyse chromosome conformation capture experiments and use the observed frequency of contacts between different parts of a chromosome in a maximum entropy approach to predict the structure and topology of the chromosomes.

I.4 Fitness landscape inference

The concept of fitness lies at the core of evolutionary biology. Fitness quantifies the average reproductive success (number of offspring) of an organism with a particular genotype, i.e., a particular DNA sequence. The dependence of fitness on the genotype can be visualized as a fitness landscape in a high-dimensional space, where fitness specifies the height of the landscape. As the number of possible sequences grows exponentially with their length, the fitness landscape requires in principle an exponentially large number of parameters to specify, and in turn those parameters need an exponentially growing amount of data to infer.

A suitable model system for the inference of a fitness landscape is HIV proteins, due to the large number of sequences stored in clinical databases and the relative ease of generating mutants and measuring the resulting fitness. In a series of papers, Chakraborty and co-workers proposed a fitness model the so-called Gag protein family (group-specific antigen) of the HIV virus dahirel2011a ; ferguson2013a ; shekhar2013a ; mann2014a . The model is based on pairwise interactions between amino acids. Retaining only the information whether the amino acid at sequence position $i$ was mutated ($s_{i}=1$) with respect to a reference sequence or not ($s_{i}=0$), Chakraborty and co-workers suggest a minimal model for the fitness landscape given by the Ising Hamiltonian (2). Again, one can view the landscape (2) as generating the maximum entropy distribution constrained by the observed one- and two-point correlations.

Adding a constant to (2) in order to make fitness (expected number of offspring) non-negative does not alter the resulting statistics. The inverse problem is to infer the couplings $J_{ij}$ and fields $h_{i}$ from frequencies of amino acids and pairs of amino acids in particular sequence positions observed in HIV sequence data. Of course it is not clear from the outset that a model with only pairwise interactions can describe the empirical fitness landscape. As a test of this approach, ferguson2013a compares the prediction of (2) for specific mutants to the results of independent experimental measurements of fitness.

Statistical models of sequences described by pairwise interactions may be useful to model a wide range of protein families with different functions hopf2015a , and have been used in other contexts as well. Santolini, Mora, and Hakim model the statistics of sequences binding transcription factors using (7), with each spin taking one of four states to characterize the nucleotides $A,C,G,T$ santolini2014a . A similar model is used in mora2010a to model the sequence diversity of the so-called IgM protein, an antibody which plays a key role in the early immune response. The model with pairwise interactions predicts non-trivial three-point correlations which compare well with those found in the data, see figure 3.

I.5 Combinatorial antibiotic treatment

Antibiotics are chemical compounds which kill specific bacteria or inhibit their growth walsh2000a ; kohanski2010a . Mutations in the bacterial DNA can lead to resistance against a particular antibiotic, which is a major hazard to public health wise1998a ; liu2015a . One strategy to slow down or eliminate the emergence of resistance is to use a combination of antibiotics either simultaneously or in rotation walsh2000a ; kohanski2010a . The key problem of this approach is to find combinations of compounds which are particularly effective against a particular strain of bacteria. Trying out all combinations experimentally is prohibitively expensive. Wood et al. use an inverse statistical approach to predict the effect of combinations of several antibiotics from data on the effect of pairs of antibiotics wood2012a . The available antibiotics are labelled $i=1,\ldots,N$; in wood2012a a distribution over continuous variables $x_{i}$ is constructed, such that $\langle x_{i}\rangle$ gives the bacterial growth rate when antibiotic $i$ is administered, $\langle x_{i}x_{j}\rangle$ gives the growth rate with both $i$ and $j$ are given, etc. for higher moments. Choosing this distribution to be a multi-variate Gaussian $P(\{x_{i}\})=\exp\left[\sum_{i\leq j}J_{ij}x_{i}x_{j}+\sum_{i}h_{i}x_{i}\right]/Z$ results in simple relationships between the different moments, which lead to predictions of the response to drug combinations that are borne out well by experiment wood2012a .

I.6 Interactions between species and between individuals

Species exist in various ecological relationships. For instance individuals of one species hunt and eat individuals of another species. Another example is microorganisms whose growth can be influenced, both positively and negatively, by the metabolic output of other microorganisms. Such relationships form a dense web of ecological interactions between species. Co-culturing and perturbation experiments (for instance species removal) lead to data which may allow the inference of these networks faust2012a ; hekstra2013 .

Interactions between organisms exist also at the level of individuals, for instance when birds form a flock, or fish form a school. This emergent collective behaviour is thought to have evolved to minimize the exposure of individuals to predators. In bialek2012a ; bialek2014a , a model with pairwise interactions between the velocities of birds in a flock is constructed. Individual birds labelled $i=1,\ldots,N$ move with velocity $\vec{v}_{i}$ in a direction specified by $\vec{s}_{i}=\vec{v}_{i}/|\vec{v}_{i}|$. The statistics of the these directions is modelled by a distribution

[TABLE]

where the couplings between the spins $\vec{s}_{i}$ need to be inferred from the experimentally observed correlations between normalized velocities. This model can be viewed as the maximum-entropy distribution constrained by pairwise correlations between normalized velocities. From the point of view of statistical physics it describes a disordered system of Heisenberg spins. As birds frequently change their neighbours in flight, the couplings are not constant in time and it makes sense to consider couplings that depend on the distance between two individuals bialek2012a . An alternative is to apply the maximum entropy principle to entire trajectories cavagna2014a .

I.7 Financial markets

Market participants exchange commodities, shares in companies, currencies, or other goods and services, usually for money. The change in prices of such goods are often correlated, as the demand for different goods can be influenced by the same events. In bury2013a ; bury2013b , Bury uses a spin model with pairwise interactions to analyse stock market data. Shares in $N$ different companies are described by binary spin variables, where spin $s_{i}=1$ indicates ‘bullish’ conditions for shares in company $i$ with prices going up at a particular time, and $s_{i}=-1$ implies ‘bearish’ conditions with decreasing prices. Couplings $J_{ij}$ describe how price changes in shares $i$ affect changes in the price of $j$, or how prices are affected jointly by external events. Bury fit stock marked data to this spin model, and found clusters in the resulting matrix of couplings bury2013a . These clusters correspond to different industries whose companies are traded on the market. In borysov2015a , a similar analysis finds that heavy tails in the distribution of inferred couplings are linked to such clusters. Slonim et al. identified clusters in stocks using an information-based metric of stock prices slonim2005a .

II Equilibrium reconstruction

The applications discussed above can be classified according to the symmetry of pairwise couplings: In network reconstruction, couplings between spins are generally asymmetric, in maximum entropy models they are symmetric. A stochastic dynamics based on symmetric couplings entails detailed balance, leading to a steady state described by the Boltzmann distribution krapivsky2010kinetic , whereas asymmetric couplings lead to a non-equilibrium steady state. This distinction shapes the structure of this review: In this section, we discuss the inverse Ising problem in equilibrium, in section III we turn to non-equilibrium scenarios.

II.0.1 Definition of the problem

We consider the Ising model with $N$ binary spin variables $\{s_{i}=\pm 1\},i=1,\ldots,N$. Pairwise couplings (or coupling strengths) $J_{ij}$ encode pairwise interactions between the spin variables, and local magnetic fields $h_{i}$ act on individual spins. The energy of a spin configuration $\boldsymbol{s}\equiv\{s_{i}\}$ is specified by the Hamiltonian

[TABLE]

The equilibrium statistics of the Ising model is described by the Boltzmann distribution

[TABLE]

where we have subsumed the temperature into couplings and fields such that $k_{B}T=1$: The statistics of spins under the Boltzmann distribution $\exp\{-\beta H\}/Z$ depends on couplings, magnetic fields, and temperature only through the products $\beta J_{ij}$ and $\beta h_{i}$. As a result, only the products $\beta J_{ij}$ and $\beta h_{i}$ can be inferred and we set $\beta$ to $1$ without loss of generality. The energy specified by the Hamiltonian (2) or its generalisation (7) is thus a dimensionless quantity.

$Z$ denotes the partition function

[TABLE]

In such a statistical description of the Ising model, each spin is represented by a random variable. Throughout, we denote a random spin variable by $\sigma$, and a particular realisation of that random variable by $s$. This distinction will become particularly useful in the context of non-equilibrium reconstruction in section III. The expectation values of spin variables and their functions then are denoted

[TABLE]

where $Q(\boldsymbol{s})$ is some function mapping a spin configuration to a number. Examples are the equilibrium magnetizations $m_{i}\equiv\left\langle\sigma_{i}\right\rangle=\sum_{\boldsymbol{s}}p(\boldsymbol{s})s_{i}$ or the pair correlations $\chi_{ij}\equiv\left\langle\sigma_{i}\sigma_{j}\right\rangle=\sum_{\boldsymbol{s}}p(\boldsymbol{s})(s_{i}s_{j})$. In statistics, the latter observable is called the pair average. We are also interested in the connected correlation $C_{ij}=\chi_{ij}-m_{i}m_{j}$, which in statistics is known as the covariance.

The equilibrium statistics of the Ising problem (11) is fully determined by the couplings between spins and the magnetic fields acting on the spins. Collectively, couplings and magnetic fields are the parameters of the Ising problem. The forward Ising problem is to compute statistical observables such as the magnetizations and correlations under the Boltzmann distribution (11); the couplings and fields are taken as given. The inverse Ising problem works in the reverse direction: The couplings and fields are unknown and are to be determined from observations of the spins. The equilibrium inverse Ising problem is to infer these parameters from spin configurations sampled independently from the Boltzmann distribution. We denote such a data set of $M$ samples by $\mathrm{D}=\{\boldsymbol{s}^{\mu}\}$ for $\mu=1,2,\ldots,M$. (This usage of the term ‘sample’ appears to differ from how it is used in the statistical mechanics of disordered systems, where a sample often refers to a random choice of the model parameters, not spin configurations. However, it is in line with the inverse nature of the problem: From the point of view of the statistical mechanics of disordered systems in an inverse statistical problem the ‘phase space variables’ are couplings and magnetic fields to be inferred, the ‘quenched disorder’ is spin configurations sampled from the Boltzmann distribution.)

Generally, neither the values of couplings nor the graph structure formed by non-zero couplings is known. Unlike in many instances of the forward problem, the couplings often do not conform to a regular, finite-dimensional lattice; there is no sense of spatial distance between spins. Instead, the couplings might be described by a fully connected graph, with all pairs of spins coupling to each other, generally all with different values of the $J_{ij}$. Alternatively, most of the couplings might be zero, and the non-zero entries of the coupling matrix might define a structure that is (at least locally) treelike. The graph formed by the couplings might also be highly heterogeneous with few highly connected nodes with many non-zero couplings and many spins coupling only to a few other spins. These distinctions can affect how well specific inference methods perform, a point we will revisit in section II.3, which compares the quality of different methods in different situations.

II.1 Maximum likelihood

The inverse Ising problem is a problem of statistical inference mackay2003 ; bishop2006 . At the heart of many methods to reconstruct the parameters of the Ising model is the maximum likelihood framework, which we discuss here.

Suppose a set of observations $x^{1},x^{2},\ldots,x^{M}$ drawn from a statistical model $p(x^{1},x^{2},\ldots,x^{M}|\theta)$. In the case of the Ising model, each observation would be a spin configuration $\boldsymbol{s}$. While the functional form of this model may be known a priori, the parameter $\theta$ is unknown to us and needs to be inferred from the observed data. Of course, with a finite amount of data, one cannot hope to determine the parameter $\theta$ exactly. The so-called maximum likelihood estimator

[TABLE]

has a number of attractive properties Cramer1961a : In the limit of a large number of samples, $\theta^{\mathrm{ML}}$ converges in probability to the value $\theta$ being estimated. This property is termed consistency. Also for large sample sizes, there is no consistent estimator with a smaller mean-squared error. For a finite number of samples, the maximum likelihood estimator may however be biased, that is, the mean of $\theta^{\mathrm{ML}}$ over many realisations of the samples does not equal $\theta$ (although the difference vanishes with the sample size). The term likelihood refers to $p(x_{1},x_{2},\ldots,x_{M}|\theta)$ viewed as a function of the parameter $\theta$ at constant values of the data $x_{1},x_{2},\ldots,x_{M}$. The same function at constant $\theta$ gives the probability of observing the data $x_{1},x_{2},\ldots,x_{M}$.

The maximum likelihood estimator (14) can also be derived using Bayes theorem mackay2003 ; bishop2006 . In Bayesian inference, one introduces a probability distribution $p(\theta)$ over the unknown parameter $\theta$. This prior distribution describes our knowledge prior to receiving the data. Upon accounting for additional information from the data, our knowledge is described by the posterior distribution given by Bayes theorem

[TABLE]

For the case where $\theta$ is a priori uniformly distributed (describing a scenario where we have no prior knowledge of the parameter value), the posterior probability distribution of the parameter conditioned on the observations $p(\theta|x_{1},x_{2},\ldots,x_{M})$ is proportional to $p(x_{1},x_{2},\ldots,x_{M}|\theta)$222This line of argument only works if the parameter space is bounded, so a uniform prior can be defined.. Then the parameter value maximizing the probability density $p(\theta|x_{1},x_{2},\ldots,x_{M})$ is given by the maximum likelihood estimator (14). Maximizing the logarithm of the likelihood function, termed the log-likelihood function, leads to the same parameter estimate, because the logarithm is a strictly monotonic function. As the likelihood scales exponentially with the number of samples, the the log-likelihood is more convenient to use. (This is simply the convenience of not having to deal with very small numbers: the logarithm is not linked to the quenched average considered in the statistical mechanics of disordered systems; there is no average involved and the likelihood depends on both the model parameters and the data.)

We now apply the principle of maximum likelihood to the inverse Ising problem. Assuming that the configurations in the dataset were sampled independently from the Boltzmann distribution (1), the log-likelihood of the model parameters given the observed configurations $\mathrm{D}=\{\boldsymbol{s}^{\mu}\}$ is derived easily.

[TABLE]

gives the log-likelihood per sample, a quantity of order zero in $M$ since the likelihood scales exponentially with the number of samples. The sample averages of spin variables and their functions are defined by

[TABLE]

Beyond the parameters of the Ising model, the log-likelihood (II.1) depends only on the correlations between pairs of spins observed in the data $\left\langle\sigma_{i}\sigma_{j}\right\rangle^{\mathrm{D}}$ and the magnetizations $\left\langle\sigma_{i}\right\rangle^{\mathrm{D}}$. To determine the maximum-likelihood estimates of the model parameters we thus only need the pair correlations and magnetizations observed in the sample (sample averages); at least in principle, further observables are superfluous. In the language of statistics, these sets of sample averages provide sufficient statistics to determine the model parameters.

The log-likelihood (II.1) has a physical interpretation: The first two terms are the sample average of the (negative of the) energy, the second term adds the free energy. Thus the log-likelihood is the (negative of the) entropy of the Ising system, based on the sample estimate of the energy. We will further discuss this connection in section II.1.5.

A second interpretation of the log-likelihood is based on the difference between the Boltzmann distribution (11) and the empirical distribution of the data in the sample $\mathrm{D}$, denoted $p^{\mathrm{D}}(\boldsymbol{s})\equiv\frac{1}{M}\sum_{\mu}\delta_{\boldsymbol{s}^{\mu},\boldsymbol{s}}$. The difference between two probability distributions $p(\boldsymbol{s})$ and $q(\boldsymbol{s})$ can be quantified by the Kullback–Leibler (KL) divergence

[TABLE]

which is non-negative and reaches zero only when the two distributions are identical Cover2006a . The KL divergence between the empirical distribution and the Boltzmann distribution is

[TABLE]

The second term (the negative empirical entropy) is independent of the model parameters; the best match between the Boltzmann distribution and the empirical distribution (minimal KL divergence) is thus achieved when the likelihood (II.1) is maximal.

Above, we derived the principle of maximum likelihood (14) from Bayes theorem under the assumption that the model parameter $\theta$ is sampled from a uniform prior distribution. Suppose we had the prior information that the parameter $\theta$ was taken from some non-uniform distribution, the posterior distribution would then acquire an additional dependence on the parameter. In the case of the inverse Ising problem, prior information might for example describe the sparsity of the coupling matrix, with a suitable prior $p_{\boldsymbol{J}}\sim\exp\left[-\gamma\sum_{i<j}|J_{ij}|\right]$ that assigns small probabilities to large entries in the coupling matrix. The resulting (log) posterior is

[TABLE]

up to terms that do not depend on the model parameters. The maximum of the posterior is now no longer achieved by maximizing the likelihood, but involves a second term that penalizes coupling matrices with large entries. Maximizing the posterior with respect to the parameters no longer makes the Boltzmann distribution as similar to the empirical distribution as possible, but strikes a balance between making these distributions similar while avoiding large values of the couplings. In the context of inference, the second term is called a regularisation term. Different regularisation terms have been used, including the absolute-value term in (20) as well as a penalty on the square values of couplings $\sum_{i<j}J_{ij}^{2}$ (called $\ell_{1}$- and $\ell_{2}$-regularisers, respectively). One standard way to determine the value of the regularisation coefficient $\gamma$ is to cross-validate with a part of the data that is initially withheld, that is to probe (as a function of $\gamma$) how well the model can predict aspects of the data not yet used to infer the model parameters Hastie2009a .

II.1.1 Exact maximization of the likelihood

The maximum likelihood estimate of couplings and magnetic fields

[TABLE]

has a simple interpretation. Since $\ln Z(\boldsymbol{J},\boldsymbol{h})$ serves as a generating function for expectation values under the Boltzmann distribution, we have

[TABLE]

At the maximum of the log-likelihood these derivatives are zero; the maximum-likelihood estimate of the parameters is reached when the expectation values of pair correlations and magnetizations under the Boltzmann statistics match their sample averages

[TABLE]

The log-likelihood (II.1) turns out to be a concave function of the model parameters, see II.1.2. Thus, in principle, it can be maximized by a convex optimization algorithm. One particular way to reach the maximum of the likelihood is a gradient-descent algorithm called Boltzmann machine learning Ackley1985a . At each step of the algorithm fields and couplings are updated according to

[TABLE]

The parameter $\eta$ is the learning rate of the algorithm, which has (23) as its fixed point.

In order to calculate the expectation values $\left\langle\sigma_{i}\right\rangle$ and $\left\langle\sigma_{i}\sigma_{j}\right\rangle$ on the left-hand side of these equations, one needs to perform thermal averages of the form (13) over all $2^{N}$ configurations, which is generally infeasible for all but the smallest system sizes. Analogously, when maximizing the log-likelihood (II.1) directly, the partition function is the sum over $2^{N}$ terms. Moreover, the expectation values or the partition function need to be evaluated many times during an iterative search for the solution of (23) or the maximum of the likelihood. As a result, also numerical sampling techniques such as Monte Carlo sampling are cumbersome, but have been used for moderate system sizes broderick2007faster . Habeck proposes a Monte Carlo sampler that draws model parameters from the posterior distribution habeck2014 . A recent algorithm uses information contained in the shape of the likelihood maximum to speed up the convergence ferrari2015a . An important development in machine learning has led to the so-called restricted Boltzmann machines, where couplings form a symmetric and bipartite graph. Variables fall into two classes, termed ‘visible’ and ‘hidden’, with couplings never linking variables of the same class. This allows fast learning algorithms fischer2012introduction at the expense of additional hidden variables.

We stress that the difficulty of maximizing the likelihood is associated with the restriction of our input to the first two moments (magnetisations and correlations) of the data. On the one hand, this restriction is natural, as the likelihood only depends on these two moments. On the other hand, computationally efficient methods have been developed that effectively use correlations in the data beyond the first two moments. An important example is pseudolikelihood, which we will discuss in section II.2. Other learning techniques that sidestep the computation of the partition function include score matching hyvarinen2005estimation and minimum probability flow sohl2011new . Also, when the number of samples is small (compared to the number of spins), the likelihood need no longer be the best quantity to optimize.

II.1.2 Uniqueness of the solution

We will show that the log-likelihood (II.1) is a strictly concave function of the model parameters (couplings and magnetic fields). As the space of parameters is convex, the maximum of the log-likelihood is unique.

We use the shorthands $\boldsymbol{\lambda}=\{\boldsymbol{J},\boldsymbol{h}\}$ and $\{Q_{k}(\boldsymbol{s})\}=\{s_{i}s_{j},s_{i}\}$ for the model parameters and the functions coupling to them and write the Boltzmann distribution as

[TABLE]

For such a general class of exponential distributions Hastie2009a , the second derivatives of the log-likelihood $L_{\mathrm{D}}$ with respect to a parameters obey

[TABLE]

This matrix of second derivatives is non-negative (has no negative eigenvalues) since $\sum_{ij}(\left\langle Q_{i}Q_{j}\right\rangle-\left\langle Q_{i}\right\rangle\left\langle Q_{j}\right\rangle)x_{i}x_{j}=\left\langle\left[\sum_{k}\left(x_{k}Q_{k}-\left\langle x_{k}Q_{k}\right\rangle\right)\right]^{2}\right\rangle\geq 0$ for all $x_{i}$. If no non-trivial linear combination of the observables $Q_{k}$ has vanishing fluctuations, the Hessian matrix is even positive-definite. For the inverse Ising problem, there are indeed no non-trivial linear combinations of the spin variables $\sigma_{i}$ and pairs of spins variables $\sigma_{i}\sigma_{j}$ that do not fluctuate under the Boltzmann measure, unless some of the couplings or fields are infinite. As a result, the maximum of the likelihood, if it exists, is unique. However, it can happen that the maximum lies at infinite values of some of the parameters (for instance when the samples contain only positive values of a particular spin, the maximum likelihood value of the corresponding magnetic field is infinite). These divergences can be avoided with the introduction of a regularisation term, see section II.1.

II.1.3 Maximum entropy modelling

The Boltzmann distribution in general and the Ising Hamiltonian (2) in particular can be derived from information theory and the principle of maximum entropy. This principle has been invoked in neural modelling Schneidman2006a , protein structure determination Weigt2009a , and DNA sequence analysis mora2010a . In this section, we discuss the statistical basis of Shannon’s entropy, the principle of maximum entropy, and their application to inverse statistical modelling.

Consider $M$ distinguishable balls, each to be placed in a box with $R$ compartments. The number of ways of placing the balls such that $n_{r}$ balls are in the $r$th compartment ($r\in 1,\ldots,R$) is

[TABLE]

with $\sum_{r=1}^{R}n_{r}=M$. For large $M$, we write $n_{r}=Mq_{r}$ and exploit Stirling’s formula $n_{r}!\approx e^{n_{r}}n_{r}^{n_{r}}$, yielding the Gibbs entropy

[TABLE]

This combinatorial result forms the basis of equilibrium statistical physics in the classic treatment due to Gibbs and can be found in standard textbooks. In the context of statistical physics, each of the $R$ compartments corresponds to a microstate of a system, and each microstate $r$ is associated with energy $E_{r}$. The $M$ balls in the compartments describe a set of copies of the system, or a so-called ensemble of replicas. The replicas may exchange energy with each other, while the ensemble of replicas itself is isolated and has a fixed total energy $ME$ (and possibly other conserved quantities). In this way, the replicas can be thought of as providing a heat-bath for each other. If we assume that each state of the ensemble of replicas with a given total energy is equally likely, the statistics of $q_{r}$ is dominated by a sharp maximum of $W$ as a function of the $q_{r}$, subject to the constraint $\sum_{r}q_{r}=1$ and $\sum_{r}E_{r}q_{r}=E$. Using Lagrange multipliers to maximize (29) subject to these constraints yields the Boltzmann distribution jaynes1957a .

This seminal line of argument can also be used to derive Shannon’s information entropy (3). The argument is due to Wallis and is recounted in Jaynes1989a . Suppose we want to find a probability distribution $p_{r}$ compatible with a certain constraint, for instance a specific expectation value $\sum_{r}p_{r}E_{r}=E$ to within some small margin of error. Consider $M$ independent casts of a fair die with $R$ faces. We denote the number of times outcome $r$ is realized in these throws as $n_{r}$. The probability of a particular set $\{n_{r}\}$ is

[TABLE]

In the limit of large $M$, the logarithm of this probability is $-M\sum_{r}q_{r}\ln q_{r}-M\ln R$ with $q_{r}=n_{r}/M$.

Each set of $M$ casts defines one instance of the $\{n_{r}\}$. In most instances, the constraint will not be realized. For those (potentially rare) instances obeying the constraint, we can ask what are the most likely values of $n_{r}$, and correspondingly $q_{r}$. Maximising the Shannon’s information entropy $-\sum_{r}q_{r}\ln q_{r}$ subject to the constraint and the normalisation $\sum_{r}q_{r}=1$ gives the so-called maximum-entropy estimate of $p_{r}$. If the underlying set of probabilities (the die with $R$ faces) differs from the uniform distribution, so outcome $r$ occurs with probability $q^{0}_{r}$, it is not the entropy but the relative entropy $-\sum_{r}q_{r}\ln\frac{q_{r}}{q^{0}_{r}}$ that is to be maximised. Up to a sign, this is the Kullback-Leibler divergence (18) between $q_{r}$ and $q_{r}^{0}$.

The maximum-entropy estimate can be used to approximate an unknown probability distribution $q_{r}$ that is under-sampled. Suppose data is sampled several times from some unknown probability distribution. With a sufficient number of samples $M$, the distribution $q_{r}$ can be easily determined from frequency counts $q_{r}=n_{r}/M$. Often this is not feasible; if $q_{r}M\ll 1$, $n_{r}$ fluctuates strongly from one set of samples to the next. This situation appears naturally when the number of possible outcomes $R$ grows exponentially with the size of the system, see e.g. section I.1. Nevertheless the data may be sufficient to pin down one or several expectation values. The maximum-entropy estimate has been proposed as the most unbiased estimate of the unknown probability distribution compatible with the observed expectation values Jaynes1989a . For a discussion of the different ways to justify the maximum entropy principle, and derivations based on robust estimates see tikochinsky1984alternative .

Many applications of the maximum- entropy estimate are in image analysis and spectral analysis gull1984a , for reviews in physics and biology see banavar2010a ; bialek2007rediscovering ; Mora2011b , and for critical discussion see aurell2016maximum ; vannimwegen2016inferring .

The connection between maximum entropy and the inverse Ising problem is simple: For a set of $N$ binary variables, the distribution with given first and second moments maximizing the information entropy is the Boltzmann distribution (1) with the Ising Hamiltonian (2). We use Lagrange multipliers to maximize the information entropy (3) subject to the normalization condition and the constraints on the first and second moments (magnetisations and pair correlations) of $p(\boldsymbol{s})$ to be $\boldsymbol{m}$ and $\boldsymbol{\chi}$. Setting the derivatives of

[TABLE]

with respect to $p(\boldsymbol{s})$ to zero yields the Ising model (1). The Lagrange multipliers $\boldsymbol{h}$ and $\boldsymbol{J}$ need to be chosen to reproduce the first and second moments (magnetisations and correlations) of the data and can be interpreted as couplings between spins and magnetic fields.

While this principle appears to provide a statistical foundation to the model (1), there is no a priori reason to disregard empirical data beyond the first two moments. Instead, the pairwise couplings result from the particular choice of making the probability distribution match the first two moments of the data. The reasons for this step may be different in different applications.

• •

Moments beyond the first and second may be poorly determined by the data. Conversely, with an increasing number of samples, the determination of higher order correlations and hence interactions between triplets of spin variables etc. becomes viable.

• •

The data may actually be generated by an equilibrium model with (at most) pairwise interactions between spin variables. This need not be obvious from observed correlations of any order, but can be tested by comparing three-point correlations predicted by a model with pairwise couplings to the corresponding correlations in the data. Examples are found in sequence analysis, where population dynamics leads to an equilibrium steady state BergWillmannLaessig2004 ; SellaHirsh2005 and the energy can often be approximated by pairwise couplings mora2010a ; santolini2014a . For a review see stein2015a . Surprisingly, also in neural data (not generated by an equilibrium model), three-point correlations are predicted well by a model with pairwise interactions tkacik2009spin ; tkavcik2014a .

• •

A model of binary variables interacting via a high-order coupling terms $J_{ijk\ldots}s_{i}s_{j}s_{k}\ldots$ can sometimes be approximated surprisingly well by pairwise interactions. This seems to be the case when the couplings are dense, so that each variable appears in several coupling terms merchan2015a .

• •

Often one seeks to describe a subset of $n$ variables $s_{1},s_{2},\ldots,s_{n}$ from a larger set of $N$ variables, for instance when only the variables in the subset can be observed. The subset of variables is characterized by effective interactions which stem from interactions between variables in the subset, and from interactions with the other variables. If the subset is sufficiently small, the resulting statistics is often described by a model with pairwise couplings roudinirenberglatham2009a .

• •

The true probability distribution underlying some data may be too complicated to calculate in practice. A more modest goal then is to describe the data using an effective statistical model such as (1), which is tractable and allows the derivation of bounds on the entropy or the free energy. Examples are the description of neural data and gene expression data using the Ising model with symmetric couplings (see I.1 and I.2).

• •

There are also useful models which are computationally tractable but do not maximize the entropy. An example is Gaussian models to generate artificial spike trains with prescribed pair correlations amari2003synchronous ; macke2009generating .

II.1.4 Information theoretic bounds on graphical model reconstruction

A particular facet of the inverse Ising problem is graphical model selection. Consider the Ising problem on a graph. A graph is a set of nodes and edges connecting these nodes, and nodes associated with spin variables. Couplings between node pairs connected by an edge are non-zero, couplings between unconnected node pairs are zero. The graphical model selection problem is to recover the underlying graph (and usually also the values of the couplings) from data sampled independently from the Boltzmann distribution. Given a particular number of samples, one can ask with which probability a given method can reconstruct the graph correctly (the reconstruction fluctuates between different realisations of the samples). Notably, there are also universal limits on graphical model selection that are independent of a particular method.

In santhanam2012information , Santhanam and Wainwright derive information-theoretic limits to graphical model selection. Key result is the dependence of the required number of samples on the smallest and on the largest coupling

[TABLE]

and on the maximum node connectivity (number of neighbours on the graph) $d$. Reconstruction of the graph, by any method, is impossible if fewer than

[TABLE]

samples are available (the precise statement is of a probabilistic nature, see santhanam2012information ). If the maximum connectivity $d$ grows with the system size, this result implies that at least $c\max\{d^{2},\alpha^{-2}\}\ln N$ samples are required (with some constant $c$) santhanam2012information . The derivation of this and other results is based on Fano’s inequality (Fano’s lemma) Cover2006a , which gives a lower bound for the probability of error of a classification function (such as the mapping from samples to the graph underlying these samples).

II.1.5 Thermodynamics of the inverse Ising problem

Calculations in statistical physics are greatly simplified by introducing thermodynamic potentials. In this section, we will discuss the method of thermodynamic potentials for the inverse Ising problem. It turns out that the maximum likelihood estimation of the fields and couplings is simply a transformation of the thermodynamic potentials.

Recall that the thermodynamic potential most useful for the forward problem, where couplings and magnetic fields are given, is the Helmholtz free energy $F(\boldsymbol{J},\boldsymbol{h})=-\ln Z(\boldsymbol{J},\boldsymbol{h})$. Derivatives of this free energy give the magnetizations, correlations, and other observables. The thermodynamic potential most useful for the inverse problem, where the pair correlations $\boldsymbol{\chi}$ and magnetizations $\boldsymbol{m}$ are given, is the Legendre transform of the Helmholtz free energy with respect to both couplings and fields Sessak2009a ; Cocco2011a ; Cocco2012a

[TABLE]

This thermodynamic potential is readily recognised as the entropy function; up to a sign, it gives the maximum likelihood (II.1) of the model parameters. The transformation (34) thus provides a link between the inference via maximum likelihood and the statistical physics of the Ising model as described by its Helmholtz free energy. The couplings and the fields are found by differentiation,

[TABLE]

where the derivatives are evaluated at the sample correlations and magnetizations. These relationships follow from the inverse transformation of (34)

[TABLE]

by setting derivatives of the term in square brackets with respect to $\boldsymbol{\chi}$ and $\boldsymbol{m}$ to zero.

In practice, performing the Legendre transformation of both $\boldsymbol{h}$ and $\boldsymbol{J}$ is often not necessary; derivatives of the Helmholtz free energy $F(\boldsymbol{J},\boldsymbol{h})$ with respect to $\boldsymbol{J}$ can also be generated by differentiating with respect to $\boldsymbol{h}$, e.g.,

[TABLE]

The thermodynamics of the inverse problem can thus be reduced to a single Legendre transform of the Helmholtz free energy, yielding the Gibbs free energy

[TABLE]

The magnetic fields are given by the first derivative of the Gibbs free energy

[TABLE]

To infer the couplings, we consider the second derivatives of Gibbs potential, which give

[TABLE]

where $\boldsymbol{C}$ is the matrix of connected correlations $C_{ij}\equiv\chi_{ij}-m_{i}m_{j}$. (40) follows from the inverse function theorem,

[TABLE]

and linear response theory

[TABLE]

which links the susceptibility of the magnetization to a small change in the magnetic field with the connected correlation stanley1987introduction .

The result (40) turns out to be central to many methods for the inverse Ising problem. The left-hand side of this expression is a function of $J_{ij}$. If the Gibbs free energy (38) can be evaluated or approximated, (40) can be solved to yield the couplings. Similarly (39) with the estimated couplings and the sample magnetisations gives the magnetic fields, completing the reconstruction of the parameters of the Ising model.

II.1.6 Variational principles

For most systems, neither the free energy $F(\boldsymbol{J},\boldsymbol{h})$ nor other thermodynamic potentials can be evaluated. However, there are many approximation schemes for $F(\boldsymbol{J},\boldsymbol{h})$ AMF , which lead to approximations for the entropy and the Gibbs free energy. Direct approximation schemes for $S(\boldsymbol{\chi},\boldsymbol{m})$ and $G(\boldsymbol{J},\boldsymbol{m})$ have also be formulated within the context of the inverse Ising problem. The key idea behind most of these approximations is the variational principle.

The variational principle for the free energy is

[TABLE]

where $q$ denotes a probability distribution over spin configurations. $U[q]\equiv\left\langle H\right\rangle_{q}$ and $S[q]\equiv\left\langle\ln q\right\rangle_{q}$ and the minimisation is taken over all distributions $q$. This principle finds its origin in information theory. Take an arbitrary trial distribution $q(\boldsymbol{s})$, the Kullback-Leibler divergence (18) quantifies the difference between $q$ and the Boltzmann distribution $p$ is positive and vanishes if and only if $q=p$ Cover2006a . One then arrives directly at (43) when rewriting $\operatorname{KL}(q|p)=U[q]-S[q]-F(\boldsymbol{J},\boldsymbol{h})$.

We will refer to $F[q]\equiv U[q]-S[q]$ as the functional Helmholtz free energy, also called the non-equilibrium free energy in the context of non-equilibrium statistical physics parrondo2015a . Another term in use is ‘Gibbs free energy’ AMF , which we have reserved for the thermodynamic potential (38).

So far, nothing has been gained as the minimum is over all possible distributions $q$, including the Boltzmann distribution itself. A practical approximation arises when a constraint is put on $q$, leading to a family of trial distributions $q$. Often the minimisation can then be carried out over that family, yielding an upper bound to the Helmholtz free energy AMF .

In the context of the inverse problem, it is useful to derive the variational principles for other thermodynamic potentials as well. Using the definition of Gibbs potential (38) and the variational principle for the Helmholtz potential (43) we obtain

[TABLE]

By means of Lagrange multipliers it is easy to show that this double extremum can be obtained by a single conditional minimisation,

[TABLE]

where the set $\mathcal{G}$ denotes all distributions $q$ with a given $\left\langle\sigma_{i}\right\rangle_{q}=m_{i}$ AMF . We will refer to the functional $G[q]=\left\{-\sum_{i<j}J_{ij}\left\langle\sigma_{i}\sigma_{j}\right\rangle_{q}-S[q]\right\}$ as the functional Gibbs free energy defined on $\mathcal{G}$.

Similarly, the variational principle can be applied to the entropy function $S(\boldsymbol{\chi},\boldsymbol{m})$, leading once again to a close relationship between statistical modelling and thermodynamics. The entropy (34) is found to be

[TABLE]

where $\mathcal{S}$ denotes distributions with $\left\langle\sigma_{i}\right\rangle_{q}=m_{i}$ and $\left\langle\sigma_{i}\sigma_{j}\right\rangle_{q}=\chi_{ij}$. This is nothing but the maximum entropy principle Jaynes1989a : the variational principle identifies the distribution with the maximum information entropy subject to the constraints on magnetisations and spin correlations, which are set equal to their sample averages (see the section on maximum entropy modelling above).

II.1.7 Mean-field theory

As a first demonstration of the variational principle, we derive the mean-field theory for the inverse Ising problem. The starting point is an ansatz for the Boltzmann distribution (11) which factorises in the sites stanley1987introduction ; AMF ; tanaka2000a

[TABLE]

thus making the different spin variables statistically independent of one another. The parameters $\tilde{m}_{i}$ of this ansatz describe the spin magnetizations; each spin has a magnetisation resulting from the effective magnetic field acting on that spin. This effective field arises from its local magnetic field $h_{i}$, as well as from couplings with other spin. The mean field giving its name to mean-field theory is the average over typical configurations of the effective field.

Using the mean-field ansatz, we now estimate the Gibbs free energy. Within the mean-field ansatz, the minimisation of the variational Gibbs potential (45) is trivial: there is only a single mean-field distribution (47) that satisfies the constraint $\mathcal{G}$ that spins have magnetisations $\boldsymbol{m}$, namely $\tilde{\boldsymbol{m}}=\boldsymbol{m}$. We can thus directly write down the mean-field Gibbs free energy

[TABLE]

The equation for the couplings $\boldsymbol{J}$ follows from the second order derivative of $G(\boldsymbol{m},\boldsymbol{J})$, cf. equation (40)

[TABLE]

Similarly, the reconstruction of the magnetic field follows from the derivative of $G(\boldsymbol{m},\boldsymbol{J})$ with respect to $m_{i}$, cf. equation (39),

[TABLE]

This result establishes a simple relationship between the observed connected correlations and the couplings between spins in terms of the inverse of the correlation matrix. The matrix inverse in (49) is of course much simpler to compute than the maximum over the likelihood (II.1) and takes only a polynomial number of steps: Gauß-Jordan elimination for inverting an $N\times N$ matrix requires ${\mathcal{O}}(N^{3})$ operations, compared to the exponentially large number of steps to compute a partition function or its derivatives.

The standard route to mean-field reconstruction proceeds somewhat differently, namely via the Helmholtz free energy rather than the Gibbs free energy. Early work to address the inverse Ising problem using the mean-field approximation was performed by Peterson and Anderson peterson1987mean . In Kappen1998b , Kappen and Rodríguez construct the Helmholtz functional free energy $F^{\mathrm{MF}}_{\tilde{\boldsymbol{m}}}(\boldsymbol{J},\boldsymbol{h})$ given by (43) under the mean-field ansatz (47). $F^{\mathrm{MF}}_{\tilde{\boldsymbol{m}}}(\boldsymbol{J},\boldsymbol{h})$ is then minimised with respect to the parameters of the mean-field ansatz $\tilde{\boldsymbol{m}}$ by setting its derivatives with respect to $\tilde{\boldsymbol{m}}$ to zero. This yields equations which determine the values $m_{i}^{\mathrm{MF}}$ of the magnetization parameters that minimize the KL divergence between the mean-field ansatz and the Boltzmann distribution, namely the well-known self-consistent equations

[TABLE]

Using $m^{\mathrm{MF}}_{i}$ as an approximation for the equilibrium magnetisations $m_{i}$ one can derive the so-called linear-response approximation for the connected correlation function $C^{\mathrm{MF}-\mathrm{LR}}_{ij}\equiv\frac{\partial m^{\mathrm{MF}}_{i}}{\partial h_{j}}(\boldsymbol{h})$. Taking derivatives of the self-consistent equations (51) with respect to local fields gives

[TABLE]

where we have used the fact that diagonal terms of the coupling matrix are zero. Identifying the result for the connected correlations $C^{\mathrm{MF}-\mathrm{LR}}_{jk}$ with the sample correlations $C_{jk}$ leads to a system of linear equations to be solved for the couplings Kappen1998b . However, to obtain (52), we have used that the diagonal elements $J_{ii}$ are zero. With these constraints, the system of equations (52) becomes over-determined and in general there is no solution. Although different procedures have been suggested to overcome this problem Kappen1998b ; Ricci2012a ; JaquinRancon2016 , there seems to be no canonical way out of this dilemma. The most common approach is to ignore the constraints on the diagonal elements altogether and invert equation (52) to get

[TABLE]

This result agrees with the reconstruction via the Gibbs free energy except for the non-zero diagonal couplings, which bear no physical meaning and are to be ignored. No diagonal couplings arise in the approach based on the Gibbs free energy since equation (40) with $j=i$ does not involve any unknown couplings $J_{ij}$.

II.1.8 The Onsager term and TAP reconstruction

The variational estimate of the Gibbs free energy (48) can be improved further. In 1977, Thouless, Anderson, and Palmer (TAP) advocated adding a term to the Gibbs free energy

[TABLE]

This term can be interpreted as describing the effect of fluctuations of a spin variable on the magnetisation of that spin via their impact on neighbouring spins Thouless1977a . It is called the Onsager term, which we will derive in section II.1.13 in the context of a systematic expansion around the mean-field ansatz. For the forward problem, adding this term modifies the self-consistent equation (51) to the so-called TAP equation

[TABLE]

In the inverse problem, the TAP free energy(54) gives an the equation for the couplings based on (40)

[TABLE]

Solving this quadratic equation gives the TAP reconstruction Kappen1998b ; tanaka1998mean

[TABLE]

where we have chosen the solution that coincides with the mean-field reconstruction when the magnetisations are zero. The magnetic fields can again be found by differentiating the Gibbs free energy

[TABLE]

II.1.9 Couplings without a loop: mapping to the minimum spanning tree

problem

The computational hardness of implementing Boltzmann machine learning (24) comes from the difficulty of computing correlations under the Boltzmann measure, which can require a computational time that scales exponentially with the system size. This scaling originates from the presence of loops in the graph of couplings between the spins. Graphs for which correlations can be computed efficiently are the acyclic graphs or trees, so it comes as no surprise that the first efficient method to solve the inverse Ising problem was developed for trees already in 1968. This was done by Chow and Liu Chow1968a in the context of a product approximation to a multi-variate probability distribution. While the method itself can be used as a crude approximation for models with loops or as reference point for more advanced methods, the exact result by Chow and Liu is of basic interest in itself as it provides a mapping of the inverse Ising problem for couplings forming a tree onto a minimum spanning tree (MST) problem. MST is a core problem in computational complexity theory, for which there are many efficient algorithms. This section on Chow and Liu’s result also provides some of the background needed in section II.1.10 on the Bethe ansatz and section II.1.11 on belief propagation.

We consider an Ising model whose pairwise couplings form a tree. The graph of couplings may consist of several parts that are not connected to each other (in any number of steps along connected node pairs), or it may form one single connected tree, but it contains no loops. We denote the set of nodes (vertices) associated with a spin of the tree $T$ by $V_{T}$ and the set of edges (couplings between nodes) by $E_{T}$. It is straightforward to show that in this case, the Boltzmann distribution for the Ising model can be written in a pairwise factorised form

[TABLE]

$\partial i$ denotes the set of neighbours of node $i$, so $|\partial i|$ is the number of nodes $i$ couples to. The distributions $p_{i}$ and $p_{ij}$ denote the one-point and two-point marginals of $p_{T}$.

The KL divergence (II.1) between the empirical distribution $p^{\mathrm{D}}(\boldsymbol{s})$ and $p_{T}$ is given by

[TABLE]

For a given tree, it is straightforward to show that the KL divergence is minimized when the marginals $p_{i}$ and $p_{ij}$ equal the empirical marginals $p^{\mathrm{D}}_{i}$ and $p^{\mathrm{D}}_{ij}$. This gives

[TABLE]

where $H=-\sum_{\{\boldsymbol{s}\}}p^{\mathrm{D}}(\boldsymbol{s})\ln p^{\mathrm{D}}(\boldsymbol{s})$ is the entropy of the empirical distribution $p^{\mathrm{D}}$, $H_{i}=-\sum_{s_{i}}p^{\mathrm{D}}_{i}(s_{i})\ln p^{\mathrm{D}}_{i}(s_{i})$ is the single site entropy and $I_{ij}$ is the mutual information between a pair of spins

[TABLE]

Assuming the graph of couplings is an (unknown) tree and that empirical estimates of all pairs of mutual informations are available, the inverse Ising problem can then be solved by minimizing the KL divergence (60) over all possible $N^{N-2}$ trees $T$. The first two terms in eq. (61) do not depend on the choice of the tree, so only the last term needs to be minimized over, which is a sum over local terms on the graph. The optimal tree topology $T_{\text{opt}}$ is thus given by

[TABLE]

This is where the mapping onto MST problem emerges: $T_{\text{opt}}$ connects all vertices of the graph and its edges are such that the total sum of their weights is minimal. In our case, each edge weight is the (negative) pairwise mutual information $I_{ij}$ between spin variables. Finding the MST does not require an infeasible exploration of the space of all possible trees. On the contrary, it can be found in a number of steps bounded by $\mathcal{O}(|V|^{2}\ln|V|)$ by greedy iterative procedures which identify the optimal edges to be added at each step ($V$ is the set of nodes in the data). The most famous algorithms for the MST problem date back to the 1950s and are known under the names of Prim’s algorithm, Kruskal’s algorithms and Boruvka’s algorithm (see e.g. Moore2011a ). In practice, one has to compute the empirical estimates of the mutual information from samples and then proceed with one of the above algorithms. An interesting observation which makes the Chow–Liu approach even easier to apply is that one may use as edge weights also the connected correlations between spins Chow1968a .

Once the optimal tree $T_{\text{opt}}$ has been identified, we still need to find the optimal values of the couplings $J_{ij}$ of the Ising model. This is, however, an easy task: the factorised form of the probability measure over the tree (59) allows one to compute the couplings using the independent-pair approximation, see subsection II.1.12.

II.1.10 The Bethe–Peierls ansatz

The factorising probability distribution (59) can also be used as an ansatz in situations where the graph of couplings is not a tree. In this context,

[TABLE]

is called the Bethe–Peierls ansatz bethe1935 ; peierls1936 . $(i,j)\in E$ runs over pairs of interacting spins, or equivalently over edges in the graph of couplings. One can parameterize the marginal distribution $p_{i}$ and $p_{ij}$ using the magnetisation parameters $\tilde{m}_{i}$ and the (connected) correlation parameters $\tilde{C}_{ij}$,

[TABLE]

subject to the constraints

[TABLE]

The Bethe–Peierls ansatz can be compared to the mean-field ansatz (47), which assigns a magnetisation (or an effective field) to each spin. The Bethe–Peierls ansatz goes one step further; it assigns to each coupled pair of spins correlations as well as magnetisations. These correlations and magnetisations are then determined self-consistently. An important feature of the Bethe–Peierls ansatz (63) is that its Shannon entropy (3) can be decomposed into spin pairs

[TABLE]

For graphs containing loops, the entropy generally contains terms involving larger sets of spins than pairs, a situation we will discuss in section II.1.12.

The Bethe–Peierls ansatz is well defined, and indeed exact, when the couplings form a tree. When the graph of couplings contains loops, the probability distribution (63) is not normalized and the ansatz is not well defined. In that case, the Bethe–Peierls ansatz is an uncontrolled approximation, although recently progress has been made in the control of the resulting error Chertkov2006a ; Parisi2006a . We start with the assumption that there are no loops.

To address the inverse Ising problem, we use the Bethe–Peierls ansatz (63) as a variational ansatz to minimise the functional Gibbs free energy (45). Again the constraint $\mathcal{G}$ in (45) implies $\tilde{\boldsymbol{m}}=\boldsymbol{m}$. The remaining minimisation is over the correlation parameters $\tilde{\boldsymbol{C}}$,

[TABLE]

and yields

[TABLE]

where $\boldsymbol{C}^{\mathrm{BP}}$ ist the optimal value of $\tilde{\boldsymbol{C}}$ and satisfies

[TABLE]

$z_{i}$ denotes the number of neighbours (interaction partners with non-zero couplings) of node $i$. From (40), the equation for the couplings can be found again by equating the second derivative of the Gibbs free energy with the inverse of the correlation matrix, that is,

[TABLE]

This quadratic equation can be solved for $C^{\mathrm{BP}}_{ij}$; inserting the solution

[TABLE]

for $(\boldsymbol{C}^{-1})_{ij}\neq 0$ and $C^{\mathrm{BP}}_{ij}=0$ for $(\boldsymbol{C}^{-1})_{ij}=0$ into (69) gives the couplings of the Bethe–Peierls reconstruction Nguyen2012a ; Ricci2012a . For the special case $m_{i}=0$, one obtains a particularly simple result

[TABLE]

In graph theory, this formula can be related to the expression for the distance in a tree whose links carry weights specified by pair correlations between spin pairs Bapat2005a . Correspondingly, the magnetic fields follow from the first derivative of the Gibbs free energy as in (39) giving

[TABLE]

In section II.1.13, we will show that the Bethe–Peierls ansatz, and hence the resulting reconstruction, is exact for couplings forming a tree. However, the reconstruction of couplings and magnetic fields based on the Bethe–Peierls ansatz can also be applied to cases where the couplings do not form a tree. Although the results then arise from an uncontrolled approximation, the quality of the reconstruction can still be rather good. For a comparison of the different approaches, see section II.3.

II.1.11 Belief propagation and susceptibility propagation

Belief propagation is a distributed algorithm to compute marginal distributions of statistical models, such as $p_{i}(s_{i})$ and $p_{ij}(s_{i},s_{j})$ of the preceding sections. Again, it is exact on trees, but also gives a good approximation when the graph of couplings is only locally treelike (so any loops present are long). The term belief propagation is used in the machine learning pearl1988probabilistic and statistical physics communities, in coding theory the approach is known as the sum-product algorithm gallager1962 . Belief propagation shares a deep conceptual link with the Bethe–Peierls ansatz; Yedidia et al. Yedidia2005a showed that belief propagation is a numerical scheme to efficiently compute the parameters of Bethe–Peierls ansatz.

A detailed exposition and many applications of belief propagation can be found in the textbook by Mézard and Montanari mezard2009information . Here, we briefly introduce the basics and discuss applications to the inverse Ising problem. We start by considering the ferromagnetic Ising model in 1D, that is, a linear chain of $N$ spins. The textbook solution of this problem considers the partitions function of the system when the last spin is constrained to take on the values $s_{N}=\pm 1$, $Z_{N}(+1)$ and $Z_{N}(-1)$. The corresponding partition functions for the linear chain with $N+1$ spins are linked to the former via the so-called transfer matrix baxter1982exactly . The partition function for a system of any size can be computed iteratively starting from a single spin and extending the 1D lattice with each multiplication of the transfer matrix. In fact the transfer matrix can also be used to solve the Ising model on a tree Eggarter1974a . Belief propagation is similar in spirit, but can be extended (as an approximation) also to graphs which are not trees. The best book-keeping device is again a restricted partition function, namely $Z_{i\to j}(s_{i})$. It is defined as the partition function of the part of the system containing $i$ when the coupling present between spins $i$ and $j$ has been deleted from the tree and spin $i$ is constrained to take on $s_{i}$. (Deleting the edge $(i,j)$ splits the tree containing $i$ and $j$ into two disconnected parts.) On a tree we obtain the recursion relation relation

[TABLE]

which can be computed recursively starting from leaves of the tree (nodes connected to a single edge only). In statistical physics, (74) is called the cavity recursion for partition functions, since deleting a link can be thought of as leaving a ‘cavity’ in the original system. The partition function for the entire tree with spin $i$ constrained to $s_{i}$ is then

[TABLE]

The marginal distribution $p_{i}(s_{i})$ can be calculated by normalizing $p_{i}(s_{i})\propto Z_{i}(s_{i})$ and the marginal distribution $p_{ij}(s_{i},s_{j})$ can be calculated by normalizing $p_{ij}(s_{i},s_{j})\propto Z_{i\to j}(s_{i})e^{J_{ij}s_{i}s_{j}}Z_{j\to i}(s_{j})$.

The power of the cavity recursion (76) lies in its application to graphs which are not trees. It is particularly effective when the graph of couplings is at least locally treelike, so any loops present are long. To extend the cavity recursion as an approximation to graphs which are not trees, it makes sense to consider normalized quantities and define $\pi_{i\to j}(s_{i})=Z_{i\to j}(s_{i})/\sum_{s}Z_{i\to j}(s)$ with the recursion

[TABLE]

leading to the estimate of the two-point marginals $\pi_{ij}(s_{i},s_{j})\propto e^{J_{ij}s_{i}s_{j}}\pi_{i\to j}(s_{i})\pi_{j\to i}(s_{j})$. (This step is necessary, as deleting the link $(i,j)$ on a tree leads to two disjoint subtrees with separate partition functions and zero connected correlations. On a locally treelike graph, correlations between $s_{i}$ and $s_{j}$ can still be small once the link $(i,j)$ has been cut, but the tree does not split into disjoint parts with separate partition functions.) Associating each edge $(i,j)$ in the graph with particular values of $\pi_{i\to j}(s_{i})$ and $\pi_{j\to i}(s_{j})$, one can update these values in an iterative scheme which replaces them with the right-hand side of (76) at each step. The fixed point of this procedure obeys (76). In this context, $\pi_{i\to j}(s_{i})$ is termed a ‘message’ that is being ‘passed’ between nodes. Belief propagation is an example of a message passing algorithm.

Belief propagation has been used to solve the inverse Ising problem in two different ways. In the first, belief propagation is used to approximately calculate the magnetisations and correlations given some parameters of the Ising model, and then fields and couplings are updated according to the Boltzmann learning rule (24). To estimate the correlations, one has to define additional messages also for the susceptibilities of each spin together with their update rules. This approach, termed susceptibility propagation, was developed by Welling and Teh Welling2003a ; Welling2004a for the forward problem. Susceptibility propagation thus solves the ‘forward’ problem, that is, it offers a computationally efficient approximation to the correlations and magnetisations, which are then used for Boltzmann machine learning.

A sophisticated variant of susceptibility propagation was developed by Mézard and Mora Mezard2009a ; Aurell2010a specifically for the inverse Ising problem, where also the couplings are updated at each step. This approach gives the same reconstruction as the Bethe–Peierls reconstruction from section II.1.10. However, the iterative equations can fail to converge, even when the analytical approach (71)–(73) gives valid approximate solutions Nguyen2012a .

II.1.12 The independent-pair approximation and the Cocco–Monasson adaptive-cluster expansion

We expect the Bethe–Peierls ansatz to work well when the couplings are locally treelike and loops are long. Conversely, we expect the ansatz to break down when the couplings generate many short loops. Cocco and Monasson developed an iterative procedure to identify clusters of spins whose couplings form short loops and evaluate their contribution to the entropy (34) Cocco2011a ; Cocco2012a . In the context of disordered systems, the expansion of the entropy in terms of clusters of connected spins is known as Kikuchi’s cluster variational method kikuchi1951 ; yedidia2001bethe .

The Cocco–Monasson adaptive-cluster expansion directly approximates the entropy potential (34). We start by considering the statistics of a single spin variable described by its magnetisation $m_{i}$, with the entropy

[TABLE]

The simplest entropy involving coupled spins is the two-spin entropy

[TABLE]

When the two spins are statistically independent, we obtain $S_{(2)}(m_{i},m_{j},\chi_{ij})=S_{(1)}(m_{i})+S_{(2)}(m_{j})$. Hence the residual entropy which accounts for the correlation between the two spins is $\Delta S_{2}(m_{i},m_{j},\chi_{ij})\equiv S_{(2)}(m_{i},m_{j},\chi_{ij})-S_{(1)}(m_{i})-S_{(1)}(m_{j})$. To make the notation uniform, we also define $\Delta S_{(1)}(m_{i})\equiv S_{(1)}(m_{i})$. A very simple approximation is based on the assumption that the $N$-spin entropy  (34) is described by pairwise terms

[TABLE]

where the pair $(i,j)$ denotes a cluster of two distinct spins. The couplings and fields can then be obtained via differentiation as in (35),

[TABLE]

This result is called the independent-pair approximation for the couplings, see roudi2009isingquality and Roudi2009b . (Expressions (8a) and (8b) in Roudi2009b differ slightly from (80) due to a typo.) When the topology of couplings is known and forms a tree, equation (79) gives the exact entropy of the system when the second sum is restricted to pairs of interacting spins (see sections II.1.9 and II.1.10). In this case, (80) gives the exact couplings.

However, in most cases the topology is not known, so second sum in (79) runs over all pairs of spins. In this case, equation (79) is only a (rather bad) approximation to the entropy. However, the independent-pair approximation (79) can serve as starting point for an expansion that includes clusters of spins of increasing size

[TABLE]

In this expansion, the contribution from clusters consisting of $3$ spins is

[TABLE]

where we have dropped the argument of $\Delta S_{(3)}$ to simplify the notation. The residual entropies of higher clusters are defined analogously.

The evaluation of $S_{(3)}$ is not as straightforward as $S_{(2)}$ and $S_{(1)}$, but still tractable. In general, the evaluation of $\Delta S_{(k)}$ requires the computation of order $2^{k}$ steps and therefore becomes quickly intractable. Cocco and Monasson argue that the contribution of $\Delta S_{k}$ decreases with $k$ and can be neglected from a certain order on Cocco2011a . This inspires an adaptive procedure to build up the library of all clusters which contribute significantly to the total entropy of the system. Starting from $1$-clusters, one constructs all $2$-clusters by merging pairs of $1$-clusters. If the entropy contribution of the new cluster is larger than some threshold, the new cluster is added to a list of clusters. One then constructs the $3$-clusters, $4$-clusters in the same way until no new cluster gives a contribution to the entropy exceeding the threshold. The total entropy and the reconstructed couplings and fields are updated each time a new cluster is added. This procedure fares well in situations when there are many short loops, like in a regular lattice in more than one dimension or when the graph of couplings is highly heterogeneous and some subsets of spins are strongly connected among each other.

The approaches introduced so far are each built on an ansatz for the Boltzmann distribution such as the mean-field ansatz (47) or the Bethe–Peierls ansatz (63): They are all uncontrolled approximations. In the next sections, we discuss controlled approximations based on either an expansion in small couplings or small correlations.

II.1.13 The Plefka expansion

In 1982, Plefka gave a systematic expansion of the Gibbs free energy of the Sherrington–Kirkpatrick model sherringtonkirkpatrick1975 in the couplings between spins Plefka1982a . Already 8 years earlier, Bogolyubov Jr. et al. had used similar ideas in the context of the ferromagnetic Ising model on a lattice bogolyubov1976high . The resulting estimate of the partition function can also be used to derive new solutions of the inverse Ising problem. Zeroth- and first-order terms of Plefka’s expansion turn out to yield mean-field theory (48), the second-order term gives the TAP free energy (54), and correspondingly the reconstructions of couplings (49) and (57).

Plefka’s expansion is a Legendre transformation of the cumulant expansion of the Helmholtz free energy $F(\boldsymbol{J},\boldsymbol{h})$. Plefka introduced a perturbation parameter $\lambda$ to the interacting part of the Hamiltonian

[TABLE]

where $H_{0}(\boldsymbol{s})=-\sum_{i}h_{i}s_{i}$ and $V(\boldsymbol{s})=-\sum_{i<j}J_{ij}s_{i}s_{j}$. The perturbation parameter $\lambda$ serves to distinguish the different orders in the strength of couplings and will be set to one at the end of the calculation. The standard cumulant expansion of the Helmholtz free energy $F^{\lambda}(\boldsymbol{J},\boldsymbol{h})$ then reads

[TABLE]

where $F^{(0)}(\boldsymbol{h})=\sum_{i}2\operatorname{cosh}h_{i}$ and $\left\langle\right\rangle_{0}$ denotes the average with respect to the Boltzmann distribution corresponding to the non-interacting part $H_{0}$ of the Hamiltonian. Next we perform the Legendre transformation of $F^{\lambda}(\boldsymbol{J},\boldsymbol{h})$ with respect to $\boldsymbol{h}$ to obtain the perturbative series for the Gibbs free energy

[TABLE]

where we suppressed the dependence of $G$ on $\boldsymbol{J}$ and $\boldsymbol{m}$ to simplify the notation. Using the definition of the Gibbs free energy (38) we solve

[TABLE]

with the local fields satisfying

[TABLE]

Plugging the perturbative series

[TABLE]

into (87) and re-expanding the right-hand side in $\lambda$, one finds successive orders of $\boldsymbol{h}$. In particular, the lowest two orders are found easily as

[TABLE]

This gives the Gibbs free energy up to first order in $\lambda$,

[TABLE]

Continuing to the second order, one obtains the Onsager term (54)

[TABLE]

A systematic way to perform this expansion has been developed by Georges and Yedidia Georges1991a leading to

[TABLE]

The mean-field and TAP reconstructions of couplings (53) and (57) then follow from (40) and the first derivatives of $G$ with respect to the magnetizations gives the corresponding reconstructions of the magnetic fields (50) and (58). Higher order terms of the Plefka expansion are discussed by Georges and Yedidia Georges1991a ; Yedidia2001a , but have not yet been applied to the inverse Ising problem. For a fully-connected ferromagnetic model with all couplings set to $J/N$ (to ensure an extensive Hamiltonian) one finds that already the second-order term $G^{(2)}$ vanishes relative to the first in the thermodynamic limit. For the Sherrington-Kirkpatrick model of a spin glass sherringtonkirkpatrick1975 , couplings are of the order of $N^{-1/2}$, again to make the energy extensive. In that case, $G^{(0)}$, $G^{(1)}$, and $G^{(2)}$ turn out to be extensive, but the third order vanishes in the thermodynamic limit. In specific instances (not independently and identically distributed couplings), higher order terms of the Plefka expansion may turn out to be important.

The terms in $J_{ij}$, $J_{ij}^{2}$, $J_{ij}^{3}$, etc. appearing in (II.1.13)-(92) resum to the results of the Bethe–Peierls ansatz Georges1991a . If the couplings form a tree, these terms are the only ones contributing to the Gibbs free energy; terms such as $J_{ij}J_{jk}J_{ki}$ in (92) quantify the coupling strengths around a closed loop of spins and are zero if couplings form a tree. This finally shows that the Bethe–Peierls ansatz is exact on a tree. However, when the couplings contain loops, these terms contribute to the Gibbs free energy. Beyond evaluating the Plefka expansion at different orders, one can also resum the contributions from specific types of loops to infinite order. We will discuss such an approach in the next section.

II.1.14 The Sessak–Monasson small-correlation expansion

The Sessak–Monasson expansion is a perturbative expansion of the entropy $S(\boldsymbol{\chi},\boldsymbol{m})$ in terms of the connected correlation $C_{ij}\equiv\chi_{ij}-m_{i}m_{j}$ Sessak2009a . For zero correlations $C_{ij}$, the couplings $J_{ij}$ should also be zero. This motivates an expansion of the free energy $F(\boldsymbol{J},\boldsymbol{h})$ in terms of the connected correlations around a non-interacting system; the Legendre transformation (34) then yields the perturbative series of the entropy function $S(\boldsymbol{\chi},\boldsymbol{m})$. Equivalently, this is an expansion of equations (23) for the fields and the couplings in the connected correlation $\boldsymbol{C}$. To make the perturbation explicit, we replace $\boldsymbol{C}$ by $\lambda\boldsymbol{C}$, where the perturbation parameter $\lambda$ is to be set to one at the end. The couplings and fields are the solutions of

[TABLE]

where the latter can be replaced by

[TABLE]

We then expand the solution $\boldsymbol{h}$ and $\boldsymbol{J}$ in a power series around the uncorrelated case $\lambda=0$

[TABLE]

At the zeroth-order in the connected correlations, spins are uncorrelated, so $\boldsymbol{J}^{(0)}=0$ and $h^{(0)}=\operatorname{artanh}(m_{i})$. The expansion of model parameters induces an expansion of the Hamiltonian

[TABLE]

where

[TABLE]

Likewise, the free energy $F(\boldsymbol{J},\boldsymbol{h})=-\ln Z(\boldsymbol{J},\boldsymbol{h})$ (or the cumulant generator) can be expanded in $\lambda$,

[TABLE]

where the subscript $\left\langle\right\rangle_{0}$ refers to the thermal average under the zeroth-order Hamiltonian $H^{(0)}(\boldsymbol{s})$. For example, the first-order term is

[TABLE]

Equations (93) and (94) then read

[TABLE]

Evaluating these equations successively gives solutions for the different orders of $h_{i}^{(k)}$ and $J^{(k)}_{ij}$, which in turn yield expressions for reconstructed fields and couplings when the magnetisations $\boldsymbol{m}$ and the connected correlations $\boldsymbol{C}$ are identified with their empirical values. Recalling the zeroth-order solution $h^{(0)}_{i}=\operatorname{artanh}(m_{i})$, the first-order contribution of the free energy (II.1.14) leads to

[TABLE]

Higher-order terms are given in Sessak2009a . Also in Sessak2009a , Sessak and Monasson give a diagrammatic framework suitable for the derivation of higher order terms in the couplings and sum specific terms to infinite order. Roudi et al. roudi2009isingquality simplified the results yielding the reconstruction

[TABLE]

for couplings with $i\neq j$, where $J^{\mathrm{IP}}_{ij}$ is the independent pair approximation (80). This result can be interpreted as follows roudi2009isingquality : The first term is the mean-field reconstruction, which is of a sub-series of the Sessak–Monasson expansion. One then adds the sub-series that constitutes the independent pair approximation (80). The last term is the overlap of the two series: the mean-field reconstruction of two spins considered independently, which needs to be subtracted. The resummation of other sub-series for special cases is also possible Sessak2009a ; in JaquinRancon2016 loop diagrams are resummed to obtain a reconstruction that is particularly robust against sampling noise.

II.2 Logistic regression and pseudolikelihood

Pseudolikelihood is an alternative to the likelihood function (II.1) and leads to the exact inference of model parameters in the limit of an infinite number of samples arnoldstrauss1991a ; hyvarinen2006consistency ; mozeika2014 . The computational complexity of this approach scales polynomially with the number of spin variables $N$, but also with the number of samples $M$. In practice, this is usually much faster than exact maximisation of the likelihood function, whose computational complexity is exponential in the system size. Pseudolikelihood inference was developed by Julian Besag in 1974 in the context of statistical inference of data with spatial dependencies Besag1974a . It is closely related to logistic regression. While pseudolikelihood and regression have been used widely in statistical inference besag1986 ; kalbfleisch1986 ; Ravikumar2010a , this approach was not well-known in the physics community until quite recently Aurell2012b ; mozeika2014 ; ekeberg2013 ; Decelle2014a .

Our derivation focuses on the physical character of regression analysis, which we then link to Besag’s pseudolikelihood. Key observation is that although the likelihood function (II.1) depends on the model parameters in a complicated way, one can simplify this dependence by separating different groups of parameters. Let us consider how the statistics of a particular spin variable $\sigma_{i}$ depends on the configuration of all other spins. We split the Hamiltonian into two parts

[TABLE]

such that the first part $H_{i}(\boldsymbol{s})$ depends on the magnetic field $h_{i}$ and the couplings of spin $i$ to other spins $J_{ij}$, while the part $H_{\setminus i}(\boldsymbol{s}\setminus s_{i})$ does not. $\boldsymbol{s}\setminus s_{i}$ denotes all spin variables except spin $s_{i}$. This splitting of variables is possible since the statistics of $\sigma_{i}$ conditioned on the remaining spins $\{s_{j}\}_{j\neq i}$ is fully captured by $h_{i}$ and $J_{ij},(j\in 1,\ldots,N)$.

The expectation values of $\sigma_{i}$ can be computed based on the partition function

[TABLE]

where only spin $i$ has been summed over. Differentiating the partition function in this form with respect to $h_{i}$ and $J_{ij}$ yields the expectation values

[TABLE]

where on both sides the thermal average is over the Boltzmann measure $e^{-H(\boldsymbol{s})}/Z$. The first equation follows from

[TABLE]

with the shorthand $\Theta_{i}=h_{i}+\sum_{k\neq i}J_{ik}\sigma_{k}$. The second equation follows analogously.

In statistical physics, equations (105) and (106) are known as Callen’s identities Callen1963a and have been used to compute coupling parameters from observables in Monte Carlo simulations for the numerical calculation of renormalisation-group trajectories Swendsen1984a . While these equations are exact, the expectation values on the right hand sides still contain the average over the $N-1$ spins other than $i$.

The crucial step and the only approximation involved is to replace the remaining averages in (106) with the sample means

[TABLE]

We now have a system of non-linear equations to be solved for $h_{i}$ and $J_{ij}$ for fixed $i$ and various $j$, $j\neq i$. Standard methods to solve such equations are Newton–Raphson or conjugated gradient descent Press2002a ; Hastie2009a .

With these steps, we have broken down the problem of estimating the magnetic fields and the full coupling matrix to $N$ separate problems of estimating one magnetic field and a single row of the coupling matrix for a specific spin $i$. Crucially, the Boltzmann average over $2^{N-1}$ states has been replaced with an average over all configurations of the samples. As a result, the computation of (108) uses not only the sample magnetizations and correlations (sufficient statistics), but all spin configurations that have been sampled. The time to evaluate (108) is thus linear in the number of samples $M$. In general, the coupling matrix inferred in this way will be asymmetric, $J_{ij}\neq J_{ji}$ due to sampling noise. A practical solution is to use the average $\frac{1}{2}(J_{ij}+J_{ji})$ as estimate of the (symmetric) coupling matrix.

The set of equations (108) can be viewed as solving a gradient-descent problem of a logistic regression Ravikumar2010a . The statistics of $\sigma_{i}$ conditioned on the values of the remaining spins $\{s_{j}\}_{j\neq i}$ can be written as

[TABLE]

From this expression for the conditional probability, the $i$-th row of couplings $J_{i\ast}$ and the field $h_{i}$ are simply the coefficients and the intercept in the multi-variate logistic regression of the response variable $\sigma_{i}$ on the variables $\{s_{j}\}_{j\neq i}$. The log-likelihood for the $i$-th row of couplings $J_{i\ast}$ and the magnetic field $h_{i}$ of this regression problem is

[TABLE]

Setting the derivatives of this likelihood function with respect to the magnetic field $h_{i}$ and the entries of $J_{i\ast}$ to zero recovers (108).

Consider all couplings and fields together, one can sum (110) over all rows of couplings to obtain the so-called (log)-pseudolikelihood besag1986

[TABLE]

which can be maximized with respect to all rows of the coupling matrix, yielding an asymmetric matrix as discussed above. The pseudolikelihood can also be maximized within the space of symmetric matrices, although the maximisation problem is harder Aurell2012b ; rather than solving $N$ independent gradient-descent problems in $N$ variables we have a single problem in $N(N+1)/2$ variables. In practice, maximising the pseudolikelihood without the symmetry constraint on the coupling matrix is preferred because of its simplicity and efficiency.

The reconstruction based on maximizing the pseudolikelihood is of a different nature than the previous methods, which approximated the Boltzmann measure. The Boltzmann measure specifies the probability of observing a particular spin configuration $\boldsymbol{s}$ in equilibrium. On the other hand, the conditional probability $\ln p(s_{i}^{\mu}|\{s_{j}^{\mu}\}_{j\neq i})$ only gives the probability of observing a given spin $s_{i}$ conditioned on the remaining spin variables and cannot be used to generate the configurations of all spins. Yet, as a function of the model parameters, the log-pseudolikelihood (110) has the same maximum as the likelihood (II.1) in the limit of $M\to\infty$, when the sample average in (108) coincides with the corresponding expectation values under the Boltzmann distribution.

Curiously, the pseudolikelihood can also be used to obtain a variant of the mean-field and TAP reconstructions. We follow Roudi and Hertz Roudi2011a and replace expressions such as $\left\langle\operatorname{th}\left(h_{i}+\sum_{j\neq i}J_{ij}\sigma_{j}\right)\right\rangle$ by $\operatorname{th}\left(h_{i}+\sum_{j\neq i}J_{ij}\left\langle\sigma_{j}\right\rangle\right)$, thus replacing the effective local field by its mean. The resulting approximation of Callen’s identity (105) is

[TABLE]

which is the mean-field equation (51) for the magnetizations. Then, replacing $\sigma_{i}$ by $\left\langle\sigma_{i}\right\rangle+(\sigma_{i}-\left\langle\sigma_{i}\right\rangle)$ in (106) and expanding the equation to the lowest orders in $(\sigma_{i}-\left\langle\sigma_{i}\right\rangle)$ gives

[TABLE]

Identifying the pseudolikelihood mean-field magnetisations $m^{\mathrm{MF}}_{i}$ and connected correlations $C^{\mathrm{PL}-\mathrm{MF}}_{ij}$ with the sample magnetisations $m_{i}$ and sampled correlations $C_{ij}$, this linear equation can be solved for $J_{ik}^{\mathrm{PL}-\mathrm{MF}}$ for a fixed $i$, giving

[TABLE]

where $\boldsymbol{C}_{\setminus i}$ is the submatrix of the correlation matrix with row and column $i$ removed. This reconstruction is closely related, but not identical to the mean-field reconstruction (53). In particular, the diagonal terms are naturally excluded. Numerical experiments show that this reconstruction gives a good correlation between the reconstructed and true couplings, but the magnitude of the couplings is systematically underestimated. Continuing the expansion to second order in $(\sigma_{i}-\left\langle\sigma_{i}\right\rangle)$ leads to a variant of the TAP reconstruction (57).

II.3 Comparison of the different approaches

Table 1 gives a classification of the different reconstruction methods based on the thermodynamic potentials they employ and the approximations used to evaluate them. Some of these approximations are exact in particular limits, and fail in others. This has consequences for how well a reconstruction method works in a given regime. For instance, the TAP equations become exact for fully connected systems (with couplings between all spin pairs drawn from a distribution with variance $1/N$) in the thermodynamic limit. Hence, we expect the TAP reconstruction to perform well when couplings are uniformly distributed across spin pairs, and couplings are sufficiently weak so there is no ergodicity breaking (see subsection II.4). Similarly, the Bethe–Peierls approximation is exact when the non-zero couplings between spin pairs form a tree. Hence, we expect the Bethe–Peierls reconstruction to work perfectly in this case, and to work well when the graph of couplings is locally treelike (so there are no short loops). The adaptive cluster expansion, on the other hand, is expected to work well even when there are short loops.

In this section, we compare the reconstruction methods discussed so far. We consider the reconstruction problem of an Ising model (2) with couplings $J^{0}_{ij}$ between pairs of spins defined on a certain graph. We explore two aspects of the model parameters: the type of the interaction graph and the coupling strength (temperature). We consider three different graphs: the fully connected graph, a random graph with fixed degree as a representative of graphs with long loops, and the $2$D square lattice as a representative of graphs with short loops. In this section, we denote the couplings of the model underlying the data by the superscript ‘[math]’ to distinguish them from the inferred couplings. For each graph, every edge is assigned a coupling $J^{0}_{ij}$ drawn from a certain distribution. We use the Gaussian distribution with zero mean and standard deviation $\beta/\sqrt{N}$ for couplings on the fully connected graph, leading to the Sherrington–Kirkpatrick (SK) model. For the tree and square lattice, the uniform distribution on the interval $[-\beta,\beta]$ is used for the couplings. By tuning $\beta$, we effectively change the coupling strength. For the SK-model, external magnetic fields field $h_{i}^{0}$ are drawn uniformly from the interval $[-0.3\beta,0.3\beta]$.

At each value of $\beta$, we generate $M$ samples (spin configurations) drawn from the Boltzmann distribution. Each sample is obtained by simulating $10^{4}N$ Monte Carlo steps using the Metropolis transition rule starting from a random configuration. Although this is not always sufficient to guarantee that the system has reached equilibrium, the results are not sensitive to increasing this breaking-in time.

Figure 4 top left compares the reconstructed couplings with the couplings of the original model at $\beta=1.3$. We chose this value as it results in couplings which are sufficiently large so that the different methods perform differently (we will see that at low $\beta$ all methods perform similarly). On the fully connected graph, the adaptive cluster expansion performs poorly as it sets too many couplings to zero. Mean-field reconstruction somewhat overestimates large couplings barton2014a . This error also affects the estimate of the magnetic fields. The TAP and Bethe–Peierls reconstructions correct this overestimate quite effectively. The reconstruction by pseudolikelihood stands out by providing an accurate reconstruction of the model parameters. One consequence is that in figure 4 the symbols indicating the results from pseudolikelihood are largely obscured by those from the exact maximization of the likelihood (II.1). We also compare the correlations and magnetizations based on the reconstructed parameters with those of the original model. The correlations and magnetisations are sampled as described above. We find significant bias in the results from all methods except for likelihood maximization and pseudolikelihood maximization (bottom panels of figure 4).

Next, we explore how the reconstruction quality depends on the coupling strength $\beta$ and the number of samples $M$. To this end, we define the relative reconstruction error

[TABLE]

which compares the reconstructed couplings $\boldsymbol{J}$ to the couplings of the original model underlying the data $\boldsymbol{J}^{0}$. Figure 5, large panel, shows the relative reconstruction error $\gamma_{J}$ as a function of coupling strength $\beta$. At low $\beta$, the connected correlations are small, so all reconstruction methods are equally limited by sampling noise: The relative reconstruction error increases with decreasing $\beta$ as the couplings become small relative to the sampling error. This is also compatible with the result that at weak couplings the relative errors of all methods decrease with an increasing number of samples, as seen in the $6$ small panels of figure 5.

We find that the pseudolikelihood reconstruction outperforms the other (non-exact) methods over the entire range of $\beta$, and correctly reconstructs the model parameters even in the glassy phase at strong couplings Aurell2012b . This is because the conditional statistics of a single spin (conditioned on the other spins) is correctly described by (109) at any coupling strength. Although at strong couplings, the error of the pseudolikelihood reconstruction grows with $\beta$, this can be compensated by increasing the number of samples (Figure 5, bottom right). This is not so for all other approximate methods (Figure 5, small panels): At high $\beta$, the approximation each method is based on breaks down, which cannot be compensated by a larger number of samples.

Figure 6 shows the reconstruction errors of different methods for a random graph of fixed degree (column A) and the square lattice (column B). Again, pseudolikelihood performs well in these cases, as does the Bethe–Peierls reconstruction. The adaptive cluster expansion shows a remarkable behaviour: For both graphs, it has a very small reconstruction error at weak coupling strength, but breaks down at strong couplings. At weak couplings, where all other methods have similar reconstruction errors arising from sampling noise, however the adaptive cluster expansion appears to avoid this source of error. The adaptive cluster expansion explicitly assumes some couplings are exactly zero, so the reconstruction is biased towards graphs which are not fully connected. It thus uses extra information beyond the data.

Our comparison of the different methods is based on the knowledge of the underlying couplings, so the reconstructed couplings can be compared to the true underlying ones. In practice, the underlying couplings are not available. However, the likelihood (II.1) can be evaluated for different reconstructions, with the better reconstruction resulting in a higher value of the likelihood. Alternatively, the correlations and magnetizations of the reconstructed model can be compared with those observed in the data as in figure 4. Whether regularizing terms improve the reconstruction can be decided based on statistical tests such as the Bayesian information criterion schwarz1978estimating or the Akaike information criterion akaike1974new .

An aspect not discussed so far is the inference of the interaction graph, this is, the distinction between non-zero and zero couplings. In practice, there is always some ambiguity between small non-zero couplings and couplings which are exactly zero, particularly at high sampling noise. A standard practice is to set a threshold and cut off small couplings below the threshold. Two exceptions are the adaptive cluster expansion and the pseudolikelihood reconstruction. The adaptive cluster expansion has a built-in procedure to set some couplings to zero. For the pseudolikelihood reconstruction, the problem is related to feature selection. Many methods for feature selection are available from statistics Hastie2009a , however there is no consensus on the best method for the inverse Ising problem. One possibility is adding an $\ell_{1}$-regularization term (known as Lasso, see section II.1) to the pseudolikelihood Ravikumar2010a . However, there is a critical coupling strength, below which the $\ell_{1}$ regularization of the pseudolikelihood fails to recover the interaction graph BentoMontanari2009 .

II.4 Reconstruction and non-ergodicity

At strong couplings, the dynamics of a system can become non-ergodic, which affects the sampling of configurations, and hence the reconstruction of parameters. Reconstruction on the basis of non-ergodic sampling is a fundamentally difficult subject where little is known to date.

At low temperatures (or strong couplings), a disordered spin system can undergo a phase transition to a state where spins are ‘frozen’ in different directions. This is the famous transition to the spin glass phase mezardparisivirasoro . At the spin glass transition, the Gibbs free energy develops multiple valleys with extensive barriers between them (thermodynamic states). These free energy barriers constrain the dynamics of the system, and the system can remain confined to one particular valley for long times. A signature of this ergodicity breaking is the emergence of non-trivial spin magnetisations: At the phase transition, the self-consistent equations (51) and (55) develop multiple solutions for the magnetisations $m_{i}$, corresponding to the different orientations the spins freeze into. The appearance of multiple free energy minima and the resulting ergodicity breaking has long been recognized as an obstacle to mean-field reconstruction marinari2010a . In fact, all methods based on self-consistent equations (mean-field, TAP, Bethe–Peierls reconstruction) fail in the glassy phase Mezard2009a .

The mean-field equations (51) describe individual thermodynamic states, not the mixture of many states that characterizes the Boltzmann measure. At large couplings (or low temperatures), mean-field and related approaches turn out not to be limited by the validity of self-consistent equations like (51) or (55), but by the identification of observed magnetisations and correlations with the corresponding quantities calculated under one of the self-consistent equations. The former may involve averages over multiple thermodynamic states, whereas each solution of the self-consistent equations describes a separate thermodynamic state. A solution to this problem is thus to consider correlations and magnetisations within a single thermodynamic state, where the mean-field result (52) is valid within the limitations of mean-field theory. The different thermodynamic states (free energy minima) can be identified from the data by searching for clusters in the sampled spin configurations. Collecting self-consistent equations from different minima and jointly solving these equations using the Moore-Penrose pseudo-inverse allows the reconstruction of couplings and fields Nguyen2012b .

The emergence of multiple states can lead to another problem: the samples need not come from the Gibbs measure in the first place, but might be taken only from a single thermodynamic state. In braunstein2011a , the reconstruction from samples from a single thermodynamic state is studied for the concrete case of the Hopfield model. In the so-called memory regime, the Hopfield model has a number of thermodynamic states (attractor states) which scales linearly with the system size. Samples generated from a single run of the model’s dynamics will generally come from a single state only. (The particular state the system is ‘attracted’ to depends on the initial conditions.) In this regime, the cavity-recursion equations (76) typically have multiple solutions, just like the TAP equations (55) in the low-temperature phase. These solutions can be identified as fixed points of belief propagation, see  II.1.10. For a system where couplings do not form a tree, each individual fixed point only approximates the marginal probability distributions inside a single non-ergodic component of the full Gibbs measure. However, this approximation can become exact in the thermodynamic limit. (A necessary condition is that loops are sufficiently long, so connected correlations measured in a particular state decay sufficiently quickly with distance, where distance is measured along the graph of non-zero couplings.)

Assuming one is able to find a fixed point of the cavity-recursion equations, it is a straightforward computation to express the Boltzmann weight (1) in terms of the full set of the belief propagation marginals and of the ratio between the cavity-recursion partition function and the true partition function of the model. For any fixed point, labelled by the index $\alpha$, the Boltzmann weights are

[TABLE]

This result also applies to more general forms of the energy function braunstein2011a . It was first derived for the zero-temperature case in the context of the so-called tree-reweighted approximation to the partition function kolmogorov2006convergent .

The probability distribution (116) can be used as the starting point for reconstruction in at least two ways. The simpler one consists in replacing $p_{ij}^{\alpha}(s_{i},s_{j})$ and $p_{i}^{\alpha}(s_{i})$ by their sample estimates inside state $\alpha$ and then solving the identity between (116) and (1) with respect to $\boldsymbol{J}$ and $\boldsymbol{h}$: $Z(\boldsymbol{J},\boldsymbol{h})$ cancels out and one is left with the independent-pair approximation of subsection II.1.12. Solving this identity can be done using belief propagation, see section II.1.11 and braunstein2011a . This approach suffers from the fact that in general there is no belief propagation fixed point for real data sets.

A more promising approach is to guide the belief propagation equations to converge to a fixed point corresponding to an appropriate ergodic component close to the empirical data. In fact, ignoring the information that the data come from a single ergodic component results in a large reconstruction error as one is effectively maximizing the wrong likelihood: the (reduced) free energy $-\ln Z$ appearing in the likelihood (II.1) needs to be replaced by a free energy restricted to a single thermodynamic state. An algorithmic implementation of this idea consists in restricting the spin configurations to a subset of configuration $\Omega_{d}(\xi)$, e.g., formed by a hypersphere of diameter $d$ centered around the centroid $\xi$ of the data samples. The system is forced to follow the measure

[TABLE]

where the indicator function $I[\boldsymbol{s}\in{\Omega}_{d}(\xi)]$ for the set of configurations $\Omega_{d}(\xi)$ is one for $\boldsymbol{s}\in\Omega_{d}(\xi)$ and zero otherwise. The BP approach can be used to enforce that both magnetizations $\{m_{i}^{(d)}\}$ and correlations $\{c_{ij}^{(d)}\}$ are computed within the subspace $\Omega_{d}(\xi)$, see braunstein2011a for algorithmic details. In this way, the reconstruction of the couplings and local fields can be done by maximising the correct likelihood, i.e., replacing the partition function in the log-likelihood (II.1) with the partition function computed over the subset of configurations $\Omega_{d}(\xi)$. In this way, the parameters describing all thermodynamic states can be inferred from configurations of the system sampled from a single state.

Inverse problems in non-ergodic systems remain a challenging topic with many open questions, for instance if the approach applied above to the Hopfield model can be extended to infer generic couplings from samples from a single thermodynamic state.

II.5 Parameter reconstruction and criticality

Empirical evidence for critical behaviour has been reported in systems as diverse as neural networks bedard2006a ; beggs2012a and financial markets gabaix2003a ; bury2013a , leading to the intriguing hypothesis that some information-processing systems may be operating at a critical point Mora2011b . A key signature of criticality is broad tails in some quantities, for instance the distribution of returns in a financial market, or avalanches of neural activity, whose sizes are distributed according to a power law. A second sign of criticality involves an inverse statistical problem: A statistical model like the Ising system with parameters chosen to match empirical data (such as neural firing patterns or financial data) shows signs of a phase transition Mora2011b ; mastromatteo2011a ; stephens2013a ; tkavcik2014a ; tkacik2015 ; mora2015dynamical . Specifically, the heat capacity of an Ising model with parameters $J_{ij}$ and $h_{i}$ reconstructed from data shows a peak in the heat capacity as a function of the temperature. Temperature is introduced by changing the couplings to $\beta J_{ij}$ and $\beta h_{i}$. Varying the inverse temperatures $\beta$, the heat capacity $C_{h}\equiv-\beta^{2}\partial\langle H\rangle/\partial\beta$ shows a pronounced maximum near or at $\beta=1$, that is, at the model parameter values inferred from the data. The implication is that the parameters of the reconstructed model occupy a very particular area in the space of all model parameters, namely one resulting in critical behaviour.

An example is shown in figure 7. It is based on recordings of $160$ neurons in the retina of a salamander taken by the Berry lab tkavcik2014a ; tkacik2015 . As described in section I.1, time is discretised into intervals of $20$ ms duration, and the spin variable $s_{i}$ takes on the values $+1$ if neuron $i$ fires during a particular interval, and $-1$ if it does not fire. Correlations and magnetisations are then computed by averaging over $297$ retinal stimulation time courses, where during each time course, the retina is subjected to the same visual stimulation of $953\times 20$ms duration. During most time intervals, a neuron is typically silent. As a result, spins variables (each representing a different neuron) have negative magnetisations, with mean and standard deviation $-0.93\pm 0.06$ over all spins. We note that the magnetisations are fairly homogeneous across spins. Connected correlations between spins are small, with a slight bias towards positive correlations; mean and standard deviation over all spin pairs are $0.006\pm 0.017$. Both points turn out to be important.

We use the pseudolikelihood method of section II.2 to reconstruct the model parameters from the firing patterns the first $N=30$ neurons. Similar results are found for different system sizes. Figure 7A shows that the reconstructed model parameters indeed occupy a particular region of the space of model parameters: External fields $h_{i}$ are linked to the sum over couplings $\sum_{j}J_{ij}$ by a linear-affine relationship. If the magnetisations $\{m_{i}\}$ are all equal to some $m$, such a relationship follows directly from the mean-field equation $h_{i}=\operatorname{artanh}m_{i}-\sum_{j}J_{ij}m_{j}$, along with slope $-m$ and zero-offset $\operatorname{artanh}m$.

Now we simulate the Ising model with the reconstructed parameters, but rescale both couplings and fields by $\beta$. At inverse temperature $\beta=1$, magnetisations and correlations are close to the magnetisations and correlations found in the original data. However, away from $\beta=1$, different fields $h_{i}=\frac{1}{\beta}\operatorname{artanh}~{}m_{i}-\frac{1}{\beta}\sum_{j}J_{ij}m_{j}$ (in the mean-field approximation) would be needed to retain the same magnetisation. However, as the fields remain fixed, the magnetisations change instead, reaching $-1(0)$ in the limit $\beta\to\infty(0)$. Between these limits, the magnetisations, and hence the energy changes rapidly, leading to a peak in the heat capacity. For small, largely positive connected correlations $C_{ij}$ we expect couplings to scale as $1/N$, resulting in a phase transition described by the Curie-Weiss model. Indeed, a finite-size analysis shows the peak getting sharper and moving closer to $\beta=1$ as the system size is increased Mora2011b ; mastromatteo2011a ; stephens2013a ; tkacik2015 ; mora2015dynamical . From this, we posit that the peak in the heat capacity indeed signals a ferromagnetic phase transition, however not necessarily one underlying the original system. Instead, we generically expect such a ferromagnetic transition whenever a system shows sufficient correlations and magnetisations that are neither plus nor minus one: The corresponding couplings and fields then lie near (not at) a critical point, which is characterized by the emergence of non-zero magnetisations. Changing the temperature will drive the system to either higher or lower magnetisations, and away from the critical point.

This effect may not be limited to a ferromagnetic transition induced by the empirical magnetisation. In mastromatteo2011a Mastromatteo and Marsili point out a link between the criticality of inferred models and information geometry balasubramanian1997a ; myung2000a : susceptibilities such as the magnetic susceptibility can be interpreted as the entries in the so-called Fisher information matrix used to calculated the covariances of maximum-likelihood estimates Cover2006a . A high susceptibility implies that different parameter values can be distinguished on the basis of limited data, whereas low susceptibilities mean that the likelihood does not differ sufficiently between different parameter values to significantly favour one parameter value over another. Thus, a critical point corresponds to a high density of models whose parameters are distinguishable on the basis of the data mastromatteo2011a .

III Non-equilibrium reconstruction

In a model of interacting magnetic spins such as (2), each pair of spins $i<j$ contributes a term $J_{ij}s_{i}s_{j}$ to the Hamiltonian. The resulting effective fields on a spin $i$, $\sum_{j}J_{ij}s_{j}+h_{i}$, involve symmetric couplings $J_{ij}=J_{ji}$; spin $i$ influences spin $j$ as much as $j$ influences $i$. A stochastic dynamics based on such local fields, such as Monte Carlo dynamics, obeys detailed balance, which allows one to determine the steady-state distribution Gardiner1985a .

In applications such as neural networks or gene regulatory networks discussed in sections I.1 and I.2 we have no reason to expect symmetric connections between neurons or between genes; a synaptic connection from neuron $i$ to neuron $j$ does not imply a link in the reverse direction. Stochastic systems generally relax to a steady state at long times, where observables no longer change with time. However, in systems with asymmetric couplings, the resulting non-equilibrium steady state (NESS) violates detailed balance and differs from the Boltzmann distribution. As a result, none of the results of section II on equilibrium reconstruction apply to systems with asymmetric couplings.

In this section we consider the inverse Ising problem in a non-equilibrium setting. We first review different types of spin dynamics, and then turn to the problem of reconstructing the parameters of a spin system from either time-series data or from samples of the non-equilibrium steady state.

III.1 Dynamics of the Ising model

The Ising model lacks a prescribed dynamics, and there are many different dynamical rules that lead to a particular steady-state distribution. One particularly simple dynamics is the so-called Glauber dynamics glauber1963a , which allows the derivation of a number of analytical results. Other dynamical rules can be used for parameter reconstruction in the same way at least in principle. What dynamical rule is suitable for particular systems is however an open question. An approach which sidesteps this issue is to apply the maximum entropy principle to stochastic trajectories, known as the principle of maximum caliber  presse2013a . In marre2009 ; nassercessac2014a such an approach is applied to analyse neural dynamics.

III.1.1 Sequential Glauber dynamics

Glauber dynamics can be based on discrete time steps, and at each step either one or all spins have new values assigned to them according to a stochastic rule. We first consider a sequential dynamics: In the transition from time $t$ to $t+1$, the label of a spin variable is picked randomly, say $i$. The value of the spin variable $\sigma_{i}$ is then updated, with $s_{i}=\pm 1$ sampled from the probability distribution

[TABLE]

where the effective local field is denoted

[TABLE]

One way to implement this dynamics is to set

[TABLE]

where $\xi_{i}(t)$ is drawn independently at each step from the distribution $p(\xi)=1-\operatorname{th}^{2}(\xi)$.

If the couplings $J_{ij}$ between spins are symmetric and there are no self-couplings $J_{ii}$, the sequential Glauber dynamics (118) obeys detailed balance, and the distribution of spin configurations at long times relaxes to a steady state described by the Boltzmann distribution (11). However, the NESS arising at long times from Glauber dynamics with non-symmetric couplings is generally not known. Nevertheless, there are several exact relations that follow from (118), and those can be exploited for inference.

Suppose that spin $i$ is picked for updating in the step from time $t$ to $t+1$. Averaging over the two possible spin configurations at time $t+1$, we obtain the magnetisation

[TABLE]

where the average is conditioned on the configuration of all other spins at time $t$ through the effective local field $\theta_{i}$. If spin $i$ is not updated in this interval, the conditioned magnetisation is trivially $m_{i}(t+1)_{\boldsymbol{s}(t)}=s_{i}(t)$. In the steady state, effective fields $\Theta_{i}=\sum_{j}J_{ij}\sigma_{j}(t)+h_{i}$ and spins configurations $\boldsymbol{\sigma}$ are random variables whose distribution no longer depends on time; their averages give

[TABLE]

Similarly, the pair correlation between spins at sites $i\neq j$ at equal times obeys in the NESS

[TABLE]

with the two terms arising from instances where the last update of $i$ was before the last update of $j$, and vice versa. Likewise, the pair correlation of spin configurations at consecutive time intervals is in the steady state

[TABLE]

These relationships are exact, but are hard to evaluate since the averages on the right-hand sides are over the statistics of spins in the NESS, which is generally unknown. Below, these equations will be used in different ways for the reconstruction of the model parameters.

The dynamics (118) defines a Markov chain of transitions between spin configurations differing by at most one spin flip. Equivalently, one can also define an asynchronous dynamics described by a Master equation, where time is continuous and the time between successive spin flips is a continuous random variable.

III.1.2 Parallel Glauber dynamics

Glauber dynamics can also be defined with parallel updates, where all spin variables can change their configurations in the time interval between $t$ and $t+1$ according to the stochastic update rule

[TABLE]

This update rule defines a Markov chain consisting of stochastic transitions between spin configurations. The resulting dynamics is not realistic for biological networks, as the synchronous update requires a central clock. It is however implemented easily in technical networks, and is widely used for its simplicity. For a symmetric coupling matrix, the steady state can still be specified in closed form using Peretto’s pseudo-Hamiltonian peretto1984a .

Magnetisations and correlations in the NESS obey simpler relationships for parallel updates than for sequential updates; with the same arguments as above, one obtains

[TABLE]

III.2 Reconstruction from time series data

The reconstruction of system parameters is surprisingly easy on the basis of time series data, which specifies the state of each spin variable at $M$ successive time points. An application where such data is widely available is the reconstruction of neural networks from temporal recordings of neural activity. Given a stochastic update rule such as (118) or (125), the (log-) likelihood given a time series $\mathrm{D}=\{\boldsymbol{s}(t)\}$ is

[TABLE]

The likelihood can be evaluated over any time interval, both in the NESS or even before the steady state has been reached. This approach is not limited to non-equilibrium systems; whenever time series data are available, it can be used equally well to reconstruct symmetric couplings.

III.2.1 Maximisation of the likelihood

For Glauber dynamics with parallel updates, the likelihood is

[TABLE]

Unlike the likelihood (II.1) arising in equilibrium statistics, the non-equilibrium likelihood can be evaluated easily, as the normalisation is already contained in the term $\operatorname{cosh}(\theta)$. Derivatives of the likelihood (128) are

[TABLE]

These derivatives can be evaluated in $MN^{2}$ computational steps and can be used to maximize the likelihood by gradient ascent Roudi2011a .

The derivatives of the likelihood can also be written as time averages over the data, which makes a conceptual connection with logistic regression apparent Roudi2011a . The derivative with respect to fields gives $\frac{\partial L_{\mathrm{D}}}{\partial h_{i}}=\langle\sigma_{i}(t+1)\rangle^{\mathrm{D}}_{T}-\langle\operatorname{th}(\sum_{j}J_{ij}\sigma_{j}(t)+h_{i})\rangle^{\mathrm{D}}_{T}$ and similarly for the derivative with respect to the couplings (the subscript refers to the temporal average, the superscript to configurations taken from the data). Parallel Glauber dynamics (125) defines a logistic regression giving the statistics of $\sigma_{i}(t+1)$ as a function of $\{s_{j}(t)\}$. In a hypothetical data set of $P$ trajectories of the system $\mathrm{D}=\{\boldsymbol{s}^{p}(t)\}_{t=1}^{M}$, each realization can be considered as $M-1$ realisations $(s^{p}_{i}(t+1),\{s^{p}_{j}(t)\})$ of such regression pairs. The discussion in section II.2 then applies directly, giving rise to regression equations for fields and couplings.

For Glauber dynamics with sequential updates, the situation is not quite so simple. In each time interval one spin is picked for updating, but if this spin is assigned the configuration it had before it is impossible to tell which spin was actually chosen zeng2013 . The solution is to sum the likelihood (127) over all spins (each of which might have been the one picked for updating with probability $1/N$). It turns out that both the intervals when a spin was actually flipped, and the intervals when no spin was flipped are required for the inference of couplings and fields.

III.2.2 Mean-field theory of the non-equilibrium steady state

As in the case of equilibrium reconstruction, mean-field theory offers an approximation to the maximum likelihood reconstruction that can be evaluated quickly. The speed-up is not quite as significant as it is in equilibrium, because the likelihood (128) can already be computed in polynomial time in $N$ and $M$.

In II.1.13, the mean-field equation (51) and the TAP equation (55) were derived in an expansion around a factorising ansatz for the equilibrium distribution. Kappen and Spanjers showed that, remarkably, exactly the same equations emerge as first- and second-order expansions around the same ansatz for the NESS as well Kappen2000a . Suppose that the (unknown) steady-state distribution of configurations $p(\boldsymbol{s})$ in the NESS is ‘close’ to another distribution with the same magnetisations, $q(\boldsymbol{s})=\prod_{i}\frac{1+m_{i}s_{i}}{2}$. According to (122), this distribution describes the NESS of a different model with fields $h_{i}^{(q)}=\operatorname{artanh}m_{i}$ and couplings $J^{(q)}_{ij}=0$. Next, we consider a small change in these fields and couplings and ask how the magnetisations change. To first order this change is given by the derivatives of (122) evaluated at $h_{i}=h_{i}^{(q)}$ and $J_{ij}=0$

[TABLE]

Setting $\Delta h_{i}=h_{i}-h_{i}^{(q)}$ and $\Delta J_{ij}=J_{ij}-J_{ij}^{(q)}=J_{ij}$ and demanding that magnetisations remain unchanged under this change of fields and couplings gives to first order $h_{i}^{(q)}=h_{i}+\sum_{j}J_{ij}m_{j}$ and thus

[TABLE]

Carrying the expansion (130) to second order in $\Delta h_{j}$ and $\Delta J_{jk}$ yields the TAP equations (55) Kappen2000a . As the equations for the magnetisations  (122) and correlations (126) are identical for sequential and parallel dynamics, the mean-field and TAP equations apply equally to both types of dynamics. Roudi and Hertz roudi2011b extended these results to time scales before the NESS is reached using a generating functional approach.

The next step is to apply the mean-field approximation to the correlations (126) for parallel Glauber updates roudi2011b ; zeng2011a . For sequential updates, analogous result can be derived from (123) and (124). We expand the effective local field $\Theta_{i}=h_{i}+\sum_{j}J_{ij}m_{j}+\sum_{j}J_{ij}\delta\sigma_{j}$ around the mean field $\theta^{\textrm{MF}}_{i}=h_{i}+\sum_{j}J_{ij}m_{j}$ writing $\delta\sigma_{i}\equiv\sigma_{i}-m_{i}$. Expanding the $\operatorname{th}\Theta_{i}$-term of the two-time pair correlation (126) in a formal expansion in powers of $\delta\sigma$ gives

[TABLE]

To first order, this equation can be read as a matrix equation in the connected correlations functions, $D_{ij}\equiv\phi_{ij}-m_{i}m_{j}=\sum_{m,l}A_{im}J_{ml}C_{lj}$ with $A_{im}=\delta_{im}(1-m_{i}^{2})$ and $C_{ij}=\chi_{ij}-m_{i}m_{j}$. Inverting this relationship leads to the mean-field reconstruction

[TABLE]

based on sample averages of magnetisations and connected correlations in the NESS. The reconstruction based on the TAP equation can be derived analogously roudi2011b ; zeng2011a ; the result is of the same form as (133) with $A_{im}=\delta_{im}(1-m_{i}^{2})(1-F_{i})$, where $F_{i}$ is the smallest root of

[TABLE]

III.2.3 The Gaussian approximation

For an asymmetric coupling matrix, with no correlation between $J_{ij}$ and $J_{ji}$, the statistics of the effective local fields in the NESS is remarkably simple Mezard2011a : In the thermodynamic limit, $\Theta_{i}$ turns out to follow a Gaussian distribution at least in some regimes. This distribution is characterized by a mean $\theta^{\textrm{MF}}_{i}$, standard deviation $\Delta_{i}$ and covariance $\epsilon_{ij}$. Using the definition of the effective local field (119), these parameters are linked to the spin observables by

[TABLE]

The key idea is that one can transform back and forth between $\{\sigma_{i}(t)\}$ and $\{\Theta_{i}(t)\}$ via

[TABLE]

and evaluate the correlation functions within the Gaussian theory. With $\theta_{i}=\theta^{\textrm{MF}}_{i}+\sqrt{\Delta_{i}}x$ and $\theta_{j}=\theta^{\textrm{MF}}_{j}+\sqrt{\Delta_{j}}y$, where $x$ and $y$ are univariate Gaussian random variables, one obtains from $D_{ik}=\langle\operatorname{th}(\Theta_{i})\sigma_{k}\rangle-\langle\operatorname{th}(\Theta_{i})\rangle\langle\sigma_{k}\rangle$

[TABLE]

In the last step we have used the fact that covariances between spin variables are small Mezard2011a . Inserting the result for the covariance (135) gives again an equation of the same form as (133), however with $A_{im}=\delta_{im}\int\textrm{d}x\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}\left[1-\operatorname{th}^{2}(\theta^{\textrm{MF}}_{i}+\sqrt{\Delta_{i}}x)\right]$. Mean-field theory neglects the fluctuations here and renders this term as $1-m^{2}_{i}$, whereas the fluctuations can be captured more accurately under the Gaussian theory. However, $A_{im}$ cannot be determined directly from the data alone, but also require the couplings; Mezard2011a gives an iterative scheme to infer the parameters of the effective local fields (135) as well as the coupling matrix and magnetic fields. The typical-case performance of the Gaussian theory in the thermodynamic limit has been analysed within the framework of statistical learning bachschmid2015learning , finding the Gaussian theory breaks down at strong couplings and a small number of samples.

The Gaussian distribution of local fields is not limited to the asymmetric Ising model. In fact, the asymmetric Ising model is one particular example from a class of models called generalized linear models. In this model class, the Gaussian approximation has been used in the context of neural network reconstruction toyoizumi2009mean already prior to its application to the asymmetric Ising model.

III.2.4 Method comparison

We compare the results of the mean-field approximation, the Gaussian approximation, and the maximization of the exact likelihood (128). As in section II.3, we draw couplings from a Gaussian distribution with mean zero and standard deviation $\beta/\sqrt{N}$, but now couplings $J^{0}_{ij}$ are statistically independent of $J^{0}_{ji}$, so the matrix of couplings is in general asymmetric. Fields are drawn from a uniform distribution on the interval $[-0.3\beta,0.3\beta]$. We then sample a time series of $M=15000$ steps by parallel updates under Glauber dynamics (125). The scatter plots in the top row of figure 8 compares couplings and fields reconstructed by different methods with the couplings and fields of the original model. It shows that couplings are significantly underestimated by the mean-field reconstruction, a bias which is avoided by the Gaussian theory. The right-hand plot shows the relative reconstruction error (115) against $\beta$. As in the equilibrium case shown in figure 5, reconstruction by any method is limited at small $\beta$ by sampling noise. At strong couplings $\beta$, mean-field theory breaks down. Also at strong couplings, the iterative algorithm for Gaussian approximation reconstruction converges very slowly and stops when the maximum number of iterations is reached (here is set $50000$ steps).

III.3 Outlook: Reconstruction from the steady state

Time series data is not available in all applications. This prompts the question how to reconstruct the parameters of a non-equilibrium model from independent samples of the steady state. It is clear that, unlike for equilibrium systems, pairwise correlations are insufficient to infer couplings: The matrix of correlations is symmetric, whereas the matrix of couplings is asymmetric for non-equilibrium systems. Hence, there are twice as many free parameters as there are observables. Similarly, using an equilibrium model like (1) (maximum-entropy model) on data generated by a non-equilibrium model would give parameters matching the two-spin observables, but entirely different from the true parameters underlying the data.

One solution uses three-spin correlations to infer the couplings dettmernguyenberg . A problem of this approach is that connected three-spin correlations are small since the effective local fields are well described by a multi-variate Gaussian distribution (see III.2.3).

A second approach uses perturbations of the non-equilibrium steady state: We measure one set of pair correlations at certain (unknown) model parameters, and then a second set of pair correlations of a perturbed version of the system. Possible perturbations include changing one or several of the couplings by a known amount, changing the magnetic fields, or fixing particular spins to a constant value. This generates two sets of coupled equations specifying two symmetric correlations matrices, which can be solved for one asymmetric coupling matrix. Conceptually, such an approach is well known in biology, where altering parts of a system and checking the consequences is a standard mode of scientific inquest. Neural stimulation can lead to a rewiring of neural connections, and the effects of this neural plasticity can be tracked in neural recordings rebesco2010rewiring . An exciting development is optogenetic tools, which allow to stimulate and monitor the activity of individual neurons in vivo tischer2014illuminating . In the context of inferring gene regulatory networks from gene expression data, perturbation-based approaches have been used both with linear tegner2003a and non-linear models of gene regulation nelander2008a ; molinelli2013a , see also section I.2. In this context, it is also fruitful to consider the genetic variation occurring in a population of cells as a source of perturbations. Such an approach has already been used in gene regulatory networks rockman2008reverse ; zhu2008integrating , and with the expansion of tools to analyze large numbers of single cells zeisel2015cell ; cadwell2015electrophysiological , it may soon spread to other types of networks.

IV Conclusions

At the end of this overview, we step back and summarize the aims and motivations behind the inverse Ising problem, discuss the efficacy of different approaches, and outline different areas of research that may involve the statistical mechanics of inverse problems in the future.

Motivation. The inverse Ising problem arises in the context of very different types of questions connected with the inference of model parameters. The first and most straightforward question appears when data actually is generated by a process obeying detailed balance and has pairwise couplings between binary spins. The problem of inferring the parameters of such a model can then be phrased in terms of the equilibrium statistics (11), the maximisation of the likelihood (II.1), and the use of approximation schemes discussed in section II.

The second question arises when data is generated by a different, and possibly entirely unknown type of process, and we seek a statistical description of the data in terms of a simpler model matching only particular aspects of the observed data. An example is models with maximum entropy given pairwise correlations (section II.1.3) used to describe neural data in section I.1.

The assumptions behind the Ising model (pairwise couplings, binary spins, $\ldots$) can be relaxed. For instance, multi-valued Potts spin variables have been used extensively in models of biological sequences  Weigt2009a ; santolini2014a . The extension of the Ising model to models with three- and four-spin couplings has not been used yet in the context of inference. Nevertheless, for many of the methods of section II, the extension beyond pairwise couplings would be straightforward. A much larger rift lies between equilibrium and non-equilibrium models. Parameter inference of a non-equilibrium model is often based on data beyond independent samples of the steady state, specifically time-series data.

Methods. The practical question how to infer the couplings and fields that parametrise an Ising model has been answered using many different approximations with different regimes of validity. Section II.3 gives an overview. For data sampled from the equilibrium distribution (1), the pseudolikelihood approach of section II.2 gives a reconstruction close to the optimal one (in fact, asymptotically close in the number of samples). However, this is paid for by a computational effort that scales with the number of samples. For the non-equilibrium regime and a time series of configurations sampled from the stochastic dynamics, the likelihood of a model can be evaluated exactly comparatively easily.

Both the pseudolikelihood and the likelihood of a time series can be understood in the language of regression. Thus a single framework links two of the most successful approaches, in both the equilibrium and non-equilibrium setting. Regression singles out one spin variable as a dependent variable and treats the remainder as independent variables. It then characterizes the statistics of the dependent variable given configurations of the independent variables.

The application of novel concepts to the inverse Ising problem and the development of new algorithms continues. An exciting new direction is interaction screening vuffray2016interaction ; lokhov2016optimal . Vuffray, Misra, Lokhov, and Chertkov introduce the objective function

[TABLE]

to be minimized with respect to the $i^{\mbox{th}}$ column of the coupling matrix and the magnetic field on spin $i$ by convex optimisation. $H_{i}(\boldsymbol{s})$ denotes the part of the Hamiltonian containing all terms in $s_{i}$, see II.2. Sparsity or other properties of the model parameters can be effected by appropriate regularisation terms. This objective function aims to find those parameters which ‘screen the interactions’ (and the magnetic field) in the data, making the sum over samples in (138) as balanced as possible. To illustrate the appeal of this objective function, we look at the infinite sampling limit, where the summation over samples in the objective function (138) can be replaced by a summation over the configurations, reweighted by the Boltzmann distribution with the true couplings $\boldsymbol{J}^{\ast}$ and fields $\boldsymbol{h}^{\ast}$,

[TABLE]

Since the last sum as a function of $J_{ij}$ and $h_{i}$ (for a fixed $i$) is convex and even when being reflected around $J^{\ast}_{ij}$ and $h_{i}^{\ast}$, so is the whole objective function. It follows immediately that this objective function has a unique global minimum at $J_{ij}=J^{\ast}_{ij}$ and $h_{i}=h^{\ast}_{i}$.

Of course, the sample average over $e^{H_{i}(\boldsymbol{s})}$ is affected by sampling fluctuations when the number of samples is small. Nevertheless, the minimum of this objective function nearly saturates the information theoretic bounds on the reconstruction of sparse Ising models of Santhanam and Wainwright santhanam2012information (see section II.1.4).

The flat histogram method in Monte Carlo simulations uses a similar rebalancing with respect to the Boltzmann measure wang2000flat . Also, there may be conceptual links between interaction screening and the fluctuations theorems such as the Jarzynski equality jarzynski1997nonequilibrium , which also take sample averages over exponentials of different thermodynamic quantities.

Another recent development concerning the sparse inverse Ising problem is the use of Bayesian model selection techniques by Bulso, Marsili and Roudi bulsomarsiliroudi2016a .

More to come. Over the last two decades, interest in inverse statistical problems has been driven by technological progress and this progress is likely to continue opening up new applications. Beyond the extrapolation of technological developments, there are several broad areas at the interface of statistical mechanics, statistics, and machine learning where inverse statistical problems such as the inverse Ising problem might play a role in the future.

• •

Stochastic control theory. Stochastic control theory seeks to steer a stochastic system towards certain desired states kappen2007a . An inverse problem arises when the parameters describing the stochastic system are only partially known. As a result, its response to changes in external control parameters (’steering’) must be predicted on the basis of its past dynamics. A recent application is the control of cell populations, specifically populations of cancer cells fischer2015value . Such cell populations evolve stochastically due to random fluctuations of cell duplications and cell deaths. Birth and death rates of at least part of the population can be controlled by therapeutic drugs.

• •

Network inference. Like regulatory connections between genes discussed in section I.2, metabolic and signalling interactions also form intricate networks. An example where such a network optimizes a specific global quantity is flux-balance analysis orth2010flux , which models the flow of metabolites through a metabolic network. Inverting the relationship between metabolic rates and metabolite concentrations allows in principle to infer the metabolic network on the basis of observed metabolite concentrations deMartino2016a .

• •

Causal analysis. The do-calculus developed by Judea Pearl seeks to establish causal relationships behind statistical dependencies pearl1995causal . do-calculus is based on interventions such as fixing a variable to a particular value, and then observing the resulting statistics of other variables. Recently, Aurell and Del Ferrano have found a link between do-calculus and the dynamic cavity method from statistical physics aurell2015causal .

• •

Maximum entropy and dynamics. The notion of a steady-state distribution with maximum entropy has been generalized to a maximum-entropy distribution over trajectories of a dynamical system presse2013a . Like the maximum entropy models of section I.1-I.3, this approach can be used to derive simple effective models of the dynamics of a system, whose parameters can be inferred from time series. In fact, Glauber dynamics with parallel updates gives the maximum entropy distribution of $\boldsymbol{\sigma}(t+1)$ given $\boldsymbol{s}(t)$. So far, applications have been in neural modelling marre2009 ; nassercessac2014a , in the effective dynamics of quantitative traits in genetics bartondevladar2009a ; bovdova2015a , and in flocking dynamics cavagna2014a .

• •

Hidden variables. Even in large-scale data sets, there will be variables that are unobserved. Yet, those hidden variables affect the statistics of other variables, and hence the inference we make about interactions between the observed variables dunn2013learning ; tkacikmarremoraetal2013a ; battistin2015belief ; rouditaylor2015a . This can lead to a signature of critical behaviour, even when the original system is not critical marsilimastromatteoroudi2013a ; schwabnemenmanMehta2014a .

• •

High-dimensional statistics and inference. In many applications, the number of systems parameters to be inferred is of the same order of magnitude or exceeds the sample size. The field of high-dimensional statistics deals with this regime, and a fruitful interaction with statistical physics has emerged over the last decade ZdeborovaKrzakala2015 , spurred by applications such as compressed sensing. Vuffray et al.  vuffray2016interaction propose the objective function (138) for the inverse Ising problem in the high-dimensional regime. Other objective functions might perform even better, and the optimal objective function may depend on the number of samples and the statistics of the underlying couplings.  berg2016 ; bachschmid2017statistical use the statistical mechanics of disordered systems to find the objective function which minimizes the difference between the reconstructed and the underlying couplings for Gaussian distributed couplings.

• •

Restricted Boltzmann machines and deep learning. Recently, (deep) feed-forward neural networks have re-established themselves as powerful learning architectures, leading to spectacular applications in the area of computer vision, speech recognition, data visualisation, and game playing lecun2015deep . This progress has also demonstrated that the most challenging step in data analysis is the extraction of features from unlabelled data. There is a wide range of methods for feature extraction, the most well-known ones being convolutional networks for images lecun1998gradient ; lecun2015deep , or auto-encoders hinton2006a ; vincent2010stacked , and restricted Boltzmann machines (RBMs) hinton2006a for less structured data. RBMs are possibly one of the most general methods, although the algorithms used for finding the optimal parameters are heuristic and approximate in many respects. For instance, in the case of multiple layers, RMBs correspond to generic Boltzmann machines (with feedback loops), for which learning is a hard computational problem. Usually, the sub-optimal solution which is adopted is used to train each layer independently. These observations are not surprising since RBMs are nothing but an inverse Ising problem with a layer of visible spins connected to a layer of hidden spins. Empirical data is available only for the visible spins. The role of the hidden variables is to compress information and identify structural features in the data. Learning consists of finding the visible-to-hidden couplings $J_{ij}=J_{ji}$ and the local fields such that summing over to the hidden variables gives back a (marginalized) probability distribution over the visible variables which is maximally consistent with the data, i.e. has minimal KL divergence (18) from the empirical distribution of the data.

• •

Learning phases of matter. One aspect of the emerging applications of machine learning in quantum physics is the identification of non-trivial quantum states from data. This is of particular interest in phases of matter where the order parameter is either unknown or hard to compute (such as the so-called entanglement entropy islam2015measuring ; levin2006detecting ). Techniques from machine learning, specifically deep learning with neural networks, have recently been used to classify quantum states without knowledge of the underlying Hamiltonian carrasquilla2017machine .

Acknowledgements

Discussion with our colleagues and students have inspired and shaped this work. In particular we would like to thank Erik Aurell, Michael Berry, Andreas Beyer, Simona Cocco, Simon Dettmer, David Gross, Jiang Yijing, David R. Jones, Bert Kappen, Alessia Marruzzo, Matteo Marsili, Marc Mézard, Rémi Monasson, Thierry Mora, Manfred Opper, Andrea Pagnani, Federico Ricci-Tersenghi, Yasser Roudi, Gašper Tkačik, Aleksandra Walczak, Martin Weigt and Pieter Rein ten Wolde. Many thanks to Gašper Tkačik and Michael Berry for making the neural recordings from tkavcik2014a ; tkacik2015 available. This work was supported by the DFG under Grant SFB 680; BMBF under Grants emed:SMOOSE and SYBACOL.