S Cocco

R Monasson

V Sessak

High-Dimensional Inference with the generalized Hopfield Model: Principal Component Analysis and Corrections

We consider the problem of inferring the interactions between a set of N binary variables from the knowledge of their frequencies and pairwise correlations. The inference framework is based on the Hopfield model, a special case of the Ising model where the interaction matrix is defined through a set of patterns in the variable space, and is of rank much smaller than N . We show that Maximum Likelihood inference is deeply related to Principal Component Analysis when the amplitude of the pattern components, ξ, is negligible compared to √ N . Using techniques from statistical mechanics, we calculate the corrections to the patterns to the first order in ξ/ √ N . We stress that it is important to generalize the Hopfield model and include both attractive and repulsive patterns, to correctly infer networks with sparse and strong interactions. We present a simple geometrical criterion to decide how many attractive and repulsive patterns should be considered as a function of the sampling noise. We moreover discuss how many sampled configurations are required for a good inference, as a function of the system size N and of the amplitude ξ. The inference approach is illustrated on synthetic and biological data.

I. INTRODUCTION

Understanding the patterns of correlations between the components of complex systems is a fundamental issue in various scientific fields, ranging from neurobiology to genomic, from finance to sociology, ... A recurrent problem is to distinguish between direct correlations, produced by physiological or functional interactions between the components, and network correlations, which are mediated by other, third-party components. Various approaches have been proposed to infer interactions from correlations, exploiting concepts related to statistical dimensional reduction [START_REF] Jolliffe | Principal Component Analysis[END_REF], causality [START_REF] Seth | [END_REF], the maximum entropy principle [START_REF] Jaynes | Proc. IEEE[END_REF], Markov random fields [START_REF] Hastie | Elements of Statistical Learning: Data Mining, Inference and Prediction[END_REF] ... A major practical and theoretical difficulty in doing so is the paucity and the quality of data: reliable analysis should be able to unveil real patterns of interactions, even if measures are affected by under-or noisy sampling. The size of the interaction network can be comparable to or larger than the number of data, a situation referred to as highdimensional inference.

The purpose of the present work is to establish a quantitative correspondence between two of those approaches, namely the inference of Boltzmann Machines (also called Ising model in statistical physics and undirected graphical models for discrete variables in statistical inference [START_REF] Hastie | Elements of Statistical Learning: Data Mining, Inference and Prediction[END_REF]) and Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF]. Inverse Boltzmann Machines (BM) are a mathematically well-founded but computationally challenging approach to infer interactions from correlations. Our scope is to find the interactions among a set of N variables σ = {σ 1 , σ 2 , . . . , σ N }. For simplicity, we consider variables σ i taking binary values ±1 only; the discussion below can be easily extended to the case of a larger number of values, e.g. to genomics where nucleotides are encoded by four-letter symbols, or to proteomics where amino-acids can take twenty values. Assume that the average values of the variables, m i = σ i , and the pairwise correlations, c ij = σ i σ j are measured, for instance, through the sampling of, say, B configurations σ b , b = 1, . . . , B. Solving the inverse BM problem consists in finding the set of interactions, J ij , and of local fields, h i , defining an Ising model, such that the equilibrium magnetizations and pairwise correlations coincide with, respectively, m i and c ij . Many procedures have been designed to tackle this inverse problem, including learning algorithms [START_REF] Ackley | [END_REF], advanced mean-field techniques [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF][START_REF] Roudi | [END_REF], message-passing procedures [8,9], cluster expansions [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF][START_REF] Cocco | [END_REF], graphical lasso [START_REF] Hastie | Elements of Statistical Learning: Data Mining, Inference and Prediction[END_REF] and its variants [12]. The performance (accuracy, running time) of those procedures depend on the structure of the underlying interaction network and on the quality of the sampling, i.e. how large B is.

Principal Component Analysis (PCA) is a widely popular tool in statistics to analyze the correlation structure of a set of variables σ = {σ 1 , σ 2 , . . . , σ N }. The principle of PCA is simple. One starts with the correlation matrix,

Γ ij = c ij -m i m j (1 -m 2 i ) (1 -m 2 j ) , (1)
which expresses the covariance between variables σ i and σ j , rescaled by the product of the expected fluctuations of the variables taken separately. Γ is then diagonalized. The projections of σ along the top eigenmodes (associated to the largest eigenvalues of Γ) identify the uncorrelated variables which contribute most to the total variance. If a few, say, p (≪ N ), eigenvalues are notably larger than the remaining ones PCA achieves an important dimensional reduction. The determination of the number p of components to be retained is a delicate issue. It may be done by comparing the spectrum of Γ to the Marcenko-Pastur (MP) spectrum for the null hypothesis, that is, for the correlation matrix calculated from the sampling of B configurations of N independent variables [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]. Generally those two spectra coincide when N is large, except for some large or small eigenvalues of Γ, retained as the relevant components.

The advantages of PCA are multiple, which explains its success. The method is very versatile and fast as it only requires to diagonalize the correlation matrix, which can be achieved in a time polynomial in the size N of the problem. In addition, PCA may be extended to incorporate prior information about the components, which is particularly helpful for processing noisy data. An illustration is sparse PCA, which looks for principal components with many vanishing entries [START_REF] Aspremont | [END_REF].

In this paper we present a conceptual and practical framework which encompasses BM and PCA in a controlled way. We show that PCA, with appropriate modifications, can be used to infer BM and discuss in detail the amount of data necessary to do so. Our framework is based on an extension of a celebrated model of statistical mechanics, the Hopfield model [15]. The Hopfield model was originally introduced to model auto-associative memories, and relies on the notion of patterns [START_REF] Amit | Modelling Brain Function: the World of Attractor Neural Networks[END_REF]. Informally speaking, a pattern ξ = (ξ 1 , . . . , ξ N ) defines an attractive direction in the N -dimensional space of the variable configurations, i.e. a direction along which σ has a tendency to align. The norm of ξ characterizes the strength of the attraction. While having only attractive patterns makes sense for auto-associative memories, it is an unnecessary assumption in the context of BM. We therefore generalize the Hopfield model by including repulsive patterns ξ, that is, directions in the N -dimensional space which σ tends to be orthogonal to [START_REF] Nokura | [END_REF]. From a technical point of view, the generalized Hopfield model with p attractive patterns and p repulsive patterns is simply a particular case of BM with an interaction matrix, J, of rank equal to p + p. If one knows a priori that the rank of the true J is indeed small, i.e. p + p ≪ N , using the generalized Hopfield model rather than a generic BM allows one to infer much less parameters and to avoid overfitting in the presence of noisy data.

We first consider the case where the components ξ i and ξi are very small compared to √ N . In this limit case we show that Maximum Likelihood (ML) inference with the generalized Hopfield model is closely related to PCA. The attractive patterns are in one-to-one correspondence with the largest components of the correlation matrix, while the repulsive patterns correspond to the smallest components, which are normally discarded by PCA. When all patterns are selected (p + p = N ) inference with the generalized Hopfield model is equivalent to the meanfield approximation [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF]. Retaining only few significative components helps, in principle, to remove noise from the data. We present a simple geometrical criterion to decide in practice how many attractive and repulsive patterns should be considered. We also address the question of how many samples (B) are required for the inference to be meaningful. We calculate the error bars over the patterns due to the the finite sampling. We then analyze the case where the data are sampled from a generalized Hopfield model, and inference amounts to learn the patterns of that model. When the system size, N , and the number of samples, B, are both sent to infinity with a fixed ratio, α = B N , there is a critical value of the ratio, α c , below which learning is not possible. The value of α c depends on the amplitude of the pattern components. This transition corresponds to the retarded learning phenomenon discovered in the context of supervised learning with continuous variables and rigorously studied in random matrix and probability theories, see [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF][START_REF] Engel | Statistical Mechanics of Learning[END_REF][START_REF] Johnstone | Proc. ICM 2006[END_REF] for reviews. We validate our findings on synthetic data generated from various Ising models with known interactions, and present applications to neurobiological and proteomic data.

In the case of a small system size, N , or of very strong components, ξ i , ξi , the ML patterns do not coincide with the components identified by PCA. We make use of techniques from the statistical mechanics of disordered systems originally intended to calculate averages over ensembles of matrices to compute the likelihood to the second order in powers of ξi √ N for a given correlation matrix. We give explicit expressions for the ML patterns in terms of non-linear combinations of the eigenvalues and eigenvectors of the correlation matrix. These corrections are validated on synthetic data. Furthermore, we discuss the issue of how many sampled configurations are necessary to improve over the leading-order ML patterns as a function of the amplitude of the pattern components and of the system size.

The plan of the paper is as follows. In Section II we define the generalized Hopfield model, the Bayesian inference framework and list our main results, that is, the expressions of the patterns without and with corrections, the criterion to decide the number of patterns, and the expressions for the error bars on the inferred patterns. Tests on synthetic data are presented in Section III. Section IV is devoted to the applications to real biological data, i.e recordings of the neocortical activity of a behaving rat and consensus multi-sequence alignment of the PDZ protein domain family. Readers interested in applying our results rather than in their derivation need not read the subsequent sections. Derivation of the loglikelihood with the generalized Hopfield model and of the main inference formulae can be found in Section V. In Section VI we study the minimal number B of samples necessary to achieve an accurate inference, and how this number depends on the number of patterns and on their amplitude. Perspectives and conclusions are given in Section VII.

II. DEFINITIONS AND MAIN RESULTS

A. Generalized Hopfield Model

We consider configurations σ = {σ 1 , , σ 2 , . . . , σ N } of N binary variables taking values σ i = ±1, drawn according to the probability

P H [σ|h, {ξ µ }, { ξµ }] = exp -E[σ, h, {ξ µ }, { ξµ }] Z[h, {ξ µ }, { ξµ }] , (2)
where the energy E is given by

E[σ, h, {ξ µ }, { ξµ }] = - N i=1 h i σ i - 1 2N p µ=1 N i=1 ξ µ i σ i 2 + 1 2N p µ=1 N i=1 ξµ i σ i 2 . ( 3
)
The partition function Z in (2) ensures the normalization of P H . The components of h = (h 1 , h 2 , .., h N ) are the local fields acting on the variables. The patterns ξ µ = {ξ µ 1 , ξ µ 2 , . . . , ξ µ N }, with µ = 1, 2, . . . , p, are attractive patterns: they define preferred directions in the configuration space σ, along which the energy E decreases (if the fields are weak enough). The patterns ξµ , with µ = 1, 2, . . . , p, are repulsive patterns: configurations σ aligned along those directions have a larger energy. The pattern components, ξ µ i , ξµ i , and the fields, h i , are realvalued. Our model is a generalized version of the original Hopfield model [15], which has only attractive patterns and corresponds to p = 0. In the following, we will assume that p + p is much smaller than N .

Energy function (3) implicitly defines the coupling J ij between the variables σ i and σ j ,

J ij = 1 N p µ=1 ξ µ i ξ µ j - 1 N p µ=1 ξµ i ξµ j . (4)
Note that any interaction matrix J ij can be written under the form (4), with p and p being, respectively, the number of positive and negative eigenvalues of J. Here, we assume that the total number of patterns, p + p, i.e. the rank of the matrix J is (much) smaller than the system size, N . The data to be analyzed consists of a set of B configurations of the N spins, σ b , b = 1, . . . , B. We assume that those configurations are drawn, independently from each other, from the distribution P H (2). The parameters defining P H , that is, the fields h and the patterns {ξ µ }, { ξµ } are unknown. Our scope is to determine the most likely values for those fields and patterns from the data. In Bayes inference framework the posterior distribution for the fields and the patterns given the data {σ b } is

P [h, {ξ µ }, { ξµ }|{σ b }] = P 0 [h, {ξ µ }, { ξµ }] P 1 [{σ b }] (5)
× B b=1 P H [σ b |h, {ξ µ }, { ξµ }] ,
where P 0 encodes some a priori information over the parameters to be inferred and P 1 is a normalization.

It is important to realize that many transformations affecting the patterns can actually leave the coupling matrix J (4) and the distribution P H unchanged. A simple example is given by an orthogonal transformation O over the attractive patterns :

ξ µ i → ξµ i = ν O µν ξ ν i .
This invariance entails that the the problem of inferring the patterns is not statistically consistent: even with an infinite number of sampled data no inference procedure can distinguish between a Hopfield model with patterns {ξ µ } and another one with patterns { ξµ }. However, the inference of the couplings is statistically consistent: two distinct matrices J define two distinct distributions over the data.

In the presence of repulsive patterns the complete invariance group is the indefinite orthogonal group O(p, p), which has 1 2 (p + p)(p + p -1) generators. To select one particular set of most likely patterns, we explicitly break the invariance through P 0 . A convenient choice we use throughout this paper is to impose that the weighted dot products of the pairs of attractive and/or repulsive patterns vanish:

i ξ µ i ξ ν i (1 -m 2 i ) = 0 1 2 p(p -1) constraints , i ξ µ i ξν i (1 -m 2 i ) = 0 pp constraints , (6)
i ξµ i ξν i (1 -m 2 i ) = 0 1 2 p(p -1) constraints .
In the following we will use the vocable Maximum Likelihood inference to refer to the case where the prior P 0 is used to break the invariance only. P 0 may also be chosen to impose specific constraints on the pattern amplitude, see Section II E devoted to regularization.

B. Maximum Likelihood Inference: lowest order

Due to the absence of three-or higher order-body interactions in E (3), P depends on the data {σ b } only through the N magnetizations, m i , and the 1 2 N (N -1) two-spin covariances, c ij , of the sampled data:

m i = 1 B b σ b i , c ij = 1 B b σ b i σ b j . (7)
We consider the correlation matrix Γ (1), and call λ 1 ≥ . . . ≥ λ k ≥ λ k+1 ≥ . . . ≥ λ N its eigenvalues. v k denotes the eigenvector attached to λ k and normalized to unity. We also introduce another notation to label the same eigenvalues and eigenvectors in the reverse order: λk ≡ λ N +1-k and vk = v N +1-k , e.g. λ1 is the smallest eigenvalue of Γ; the motivation for doing so will be transparent below. Note that Γ is, by construction, a semi-definite positive matrix: all its eigenvalues are positive. In addition, the sum of the eigenvalues is equal to N since Γ ii = 1, ∀i. Hence the largest and smallest eigenvalues are guaranteed to be, respectively, larger and smaller than unity.

In the following Greek indices, i.e. µ, ν, ρ, correspond to integers comprised between 1 and p or p, while roman letters, i.e. i, j, k denote integers ranging from 1 to N .

Finding the patterns and fields maximizing P (5) is a very hard computational task. We introduce an approximation scheme for those parameters

ξ µ i = (ξ 0 ) µ i + (ξ 1 ) µ i + . . . , ξµ i = ( ξ0 ) µ i + ( ξ1 ) µ i + . . . , h i = (h 0 ) i + (h 1 ) i + . . . . (8)
The derivation of this systematic approximation scheme and the discussion of how smaller the contributions get with the order of the approximation can be found in Section V A. To the lowest order the patterns are given by

(ξ 0 ) µ i = N 1 - 1 λ µ v µ i 1 -m 2 i (1 ≤ µ ≤ p) (9) ( ξ0 ) µ i = N 1 λµ -1 vµ i 1 -m 2 i (1 ≤ µ ≤ p) .
The above expressions require that λ µ > 1 for an attractive pattern and λµ < 1 for a repulsive pattern. Once the patterns are computed the interactions, (J 0 ) ij , can be calculated from ( 4),

(J 0 ) ij = 1 (1 -m 2 i )(1 -m 2 j ) p µ=1 1 - 1 λ µ v µ i v µ j - p µ=1 1 λµ -1 vµ i vµ j . ( 10
)
The values of the local fields are then obtained from

(h 0 ) i = tanh -1 m i - j (J 0 ) ij m j , (11)
which has a straightforward mean-field interpretation.

The above results are reminiscent of PCA, but differ in several significative aspects. First, the patterns do not coincide with the eigenvectors due to the presence of m i -dependent terms. Secondly, the presence of the λ µ -dependent factor in (9) discounts the patterns corresponding to eigenvalues close to unity. This effect is easy to understand in the limit case of independent spins and perfect sampling (B → ∞): Γ is the identity matrix, which gives λ µ = 1, ∀µ, and the patterns rightly vanish. Thirdly, and most importantly, not only the largest but also the smallest eigenmodes must be taken into account to calculate the interactions.

The couplings J 0 (10) calculated from the lowest-order approximation for the patterns are closely related to the mean-field (MF) interactions [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF],

J MF ij = - (Γ -1 ) ij (1 -m 2 i )(1 -m 2 j ) , (12)
where Γ -1 denotes the inverse matrix of Γ (1). However, while all the eigenmodes of Γ are taken into account in the MF interactions (12), our lowest-order interactions [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF] include contributions from the p largest and the p smallest eigenmodes only. As the values of p, p can be chosen depending on the number of available data, the generalized Hopfield interactions (10) is a priori less sensitive to overfitting. In particular, it is possible to avoid considering the bulk part of the spectrum of Γ, which is essentially due to undersampling ( [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF] and Section VI B 2).

C. Sampling error bars on the patterns

The posterior distribution P can locally be approximated with a Gaussian distribution centered in the most likely values for the patterns, {(ξ 0 ) µ }, {( ξ0 ) µ }, and the fields, h 0 . We obtain the covariance matrix of the fluctuations of the patterns around their most likely values,

∆ξ µ i ∆ξ ν j = N M ξξ µν ij B (1 -m 2 i )(1 -m 2 j ) . ( 13
)
and identical expressions for ∆ξ µ i ∆ ξν j and ∆ ξµ i ∆ ξν

M ξξ µν ij = δ µν   N -p k=p+1 v k i v k j |λ k -λµ | + p ρ=1 |λ µ -1|λ ρ v ρ i v ρ j G 1 (λ ρ , λ µ ) + p ρ=1 |λ µ -1| λρ vρ i vρ j G 1 ( λρ , λ µ ) + G 2 (λ µ , λ ν ) G 1 (λ µ , λ ν ) v µ j v ν i , M ξ ξ µν ij = G 2 (λ µ , λν ) G 1 (λ µ , λν ) v µ j vν i , (14)
and M ξ ξ µν ij is obtained from M ξξ µν ij upon substitution of λ µ , λ ν , v µ i , v ν i with, respectively, λµ , λν , vµ i , vν i . Functions G 1 and G 2 are defined through

G 1 (x, y) = (x |y -1| + y |x -1|) 2 , G 2 (x, y) = x y |x -1| |y -1| . ( 15
)
The covariance matrix of the fluctuations of the fields is given in Section V D. Error bars on the couplings (4) can be calculated from the ones on the patterns. Formula (13) tells us how significative are the inferred values of the patterns in the presence of finite sampling. For instance, if the error bar (∆ξ µ i ) 2 1/2 is larger than, or comparable with the pattern component (ξ 0 ) µ i calculated from (9) then this component is statistically compatible with zero. According to formula [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF] we expect error bars of the order of 1 √ α over the pattern components, where α = B N . [START_REF] Amit | Modelling Brain Function: the World of Attractor Neural Networks[END_REF], expressing the rescaled pattern ξ ′ as a linear combination of the eigenvector v and of the orthogonal fluctuations β. The most likely rescaled pattern, (ξ ′ ) 0 , corresponds to a = 1 -1 λ , β = 0.The dashed arc has radius 1 -1 λ . The subscript µ has been dropped to lighten notations.

λ θ β ξ ' v (ξ ) = 0 ' a v 1-1 FIG. 1: Geometrical representation of identity

D. Optimal numbers of attractive and repulsive patterns

We now determine the numbers of patterns, p and p, based on a simple geometric criterion; the reader is referred to Section V E for detailed calculations. To each attractive pattern ξ µ we associate the rescaled pattern (ξ µ ) ′ , whose components are (

ξ µ i ) ′ = ξ µ i 1 -m 2 i / √ N . We write (ξ µ ) ′ = √ a µ v µ + β µ , (16)
where a µ is a positive coefficient, and β µ is a vector orthogonal to all rescaled patterns by virtue of (6) (Fig. 1).

Our lowest order formula (9) for the Maximum Likelihood estimators gives a µ = 1 -1 λ µ and β µ = 0, see Fig. 1. This result is, to some extent, misleading. While the most likely value for the vector β µ is indeed zero, its norm is almost surely not vanishing! The statement may appear paradoxical but is well-known to hold for stochastic variables: while the average or typical value of the location of an isotropic random walk vanishes, its average squared displacement does not. Here, β µ represents the stochastic difference between the pattern to be inferred and the direction of one of the largest eigenvectors of Γ. We expect the squared norm (β µ ) 2 to have a non-zero value in the N, B → ∞ limit at fixed ratio α = B N > 0. Its average value can be straightforwardly computed from formula [START_REF] Aspremont | [END_REF],

(β µ ) 2 = 1 B i M ξξ µµ ii = 1 B N -p k=p+1 1 λ µ -λ k , ( 17
)
where µ is the index of the pattern. We define the angle θ µ between the eigenvector v µ and the rescaled pattern

(ξ µ ) ′ through θ µ = sin -1 (β µ ) 2 1 -1 λ µ , (18)
see Fig. 1. Small values of θ µ correspond to reliable patterns, while large θ µ indicate that the Maximum Likelihood estimator of the µ th pattern is plagued by noise. The value of p such that θ p is, say, about π 4 is our estimate for the number of attractive patterns.

The above approach can be easily repeated in the case of repulsive patterns. We obtain, with obvious notations,

( βµ ) 2 = 1 B i M ξ ξ µµ ii = 1 B N -p k=p+1 1 λ k -λµ , (19)
and

θµ = sin -1 ( βµ ) 2 1 λµ -1 . ( 20
)
The value of p such that θp is, say, about π 4 is our estimate for the number of repulsive patterns.

E. Regularization

So far we have considered that the prior probability P 0 over the patterns was uniform, and was used to break the invariance through the conditions [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF]. The prior probability can be used to constrain the amplitude of the patterns. For instance, we can introduce a Gaussian prior on the patterns,

P 0 ∝ exp - γ 2 N i=1 (1 -m 2 i ) p µ=1 (ξ µ i ) 2 + p µ=1 ( ξµ i ) 2 , (21)
which penalizes large pattern components [START_REF] Cocco | [END_REF]. The presence of the (1 -m 2 i ) factor entails that the effective strength of the regularization term, γ(1 -m 2 i ), depends on the site magnetization. Regularization is particularly useful in the case of severe undersampling. With regularization (21) the lowest order expression for the pattern is still given by ( 9), after carrying out the following transformation on the eigenvalues,

λ µ → λ µ -γ , (µ = 1, . . . , p) , λ k → λ k , (k = p + 1, . . . , N -p) , λµ → λµ + γ , (µ = 1, . . . , p) . ( 22
)
The values of p and p must be such that the transformed λ p and λp are, respectively, larger and smaller than unity. Regularization (21) ensures that the couplings do not blow up, even in the presence of zero eigenvalues in Γ.

Applications will be presented in Sections III and IV. The value of the regularization strength γ can be chosen based on a Bayesian criterion [START_REF] Mackay | [END_REF].

F. Maximum likelihood inference: first corrections

We now give the expression for the first-order correction to the attractive patterns,

(ξ 1 ) µ i = N 1 -m 2 i N k=1 A kµ B kµ v k i , (23)
where

A kµ = C k C µ +   p ρ=1 + N ρ=N +1- p  (λ ρ -1) × i v k i v µ i (v ρ i ) 2 + 2 m i C ρ v ρ i 1 -m 2 i ( 24
)
and

B kµ =        1 2 λ µ λ µ -1 if k ≤ p , √ λ µ (λ µ -1) λ µ -λ k if k ≥ p + 1 . (25)
and

C k = i m i v k i 1 -m 2 i   p ρ=1 + N ρ=N +1- p  (λ ρ -1) (v ρ i ) 2 .
(26) Similarly, the first corrections to the repulsive patterns are ( ξ1

) µ i = N 1 -m 2 i N k=1 Âkµ Bkµ v k i . (27)
The definition of Âkµ is identical to (24), with C µ and v µ i replaced with, respectively, C N +1-µ and vµ i . Finally,

Bkµ =        1 2 λµ 1-λµ if k ≥ N -p + 1 , √ λµ (1-λµ ) λµ -λ k if k ≤ N -p . (28)
The first order corrections to the fields h i can be found in Section V F. It is interesting to note that the corrections to the pattern ξ µ involve non-linear interactions between the eigenmodes of Γ. Formula (24) for A kµ shows that the modes µ and k interact through a multi-body overlap with mode ρ (provided λ ρ = 1). In addition, A kµ does not a priori vanish for k ≥ p + 1: corrections to the patterns have non-zero projections over the 'noisy' modes of Γ. In other words, valuable information over the true values of the patterns can be extracted from the eigenmodes of Γ associated to bulk eigenvalues.

G. Quality of the inference vs. size of the data set

The accuracy ǫ on the inferred pattern is limited both by the sampling error resulting from the finite number of data and the intrinsic error due to the expansion (8). According to Section II C, the sampling error on the pattern component is expected to decrease as ∼ N B . The intrinsic error depends on the order of the expansion, on the size N and on the amplitude of the patterns.

No inference is possible unless the ratio α = B N exceeds a critical value, referred to as α c in the following (Section VI A 2). This phenomenon is similar to the retarded learning phase transition discovered in the context of unsupervised learning [START_REF] Engel | Statistical Mechanics of Learning[END_REF].

Assume that the pattern components ξ i are of the order of one (compared to N ), that is, that the couplings are almost all non zero and of the order of 1 N . Then, the intrinsic error is of the order of 1 N with the lowest order formula (9), and of the order of 1 N 2 when corrections ( 23) are taken into account; for a more precise statement see Section V A and formula (49). The corresponding values of B at which saturation takes place are, respectively, of the order of N 3 and N 5 . The behaviour of the relative error between the true and inferred patterns, ǫ (32), is summarized in Fig. 2. In general we expect that B ∼ N 1+2a samples at least are required to have a more accurate inference with a th -order patterns than with (a -1) th -order patterns. Furthermore there is no need to sample more than N 3+2a configurations when using the a th -order expression for the patterns.

If the system has O(N ) non vanishing couplings J ij of the order of J, then patterns have few large components, of the order of √ J. In this case the intrisic error over the patterns will be of the order of J with the lowest order inference formulae, and of the order of J 2 with the first corrections. The numbers of sampled configurations, B, required to reach those minimal errors will be, respectively, of the order of N J 2 and N J 4 .

III. TESTS ON SYNTHETIC DATA

In this Section we test the formulae of Section II for the patterns and fields against synthetic data generated from various Ising models with known interactions. We consider four models:

• Model A is a Hopfield model with N = 100 spins, p (= 1 or 3) attractive patterns and no repulsive pattern (p = 0). The components of the patterns are Gaussian random variables with zero mean and standard deviation ξ, specified later. The local fields h i are set equal to zero.

• Model B: Model B consists of N spins, grouped into four blocks of N 4 spins each. The p = 3 patterns have uniform components over the blocks:

ξ 1 = 2 √ 3 5 (0, 1, 1, 1), ξ 2 = 2 5 ( √ 3, 1, -2, 1), ξ 3 = 2 5 ( √ 3, -2, 1, 1
). The fields are set to zero. Those choices ensure that the pattern are orthogonal to each other, and have a weak intensity: on average, |ξ| 2 = 9 25 < 1. • Model C is a very simple Ising model where all fields and couplings vanish, except coupling J 12 ≡ J between the first two spins.

• Model D is an Ising model with N = 50 spins, on an Erdos-Renyi random graph with average connectivity (number of neighbors for each spin) equal to d = 5 and coupling values J distributed uniformly between -1 and 1. Model D is an instance of the Viana-Bray model [21]. In the thermodynamic limit N → ∞ this model is in the spin glass phase since d tanh 2 (J) J > 1 [21].

For each one of the models above, the magnetizations and pairwise correlations can be estimated through the sampling of B configurations at equilibrium using Monte Carlo simulations. This allows us to estimate the consequence of sampling noise on the inference quality by varying the value of B. Furthermore, for models B and C, it is possible to obtained the exact Gibbs values for m i and c ij (corresponding to a perfect sampling, B = ∞) [40]. This allows us to study the systematic error resulting from formulae (9,23,27), irrespectively of the sampling noise.

Model A is used to test the lower order formula for the patterns, and how the quality of inference depends on the amplitude of the patterns. Models C and D are highly diluted networks with strong J = O(1) interactions, while models A and B correspond to dense networks with weak J = O( 1 N ) couplings. Models C and D are, therefore, harder benchmarks for the generalized Hopfield model. In addition, the couplings implicitly define, through (4), both attractive and repulsive patterns. The standard deviation of the pattern components is ξ = .7. Left: comparison of the true and inferred couplings for each pair (i, j). Right: comparison of the true and inferred components ξi of the pattern, with the error bars calculated from [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]. The dashed lines have slope unity. Inference is done with p = 1 attractive pattern and no repulsive pattern.

Those models can thus be used to determine how much repulsive patterns are required for an accurate inference of general Ising models.

A. Dominant order formula for the patterns

We start with Model A with p = 1 pattern. In this case, no ambiguity over the inferred pattern is possible since the energy E is not invariant under continuous transformations, see Section II A. We may therefore directly compare the true and the inferred patterns. Figures 3 and 4 show the accuracy of the lowest order formula for the patterns, eqn (9). If the pattern components are weak, each sampled configuration σ is weakly aligned along the pattern ξ. If the number B of sampled configurations is small, the largest eigenvector of Γ is uncorrelated with the pattern direction (Fig. 3). When the size of the data set is sufficiently large, i.e. B > α c N (Section VI A 2), formula (9) captures the right direction of the pattern, and the inferred couplings are representative of the true interactions. Conversely, if the amplitudes of the components of the pattern ξ are strong enough, each sampled configuration σ is likely to be aligned along the pattern. A small number B (compared to N ) of those configurations suffice to determine the pattern (Fig. 4). In the latter case, we see that the largest components ξ i are systematically underestimated. A systematic study of how large B should be for the inference to be reliable can be found in Section VI. We now use model B to generate the data. As model B includes more than one pattern, the inferred patterns cannot be compared to the true one easily due to the invariance of Section II A. We therefore compare in Fig. 5 the true couplings and the interactions found using (9) for three sizes, N = 52, 100 and 200. The size N sets also the amplitude of the couplings, which decreases as 1 N from (4). As the patterns are uniform among each one of the four blocks there are ten possibles values for the couplings J ij , depending on the labels a and b of the blocks to which i and j belong, with 1 ≤ a ≤ b ≤ 4. For N = 100 spins, the relative errors range between 3 and 5.5%. When the number of spins is doubled (respectively, halved) the relative errors are about twice smaller (respectively, larger). This result confirms that formula (9) is exact in the infinite N limit only, and that corrections of the order of O( 1N ) are expected for finite system sizes (Inset of Fig. 5). This scaling was expected from Section II G.

We now consider model C. For perfect sampling (B = ∞) the correlation matrix ( 1) is

Γ =      1 tanh J 0 . . . 0 tanh J 1 0 . . . 0 0 0 1 . . . 0 0 . . . 0 1 0 0 . . . 0 0 1      . ( 29
)
The top eigenvalue, λ 1 = 1+tanh J > 1, and the smallest eigenvalue, λ1 = λ N = 1 -tanh J < 1, are attached to the eigenvectors

v 1 = 1 √ 2       1 1 0 . . . 0       , v1 = 1 √ 2       1 -1 0 . . . 0       . (30)
The remaining N -2 eigenvalues are equal to 1. Using formula [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF] for the lowest order coupling, J 0 , we find that those eigenmodes do not contribute and that the interaction can take three values, depending on the choices for p and p:

(J 0 ) p=1, p=0 = tanh J 2 (1 + tanh J) ≃ J 2 - J 2 2 + J 3 3 + . . . , (J 0 ) p=0, p=1 = tanh J 2 (1 -tanh J) ≃ J 2 + J 2 2 + J 3 3 + . . . , (J 0 ) p=1, p=1 = tanh J 1 -tanh 2 J ≃ J + 2 J 3 3 + . . . . (31)
Those expressions are plotted in Fig. 6. The coupling (J 0 ) 1,0 (dashed line), corresponding to the standard Hopfield model, saturates at the value 1 4 and does not diverge with J. Even the small J behavior, (J 0 ) 1,0 ≃ J 2 , is erroneous. Adding the repulsive pattern leads to a visible improvement, as fluctuations of the spin configurations along the eigenvector v1 (one spin up and the other down) are penalized. The inferred coupling, (J 0 ) 1,1 (bold line), is now correct for small J, (J 0 ) 1,1 ≃ J, and diverges for large values of J. 7: Inferred vs. true couplings for Model D, with B = 4500 sampled configurations. Left: Hopfield model with p = 9 (corresponding to the optimal number of patterns selected by the geometrical criterion); no repulsive pattern is considered (p = 0). Right: Generalized Hopfield model with (p, p) = (9, 35) (optimal numbers).

We now turn to Model D. An illustration of formula [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF] for the error bars is shown in Fig. 3, where we compare the components of the true pattern used to generate data in Model A with the inferred one, (ξ 0 ) i , and the error bar, (∆ξ i ) 2 . For small α = B N the inferred pattern components are uncorrelated with the true pattern and compatible with zero within the error bars. For larger values of α, the discrepancy between the inferred and the true components are stochastic quantities of the order of the calculated error bars.

We report in Fig. 8 the tests of the criterion for determining p and p on artificially generated data from an extension of model A with p = 3 patterns. For very poor sampling (Fig. 8, top) the angle θ 1 is close to π 4 : even the first pattern cannot be inferred correctly. This prediction is confirmed by the very poor comparison of the true interactions and the inferred couplings calculated from the first inferred pattern. For moderately accurate sampling (Fig. 8, middle) the strongest pattern can be inferred; the accuracy on the inferred couplings worsens when the second pattern is added. Excellent sampling allows for a good inference of the structure of the underlying model: the angle θ µ is small for µ = 1, 2, 3 (Fig. 8, bottom), and larger than π 4 for µ ≥ 4 (not shown). Not surprisingly large couplings are systematically affected by errors. Those errors can be corrected by taking into account O( ξ √ N ) corrections to the patterns if the number of data, B , is large enough (Section VI).

Figure 9 compares the inferred and true couplings for B = 4500 sampled configurations of Model D. The optimal number of patterns given by the geometrical criterion is (p = 9, p = 35), see Fig. 7. Hence most of the components of Γ are retained and the interactions inferred with the generalized Hopfield model do not differ much from the MF couplings.

C. Corrections to the patterns

Formula (23) for the corrections to the patterns was tested on model B in the case of perfect sampling. Results are reported in Fig. 10 and show that the errors in the inferred couplings are much smaller than in Fig. 5. Inset of Fig. 10 shows that the relative errors are of the order of 1 N 2 only. This scaling was expected from Section II G. Pushing our expansion of ξ to the next order in powers of 1 N could in principle give explicit expressions for those corrections. We have also tested our higher order formula when the fields h i are non-zero. For instance we have considered the same Hopfield model with p = 3 patterns as above, and with block pseudo-magnetizations t = 1 15 (2 √ 3, 2, 2, -4). Hence, t was orthogonal to the patterns, and the field components were simply given by h i = tanh -1 t i , according to [START_REF] Weigt | [END_REF] 16) and [START_REF] Engel | Statistical Mechanics of Learning[END_REF]. For each value of B one data set was generated from Model A with p = 3 patterns, and standard deviations ξ 1 = .95, ξ 2 = .83, and ξ 3 = .77. aged over the four blocks a) was ∆ta ta ≃ .0301 with the large-N formula (9) and ∆ta ta ≃ 0.0029 with the finite-N formulae (23) and (78).

Corrections to the PCA were also tested when data are corrupted by sampling noise. We compare in Fig. 11 the components of the pattern of Model A found with the lowest order approximation (9) and with our first order formulae (23) (case of strong pattern). A clear improvement in the quality of the inference is observed, even when the sampling noise is strong. Our second example is Model B. We show in Fig. 12 the relative errors

ǫ J = 2 N (N -1) i<j ∆J ij J ij ( 32
)
between the true and the inferred couplings, with formulas ( 9) and ( 23), as a function of the number of sampled configurations, B, and for N = 52 spins. As B increases the relative error with the lowest order patterns (PCA) first decreases as B -1/2 , then saturates to the value ≃ .0794, as expected from Fig. 5. The relative error with the correction to the patterns also decreases as B -1/2 , and is expected to saturate to the lower value ≃ .00374 (Fig. 10). We remark that the gain in accuracy over the inferred couplings resulting from the corrections (23) to the patterns is obtained only when B is very large. B ∼ N 3 configurations at least should be sampled to obtain an improvement over the lowest order formula (9). This scaling holds when the couplings are weak, and decrease as 1 N . If the interaction network is diluted and carries couplings J = O(1), we expect that B ∼ N/J 2 configurations have to be sampled to make the first-corrections to the patterns effective.

We have applied our formula (23) to calculate the first correction to the couplings (31) for Models C and D. As for Model C, we find that the correction to the coupling (J 0 ) 1,1 vanishes; this result is due to the fact that (J 0 ) 1,1 is already correct to the second order in J, and that higher order corrections would be needed. The corrections to the coupling (J 0 ) 1,0 are equal to

(J 1 ) 1,0 = tanh J 2 √ 2 + tanh J(1 + tanh J) 16 = 1 16 + 1 2 √ 2 J + J 2 16 + . . . ... . ( 33
)
The resulting coupling, (J 0 + J 1 ) 10 , is plotted as a function of J in Fig. 6, and qualitatively improves over the lowest order result (31). In particular, for small J, the inferred coupling is now (J 0 + J 1 ) 10 ≃ .916 J -.438 J 2 , which is definitely closer to J than (31). In the case of Model D, the first-order corrections improve only slightly the estimates for the large couplings.

IV. APPLICATION TO BIOLOGICAL DATA

In this Section we show how the inference approach can be applied to real biological data, and compared to other Boltzmann Machine learning procedures.

A. Cortical activity of the rat

We have first analyzed data coming from the recording of 37 neurons in the prefrontal cortex of rats. The experiment, done by A. Peyrache, F. Battaglia and their collaborators, consists in recording the neural activity during a task and during the Slow Wave sleep preceding and following the learning of the task [22]. PCA allowed Peyrache et al. to identify patterns in the activity, which are generated when the rat learns a task and are replayed during the sleep [22].

We have analyzed with the generalized Hopfield model the data corresponding to a 20 minute-long recording of the activity of a rat during the task (data shown in Fig. 1 of [22]). The raster plot was binned with a 10 msec window to obtain binary configurations of the neurons (active or silent in the time-bin). We have then calcu- lated the average frequencies, m i , and the pairwise correlations, c ij . We calculate the couplings with p attractive and p repulsive patterns according to ( 9) and [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF]. The numbers p and p are calculated according to the geometrical criteria [START_REF] Engel | Statistical Mechanics of Learning[END_REF] and [START_REF] Mackay | [END_REF]. Hereafter, we compare the couplings obtained this way to the ones found with the adaptive cluster expansion (ACE) of [START_REF] Cocco | [END_REF], which is not based on the expansion of the loglikelihood used in the present work.

In Fig. 13 (top) we compare the Hopfield (p = 0) couplings with p = 4, 8, 17 selected patterns to the ACE couplings. The agreement is quite good for p * = 8. In [22] p = 6 patterns were kept in the PCA; this value is close to the optimal value, p = 8, we find using the geometrical criterion. Addition of repulsive patterns (bottom of Fig. 13) slightly improves the similarity with the ACE couplings. We find, indeed, that the couplings J ij are rather weak, and that repulsive patterns do not play an important role. Calculating the couplings with all eigenmodes (p = 17, p = 20) is equivalent ot the mean-field (MF) approximation. A clear discrepancy between the Hopfield and the ACE couplings is found for the largest (in absolute value) interactions. We have checked that this discrepancy is not reduced when the first order corrections to the patterns are included, presumably because the number of data is not sufficient. Couplings are not significatively changed in the presence of the regularization (21) for sensible values of γ.

B. Protein-domain families

We have next analyzed the alignement of a family of 240 sequences of PDZ, a commonly encountered domain binding the C-terminus of proteins, with 92 amino-acids [24]. R. Ranganathan and collaborators have elaborated an approach, called Statistical Coupling Analysis(SCA), to extract interactions between residues by using evolutionary data for the protein, i.e. by sampling the singlesite and pairwise frequencies from multi-sequence alignments of the family [23]. Briefly speaking, SCA consists in doing a PCA analysis of a weighted correlation matrix, D i Γ ij D j , where the weight D i on site i is small for poorly conserved residues [24].

We have taken the binary data representation of the 240 PDZ sequences in the alignement given in [25] (Supplementary Material). This consensus approximation amounts to replace the amino-acid on each site (20 possible types) with a binary variable σ b i , equal to +1 if the amino-acid i in the b th sequence is the most common amino-acid at that position in the alignment, to -1 otherwise. The consensus representation does not allow to keep track of all the information contained in the alignment but is indicative of the conservation pattern in the family.

The inferred couplings, denoted by J 92 , are shown in Fig. 14. As in the case of Model D in Section III we find that proteomic data are better accounted with the generalized Hopfield model than with the standard Hopfield model: repulsive patterns seems necessary to recover the couplings found with the ACE method. The couplings found with attractive patterns only are not correlated with the ACE couplings (top of Fig. 14), while the agree- ment is quite good when taking into account attractive and repulsive patterns; the optimal numbers of patterns are p = 4 and p = 10.

-2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 p=3 p=16 p=30 p=3 p=3 p=4 p=10 p=16 p=28 -2 -1 0 1 2 -2 -1 0 1 2 J 44 (Hopfield) -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 J 44 (Adaptive Cluster Expansion) -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0
We have also calculated the couplings when discarding all but the most weighted sites. More precisely, we have recalculated the distribution of the weights D i as in [24,25], and found a bimodal distribution, which suggests a natural cut-off between large and small weights. We have redone the previous inference when keeping only the 44 residues (out of 92) with the largest weights, corresponding to the red sites in Fig. C of [24]. The resulting interactions, denoted by J 44 , are shown in Fig. 15. Again we compare the couplings found with the Hopfield model and with the ACE. The agreement is not good with attractive patterns only (as done in usual PCA), and is very good when repulsive patterns are included.

An interesting question is whether the couplings obtained between the 44 most conserved residues are strongly affected by the presence or the absence of the remaining 48 residues in the inference. The interactions in the 44-site model are effective and a priori differ from their values in the 92-site model, in that they account for chains of interactions going through the remaining 48 sites. Nevertheless, we find that the couplings calculated with all 92 residues and the couplings obtained from the subset of 44 sites with large weights are similar, see Fig. 16. This result suggests that the 48 residues removed from our second analysis are not strongly interacting with the 44 retained sites.

V. EXPANSION OF THE CROSS ENTROPY AND MAXIMUM LIKELIHOOD INFERENCE

This Section is intended to provide the derivations of the results announced in Section II. Maximizing the posterior probability [START_REF] Ackley | [END_REF] with respect to the patterns and the fields is equivalent to minimizing the cross entropy of the Hopfield model given the data,

Φ[h, {ξ µ }, {σ b }] = log Z[h, {ξ µ }] + U [h, {ξ µ }, {σ b }] , ( 34
)
where Z is the partition function appearing in [START_REF] Seth | [END_REF],

Z[h, {ξ µ }] = σ exp -E[σ, h, {ξ µ }] , (35)
and U is the average value of the energy E (3) over the sampled configurations:

U [h, {ξ µ }, {σ b }] = - N i=1 h i m i - 1 2 i,j J ij c ij , (36)
where the couplings J ij are calculated from the patterns according to (4). The calculation of the partition function, which is defined as a sum over 2 N configurations, cannot generally be done in a reasonable time for large sizes N . In the next section we show how the use of statistical mechanics techniques allows one to obtain a systematic expansion of Z, and, thus, of the cross entropy

Φ = Φ 0 + Φ 1 + . . . , (37)
in powers on ξi √ N and ξi √ N .

A. Expansion of the free energy of the Hopfield model in powers of

ξ i √ N , ξi √ N
To lighten notations calculations are presented for the case of attractive patterns only. We explain at the end of the Section how formulae are modified in the presence of repulsive patterns.

For technical reasons to be made clear below it results convenient to make the change of variables h → t described by

h i = tanh -1 t i - 1 N µ j ξ µ i ξ µ j t j , (38)
where the t i , hereafter called pseudo-magnetizations, are real-valued numbers comprised between -1 and 1. Hereafter, we will infer the most likely values for t, and will recover the fields h through [START_REF] Weigt | [END_REF]. The change h → t amounts to consider the energy function

E = - N i=1 σ i tanh -1 t i - 1 2N p µ=1 N i=1 ξ µ i σ i -t i 2 , (39)
instead of the original expression for E (3) (with p = 0). Obviously, when the identities (38) are fulfilled, both energies are equal (up to a σ-independent additive term) and define the same likelihood function [START_REF] Seth | [END_REF].

We unravel the squared terms in the partition function (35) through a set of p auxiliary Gaussian variables x = (x 1 , . . . , x p ), and carry out the summation over the spin configurations. We obtain

Z = µ dx µ √ 2π exp - 1 2 µ (x µ ) 2 - i,µ x µ ξ µ i t i √ N + i log 2 cosh tanh -1 t i + µ x µ ξ µ i √ N . ( 40
)
If N is large enough the dominant contribution to the integral will come from x * , the value of x maximizing the argument of the exponential above. We obtain the following saddle point equation for x,

(x µ ) * = 1 √ N i ξ µ i (T i -t i ) , (41)
where

T i ≡ tanh tanh -1 t i + µ (x µ ) * ξ µ i √ N ( 42
)
We then write x µ = (x µ ) * +y µ and expand the hyperbolic cosine function in powers of y µ . The change of variable [START_REF] Weigt | [END_REF] is such that the linear term in y µ in the expansion of the hyperbolic cosine function cancels out with the linear term in the exponential,i,µ

y µ ξ µ i t i √ N , indepen-
dently of the value of (x µ ) * . Expanding the hyperbolic cosine up to the second order in y µ we find our lowest order approximation to the partition function,

Z 0 = e F * µ dy µ √ 2π exp - 1 2 µ (y µ ) 2 (43) + 1 2 N i µ,ν ξ µ i ξ ν i y µ y ν 1 -T 2 i = e F * √ det A
where F * is the the argument of the exponential in (40) calculated in x µ * ,

F * = N log 2 + 1 2 i log(1 -T 2 i ) - µ,ij ξ µ i ξ µ j (T i T j -t i t j ) ,
(44) and A is the p × p matrix with entries,

A µν = δ µν - 1 N i ξ µ i ξ ν i 1 -T 2 i . ( 45
)
We then compute the average energy U (36),

U = - i m i tanh -1 t i (46) - 1 2N µ,ij ξ µ i ξ µ j (c ij -m i t j -t i m j + t i t j ) .
Our lowest order approximation for the cross entropy is, according to ( 34), ( 44) and ( 46):

Φ 0 = - N i=1 m i tanh -1 T i + N log 2 + 1 2 i log(1 -T 2 i ) - 1 2N µ,ij ξ µ i c ij -m i m j ξ µ j - 1 2 log det A + 1 2N µ i ξ µ i T i -m i 2 . ( 47
)
The first order contribution to the cross entropy, Φ 1 in (37), is obtained by retaining the fourth order in y µ in the expansion of the hyperbolic cosine function in (40),

Φ 1 = 1 4N 2 i (1 -4T 2 i + 3T 4 i ) µ,ν ξ µ i (A -1 ) µν ξ ν i 2 .
(48) We expect the differences Φ -Φ 0 and Φ -(Φ 0 + Φ 1 ) between, respectively, the true and the lowest order cross entropies and the true and the first order cross entropies to be of the order of, respectively, R 2 and R 3 , where

R = p N ξ 2 (1 -m 2 ) Λ . ( 49
)
Here, ξ 2 is the order of magnitude of the pattern components, which can range from 1 if the patterns are extended over the whole system to ∼ √ N for highly sparse patterns, m is the typical value of the local magnetization, and Λ is the order of magnitude of the eigenvalues of A -1 , which can range from 1 to N . The value of R fixes the instrinsic error ǫ on the inferred patterns discussed in Section II G, ǫ ∼ R for the lowest order approximation and ǫ ∼ R 2 with the first order corrections.

The above calculation can be straightforwardly extended to the case of the generalized Hopfield model by considering the p repulsive patterns as patterns with purely imaginary components, ξ µ = i ξµ , with i 2 = -1. For instance the general lowest order expression for the cross entropy is

Φ 0 = - N i=1 m i tanh -1 T i + N log 2 + 1 2 i log(1 -T 2 i ) - 1 2N ij c ij -m i m j p µ=1 ξ µ i ξ µ j - p µ=1 ξµ i ξµ j + 1 2N p µ=1 i ξ µ i T i -m i 2 - 1 2N p µ=1 i ξµ i T i -m i 2 - 1 2 log det A i  -i ÂT  , (50)
where

T i = tanh tanh -1 t i + p µ=1 (x µ ) * ξ µ i √ N - p µ=1 (x µ ) * ξµ i √ N , (x µ ) * = 1 √ N i ξµ i (T i -t i ) , µν = 1 N i ξ µ i ξν i 1 -T 2 i , µν = δ µν + 1 N i ξµ i ξν i 1 -T 2 i . ( 51
)
The first order correction (48) can be easily written for the case of repulsive patterns, too.

B. Are the physical properties of the system relevant for the inference?

The Hopfield model was first introduced as a model for which a set of p desired ground states ξ µ (or fixed points of the zero temperature Glauber dynamics) could be programmed through an adequate choice of the interactions. Each fixed point has a basin of attraction in the configuration space, corresponding to a phase of the system. The order parameters are the overlaps

q µ = σ P H [σ|h, ξ] 1 N i ξ µ i σ i , (52)
which quantify how much the configurations are on average aligned along each pattern. The amplitudes and directions of the pattern and the field vectors determine if spin configurations tend to be aligned along the field, or along one or more patterns. In the infinite size limit (N → ∞) the overlaps are the roots of p coupled and self-consistent equations,

q µ = lim N →∞ 1 N i ξ µ i tanh h i + ρ q ρ ξ ρ i . (53)
Using [START_REF] Weigt | [END_REF] and the saddle point equation ( 41) it is easy to check that the overlaps

q µ = 1 N i ξ µ i T i ( 54
)
are solutions to the set of equations ( 53). Solutions are in one-to-one correspondance with the saddle points (x µ ) * . The saddle-point solution x * = 0 corresponds to T i = t i . The average interaction term in the energy function [START_REF] Tibshirani | [END_REF] vanishes, meaning that configurations tend to be mainly determined by the fields. Such a behaviour corresponds to the paramagnetic phase. The solution x = 0 is locally stable if the eigenvalues of the matrix A are all positive and, thus, if the patterns are weak enough. Solutions with x * = 0 correspond to stronger patterns and interaction terms in [START_REF] Tibshirani | [END_REF] having non zero values on average: they correspond to magnetized phases.

The cross entropy Φ depends on the solution x * through the variables T i only. Once the T i 's and the patterns ξ µ 's are inferred, it is easy to calculate the value of the fields h i based on equations [START_REF] Weigt | [END_REF], ( 41) and (42). One finds that h i is given by [START_REF] Weigt | [END_REF] where t i is substituted with T i . Hence, the inferred parameters do not explicitely depend on the value of x * . The procedure followed to infer the patterns and the fields is not affected by the physical phase (paramagnetic or magnetized) of the system, though the values of the data m i and c ij obviously depend on those physical properties.

It may accidentally happen that equations (41) have different solutions with equal or almost equal contributions to the partition function Z. The most natural illustration is the case of zero field (t i = 0) and one strong pattern, where two ferromagnetic states with opposite overlaps, (x 1 ) * and -(x 1 ) * , coexist. In this latter case both states give equal contributions to the partition function.

C. Maximum Likelihood inference: lowest order

We first infer the patterns and the pseudomagnetizations from Φ 0 . Minimization of Φ 0 (47) over T immediately shows that, up to O(R) corrections, pseudoand true magnetizations coincide:

(T i ) * = m i . ( 55
)
Without loss of generality we may write the patterns to infer as

ξ 0 µ i = √ N a µ v µ i + √ N β µ i 1 -m 2 i , ξ0 µ i = √ N âµ vµ i + √ N βµ i 1 -m 2 i , (56)
where √ a µ , √ âµ are real-valued coefficients, and v µ and vµ are eigenvectors of Γ. According to identity (55) the conditions [START_REF]Advanced Mean Field Methods: Theory and Practice[END_REF] are fulfilled in the large N limit if the (p+p) vectors β µ and βν are orthogonal to each other, and to all the patterns (ξ 0 ) ν and ( ξ0 ) ν . The matrices A (45) and  (51) are then diagonal, while  vanishes. We rewrite the cross entropy (50) as

Φ 0 = - i σ=±1 1 + σ m i 2 log 1 + σ m i 2 - 1 2 µ λ µ a µ - 1 2 ij,µ β µ i Γ (r) ij β µ j , + 1 2 µ λµ âµ + 1 2 ij,µ βµ i Γ (r) ij βµ j , - 1 2 µ log 1 -a µ - i (β µ i ) 2 - 1 2 µ log 1 + âµ + i ( βµ i ) 2 (57)
where Γ (r) is the restriction of Γ to the (N -p -p)dimensional subspace orthogonal to the p largest and p smallest eigenvectors:

Γ (r) ij = N -p k=p+1 λ k v k i v k j . ( 58
)
Minimizing Φ 0 over the coefficients a µ and the vectors β µ gives the coupled set of equations

λ µ = 1 1 -a µ -b µ , ( 59
) j Γ (r) ij β µ j = β µ i 1 -a µ -b µ , ( 60
)
where b µ = (β µ ) 2 is the squared norm of β µ . If the vector β µ were non zero, it would be an eigenvector of Γ with eigenvalue λ µ according to (60). This cannot be true as the largest eigenvalue of Γ (r) is smaller than λ p . Hence, β µ = b µ = 0. From (59) we obtain

a µ = 1 - 1 λ µ . (61)
We conclude that the maximum likelihood values for the p attractive patterns are given by (9). The minimization of Φ 0 over the coefficients âµ and the vectors βµ can be done along the same lines. We find

âµ = 1 λµ -1 . (62)
and βµ = 0. The maximum likelihood estimators for the p repulsive patterns are given by ( 9) again. Once the patterns are computed the values of the local fields h i are obtained from [START_REF] Cocco | [END_REF].

Notice that v µ i , vµ i are typically of the order of N -1 2 , which entails that the components of the patterns are of the order of unity. Though keeping each ξ i , ξi of the order of unity is a natural scaling in the infinite size limit N → ∞, other scalings are possible. Consider a pair of strongly coupled spins, i.e. such that the correlation Γ ij is sizeably larger than 1 N . According to expression (4) for the coupling J ij induced by the patterns between spins i and j, we expect the pattern components to be of the order of √ N . There is thus no compelling reason to assume that ξi √ N , ξi √ N is vanishingly small for all components i. To end with we compute the decrease in cross entropy when adding a pattern attached to the eigenvalue λ (= λ µ or λµ ). Inserting expressions (61,62) for a µ , âµ in (57) we obtain

∆Φ = - 1 2 λ -1 -log λ , (63)
a quantity which is strictly negative for λ = 1. Not surprisingly, adding more parameters to the model allows for a better fit of the data. We will see in Section V E how the values of p and p can be determined.

D. Error bars on the patterns and fields

When the sample size B is large the posterior distribution P tends to a Gaussian law centered in the most likely values for the patterns, {ξ µ }, { ξµ }, and the pseudomagnetizations, T. For the sake of simplicty we consider below the case of attractive patterns only; repulsive patterns can formally be seen as purely imaginary attractive patterns, see Section V A. Let H denote the Hessian matrix of Φ 0 . We find, to the leading orders,

(H tt ) ij ≡ ∂ 2 Φ 0 ∂T 1 i ∂T 1 j = δ ij 1 -m 2 i -(J 0 ) ij , (H ξξ ) µν ij ≡ ∂ 2 Φ 0 ∂(ξ 0 ) µ i ∂(ξ 0 ) ν j = δ µν N m i m j -c ij + (1 -m 2 i ) × λ µ δ ij + (1 -m 2 j ) ρ λ ρ N (ξ 0 ) ρ i (ξ 0 ) ρ j + λ µ λ ν N 2 (1 -m 2 i )(1 -m 2 j )(ξ 0 ) ν i (ξ 0 ) µ j , (64)
(H tξ ) ν ij ≡ ∂ 2 Φ 0 ∂T i ∂(ξ 0 ) ν j ≃ 0 . (65)
Here, δ denotes the Kronecker function and the expression of the lowest order coupling matrix, J 0 , is given in [START_REF] Cocco | Proc. Nat. Acad. Sci[END_REF]. The sum over ρ runs over all pattern indices. The cross second derivative, H tξ , of the order of |ξ| N , is much smaller than the expected order, |ξ| √ N , and can be neglected.

The covariance matrix of the Gaussian posterior probability P is the inverse matrix of B H. The inverse is properly defined in the subspace of dimension N (p + p + 1) -1 2 (p+ p)(p+ p-1), orthogonal to the modes generating the invariance over the patterns, see Section II A. We write H = DHD, where D is a diagonal matrix with elements:

D i = 1 -m 2
i in the T-sector, and

D µ i = N 1-m 2 i
in the ξ µ -sector. Matrix H has a particularly simple expression in the eigenbasis of the correlation matrix Γ, and can be diagonalized exactly after some simple algebra. We obtain the following expression for the covariance matrix of the fluctuations:

∆T i ∆T j = (1 -m 2 i )(1 -m 2 j ) B M tt ij , (66)
where

M tt ij = δ ij + p ρ=1 (λ ρ -1) v ρ i v ρ j + p ρ=1 ( λρ -1) vρ i vρ j . (67)
The expressions for the fluctuations of the pattern components are reported in [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]. Note that the cross-term ∆T i ∆ξ ν j vanishes at the expected order of √ N B , and is actually of the order of 1 B only. Using formula [START_REF] Weigt | [END_REF] we find that the error over the fields h i is of the order of p √ α , where α = B N .

E. Optimal number of patterns

So far we have assumed that the number of patterns, p, was known. In practice p is often determined based on simple criteria, such as how many eigenvalues 'come out' from the spectrum of the correlation matrix (Section VI B 2). Alternative approaches exist, e.g. Bayesian Information Criterion (BIC) [26]. In the BIC the decrease B∆Φ (63) in cross entropy obtained with a new pattern is added a 'cost' N log B, equal to the number of new parameters times the logarithm of the number of data. As the index µ increases the selected eigenvalue λ µ or λµ gets closer to one; B|∆Φ| (63) decreases in absolute value, and, eventually, is counterbalanced by the cost term N log B. The value of µ for which the two terms balance each other depends on the size of the data set: the higher B, the more significative are the correlations and the more patterns we need to represent the interactions. However BIC is mathematically justified when B is large compared to N , which is not always the case in real data sets.

Hereafter, we propose a different approach based on Bayesian and geometric considerations. Based on the discussion in Section II D we expect the squared norm b µ of the transerve fluctuations β µ to be non vanishing in the B, N → ∞ limits. Let us call a µ the squared projection of the µ th rescaled pattern onto v µ (16). The same quantities, âν and bν , can be defined for repulsive patterns. We define the marginal probability P M of the squared projections a µ , âν and of the squared norms b µ , bν through

P M = µ,i dβ µ i 1 -m 2 i ν,i d βν i 1 -m 2 i µ dΩ µ πiαN/2 × ν d Ων πiαN/2 exp - α 2 µ Ω µ (β µ ) 2 -N b µ × exp - α 2 ν Ων ( βν ) 2 -N bν (68)
× P {T 0 i , √ N a µ v µ i + √ N β µ i 1 -m 2 i , √ N âν vν i + √ N βν i 1 -m 2 i } ,
where P is the posterior probability (5), and the sums over µ and ν run from 1 to, respectively, p and p. After carrying out the integrals over the fluctuations β µ and βν we obtain

P M = 1 Z 1 µ dΩ µ ν d Ων (69) × exp - B 2 µ ∆Φ M (Ω µ ) - B 2 ν ∆ ΦM ( Ων )
where Z 1 is a normalization constant and

∆Φ M (Ω µ ) = λ µ a µ + Ω µ b µ + log (1 -a µ -b µ ) (70) - 1 B log det Ω µ 1 -Γ (r) + O log N N ,
∆ ΦM ( Ων ) = -λν âν + Ων bν + log 1 + âν + bν (71)

- 1 B log det Ων 1 + Γ (r) + O log N N ,
Here 1 denotes the N -dimensional identity matrix. When B is large the integrals in (69) are dominated by the contributions coming from the vicinity of the roots of

∂∆Φ M ∂Ω µ = ∂∆ ΦM ∂ Ων = 0 . (72)
Maximimization of ∆Φ M with respect to the a µ , b µ 's gives equations (59) and

Ω µ = λ µ , (73)
for each µ = 1, . . . , p. We then compute the squared norm b µ from the extremization condition (72) and obtain

b µ = 1 B N -p k=p+1 1 λ µ -λ k , (74)
a µ = 1 - 1 λ µ -b µ . (75)
Repeating the same procedure to maximize ∆ ΦM gives bν

= 1 B N -p k=p+1 1 λ k -λν , âν = 1 λν -1 -bν . (76)
The difference between expressions (61) and (75) for the coefficients a µ must be emphasized. P defined in ( 5) is a probability density over pN pattern components, once the pseudo-magnetizations have been inferred. Maximization of P , or, equivalently, of Φ over this large-dimensional space gives expression (61) for the projection a µ of the pattern ξ µ onto the µ th largest eigenvector of Γ, v µ . Instead of directly maximizing P , we may first integrate out the orthogonal fluctuations to v µ in P , and obtain the marginal probability density P M for 2p parameters only, namely the squared projections on the eigenvectors, a µ , and the squared norms of the orthogonal fluctuations, b µ . Maximizing the marginal probability density P M or, equivalently, minimizing Φ M shows that b µ (75) does not vanish, and that the value of the squared projection a µ (75) is smaller than (61). Figure 1 sketches the geometrical meaning of the coefficient √ a µ and the fluctuations β µ , see [START_REF] Amit | Modelling Brain Function: the World of Attractor Neural Networks[END_REF]. Small values of the angle θ µ are expected for reliable patterns. A similar picture can be drawn for repulsive patterns. We will see how expression (75) for the squared norm b µ naturally arises in the context of random matrix theory.

F. Maximum likelihood inference: first corrections

We now look for the corrections to the lowest order expressions of the patterns and the fields (9,55), encoded in expressions (8) and T i = T 0 i + T 1 i . The first order contribution to the cross entropy, Φ 1 , can be seen as a perturbation to the lowest order cross entropy, Φ 0 , according to [START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF]. Within linear response theory this perturbation will shift the maximum likelihood estimators by

  T 1 {(ξ 1 ) µ } { (ξ 1 ) µ }   = -H -1    ∂Φ 1 ∂T { ∂Φ 1 ∂ξ µ } { ∂Φ 1 ∂ ξµ }    , (77)
where the inverse of the Hessian matrix of Φ 0 , H, was given in Section V D. The calculation of the gradient of Φ 1 does not present any particular difficulty. The resulting corrections to the patterns are given in eqn (23). The expression for the shift in the pseudo-magnetization is

T 1 i = p µ=1 (λ µ -1) C µ v µ i 1 -m 2 i + m i (v µ i ) 2 (78)
+ p µ=1 ( λµ -1) C N +1-µ vµ i 1 -m 2 i + m i (v µ i ) 2 .
where C k is given in (26). Notice that, if the magnetizations m i vanish, so do the dominant and first-order contributions to the pseudo-magnetizations.

VI. RELIABILITY OF THE INFERENCE

An important issue is to determine how many configurations should be sampled in order to ensure that the inference of the patterns is accurate. To do so, we assume that the examples σ b are drawn independently and at random from the equilibrium probability P H (2) of a Hopfield model, with fixed fields h and patterns ξ. We call S[{σ b }] the entropy of the posterior distribution P (5) for the fields h and patterns ξ. In the large N limit, we expect this entropy to be self-averaging, that is, to depend on the set of examples only through their number B. We want to determine how fast S decays with B. To do so it is instructive to first consider the simple case where the local fields are known, and only one pattern has to be inferred. This specific situation is treated in great analytical details in Section VI A. The general (and harder) case where both fields and patterns have to inferred is treated in Section VI B.

A. Case of one unknown pattern and known fields

Throughout this Section, we assume that the local fields vanish, h = 0 and that the number of patterns to be inferred is p = 1. The posterior entropy, [27], the examples σ b have vanishingly small overlap (52) with the pattern. We expect that a large number B (diverging with N ) of examples is necessary to convey reliable information about the pattern. Conversely, few configurations sampled in a ferromagnetic state around a strong pattern (or its opposite) should be sufficient to reconstruct the pattern.

We now make this scenario quantitative in various cases. An important simplication arises when the pattern is restricted to have binary components, ξ = { ξi = ± ξ}, with ξ > 0. Hamiltonian (3) with p = 1 pattern is invariant under the exchange of the spin configuration and the pattern: patterns. In particular, the posterior entropy S is equal to the entropy of the dual Hopfield model at inverse temperature

E[σ, 0, ξ] = E[ξ, 0, σ].
β = ξ2 . ( 80
)
The duality property allows us to exploit the wellunderstood physics of the Hopfield model [27] to simplify the study of our inference problem.

Strong components

In the ferromagnetic regime ( ξ > 1), the dual spin configuration is strongly magnetized along the dual pat-terns. Going back to the inference problem, we find that the overlap between the inferred pattern and a sampled configuration,

q b = {σ b },ξ P [0, ξ|{σ b }] b P H [σ b , ξ] 1 N i ξ i σ 1 i , (81)
may take values +q or -q, where q is the positive root of q = tanh(q ξ2 ). The sign of the overlap q b is random, depending on which one of the two states with opposite magnetizations the configuration σ b in sampled in; it is equal to + or -with equal probabilities 1 2 . These statements hold if the thermodynamical limit, N → ∞, is taken while B is kept fixed. We find that S is equal to the entropy of a single spin at inverse temperature β, interacting with B other spins of magnetization q, S =

B b=0 B b 1 + q 2 b 1 -q 2 B-b S (B -2b)q ξ2 , (82)
where S(u) = log(2 cosh u) -u tanh u. Figure 17A shows that the entropy is almost a pure exponential: log S ≃ -B/B c where the decay constant, B c = 1/ log cosh(q ξ2 ), is finite (compared to N ). In the ferromagnetic regime few sampled configurations are sufficient to determine ξ accurately.

This result also applies to the case of a single ferromagnetic state. If the field h does not strictly vanish and explicitely breaks the reversal symmetry between the two states, all configurations are sampled from the same state, with probability 1 -exp(-O(N )). Remarkably, expression (82) for the entropy still holds. Again we find that B = O(1) configurations are sufficient to infer the pattern. We will discuss in more details the inference in the ferromagnetic regime in Sections VI B 1 and VI B 3.

Weak components

In the paramagnetic phase ( ξ < 1), the overlap (81) between the inferred pattern and an example is typically very small, q ∼ N -1/2 . No inference is possible unless the number of examples, B, scales linearly with N ; we denote α = B/N . In this regime, we expect the entropy to be self-averaging: S[{σ b }] does not depend on the detailed composition of the data set and is a function of the value of the macroscopic parameters, e.g. the ratio α, only. To calculate this function S we use the replica method [START_REF] Amit | Modelling Brain Function: the World of Attractor Neural Networks[END_REF]27]. We report below the results of the replica symmetric calculation; technical details can be found in Appendix. The order parameter is the average overlap r between the inferred and the true patterns,

r = {σ b },ξ P [0, ξ|{σ b }] b P H [σ b , ξ] 1 N i ξ i ξi . ( 83
)
which is solution of the self-consistent equation

r = ∞ -∞ Dz tanh(z √ γ + γ) , (84)
where Dz = dz √ 2π e -z2 /2 is the Gaussian measure, and

γ = αβ 2 r (1 -β)(1 -β + βr) . (85)
The posterior entropy is equal to

S = ∞ -∞ Dz log 2 cosh(z √ γ + γ) - α 2 log(1 -β + βr) - αβ(1 -β -r + 3βr) 2(1 -β)(1 -β + βr) , (86)
and is plotted in Fig. 17B. To check this analytical prediction we have run extensive numerical simulations on small-size systems (N = 10, 20). The numerical procedure follows three steps: 1. evaluate the partition function Z in (2) through an exact enumeration; 2. generate a data set of B = αN configurations {σ b i } according to the Hopfield measure P H by rejection sampling; 3. evaluate P 1 in (5) and S in (79) through exact enumerations. The resulting entropy, averaged over one hundred data sets, is compatible with the analytical prediction and the existence of 1 N finite-size effects. Inset of Fig. 17B shows that the overlap r remains null until α reaches the critical value

α c = 1 ξ2 -1 2 . ( 87
)
Hence, in the range [0; α c ], the posterior probability becomes more concentrated (S decreases), but not around the true pattern ξ. The existence of a lagging phase before any meaningful inference is possible is similar to the 'retarded learning' phenomenon discovered in the field of unsupervised learning, where the variables to be learned are real-valued [28][29][30]. In the present case of binary spins we expect the replica symmetric assumption to break down at large α. The entropy (86) indeed becomes negative when α > α 0 ≃ 42 for the case studied in Fig. 17B. Nevertheless we may conjecture that the entropy decays as S ∼ 1 α when α → ∞. The dual Hopfield model has random couplings J ij , with second moment equal to J 2 ij -J ij 2 = α N . Hence T = 1 √ α sets the temperature scale of the dual model. The low temperature scaling of the entropy of the Sherrington-Kirkpatrick (SK) model suggests that S ∝ T 2 [31]; this scaling is compatible with the small-N results of Fig. 17B. However the dual and SK models are not strictly identical when α → ∞: the coupling matrix J of the dual model is guaranteed to be semidefinite positive, while the entries of J are independent in the SK model. A complete calculation of the entropy valid for any (large) α would require a replica symmetry broken Ansatz for the order parameters [32], and is beyond the scope of this article.

Note that the calculations above can be extended to real patterns; β in (80) is then replaced with ξ 2 , where the average is taken over the pattern components. The entropy is not constrained to be positive as in the binary case. The distinction between the strong-and weakcomponent regimes remains qualitatively unchanged, and so does the value of the critical ratio α c (87), which does not depend on the third and higher moments of ξi .

B. General case of unknown patterns and fields

In this Section, we first interpret the above results. We show that, while B = O(1) configurations can be sufficient in a particular context, B = O(N ) data are generally necessary for the inference to be sucessful. The connection between the results of Section VI A and random matrix theory are emphasized.

Inference from the magnetizations

Consider first the case where a single state exists, i.e. equations (53) admit a single solution {q µ }; the case where states coexist will be discussed in Section VI B 3. For large N , the average value of spin i with the measure P H (2) is

m i = tanh h i + µ q µ ξ µ i . ( 88
)
As the error on the estimate of m i decreases as

∼ 1-m 2 i B with B, O (1)
configurations are sufficient to sample the magnetizations accurately. Few sampled configurations therefore give access to the knowledge of a linear combination of the field vector and pattern vectors with non zero-overlaps q µ . This linear combination is simply T 0 i , and equation (88) coincides with (55).

When the fields h i are known and the model consists of a single strong pattern (p = 1) the pattern components ξ 1 i can be readily calculated from the magnetizations (88) through

ξ 1 i = 1 q tanh -1 m i where q 2 = 1 N j m j tanh -1 m j .
(89) This particular case was encountered at the end of Section VI A 1, when the fields h i are sent to zero after having broken the reversal symmetry of the system to avoid state coexistence. In the generic situation of unknown fields and patterns, knowledge of the magnetizations does not suffice to determine the field and the patterns, and must be supplemented with the information coming from the correlation matrix Γ ij . patterns in the energy (3) the spins are correlated. The entries of the correlation matrix are, for large N [42],

Γ ij = δ ij + 1 N ξ i ξ j (1 -m 2 i )(1 -m 2 j ) 1 -1 N k ξ 2 k (1 -m k ) 2 (90)
where we have considered the case of a single pattern (p = 1, p = 0) to lighten notations. Though the pattern affects each correlation Γ ij by O( 1 N ) only, these small contributions add up to boost the largest eigenvalue from one (in the absence of pattern) to

L = 1 1 -1 N k ξ 2 k (1 -m k ) 2 . (91)
The eigenvector attached to L has components v i ∝ ξ i 1 -m 2 i and ML inference perfectly recovers the pattern.

In the presence of sampling noise (finite B), each correlation (90) is corrupted by a stochastic term of the order of x = 1 √ B . This stochastic term will, in turn, produce an overall contribution of the order of x √ N = 1 √ α to the largest eigenvalue. Intuitively, whether α is large or small compared to L -2 should tell us how hard or easy it is to extract the pattern ξ from Γ. Several studies in the physics [33,34] and in the mathematics [35] literatures have indeed found that an abrupt phase transition takes place at the critical ratio

α c = 1 (L -1) 2 . (92)
It is a simple check that α c coincides with the ratio (87) for the retarded learning transition calculated in Sections VI A 2.

In the strong noise regime (α < α c ) the largest eigenvector v 1 of Γ is uncorrelated with (orthogonal to) the pattern ξ, and the spectrum of Γ is identical to the one of the sample correlation matrix of independent spins, whose density of eigenvalues is given by the Marcenko-Pastur (MP) law, [START_REF] Johnstone | Proc. ICM 2006[END_REF]. The edges of the continuous component of the MP spectrum are given by

ρ MP (λ ′ ) = v(1-α) δ(λ ′ )+ α 2πλ ′ v (λ + -λ ′ )(λ ′ -λ -) (93) with v(u) = max(u, 0)
λ ± = 1 - 1 √ α 2 . ( 94
)
The largest eigenvalue of Γ, λ + , is not related to the value of L.

In the weak noise regime (α > α c ) the largest eigenvalue of Γ is [35]

λ 1 = L 1 + 1 α (L -1) . (95)
It exceeds L for any finite α, and converges to L when α → ∞. The rest of the spectrum is described by the MP density (93). Expression (74) for the squared norm b 1 of the orthogonal fluctuations leads to the analytical formula

b 1 = 1 α λ+ λ- dλ ′ ρ MP (λ ′ ) λ 1 -λ ′ = λ 1 -L λ 1 , (96)
where we have used the analytical expression of the Stieltjes transform of ρ MP [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]. Using (75) we deduce the value of the squared projection of the inferred rescaled pattern (ξ 1 ) ′ onto v 1 ,

a 1 = L -1 λ 1 . (97)
Identities ( 96) and (97) are graphically interpreted in Fig. 1: b 1 is the squared norm of the orthogonal fluctuations β, while a 1 is the squared projection of the rescaled pattern ξ onto v 1 .

The above discussion is illustrated on the simple case of a Hopfield model with p = 1, p = 0 patterns in Fig. 18, see caption for the description of the model. Using formula (91) we compute the largest eigenvalue of the correlation matrix for perfect sampling, L = 2. Figure 18 shows that a large eigenvalue clearly pulls out from the bulk spectrum for the ratio α = 4 (top spectrum), larger than the critical ratio α c = 1 according to (92) (bottom). For α = 4, the infinite-N predicted values for the largest eigenvalue, λ 1 = 2.5 (95), and for the edges of the MP spectrum, λ -= .25, λ + = 2.25 (94), are in good agreement with the numerical results for N = 100.

Formulae (96) and (97) hold for each pattern µ when p ≥ 2 patterns are present, provided that p remains finite when N → ∞. The case of p = 2 patterns, where one pattern is strong and has overlap q > 0 (81) with the sampled configurations, and the second pattern has weak components, is of particular interest. Again, we assume that the fields vanish. Repeating the calculation of Section VI A 2 and Appendix A we find that the entropy S/N quickly decreases with B from 2 bits down to 1 for B = O(1). When B ∝ N , the entropy decreases from 1 down to 0; the expression of S coincides with (86) where β is replaced with β(1 -q 2 ). Hence we have a two-step behaviour: the strong pattern is determined with O(1) examples, the weak pattern requires O(N ) sampled configurations. Learning of the weak pattern is possible if

α ≥ 1 ξ2 (1 -q 2 ) -1 2 , (98)
according to (87). The two-step behaviour agrees with the discussion of Section VI B 1.

Coexistence of ferromagnetic states

Consider now the case of the coexistence of two ferromagnetic states exposed in Section V B. Data are generated from a Hopfield model, with zero fields and one strong pattern ξ, as in Fig. 4. In the up-state the spins are magnetized with m + i = tanh(q ξ i ). In the down-state the local magnetization is m - i = -m + i . On the overall the local magnetization is m i = 1 2 m + i + 1 2 m - i = 0, up to O( 1 √ B ) fluctuations. The discrepancy between the Gibbs magnetizations, m i = 0, and the state magnetizations, m ± i , results in a O(1) contribution m + i m + j (= m - i m - j ) to the correlation matrix entry Γ ij , dominating the O( 1 N ) contributions due to the interactions between spins. The largest eigenvalue of Γ,

λ 1 = i (m + i ) 2 , (99)
is of the order of N ; the corresponding eigenvector is v 1 = (m + 1 , m + 2 , . . . , m + N )/ √ λ 1 . Informally speaking, the information about the state magnetizations is not conveyed by the Gibbs magnetizations (as in Section VI B 1) but by the correlation matrix [36]. According to formula (55) the pseudo-magnetization T i vanishes; hence we correctly infer that the fields h i have zero values. Using formula (9) we obtain

(ξ 0 ) i ≃ N λ 1 m + i . (100)
Therefore, the inferred pattern component is not equal to the true pattern component, but is proportional to its hyperbolic tangent. This non linear transform is clearly seen in Fig. 4. The discrepancy between the true and inferred components is a nice illustration of the claimed scaling for the higher order corrections in (49) (recall that the eigenvalues of A -1 are the p largest eigenvalues of Γ). In the presence of coexistent states, while ξ 2 is small compared to N , λ 1 is of the order of N , making the ratio

λ 1 ξ 2
N of the order of unity. Corrections are required and shown to improve the quality of the inferred pattern in Fig. 11.

VII. CONCLUSION

In this paper we have studied how to infer a small-rank interaction matrix between N binary variables given the average values and pairwise correlations of those variables. We have seen that the generalized Hopfield model, where the interactions are encoded into a set of attractive and repulsive patterns ξ, is a natural framework for Maximum Likelihood (ML) inference. Using techniques from the statistical physics of disordered systems, we have presented a systematic expansion of the log-likelihood in powers of λ ξ 2 N , where λ is the largest eigenvalue of the correlation matrix Γ (1). We have then calculated the ML estimators for the patterns and the fields to the lowest and first order in this expansion in a variety of physical regimes. The lowest order is a simple extension of Principal Component Analysis, where not only the largest but also the smallest eigenmodes build in the interactions. First order corrections involve non-linear combinations of the eigenvalues and eigenvectors of Γ. We have validated our ML expressions for the patterns on synthetic data generated by Hopfield models with known patterns and fields, and by Ising models with sparse interactions. We have also presented a simple geometrical criterion for deciding the number of patterns. Those results have been discussed and compared to previous studies in the unsupervised learning and random matrix literatures.

The quality of the inference strongly depends on the number of sampled configurations, B. The sampling error on each magnetization, m i , and pairwise correlation, c ij , is of the order of B -1/2 . Elementary insights from random matrix theory suggest that the resulting errors on the eigenvectors of the matrix Γ are √ N times larger. The error on the inferred patterns, ǫ, picks up a contribution ∼ N B 1/2 due to finite sampling, as found in Section II C. This scaling has several important consequences. First, inference is retarded: no information about the true couplings can be obtained unless the ratio B N exceeds a critical value (Sections VI A 2 and VI B 2). Secondly, for larger B, ǫ decreases as B -1/2 , which is confirmed by the simulations presented in Fig. 12, and then saturates to the intrinsic error resulting from our approximate expressions for the patterns. The intrinsic error depends on the order in the expansion used for the calculation of the cross-entropy in Section V. Note that other inference methods, looking for the local structure of the interaction network [START_REF] Cocco | [END_REF]12], may unveil strong cou-plings J = O(1) from a much smaller number of sampled configurations, B = O(log N ), and do not suffer from the retarded learning transition.

Our study could be extended in several directions. It would be particularly interesting to consider the case of spins taking Q > 2 values (Potts model), e.g. for applications to the study of coevolution between residues in protein sequences [23,25,[START_REF] Weigt | Proc. Nat. Acad. Sci[END_REF]. Mean-field inference methods provide a simple and efficient way to get interactions from correlations [START_REF] Weigt | [END_REF]. Knowing how MF interactions are modified when some eigenmodes are rejected (using the criterion of Section II D) or first-order corrections are taken into account would be of interest. However the linear increase in the number of possible symbols with Q (= 20 for amino-acids) may make the effective size of the problem, N × Q, larger than the number of configurations, B, in practical applications. A large number of vanishing eigenvalues is expected in those cases, and extracting repulsive patterns may become a difficult task.

Appropriate priors P 0 could also be used to force many pattern components to identically vanish, instead of acquiring small values as in Section II E. This can be particularly useful when the true patterns are known to be highly sparse and few data are available. Inspired by the so-called Lasso regression method [START_REF] Tibshirani | [END_REF], a natural prior is

P 0 ∝ exp -γ N i=1 1 -m 2 i p µ=1 |ξ µ i | + p µ=1 | ξµ i | .
(101) Contrary to the case of the quadratic penalty (21) the most likely values for the patterns cannot be expressed by means of simple analytical formulae. However, they could be efficiently obtained using convex optimization algorithms minimizing the sum of the cross entropy and of the penalty term (101).

Last of all, we have considered in this work that the configurations were sampled at equilibrium. In practice, when more than one state exist, the equilibration time may be prohibitive and a reasonable assumption would be to sample from one state only. To what extent ergodicity breaking in the sampling affects the quality of inference is an interesting question.

  entries of the M matrices are

FIG. 2 :

 2 FIG.2: Schematic behaviour of the error ǫ on the inferred patterns as a function of the number B of sampled configurations and for a problem size equal to N , when the pattern components are of the order of unity compared to N . See main text for the case of few large pattern components, of the order of √ N , i.e. couplings J of the order of 1.

FIG. 3 :

 3 FIG.3: Application of formula(9) to two sets of B = 40 (top) and 400 (bottom) configurations, randomly generated from the distribution PH (2) for model A with p = 1 pattern. The standard deviation of the pattern components is ξ = .7. Left: comparison of the true and inferred couplings for each pair (i, j). Right: comparison of the true and inferred components ξi of the pattern, with the error bars calculated from[START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]. The dashed lines have slope unity. Inference is done with p = 1 attractive pattern and no repulsive pattern.

FIG. 4 :

 4 FIG.4: Same as Fig.3, but with a standard deviation ξ = 1.3 instead of ξ = .7. The amplitude is strong enough to magnetize the configurations along the pattern, see Sections VI A 1 and VI B 3.

FIG. 5 :

 5 FIG.5: Relative differences between the true and the inferred couplings, ∆J ab /J ab for three system sizes, N . The inference was done using the lowest order ML formulae(9) for the patterns. Data were generated from Model B (perfect sampling); there are a priori ten distinct values of the couplings, one for each pair of blocks a and b. Inset: average value of ∆J ab /J ab as a function of 1 N . Circles, squares and diamonds correspond to, respectively, N = 52, 100 and 200 spins.

6 :

 6 FIG.

  FIG.

  Figure 7 compares the inferred and true couplings for B = 4500 sampled configurations. The generalized Hopfield model outperforms the standard Hopfield model (p = 0), showing the importance of repulsive patterns in the inference of sparse networks with strong interactions. Large couplings, either positive or negative, are overestimated by the lowest order ML estimators for the patterns.

  B. Error bars and criterion for p, p

FIG. 8 :

 8 FIG.8: Criterion to decide the number p of patterns and performance of the ML inference procedure for three different sizes of the data set, B. Left: inferred vs. true interactions with p = 1, 2 or 3 patterns; the dashed line has slope unity. Right: coefficients a µ = (ρ µ ) 2 and b µ = (β µ ) 2 vs. pattern index µ, and angles θ µ , divided by π 2 , see definitions (16) and[START_REF] Engel | Statistical Mechanics of Learning[END_REF]. For each value of B one data set was generated from Model A with p = 3 patterns, and standard deviations ξ 1 = .95, ξ 2 = .83, and ξ 3 = .77.

FIG. 9 :FIG. 10 :

 910 FIG. 9: Inferred vs. true couplings for Model D, with B = 4500 sampled configurations. Left: Generalized Hopfield model with (p, p) = (9, 10) and[START_REF] Cocco | [END_REF][START_REF] Tibshirani | [END_REF] (corresponding to the numbers of eigenvalues, respectively, larger and smaller than unity). Right: angles θ µ and θµ for, respectively, attractive (triangle) and repulsive (diamond) patterns.

FIG. 11 :

 11 FIG.11: True vs. inferred components of the patterns, ξi, for the model with N = 100 spins described in Fig.4. Full circles are the result of the lowest order inference formula (9), while empty circles show the outcome of the first order formulae(23).

FIG. 12 :

 12 FIG.12: Relative error between the inferred and true couplings for Model B (with N = 52 spins) vs. number of sampled configurations, B. The two curves correspond to the inference done with the 0 th order formula (9) (black circles) and the 1 st order formula (23) (squares). Each data point is the average over 10 samples; relative error bars are about 1%, and are much smaller than the symbol size. The asymptotic value of the errors, corresponding to perfect sampling (B = ∞), are extracted from Figs. 5 and 10.

1 FIG. 13 :

 113 FIG.13: Couplings calculated with the generalized Hopfield model vs. couplings calculated with the adaptive cluster expansion of[START_REF] Cocco | [END_REF] for 37 cells recorded in the prefrontal cortex of a behaving rat. Top: Hopfield model with p = 4, 8 (corresponding to the optimal number of patterns selected by the geometrical criterion) and 17; no repulsive pattern is considered (p = 0). Bottom. Generalized Hopfield model with (p, p) = (4, 4),(8,8) (optimal numbers) and[START_REF] Nokura | [END_REF][START_REF] Mackay | [END_REF] (corresponding to the numbers of eigenvalues, respectively, larger and smaller than unity).

IG. 14 :

 14 Couplings calculated with the Generalized Hopfield model versus coupling calculated with the adaptive cluster expansion for 92 amino-acids in the PDZ domain.The values of p, p are given in the Figure.Note that p = 0 for the top panels. The middle panels correspond to the optimal values for the number of patterns.

1 2FIG. 15 :

 115 FIG.15: Same as Fig.14when retaining the 44 residues with the largest weights Di only[24]. The values of p, p are given in the Figure.Note that p = 0 for the top panels. The middle panels correspond to the optimal values for the number of patterns.

FIG. 16 :

 16 FIG.16: Comparison between the couplings Jij calculated with all 92 residues and with the 44 most weighted residues only, for each one of the 44 × 43/2 pairs (i, j) of residues.

S

  [{σ b }] = -{ξi=± ξ} P [0, ξ|{σ b }] log P [0, ξ|{σ b }] , (79) therefore measures the uncertainty about this unique pattern given a set B sampled configurations. Intuitively, the dependence of S on B is closely related to the physics of the Hopfield model (with pattern ξ and zero fields) used to generate the examples. If the model is in the paramagnetic phase, i.e. if the components of the pattern are weak

FIG. 17 :

 17 FIG. 17: Entropy of the posterior distribution for the patterns, S (in bits and per component), as a function of the number of sampled configurations, B, when the local fields hi are known to vanish. (a). Ferromagnetic regime ( ξ2 = 1.1): the entropy decays exponentially with B. Inset: comparison with the theoretical prediction exp(-B/Bc) (dashed line), with Bc ≃ 6.85, in semi-log scale. (b). Paramagnetic regime ( ξ2 = .5): S (86) is a decreasing function of α = B/N . The entropies calculated from numerical calculations are shown for N = 10 and N = 20. Inset: the overlap r (83) between the inferred and true patterns is positive when α exceeds αc = 1 (87).

FIG. 18 :

 18 FIG.18: Spectrum of the correlation matrix for a Hopfield model with p = 1 pattern, N = 100 spins, and for B = 100 (bottom) and 400 (top) randomly sampled configurations at equilibrium. The bulk parts of the spectra coincide with the Marcenko-Pastur law for random correlation matrices. When B is large the top eigenvalue clearly comes out from the noisy bulk and the corresponding eigenvector approximately corresponds to the pattern. The pattern components are i.i.d. Gaussian variables, of zero mean and variance ξ 2 = .5; local fields hi have zero values.

Inference from the correlations: relationship with random matrix theoryWhat is the order of magnitude of Γ ij ? We first consider the ideal case of perfect sampling (B → ∞ while N is large but finite). As a result of the presence of the

Acknowledgments: We thank S. Leibler for numerous discussions. V.S. thanks the Simons Center for Systems Biology for its hospitality. This work was partially funded by the ANR contract 06-JC-JC-051.

Appendix A: Replica calculation of the entropy S for weak patterns

When the pattern has binary components ξi = ± ξ we make the change of variables σ ′ i = ξ i σ i to rewrite the partition function (35) of the Hopfield model through

where the inverse temperature β is defined in (80). The partition function is thus independent of the pattern direction, which makes the calculation considerably simpler. The posterior entropy (79) can be written as

where Equation (A2) gives the entropy of the system for a particular set of measures {σ b }. It is natural to expect the entropy to be reproducible across different sets of measurements. In this context, we are interested in evaluating the average of the entropy with respect to all possible measurements. Assuming that the configurations {σ b } are sampled from the equilibrium measure of a Hopfield model with one pattern ξ, we write the average entropy as

where

where we have introduced a new variable β since we should not take the derivative only with respect to β in (A4).

To calculate the average value of the logarithm of Ñ in (A5) we use the replica trick [27] and estimate the n th moment of Ñ ,

We introduce auxiliary Gaussian variables, denoted by mb , to linearize the quadratic term in the spins σ b i . We obtain, after summation over the spins,

In the paramagnetic phase we expect the variables m ρ b and mb to be of the order of 1 √ N . Expanding the hyperbolic cosine to the second order in those variables and carrying out the resulting Gaussian integral we obtain

Here, M is the (n + 1) × (n + 1) matrix with elements

with the overlaps defined through r ρσ = 1 N i ξ ρ i ξ σ i and t ρ = 1 N i ξ ρ i ξi . We now enforce the definitions of the overlaps using conjugated Lagrange multipliers, rρσ and tρ , and obtain

where Ξ is given by

We look for a replica-symmetric saddle point of Ξ: r ρσ = r, t ρ = t, rρσ = r and tρ = t. We obtain, after some elementary algebra,

where Dz = dz e -z 2 /2 / √ 2π is the Gaussian measure and

We now send n to zero. The saddle-point equations show that t = r; this result was expected from the fact that, if β = β, the true pattern ξ plays the role of an extra replicated pattern ξ. In addition, t = r ≡ γ, where γ is defined in (85). The self-consistent equations for r and the entropy S are given by, respectively eqns (84) and (86). on the N -spin configuration through the four block magnetizations (sums of the N 4 spins in each block) only. Hence, the correlations cij and magnetizations mi can be calculated in a time growing as N 4 (instead of 2 N ), which allows us to reach sizes equal to a few hundreds easily.

[41] The corresponding magnetizations were ≃ (-.26, .13, .13, .23) for N = 52 spins.

[42] Formula (90) can be found by inverting identity (12), with Jij = 1 N ξiξj.