S Cocco

R Monasson

Adaptive cluster expansion for the inverse Ising problem: convergence, algorithm and tests

We present a procedure to solve the inverse Ising problem, that is to find the interactions between a set of binary variables from the measure of their equilibrium correlations. The method consists in constructing and selecting specific clusters of variables, based on their contributions to the crossentropy of the Ising model. Small contributions are discarded to avoid overfitting and to make the computation tractable. The properties of the cluster expansion and its performances on synthetic data are studied. To make the implementation easier we give the pseudo-code of the algorithm.

• A: find the interaction network from a set of spin configurations σ τ . It is generally assumed in the graphical model community that the Ising model is exact, that is, that the underlying distribution of the data is truly an Ising model with unknown interaction parameters J. The question is to find which interactions J ij are non zero (or larger than some J min is absolute value), and how many configurations (value of B) should be sampled to achieve this goal with acceptable probability.

• B: find the interactions J ij and the fields h i from the frequencies p i , p ij only. Those frequencies should be reproduced within a prescribed accuracy, ǫ, not too small (compared to the error on the data) to avoid overfitting. Note that in general the Ising model is not the true underlying model for the data here; it is only the model with maximal entropy given the constraints on 1-and 2-point correlations.

.

. The purpose of the present paper is to discuss its performances and limitations.

It is essential to be aware of the presence of noise in the data, e.g. due to the imperfect sampling (finite number B of configurations). A potential risk is overfitting: the network of interactions we find at the end of the inference process could reproduce the mere noisy data, rather than the 'true' interactions. How can one disentangle noise from signal in the data? A popular approach in the statistics community is to require that the inferred interaction network be sparse. The rationale for imposing sparsity is two-fold. First, physical lattices are very sparse, and connect only close sites in the space; it is possible but not at all obvious that networks modeling other e.g. biological data enjoy a similar property. Secondly, an Ising model with a sparse interaction network reproducing a set of correlations is a sparing representation of the statistics of the data, and, in much the same spirit as the minimal message length approach [14], should be preferred to models with denser networks. The appeal of the approach is largely due to the fact that imposing sparsity is computationally tractable.

The criterion required by our procedure is not that the interaction network should be sparse, but that the inverse Ising problem should be well-conditioned. To illustrate this notion, consider a set of data, i.e. of frequencies p i , p kl , and assume one has found the solution h i , J kl to the corresponding inverse Ising problem. Let us now slightly modify one or a few frequencies, say, p 12 → p ′ 12 = p 12 + δp 12 , and solve again the corresponding inverse Ising problem, with the results h ′ i , J ′ kl . Let δJ kl = J ′ kl -J kl and δh i = h ′ ih i measure the response of the interaction parameters to the small modification of p 12 alone. Two extreme cases are:

I. INTRODUCTION

The Ising model is a paradigm of statistical physics, and has been extensively studied to understand the equilibrium properties and the nature of the phase transitions in various systems in condensed matter [1]. In its usual formulation, the Ising model is defined over a set of N binary variables σ i , with i = 1, 2, . . . , N . The variables, called spins, are submitted to a set of N local fields, h i , and of 1 2 N (N -1) pairwise couplings, J ij . The observables of the model, such as the average values of the spins or of the spin-spin correlations over the Gibbs measure,

σ i , σ k σ l , (1)
are well-defined and can be calculated from the knowledge of those interaction parameters. We will refer to the task of calculating (1) given the interaction parameters as to the direct Ising problem. In many experimental cases, the interaction parameters are unknown, while the values of observables can be estimated from measurements. A natural question is to know if and how the interaction parameters can be deduced from the data ( [2][3][4][START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF][START_REF] Balakrishnan | [END_REF][7]). When the coupling matrix is known a priori to have a specific and simple structure, this question can be answered with an ordinary fit. For instance, in a two-dimensional and uniform ferromagnet, all couplings vanish but between neighbors on the lattice, and J ij = J for contiguous sites i and j. In such a case, the observable such as the average correlation between neighboring spins, c, depends on a single parameter, J. The measurement of c gives a direct access to a value of J. However, data coming from complex systems arising in biology, sociology, finance, ... can generally not be interpreted with such a simple Ising model, and the fit procedure is much more complicated for two reasons. First, in the absence of any prior knowledge about the interaction network, the number of interaction parameters J ij to be inferred scales quadratically with the system size N , and can be very large. Secondly, the quality of the data is a crucial issue. Experimental data are plagued by noise, coming either from the measurement apparatus or from imperfect sampling. The task of fitting a very large number of interaction parameters from 'noisy' data has received much attention in the statistics community, under the name of high-dimensional inference [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF].

To be more specific, the inverse Ising problem is defined as follows. Assume that a set of B configurations σ τ = {σ τ 1 , σ τ 2 , . . . , σ τ N }, with τ = 1, 2, . . . , B are available from measurements. We compute the empirical 1-and 2-point averages through

p i = 1 B B τ =1 σ τ i , p kl = 1 B B τ =1 σ τ k σ τ l . (2)
The inverse Ising problem consists in finding the values of the N local fields, h i , and of the 1 2 N (N -1) interactions, J ij , such that the individual and pairwise frequencies of the spins (1) defined from the Gibbs measure coincide with their empirical counterparts, p i and p kl . While the Gibbs measure corresponding to the Ising model is by no means the unique measure allowing one to reproduce the data p i and p kl , it is the distribution with the largest entropy doing so [START_REF] Jaynes | Proc. IEEE[END_REF]. In other words, the Ising model is the least constrained model capable of matching the empirical values of the 1-and 2-point observables. This property explains the recent surge of interest in defining and solving the inverse Ising problem in the context of the analysis of biological, e.g. neurobiological [2][3][4][START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF] and proteomic [START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF][START_REF] Balakrishnan | [END_REF] data.

As a result of its generality, the inverse Ising problem has been studied in various fields under different names, such as Boltzmann machine learning in learning theory [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF][START_REF] Ackley | [END_REF] or graphical model selection in statistical inference [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF][START_REF] Ravikumar | [END_REF]16]. While the research field is currently very active, the diversity of the tools and, sometimes, of the goals make somewhat difficult to compare the results obtained across the disciplines. Several variants of the inverse Ising problem can be defined:

• Localized response: the response is restricted to the parameters involving spins 1 and 2 only, i.e. δh 1 , δh 2 , δJ 12 = 0; it vanishes for all the other parameters.

• Extended response: the response spreads all over the spin system, and all the quantities δh i , δJ kl are non-zero.

Intermediate cases will generically be encountered, and are symbolized in Fig. 1(a)&(b). For instance, if the response is non-zero over a small number of parameters only, which define a 'neighborhood' of the spins 1, 2, we will consider it is localized. Obviously, the notion of 'smallness' cannot be rigorously defined here, unless the system size N can be made arbitrarily large and sent to infinity. Drawing our inspiration from the vocabulary of numerical analysis, we will say that the inverse Ising problem is well-conditioned if the response is localized. For a well-conditioned problem, a small change of one or a few variables essentially affects one or a few interaction parameters. On the contrary, most if not all interaction parameters of a ill-conditioned inverse Ising problem are affected by an elementary modification of the data. This notion must be distinguished from the concept of ill-posed problem. As we will see in Section II, the inverse Ising problem is always well-posed, once an appropriate regularization is introduced: given the frequencies, there exists a unique set of interaction parameters reproducing those data, regardless of how hard it is to compute. Not all inverse Ising problems are well-conditioned. However, it is our opinion that only those ones should be solved. The reason is that, in generic experimental situations, only a (small) region of the system is accessible. Solving the inverse problem attached to this sub-system makes sense only if the problem is well-conditioned. If it is ill-conditioned, extending even by a bit the sub-system would considerably affect the values of most of the inferred parameters (Fig. 1(c)). Hence, the interaction parameters would be very much dependent on the part of the system which is not measured! Such a possibility simply means that the inverse problem, though mathematically well-posed, is not meaningful. FIG. 1: Schematic representation of a well-conditioned (a) and an ill-conditioned (b) inverse Ising problems. The gray areas symbolize the set of spins (full dots) whose interactions (links) and fields are affected by a change of the frequency p12 of spins 1 and 2. The response is localized in (a) and extended in (b). Experiments usually measure a restricted part of the system (dashed contour) only (c). Increasing the size of the measured sub-system, e.g. by including the frequencies of the extra-variables 3 and 4, will modify most of the inferred interaction parameters if the problem is ill-conditioned.

Interestingly, the response of the interactions to a change of a few correlations can be localized, while the response of the correlations to a change of a few interactions is extended. An example is given by 'critical' Ising models, where correlations extend over the whole system. However, the corresponding inverse Ising problem may be well-conditioned.

The presence of noise in the data considerably affects the status of the inverse Ising problem. As we will see later, even well-conditioned problems in the limit of perfect sampling (B → ∞) become ill-conditioned as soon as sampling is imperfect (finite B). The same statement holds for the sparsity-based criterion mentioned above: when data are generated by a sparse interaction network, the solution to the inverse Ising model is not sparse as a consequence of imperfect sampling. Only the presence of an explicit and additional regularization forces the solution to be sparse. In much the same way, the procedure we present hereafter builds a well-conditioned inverse Ising problem, which prevents overfitting of the noise. This procedure is based on the expansion of the entropy at fixed frequencies in clusters of spins, a notion closely related to the neighborhoods appearing in the localized responses.

The plan of the article is as follows. In Section II we give the notations and precise definitions of the inverse Ising problem, and briefly review some of the resolution procedures in the literature. In Section III, we explain how the entropy can be expanded as a sum of contributions, one for each cluster (or sub-set) of spins, and review the properties of those entropic contributions. The procedure to truncate the expansion and keep only relevant clusters is discussed in Section IV. The pseudo-codes and details necessary for the implementation of the algorithm can be found in Section V. Applications to artificial data are discussed at length in Section VI. Finally, Section VII presents some perspectives and conclusions. To improve the readability of the paper most technical details have been relegated to technical appendices.

II. THE INVERSE ISING PROBLEM: FORMULATIONS AND ISSUES A. Maximum Entropy Principle formulation

We consider a system of N binary variables, σ i = 0, 1, where i = 1, 2, . . . , N . The average values of the variables, p i , and of their correlations, p kl , are measured, for instance through the empirical average over B sampled configurations of the system, see equations (2). As the correlations p kl are obtained from the empirical measure, the problem is realizable [START_REF] Cagliotti | [END_REF]. Let p = {p i , p kl } denote the data. The Maximum Entropy Principle (MEP) [START_REF] Jaynes | Proc. IEEE[END_REF] postulates that the probabilistic model P (σ) should maximize the entropy S of the distribution P under the constraints σ P (σ) = 1 , σ P (σ) σ i = p i , σ P (σ) σ k σ l = p kl .

(

In practice these constraints are enforced by the Lagrange multipliers λ and J = {h i , J kl }. The maximal entropy is

The maximization condition over P shows that the MEP probability corresponds to the Gibbs measure P J of the celebrated Ising model,

P J [σ] = e -HIsing [σ|J] Z[J] (5)
where the energy function is

H Ising [σ|J] = - i h i σ i - k<l J kl σ k σ l (6)
and Z[J] = σ exp(-H Ising [σ|J]) denotes the partition function. The values of the couplings and fields [51] are then found through the minimization of

S Ising [J|p] = log Z[J] - i h i p i - k<l J kl p kl . (7)
over J. The minimal value of S Ising coincides with S defined in (4). The cross-entropy S Ising has a simple interpretation in terms of the Kullback-Leibler divergence between the Ising distribution P J [σ] and the empirical measure over the observed configurations, P obs [σ]. Assume B configurations of the N variables, σ τ , with τ = 1, 2, . . . , B, are sampled. We define the empirical distribution through

P obs [σ] = 1 B B τ =1 δ σ,σ τ , (8)
where δ denotes the N -dimensional Kronecker delta function. It is easy to check from (7) that

S Ising [J|p] = - σ P obs [σ] log P J [σ] = - σ P obs [σ] log P obs [σ] + D P obs ||P J , (9)
where D denotes the KL-divergence. Hence, the minimization procedure over J ensures that the 'best' Ising measure (as close as possible to the empirical measure) is found.

B. Regularization and Bayesian formulation

We consider the Hessian of the cross-entropy S Ising , also called Fisher information matrix, which is a matrix of dimension 1 2 N (N + 1), defined through

χ = ∂ 2 S Ising ∂J∂J = χ i,i ′ χ i,k ′ l ′ χ kl,i ′ χ kl,k ′ l ′ . ( 10
)
The entries of χ are obtained upon repeated differentiations of the partition function Z[J], and can be expressed in terms of averages over the Ising Gibbs measure • J ,

χ i,i ′ = σ i σ i ′ J -σ i J σ i ′ J , (11)
χ i,k ′ l ′ = σ i σ k ′ σ l ′ J -σ i J σ k ′ σ l ′ J , χ kl,k ′ l ′ = σ k σ l σ k ′ σ l ′ J -σ k σ l J σ k ′ σ l ′ J .
Consider now an arbitrary 1 2 N (N + 1)-dimensional vector x = {x i , x kl }. The quadratic form

x † • χ • x = i x i σ i -σ i J + k<l x kl σ k σ l -σ k σ l J 2 J ( 12
)
is semi-definite positive. Hence, S Ising is a convex function. However the minimum is not guaranteed to be unique if χ has zero modes, nor to be finite. To circumvent those difficulties, one can 'regularize' the cross-entropy S Ising by adding a quadratic term in the interaction parameters, which forces χ to become definite positive, and ensures the uniqueness and finiteness of the minimum of S Ising . In many applications, no regularization is needed for the fields h i . The reason can be understood intuitively as follows. Consider a data set where all variables are independent, with small but strictly positive means p i . Then, the empirical average products, p kl , may vanish if the number B of sampled configurations is not much larger than (p k p l ) -1 . This condition is often violated in practical applications, e.g. the analysis of neurobiological or protein data [2,[START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF][START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF]. Hence, poor sampling may produce infinite negative couplings. We therefore add the following regularization term to S Ising ,

γ k<l J 2 kl p k (1 -p k )p l (1 -p l ) . (13)
The precise expression of the regularization term is somewhat arbitrary, and is a matter of convenience. The dependence on the p i 's in (13) will be explained in Section III B. Other regularization schemes, based on the L 1 norm rather than on the L 2 norm are possible, such as

γ k<l |J kl | p k (1 -p k )p l (1 -p l ) . (14)
The above regularization is especially popular among the graphical model selection community [START_REF] Ravikumar | [END_REF], and favors sparse coupling networks, i.e. with many zero interactions. The introduction of a regularization is natural in the context of Bayesian inference. The Gibbs probability P J [σ] defines the likelihood of a configuration σ. The likelihood of a set of B independently drawn configurations σ τ is given by the product of the likelihoods of each configuration. The posterior probability of the parameters (fields and couplings) J given the configurations σ τ , τ = 1, 2, . . . , B, is, according to Bayes' rule,

P post [J|{σ τ }] ∝ B τ =1 P J [σ τ ] P 0 [J] , (15)
up to an irrelevant J-independent multiplicative factor. In the equation above, P 0 is a prior probability over the couplings and fields, encoding the knowledge about their values in the absence of any data. Taking the logarithm of (15), we obtain, up to an additive J-independent constant,

log P post [J|{σ τ }] = -B S Ising [J|p] + log P 0 [J] . (16)
Hence, the most likely value for the parameters J is the one minimizing S Ising [J|p] -1 B log P 0 [J]. The regularization terms ( 13) and ( 14) then correspond to, respectively, Gaussian and exponential priors over the parameters. In addition, as the prior is independent of the number B of configurations, we expect the strength γ to scale as 1 B . The optimal value of γ can be also determined based on Bayesian criteria [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF][START_REF] Mackay | [END_REF] (Appendix A).

We emphasize that the Bayesian framework changes the scope of the inference. While the MEP aims to reproduce the data, the presence of a regularization term leads to a compromise between two different objectives: finding an Ising model whose observables (one-and two-point functions) are close to the empirical values and ensuring that the interaction parameters J have a large prior probability P 0 . In other words, a compromise is sought between the faithfulness to the data and the prior knowledge about the solution. The latter is especially important in the case of poor sampling (small value of B or data corrupted by noise). For instance, the regularization term based on the L 1 -norm [START_REF] Wallace | Statistical and Inductive Inference by Minimum Message Length[END_REF] generally produces more couplings equal to zero than its L 2 -norm counterpart [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. This property is desirable if one a priori knows that the interaction graph is sparse. Hence, the introduction of a regularization term can be interpreted as an attempt to approximately solve the inverse Ising problem while fulfilling an important constraint about the structure of the solution. We will discuss the nature of the structural constraints corresponding to our adaptive cluster algorithm in Section IV C.

Knowledge of the inverse of the Fisher information matrix, χ -1 , allows for the computation of the statistical fluctuations of the inferred fields and couplings due to a finite number B of sampled configurations. According to the asymptotic theory of inference, the posterior probability P post [J|{σ τ }] over the fields and couplings becomes, as B gets very large, a normal law centered in the minimum of S Ising [J|p]. The covariance matrix of this normal law is simply given by 1 B χ -1 . Consequently the standard deviations of the fields h i and of the couplings J kl are, respectively,

δh i = 1 B (χ -1 ) i,i , δJ kl = 1 B (χ -1 ) kl,kl . (17)
In order to remove the zero modes of χ and have a well-defined inverse matrix χ -1 , the Ising model entropy S Ising (7) can be added a regularization term, e.g. [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF], which guarantees that χ is positively defined. The Fisher information matrix, χ, can also be used to estimate the statistical deviations of the observables coming from the finite sampling. If the data were generated by an Ising model with parameters J, we would expect, again in the large B setting, that the frequencies p i , p kl would be normally distributed with a covariance matrix equal to 1 B χ. Hence, the typical uncertainties over the 1-and 2-point frequencies are given by

δp i = 1 B χ i,i = σ i J (1 -σ i J ) B , δp kl = 1 B χ kl,kl = σ k σ l J (1 -σ k σ l J ) B . (18)
In practice, we can replace the Gibbs averages above with the empirical averages p i and p kl to obtain estimates for the expected deviations. These estimates will be used to decide whether the inference procedure is reliable, or leads to an overfitting of the data in Section VI.

C. Methods

The inverse Ising problem has been studied in statistics, under the name of graphical model selection, in the machine learning community under the name of (inverse) Boltzmann machine learning, and in the statistical physics literature. Different methods have been developed, with various applications. Some of the methods are briefly discussed below.

A direct calculation of the partition function Z[J] generally requires a time growing exponentially with the number N of variables, and is not feasible when N exceeds a few tens. Inference procedures therefore tend to avoid the computation of Z[J]:

• A popular algorithm is the Boltzmann learning procedure, where the fields and couplings are iteratively updated until the averages σ i J 's and σ k σ l J 's, calculated from Monte Carlo simulations, match the imposed values [START_REF] Ackley | [END_REF]. The number of updatings can be very large in the absence of a good initial guess for the parameters J. Furthermore, for each set of parameters, thermalization may require prohibitive computational efforts for large system sizes N , and problems with more than a few tens of spins can hardly be tackled. Finally, learning data exactly leads to overfitting in the case of poor sampling.

• the Pseudo-Likelihood-based algorithm by Ravikumar et al. [START_REF] Ravikumar | [END_REF]16] is an extension to the binary variable case of Meinshausen and Bühlmann's algorithm [20] and is related to a renormalisation approach introduced by Swendsen [19]. The procedure requires the complete knowledge of the configurations {σ τ } (and not only of the one-and two-point functions p). The starting point is given by well-known Callen's identities for the Ising model,

σ i J = 1 1 + exp -j J ij σ j -h i J ≃ 1 B B τ =1 1 1 + exp -j J ij σ τ j -h i (19)
where the last approximation consists in replacing the Gibbs average with the empirical average over the sampled configurations. Imposing that the Gibbs average σ i J coincides with p i is equivalent to minimizing the following pseudo-likelihood over the field h i ,

S i,P L [h i , {J ij , j = i}] = 1 B B τ =1 log 1 + exp j J ij σ τ j + h i -h i p i - j( =i) J ij p ij . (20)
The minimization equations over the couplings J ij , with j = i (and fixed i), correspond to Callen identities for two-point functions. Informally speaking, the pseudo-likelihood approach simplifies the original N -body problem into N independent 1-body problem, each one in a bath of N -1 quenched variables. Note that the couplings J ij and J ji (found by minimizing S j,P L ) will generally not be equal. However, as far as graphical model selection is concerned, what matters is whether J ij and J ji are both different from zero.

The pseudo-entropy S i,P L is convex, and can be minimized after addition of a L 1 -norm regularization term [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF][START_REF] Ravikumar | [END_REF]21]. The procedure is guaranteed to find strong enough couplings [START_REF]The minimal strength of the couplings which can be 'detected' depend on the quality of the sampling, and scales as log N/B[END_REF] in a polynomial time in N , provided that the data were generated by an Ising model (which is usually not the case in practical applications) and that a quantity closely related to the susceptibility χ (10) is small enough. The latter condition holds for weak couplings and may break down for strong couplings [16]. For a review of the literature in the statistics community, see [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF].

In specific cases, however, the partition function can be obtained in polynomial time. Two tractable examples are:

• Mean-field models, which are characterized by dense but weak interactions. An example is the Sherrington-Kirkpatrick model where every pair of spins interact through couplings of the order of N -1/2 [27]. The entropy S[p] coincides asymptotically with

S MF (p) = 1 2 log det M (p), where M ij (p) = p ij -p i p j p i (1 -p i )p j (1 -p j ) , (21)
which can be calculated in O(N 3 ) time [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF]30]. Expression (21) has been obtained from the high temperature expansion [23][24][25] of the Legendre transform of the free energy , and is consistent with the so-called TAP equations [26]. The derivative of S MF with respect to p gives the value of the couplings and the fields,

(J MF ) kl = - ∂S MF ∂p kl = - (M -1 ) kl p k (1 -p k )p l (1 -p l ) , (h MF ) i = - ∂S MF ∂p i = j( =i) (J MF ) ij c ij p i -1 2 p i (1 -p i ) -p j , (22)
where c ij = p ijp i p j is the connected correlation. From a practical point of view, expression ( 21) is a good approximation for solving the inverse Ising problem [28][29][30] on dense and weak interaction networks , but fails to reproduce dilute graphs with strong interactions.

• Ising models on tree-like structures, i.e. with no or few interaction loops. Message passing methods are guaranteed to solve the associated inverse Ising problems. For trees, the partition functions can be calculated in a time linear in N . Sparse networks of strong interactions with long-range loops, such as Erdös-Renyi random graphs, can also be successfully treated in polynomial time by message-passing procedures [START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF]31,32]. However, these methods generally break down in the presence of strongly interacting groups (clusters) of spins.

When an exact calculation of the partition function is out-of-reach, accurate estimates can be obtained through cluster expansions. Expansions have a rich history in statistical mechanics, e.g. the virial expansion in the theory of liquids [33,[START_REF] Hansen | McDonald Theory of Simple Liquids[END_REF]. However, cluster expansions suffer from several drawbacks. First, in cluster variational methods [31,[START_REF]The Equilibrium Theory of Classical Fluids (A lecture note and reprint[END_REF], the calculation of the contributions coming from each cluster generally involves the resolution of non trivial and self-consistent equations for the local fields, which seriously limits the maximal size of clusters considered in the expansion. Secondly, the composition and the size of the clusters is usually fixed a priori, and does not adapt to the specificity of the data [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF]. The combinatorial growth of the number of clusters with their size entails strong limits upon the maximal sizes of the network, N , and of the clusters, K. Last of all, cluster expansions generally ignore the issue of overfitting.

Recently, we have proposed a new cluster expansion, where clusters are built recursively, and are selected or discarded, according to their contribution to the cross-entropy S [8]. This selection procedure allows us to fully account for the complex interaction patterns present in experimental systems, while preventing a blow-up of the computational time. The purpose of this paper is to illustrate this method and discuss its advantages and limitations.

III. CLUSTER EXPANSION OF THE CROSS-ENTROPY

A. Principle of the expansion

In this Section, we propose a cluster expansion for the entropy S(p). A cluster, Γ, is defined here as a non-empty subset of (1, 2, . . . , N ). To illustrate how the expansion is built we start with the simple cases of systems with a few variables (N = 1, 2), in the absence of the regularization term [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. Consider first the case of a single variable, N = 1, with average value p 1 . The entropy S(p 1 ) can be easily computed according to the definitions given in Section II, with the result

S (1) (p 1 ) = min h1 S Ising [h 1 |p 1 ] = min h1 log 1 + e h1 -h 1 p 1 = -p 1 log p 1 -(1 -p 1 ) log(1 -p 1 ) . ( 23
)
We recognize the well-known expression for the entropy of a 0-1 variable with mean value p 1 . For reasons which will be obvious in the next paragraph, we will hereafter use the notation ∆S (1) (p 1 ) to denote the same quantity as S (1) (p 1 ). The subscript (1) refers to the index of the (unique) variable in the system. Consider next a system with two variables, with mean values p 1 , p 2 and two-point average p 12 . The entropy S (1,2) (p 1 , p 2 , p 12 ) can be explicitly computed:

S (1,2) (p 1 , p 2 , p 12 ) = min h1,h2,J12 S Ising [h 1 , h 2 , J 12 |p 1 , p 2 , p 12 ] = min h1,h2,J12 log 1 + e h1 + e h2 + e h1+h2+J12 -h 1 p 1 -h 2 p 2 -J 12 p 12 = -(p 1 -p 12 ) log (p 1 -p 12 ) -(p 2 -p 12 ) log(p 2 -p 12 ) -p 12 log p 12 -(1 -p 1 -p 2 + p 12 ) log (1 -p 1 -p 2 + p 12 ) . ( 24
)
We now define the entropy ∆S (1,2) of the cluster of the two variables 1, 2 as the difference between the entropy S (1,2) (p 1 , p 2 , p 12 ) calculated above and the two single-variable contributions ∆S (1) (p 1 ) and ∆S (1) (p 2 ) coming from the variables 1 and 2 taken separately:

∆S (1,2) (p 1 , p 2 , p 12 ) = S (1,2) (p 1 , p 2 , p 12 ) -∆S (1) (p 1 ) -∆S (2) (p 2 ) . (25)
In other words, ∆S (1,2) measures the loss of entropy between the system of two isolated variables, constrained to have means equal to, respectively, p 1 and p 2 , and the same system when, in addition, the average product of the variables is constrained to take value p 12 . Using expressions (23) and (24), we find

∆S (1,2) (p 1 , p 2 , p 12 ) = -(p 1 -p 12 ) log p 1 -p 12 p 1 -p 1 p 2 -(p 2 -p 12 ) log p 2 -p 12 p 2 -p 1 p 2 -p 12 log p 12 p 1 p 2 -(1 -p 1 -p 2 + p 12 ) log 1 -p 1 -p 2 + p 12 1 -p 1 -p 2 + p 1 p 2 . ( 26
)
The entropy of the cluster (1, 2) is therefore equal to the Kullback-Leibler divergence between the true distribution of probability over the two spins and the one corresponding to two independent spins with averages p 1 and p 2 . It vanishes for p 12 = p 1 p 2 . Formula (25) can be generalized to define the entropies of clusters with larger sizes N ≥ 3. Again p = {p i , p kl } denotes the data. For any non-empty subset Γ including 1 ≤ K ≤ N variables, we define two entropies:

• the subset-entropy S Γ (p), which is the entropy of the subset of the K variables for fixed data. It is defined as the right hand side of (4), when the variable indices, i, k, l are restricted to Γ. Note that, when the subset Γ includes all N variables, S Γ (p) coincides with S(p).

• the cluster-entropy ∆S Γ (p), which is the remaining contribution to the subset-entropy S Γ (p), once all other cluster-entropies of smaller clusters have been substracted. The cluster entropies are then implicitly defined through the identity

S Γ (p) = Γ ′ ⊂Γ ∆S Γ ′ (p) , (27)
where the sums runs over all 2 K -1 non-empty clusters Γ ′ of variables in Γ.

Identity (27) states that the entropy of a system (for fixed data) is equal to the sum of the entropies of all its clusters. Figure 2 sketches the cluster decomposition of the entropy for a system of N = 4 variables. For Γ = (1), equation (27) simply expresses that S (1) (p 1 ) = ∆S (1) (p 1 ). For Γ = (1, 2), equation (27) coincides with (25). For Γ = (1, 2, 3), we obtain the definition of the entropy of a cluster made of a triplet of variables:

∆S (1,2,3) (p 1 , p 2 , p 3 , p 12 , p 13 , p 23 ) = S (1,2,3 )(p 1 , p 2 , p 3 , p 12 , p 13 , p 23 ) -∆S (1) (p 1 ) -∆S (2) (p 2 ) -∆S (3) (p 3 ) -∆S (1,2) (p 1 , p 2 , p 12 ) -∆S (1,3) (p 1 , p 3 , p 13 ) -∆S (2,3) (p 2 , p 3 , p 23 ) . (28)
The analytical expression of the cluster-entropy ∆S (1,2,3) is given in Appendix B. The examples above illustrate three general properties of cluster-entropies:

• the entropy of the cluster Γ, ∆S Γ , depends only on the frequencies p i , p ij of the variables i, j in the cluster Γ (and not on all the data in p).

• the entropy of a cluster with, say, K variables, can be recursively calculated from the knowledge of the subsetentropies S Γ ′ (p) of all the subsets Γ ′ ∈ Γ with K ′ ≤ K variables. According to Möbius inversion formula,

∆S Γ (p) = Γ ′ ⊂Γ (-1) K ′ -K S Γ ′ (p) . (29)
• the sum of the entropies of all 2 N -1 clusters of a system of N spins is the exact entropy of the system, see (27) with Γ = (1, 2, . . . , N ).

In practice, to calculate S(p), one first computes the partition function Z[J] by summing over the 2 K configurations σ and, then, minimizes S Ising [J|p] (7) over the interaction parameters J. The minimization of a convex function of 1 2 K(K + 1) variables can be done in time growing polynomially with K. Moreover the addition of the regularization term [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF] can be easily handled. The limiting step is therefore the calculation of Z, which can be done exactly for clusters with less than, say, K = 20 spins.

Hence, only a small number of the 2 N -1 terms in ( 27) can be calculated. In the present work we claim that, in a wide set of circumstances, a good approximation to the entropy S(p) can be already obtained from the contributions of well-chosen clusters of small sizes,

S(p) ≃ Γ∈L ∆S Γ (p) , (30)
We will explain in Section IV how the list of selected clusters, L, is established.

B. The reference entropy S0

So far we have explained how the entropy S(p) can be expanded as a sum of contributions ∆S Γ (p) attached to the clusters Γ. In this Section we present the expansion against a reference entropy, S 0 (p), and two possible choices for the reference entropy.

The idea underlying the introduction of a reference entropy is the following. Assume one can calculate a (rough) approximation S 0 (p) to the true entropy S(p). Then, the difference S(p) -S 0 (p) is smaller than S(p), and it makes sense to expand the former rather than the latter. We expect, indeed, the cluster-entropies to be smaller when the reference entropy S 0 (p) is substracted from the true entropy. We substitute the original definition (27) with the new definition

S(p) = S 0 (p) + Γ⊂(1,2,...,N ) ∆S Γ (p) . (31)
With this new definition, the values of the cluster-entropies ∆S Γ depend on the choice of S 0 ; the previous definition (27) is found back when S 0 = 0. The procedure for the calculation of the cluster-entropies ∆S Γ (p) is the same as in Section III A, upon replacement of S(p) with S(p) -S 0 (p). The three properties of the cluster expansion listed above still hold.

Our final estimate for the entropy will be, compare to (30),

S(p) ≃ S 0 (p) + Γ∈L ∆S Γ (p) . (32)
Hence, the cluster expansion is a way to calculate a correction to the approximation S 0 to the true entropy S.

Obviously, the introduction of a reference entropy is useful in practice only if S 0 (p) can be quickly calculated for the entire system of size N . In other words, the computational effort required to obtain S 0 should scale only polynomially with N . A natural choice for the reference entropy is S 0 = S MF (21), the mean-field entropy discussed in Section II C. As the calculation of S MF requires the one of the determinant of the matrix M (p), it can be performed in a time scaling as N 3 only. In addition, we expect S MF to be a sensible approximation to S for systems with rather weak interactions. Corrections coming from the strongest interactions will be taken care of by the cluster expansion. Regularized versions of the Mean Field entropy can be derived as follows. First, we use the MF expression for the cross-entropy at fixed couplings J kl and frequencies p i , see (7) and [24], to rewrite

S Ising ({p i }, {J kl }) = - 1 2 log det Id -J ′ - k<l J kl (p kl -p k p l ) , where J ′ kl = J kl p k (1 -p k )p l (1 -p l ) , (33)
and Id denotes the N -dimensional identity matrix. We consider the L 2 -norm regularization [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. The entropy at fixed data p is

S 0 (p) = min {J kl } S Ising ({p i }, {J kl }) + γ k<l J 2 kl p k (1 -p k )p l (1 -p l ) = min {J ′ kl } - 1 2 log det Id -J ′ - 1 2 Trace J ′ • M (p) + γ 2 Trace J ′ ) 2 , ( 34
)
where M (p) is defined in (21). The optimal interaction matrix J ′ is the root of the equation

Id -J ′ -1 -M (p) + γ J ′ = 0 . (35)
Hence, J ′ has the same eigenvectors as M (p), a consequence of the dependence on p i we have chosen for the quadratic regularization term in [START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF]. Let j q denote its q th eigenvalue, and mq = (1j q ) -1 . Then,

S 0 (p, γ) = 1 2 N q=1 log mq + 1 -mq , ( 36
)
where mq is the largest root of m2 qmq (m qγ) = γ, and m q is the q th eigenvalue of M (p). Note that mq = m q when γ = 0, as expected.

C. Properties of the cluster entropies ∆SΓ

Diagrammatic expansion in powers of the connected correlations

A better understanding of the cluster expansion and of the role of the reference entropy S 0 can be gained through the diagrammatic expansion of the entropy S(p) in powers of the connected correlations (high-temperature expansion), Note that the entry M ij of the matrix M defined in (21) vanishes linearly with c ij . Thus, an expansion in powers of c ij is equivalent to an expansion in powers of M ij . A procedure to derive in a systematic way the diagrammatic expansion of S(p) is proposed in [30]. The diagrammatic expansion provides a simple representation of the cluster-entropies, in which the entropy S(p) can be represented as a sum of connected diagrams (Fig. 3). Each diagram is made of sites, connected or not by one or more edges. Each point symbolizes a variable, and carries a factor p i . The presence of n(≥ 0) edges between the sites k and l results in a multiplicative factor (c kl ) n . The contribution of a diagram to the entropy is the product of the previous factors, times a function of the p i specific to the topology of the diagram, see [30]. Diagrams of interest include (Fig. 3):

c ij = p ij -p i p j . (37)
• the N single-point diagrams, whose contributions are ∆S (i) (p i );

• the 'loop' diagrams, which consist of a circuit with K edges going through K sites i 1 → i 2 → . . . → i K → i 1 , whose contributions to the entropy are

S loop (p|i 1 , i 2 , . . . , i K ) = (-1) K-1 M i1,i2 M i2,i3 . . . M iK-1,iK M iK ,i1 ; (38)
• the Eulerian circuit diagrams, for which there exists a closed path visiting each edge exactly once;

• the non-Eulerian diagrams, with the lowest number of links (smallest power in M ).

The entropy for two variables i, j, S(p i , p j , p ij ) (24), is the sum of the two single-point diagrams i and j, plus the sum of all connected diagrams made of the two sites i and j with an arbirtrary large number of edges (n ≥ 2) in between (first two columns in Fig. 3). According to (25), the cluster-entropy ∆S (i,j) (p i , p j , p ij ) is equal to the latter sum (second column in Fig. 3). More generally, the entropy of a cluster ∆S Γ (p) is the infinite sum of all diagrams whose sites are the indices in Γ. We now interpret the Mean Field expression for the entropy, S MF , in the diagrammatic framework. We start from identity (21), and rewrite, Using the fact that the diagonal elements of M are equal to unity, the term corresponding to K = 1 above vanishes. For K ≥ 2, we have

S MF (p) = 1 2 Trace log M = 1 2 Trace log Id -Id -M = K≥1 -Trace Id -M K 2K . (39)
-Trace M -Id K = - i1,i2,...,iK δ i1,i2 -M i1,i2 δ i2,i3 -M i2,i3 . . . δ iK-1,iK -M iK-1,iK δ iK ,i1 -M iK ,i1 = i1,i2,...,iK (-1) K-1 Mi1,i2 Mi2,i3 . . . MiK-1,iK MiK,i1 , (40)
where the matrix M has the same off-diagonal elements as M , and has zero diagonal elements. Each term in the above sum corresponds to an Eulerian circuit over K ′ ≤ K sites, where K ′ is the number of distinct indices in (i 1 , i 2 , . . . , i K ). Note that the same circuit can be obtained from different K-uplets of indices. Consider for instance the longest circuits, obtained for K ′ = K, i.e. all distinct indices. 2K different K-uplets (i 1 , i 2 , . . . , i K ) correspond to the same circuit, as neither the starting site nor the orientation of the loop matter. For instance,

i 1 → i 2 → i 3 → i 1 , i 2 → i 3 → i 1 → i 2 , i 1 → i 3 → i 2 → i 1 , .
.. are all equivalent. This multiplicity factor 2K precisely cancels the 2K at the denominator in [START_REF] Mackay | [END_REF]. The contribution corresponding to a circuit therefore coincides with expression [START_REF] Gumbel | Statistics of Extremes[END_REF] for the loop entropy. We conclude that

• S MF (p) sums up all loop diagrams exactly;

• S MF (p), in addition, sums up Eulerian circuit diagrams, but with weights a priori different from their values in the cross-entropy S(p) [53]. An exception is the three-variable Eulerian diagram shown in Fig. 3, whose weights in S MF and S coincide.

• no non-Eulerian diagram is taken into account in S MF (p).

As a conclusion, the diagrammatic expansion provides a natural justification for the choice of the reference entropy S 0 (p) = S MF (p). In addition, it provides us with the dominant contribution to the cluster-entropies once the Mean-Field entropy is substracted, see Fig. 3. A detailed study of those dominant contributions is presented in Appendix C.

Dependence on the cluster size and on the interaction path length

To reach a better understanding of what the cluster-entropy means, we consider the case of finite-dimensional Ising model, e.g. with coupling J > 0 between nearest-neighbors on a D-dimensional lattice. We call ξ the correlation length: the connected correlation c between two sites at large distance d decays as ∼ exp(-d/ξ). We want to characterize the behavior of the cluster-entropy ∆S Γ when the K sites in the cluster Γ are far apart on the lattice. We first choose no reference entropy (S 0 = 0). According to the above diagrammatic expansion, the lowest order diagram (in powers of c) with K sites has the loop topology. We look for the shortest closed path joining all the sites in Γ; let L(Γ) be this contour length, that is, the sum of the distances between neighboring sites along the path (Fig. 4). Then, according to [START_REF] Gumbel | Statistics of Extremes[END_REF], the largest contribution (in absolute value) to the cluster entropy is

∆S Γ ≃ A(p, K) (-1) K-1 exp -L(Γ)/ξ , ( 41
)
where A is a positive function K and of p, the representative value of the frequencies p i of the variables in Γ. We conclude that the sign of the cluster-entropy depends on the parity of the number of sites. Furthermore, ∆S Γ decreases exponentially fast (in absolute value) with the length of the shortest path joining the sites in the cluster. As soon as one site is very far away from the remaining K -1 ones, the cluster-entropy is small. As a consequence, the sum (27) is alternate, and we expect cancellation between contributions coming from clusters sharing the same shortest path, but with different sizes. This crucial point is perfectly illustrated by the onedimensional Ising model. The correlation between two sites at distance

d ij = j -i is, in one dimension, c ij = p i (1 -p i )p j (1 -p j ) exp(-d ij /ξ) (Appendix F).
The matrix M defined in ( 21) has elements

M ij = e -dij/ξ (42)
Then, according to [START_REF] Gumbel | Statistics of Extremes[END_REF], the largest contribution (in absolute value) to the cluster entropy of a cluster containing the K spins i 1 < i 2 < . . . < i K is given by ( 41) with

L Γ = (i 1 , i 2 , . . . , i k ) = 2(i k -i 1 ) , (43)
and

A(p, K) = 1 2 . An exact calculation, reported in Appendix F, shows that ∆S (i1,i2,...,iK ) = (-1) K F exp - i k -i 1 ξ , ( 44
)
where F is a smooth function given in (F11), such that

F (0) = F ′ (0) = 0, F ′′ (0) = -1.
This identity is in agreement with [START_REF] Johnstone | Proc. ICM 2006[END_REF], since the shortest path encircling all sites has length L(Γ) = 2(i Ki 1 ). Hence, all clusters sharing the same 'extremities, i.e. the same values of i 1 and i K , have the same entropies in absolute value. The sign is determined by the parity of K as mentioned above.

Let i K -i 1 ≡ d. Γ = (i 1 , i K ) is the unique cluster of size K = 2 having its 'extremities' equal to i 1 and i K ; its entropy is ∆S * (i1,iK ) = F (exp(-d/ξ)).
There is (d -1) clusters of size K = 3 with the same extremities, each having an entropy equal to -∆S * (i1,iK ) . More generally, there are d-1 K-2 clusters of size K with the same extremities, each having an entropy equal to (-1) K-2 ∆S * (i1,iK ) . The total contribution to the entropy of all those clusters (at fixed extremities i

1 , i k ) is ∆S fixed i1,i k = d+1 K=2 (-1) K-2 d -1 K -2 ∆S * (i1,iK ) = (1 -1) d-1 ∆S * (i1,iK ) = ∆S * (i1,iK ) if d = 1 , 0 if d ≥ 2 . (45)
The above calculation nicely exemplifies the cancellation of cluster-entropies. The contributions of all clusters sharing the same extremities exactly compensate each other, unless those extremities are nearest-neighbors on the lattice. We show in Appendix F that this exact cancellation is a consequence of the existence of a unique interaction path along the unidimensional chain. As a result, in dimension D = 1, the cross-entropy S is simply the sum of the entropies of the clusters made of nearest neighbours.

In the presence of a reference entropy, S 0 = S MF , the asymptotic scaling of the cluster-entropy with its contour length L changes, as the dominant contribution coming from loop diagrams is removed from the cluster expansion and absorbed into S 0 . The subleading contribution to the cluster-entropies is depicted in bold in Fig. 3 and derived in Appendix C. In dimension D = 1, formula ( 41) is replaced with

∆S (i1,i2,...,iK ) = A ′ (p, K) (-1) K-1 exp -3(i k -i 1 )/ξ . (46)
Note the sharper asymptotics decay with the distance between the extremities of Γ than in the absence of reference entropy. As expected, the terms in the expansion of S -S 0 are smaller than the one in the expansion of S alone.

Remarkably, the exact cancellation property studied above also holds when the reference entropy is non-zero, as proven in Appendix F.

In dimension D = 2 or higher, more than one interaction path connect any two spins, and cluster-entropies with the same contour path do not cancel exactly as in the D = 1 case. However, partial cancellations are present. Figure 5 shows the values of the cluster-entropies versus the length of the shortest path, L(Γ), for a small bidimensional 3×3 grid. For such a small system all data p and cluster-entropies ∆S Γ (with up to K = 9 spins) can be calculated by exact enumeration methods. We observe that:

• |∆S Γ | is sensitive to the value of L Γ more than to the size K of the cluster;

• |∆S Γ | rapidly decreases with L Γ ;
• the values of the cluster-entropies reflect the structural properties of the lattice, e.g. clusters made of central sites, such as 4-5, have a larger entropy than the clusters including pairs of edge spins, such as 1-2; [START_REF]0,1, the couplings (Jij ) and fields (hi) are[END_REF] are chosen to make the system critical in the infinite grid size limit, see Section VI C.

= -1 2 j( =i) Jij
• the sign of ∆S Γ changes with the parity of the size of the cluster.

As a result, the contributions to the entropies coming from the clusters sharing the same path, of length L, partially cancel each other. Consider for example the path 1-2-4-5 of length L = 4; all 7 clusters that share this path have similar |∆S|, ranging between 0.0024 and 0.0046, and so does their sum, ∆S (2,4)

+ ∆S (1,5) + ∆S (2,4,5) + ∆S (1,4,5) + ∆S (1,2,5) + ∆S (1,2,4) + ∆S (1,2,4,5) = -0.00242 [54].
The sum of the entropies of the clusters sharing the same path is generally of the same order of magnitude as, or even smaller than the single contributions. Figure 6 shows that the sum of the 12 clusters of contour length L = 2 and of the 4 square-path contributions (|∆S| ≥ .0024) approximates the entropy within 10 -6 .

IV. TRUNCATION OF THE CLUSTER EXPANSION

In this Section we present a truncation scheme for the cluster expansion, which consists in discarding all clusters with entropies smaller than a threshold Θ. We explain why this scheme is efficient, in particular in the presence of sampling noise, and robust against strong correlations in the data (large correlation length). The behavior of the expansion as a function of the threshold is discussed.

A. Schemes for truncating the expansion in the noiseless case

Expansion (27) for S(p) includes 2 N -1 terms, and is useless unless an accurate truncation scheme is available. A naive truncation consists in keeping the contributions from the clusters with ≤ K spins, where K is an arbitrary size. This procedure was applied to neurobiological data (with N ≤ 40, K = 7) in [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF], which are characterized by large negative fields. However it suffers from two drawbacks. First, the combinatorial growth of the number of clusters with N and K impedes its application to very large systems. Secondly, the truncation does not converge properly with increasing K if the correlation length of the system is large.

As an illustration, consider again the 1D-ferromagnetic Ising model, with correlation length ξ. The sign of ∆S Γ alternates with the parity of the size K of Γ; its modulus decays asymptotically as exp(-Ωd/ξ), where d is the maximal distance between any two spins in Γ (Section III C 2), and Ω = 2 if there is no reference entropy (S 0 = 0), Ω = 3 if S 0 = S MF . Let ∆S(K) be the sum of ∆S Γ over all the clusters Γ with K spins. In the thermodynamic limit

(N → ∞), 1 N ∆S(K) ∼ (-1) K-1 d≥K-1 d -1 K -2 exp -Ω d/ξ = (-1) K-1 exp(Ω/ξ) -1 K-1 . (47)
Consider then the series summing all 1 N ∆S(K) with K ≥ 2. The series is convergent if ξ < ξ c = Ω log 2 , and divergent when ξ > ξ c . In the latter case, for a finite-N system, the maximum of |∆S(K)| is exponentially large in N , and is reached in K ≃ N 2 . As a consequence, for ξ > ξ c , the sum (27) can not be truncated according to the size of the clusters. This result is not specific to the dimension unity, and holds for other interaction networks. The expansion of S(p) defines an alternate series, and the order of its terms matters for its convergence in the N → ∞ limit. For an Ising model on a generic lattice with fixed degree (number of neighbours) v, the largest value of ξ such that the series (27) (after division by N ) is absolutely convergent in the N → ∞ limit is ξ c = Ω log v (Appendix D). A better truncation scheme consists in keeping cluster-entropies larger than a threshold Θ only. Let us define

S(p, Θ) = Γ⊂(1,2,...,N ) |∆SΓ(p)|>Θ ∆S Γ (p) . ( 48
)
The rationale is that, due to the properties of the cluster entropies and to the cancellation mechanism exposed in Section III C 2, summing large cluster-entropies may provide a good approximation to the true value of S(p). In the D = 1 Ising model case, the exact value of S(p) is, indeed, obtained as soon as Θ < ∆S (1,2) . We show in Fig. 6 the residual error in the cross-entropy due to the truncation as a function of the threshold Θ for the same small D = 2 grid as in Fig. 5. The error S(p, Θ) -S(p) is very small, and equal to 10 -6 when all clusters with contour length smaller than 4 are taken into account. As Θ is made smaller, clusters with larger contour lengths are summed up, and the error reaches the numerical accuracy ∼ 10 -14 . On top of this trend, positive fluctuations, corresponding to larger errors, arise when not all the clusters with the same interaction path (and length L) are summed up, and the cancellation of those contributions is not effective (Fig. 6 andcaption). We will study in more details this phenomenon in Section IV D. We now explain why the presence of noise in the data provides a compelling argument supporting the introduction of the cut-off Θ.

B. Distribution of small cluster-entropies in the presence of noisy data

In this Section, we investigate how limited sampling affects the values of the cluster-entropies. We assume that B configurations σ τ are sampled from the Gibbs distribution of an Ising model with interaction parameters J using a Monte Carlo procedure to generate the data p.

Universality at small |∆S|: numerical evidence

The empirical correlations, c ij = p ijp i p j , differ from the Gibbs correlations, σ i σ j Jσ i J σ j J , by random fluctuations of amplitude

c B ≃ p(1 -p) √ B , ( 49
)
where p is the typical value of the p i . For pairs i, j with weak Gibbs correlations (< c B in absolute value), the experimental correlations are dominated by the noise. As a consequence, the distribution of the cluster-entropies is universal for small ∆S. Its structure is a consequence of the noise in the data, and not of the interaction network of the model used to generate the data. Figure 7 shows the histograms H (full distributions) of the entropies ∆S (i,j,k) for the K = 3-clusters for a onedimensional Ising model, for three values of the numbers B of sampled configurations. The histograms are made of two components: a bell-shaped distribution at small |∆S|, and isolated peaks at larger |∆S|. The cluster-entropies corresponding to the isolated peaks have the same values as in the perfect sampling case (B = ∞, impulses). When B increases, the bell shapes move towards smaller entropies (in the log-scale of Fig. 7), and more peaks are unveiled in H.

We show also in Fig. 7 the histograms H IS for a system of Independent Spins (IS), with the same p i 's as the original system, and the same number B of sampled configurations. Contrary to H, H IS does not exhibit isolated Length FIG. 6: Effect of truncation on the bidimensional 3×3 grid with J = 1.777 (units of kBT ) and fields hi =j( =i) Jij [START_REF]0,1, the couplings (Jij ) and fields (hi) are[END_REF], for these parameters values the model has a phase transition between the paramagnetic and ferromagnetic phase and therefore the correlation length is proportional to the linear size of the system. Top: difference between the truncated and the true cross-entropies as a function of the cut-off on the absolute cluster-entropies, Θ. Bottom: contour lengths L(Γ) vs. Θ. The fluctuations of S(Θ) -S reflect the cancellation phenomenon. Summation of the 12 clusters of nearest-neighbours with K = 2 and L = 2 gives S(Θ = 0.1) -S ≃ 0.01, of the 21 clusters contributions corresponding to squared paths, e.g. 1-2-4-5), with K = 2, 3, 4 and L = 4 gives S(Θ = 0.002) -S ≃ 10 -6 . Fluctuations arise if only a part of the clusters that share the same interaction path are summed up, and cancellation is incomplete. For instance, fixing Θ = 0.0025 discards (1, 2, 4, 5), which has the same interaction path as (2,4).

peaks at well-defined, B-independent cluster-entropies. The histograms H IS concentrate around smaller |∆S| as the number B of configurations increases. Note that the histograms H IS roughly correspond to the bell-shape parts of the distributions H for the same value of B. We have checked that these features are largely independent of the particular sample and of the cluster size, K.

The histograms H IS depend on B through their standard deviation, σ IS (B). The calculation of σ IS (B) from the dominant contribution (C7) in the diagrammatic expansion of the cluster entropies (Section III C 1) is presented in Appendix E. We obtain that, for clusters of size K and in the case of uniform averages p i = p different from 0, 1 2 , and 1 [55],

σ IS (B) ≃ 3 K K! 8 (2p -1) 2 p(1 -p) 1 B K-1 2 . ( 50
)
Figure 8 shows how the small-entropy regions of the histograms H obtained for different B collapse onto each other once rescaled by σ IS (B). As expected, the rescaled H coincide with H IS in the small |∆S| ≤ σ IS (B) region, which concentrates most of the distribution (Fig. 8). The universality of the distribution at small ∆S is not specific to the one-dimensional Ising model, but holds, in the thermodynamic limit, for all interacting spin systems when the measured connected correlations are corrupted by noise. For a finite system in dimension D with correlation length ξ, we expect that the small-∆S is universal when N > ℓ D , where ℓ = ξ log(1/c B ). Indeed, the number of large c ij coming out of the noisy background is ≈ N ℓ D , while, for most of the N 2 pairs of spins i, j, the connected correlations have random values of amplitude c B . The full distribution H IS can be characterized analytically in the N → ∞ limit. Details can be found in Appendix E. We find the following scalings, depending on the value of the cluster size, K:

H IS (∆S) ∼ 1 |∆S| 2/3 for K = 2 , (-log ∆S) K-2 |∆S| for K ≥ 3 (|∆S| → 0) , ∼ exp -C 2 (B) |∆S| 2/(2K-1) (1 + o(1)) for every K ≥ 2 (large |∆S|) , (51)
where

C 2 (B) = 2×3 (K-1)/(2K-1) K (2K-3)/(2K-1) /(2K -1) 2 (σ ∆S ) -2/(2K-1
) is proportional to B, see Appendix E and equation (E21). The distribution is therefore characterized by a divergence at the origin, and stretched exponential tails. The scalings above were derived with the choice S 0 = S MF ; in the absence of the reference entropy, the stretched exponential has exponent 2 K instead of 2 2K-1 .

Finite-N effects and lower bound to the threshold Θ

The discussion about the localized peaks and the bell-shape distribution in H IS in the previous Section is an oversimplification. In reality, for finite systems, large fluctuations of the sampled correlations take place, and no clearcut boundary exist between cluster-entropies due to the noise and the ones deriving from the interaction network. From extreme value theory [START_REF] Gumbel | Statistics of Extremes[END_REF], the largest value of the correlations are of the order of c MAX ij = c B √ 4 log N . Therefore, the largest cluster-entropy is, according to (C7), of the order of

∆S max ≈ (4 log N ) (2K-1)/2 σ IS . (52)
A more detailed calculation to estimate where this fuzzy boundary between the signal and the noise in the entropy distribution takes place is presented below. Let M K (Θ) be the average number of clusters of size K with entropies |∆S| > Θ. According to (51),

M (Θ) = N K Θ d(∆S) H IS (∆S) ≃ exp -C 2 (B) Θ 2/(2K-1) (2K -1) N K Θ (2K-3)/(2K-1) 2 K! C 2 (B) , (53)
for large Θ and N (compared to K). The value of the threshold Θ such that M (Θ) = N α , with α < K, is, to the leading order in N , In particular, using the formula above for α = 0, it is likely that no cluster have entropy larger than Θ(0), in agreement with [START_REF]The minimal strength of the couplings which can be 'detected' depend on the quality of the sampling, and scales as log N/B[END_REF].

Θ(α) ≃ K -α C 2 (B) log N K-1 2 . ( 54
We have tested formula (54) through a computation based of a system of N = 40 Independent Spins, with uniform mean p = .0248; these parameters were chosen to mimick real data described in [2]. Figure 9(a) shows the number of clusters with entropies larger than Θ in absolute value, for a random set of B = 10 6 configurations (K = 3). The theoretical predictions based on (54) are in very good agreement with the simulations. The vast majority of clusters have entropies smaller than, say, Θ(2). On a smaller entropy scale, the histogram H IS of the small cluster entropies is strongly concentratred around zero as predicted in 51 (Fig. 9(b)).

As a conclusion, due to the sampling noise, most small cluster-entropies are random quantities, and provide no information about the underlying interactions parameters. Imposing a threshold Θ allows one to remove these artifact contributions. A lower bound to the value of Θ is given by (54), with, say α = 1 or 2. In practice, we will see that higher values of Θ may be sufficient for an accurate solution of the inverse Ising problem.

C. Properties of the susceptibility matrix and of its inverse

We now present a theoretical argument suggesting that the truncation scheme we have introduced is robust against an increase of the correlation length of the system. More precisely, the maximal size of the clusters to be summed up to reach an accurate solution of the inverse problem is not directly related to the correlation length, but rather depends on the structure of the interaction graph.

The susceptibility matrix χ (10) characterizes how the observables of the Ising model, such as the averages and correlations in p, are modified in response to an infinitesimal change in one or more interaction parameters in J. As far as the inverse Ising problem is concerned, it is more natural to ask the following question. Assume the inverse problem has been solved for a set of data p and the corresponding interations J have been found. Now imagine that the data are slightly changed, p → p + δp. How large will be the resulting change δJ in the interactions? The response function characterizing the inverse problem,

δJ δp = - ∂ 2 S(p) ∂p∂p = χ -1 , (55)
is simply the inverse of the susceptibility matrix χ. Whether the inverse problem is well-behaved or not will therefore depend on the properties of χ -1 . In particular, it will depend on the largest eigenvalues of χ -1 and on the structure of the corresponding eigenvectors.

A quantity which is closely related to (55) in liquid theory is the Ornstein-Zernike direct correlation function. The direct correlation is widely believed to be short-ranged, as the interaction potential [46]. This property is used in closure schemes such as the Percus-Yevick scheme to obtain the equation of state [START_REF] Percus | [END_REF]. We discuss below in details the property of the inverse susceptibility matrix in the case of the spherical model and of the unidimensional Ising model.

Case of perfect sampling: properties of χ and χ -1

Consider first the O(m) model, where each site i = 1, 2, . . . , N carries a m-dimensional real-valued spin vector

σ i = (σ 1 i , σ 2 i , . . . , σ m i ), of norm √ m.
As usual, two spins, say, i and j, are coupled through the dot product of their spin vectors, -J ij σ i • σ i (units of k B T ). Hence the interaction J ij couples the same component (α) of the spins in the pair i, j. The fields h α i , with α = 1, 2, . . . , m, are chosen to vanish for simplicity. In the large-m limit the model can be exactly solved [40]. The cross-entropy is equal to

S(p) = m 2 log det p + O log m , ( 56
)
where p is the N × N matrix with diagonal elements pii = 1 and off-diagonal elements pij , equal to the average of the product of the components α of spins i and j. The elements of the inverse susceptibility matrix are obtained by differentiating S(p) twice with respect to p,

χ -1 kl,k ′ l ′ = 1 2 J k,k ′ J l,l ′ + J k,l ′ J l,k ′ . (57)
Hence, the inverse susceptibility has the same structure as the interaction graph. In particular, if the coupling matrix J is sparse (has many zero elements), so is χ -1 . On the contrary, the susceptibility matrix χ is generally not sparse. The observation above is not specific to spherical spins. Consider now the D-dimensional Ising model with σ i = 0, 1 spins on a hypercubic lattice, with nearest neighbour interactions J ij . In the D = 1 case the susceptibility matrix (top right corner in matrix [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF]) is non zero for all i, i ′ : χ i,i ∝ x |i-i ′ | , where the proportionality constant does not depend on i, i ′ , and x = exp(-1/ξ). The inverse susceptibility matrix is a tridiagonal matrix [43]: the only non-zero elements are

(χ -1 ) ii = 1 + x 2 1 -x 2 and (χ -1 ) i,i±1 = - x 1 -x 2 . ( 58
)
As in the spherical model, the structure of the inverse susceptibility matrix is the same as the one of the interaction matrix. In dimension D ≥ 2 the inverse susceptibility matrix is not, strictly speaking, sparse. However it exhibits a much faster decay with the distance r = |ii ′ | than the susceptibility itself [56]. At the critical point, the latter decays as χ(r) ∼ r -(D-2+η) , where the critical exponent attached to the decay of the spin-spin correlation, η, vanishes in dimension D ≥ 4, and is positive and small in dimension D ≤ 3, i.e. η = 1 4 for D = 2. The inverse susceptibility is the Laplacian in dimension D ≥ 4, a purely local operator, and decays as χ -1 (r) ∼ r -(D+2-η) for D ≤ 3. While both quantities decrease as power laws in r, the inverse susceptibility has a much sharper decay than the susceptibility itself. In particular, the integrated contribution to the susceptibility coming from distances larger than R,

∞ R d D r χ(r) = R 1-η , (59)
diverges when R → ∞, while the same quantity calculated for the inverse susceptibility,

∞ R d D r χ -1 (r) = 1 R 3-η , (60)
tends to zero as R → ∞. This fact is a good news for the inverse problem. According to (55) the error on the field h i done when discarding all the spins at distance R > ǫ -1/(3-η) is of the order of ǫ only. In this regard, the inverse Ising problem remains local even at the critical point. While the discussion above is related to the response of a field h i to a change in the average p i ′ of spin i ′ , the response of a coupling J kl following a modification of the 2-point average p k ′ l ′ , see [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF], is also of interest. Unfortunately, to our best knowledge, this quantity has not been studied in the case of the Ising model so far. As a first step, we focus here on the D = 1-Ising model with uniform nearest-neighbour interactions, and in the thermodynamical limit (N → ∞). The four-spin connected correlation function is, up to a p-dependent multiplicative constant, equal to

χ ij,kl = x i4-i3+i2-i1 -x j-i+k-l , (61)
where i 1 ≤ i 2 ≤ i 3 ≤ i 4 are the same indices as i, j, k, l but sorted in increasing order, and x = exp(-1/ξ) < 1. We show in Appendix G that the inverse susceptibility matrix is given by

χ -1 ij,kl =                (1+x 2 ) 2 (1-x 2 ) 2 if i = k, j = l and j ≥ i + 2 , 1+x 2 +x 4 (1-x 2 ) 2 if i = k, j = l and j = i + 1 , -x(1+x 2 ) (1-x 2 ) 2 if i = k ± 1, j = l or i = k, j = l ± 1 , x 2 (1-x 2 ) 2 if i = k ± 1, j = l ± 1 , 0 otherwise . (62)
Hence, the inverse susceptibility matrix is sparse, with at most 9 non-zero elements per line, while the dimension of the matrix is 1 2 N (N -1) → ∞. In dimension D ≥ 2, we do not expect χ -1 to be sparse. However we conjecture that χ -1 ij,kl decays quickly with the minimal distance between the four points i, j, k, l (each index, e.g. j, is now a D-dimensional vector).

Influence of the sampling noise on the norms of χ and χ -1

To corroborate this statement we have carried out exact numerical analysis of small bidimensional grids (Section III C 2). We show in Fig. 10(a) the fraction of elements χ ij,kl of the susceptibility matrix larger than ǫ = 10 -7 in absolute value (the largest elements have magnitude ∼ 1). This fraction is closed to 1 for all the values J of the coupling we have studied. As expected, the inverse susceptibility matrix has many more small elements (Fig. 10(b)). In addition, the fraction of entries in χ -1 smaller than ǫ seem to increase with the size N of the grid.

In the presence of noise in the sampling process the inverse matrix χ -1 loses its quasi-sparse structure. More precisely, for the number B of sampled configurations chosen in Fig. 10(b), all the elements (χ -1 ) ij,kl are larger than ǫ in absolute value. Indeed, the quasi-sparsity χ -1 in the perfect sampling case reflects the sparse structure of the underlying interaction matrix. When data are corrupted by noise, the Ising model (over)fitting the data has no reason to be sparse anymore, and neither has the inverse susceptibility.

The influence of the sampling noise on the susceptibility matrix and on its inverse can be measured through the largest and smallest eigenvalue of χ, denoted by, respectively, λ max and λ min . According to Figs. 11(a,b), we have that:

• λ max increases with the size of the system (we expect λ max to diverge at the critical coupling J ≃ 1.778 in the thermodynamical limit), but is not affected by the sampling noise (the black and red/gray curves associated to the same size are nearly indistinguishable in Fig. 11(a)).

• λ min is not strongly affected by the system size in the case of perfect sampling. In case of noisy sampling, λ min acquires a smaller value. The effect of the noise increases with the system size (Fig. 11(b)). Those facts are observed from the study of the norms of the two matrices χ and χ -1 . Here, we define the norm of the matrix A through

||A|| = max i j |A i,j | . ( 63
)
Figures 11(c,d) show that the behaviours of the norms ||χ|| and ||χ -1 || are, from a qualitative point of view, similar to the ones of, respectively, λ max and 1/λ min . However the norms are directly related to the magnitudes of the elements of the matrices, according to (63). The independence of ||χ -1 || from the size N , contrary to the strong increase of ||χ||, supports the notion that most elements of χ -1 are very small (or even zero) in the case of perfect sampling. This property is lost when the sampling is not perfect: the presence of noise in the correlation makes the norm ||χ -1 || increases with N (Fig. 11(d)). Hereafter, we study how the error on the entropy S(Θ) resulting from the truncation varies with the threshold Θ and we discuss the fluctuations of S(Θ) -S observed in Fig. 6. We start by sorting the absolute values of the cluster-entropies |∆S Γ | in decreasing order:

∆S 1 ≥ ∆S 2 ≥ ∆S 3 ≥ . . . ≥ ∆S n ≥ . . . ≥ 0 . ( 64
)
We call η n = ±1 the sign of the cluster-entropy ∆S Γ attached (equal in absolute value) to ∆S n . Given the threshold Θ, we define n * (Θ) as the index of the smallest cluster-entropy larger than Θ: ∆S n * (Θ) ≥ Θ > ∆S n * (Θ)+1 . The truncated entropy (48) can be rewritten as S(p, Θ) = Σ(Θ), where

Σ(Θ) = n * (Θ) n=1 η n ∆S n . (65)
We want to study how Σ(Θ) behaves when Θ is made small. In particular, how does the difference ǫ s (Θ) = Σ(Θ)-Σ(0) behave as a function of Θ? Is it a smooth function, or does it exhibit large and irregular fluctuations? From a mathematical point of view, it is convenient to imagine that N → ∞. The above question can be formalized as whether 1 N Σ(Θ) converges to some limit value; the normalization factor comes from the fact that we expect the cross-entropy to be extensive in the system size N . Depending on the system under consideration, different situations can be encountered.

The most favorable case is when lim

N →∞ 1 N n≥1 ∆S n < ∞ . (66)
If this condition holds, the difference ǫ s (Θ) can be made arbitrarily small if Θ is small enough. An illustration is provided by the one-dimensional Ising model with small correlation length ξ and perfect sampling (B = ∞). For this model, the sequence of ∆S n is highly degenerate, and its distinct values are in one-to-one correspondence with the integer distances d ≥ 1 between the extremities of the clusters (Fig. 7). The cluster-entropy ∆S(d) asymptotically decays as exp(-3d/ξ), and has multiplity 2 d-1 , since each point between the extremities may or may not belong to the cluster. We find

1 N n≥1 ∆S n ≃ d≥1 2 d-1 exp(-3d/ξ) , (67)
which converges if ξ < ξ c = 3 log 2 . The calculation above is very similar to the one of Section IV A. Indeed, when the series with general term ∆S(K) is absolutely convergent, any ordering of the cluster-entropies is possible. In particular, one is allowed to sum all the clusters of a given size K as proposed at the beginning of Section IV A.

What happens when condition (66) is violated? Again consider the one-dimensional Ising model. For perfect sampling, the cancellation property discussed in Section III C 2 ensures that 1 N Σ(Θ) has reached its limit ∆S(d = 1) as soon as Θ < ∆S(1). In the case of noisy sampling (finite B), the situation is more complex. In the presence of noise in the correlations c kl the cluster-entropies with the same distance d between extremities are not degenerate any longer. We show in Fig. 12(a) the value of Σ as function of Θ for a large correlation length ξ compared to ξ c , and B = 10 5 sampled configurations. We observe the appearance of 'packets' of cluster-entropies, located around the noiseless values ∆S(d ≥ 2). The width of a packet depends on the amount of noise due to the sampling, i.e. on the number B of sampled configurations. The values of Σ at the two edges of the packet are very close to one another due to the cancellation property. As Θ spans the range of cluster-entropies in the packet, Σ fluctuates. The maximal amplitude of the fluctuations seems to weakly increase as we look at packets with smaller and smaller entropies (Fig. 12(a)).

We have analyzed the statistics of the signs ǫ n of the cluters -entropies in (65). Writing the sequence of signs η = (η 1 , η 2 , η 3 , . . .), we consider the blocks j of contiguous and equal signs, and defines their lengths ℓ j . For instance, the block lengths corresponding to η = (+, +, -, +, +, +, +, -, -, -, -, -, +, -, ....) are ℓ 1 = 2, ℓ 2 = 1, ℓ 3 = 4, ℓ 4 = 5, ℓ 5 = 1, .... The histogram of the block-lengths is shown in Fig. 13. The two main features are:

• The frequence of ℓ decreases exponentially when ℓ ≪ N , and is in very good agreement with the exponential law 1 2 ℓ .

• A large 'structural' block of length ℓ ≃ N is present. This block corresponds to the N clusters of size K = 2 (having all sign +), and the cluster of size K = 3 with largest entropy (which has the same sign +).

We have calculated the correlation between successive block lengths, normalized by the variance of the block length,

ρ = l j l j+1 -l j 2 l 2 j -l j 2 , (68)
where (•) denotes the average over the blocks j. For the model and the data shown in Figs. 12 and 13, we find ρ ≃ .012. Changing the set of sampled configurations does not affect the amplitude of the ratio ρ, which is always found to be about 1%. This ratio coincides with the inverse of the square root of the number of blocks, equal to a few thousands. Hence, the analysis is compatible with the absence of any correlation between the lengths of successive blocks. The same conclusion is reached with experimental data, e.g. multi-electrode recordings of the activity of a neural population [2,37] (not shown).

The simple statistics sets above suggests the following 'random sign' model, allowing us to deepen our theoretical understanding of the behavior of Σ(Θ). In the random sign model, the signs η n are replaced with random variables, equal to ± with probabilities 1 2 , and independent from each other. We emphasize that, in Σ(Θ) defined in (65), the signs are deterministic (for given data p). The random sign model is therefore an approximation motivated by the statistical analysis above. Assume now that the value chosen for the threshold Θ falls within a packet p including N p clusters. Fluctuations of the order of ∆Σ ∼ ±Θ N p (69) are expected on the entropy. As Θ decreases, the size of the packets, N p , tends to be bigger. Loosely speaking, smaller entropies correspond to longer interaction paths, shared by many more clusters. In the case of the one-dimensional Ising model, as Θ decreases, the distance between the extremities of the clusters involved in a packet, d, increases. We have N p = 2 d-1 ; hence, ∆Σ ∼ exp(-3d/ξ) √ 2 d . We conclude that the error on the entropy tends to zero if ξ < 2 ξ c . From the above discussion, it appears that a general, sufficient condition for the amplitude of the fluctuations to vanish as Θ → 0 is lim

N →∞ 1 N n (∆S n ) 2 < ∞ . (70)
Indeed, if condition (70) is fulfilled, the sum of the fluctuations due to all packets corresponding to cluster-entropies smaller than Θ is guaranteed to vanish with Θ. Hence condition (70) not only ensures that the fluctuations ∆Σ attached to the packet 'cut' by Θ vanishes, but also that the error on the entropy, ǫ s (Θ), tends to zero when Θ → 0.

It is important to realize that the guarantee is of probabilistic nature. Arbitrary large fluctuations are possible (in the N → ∞ limit), but are very unlikely. More precisely, within the random sign model, the error is a normal variable,

ǫ s (Θ) = N 0, 1 N n>n * (Θ) (∆S n ) 2 , (71)
with a variance vanishing with Θ according to (70). The true error is expected to be even smaller than the random sign estimate (71). Indeed, packets need not be isolated from each other as in Fig. 12. In the presence of a strong sampling noise, or in higher dimension than D = 1, packets will overlap. As a consequence, the number of packets 'cut' by the threshold and their size will determine the amplitude of ∆Σ. Further investigations of those points are needed.

V. ADAPTIVE ALGORITHM FOR THE INVERSE ISING PROBLEM

A. Procedure to construct and select clusters

As explained above discarding the cluster-entropies smaller than a threshold Θ is an efficient step against overfitting of the sampling noise. In addition, for systems with dilute and strong interactions, we expect that only the clusters of neighboring sites on the interaction network will have substantial entropies. These arguments provide a heuristic basis for the threshold-based truncation of the expansion (27).

How can we implement the truncation scheme in practice? The combinatorial explosion of the number of clusters of size K among N sites impedes any brute force computation approach, as soon as N is larger than a few tens. Even for small-N system for which it is feasible, computing ∼ 2 N cluster-entropies and, then, discarding most of them does not sound like an efficient procedure.

We propose below an alternative approach, based on a recursive and selective construction of relevant clusters. The approach is based on the principle that clusters with large entropies should be compatible with the interaction network to be inferred. Suppose that two clusters Γ and Γ ′ have both large entropies, and share most of their spins. Then, the union Γ ∪ Γ ′ is a good candidate for a bigger cluster. If the entropy of the union cluster is large, a new part of the interaction network will be unveiled. Conversely, if it is small, no new interaction path with respect to the one discovered from Γ and Γ ′ separately exists. Hence, combining strongly overlapping clusters should allow us to progressively deepen our knowledge of the local structure of the interaction graph.

The above heuristics is formalized as follows:

A1. Initial step: build the list of all clusters of size one: L 1 = {(i) : i = 1, 2, . . . , N }. All the other lists L K for K ≥ 2 are empty.

A2. Iteration: assume the current size of clusters is K ≥ 1, i.e. L K is not empty while L K+1 is empty. For every pair Γ 1 , Γ 2 in L K :

A21. Construction: build Γ = Γ 1 ∪ Γ 2 A22. Selection: if Γ is of size K + 1 and if |∆S Γ (p)| ≥ Θ, then select Γ and add it to L K+1 .
A3. Recursion: if at least one cluster has been selected, then add 1 to K, and go to step 2 to pursue the construction process. If no cluster has been selected, the construction process is over.

The first condition in A22 is about the size of Γ. The union of two clusters of size K has size K + 1 if and only if they have exactly K -1 common spins. Γ 1 = (i 1 , i 2 . . . i K-1 , x) and Γ 2 = (i 1 , i 2 , . . . , i K-1 , y) can be merged into Γ = (i 1 , i 2 , . . . , i K-1 , x, y); the ordering of x, y, and of the i l 's is irrelevant here.

B. Calculation of ∆SΓ(p)

Step A22 requires the calculation of the cluster-entropy ∆S Γ (p) for each selected cluster Γ (of size K). In order to do so we make use of the formula

∆S Γ (p) = S Γ (p) -(S 0 ) Γ (p) - Γ ′ ⊂Γ (Γ ′ =Γ) ∆S Γ ′ (p) , (72)
which can be easily deduced from (27). The procedure is as follows:

B1. calculate the subset-entropy S Γ (p) through the minimization of S Ising (J|p) (7) with respect to the fields and couplings. The partition function Z[J] is computed as the sum over the 2 K configurations of the spins in Γ.

B2. substract the reference entropy (S 0 ) Γ (p). For the mean-field reference entropy, (S 0 ) Γ (p) = 1 2 log detM Γ (p), according to formula (21); M Γ (p) is the K × K restriction of matrix M (p) to the indices i 1 , i 2 , . . . , i K in Γ. In presence of a regularization term ( 13) equation ( 36) has to be used instead of (21) to calculate (S 0 ) Γ (p).

B3. Substract the entropies ∆S

Γ ′ (p) of all the sub-clusters Γ ′ of size K ′ < K, included in Γ.
The last step (B3) assumes that the entropies of all the sub-clusters of Γ are known, i.e. have been computed at a previous step in the algorithm. This is true for K ′ = 2, but not necessarily so for K ′ ≥ 3. To circumvent this difficulty we maintain at all times during the execution of the algorithm the list L all of all the clusters and of their entropies calculated so far; L all is a larger list than the one of the selected clusters (union of all L K ). The procedure to compute ∆S Γ (p) is then: B0. build the list L of all the sub-clusters Γ ′ in Γ not already present in L all . For each Γ ′ ∈ L, starting from the smallest sub-cluster and ending up with the largest one, run steps B1, B2, B3 to obtain ∆S Γ ′ (p), and add Γ ′ and its entropy to the list L all .

The ordering of L ensures that all the sub-clusters of Γ ′ required to calculated its entropy are in L all when step B3 is executed.

C. Calculation of the cross-entropy, couplings and fields

Once the construction process is finished, the list L sel = L 1 ∪L 2 ∪L 3 ∪. . .∪L Kmax of all selected clusters is available. Here, K max is the size of the largest cluster selected by the construction procedure. We then C1. estimate the cross-entropy through

S(p) = S 0 (p) + Γ∈L sel ∆S Γ (p) . ( 73
)
Next we need to estimate the values of the fields and of the couplings, solution to the inverse Ising problem. One possibility would be to use recursion relations similar to (72) for ∆h i,Γ (p) and ∆J ij,Γ (p), that is, the contributions to, respectively, the field h i and the coupling J ij coming from the cluster Γ. Next we could sum up those contributions over the clusters included in L sel . However, to save memory space, it is possible to resort to the following, alternative procedure:

C2. define the 'multiplicities' m Γ of the subsets Γ through:

C21. let L sub be the list of all clusters in L sel and of all their subsets. Initialize m Γ = 0 for every Γ ∈ L sub .

C22. For each Γ ∈ L sel , and for each Γ ′ ⊂ Γ, add (-1) K-K ′ (see (29)) to m Γ ′ , where K, K ′ are the sizes of, respectively, Γ, Γ ′ . The sub-clusters Γ ′ = Γ must be taken into account in the addition process.

C3. estimate the fields and the couplings through

h i (p) = (h 0 ) i (p) + Γ∈L sub :(i)⊂Γ m Γ h i,Γ (p) -(h 0 ) i,Γ (p) , J ij (p) = (J 0 ) ij (p) + Γ∈L sub :(i,j)⊂Γ m Γ J ij,Γ (p) -(J 0 ) ij,Γ (p) . (74)
The fields h i,Γ and the couplings J ij,Γ in step C3 above are the ones obtained through the minimization of S Ising (J|p) over J = {h i,Γ , J ij,Γ } in step B1. The fields (h 0 ) i and the couplings (J 0 ) ij are (minus) the derivatives of the reference entropy S 0 (p) with respect to p i and p ij , see formulas (22). The fields (h 0 ) i,Γ and the couplings (J 0 ) ij,Γ are their counterparts for the subset Γ only, i.e. the derivatives of (S 0 ) Γ (p); their expressions are given by ( 22) again, upon substitution of the N × N matrix M (p) with the K × K matrix M Γ (p) restricted to the K elements of Γ only.

D. Pseudo-code of the algorithm

We now give the pseudo-code useful for the implementation of the procedures above. To improve the readability the code is broken into several parts.

We start with Algorithm 1, which computes the cross-entropy and the reference entropy for a subset Γ. The energy function H Ising is defined in [START_REF] Balakrishnan | [END_REF]. The minimization over J can be done using standard numerical algorithms for convex functions. A speed-up is generally obtained when we start with J MF , the value of the interaction parameters obtained from the MF approximation (22), as an initial guess for the value of J [45]. In the absence of regularization, the parameter γ is set to 0. It is straightforward to change the pseudo-code to introduce the L 1 -regularization instead of the L 2 -norm, see formulas [START_REF] Wallace | Statistical and Inductive Inference by Minimum Message Length[END_REF].

Algorithm 1 Computation of entropy S Γ (p) -(S 0 ) Γ (p)
Require: Γ (of size K), data p, regularization parameter γ Computation of SΓ:

Define SIsing Γ [J|p] ← log   σ∈{0,1} K exp -HIsing[σ|J]   - i∈Γ hi pi - i<j∈Γ Jij pij + γ i<j J 2 ij pi (1 -pi) pj (1 -pj) , where J = {hi, Jij } is of dimension 1 2 K(K + 1). SΓ(p) ← min J SIsing Γ [J|p] Computation of (S0)Γ: MΓ ← K × K matrix with elements (MΓ)ij = p ij -p i p j √ p i (1-p i )p j (1-p j ) with i, j ∈ Γ (S0)Γ(p) ← 1 2 log det MΓ if γ = 0, or use formula (36) if γ > 0. Output: SΓ(p) -(S0)Γ(p)
Algorithm 2 calculates the entropy ∆S Γ of the cluster Γ and maintains the list L all of all cluster-entropies computed so far. It calls Algorithm 1 as a subroutine.

Algorithm 2 Computation of cluster-entropy ∆S Γ (p)

Require: Γ (of size K), data p, list L all = {Γ ′ , ∆S Γ ′ (p)} of known cluster-entropies.

L ← {Γ ′ : Γ ′ ⊂ Γ and Γ ′ / ∈ L all } (ordered in increasing sizes) for Γ ′ ∈ L do ∆S Γ ′ (p) ← SΓ(p) -(S0)Γ(p) - Γ ′′ ⊂Γ ′ (Γ ′′ =Γ ′ )
∆S Γ ′′ (p) using list L all of cluster-entropies calculated so far update L all ← L all U {Γ ′ , ∆S Γ ′ (p)} end for

Output: ∆SΓ(p) and L all

We can now give the core part of the procedure, which produces the list of selected clusters: Algorithm 4 calculates the estimates for the total cross-entropy, and for the interaction parameters once the list of selected clusters L sel has been obtained. It requires Algorithms 1 and 2; function (S Ising ) Γ and matrix M Γ are defined in the pseudo-code of Algorithm 1. compute S0(p) using formula ( 21) or ( 13)

Algorithm 3 Adaptive cluster algorithm for the inverse Ising problem

Require: data p, threshold Θ L1 ← {(i) : i = 1, 2, . . . , N } L sel ← ∅ K ← 1 while LK is not empty L sel ← L sel U LK K ← K + 1 LK ← ∅ for Γ1, Γ2 ∈ LK-1 do Γ ← Γ1 ∪ Γ2 if Γ is of size K -
SΓ(p) ← S0(p) + Γ∈L sel ∆SΓ(p)
Computation of fields and couplings J = {hi, Jij }: compute {(h0)i(p), (J0)ij(p)} using formula ( 22)

L sub ← {Γ ′ ⊂ Γ : Γ ∈ L sel } for Γ ′ ∈ L sub do m Γ ′ = Γ∈L sel :Γ ′ ⊂Γ (-1) |Γ|-|Γ ′ | ,

VI. APPLICATIONS

In this Section we report the results of our algorithm when applied to data generated from Ising models with diverse interaction structures and various numbers B of sampled configurations. We define:

• the number N clu of clusters generated by the algorithm and the size K max of the largest clusters.

• the average error on the inferred couplings and fields:

ǫ h = 1 N i (h inf i -h i ) 2 1 2 , ǫ J =   2 N (N -1) i<j (J inf ij -J ij ) 2   1 2 . ( 75
)
Here, J inf ij and h inf i denote the values of, respectively, the inferred couplings and fields, while J ij and h i are the values of the couplings and fields in the model used to generate the data.

• The error bars δh i and δJ kl on the inferred couplings and fields, resulting from the finite sampling. Those statistical fluctuations are asymptotically given by the inverse of the susceptibility matrix of the cross-entropy S Ising , see equation [START_REF] Jerrum | The Markov chain Monte Carlo method: an approach to approximate counting and integration in Approximation Algorithms for NP-hard Problems[END_REF]. The entries of χ can be calculated from a Monte Carlo simulation, to estimate the multi-spins correlations. In practice, a good approximation of χ can already be obtained from the empirical average over the B configurations in the sampling set. This procedure avoids the use of the Monte Carlo. In the presence of a L 2 -regularization ( 13), γ p k (1p k )p l (1p l ) is added to the diagonal element χ kl,kl of the susceptibility matrix, before the inversion is performed. Hence, all the eigenvalues are strictly positive and the inverse is well defined. The inversion of χ can be done with standard linear algebra routines.

Inferred couplings are called 'reliable' when their absolute value is larger than three times their statistical error-bar: |J kl | > 3 δJ kl .

• The reconstructed observables, p rec i and c rec ij , which we compare to the data, p i and c ij . Those reconstructed averages are obtained using Monte Carlo simulations of the Ising model with the inferred fields, h inf i , and couplings, J inf ij . For those simulations the number of sampled configurations is chosen to be much larger than B, e.g. 100 B, to minimize the uncertainty on the reconstructed averages.

• the relative errors on the reconstructed averages p i and connected correlations c ij , with respect to their statistical fluctuations due to finite sampling:

ǫ p = 1 N i (p rec i -p i ) 2 (δp i ) 2 1 2 , ǫ c = 2 N (N -1) k<l (c rec kl -c kl ) 2 (δc kl ) 2 1 2 . ( 76
)
where the denominators in (76) measure the typical fluctuations of the data expected at thermal equilibrium, see [START_REF] Cagliotti | [END_REF], and

δc kl = δp kl + p k δp l + p l δp k . ( 77
)
If not explicitly stated otherwise we start from the value Θ = 1 for the threshold and run the algorithm several times, dividing the threshold by 1.01 after each execution. The algorithm is stopped when both errors ǫ p and ǫ c are close to 1. We call Θ * the final value of the threshold corresponding to this criterion. Unless explicitly stated otherwise a L 2 -regularization term ( 13) is present, with γ = 1/(10 B p 2 (1-p) 2 ), where p is the average value of the p i 's (Appendix A). As explained in Section II B the regularization term is important in case of undersampling and guarantees the convergence of the numerical minimization of S Ising .

A. Independent Spin Model

It is instructive to run first the algorithm on the Independent Spins model, where each spin has a probability p i to be 1, 1p i to be 0, independently of the other variables. Due to the noise in the sampling (finite value of B), the connected correlations c ij are not equal to zero. Figure 14 shows the outcome for a system of N = 40, as a function of the threshold Θ. The errors of reconstruction, ǫ p and ǫ c , are already smaller than one for the initial threshold value Θ * = 1. For this value of the threshold, cluster of size one only are selected. In other words, the interaction network J 0 , calculated from the reference entropy S 0 = S MF alone, is already overfitting the data as it attempts to reproduce the correlations due to statistical fluctuations. For smaller thresholds Θ contributions from clusters of size K ≥ 2 allow for an even more precise reproduction of the data .

The histogram of the inferred couplings, J inf ij , is shown in Fig. 15. It is centered around zero, and is approximately Gaussian. The standard deviation of the distribution is compatible with the statistical error bar on couplings δJ ij [START_REF] Jerrum | The Markov chain Monte Carlo method: an approach to approximate counting and integration in Approximation Algorithms for NP-hard Problems[END_REF] averaged on all the couplings. For the particular case of Fig. 14, 815 of the 820 inferred couplings are away from zero by less than three error bars, and are, therefore, classified as unreliable. This result is compatible with the fact that the non-zero couplings are the consequence of overfitting and do not reflect any real interactions.

Another possibility to avoid overfitting in this case is to apply the cluster expansion to the entropy S in the absence of reference entropy (S 0 = 0). We find that, for Θ * = 1, the reconstruction errors are ǫ p = 0.07 and ǫ c = 0.8. Therefore only one-spin clusters are taken into account, and all couplings are equal to zero exactly (since J 0 = 0).

B. Unidimensional Ising model

We now test the algorithm on the unidimensional Ising model with first-neighbor interactions, J i,i+1 = J, and uniform fields, h i = h. The model is placed on a ring with N sites (periodic boundary conditions). Data p can be computed exactly (Appendix F) or through an average over B configurations, sampled by Monte Carlo simulations. We compare the performance of the inference procedure for various values of B and two values of J, h, corresponding 1. Accuracy of the cluster expansion as a function of the threshold: errors on the entropy, couplings and fields

We start with the small-ξ case. Figure 16(top) shows ǫ S , the absolute value of the difference between the crossentropy S(p, Θ) (48) and the entropy of the model for perfect sampling, and the errors ǫ J and ǫ h , for various values of B. We observe that ǫ S sharply decreases around Θ 1 = 0.05, that is, the entropy of nearest-neighbor clusters ∆S (i,i+1) ; discarding all entropies smaller than this value would be exact in the perfect sampling case (Section III C 2). As Θ is decreased, ǫ S exhibits fluctuations centered around a discrete sequence of threshold values, Θ d , with d ≥ 2. As explained in Section III C 2, these values correspond to the cluster-entropies ∆S (i,i+d) in the absence of noise, see identity (45). Fluctuations spread over a small window around Θ d . They correspond to imperfect cancellations of 'packets' of entropies whose associated clusters share the same interaction path with length L = 2d (Section IV D). Since the correlation length ξ is small, and, therefore, the cross entropy expansion is absolutely convergent, the magnitude of the fluctuations quickly decreases with d. For B = 10 5 two bursts of fluctuations are visible (corresponding to d = 2, 3). For B = 10 7 , fluctuations are smaller and spread over narrower windows; only the d = 2 burst can be observed. In between two bursts of fluctuations, ǫ S reaches a plateau. Note that the value of ǫ S on the plateau is not zero due to the sampling noise (finite B).

The errors on the inferred couplings and fields have the same behaviour as ǫ S (Fig. 16(top)). Their magnitude are comparable to the expected statistical fluctuations calculated from the 4-spin correlations (inverse susceptibility in ( 17)), which decrease as B -1/2 .

We now test our algorithm on the unidimensional Ising model with a large correlation length, ξ = 9. The errors ǫ S , ǫ J , ǫ h are shown in Fig. 16(bottom) as functions of Θ. The correlation length of the model is larger than 2ξ c , at which the series of the squared cluster-entropies is divergent. Therefore, we expect large fluctuations of ǫ S , corresponding to packets of clusters with the same interaction path (Section IV D). Figure 16 shows, indeed, that fluctuations are much larger for ξ = 9 (bottom) than for ξ = 1 (top). Furthermore, fluctuations do not decrease much in amplitude as Θ is decreased. In the case of severe undersampling (B = 1000) the notion of bursts of fluctuations separated by plateaus is blurred out. The distributions of cluster-entropies in a packet is so wide that it overlaps with the distributions of entropies associated to the neighbouring packets. As for the previous case the magnitude of ǫ J , ǫ h are comparable to the expected statistical fluctuations calculated from [START_REF] Jerrum | The Markov chain Monte Carlo method: an approach to approximate counting and integration in Approximation Algorithms for NP-hard Problems[END_REF].

Quality of reconstruction and choice of the threshold Θ *

We now study the reconstruction errors ǫ c and ǫ p as functions of Θ, for the small and large correlation lengths and for the weak and strong sampling noises. We show in Fig. 17 the errors ǫ c and ǫ p , as well as the maximal size, K max , and the number, N clu , of clusters vs. the threshold Θ in the case of a small correlation length, ξ = 1 . The threshold at which both ǫ c and ǫ p are close to 1 can be chosen in the range 0.0007 < Θ * < 0.003, for which all the clusters of lengths L ≤ 4 are processed. It is possible to check in Fig. 16(top) that this threshold value gives the minimum values of ǫ S , ǫ J , ǫ h . Contrary to the case of perfect sampling, it is not sufficient to take into account clusters with contour length L = 2 only. The selected clusters correspond to three groups of N = 30 clusters each; the first group gathers the clusters (i, i + 1) (L = 2), the second one, the clusters (i, i + 2) (L = 4) and the third one, the clusters (i, i + 1, i + 2) (L = 4). In Fig. 17 we show the reconstructed averages p i and correlations c ij , at Θ * and for the largest sampling noise case B = 10 5 , vs. their values in the data . The agreement is very good, and falls within the statistical fluctuations expected at equilibrium for the Ising model, given in (77). We stress that the optimal value of the threshold and the maximal size of selected clusters depend on the particular realization of the data p, and can vary from sample to sample. By further decreasing the threshold Θ below Θ * , the reconstruction errors ǫ p and ǫ c decrease to values smaller than one (Fig. 17). This regime corresponds to an overffiting of the data, as the errors on the couplings, ǫ J , and on the fields, ǫ h , cease to decrease (Fig. 16(top)).

Results for the case of a larger correlation length (ξ = 9) with B = 10 5 and B = 10 3 sampling configurations are reported in Fig. 18. For the poor sampling case (B = 10 3 ) plateaus are not present any longer, but a good inference is still obtained for a value Θ * of the threshold, which, as in the ξ = 1 case, corresponds to the summation of all the clusters with contour length L = 4. This finding supports the discussion of Section IV C: the contour length required for a good inference is largely independent of the correlation length. However, in the poor sampling case, finding the right value for Θ * is harder for larger ξ due to the mixing of packets.

Quality of the inference: histograms of couplings

To better understand the quality of the inference we plot in Fig. 19 (up and middle panels) the histogram of the inferred couplings J inf ij at the threshold Θ * . The distribution is bimodal: a Gaussian-like peak centered in J inf = 0 and a smaller distribution around J inf = 4. The two sub-distributions are separated by a wide gap. The inference algorithm makes no classification error: the sub-distribution centered around J inf = 0 contains all the pairs (i, j) such that J ij = 0, and the one around J inf = 4 includes all the pairs of nearest neighbours (i, i + 1). All the couplings centered around zero are unreliable. Moreover the standard deviation of the distribution of the couplings around the zero value (equal to the minimal value of ǫ J reached on the plateau in Fig. 16 ) agrees with the statistical fluctuations [START_REF] Jerrum | The Markov chain Monte Carlo method: an approach to approximate counting and integration in Approximation Algorithms for NP-hard Problems[END_REF]; all the couplings around zero are therefore unreliable. The structure of the interaction network is perfectly recovered.

We show the histogram of couplings for Θ > Θ * , i.e. ǫ c > 1, in Fig. 19 (bottom); the structure of the interaction network is still perfectly recovered but the values of the positive inferred couplings is less accurate. The histogram of inferred couplings for the large correlation length , ξ = 9, is shown in Fig. 20. Even when the sampling noise is large, a good separation of the two sub-distributions corresponding to interacting and non-interacting pairs of spins is achieved for large values of the threshold Θ * , and the couplings are correctly inferred (up to the expected statistical fluctuations).

C. Regular bidimensional grid

We now analyze the performance of the algorithm on bidimensional grids of different sizes, N = N ′ × N ′ . Nearest neighbors on the grid interact through the coupling J. The value of J ≃ 1.778 is chosen to make the grid critical in the thermodynamical limit, N ′ → ∞ [47,48]. Hence, the system is at the paramagnetic/ferromagnetic critical point, and the correlation length ξ ′ diverges with N ′ .

We have described in Section IV A and Fig. 6 the partial cancellation of the cluster-entropies for a small bidimensional grid (N ′ = 3), and no sampling noise. Due to this cancellation property, taking into account clusters of contour length L ≤ 4 was sufficient to obtain a very accurate approximation to the cross-entropy. Hereafter, we show that this result is not affected by the presence of sampling noise. Furthermore, we will see that the size of the clusters necessary for a good inference of the interactions remains rather constant when the grid size is increased from N ′ = 7 to N ′ = 20, and is thus, as in Section VI B, largely independent of the correlation length ξ.

1. The small 3 × 3 grid revisited: influence of the sampling noise

We start with the 3 × 3 grid of Section IV A, for which the summations of all clusters up to size K = 9 gives the exact solution of the inverse problem, and all 1-and 2-point averages can be calculated exactly in the perfect sampling case. The reader is kindly referred to Fig. 6 and to the related discussion. Figure 21 shows the errors on the entropy, the couplings, the fields, the reconstructed 2-and 1-point averages, and the size and the number of clusters as functions of Θ. For a given amount of sampling noise (set by the value of B), the errors ǫ S , ǫ J , ǫ h follow their perfect-sampling counterparts, until a threshold value Θ sat , and they saturate for Θ < Θ sat . The saturations are interrupted by fluctuations due to the imperfect cancellations of clusters within a packet (Section IV D). The saturation values of ǫ J and ǫ h decrease with increasing B, and are compatible with the expected statistical fluctuations δJ ij given by [START_REF] Jerrum | The Markov chain Monte Carlo method: an approach to approximate counting and integration in Approximation Algorithms for NP-hard Problems[END_REF]. The value of Θ sat approximately coincides with the threshold Θ * , at which both ǫ c and ǫ p are ≃ 1 (Fig. 21). For B = 4500, only clusters of size K = 2 are taken into account at Θ * . For B = 10 7 , clusters made of spins on the elementary squares of Fig. 5 and of size up to K = 4 are selected, e.g. 1, 5 or 1, 2, 4, 5 with L = 4 in Fig. 5, while clusters such as (1, 2, 3) or (1, 3), which have the same contour length but are not on elementary squares, are discarded. At Θ * the histogram of the inferred couplings is made of two far apart sub-distributions (not shown): the first one corresponds to the 12 non zero couplings and is centered around J inf ≃ 1.8, and the second one, peaked around J inf = 0, contains the remaining 24 pairs of sites. Hence, the structure of the interaction graph is correctly found back.

Study of larger critical grids

We now turn to larger grids N ′ × N ', where N ′ ranges between 7 and 20. Data are calculated from B = 4500 configurations sampled through a Monte Carlo simulation. The error ǫ J on the couplings is shown in Fig. 22. As Θ decreases, ǫ J saturates to a value close to the average of the expected statistical error, δJ ij , which lies between .1 and .2. Saturation begins at large values of the threshold, even when the linear size N ′ of the grid is increased. The asymptotic values depends strongly on N due to the non periodic boundary conditions [57].

The reconstructions errors ǫ p and ǫ c are shown in Fig. 23. The maximal size of the clusters at the optimal threshold Θ * is bounded (K max ≤ 4), as the correlation length ξ diverges with N ′ . As a consequence, the running time of the algorithm increases linearly with N . For threshold values smaller than Θ * , ǫ p and ǫ c decrease. The 7 × 7 grid in Fig. 23 provides a clear illustration of data overfitting. Keeping B fixed while N ′ and N increase make the data effectively more and more noisy. This is the reason why the value of Θ * slightly increases with N .

-1 0 1 2 3 0 1 2 3 -1 0 1 2 3 0 1 2 3 -1 0 1 2 3 inferred couplings J
The histograms of the inferred couplings at the threshold Θ * are shown in Fig. 24. The structure of the grid is perfectly reconstructed for all sizes N ′ . We find that all the inferred couplings in the sub-distribution centered around zero are smaller (in absolute value) than three times their standard deviation (corresponding to the asymptotic value of Fig. 22) and are, therefore, unreliable.

D. Randomly diluted bidimensional grid

We now remove a fraction 1-ρ of the couplings on the grid, independently and at random. Our goal is to study how the algorithm performs on such disordered systems, at the phase transition and in the low temperature phase. We will compare the performance with another low complexity algorithm, the regularized logistic regression algorithm, guaranteed to perform well at high temperature and to fail at low temperature [START_REF] Ravikumar | [END_REF]. To compare with the numerical experiments of [16], we have generated 7 × 7 bidimensional grids, and keep each bond with probability ρ = .7. The remaining bonds are all equal to J [58]. We generate, for each value of J ranging from 0.4 to 4.4, eight randomly diluted grids. For each grid we calculate the data p by sampling over B = 4500 configurations generated by a Monte Carlo dynamics.

1. Inference of the network structure from the mean field entropy SMF Our first task consists, as in [16], in reconstructing the structure of the interaction graph only. We do not want to accurately determine the value of the coupling constants J ij , but only if it is positive or null. This task is easier than the precise inference of the couplings, and we will first handle by approximating the cross-entropy S with the reference entropy S 0 = S MF only. Equivalently, we choose Θ * to be large enough that no cluster is selected by our algorithm. We compute the mean-field couplings, (J 0 ) ij , and, for each pair (i, j), and decide that a bond is present if (J 0 ) ij > J/2, absent otherwise. The performances of this simple, Mean-Field algorithm are shown in Fig. 25. We say that the neighborhood of a vertex i is reconstructed if the sets of its neighbors j (J ij = 0) and of its non-neighbors (J ij = 0) are correctly inferred. In Fig. 25 (top panel) we report the fraction Q succ of the neighborhoods which are reconstructed (straight line), as a function of the coupling strength J. We compare the performance of the mean-field FIG. 25: Probabilities that a neighborhoods is reconstructed (top) and that a bond is inferred (bottom) as functions of the coupling J for a bidimensional random 7 × 7 grid density ρ = .7 and for various values of the numbers B of sampled configurations. Error bars are calculated from the standard deviations over eight random grids. Probabilities were obtained from the simple Mean-Field algorithm (S0 = SMF , Θ = 1). The value of Qsucc (top) is compared to the performances of the pseudo-likelihood algorithm (Rlr) of [START_REF] Ravikumar | [END_REF]16] in the top panel.

algorithm with the pseudo-likelihood algorithm of [START_REF] Ravikumar | [END_REF] in Fig. 25(right). Contrary to the pseudo-likelihood algorithm case, the neighborhoods are perfectly reconstructed at the phase transition. Furthermore, Q succ remains large in the ferromagnetic phase: for instance, for J = 3.2, more than 80% of neighborhoods are perfectly inferred. For very large J (low temperatures) the average p i are too close to 0 in the down state and to 1 in the up state. Most of the sampled configurations coincide with one of the two ground states, and the inference is difficult. Fig. 25(top panel) shows the increase of Q succ resulting from a ten-fold increase of B.

Another measure of the performances is shown in Fig. 25(bottom panel). We plot P succ , the fraction of bonds in the grid correctly predicted to exist, averaged over the data realizations, as a function of the coupling strength J. For J = 3.2 more than 99.78% (respectively 99.96% ) of the bonds are correctly predicted with B = 4500 (respectively B = 45000) configurations. For even larger values of the coupling constant, J = 4, more than 95 % (respectively 99%) of the bonds are correctly predicted with B = 4500 (respectively B = 45000) with the Mean-Field algorithm.

In Section VI D 3 we show that the probability of success increases in the ferromagnetic phase, when using our algorithm with a well-chosen threshold Θ rather than the simple Mean-Field procedure (Θ = 1), e.g. all neighborhoods are perfectly reconstructed for J = 3.2 (Q succ = 1).

Is thermalization relevant to the inference at low temperatures?

The diluted bidimensional grid, with a fraction ρ = .7 of non-zero bonds, undergoes a transition from a paramagnetic to a ferromagnetic phase at the value J crit (ρ = .7) ≃ 2.8 (vertical line in Fig. 25) in the infinite size limit [47,48]. In the ferromagnetic phase, J > J crit (ρ = 0.7), and for small bidimensional grids, two competing 'states' coexist: the 'down' state, where most spins are 0, and the up state, where most spins are equal to 1. The system 'jumps' from one state to the other, as shown by the time-dependence of the average activity,

µ(t) = 1 N N i=1 σ i (t) , ( 78
)
where t is the Monte Carlo time. Figure 26(a) shows that the two states are equally sampled on a 9 × 9 grid, with N A = 10, 000 single spin-flip attempts with the Metropolis rule in between two sampled configurations (the results of Fig. 25 on the 7 × 7 grid were obtained with the same value of N A ). To investigate the performance of the algorithm when the two states are not well sampled we have studied a 9×9 grid, with N A = 100 and N A = 1, 000. For N A = 100, fig. 26(b) shows that few transitions occur, and that the two states will likely not be weighted equally. On larger 12 × 12 grids no jump occurs, even with N A = 1, 000 spin-flip attempts (Fig. 26(c)). The values of the spin averages Left: Relative errors ǫp, ǫc as functions of the threshold (top) and comparison of reconstructed and data correlations at Θ * (bottom). Middle and Right: comparison of the inferred and true fields (top) and histograms of inferred couplings (bottom) for our cluster algorithm and the mean field procedure. Color code for the histograms: brown/gray: unreliable couplings (which also correspond to zero couplings in the true network), black : reliable couplings. The two sub-distributions are normalized separately.

p i will strongly vary between the three cases above: p i ≃ .5 in the mixed case (a), p i ≃ .7 in the particular realization (b) of the partially mixed case, and p i close to zero in case (c). Remarkably the probability of success of the algorithm is not sensitive to the nature of the mixing. Fig. 26(d) shows, indeed, that the reconstruction performances do not significantly decrease even in the partially and non mixed cases compared to the fully thermalized case.

Inference of the couplings and reconstruction of 1-and 2-point averages with the cluster expansion

It is harder to determine the values of the fields and of the couplings and to reconstruct the frequencies than to infer the structure of the interaction graph alone. To this aim the minimization of the MF entropy S MF is generally not sufficient, and the cluster expansion of S -S MF has to be carried out. In Fig. 27 we show the relative errors ǫ p and ǫ c in the reconstruction of one-and two-site frequencies as a function of the threshold Θ for the same 7 × 7 randomly diluted grid as in Section VI D 1. We also compare the fields and couplings inferred with our cluster algorithm (middle panels in Fig. 27, Θ * = .05) to the ones found with the simple MF procedure (right panels in Fig. 27, large value of Θ). Note that the small error done in the graph learning for J = 3.2 (P succ = 0.997) can be avoided when the threshold value is optimized for each data realization.

We show in Fig. 28 the performances of the algorithm in the case of poor mixing, when the two states are not equally sampled. For the particular realization corresponding to Fig. 28, the frequencies are p i ≃ 0.3 instead of p i ≃ 0.5. In spite of the poor mixing the inference of the fields and couplings is as accurate as in the case of well-mixed sampling. The difference between the fields corresponding to the apparent frequencies p i ≃ 0.3 and the true one (p i = .5) are, indeed, smaller than the statistical uncertainty on the fields due to the limited sampling (finite value of B). The reason is that, near a critical point, a small variation in the field is sufficient to produce a large change in the average values of the spins.

E. Erdös-Renyi random graphs

In this Section we report the results of our inference algorithm when applied to disordered Ising models on random graphs. The random networks are generated from the Erdös-Renyi ensemble, where M = d 2 N edges are drawn, uniformly and at random, between N points. Parameter d is the average degree of a vertex on the network.

Figure 29 shows the outcome of the algorithm when data are generated from an Erdös-Renyi model of connectivity d = 10. On the selected bonds (i, j) the couplings J ij were chosen uniformly at random in [-3; 3]. All other couplings J ij were set to zero, and the fields were h i = -1. Values of the parameters are such that the system is in the paramagnetic phase (in the thermodynamic limit). Panel A shows the inference with good sampling (the data are obtained by averaging over B = 10 6 Monte Carlo configurations), while Panel B shows the inference with poor sampling (B = 10 3 ). At Θ * the data are reconstructed within the expected statistical fluctuations and, correspondingly, couplings are found back within the statistical error bars δJ ij . In the case of a large sampling noise case B = 10 3 , the statistical fluctuations δJ ij are so large that most of the inferred couplings are unreliable. The inference of the complete network is thus not possible. The maximal size of clusters at θ * increases with the average degree (Section VI F); we find K max = 9 and K max = 7 for, respectively, B = 10 6 and B = 10 algorithm is able to reach ǫ c = 1 at large thresholds, with a small number of selected clusters, e.g. N clu < 1000 and K max = 7 for N = 100. The threshold Θ * for which ǫ c = 1 corresponds to the beginning of the plateau for ǫ J . For smaller thresholds ǫ c decreases and data are overfitted. The height of the plateau for ǫ J coincides with the calculated statistical error δJ; it scales as 1/ √ B and does not strongly depend on N .

F. Computational Time

For a given value Θ of the threshold the computational time can be estimated through time

≃ Kmax K=1 N clu (K) 2 K , (79)
where N clu (K) is the number of selected clusters of size K, and 2 K is the number of operations necessary to calculate exactly the partition function of a sub-system of size K. As the number of selected clusters depends on the interaction graph, the computational time is sensitive to the structure of the interaction graph, while it does not depend too much on the correlation length of the system. For instance, the number of processed clusters and the computational time for Erdös-Renyi graphs is larger for the connectivity d = 10 than for d = 5. Moreover, as the sampling noise increases (the value of B is made smaller), so does the threshold value Θ . As less precision is needed in the reconstructed frequencies and correlations, the size of the selected clusters is reduced. As a consequence, the computational time is reduced. Figure 31 illustrates this statement for Erdös-Renyi random graphs: the running time increases with the quality of the sampling, i.e. increases with the value of B. In some very noisy cases, however, large size clusters which are due only to the noise and do not reflect the interaction network can be processed. As an example, the number of clusters for the unidimensional Ising model with ξ = 9, B = 1000 is larger than the one for B = 4500. Another illustration is given by the diluted 7 × 7 grid with J = 3.2, which requires the processing of many clusters of large size (K max = 8).

VII. CONCLUSION AND PERSPECTIVES

In this paper, we have presented an adaptive cluster expansion to infer the interactions between a set of Ising variables from the measure of their equilibrium correlations. We have discussed the statistical mechanics of this expansion, and shown applications of the algorithm to artificial data generated using Ising models on unidimensional and bidimensional lattices, as well as on Erdos-Renyi random graphs.

We have in particular underlined the important conditions on the inverse problem that should be fulfilled for our algorithm to be efficient. The essential condition is that the inverse susceptibility, which determines the change of a coupling (or a field) resulting from a change in the data (1-or 2-spin frequencies) should be well-conditioned. We stress that this property is not incompatible with the presence of a long-range susceptibility. Hence, the inverse problem can be easy to handle even in the presence of long-range correlations. As far as our algorithm is concerned, this condition entails that the maximal size K clu of the clusters which need to be taken into account remains small even if the correlation length of the system is large.

The origin of this condition is that our algorithm builds up, by definition, an interaction network defining a wellconditioned Ising model. Indeed, in the absence of reference entropy (S 0 = 0), the cross-entropy S(p) is approximated through a sum of a restricted number of cluster-entropies, see (30). For sufficiently large thresholds Θ, most quadruplets of variables, say, i, j, k, l, do not appear in any selected cluster (of size K ≥ 4); hence, most of the entries (χ -1 ) ij,kl of the inverse susceptibility matrix entries vanish according to (55). In the presence of the reference entropy S 0 = S MF , χ -1 is not guaranteed to be sparse any more due to the contribution χ -1 0 = -∂S0 ∂p∂p . However, when a regularization is introduced, e.g. based on the norm L 1 [START_REF] Wallace | Statistical and Inductive Inference by Minimum Message Length[END_REF], the network of interactions (J 0 ) ij is highly diluted, and we expect χ -1 0 to be well-conditioned, too. Further investigations of this point would be very useful. According to the discussion of Section IV C inverse problems corresponding to Ising models on finite-dimensional lattices are well-conditioned in the perfect sampling limit. The introduction of a threshold over the minimal values of cluster-entropies allows us to force the inverse problem to be well-conditioned even in the presence of sampling noise. We have checked this statement on inverse problems corresponding to 'critical' Ising models. While the correlation length increases with the size of the system, the maximal size of the clusters, K clu , remains roughly constant. Therefore, the computational complexity of the algorithm increases only linearly with the system size.

An essential feature of inverse problems is that data are generally obtained from a finite sampling and, therefore, frequencies and pairwise correlations are plagued by sampling noise. Avoiding overfitting is a primary goal for an inference algorithm. This goal is achieved, in our algorithm, by the introduction of the threshold Θ. As a result most of the clusters are discarded, and in particular, those whose contributions would convey very little information about the true nature of the underlying interaction network. Fixing the threshold value such that the relative reconstruction errors ǫ p and ǫ c are of the order of one corresponds to the maximal accuracy allowed by the quality of the data. The cluster expansion introduced here differs from other classical cluster expansions, developed in the contexts of the theory of liquids and of computational physics. In particular we do not impose consistency equations for the marginal probabilities over the clusters. Our expansion scheme is simpler, and requires only the knowledge of the individual and pairwise frequencies of the variables in the cluster. Moreover, the cluster construction and selection rules prevents any combinatorial explosion of the computational time.

Several points would deserve further investigations. Among them the discussion of the convergence properties of the expansion, started in Section IV D, should be expanded and improved. A natural and interesting question is to ask how the series behaves when the packets of Fig. 12 start mixing, i.e. in the presence of a strong sampling noise. Another aspect which should be better understood is the influence of the construction rule. Our heuristic consists in merging two almost completely overlapping clusters of size K to build a new cluster of size K + 1 (provided its entropy is larger than Θ). This rule has a simple intuitive interpretation, compatible with the notion of interaction path, and attempts with other rules have been less fruitful. However, a deeper theoretical understanding and justification is clearly needed. Last of all, the a posteriori validation of the method relies on the use of a Monte Carlo simulation to calculate ǫ c and ǫ p . We have tested another procedure to avoid the use of a Monte Carlo calculation, based on a partial resummation of the cluster contributions corresponding to the free-energy (at fixed couplings and fields). This procedure, whose applicability goes beyond the inverse problem, will be detailed in a further publication.

where the denominator P(p, B, σ 2 ) (marginal likelihood) is simply the integral of the numerator over all real-valued fields h. Given p and B we plot I =log P(p, B, σ 2 )/B as a function of σ 2 . The general shape of I is shown in Fig. 32. The value of σ 2 minimizing I is the most likely to have generated the data, and should be chosen on Bayesian grounds. For more than one spin calculating the marginal likelihood would be difficult. We thus need an alternative way of obtaining the best value for σ 2 . The idea is to calculate I through a saddle-point method, and include the Gaussian corrections which turn out to be crucial. This approach is correct when the size of the data set is large. A straightforward calculation leads to

I ≃ log 1 + exp(h * ) -p h * + Γ 2 (h * ) 2 + 1 2B log 1 + 1 Γ exp(-h * ) (1 + exp(-h * )) 2 (A3)
where Γ = 1/(Bσ 2 ) and h * denotes the root of (1 + exp(-h * )) -1p + Γ h * = 0. I decreases from I(σ 2 = 0) = log 2 with a strong negative slope, dI/dσ 2 (0) ≃ -B p 2 , and increases as log σ 2 /(2B) for large values of the variance. Expression (A3) cannot be distinguished from the logarithm of the true marginal likelihood I shown in Fig. 32.

f) i j k i j k i j k j i k i k j j k i a) d) c) b) e)
FIG. 33: Leading diagrams to the order 5 (top) and 6 (bottom) in the connected correlation for the entropy of 3-clusters.

Case K = 3

For K = 3 the leading term to ∆S (i,j,k ) in powers of M ij was not derived analytically in [30]. Based on the studies of the unidimensional Ising model and the independent spin models (Appendix F), we find that the leading diagrams are diagrams (a), (b), (c) in Fig. 33 (bold diagrams in Fig. 3), whose sum is given by ∆S

(5) i,j,k = α (5) ijk (c ij ) 2 (c jk ) 2 c ki + α (5) jik (c ij ) 2 c jk (c ki ) 2 + α (5) jki c ij (c jk ) 2 (c ki ) 2 , (C3) with α (5) ijk = - (2p i -1) (2p k -1) 2 (p i ) 2 (1 -p i ) 2 (p j ) 2 (1 -p j ) 2 (p k ) 2 (1 -p k ) 2 . (C4) Note that α (5)
ijk differs from α

jik . We have also found the coefficients of the subsequent diagrams, of the order of M 6 . These diagrams are labelled by (d), (e), (f) in Fig. 33. Their total contribution to the cluster-entropy is

∆S (6) i,j,k = α (6) ijk (c ij ) 3 (c jk ) 3 + α (6) jik (c ij ) 3 (c ki ) 3 + α (6) jki (c jk ) 3 (c ki ) 3 (C5) with α (6) ijk = (2p i -1) (2p k -1) 3 (p i ) 2 (1 -p i ) 2 (p j ) 3 (1 -p j ) 3 (p k ) 2 (1 -p k ) 2 . (C6) 3. Generic case K ≥ 4
The above results for K = 3 are easily generalized to any value of the cluster size K ≥ 4. The diagrammatic expansion of a K-spin cluster includes all circuits where pairs of spins are linked together. Each diagram with (one or two) links between i l and i l+1 (l = 1, . . . , K -1) and (one or two) links between i 1 and i K gives

∆S (2K-1) i1,...i k = (-1) K 2 K l=1 (p i l ) 2 (1 -p i l ) 2 (2p i k-1 -1) (2p i k -1) (c i1,i2 ) 2 (c i2,i3 ) 2 . . . (c i k-2 ,i k-1 ) 2 c i k-1 ,i k + (2p i k-2 -1) (2p i k-1 -1) (c i1,i2 ) 2 (c i2,i3 ) 2 . . . c i k-2 ,i k-1 (c i k-1 ,i k ) 2 + . . . + (2p i1 -1) (2p i2 -1) c i1,i2 (c i2,i3 ) 2 . . . (c i k-2 ,i k-1 ) 2 (c i k-1 ,i k ) 2 . ( C7
)
At the next order in power of M ij , each diagram with three links between i l and i l+1 (l = 1, . . . , K -1) gives a contribution ∆S 1. Distribution of cluster-entropies for K = 2

(3K-3) i1,...i k = (-1) K-1 (2p i1 -1) (2p i k -1) 3(p i1 ) 2 (1 -p i1 ) 2 (p i k ) 2 (1 -p i k ) 2 K-1 l=2 (p i l ) 3 (1 -p i l ) 3 K-1 l=1 (c i l ,i l+1 ) 3 . ( C8
The distribution of the entropy of K = 2-clusters for a set of B = 10 6 configurations is shown in Fig. 35. To derive the analytical expression of the distribution in the N → ∞ limit, we use the small-correlation formula (C1) for ∆S (1,2) , and the fact that the distribution of the connected correlation is Gaussian. As a result, approximating α 1,2 with its average value α obtained by substituting p 1 and p 2 with p in (C2), we obtain

H IS (∆S (1,2) ) = exp - (∆S (1,2) ) 2/3 2 (cB ) 2 α 2/3 3 α 1/3 2π(c B ) 2 (∆S (1,2) ) 2/3 (E2)
This distribution is a stretched exponential at infinity, and diverges in zero. Its standard deviation is approximately equal to 1.2 10 -6 for N = 40. This value is about twice the largest cluster-entropy observed in Fig. 35 for one particular realization of the sampled configurations.

2. Distribution of the cluster-entropies for K = 3

The leading order contribution to the entropy of a 3-cluster is given by (C3). We want to calculate the distribution of ∆S (1,2,3) when the connected correlations c ij are random Gaussian variables, of zero mean and variance (c B ) 2 . We neglect the correlations between c 12 , c 13 , c 23 , which is legitimate for large B. Let us call x = ∆S (1,2,3) /(α(c B ) 5 ), with α given by (C4), and let P (x) be the probability density of x. Though we have not been able to find a closed expression for P (x), the asymptotics behavior of P for large or small arguments can characterized analytically. We now use the saddle-point method again, this time to estimate the integral on the left hand side of (E5). We obtain (E9) at large x. The distribution of the cluster entropies ∆S (1,2,3) thus follows a stretched exponential with exponent 2 5 . This decay is much slower than an exponential, and leads to large tails as can be seen from Fig. 9.

b. Small x behavior

In order for the rescaled entropy x to be small, at least one among the three correlations should be small according to (C3). Without restriction, we may assume that c 12 is the smallest of the three correlations. As c 12 appears once with power one, and twice with power two in (C3), we approximate x ≃ c 12 c 2 13 c 2 23 . The Mellin transform of P is, for negative λ,

(E10)

The largest pole is located in λ = -1 2 , and is of order 2. According to standard results on the inversion of Mellin transforms [START_REF] Flajolet | [END_REF], we obtain a precise characterization of the divergence of the probability density at small x,

P (x) ≃ C (-log x) √ x , ( E11
)
where C is a constant.

c. Typical value of x

The typical value of the x is defined through As z is a sum of independent random variables its Laplace transform is the product of their Laplace transforms,

x typ = exp
e -u z z = ∞ 0 dc 2 √ 2π e -c 2 /2-u/c 3 = λ 2π √ 2 G 30 03 λ 2 8 - -1 2 , 0, 0 3 , ( E17
)
where G is the Meijer-G function. we have calculated the integral (E16) using the Mathematica software. Some care must be taken for the numerical accuracy when z → 0. The outcome is log z z ≃ 2.09643. Putting all contributions together we obtain x typ ≃ 0.18. The corresponding values of ∆S (1,2,3) are 3.5 × 10 -15 for B = 10 6 , and 1.1 × 10 -12 for B = 10 5 , in good agreement with the numerical value, respectively, 3 × 10 -15 and 9 × 10 -13 .

d. Standard deviation of x

We can easily evaluate the variance of each of the three terms of the sum in (C3) as the product of the variances of the three terms in the product, based on the approximation that the connected correlations c ij are independent stochastic variables. We obtain

σ ∆S ijk = 3 √ 3(2p -1) 2 (c B ) 5 2 p 6 (1 -p) 6 = 3 √ 3 (2 p -1) 2 2 p (1 -p) B 5/2 . (E18)
With the values of N and p chosen in Fig. 35, we find that the standard deviation is of the order of 10 -13 for B = 10 6 , and 2 10 -11 for B = 10 5 , see Fig. 9.

3. Distribution of cluster-entropies for generic K ≥ 4

In general, for K ≥ 3, the leading contribution to ∆S (i1,i2,...iK ) (C7) contains the sum of K × (K -1)!/2 terms, each one being the product of K random variables, among which (K -1) are elevated to power two, and 1 is elevated to power 1. The factor K comes from the fact that there are K way of choosing the single link in the circuits with K spins. The factor (K -1)!/2 is the number of non equivalent circuits going through K spins. We define the rescaled entropy x through Consider now the Gaussian model over N real-values variables ϕ i , whose energy function is given by

E[ϕ] = 1 2 i,j G ij ϕ i ϕ j . ( F14
)
For this Gaussian model, the logarithm of the partition function is (up an irrelevant additional constant), log Z[G] = -1 2 log det G. By construction, model (F14) is the solution of the inverse Gaussian problem, with data: ϕ i = 0, ϕ i ϕ j = M ij . Hence, S 0 can be interpreted as the cross-entropy of Gaussian model (F14) under those data. A key feature of the Gaussian model above is that its interaction matrix G ij is tridiagonal. Only nearest neighbour variables are coupled to each other according to (F13). We conclude that the Gaussian model is a one-dimensional model. Consequently, it obeys a chain rule similar to (F5). This is the only requirement for the main conclusion of Section F 1 a to hold: in the cluster expansion of S 0 , the entropy of a K-cluster is simply equal to the entropy of the 2-cluster with the same extremities, multiplied by (-1) K-2 . As both the expansions of S and the one of S 0 enjoy this property, so does the expansion of S -S 0 .

We conclude this Section by the expression of the 2-cluster entropy ∆S (i,i+d) . In the presence of the reference entropy S 0 = S MF , we substract the following contribution to expression (F10), see (21),

(∆S 0 (i,i+d) = 1 2 log det 1 M i,i+d M i,i+d 1 . (F15)
Hence, function F (u) defined in (F11) should be substracted 1 2 log(1u 2 ). It is a simple check that F (u) -1 2 log(1u 2 ) = O(u 3 ), in agreement with scaling (46).

Examples and calculation of diagrammatic coefficients

We now show the histograms of cluster-entropies for K = 2 and K = 3 for specific choices of J, h. The averages p i and p ij were calculated exactly through formulas (F3) and (F4) (perfect sampling). Figure 36(a) shows the histogram of entropies for clusters of the type (i, i + d). Entropy values are discrete and labelled by the distance d. They range from 10 -2 (for nearest neighbours, distance d = 1) to values smaller than 10 -15 for 6 < d < 15. All entropies smaller than the numerical accuracy ≃ 10 -15 are put in the peak at the origin. Expanding F (u) to the lowest order in u (for S 0 = S MF ) we find the asymptotic formula for the 2-cluster entropy: ∆S i,i+d ≃ (2p -1) 2 6 p (1p) e -3 d/ξ , (F16) in agreement with (C1). We have verified that this formula is in very good agreement with the numerics as soon as d ≥ 4 for the parameters of Fig. 36. Figure 36(b) show the histogram of the entropies of 3-clusters (i, j, k). Let d = ki be the distance between the extremities. We observe that the entropies are gathered into peaks, and are exactly the opposite of the ones found in Fig. 36(a) as expected. Two differences are: for various values of (m, n). Let d = ji. For m = i, n = j, constraint (G1) gives

γ(1 -x 2d ) + 2β 2x -x 2d-1 -x 2d+1 + α 4x 2 -x 2d-2 -x 2d -x 2d+2 = 1 , (G2)
which should hold for all d ≥ 3. We deduce two coupled equations for the three unknown variables:

γ + 2 x + 1 x β + 4 x + 1 x 2 α = 0 , (G3) γ + 4 x β + 4 x 2 α = 1 . ( G4
)
For m = i + 1, n = j, constraint (G1) is equivalent to γ(xx 2d-1 ) + β 1 + 3x 2x 2d-2 -3x 2d + α 2x + 2x 3x 2d-3 -2x 2d-1x 2d+1 = 0 . (G5)

The d-dependent term in the equation above cancels by virtue of (G3). We are left with an additional equation over α, β, γ:

γ x + β (1 + 3x 2 ) + 2 x (1 + x 2 ) α = 0 . (G6)
By symmetry of the matrices χ, χ -1 , no new constraint is obtained when the values of m, n are further varied. Solving (G3), (G4), (G6) we obtain

α = x 2 (1 -x 2 ) 2 , β = - x(1 + x 2 ) (1 -x 2 ) 2 , γ = (1 + x 2 ) 2 (1 -x 2 ) 2 .
(G7)

The analysis of the other cases j = i + 2 (Fig. 37(c)) and j = i + 1 (Fig. 37(d)) can be done along the same lines. We do not write the calculations in details, and simply report the results. The case j = i + 2 is very similar to the previous case. There are 8 coefficients to be calculated, with three independent values, α ′ , β ′ , γ ′ . It turns out that

α ′ = α , β ′ = β , γ ′ = γ . (G8)
As for the last case, j = i+1, we call α ′′ the values of the entries of χ -1 with (k, l) = (i-1, j-1), (i-1, j+1), (i+1, j+1), β ′′ the values of the entries with (k, l) = (i-1, j), (i, j +1), and γ ′′ the diagonal element corresponding to (k, l) = (i, j).

After some elementary algebra, we find

α ′′ = α , β ′′ = β , γ ′′ = 1 + x 2 + x 4 (1 -x 2 ) 2 . ( G9
)
All those results are reported in (62).

[1] S.G. Brush, Rev. Mod. Phys. 39, 883 (1967).

P

  (σ) σ ip i + k<l J kl σ P (σ) σ k σ lp kl .

FIG. 2 :

 2 FIG.2: Decomposition of the cross-entropy S(p) for a system of 4 spins, indicated with different colors, as the sum of cluster contributions. Each cluster-entropy ∆SΓ(p) depends only on the one-and two-point frequencies of the variables in the cluster: it can be calculated in a recursive way, see main text. Dotted clusters are decomposed into a diagrammatic expansion in Fig.3.

FIG. 4 :

 4 FIG.4: Examples of contour paths for three different graphs. Spins are labelled by 1, 2, 3, 4 and first-neighbor interactions are represented by bold lines. The contour path is depicted with a dotted line. The contour length L, which can be calculated as the sum of distances along the contour path (dotted arcs have length 2) is indicated above each graph. Different clusters may have the same contour path and contour length. Left: (1, 3) and (1, 2, 3) have contour length L = 4, while (1, 2) has L = 2. Middle: clusters (1, 3), (1, 2, 3), (1, 2, 3, 4), (1, 3, 4) have the same contour path. Right: (1, 2, 3, 4) has length L = 6.

5 FIG. 5 :

 55 FIG.5: Cluster entropy contribution ∆SΓ for a 3 × 3 grid (top-right) with nearest-neighbour couplings J = 1.777 (in units of kBT ) as a function of the contour length LΓ of the shortest closed path on the lattice joining the spins in Γ. To illustrate cancellation effects, some ∆SΓ are labelled with the indices in Γ, see main text. The values of J and of the fields hi = -1

6 FIG. 7 :

 67 FIG.7: Histograms of ∆S (i,j,k) for the 1D-Ising (full distributions, p = .02, ξ = 1) and Independent Spin (dotted distribution) models, with N = 50 spins. Each histogram correspond to one random sample of B configurations. Impulses show the histogram for perfect sampling (B = ∞), with Dirac peaks located at log |∆SΓ| = -3d Γ ξ + Cst. The cut-off at small entropies comes from the finite value of N .

6 FIG. 8 :FIG. 9 :

 689 FIG.8: Same as Fig.7, after rescaling by the standard deviation σIS(B)(50) of the Independent Spin model. Note the linear scale of the x-axis. As a result of the presence of the interaction network, the Ising histograms H are asymmetric in ∆S → -∆S for large values of ∆S, while the IS distributions HIS are obviously symmetric when averaged over the realizations of B configurations (not shown).

FIG. 10 :FIG. 11 :

 1011 FIG.10: Fraction of elements larger than 10 -7 (in absolute value) for the susceptibility χ (a) and the inverse susceptibility χ -1 (b) matrices vs. strength J of the nearest-neighbour coupling. The sizes N of the grids are indicated. Data were obtained from exact numerations for sizes 3 × 3, 4 × 4, 5 × 5 (perfect sampling, black) and from Monte Carlo simulations for all sizes (one realization of B = 4500 configurations, red/gray). Periodic boundary conditions were used.

FIG. 12 :

 12 FIG.12:(a) Sum Σ of the cluster-entropies larger than Θ (in absolute value) for the nearest neighbour one-dimensional Ising model with ξ ≃ 8.97 ≫ ξc ≃ 4.33, and B = 10 5 configurations. The initial increase from zero, taking place at small Θ ≃ ∆S(1) ≃ 0.41, is not shown. (b) magnification of (a) in the range .029 < Θ < .033. Within the random sign model, the behavior of Σ within a packet is similar to a Brownian bridge, see main text.

FIG. 13 :

 13 FIG. 13: Frequences of the block length ℓ for the 1D-Ising model, with N = 30 spins, ξ = 8.96, and one set of B = 10 5 sampled configurations. The statistics takes into account the cluster-entropies used to draw Fig. 12(a) only.

  2 and if |∆SΓ(p)| < Θ, then LK ← LK U Γ end for end while Output: list L sel of selected clusters Algorithm 4 Estimates for the cross-entropy and for the interaction parameters Require: data p, list L sel of selected clusters Computation of cross-entropy S:

  where |Γ|, |Γ ′ | are the sizes of Γ, Γ ′ . {hi(p), Jij (p)} ← arg min J SIsing Γ [J|p] compute {(h0)i,Γ(p), (J0)ij,Γ(p)} using formula (22), with MΓ replacing M . end for compute {hi(p), Jij (p)} using formula (74) Output: S(p), {hi(p), Jij (p)}

FIG. 14 :FIG. 15 :

 1415 FIG. 14: Performance of the algorithm as a function of Θ for the Independent Spin model: (a) errors ǫp and ǫc; (b) crossentropy S The entropy of the model in absence of sampling noise is ≃ 4.182071; (c) size Kmax of largest clusters; (d) number N clu of clusters. Panels (e) and (f ) show the reconstructed pi and cij vs. their values in the data. Error bars on pi and cij are calculated through (77). Data were obtained by sampling B = 10 5 configurations of the N = 40 spins, with spin dependent means pi equal to the average activity of the neurons in [2]. The optimal threshold Θ * is already obtained with clusters of size K = 1.

FIG. 16 :

 16 FIG.16: Errors on the entropy (ǫS), the couplings (ǫJ ), and the fields (ǫ h ) vs. threshold Θ for the unidimensional Ising model with ξ = 1 (top) and ξ = 9 (bottom). Full lines correspond to large sampling noise (B = 10 5 for ξ = 1, B = 10 3 for ξ = 9), dashed line to weak sampling noise (B = 10 7 for ξ = 1, B = 10 5 for ξ = 9). The size is N = 30.

FIG. 17 :

 17 FIG.17: Performance of the algorithm on the unidimensional Ising model with ξ = 1, B = 10 5 (full-bold line) B = 10 7 (dashed line). For both values of B the optimal threshold 0.0007 < Θ * < 0.003 is reached with clusters of size K = 3 and length L = 6. The reconstruction of data pi and cij is shown only for the B = 10 5 case (large sampling noise). Error bars on pi and cij are computed from (77).

FIG. 18 :

 18 FIG. 18: Performance of the algorithm on the unidimensional Ising model with ξ = 9, B = 10 3 (full-bold line) and B = 10 5 (dashed line). For the largest sampling noise case (B = 10 3 ) at the threshold Θ * ≃ .104, 89 clusters of size K = 2, 92 clusters of size K = 3, 35 clusters of size K = 4, and 2 clusters of size K = 5 are selected. The reconstruction of data pi and cij is shown for the B = 10 3 case. Error bars on pi and cij are computed from (77).

FIG. 19 :

 19 FIG. 19: Histograms of the inferred couplings J inf ij for the unidimensional Ising model with ξ = 1. Couplings equal to J = 4 and J = 0 in the model are shown in, respectively, black and gray. Gray couplings are unreliable, as they differ from zero by less than three standard deviations. Values of B and Θ are indicated on the figure. The bin width is ∆J = .08.

FIG. 20 :

 20 FIG.20: Same as Fig.19, but with ξ = 9 instead of ξ = 1. Couplings equal to J = 6 and J = 0 in the model are shown in, respectively, black and gray. Gray couplings are unreliable, as they differ from zero by less than three standard deviations.

FIG. 21 :FIG. 22 :FIG. 23 :

 212223 FIG.21: Errors ǫs, ǫJ , ǫ h , ǫc, ǫp, size Kmax, and number K clu of clusters vs. Θ for a 3 × 3 grid with J = 1.778. The dashed line corresponds to B = 4500 sampled configurations, the full line to B = 10 7 ; the perfect sampling curves are shown with dotted lines in the top three panels (same data as Fig.6). The values of Θ * are shown with vertical lines: Θ * = .0017 (full, B = 10 7 ), and .007 (dashed, B = 4500). At Θ * 21 clusters (B = 4500) and 49 clusters (B = 10 7 ) are selected .

FIG. 24 :

 24 FIG. 24: Histograms of inferred couplings for bidimensional grids of different sizes. Couplings equal to J = 1.778 and 0 in the model are shown in, respectively, black and gray; the integral of each sub-distribution is normalized to one by hand. The values of Θ * for each size are shown in Fig. 23.

FIG. 26 :FIG. 27 :

 2627 FIG.26: (a,b,c). Average activity µ (78) of the B = 4500 sampled configurations for three grid sizes and numbers of spin-flip attempts in between two samplings. (d) Probability of success for J > 2.8, for the same data as in (a,b,c) and for the 7 × 7 grid of Fig.25.

FIG. 28 :

 28 FIG. 28: Performance of the inference algorithm in the case of poor mixing, for a 7 × 7 randomly diluted grid, with J = 3.2, B = 4500, NA = 10 2 . Top & left: relative errors ǫp, ǫc as functions of the threshold Θ. Bottom & left: reconstructed pi and cij versus their data values. Right: Inferred fields h inf i vs. their true values for Θ * = 0.05 (top) and histogram of the inferred couplings J inf (bottom). As usual, unreliable couplings, which correspond to zero couplings in the true network, are depicted in gray, and reliable couplings in black.

3 .FIG. 29 :

 329 FIG. 29: Outcome of the inference algorithm for an Erdös-Renyi random graph with N = 50 spins, connectivity d = 10, and with B = 10 6 (a) and B = 10 3 (b) sampled configurations. For each value of B we show the errors ǫp, ǫc vs. θ, the inferred vs. data values of the correlations cij , and of the couplings Jij . A few large error bars over Jij (calculated from χ -1 ) are shown. Values of Jij were chosen uniformly at random in [-3; 3], and fields were set to hi = -1.

FIG. 30 : 2 j

 302 FIG. 30: Outcome of the inference algorithm for an Erdös-Renyi random graph with N spins, and connectivity d = 5. Values of Jij were chosen uniformly at random in [-4; 4], and fields were set to hi = -1 2 j( =i) Jij . Left: error ǫJ on the inferred couplings as a function of Θ for N = 50 and B = 1000, 4500, 10000 configurations. Middle: error ǫJ vs. Θ for N = 50, 100, 200 and for B = 4500. Right: ǫc (top) and K clu (middle) as functions of Θ; the inferred couplings J inf are compared to their true values in the bottom panel for N = 100 and threshold Θ * .

6 IFIG. 32 :

 632 FIG. 32: Logarithm of the marginal likelihood (with a minus sign, and divided by the size B of the data set) versus variance σ 2 of the prior distribution of the field. Parameters are B = 100, p = .02.

FIG. 35 :

 35 FIG.35: Theoretical (red) and numerical (black) histograms HIS for the entropies ∆S of 2-spin clusters in a system of independent spins and B = 10 6 configurations. Simulations were done with N = 40 spins, with heterogeneous pi ≃ .0238, see caption of Fig.34.

σ

  ∆S (1,2) = √ 15 α (c B ) 3 = √ 15 (2p -1) 2 6 p(1p) B 3/2 (E3)For B = 10 6 and p = 0.0238 we obtain that the standard deviation is ≃ 2.7 10 -8 . Distribution (E2) is compared to the histogram obtains from numerics in Fig.35. The standard deviation and the distribution at small entropies are in good agreement. Large values of the correlations (E1) give rise to isolated values of ∆S (1,2) , of the order of ∆S MAX (1,2) ≃ (4 log N ) 3/2 (∆S (i,j) ) 2 1/2 , (E4)

dc 12

 12 dc 13 dc 23 e -F (c12,c13,c23) (E5) where F (c 12 , c 13 , c 23 ) = -P (x) at large x can be studied by considering large values of λ. We expect the dominant contribution to the multiple integral on the right hand side of (E5) to come from large correlations. The location of the main contribution to the integral is the value of (c 12 , c 13 , c 23 ) which maximizes F . As F is invariant under any permutation of its arguments, we look for a maximum where c 12 = c 13 = c 23 ≡ c * . A straigthforward calculation shows that c * (λ

  max x log P (x) + λ log x = F * (λ) ,(E8)which is true when λ is very large. Hence, F * (λ) is the Legendre transform of log P (x). Solving (E8) gives log P (x) ≃ -

= 6 ∞ 0 dc 2 log c √ 2π e -c 2 / 2 +

 6022 less sensitive than the average value of x to the presence of the long tails in P (x) at large x. We write x = (c 12 c 13 c 23 ) the logarithm, and averaging over the correlation, we obtain the following expression for the average value of the logarithm of x, log x log z z .(E14)The integral over c in the above equation can be calculated numerically, with a value ≃ -.63518. To calculate the average value of log z, we first use the identity log z =

xFIG. 36 :

 36 FIG. 36: Histograms of the entropies for clusters of size 2 (a) and 3 (b); in the latter case, entropies are negative. Data are generated from the unidimensional Ising model (F1) with N = 30 spins, and parameters J = 4 and h = -6. Each peak is labelled by the distance d between the extremities of the clusters. The reference entropy is S0 = SMF .

FIG. 37 :

 37 FIG.37:Half grid representing the index (i, j) of the entries of the inverse susceptibility matrix, with i < j (a). Black circles locate the nearest-neighbours and the diagonally opposed sites (k, l) of (i, j) (cross), with i = 4 and j = 9 (b), 6 (c), 5 (d).

  Computational time of the cluster algorithm at threshold Θ * for the different examples shown Section VI. The computational time grows with the number of processed clusters, which depends on the structure of the interaction graph and on the number of sampled configurations more than the sole number of variables, N . For a fixed maximal size of clusters, Kmax, the running time is roughly proportional to the number of clusters. Unless explicitly stated otherwise, the sampling is realized in good mixing conditions. Times were measured on one core of a 2.8 GHz Intel Core 2 Quad desktop computer.

	Grid 7x7, ρ=0.7, J=3.2, poor mixing, B=4500 Grid 7x7, ρ=0.7, J=3.2, B=4500 Grid 20x20, ρ=1, J=1.78, B=4500
	Grid 12x12, J=1.78, B=4500 Grid 7x7, ρ=1, J=1.78, B=4500
	Ising 1d, N=30, ξ=9, B=10 5
	Ising 1d, N=30, ξ=9, B=1000
	Erdos Renyi, N=50, d=10, B=10	6
	Erdos Renyi, N=50, d=10, B=10	3
	Erdos Renyi, N=50, d=5, B=4500
	Erdos Renyi, N=100, d=5, B=4500
	FIG. 31:

j( =i) J ij , in such a way that the corresponding fields in spin variable ±1 vanish. This is an example of a smaller connectivity system in the spin-glass phase. We have studied the performance of the algorithm as a function of the threshold Θ, by varying the system size N from 50 to 200 and the number of sampled configurations from B = 1000 to 10000. The

j =i Jij .

Acknowledgements: We are grateful to J. Barton, J. Lebowitz, E. Speer for very useful and stimulating discussions, in particular regarding the correspondence between the inverse susceptibility and the direct correlation functions and the practical implementation of the inference algorithm. We thank E. Aurell for pointing to us the difference between P and Q, see Section VI D 1.

Appendix A: Optimal choice for the regularization parameter γ

In this Appendix we discuss how the optimal value for the parameter γ in the L 2 -regularization (13) can be determined. As explained in Section II B, the regularization term can be interpreted as a Gaussian prior P 0 over the couplings. Let us call σ 2 the variance of this prior. Parameters γ and 1/σ 2 are related through

where we have assumed that the 1-site frequencies p i are uniformly equal to p. To calculate the optimal value for γ, or, equivalently, for σ 2 , we start with the case of a single spin for the sake of simplicity, and then turn to the general case of more than one spin.

a. Case of N = 1 spin

For a unique spin subjected to a field h the likelihood of the set of sampled spin values, {σ τ }, is P h [σ] = exp(B p h)/(1 + e h ) B . Here p denotes the average value of the spin over the sampled configurations (2). We obtain the a posteriori probability [START_REF] Ravikumar | [END_REF] for the field h given the frequency p,

The above saddle-point approach can be generalized to any number N of spins, with the result

where S Ising was defined in (7) and χ is the N + 1 2 N (N -1) = 1 2 N (N + 1)-dimensional Hessian matrix composed of the second derivatives of S Ising with respect to the couplings and fields [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF]. In principle χ could be diagonalized and the expression (A4) calculated. However this task would be time-consuming. As we have seen in the previous subsection we expect I not to increase too quickly with σ 2 (for not too small variances) and approximate calculations of I can be done under some data-dependent hypothesis. We now give an example of such an approximation, valid in the case of multi-electrode recordings of neural cell populations.

A simplification arises when the number B of configurations and the frequency p are such that: (a) each spin i is active (= 1) in a number of configurations much larger than 1 and much smaller than B i.e. 1 ≪ B × p ≪ B; (b) the number n 2 of pairs of spins that are never active together is much larger than one and much smaller than N (N -1) 2 . These assumptions are generically true for the applications to neurobiological data. For instance, the recording of the activity of N = 40 salamander retinal ganglion cells in [2] fulfills conditions (a) and (b) for a binning time ∆t = 5 msec: a cell i firing at least once in a time-bin corresponds to σ i = 1, while a silent cell is indicated by σ i = 0. More precisely: (a) the least and most active neurons respectively fire 891 and 17,163 times (among B = 636, 000 configurations); (b) n 2 = 34 pairs of cells (among 780 pairs) are never active together.

Condition (a) allows us to omit the presence of Γ in the calculation of the fields, h i ≃ log p i , to the first order of a large (negative) field expansion. Condition (b) forces us to introduce a non-zero Γ to calculate the couplings, with the result that interactions between pairs i, j of cells not active together are equal to

Finally we obtain the asymptotic scaling of the entropy when Γ → 0,

We are now left with the calculation of the determinant in (A4). From assumption (b) the number of pairs of neurons not spiking together is small with respect to N 2 , meaning that most of the eigenvalues λ a of the Hessian matrix of S Ising are non zero. Hence,

Putting both contributions to I together we get

The optimal value for the variance σ 2 is the root of

We finally deduce the optimal variance

For the data described above we find σ 2 ≃ 8.

Appendix B: Expression of the entropy of clusters with size K = 3

In this Appendix, we give the analytical expression for the entropy of a cluster with K = 3 spins. Using this expression instead of minimizing the cross-entropy (7) offers a valuable computational speed-up as there are O(N 3 ) clusters of size K = 3. We start with the definition of the entropy P (σ):

We then replace the probabilities P (σ 1 , σ 2 , σ 3 ) of the eight configurations of the three spins above with their expressions in terms of the probabilities {p i , p kl } in the data, and of the probability p 123 that the three spins are equal to 1:

Upon substitution of p 123 in (B1) we obtain the desired cross-entropy S 3 , as a function of the three average values p i and the three two-point averages p kl . The expression of the cluster-entropy is given by,

according to (28). The expressions of the cluster-entropies for one and two spins are given by, respectively, ( 23) and (25). Similarly, one obtains the expressions for the contributions of the 3-spin cluster to the values of the interactions parameters by differentiating ∆S with respect to the p i 's and the p kl 's.

Appendix C: Leading diagrammatic contributions to small cluster-entropies

We analyze the dominant diagrams contributing to the cluster-entropies for the various values of the cluster sizes, K, in the limit of small connected correlations c kl .

The entropy ∆S (i,j) of a 2-spin cluster is the sum of all diagrammatic contributions containing two spins and an arbitrary number of links between them, corresponding to the power of the expansion parameter 3) and Section III C 1. For small values of M ij the largest contribution to ∆S (i,j) is the one with three links (cubic power of M ij ), if the reference entropy S 0 = S MF removes the two-link loop diagram. The entropy contribution of this diagram was computed in [30], with the result

(C2)

The superscript 3 refers to the power of the connected correlation. The system is made of N = 40 independent spins, with the same pi as in the neural data of Ref. [2]; the average value of the pi's is p = 0.0238. Theoretical predictions correspond to (C3) and (C5). The number of sampled configurations is B = 10 6 .

Appendix D: Critical correlation length ξc for the absolute convergence

In this Appendix, we briefly explain why the cluster-entropy series is absolutely convergent if and only if the correlation length ξ is smaller than

Here, Ω = 2 when the reference entropy is S 0 = 0, and Ω = 3 when S 0 = S MF . Parameter v denotes the number of neighbours of a site on the lattice, supposed to be uniform. For instance, v = 2D on a hypercubic lattice in dimension D ≥ 1. Consider a set of K distinct points on the lattice. Let N (L) be the number of closed paths of length L visiting all K points. We obviously have N (L) ≤ v L . Hence, the series

is convergent if ξ < ξ c . Reciprocally, let L 0 be the length of the shortest closed path C 0 encircling the K points. A closed path of length L 1 + L 0 can be built from C 0 by attaching a closed loop of length L 1 to any one of the sites in C 0 . Hence, for L ≥ L 0 + 2, N (L) ≥ L 0 v L-L0 . We deduce that the series (D2) is divergent if ξ > ξ c .

Appendix E: Distribution of cluster-entropies for the Independent Spin model

We generate B configurations of N independent spins σ i . Spin i is equal to 1 with probability p and to zero with probability 1p (for simplicity we assume here that all the frequencies p i are equal to the same value p). The empirical connected correlations c ij computed from the B sampled configurations of spins are generally non zero. The marginal distribution of c ij is a normal law, with zero mean and standard deviation [START_REF]0,1, the couplings (Jij ) and fields (hi) are[END_REF]. The largest values of the correlations are, for a system with N spins, of the order of

according to extreme value theory. We compare in Fig. 34 formulas (C3) and (C5) for, respectively, the cluster-entropies ∆S ij and ∆S ijk with numerics carried out from randomly sampled configurations. Each pair (i, j) (Fig. 34(a)) and triplet (i, j, k) (Fig. 34(b)) define a point, whose coordinates are the numerical and theoretical values of the entropy corresponding to the pair-or triplet-cluster. The agreement, for B = 10 6 sampled configurations, is excellent due to the small value of c B ≃ 2 10 -5 .

The approach followed in Section E 2 to calculate the asymptotic behaviour of the probability density P of x for K = 3 can be extended without difficulty to any value of K > 3. We find that P (x) diverges when x → 0, with

where C is a constant. Hence the shape of the distribution of x is, up to logarithmic terms, independent of K. On the contrary, the tail of the distribution for large x is very sensitive to K,

As in the K = 3 case, the distribution of the cluster entropies ∆S follows a stretched exponential. The exponent of the stretched exponential decreases with K. The variance of the distribution can be easily evaluated, with the result

Appendix F: Properties of the cluster-entropies of the one-dimensional Ising model

Consider the one-dimensional Ising model with nearest-neighbour couplings and periodic boundary conditions. The Hamiltonian of the model is

where the spins σ i take 0,1 values. The parameters of the model are the N identical fields h i = h, the N couplings J i,i+1 = J between neighbours and the remaining N × (N -3)/2 zero couplings J i,j = 0 between non neighbours. We recall a few elementary facts about the model. The transfer matrix is

The eigenvalues are λ ± = 1 2 e J+h + 1 ± (e J+h -1) 2 + 4 e h , and the two components of the eigenvectors are, respectively, v ± (1) = -(1λ ± )/ e h + (1λ ± ) 2 and v ± (2) = e h/2 / e h + (1λ ± ) 2 . The probability that a spin is up is given by, in the N → ∞ limit,

and the connected correlation at distance d is

where the correlation length is given by ξ = -1/ log(λ -/λ + ).

Calculation of the cluster-entropies and cancellation property

In this Section, we show the exact cancellation property between the entropies of clusters with different sizes discussed in Section III C 2. We will see that this property is a direct consequence of the existence of a unique interaction path along the unidimensional chain. We first consider the case where the reference entropy is zero. Let Γ = (i 1 , i 2 , . . . , i K ) be a cluster of size K, with i 1 < i 2 < . . . < i K . Due to the unidimensional nature of the interactions, the Gibbs distribution over the K-spin configurations σ obeys the chain rule,

where P (•, •) and P (•|•) denote, respectively, joint and conditional probabilites. Inserting the above formula into expression [START_REF] Jaynes | Proc. IEEE[END_REF] for the cross-entropy, we obtain

Each variable σ i l , with l = 2, . . . , K -1, appears three times in (F6), which explains the presence of three fields h with the same index i l . After optimization over J = ({J i l ,i l+1 }, {h → i l }, {h ← i l }, {h 0 i l }) all these fields are equal, and we obtain

Hence the cross-entroy S(p) is the sum of the 1-cluster entropies and of the entropies of the 2-clusters made of adjacent sites. None of the other cluster-entropies appear, which proves that they cancel each other. To illustrate the cancellation mechanism, consider the case K = 3. According to (F7),

Comparing with ( 28) we obtain

which shows that the entropy of a 3-cluster and the one of a 2-cluster with the same extremities i 1 , i 3 are opposite to one another. By a recursive applications of (F7) this result can be immediately generalized to higher values of K.

The entropy of a K-cluster is simply the entropy of the 2-cluster with the same extremities, multiplied by (-1) K-2 . Hence, identity ( 44) is established. According to formula (F4) for the connected correlation, the entropy of a two-site cluster is a function of the distance d between the two sites:

where

To obtain the expression (F11) for F , we have used formula (26) for the 2-spin cluster-entropy, with p 1 = p 2 = p and p 12 = p 2 + c 12 , where the correlation c 12 is given by (F4). Note that F (u) = O(u 2 ) for small u, in agreement with scaling [START_REF] Johnstone | Proc. ICM 2006[END_REF].

b. Case S0 = SMF

We now introduce the reference entropy S 0 = S MF . The matrix M defined in (21) has elements

The inverse of M , G = M -1 , is a tridiagonal matrix, whose non zero elements are

• The peak in d = 1 is not present because the minimal distance between three spins is d = 2. The largest 3-spins entropy thus corresponds to triplets of the type (i, i + 1, i + 2).

• The height of the peak (number of clusters) corresponding to distance d is (d -1) N . The degeneracy (d -1) is the number of ways of choosing the location of site i 2 in between i 1 and i 3 .

We now show how the value of the cluster entropy can be found back from the leading terms in the diagrammatic expansion calculated in Section C 2. Let us call d ′ = ji < d the distance between the first two sites in the cluster. For each diagram in Fig. 33 we give in Table F 2 the sum of the distances of its links, i.e. the power of exp(-1/ξ).

Interestingly, the lowest total distances are found in diagrams a) and d), while the latter diagram is of a higher power [START_REF] Balakrishnan | [END_REF] in terms of the correlated function than the former [START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF]. Hence, contrary to the case of independent spins (Section E), diagrams a) and d) give the dominant contributions to the entropy. Summing the contributions of a) and d) we find

To derive the coefficients α (5) and α (6) , we impose that ∆S (i,j,k) is the opposite of (F16). We deduce that α (5) and α (6) are given by, respectively, (C4) and (C6). The exact cancellation property discussed above has important consequences for the inferred fields and couplings. Consider for instance the coupling J i,i+2 , which vanishes in the 1D-Ising model with nearest-neighbour interactions (F1). As the connected correlation c i,i+2 is not equal to zero, a contribution to the coupling will be collected from the cluster (i, i + 2) itself, equal to

Other contributions will come from larger clusters. For instance the cluster (i, i + 1, i + 2) will give an additional

The sum of the two contributions above vanishes due to the cancellation property. It can be checked that the contributions coming from all the other clusters vanish, too, which makes the coupling J i,i+2 = 0 as it should.

Appendix G: Inverse susceptibility matrix for the unidimensional Ising model

Hereafter, we want to invert the matrix χ, whose elements are given in (61). The matrix is of dimension 1 2 N (N -1), and each element is labelled by two indices (i, j) and (k, l), with i < j and k < l. Each index (i, j) can be represented by a site of coordinates i and j on the half-grid of Fig. 37(a). We now show that the non-zero entries of the inverse susceptibility matrix, χ -1 ij,kl , are in one-to-one correspondence with the sites (i, j) and (k, l) that are either identical, or nearest neighbours, or diagonally opposed on the elementary mesh of the half-grid (Fig. 37(b,c,d)). Depending on the value of the difference ji, the number of those sites is equal to 9, 8, or 6.

We start with the case ji ≥ 3 (Fig. 37(b)). By symmetry, the nine unknown matrix elements χ -1 ij,kl take only three independent values, denoted by γ for (k, l) = (i, j), β for (k, l) and (i, j) nearest neighbours, and α for (k, l) = (i ± 1, j ± 1). We now write the matrix inversion identity, k<l χ -1 ij,kl χ kl,mn = δ i,m δ j,n , (G1)

[53] In mean-field spin-glasses, as the couplings scale as the inverse square root of the number N of spins, only loop diagrams have non-zero weights in the thermodynamical limit.

[54] A further theoretical argument supporting the existence of the cancellation property is, in the case of perfect sampling, the fact that the entropy S must be extensive in N . As S is the sum of ∼ 2 N cluster entropies, those contributions must compensate each other.

[55] If p = 1 2 , ∆SIS is of the order of B -K [56] This statement is widely believed to be true in the theory of liquids literature. The fast decay of the inverse susceptibility, χ -1 (r), or, equivalently, of the direct pair correlation, g(r), is used to approximately close the hierarchy of correlation functions. The Percus-Yevik closure scheme, which gives an accurate equation of state for liquids of hard spheres, assumes that the inverse susceptibility vanishes above the interaction range of the potential (diameter of a particle).

[57] We have verified that the dependence on N is weaker in the case of periodic boundary conditions.

[58] The numerical experiments of [16] were done with ±1 spins and with coupling parameter θ and field ν = 0; in the present work where spins are equal to 0,1, the corresponding couplings and fields are: J = 4θ and hi = -1