{'fraction_non_alphanumeric': 0.036658227066415534, 'fraction_numerical': 0.009664371507086177, 'mean_word_length': 4.4742869216141985, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the typical properties of inverse problems in statistical mechanics

In this work we consider the problem of extracting a set of interaction parameters from an high-dimensional dataset describing T independent configurations of a complex system composed of N binary units. This problem is formulated in the language of statistical mechanics as the problem of finding a family of couplings compatible with a corresponding set of empirical observables in the limit of large N. We focus on the typical properties of its solutions and highlight the possible spurious features which are associated with this regime (model condensation, degenerate representations of data, criticality of the inferred model). We present a class of models (complete models) for which the analytical solution of this inverse problem can be obtained, allowing us to characterize in this context the notion of stability and locality. We clarify the geometric interpretation of some of those aspects by using results of differential geometry, which provides means to quantify consistency, stability and criticality in the inverse problem. In order to provide simple illustrative examples of these concepts we finally apply these ideas to datasets describing two stochastic processes (simulated realizations of a Hawkes point-process and a set of time-series describing financial transactions in a real market).


Introduction
The generations living during the last twenty or thirty years witnessed a huge scientific revolution which has been, essentially, technology driven. An impressive amount of computational power became cheaply available for people and institutions, while at the same time the quantity of data describing many aspects of our world started to grow in a seemingly unbound fashion: the human genoma can be efficiently sequenced in some days [75,88], the interactions among proteins in a human body can in principle be enumerated one-by-one [69], financial transactions are recorded with resolutions well below one second [1], the dynamics of networks of all kinds (social, economics, neural, biological) can be tracked in real-time. Parallel to this, the widely accepted scientific paradigm according to which it is necessary to ground reliable models on solid first principles started to crumble: promising results evidenced that it is possible to extract accurate statistical models from empirical datasets without even trying to guess what is their underlying structure, nor to characterize which input-output relations govern their behavior. Large datasets can be automatically and faithfully compressed in small sets of coefficients [31], their features can be described accurately with unsupervised algorithms, new data can be predicted with a given degree of accuracy on the basis of the older one (see for example [3]). Google uses pattern recognition and Bayesian techniques to translate from one language to the other regardless of the formal rules of the grammar [2], and Netflix can predict how much you will rate a movie (one to five) with an error around 0.85 without knowing anything about you but a few of your former preferences [4]. The embarrassing success of this approach compels a basic epistemological question about modeling: does an approach based solely on statistical learning lead to any actual understanding? What does one learn about a system when processing data in this way?
This problem is particularly relevant when dealing with the task of high-dimensional inference, in which a typically large set of parameters is extracted from an even larger dataset of empirical observations. What meaning has to be associated with each of the many parameters extracted from data? Are there combinations of such numbers describing global, macroscopic features of the system? A prototypical example is provided by the study of networks of neurons, in which one would like to understand how the brain works (e.g., the presence of collective states of the network, the possibility to store and retrieve informations) by processing data describing the behavior of a huge set of elementary units (the neurons). This task can be thought of as a seemingly hopeless one: in a way it is similar to reverse-engineering how a laptop works by probing the electric signal propagating through its circuitry. A modern answer to this type of arguments is the idea that if data is sufficient and the inference algorithm is good enough, some of the actual features of the system will eventually be detected. In the case of a laptop, one can think to extract from data not only the wiring pattern a set of cables, but to detect collective features such as the fact that a computer is an essentially deterministic object (in contrast to biological networks, where fluctuations are essential), or that it possesses multiple collective states (say, switched-on, switched-off or sleepy).
Physics, and in particular statistical mechanics, has much to do with all of this picture for two main reasons. The first one is technical: while the high-dimensional limit is a relatively new regime in the field of statistical inference, statistical mechanics has since long developed mathematical descriptions of systems composed by a very large (or better, infinite) number of interacting components [40]. Hence, mapping problems of statistical inference onto problems of statistical mechanics opens the way to a remarkable amount of mathematical machinery which can be used to solve quickly and accurately problems which become very complicated for large systems [45,79]. This is even more true since the study of heterogeneous and glassy materials produced sophisticated tools (replica trick, cavity methods) suitable to study systems in which no apparent symmetry or regularity is present, as often found in data describing complex systems [56]. The second, and more philosophical, reason is that statistical mechanics is naturally built to explain collective behaviors on the basis of individual interactions. Just as the ideal gas can be understood by studying the aggregate behavior of many non-interacting particles, or the emergence of spontaneous magnetization can be derived by studying the interactions of single spins, statistical mechanics can be used to predict the collective behavior of biological, social and economic systems starting from a given set of rules describing the interaction of some fundamental units [30]. In 1904 Ludwig Boltzmann, almost a century before anyone could take him literally, anticipated that "The wide perspectives opening up if we think of applying this science to the statistics of living beings, human society, sociology and so on, instead of only to mechanical bodies, can here only be hinted at in a few words." Hence, from the perspective of (large-scale) statistical learning, it is natural to use statistical mechanics methods to study the emergence of collective properties of a system once the microscopic interactions of the fundamental units have been reconstructed through a careful analysis of empirical data.
Unfortunately, even if one is able to do that, it is not always easy to understand how much of the inferred model faithfully describes the system: it is possible, and it is often the case, that the procedure which is used to perform the data analysis influences so much the outcome that the actual properties of the system get lost along the way, and the inferred model shows a spurious behavior determined just by the fitting procedure. For example, models with binary interactions may describe very well systems in which the interaction is actually multi-body [39], just as critical models (strongly fluctuating statistical systems) may fit random set of observables much better than ordinary ones [53]. Noise itself may be fitted very well by sophisticated models, while non-stationary systems might be accurately described by using equilibrium distributions [84]. In all of these cases, it is important to develop quantitative tools which allow to distinguish between spurious features of the inferred model and genuine ones.
The purpose of this work is precisely to inquire some of those aspects in the simpler setting in which we consider a statistical system consisting in a string of N binary variables, used to model T independently drawn configurations. We will show that, while the small N regime the problem of inference can be completely controlled (chapter 2), the large N regime becomes computationally intractable and non-trivial collective properties may emerge (chapter 3). Such features be observed independently of the data, and have to be associated uniquely with the properties of the model which is used to perform the inference, regardless of the system which one is trying to describe.
In chapter 4 we will show under which conditions the problem of inferring a model is easy, showing in some cases its explicit solution. We will also evidence the limits of non-parametric inference, highlighting that for under-sampled systems correlations might be confused with genuine interactions. In chapter 5 we will provide a geometric interpretation for the problem of inference, showing a metric which can be used to meaningfully assess the collective phase of an inferred system. We will apply these ideas to two datasets, describing extensively the results of their analysis in the light of our approach.
In this chapter we will describe the problem of extracting information from empirical datasets describing a stationary system composed of a large number of interacting units. Interestingly, this problem has almost simultaneously received a great deal of attention from the literature of diverse communities (biology [77,87], genetics [22], neuroscience [72,76,24], economy, finance [48,59,29], sociology). This can be traced back to two main reasons: first, it is now possible across many fields to analyze the synchronous activity of the components of a complex system (e.g., proteins in a cell, neurons in the brain, traders in a financial market) due to technological advantages either in the data acquisition procedures or in the experimental techniques used to probe the system. Secondly, data highly resolved in time is often available, which (beyond implying that finer time-scales can be explored) provides researchers with a large number of observations of the system. Defining as N the number of components of the system and as T the number of available samples, these last observations can be summarized by asserting that the limit of large N and large T can be accessed for a large number of complex systems. In this work we will restrict ourselves to the more specific case in which such systems are described by binary units, reminding to the reader that (i) most of what will be shown can be generalized to the case of non-binary (Potts) or continuous variables [85] and (ii) the binary case already allows to describe in detail several systems [72,76,24]. In section 2.1 we describe the models that we consider, which usually go under the name of exponential families and are justified on the basis of the maximum entropy principle (appendix A.1), and state the direct problem, alias the calculation of the observables given the model. In section 2.2 we present the problem of inferring a model from data (the inverse problem) and characterize it as the Legendre-conjugated of the direct one. In section 2.3 we present the regularization techniques which can be used to cure the pathological behavior of some inverse problems and improve their generalizability. Although the results presented in this chapter are far from being original, we aim to show as transparently as possible the deep connections between information theory and statistical mechanics, emphasizing the strong analogy between direct and inverse problems.

The direct problem
We introduce in this section the direct problem -which deals with finding the observables associated with a given statistical model -as a preliminary step towards the formulation of an inference problem. This is the problem typically considered by statistical mechanics, hence we will adopt most of the terminology and the notation from this field. The main results that we will present are associated with the free energy -which we use in order to generate the averages and the covariances of the model -and to its relations with the notion of Kullback-Leibler divergence and the one of Shannon entropy. Finally, we will characterize the large and small deviation properties of the empirical averages of the observables under the model.

Statistical model
We consider a system of N binary spins s = (s 1 , . . . , s N ) ∈ {−1, 1} N = Ω, indexed by i ∈ V = {1, . . . , N }. A probability density p is defined as any positive function p : Ω → R such that s p(s) = 1, while the space of all possible probability densities on Ω is denoted as M(Ω). We also consider a families of real-valued functions φ : Ω → R |φ| with components φ(s) = (φ 1 (s), . . . , φ |φ| (s)), which will be referred as binary operators, and are more commonly known in the literature of statistical learning as sufficient statistics or potential functions [85], and will be used in order to construct a probability density on the configuration space of the system. g µ φ µ (s) (2.2) and is referred as the partition function of the model. The free energy F (g) is defined as F (g) = − log Z(g).
For conciseness, the identity operator will always be labeled as the zero operator φ 0 (s) = 1, in order to reabsorb the normalization constant Z(g) into its conjugated coupling g 0 . The probability density will be written as p(s) = p s , so that (2.1) will be compactly written as With these definitions, the coupling g 0 results equal to the free energy g 0 = − log Z(g) = F (g). Given a family of operators φ, we also denote as M(φ) the set of all the statistical models of the form (2.1) obtained by varying the coupling vector g. Given the probability density (2.1) and a generic subset Γ ⊆ V (which we call a cluster ), we also define the marginal p Γ (s Γ ) as which expresses the probability to find spins belonging to the Γ in a given configuration once the degrees of freedom associated with spins outside such cluster have been integrated out (whereas p ∅ = 1 and p V (s) = p(s)).
This construction will be used to study inference problems in which the M operators φ(s) are a priori known. We will disregard for the moment the issue of optimally selecting the most appropriate operators in order to describe a given set of data, an important problem known as model selection. Let us indeed remind that models of the form (2.1) can be justified on the basis of the maximum entropy principle, which will be stated in appendix A.1. The next notions which will be defined are the one of ensemble average and the one of susceptibility which will be extensively used throughout our discussion. while the generalized susceptivity matrixχ is defined as the covariance matrix whose elements are given by Beyond describing fluctuations around the ensemble average of the φ operators, the generalized susceptibilityχ is a fundamental object in the field of information theory [27], in whose context is more often referred as Fisher information, and is more commonly defined as Its relevance in the field of information theory and statistical learning will later be elucidated by properties (2.23) and (2.24) which concern with the direct problem.
Proposition 2.1. The free energy function enjoys the properties and χ µ,ν = − ∂ 2 F ∂g µ ∂g ν , (2.9) thus it is the generating function of the averages and of the fluctuations of the operators φ µ contained in the model.
Equation (2.9) implies that covariances χ µ,ν are related to the response of the ensemble averages with respect to changes of the couplings through a relation known as fluctuation-dissipation relation, which is a direct consequence of the stationary nature of the probability distribution (2.1). Another fundamental property of the free energy function F (g) is its concavity, which will later allow us to relate the field of statistical inference with the one of convex optimization (appendix C). It can be shown (appendix A.2) that: • The susceptibility matrixχ is a positive semidefinite matrix, thus the free energy F (g) is a concave function.
• If the family of operators φ is minimal (i.e. it doesn't exist a non-zero vector x such that µ x µ φ µ,s is constant in s), then the susceptibility matrixχ is strictly positive definite and the free energy F (g) is strictly concave.

Entropy and Kullback-Leibler divergence
In this section we will define the concept of Shannon entropy, which will be used as an information theoretic measure of the information content of a distribution. The quantity S(p) measures the amount of disorder associated with the random variable s, and satisfies the following properties: • 0 ≤ S(p) ≤ log |Ω|. In particular S(p) = 0 for p(s) = δ s,s (when the variable s is maximally informative), while S(p) = log |Ω| for the flat case p(s) = 1/|Ω| (in which s is maximally undetermined).
• The function S(p) is concave in p.
They can be proven straightforwardly, as for example in [27]. Another informationtheoretic notion which will be extensively used is the Kullback-Leibler divergence D KL (p|q), which characterizes the distance between two probability distributions.
Although it doesn't satisfy the symmetry condition nor the triangular inequality required to define a proper measure of distance, in chapter 5 we will show that indeed a rigorous concept of distance can be extracted by means of the Kullback-Leibler divergence. Such quantity enjoys the following properties: • D KL (p||q) ≥ 0 for any pair of probability densities p, q.
• D KL (p||q) = 0 if and only if p = q.
• D KL (p||q) is a convex function in both p and q.
These property justify the role played by the Kullback-Leibler divergence in information theory, and can be proven straightforwardly (see [27]). Notice indeed that given two statistical models (φ, g) and (φ, g ) respectively associated with densities p and p , the entropy and the Kullback-Leibler divergence can be written as so that the concavity properties of S(p) and D KL (p||q) can be related to the ones of the free energy F (g). These quantities will be relevant in order to characterize the large deviation properties both for the direct and of the inverse problem.

Observables
Throughout all our discussion, we will focus on the case in which T independent, φ µ (s (t) ) (2. 16) are called empirical averages. It is worth remarking that P T (ŝ) depend on the observed configurations just through the empirical averagesφ. We will denote averages over the measure P T (ŝ|g) with the notation . . . T . We also define the empirical frequencies (also known as type)p as the vector with components which enjoys the following properties: • It is positive and normalized ( sp s = 1), thus it defines a probability density on Ω (i.e.,p ∈ M(Ω)).
• If the datasetŝ is generated by a probability distribution p, thenp is distributed according to the multinomial distribution where T s = Tp s . Its first and second momenta are Finally, given a collection of operators φ we will denote the set of all empirical averages φ that are compatible with at least one probability density in Ω with which is called in the literature marginal polytope [85]. It can be proven that (see for example [85]): • G(φ) is a convex set (i.e., givenφ,φ ∈ G(φ), for any α ∈ [0, 1] also

Small and large deviations
In the case of the direct problem it is natural to formulate the following questions: 1. What are the most likely values for the empirical averagesφ?
2. How probable it is to find rare instancesŝ?
The first question is relatively easy to answer, and characterizes the role of the generalized susceptivity in the direct problem as ruling the convergence of the empirical averages to the ensemble averages 1 , as shown in the following and proven in appendix Proposition 2.3. Given a statistical model (φ, g), the empirical averagesφ satisfy the relations The explicit form of the likelihood function (2.15) allows to answer exhaustively also to the second question.
Proposition 2.4. Given a probability density p defined by a statistical model (φ, g), This implies that the probability of observing dataset a genericŝ decays exponentially in T , with a non-trivial rate function I p (p) determined by the empirical 1 In the framework that we are considering (i.i.d. sampling of configurations drawn by the same distribution) empirical averages always converge to ensemble averages with an error scaling as 1/ √ T . Indeed it makes sense to model the case in which the probability measure p breaks into states, so that for any finite length experiment, just samples belonging to the same state are observed. This is meant to model the phenomenon of ergodicity breaking, which we will comment about in section 3.4. averagesφ only. Also notice that the large deviation function can be expressed entirely in terms of the entropy and the Kullback-Leibler divergence as (2.25)

The inverse problem
In this section we introduce the inverse problem of extracting a coupling vector g given a set of operators φ and a vector of empirical averagesφ. We will present this problem as dual with respect to the direct one, showing that just as the knowledge of the free energy F (g) completely solves the direct problem, the Legendre transform of F (g) denoted as S(φ) and characterized as the Shannon entropy, analogously controls the inverse one.

Bayesian formulation
We will be interested in calculating the set of couplings g which best describes a given set of dataŝ of length T within the statistical model (φ, g). Bayes theorem provides a mathematical framework in which the problem can be rigorously stated, by connecting the likelihood function P T (ŝ|g) described in section 2.1.3 to the posterior of the model P T (g|ŝ), which specifies the probability that the dataŝ has been generated by model g. Bayes theorem states in fact that where P 0 (g) is known as the prior, and quantifies the amount of information which is a priori available about the model by penalizing or enhancing the probability of models specified by g by an amount P 0 (g). Bayes theorem also links the concept of prior to the one of regularization which will be discussed in section 2.3, but for the moment we will consider the prior P 0 (g) to be uniform (i.e. a g-independent constant), so that it can be reabsorbed into the pre factor of equation (2.26). In this case finding the best model to describe the empirical averages may mean: • Finding the point in the space of couplings g in which the function P T (ŝ|g) is maximum (maximum likelihood approach).
• Finding the region of the space of couplings in which such probability is high (Bayesian approach).
These two approaches lead to very similar results in the case in which the likelihood function is strictly concave, as one can prove by means of large deviation theory (see section 2.2.4 and appendix A.6). Roughly speaking, when the number of observations T is large, the posterior P T (g|ŝ) concentrates around the maximum likelihood parameter, being the rate of convergence fixed by the stability matrix of the maximum and the number of samples T . Hence we will later define as the inverse problem the characterization of the maximum likelihood parameters and of their linear stability, disregarding the detailed shape of the function P T (g|ŝ).

Maximum likelihood criteria
The maximum likelihood criteria requires to find the maximum of the likelihood function P T (ŝ|g), whose solution is obtained by differentiation of equation (2.15) with respect to the couplings g µ , and reads for each µ

Statement of the inverse problem
The concavity properties of the likelihood function (or equivalently, of the free energy F (g)), allow for a characterization of the problem of inferring the maximum likelihood parameters g given dataŝ in terms of a Legendre transform of F (g).
Definition 2.6. Given a minimal set of operators φ and a set of empirical averages φ, the function S(φ) is defined as the Legendre transform We denote with g the (only) value of g maximizing equation (2.28). Such quantity By construction the statistical model (φ, g ) verifies the matching condition (2.27).
By considering the Shannon entropy S(p) = − s p s log p s and by plugging probability density p inside its definition, one can see that it holds so that we can prove the following proposition.
Proposition 2.6. Given a vector of empirical frequenciesp, the large deviation function for the inverse problem Ip(p) ∝ − 1 T log P T (p|p) is given by the Kullback-Leibler divergence (2. 34) This implies that the probability for datap to be generated by any model p decays exponentially fast in T with a rate function given by the large deviation function D KL (p||p). This result can be seen as a particular case of a more general theorem, which is known as Sanov theorem and whose proof can be found in appendix A.4. 2 Theorem 2.1. Consider a statistical model defined by a probability distribution p, and a (compact) set of probability densities M ⊆ M(Ω). Then ifp is a vector of 2 We won't adopt the informal version of the theorem often found in literature (see for example [54]), which doesn't require the introduction of the set M . In such form the theorem is not valid when, for any value of T , M has empty intersection with the set of realizable empirical frequencies, as the probability for any point in M to be realized is strictly zero regardless of T . empirical frequencies sampled from the distribution P T (ŝ|p), it holds that Building on these results, we can provide an answer for our first question and find out what are the most likely distributions p having generated datap. In particular, it is possible to expand the Kullback-Leibler divergence around its minimum and perform a saddle-point estimation, obtaining the following result.
Proposition 2.7. Consider a generic datasetŝ defining the empirical distribution p ∈ M(Ω). Then, given a family of operators φ, the posterior probability (with uniform prior) P T (g|p) ∝ P T (p|g) defines a probability measure on space M(φ), parametrized by the coupling vector g which defines the statistical model p. The averages and the covariances under this measure are given in the large T limit by where g is the maximum likelihood estimator of g andχ −1 is the inverse of the Fisher information matrix calculated in g .
This result (proved in appendix A.6) characterizes the inverse of the generalized susceptibility as the matrix quantifying the speed in T at which the probability measure on the inferred couplings concentrates around the maximum likelihood estimate.
The centrality of this matrix in the inverse problem is also provided by a rigorous bound that can be proven for the covariance of any unbiased estimator, and known as Cramér-Rao bound. From this perspective, χ −1 µ,ν , can be seen as establishing a bound to the maximum rate of convergence for the estimator of a coupling.
Theorem 2.2. Consider a statistical model (φ, g) with F (g) strictly concave and an unbiased estimator of the couplings g (i.e., such that g µ T = g µ ). Then the covariance matrix of g under the measure . . . T satisfies where withX Ŷ we indicate that the matrixX −Ŷ is positive semidefinite.
The proof of this theorem is presented in the appendix A.5.

Examples Independent spins model
The simplest model of the form (2.1) which can be considered is of the form and will be called independent spin model. The model contains N operators of the form {φ {i} (s) = s i } i∈V (called in the following magnetizations), whose conjugated couplings are denoted as g {i} = h i (and referred as external fields). The empirical magnetizations will be denoted ass i = m i . The direct problem can be solved by evaluating the partition function of the model, so that the free energy F (h) results The ensemble averages and generalized susceptibilities can be obtained by differentiation, and are given by The inverse problem is also easily solvable, as the Legendre transformation of F (h) can explicitly be computed, and the entropy results while by differentiation one finds The additivity both of the entropy and of the free energy, which are crucial in order to solve the model, descend directly by the independence of p(s), which can be written as a product of single spin marginals (2.46) Notice that the existence of the solution is guaranteed for any m in the hypercube [−1, 1] N , while its uniqueness is enforced by the minimality of the operator set (which is additionally an orthogonal set in the sense that will be defined in (4.1)). As expected, for m i = ±1, the estimator h i is divergent, so that h i (m i = ±1) = ±∞.

The pairwise model
The next model that will be presented is known in a large variety of fields with different names (Ising model in physics, graphical model in the field of statistical learning), and is defined by the probability density where E is a given set of edges, that is, a given subset of While in statistical mechanics it has been extensively used since 1925 as a prototypical model to study magnetic materials [41,15], it has deserved a special interest in the field of statistical learning as it is the simplest model which is able to capture the correlation structure of a given dataset 3 . The operator content of this model is a set of N magnetizations, conjugated to their corresponding external fields (as in section 2.2.5), and a set of |E| ≤ N (N −1) 2 operators {φ {i,j} (s i , s j ) = s i s j } (i,j)∈E conjugated to a set of pairwise couplings g {i,j} = J ij . We will call empirical correlations the averages Remark 2.3. This direct problem for the pairwise model is hard to solve in the general case for even moderate values of N , in the sense that the calculation of the partition function Z(g) is a problem which is known to belong to the #P-complete class [43,42]. Only for some subclasses of this general problem an exact, analytical solution for the partition function can be obtained (e.g., regular lattices, trees) and evaluated in polynomial time, while in general just approximate solutions can be obtained in polynomial time [42]. Another possible approach consists in finding approximated expressions for the partition function Z(h,Ĵ ) which are proven to converge in the limit of large system size or weak interaction to the exact result for the free energy of the model (mean-field approximations).
In the next sections we will introduce specific versions of model (2.47) for which we will be able to solve the inverse problem, namely the fully connected ferromagnet (section 3.3) and the pairwise tree (section 4.2.3).

The regularized inverse problem
The inverse problem described in section 2.2 may appear extremely easy to solve due to the concavity of the free energy F (g). The optimization of concave functions is usually very easy because fast algorithms such as gradient ascents can find in short time a maximizer (if any) for F (g) (appendix C). Despite that, there are several cases in which this procedure may be problematic, so that the function F (g) is often replaced by a modified function F (g) − H 0 (g) which enforces a better behavior for the inverse problem. In this case the function H 0 (g) is called a regularizer. In a Bayesian setting, regularization can be understood as an injection of a priori information about a statistical model. Indeed the issue of regularization is a topic of fundamental importance in the field of statistical inference well beyond the need of enforcing mathematical tractability of the model. In particular it can be used to deal with these cases: • Divergencies: Regularization can cure divergencies, by removing infinite couplings. A solution to any inverse problem is guaranteed to exist for any set of empirical averagesφ ∈ G(φ), but such solution may be located at the boundary of the coupling space, in whose case one or more couplings are divergent.
Penalizing large couplings with a regularizer ensures that the inferred couplings attain a finite value. This is often the case for neurobiological or protein data  [7] or the Bayesian information criterion [73]), in order to tune the regularizer to a correct value (see also section 5.1.3).
• Model selection: Finally, regularization can be used as a tool to perform model selection. In the case in which data are distributed according to a specific, unknown, statistical model, it is possible to perform inference by using a more general distribution which is likely to contain (or to be very close) to the true one. By adding a suitable regularizing term (such as an L-1 or L-0 norm) it is sometimes possible to recover the original model as a particular sub-class of a more general distribution. For example, this has been used in the context of graph reconstruction, where models defined by specific topologies have been successfully selected by a regularizer out of the space of all possible graph structures [63,86].

Bayesian formulation
Consider an empirical datasetŝ and a model defined by a set of operators φ. Then the posterior of the model can be written as in (2.26), in which it is P T (g|ŝ) ∝ P T (ŝ|g)P 0 (g), so that the problem of inference can be reformulated as the minimization of the function Definition 2.8. Given a statistical model (φ, g) and a positive prior function P 0 (g) we define a regularizer as the function H 0 (g) = − log P 0 (g).
Notice that due to convexity of theχ matrix, if the regularizer H 0 (g) is (strictly) convex, also H(g|ŝ) is (strictly) convex. Hence, the introduction of a strictly convex prior can be used to remove zero modes from theχ matrix thus enforcing a unique solution for the inverse problem. In our analysis we will restrict to the case of convex regularizers. Also notice that if H 0 (g) = +∞ when any component of g is divergent, the solution of the inverse problem is confined to a finite region of the coupling space. Given a statistical model (φ, g), a set of empirical frequenciesφ and a vector β such that it is component-wise β µ > 0, we consider the minimization problem

Two popular regularizers
which we call the L-2 regularized inverse problem. This choice for H 0 (g) enforces strict concavity of the problem and finiteness of the values of g , which should satisfy the set of equationsφ This regularization corresponds to the Gaussian prior P 0 (g) ∝ exp − µ βµ 2 g 2 µ .
Notice also that the regularizer is differentiable, so that a solution of this problem can be addressed efficiently by using techniques such as the ones described in the first part of appendix C. As in the non-regularized case, the main computational limitation consists in calculating the gradient of the minus log-likelihood function − log P T (φ|g), which requires the knowledge of the averages φ as functions of the coupling vector g. This regularization procedure is typically used to remove infinite values arising in the solution of the non-regularized inverse problem.

L-1 regularization
We also present the L-1 regularized inverse problem, which is defined by the minimization problem corresponding to the choice of an exponentially decaying prior exp − µ β µ |g µ | ∝ P 0 (g) for the coupling vector g. Analogously to the L-2 case, this regularizer is convex and enforces a finite value for the inferred couplings g . Unlike that, this regularizer is non-differentiable. This introduces some difficulties in the solution of the minimization problem, as shown in the second part of appendix C, where it is shown that the inferred coupling vector should satisfy the equation where sgn(x) is the set-valued function defined by equation (C.11). The main interest in this regularizer arises from its efficacy as a feature-selector, as it is able to provide sparse solutions for g , i.e., to put exactly to zero some components of the inferred couplings vector. Despite being first used in the field of compressed-sensing, in which the use of the L-1 regularizer has been exploited to solve underconstrained sets of linear equations [31], this regularized has been successfully applied in the field of binary inference (also called logistic regression), in which it has been useful to reconstruct the structure of an interaction network of the form (2.47) [63,86] and even in more general cases dealing with non-binary interaction [71].
Remark 2.4. The two regularizers presented so far are special cases of the L-p regularization scheme, which is associated with the choice H 0 (g) ∝ µ β µ ||g µ || p , where ||x|| p = |x| p is the L-p norm of x. Notice that the L-p regularizer is convex (hence leading to computationally tractable minimization problems) for p ≥ 1, and is strictly so for p > 1. In particular, the L-1 regularizer can be seen as the simple ( alias, convex) regularizer that is closer to the L-0 one, which is associated with the problem of minimizing the number of inferred parameters for a fixed value of the posterior, a criterium which one would think to use in order to minimize the complexity of the inferred model.

Examples Independent spins model
Consider the model defined by the probability density (2.39). Then we can consider the regularized inverse problems in which one tries to minimize to the L-2 and L-1 norm. In the first case the (decoupled) set of equations which has to be solved in order to find the vector h is show that the inferred couplings h i attain a finite value for any of −1 ≤ m i ≤ 1 and 0 ≤ β i < ∞. In the case of the L-1 norm, one has to solve the decoupled set of equations whose solution is The solution for h i in the two cases for a specific value of β i is plotted in figure 2.2.
Notice that: • Both regularizations schemes produce a finite h i in the case |m i | = 1.
• The zero-field solutions of the L-1 regularized problem can be seen as arising from a complexity related criteria, stating that operators which do not add enough descriptive power to the model should be suppressed by the assignment of zero weight to their conjugated coupling. In this example the notion of "enough descriptive power " is quantified through the comparison of β i against the direc- Despite its trivial solution, we have chosen to present this problem as it shows with simplicity some basic features of the L-1 and L-2 regularizers which are retained even in more complicated scenarios.

High-dimensional inference
Rigorous results in information theory -such as the ones presented in section 2.2.4 -are able to provide both qualitative and quantitative understanding of the inverse problem in the regime of finite N and large T , the case most of the literature on statistical learning deals with, while computational techniques such as the ones described in appendix C provide efficient means to find its solution. Nevertheless, recent technological advances in several fields (such as biology, neuroscience, finance, economy) are pushing the fields of statistical learning towards a less trivial regime, in which both N and T are large, with a given relation among system size and number of samples keeping fixed their scaling. The reason for this change of perspective is that it is now possible for several complex systems to record a large number of data samples describing simultaneous the activity of the many microscopic constituents [22,72,76,24,48,59,29]. The question that naturally arises in this case is whether it makes sense to consider a model with a large (possibly very large) number of parameters, if the data available is also very large. The answer is non-trivial, and requires the addition of some degree of complexity to the problem of inference. The first problem which has to be addressed (section 3.1) is of purely technical nature, and deals with the problem of finding the minimum of a convex function when its gradient is computationally intractable. Then, we will describe some interesting conceptual problems which arise when considering the large N limit. For simplicity, we will consider initially the problem in which both N and T are large, but the number of inferred parameters M is finite (section 3.2). Discussing the case in which M scales with N as well will require the introduction of the notion of disorder, which we will briefly comment about in section 3.5.

Computational limitations and approximate inference schemes
In appendix C we show how it is possible to construct algorithms which are guaranteed to find a minimum (if any) for a convex function. Then the solution of the inverse problem can be written as a minimization problem over a convex function H(g) of the form that problem is in principle solved. Indeed, the problem which often arises in many practical cases is that the naive minimization of this function can be extremely slow, and ad-hoc techniques have to be implemented in order to overcome this problem.

Boltzmann Learning
One of the most intuitive algorithms to solve the inverse problem is provided by the Boltzmann learning procedure [6], which consists in the application of algorithm C.1.1 to the inverse problem described in section 2.2. In that case, the minimization procedure of H(g) consists in constructing a succession {g (k) } K k=1 of the form where { k } K k=1 is a schedule satisfying the set of conditions (C.4) which enforce the convergence of g (k) to the minimum (if any) g . Indeed the computation of each of the g (k) requires the evaluation of H(g (k) ) and the calculation of a gradient of the form The calculation of the gradient (or the sub-gradient) of H(g) requires evaluating the ensemble averages of the operators φ, which is a computationally challenging task if N is even moderately large. This is true even when the function H and the ensemble averages φ are not computed via direct enumeration (which would in principle entail a summation over 2 N states for each of the M operators plus the identity), and are instead calculated with Monte Carlo methods. The number of iterations required to calculate each of the gradients and the function H with a controlled precision is in fact typically fast growing in N , being the quality of the approximation and the time computational power required to obtain it dependent on the algorithm which is adopted to compute the averages (see for example [6,54,50,47]). Summarizing: • Boltzmann learning is able to solve with arbitrary precision any inverse problem.
• The computational power required to solve the inverse problem through the Boltzmann learning procedure with a given degree of accuracy (i.e. H(g (K) ) − H(g ) smaller than a fixed ) grows fast in N .

Mean field approaches for pairwise models
An alternative approach to the Boltzmann learning procedure can be constructed by adopting so-called mean-field techniques, which allow to obtain efficient approximations for the free energy F (g) and the averages φ of a statistical model. Such techniques are suitable for systems whose partition function can be quickly, although approximately, evaluated with a precision which either increases with the system size CHAPTER 3. HIGH-DIMENSIONAL INFERENCE N or decreases with the magnitude of the interactions, so that in many practical applications the difference between the approximated observables and the exact ones is very small [68,67]. For pairwise models of the form (2.47), mean-field approximations are well-known since long time in statistical physics. In particular we will consider approaches in which the free energy of the model (2.47) is expanded in a series around a non-interacting or a weakly correlated model (naive mean field, TAP approximation, Sessak-Monasson approximation), or obtained by assuming a factorization property of the probability distribution in terms of one and two body marginals (Bethe approximation). We will briefly describe these approximate inference schemes without providing explicit derivations, supplying the interested reader with the necessary ref-

erences.
In order to motivate the mean-field approach, we first state the result [62].
Proposition 3.1. Consider a pairwise model of the form where β > 0 is an expansion parameter. Then its free energy can be written as where the terms ∂ n F β ∂β n are functions such that: (i) depend only on the couplings J ij and the ensemble magnetizations s i (ii) for n ≥ 1 the n-th term involves n-th powers of J ij (iii) the ensemble magnetizations satisfy the self-consistency equations .

CHAPTER 3. HIGH-DIMENSIONAL INFERENCE
Leaving aside the problem of convergence of the series (3.5), the free energy for a generic pairwise model can in principle be obtained by setting β = 1 in the above expansion.
• Naive mean field: The naive mean field approximation can be obtained by truncating the series (3.5) for n = 2, thus obtaining the expression while the self-consistency equations become The solution of the inverse problem within this inference scheme can be obtained by inserting the momentum matching condition s i = m i in the previous expression, yielding a first set of relations among h ,Ĵ and m. Matching the correlations c ij with the ensemble averages s i s j requires instead the use of linear response theory 1 [45], which can be used to to prove that Putting those informations together, one finds that • TAP approximation: The Thouless-Anderson-Palmer (TAP) approximation can be obtained by considering an additional term in the expansion (3.5), often denoted as Onsager reaction [81], leading to the expression for the free energy and the self-consistency relation 2 (3. 13) Also in this case, in order to apply this approximation to the inverse problem [79], one has to use the momentum matching conditions together with linear response theory, leading to the expression [64] (ĉ − mm T ) −1 While the expansion (3.5) is a series for F (h,Ĵ ), and is hence associated with the direct problem, it is also possible to find an analogous expansion for the entropy S(m,ĉ) due to Sessak and Monasson which is more naturally associated with the inverse problem [74].  analytically resumed in [74], where it is found that their contribution leads to which is commonly referred as the Sessak-Monasson approximation.
Notice that the expansion (3.16) automatically leads to a series expansion for the external fields and the couplings by using relation (2.31) and exploiting the linearity of the derivative, without the need of resorting to linear response theory.
A different type of approximation is the so-called Bethe approximation, in which the free energy is written as Notice that for generic systems, the self-consistence equations are not guaranteed to yield a unique, stable solution, being the solutions to the minimization conditions associated with fixed points of the so-called Belief-Propagation (BP) algorithm for constraint satisfaction problems [54]. The expression for the averages obtained by using the free-energy (3.18) is given by [64] (3.20) • Bethe approximation The use of linear response theory together with equation (3.20) allows to find a solution of the inverse problem in Bethe approxima-tion, yielding where δc =ĉ − mm T . Notice that this equation describes the fixed point solution of the susceptibility propagation algorithm (SuscProp) [55] without the need of numerically iterating the algorithm itself [64].
Remark 3.1. The techniques described above have been extensively used in order to solve the inverse problem for the pairwise model. Indeed no general result for the quality of these approximations is rigorously known, thus it is worth remarking that (i) several approximations have been tested on synthetic and experimental data (see for example [24,68,67,52,64,12,26]) in order to check their performance and (ii) those approximations describe the correct expression of the free energy for some specific models. In particular the free energy (3.7) is the exact free energy for the (either homogeneous or heterogeneous) Curie-Weiss model in the limit of large N , (3.12) is the correct free energy for the Sherrington-Kirkpatrick model [62] and the Bethe approximation is exact for loop-less graphs (appendix D.2).

The large N , finite M regime
We will be interested in sketching some features of the inverse problem which arise for large values of N (a regime known in statistical mechanics as the thermodynamic limit), and in commenting about their role in the solution of an inference problem such as the one described in section 2.2. In particular we will consider the following issues: • Loss of concavity: A model defined by a strictly concave free-energy F (g) may develop null-modes associated with the matrixχ. This implies that the solution of the inverse problem may lose its uniqueness or, more precisely, large regions of the space M(φ) might be associated with similar sets of empirical averagesφ.
• Model condensation: Models undergoing a so-called second order phase transition display a divergence of one or more components of the generalized susceptibility matrixχ. This indicates that large portions of the marginal polytope G(φ) can be described by slightly shifting the values of g around the critical point in whichχ diverges. More generally, even for non-critical points finite regions of the space of the empirical averages can be mapped by the inverse problem onto sets of apparently vanishing measure of the space M(φ). We call this behavior model condensation, a phenomenon which will be discussed in great detail in chapter 5.
• Ergodicity breaking: The probability measure p may break in a set of P states, each of them characterized by a different probability density p (α) (with p = P α=1 q α p (α) and P α=1 q α = 1). If this is the case, empirical averages produced with a finite amount of data T |Ω| by any realistic dynamics concentrate according to the measure p (α) rather than the full measure p. Then, equation (2.23) fails to hold and the sampled averages are no longer representative of the global probability measure. Hence, the notion of ergodicity breaking deals with the direct problem more than with the inverse one, as it relates to the problem of the convergence of the averagesφ to the empirical ensemble averages φ . As the discussion of this phenomenon will require the addition of some structure to the direct problem, we will briefly comment its role in section 3.4.
Those features are expected to be universal, i.e., present in several models in the limit N → ∞ limit. Nevertheless, we will just study a single model known as the fully connected ferromagnet, and try to underline the characteristics which are expected to generalize also to other type of models.

Fully-connected ferromagnet
We want to illustrate some of the features described above by discussing a completely solvable model. Such model is a particular case of the pairwise model (2.47), and is also known as the Curie-Weiss model of magnetism. It has been used as a prototypical model to study the emergence of a spontaneous magnetization in ferromagnetic materials, as it is one of the simplest statistical models which are able to describe a thermodynamic phase transition between a non-ordered phase and an ordered one.
Definition 3.1. Consider the pair of operators φ = i s i , 1 N i<j s i s j , and the statistical model (φ, g) defined by g = (h, J), so that its associated probability density is given by We call this model a fully connected ferromagnet. As for the pairwise model, we will Due to symmetry, we will consider without loss of generality the model in the region h ≥ 0. The free energy of the model F (h, J) can be calculated in the large N limit using a saddle-point approximation, and can be written as

The mean-field solution
The solution of the direct problem considering only F 0 (h, J) will be called mean-field solution. Notice that due to the scaling F 0 (h, J) ∝ N , for large values of N this contribution dominates the free energy F (h, J).
Proposition 3.3. For all i = j the mean-field solution for the fully connected ferromagnet is : while the susceptibility matrix is given by This implies the following fact concerning the inverse problem.
while the couplings belong to the space restricted to the region in which sign(m) = sign(h). Finally, the inverse susceptibility matrix is divergent.
This last fact can be understood by checking that the matrix χ 0 has eigenvalue . In particular, the null eigenvalue has eigenvector (−m, 1) which indicates that the mean field solution of the direct problem is invariant under the change of couplings Thus, the inverse problem maps all the points belonging to the one-dimensional region (m, m 2 ) on the two-dimensional plane (h, J). This apparently contradicts the remark in section 2.2 about the existence of solutions to the inverse problem for any point belonging to the marginal polytope G(φ). Indeed, we will show in the next section that keeping properly into account the presence of the h = 0, J > 1 line allows to understand this discrepancy. Interestingly, the two-dimensional region

Finite N corrections
Keeping into account the terms F f luct and F trans allows to describe the transition from the finite N regime to the mean-field one. In particular, the Gaussian fluctuations around the mean-field solution extend the region in which the inverse problem is solvable to a strip of finite width in the space G(φ).
Proposition 3.5. Given (N m, N −1 2 c) ∈ G(φ), the inverse problem for a fully connected ferromagnet described by the terms F 0 and F f luct of equation (3.24) has solution if and only if c = m 2 + δc N with δc finite, and reads 3 . The null eigenvalue of the matrixχ 0 is lifted to a finite value, as one can see that Proposition 3.6. The inverse problem for the fully connected ferromagnet described by the terms F 0 + F trans has solution for any point (m, c) ∈ G(φ) excluding the region 3 In the literature concerning the so-called inverse Ising model, this result is typically derived by differentiating the relation with respect to m j , and by recognizing that through linear response one can write (∂h/∂m Also in this case one can show that in the limit h N , the null mode ofχ 0 is Finally, one can draw the following conclusion, which despite being a trivial consequence of what shown above, shows that the N → ∞ limit can lead to counterintuitive results. This simple example shows some of the features discussed above concerning the limit of large N , namely: 1. The free energy loses (strict) concavity, as one has detχ −−−→ N →∞ detχ 0 = 0. This indicates that some directions in the coupling space cannot be discriminated. In this example, when N is large, interactions are no longer distinguishable from external fields due to the presence of an eigenvector (−m, 1) associated with the null eigenvalue.
2. Model condensation takes place, as all the region G(φ) but a set of null measure is mapped on a one-dimensional strip. This will be better elucidated in chapter 5, where we will be able to quantify the density of models contained in a finite region of the space (h, J).

Saddle-point approach to mean-field systems
In this section we generalize the procedure employed in the case of the fully connected ferromagnet to the case in which a saddle-point approach is used to solve the direct problem for a generic system. In particular, we consider a statistical model (φ, g) with partition function 42) and suppose that the operators φ µ,s can be written as functions of a small set of Then it is possible to write where e Σ(ψ) = s δ (ψ − ψ(s)), and Σ(ψ) is often referred as entropy for the value of the order parameter ψ. For many statistical models, one has that the limit is finite, and f (ψ) is often called (intensive) free-energy for the value of the order parameter ψ. In this case, one can exploit a saddle-point approximation to evaluate the partition function Z(g) at large N . It results .
(3. 45) where we use the notation f (n) (ψ) for the tensor with components f (n) , and ψ is the global minimum of the function f (ψ), which in Besides providing us with a mean to calculate the free energy F (g), the ensemble averages φ and the susceptibilitiesχ, the notions defined above allow us to introduce the concept of state, which we will use to characterize the phenomenon of ergodicity breaking.
Definition 3.2. Consider a statistical model (φ, g) which can be described by a set of order parameters ψ, and such that at large N its partition function can be approximated by (3.45). Then we call a state any local minima of the saddle-point equations (3.46).
We will label any of those minima as ψ (α) with α = 1, . . . , P , and use the a superscript α to identify quantities associated with the state α, as for example In particular we informally remind that for large statistical models (φ, g) endowed with a realistic dynamics (e.g., Metropolis-Hastings [80,36]) leading in the limit of exponentially large T to the stationary distribution p associated with (φ, g), states naturally emerge when observing a finite amount of configurations. In fact, the iteration of a dynamics for T 2 N time steps typically produces configurations belonging to the same state as the initial one, while in the opposite limit of large T the probability of observing a state belonging to a configuration α is proportional to e N f (ψ (α) ) .
Hence, unless data obtained from an experiment are exponentially large in the size of the system (which isn't typically the case in real world applications of the inverse problem), one expects empirical averages to concentrate around averages which are in principle different from the ensemble ones, and that are associated with a specific state α. Accordingly, we define the notion of state average φ (α) , which is expected in the regime of T 2 N to model the averages obtained by experiments of finite length as follows: Definition 3.3. Given a system (φ, g) whose partition function can be approximated by the partition function (3.45), we define the state averages and the state susceptibilities The correctness of above construction has been verified for several statistical models subject to different dynamics [57,89], nevertheless to the best of our knowledge no fully general, rigorous result concerning this phenomenon is available yet. In particular, in order to rigorously motivate the notion of state average, it would be necessary to show that for a generic, local dynamics a decomposition property of the where P α=1 q α = 1 and p (α) ∈ M(Ω) holds for the Gibbs measure, which again is known to be correct just for specific models.
In that case, the state averages and the susceptibilities can be explicitly computed by explicitly deriving the above free-energy, allowing to prove the following result.
Proposition 3.7. The direct problem for a statistical model (φ, g) which can be described with order parameters ψ and an order parameter free-energy f (ψ) can be solved in saddle-point approximation in any state α, leading to (1) (1) f φ (2) µ;e,f f ) and by convention repeated index are summed.
This result allows us to characterize the behavior of the inverse problem in the large N limit. In fact one can see that at leading order in N , the momentum matching condition (2.27) becomes where we remark that the averages in the state α do not depend explicitly on g, being their dependence contained in the order parameter ψ (α) . This implies that, given two statistical models (φ, g) and (φ, g ) with g = g such that there exist a couple of states (respectively α and α ) solving the saddle point equations with ψ (α) = ψ (α ) , in the large N limit those models cannot be discriminated.
Remark 3.4. Consider an empirical datasetφ (α) generated by a system in state α.
Unless one doesn't consider a matching condition in which the state average contains the corrections of order 1/N indicated in the right term of formula (3.51), it is not generally guaranteed that it is possible to reconstruct the state α which generated the empirical averages. pairwise couplings. In this case one may have several problems in studying the features introduced in section 3.2 as we did above. In particular: • The averages m,ĉ and the generalized susceptibilityχ are hard to compute for a generic value ofĴ and h. Therefore, it is not possible to understand which points of marginal polytope G(φ) are associated with zero modes inχ.
Moreover, the limit N → ∞ is ambiguously defined if no prescription is provided for how should the empirical averages scale with N .
• For the same reason, it is not possible to find in which points one expects model concentration to occur, as this would require knowing which eigenvalues ofχ are divergent in the thermodynamic limit N → ∞ for generic points (m,ĉ) ∈ G(φ).
• No saddle-point approach is justified for generic empirical averages m,ĉ. Thus, an approach analogous to the one in 3.4 cannot be considered, and the notion of state cannot be described in such terms.
These difficulties could be overcome by resorting to the notion of disorder, which is commonly used in the field of statistical mechanics of heterogeneous systems. In particular we want to show, as a possible outlook of this work, an approach to the analysis of the large M limit borrowed from that field [56] which could be applied to this problem.

Self-averaging properties and inverse problem
Given an operator set φ, consider a set of statistical models M(φ) and a prior P 0 (g) on this space. Then, suppose that a statistical model (φ, g) is sampled according to P 0 (g), and successively a set of empirical data of length T is drawn by such distribution. Several functions of the estimator g (φ) can be built in order to analyze the properties of an instance of the inverse problem, such as the quantities which quantifies the average error in the inferred coupling and whose divergence signals critical properties of the generalized susceptibility matrixχ.
If these of quantities are self-averaging for large N and T (i.e., they concentrate around an average value determined by P 0 (g)), then one expects that specific instances of inverse problems drawn by the same prior P 0 (g) to share the same collective features.
As an example, if one considers a Gaussian prior for the ferromagnetic model of the Remark 3.5. The idea of disorder in the context of the inverse problem is obviously linked to the existence of a prior P 0 (g) on the space M(φ), so that in principle the case of a flat prior cannot be treated with these techniques. Nevertheless, fixing implicitly a specific class of models through P 0 (g) is the price to pay to answer to very interesting questions, which wouldn't otherwise be well-posed namely: (i) can a specific model be learnt with high probability according to a given inference prescription? (ii) Are the global properties of an heterogeneous system equivalent the the ones of an homogeneous one? (iii) Is it possible to understand the generic properties ofχ?

Complete representations
In this chapter we will introduce the notion of complete family of operators, which can be used to gain some insight about the inverse problem. Although in general this approach may introduce a high degree of over fitting, dealing with complete families allows to discuss very transparently some features of inference which are related to algorithmic complexity (section 4.2.1). Moreover, completeness allows for an explicit reparametrization of the probability distribution (2.1) in terms of state probabilities, allowing for a complete understanding of properties of the inverse problem which are less clear by using the Gibbs form for the probability density. More interestingly, in this language we will be able to differentiate local features of the direct and of the inverse problem, which in turn rely on the locality of the marginals. In this chapter the inverse problem for some models will be exactly and explicitly solved, while some ideas will be presented in order to generalize this methods to more relevant problems

Orthogonality and completeness
We define in this section the notion of orthogonality and completeness for families of operators. While the orthogonality condition is related to the one of minimality, the one of completeness will allow to formally invert the relation among ensemble averages φ and couplings g.
Definition 4.1. Given a family of operators φ, we call it orthogonal if it satisfies Property (4.1) can be seen as expressing the fact that in an orthogonal family any pair of operators decorrelate when averaged with respect to a uniform probability density (at infinite temperature in the language of statistical mechanics). Additionally, if φ 0 ∈ φ, one can see that in an orthogonal family, for µ = 0 i.e., φ µ has zero mean at infinite temperature for any µ = 0. Finally, if φ is an orthogonal family, then it is easy to see that φ\{φ 0 } is minimal. The main result that derives instead from equation (4.2) is the explicit one-to-one mapping between couplings g, state probabilities p and averages φ , as clarified by the next proposition.
Proposition 4.1. Given a family φ satisfying (4.1) and (4.2), the statistical model (φ\{φ 0 }, g) associated with the probability density p satisfies Additionally, state probabilities can be expressed as Proof. These relations are a direct consequence of the axioms (4.1) and (4.2) and can be checked by direct substitution.

Monomials
Throughout most of the following discussion, we will focus on families of operators φ formed by monomials, for which axiom (4.1) trivially applies. More precisely, given a cluster of spins Γ, we define the monomial φ Γ (s) as while the identity is associated with the empty cluster φ 0 (s) = φ ∅ (s) = 1. It is easy to show the following: • φ is an orthogonal family; • φ is complete if and only if it contains all possible monomials, whose number Moreover, monomials satisfy a very important relation which will be used extensively in the following.
Proposition 4.3. Consider a complete family of monomials φ. Then the marginals of the probability density p associated with the model (φ\{φ 0 }, g) can be expressed as Proof. This can be checked by using equation (4.6) and showing that for each mono- This property expresses the locality of marginals once they are expressed in terms of ensemble averages. This should be compared with the expression of a marginal written as a function of the couplings (2.4), in whose form the locality properties are hidden by the interaction structure.

Inference on complete models 4.2.1 The complete inverse problem
The techniques shown in the above section can be used to write a formal solution of the inverse problem in full generality. The main drawback of this procedure is the overfitting issue which has to be associated with the presence of an exponential number of couplings, which in practical cases makes this approach unfeasible unless the system has small size (typically N ∼ 10 1 ). Indeed, as the solution of the complete inverse problem illustrates with simplicity some very general features of many inverse problem, we choose to present its solution.
Definition 4.2. The complete inverse problem is the inverse problem associated with the statistical model defined by the complete family of monomials {φ Γ (s)} Γ⊆V . Its probability density can be written as It is easy to write the formal solution for the entropy by using the relation (4.2), while its differentiation (or the direct use of the relation (4.5)) leads to an exact expression for the couplings and the susceptibility matrix.
Proposition 4.4. The expression for the entropy of the complete inverse problem reads while the inferred couplings g and the inverse susceptibility matrixχ −1 result (4. 13) This solution has a simple interpretation in terms of empirical frequencies, once one rewrites above expression using the relationp In this form it is possible to appreciate that the solution simply corresponds to a matching of state probabilities with empirical probabilities. which can be used to express the log-likelihood function as a function of the probabilities p instead of the coupling vector g. Its maximization can be seen equivalently as performed over the state probabilities p. In this case, the obvious solution is p =p, so that the expression (4.14) is describing an approach in which the state probabilities are matched with the empirical ones one-by-one. In particular, if a configuration is not observed, the inferred probability for that configuration is strictly zero.

Divergencies
The formal solution (4.14) shows that the inferred couplings can be infinite if there are states which are never sampled in the dataŝ. In particular if data are generated by an actual probability distribution p assigning zero weight to some configuration, the Ω space splits into an accessible and a non-accessible sector, and divergencies can be seen as are required to implement an hard constraint on the set of accessible configurations. Couplings obtained by using this scheme are finite either when all states are measured or when divergencies cancel out for a given region of the coupling space. Indeed, the presence of an unaccessible sector has to be considered a spurious result unless p s ≈p s , which is expected to hold just in the large T limit. In particular for T < |Ω| ,p s = 0 for at least |Ω|−T > 0 configurations, regardless of the presence or absence of a forbidden sector. Therefore it is not possible to distinguish if divergencies are due to the presence of an unaccessible sector or to poor sampling. In this case, regularization schemes such as the use of Laplacian smoothing or an L-2 norm can be used to obtain finite results. This basically corresponds to lift the probability for non-measured configurations from zero to some finite value. For example, Laplacian smoothing procedure [70] corresponds to the choice: Finally, we remark that the same type of divergence arises in all the cases that will be analyzed (see sections 4.2.3 and 4.2.4), and is a very general characteristic of inverse problems, which typically relates to under sampling. This is the simplest setting in which this problem can be analyzed in full generality.

Observed sector
The expression for the inferred couplings (4.14) involves a summation over all the configuration space Ω, so that a summation over |Ω| = 2 N terms seems to be required to calculate any of them. Indeed, those expressions may be rewritten exploiting the orthogonality relation (4.1), which implies that whereĪ = {s ∈ Ω |p s > 0} is the set of observed configuration. Then, one can rewrite (4.14) as where the term proportional top 0 = 0 account for the divergencies, and the sum over states runs over a number |Ī| ≤ T terms. In the case of the regularized complete inverse problem (section 4.2.2), we will see that it will be possible to write an analogous expression for the couplings, in which the weight assigned to non-observed configuration will be finite.

Rate of convergence
Given an underlying statistical model p for the complete inverse problem, large deviation theory (as described in section 2.2.4) states that for large T the variance of the inferred couplings g with respect to the measure given by P T (p|p) ∝ P T (p|p) is Incidentally, the same quantity can also be obtained by averaging with respect to the . . . T measure, a result which allows to express the rate of convergence for the complete inverse problem (appendix D.1). While the 1/T pre factor expresses the expected scaling for the error on the inferred coupling, the χ −1 Γ,Γ term is non trivial. In particular we observe that: 1. The fluctuations of the inferred couplings are identical for all the operators.
2. The value of the fluctuations is bound by the inequality: where p min = min s p s .
3. The speed of convergence is limited by the presence of rare configurations. In particular if p min = 0, the variance diverges.
The generalization to the case in which the sector of observable states I = {s ∈ Ω | p s > 0} is smaller than the entire phase space I ⊂ Ω is straightforward (appendix D.1). Indeed, it is necessary to define a set of regular operators φ reg such that

Regularization of the complete inverse problem
The generality of the complete inverse problem renders its regularization relevant for a strong theoretical reason. In fact, the complete inverse problem is totally nonparametric in the sense that the probability distribution (4.10) contains all possible statistical models describing a set of N binary variables. Then one could think of selecting the most appropriate statistical model to describe a dataset of binary data simply by applying a suitable regularizer to this general problem, and let the regularization term itself perform the task of model selection (an approach successfully adopted in [63,86] in a less general scenario). We present in the following the results obtained by using different regularization terms, and comment about the interpretation of the solutions of the regularized inverse problem. Finally, we will characterize a symmetry property of regularizers which can be used to study their suitability in the field of high-dimensional inference (i.e., for large values of N ).

L-2 regularization
The simplest regularized version of the complete inverse problem is the one defined by the function which implements the Gaussian prior over the L-2 norm of the coupling vector described in section 2.3.2. In terms of state probabilities, equation (4.24) can be written as and its minimization with respect to p s (constrained to s p s = 1) leads to the set of implicit equations Its solution determines the value of the couplings g through the relation We observe that: 1. The summation over the configuration space requires considering in principle an exponential number of terms, but this issue can be avoided as explained in section 4.2.1.
2. The expression for g µ is always finite, as the presence of infinite couplings is suppressed by the cost associated with the L-2 norm.
3. The parameter β controls the total value of the L-2 norm of the coupling vector g and the entropy of the inferred distribution. In particular the total L-2 norm can be expressed as This considered, the expression for the couplings obtained using this regularization scheme is where the p s and the p 0 satisfy the set of implicit equations: (4. 30) We remark that the calculation of the regularized couplings can be performed in polynomial time in T .

Entropy regularization
Another choice for the regularization is motivated by the following argument. If a datasetŝ of length T is associated with an entropy S(p) ∼ log T , with log T N , it is likely for the model to be in the under sampled regime, as the entropy per variable is expected to be finite (i.e., S(p) ∼ N ) for well-behaved models. Then, it is possible to consider a regularizing term which penalizes low entropy distribution, so that where as usual p is the density associated with the statistical model (φ, g), so that S(p) = − M µ=0 g µ φ µ . The minimization of above expression with respect to g µ leads toφ (4.32) After some manipulation and after using the completeness relation (4.2) one finds Finally, by writing s s = −p s log p s , one is led to a set of implicit equations (4. 34) which is analogous to the one described in the L-2 case. Also in this case the system has to be solved numerically, by exploiting the fact that the probabilities p s depend on the s index through the empirical frequencyp s (i.e., states visited the same number of times are associated with the same inferred probability). Equation (4.2) can finally be used to extract the inferred couplings from the probability density of p s .

Susceptibility regularization
The inverse generalized susceptibility of a modelχ −1 provides an indication of the generalizability of an inference procedure through equation (2.32), which implies that the response of the inferred couplings g to a shift of the empirical averagesφ is (4. 35) Then one could think to favor generalizability in an inference procedure by introducing a regularization term of the form (4.36) By employing equation (4.13) it is easy to see that the inverse susceptibility matrix can be written as a function of the coupling vector g as so that the total energy can be written as Its minimization leads tō whose solution requires solving a set of implicit equations analogous to (4.30) and (4.34) of the form (4.40) By using equation (4.6) the solution p can be used to explicitly express g . The minimization of above expression with respect to g leads to where we define the set valued function sgn(x) as in appendix C. We remark several issues concerning this regularizer: 1. Unlike the L-2 case, the completeness relation does not allow to switch from a summation on operators to a summation over configurations, hence algebraic properties cannot be fully exploited to manipulate the above equation.

Explicit selection of couplings
An interesting case is the one in which a specific inverse problem -such as the inverse pairwise model -is seen as a regularized version of the complete inverse problem.
This implicitly implies that unlike with the previous regularizers, in this particular example we are not interested in the problem of model selection, but we mean to offer a different perspective on a problem which is known to be hard, in order to characterize it from a different point of view. In particular, we consider the regularized minus-log-likelihood 43) in which θ ∈ {0, 1} M determines the couplings that are penalized by the L-2 norm, and we consider the limit of large, positive β, so that g µ ≈ 0 if θ = 1. The minimization of (4.43) leads to the set of equations which in the parametrization of states becomes The last term in (4.45) is finite in the limit of large β, and encodes the constraint specified by θ. Within this formulation the intrinsic difficulty of an inverse, noncomplete problem emerges as the fact that the probabilities p s can be different for states visited with the same frequency. This is associated with the dependence of the second term of equation (4.45) upon the index s associated with the operators φ µ,s , and is analogous to the case of the L-1 norm described above.
Remark 4.3. A formal solution for this problem can be written by studying the limit β → ∞, which is associated with couplings g µ = 0 for θ µ = 0. The equation g µ = 0 can be expressed in term of operator averages by using equation (4.5) as follows: where φ ν indicates the ensemble average of the operator φ µ under the distribution p s .
This result expresses a relation among observables which must hold whenever couplings are zero, which is typically used to express higher order correlations in terms of low order ones, which can be expressed as roots of polynomial equations. Then equation (4.5) can be used to write the remaining couplings, and the roots of equation (4. 46) can in principle be used to obtain an expression for the non-zero components of g .
Symmetry properties of the regularizers   Proof. To prove the above proposition it is sufficient to notice that by symmetry the coupling vector g depends on the two values p 0 and p 1/T associated with the states sampled zero (p s = 0) and once (p s = 1/T ). Then equation (4.6) implies that Remark 4.4. The symmetry broken by the L-1 regularizer and by the explicit coupling selection is associated with the following consideration: in principle, unless there is an explicit information that allows to distinguish between states s and s that are observed the same number of times, then inference should assign the same weight to those states. In the first case such symmetry is spontaneously broken (the information injected by the prior doesn't specifically favor any state), while in the second it is explicitly broken.

Pairwise model on trees
One of the simplest cases in which the pairwise model defined by equation ( 50) in the case in which the set of edges E does not contain any cycle. Then this model is called a tree (see appendix D.2 for a more precise definition).
For such models the inverse problem is easy to solve due to the factorization property shown in appendix D.2, which allows to write the probability density as Hence, the entropy can be written as and the inverse problem can be solved, as shown in the next proposition.
while the fields h and the couplings J result and the inverse susceptibility matrixχ −1 is given by The structure of this solution is reminiscent of the one shown in the case of the complete inverse problem described in section 4.2.1, and can intuitively be understood as follows. To solve an inverse problem it is necessary to find the clusters which allow to express the entropy (in that case all clusters had to be included, while in this case single spin and two spins clusters alone are sufficient to write the full entropy). The presence of a large number of delta functions is due to the fact that the entropy is built by a small number of cluster contributions, so that the response of the couplings to a shift in the value of the conjugated average is strongly localized: either the perturbation is applied to a neighbor, in whose case the response is finite, or it is zero. This has to be compared with the direct problem, in which a perturbation in the couplings changes the average of a finite number of operators in general. In that case, one roughly expects that 56) where ξ is the correlation length of the system. This was first noted in [25,26], where it is shown that for a large number of statistical models that the structure ofχ is dense, while the one ofχ −1 tends to be sparse.  where T is a set of translation operators characterized by a periodicity parameter ρ and φ is a set of M operators of range R.

One
We are interested in solving the inverse problem for this type of system, which means to calculate the entropy S(Tφ) as a function of the empirical averages of the operators T φ = N/ρ n=0 T n φ. In order for the entropy to be well-behaved, and in order to exploit the property of completeness (4.2), we need to require a specific choice for the set φ. This relation, together with equation (4.9) which expresses the locality of marginals, allows to explicitly find the expression of the entropy of a one-dimensional chain.
Proposition 4.7. The inverse problem for an orthogonal, complete one-dimensional chain of monomials has the following solution. The entropy can be expressed as 1 where c µ,Γ = 1 if ∃n such that T n φ µ = φ Γ and c µ,Γ = 0 otherwise. The couplings while the inverse susceptibilities are given by . Remark 4.6. In order to apply these ideas to empirical data, the information about the one-dimensional nature of the problem should be a priori known. Indeed, the exact nature of the interactions needs not to be known, provided that the R parameter is larger than the actual range of the interactions.

Complete inverse problem
The techniques shown in section 4.2.1 have been tested on synthetic datasets in order to check their performance. As expected, they are suitable for systems in the small N regime, due to the slow convergence in T of the inferred coupling vector g to the true coupling vector g, which can be seen as a consequence of the over fitting problem associated with the presence of an exponential number of couplings. We have considered for simplicity a system of N = 8 spins, with couplings corresponding to several models, namely:    Finally, we show how this procedure might be employed in the case in which one or more couplings are infinite. We consider complete models in which all couplings g Γ are put to zero, but a random set which are set to g Γ = ∞.

Pairwise model on a tree
We tested the results shown in section 4.2.3 providing a solution for the inverse problem for pairwise models with tree-like structure. We considered trees of size  We also studied how the quality of the reconstruction of the couplings degrades by raising the β parameter. We find that within this inference scheme it is possible to reconstruct accurately the couplings as long as local fluctuations are sampled. More

One-dimensional periodic chain
We studied the performance of the inference procedure described in  on βg is minimum, while for the parameters g the reconstruction error is minimal for a finite value of β which optimize the signal-to-noise ratio. We remark that also in this case the quality of the reconstruction of the couplings is determined by the sampling of the configurations belonging to clusters of R spins. If local fluctuations are not sampled well-enough, the error on the inferred couplings is large as predicted by equation (4.63). As observed above, it is not necessary to probe global fluctuations of the system in order to accurately reconstruct the couplings.

Chapter 5
Information geometry and criticality In this chapter we will be interested in studying the natural structure of Riemannian manifold which characterizes the space of probability distributions [9]. This structure provides a mean to rigorously define a distance between statistical models, which can be used to characterize the consistency of the solution of the inverse problem through the notion of distinguishable distribution [61]. The metric structure of the coupling space becomes especially interesting in the case of models displaying a critical behavior at large N , as it allows for a characterization of (second-order) criticality from the point of view of information theory. In this scenario critical points can be seen as regions of the space of statistical models which are infinitely descriptive, in the sense that any finite region of the coupling space around a critical point can encode an anomalously high number of distinguishable statistical models. We call this phenomenon model condensation. An illustrative example is presented by discussing the thermodynamic limit of a fully connected ferromagnet. Finally, we will introduce a model of a stochastic point-process known as Hawkes process which we will use as a toy model to study the features of the inverse problem when applied to a realistic dataset, and compare the results to the ones obtained by studying real data describing financial transactions in a stock market. This will allow to distinguish among spurious and the genuine collective features which emerge from the analysis of empirical data similar to the one considered in [72,76,24] in the context of neurobiology. in the field of information geometry, in which the geometric properties of the space of probability distributions are inquired by using methods of differential geometry (see [9,10] and [11] for a pedagogical review), which we will briefly present in the following sections. We will be interested in using these methods to answer several M µ=1 X µ ∂ µ and Y = M µ=1 Y µ ∂ µ and a point p ∈ T (φ) one can define the scalar product ·, · p : T (φ) × T (φ) → R as: It can be shown (appendix A.2) that for any X, Y ∈ T (φ) one has Hence, ·, · p is a metrics which we define the Fisher metrics associated with M(φ).
Notice that the scalar product X, Y p is independent of the parametrization used to describe the distribution p due to the transformation law of χ µ,ν = ∂ µ log p(s)∂ ν log p(s) . This fact, and the choice of this metric itself, will be intuitively justified in the next section where the notion of distinguishable distribution will be introduced. The Fisher metrics allows to define the length of a curve in the space M(φ).
Definition 5.2. Given a curve γ, i.e., a one-to-one function γ : with components γ = (γ 1 , . . . γ M ), we define its length as It is easy to show that the length of a curve (i) is independent of the parametrization of γ, (ii) is independent of the parametrization of M(φ) (iii) is additive, i.e., 1 It is customary in literature to use superscripts for contravariant tensors and superscripts for covariant ones. We will disregard for simplicity this distinction and use lower indices for any tensor or vector field, as their use will be unambiguous. is called a geodesics, and its coordinates γ = (γ 1 , . . . , γ M ) satisfy the linear differential equation where the Christoffel symbols Γ µ ν,ρ are given by In appendix (E.1) we prove this well-known result by explicitly varying the length functional (γ).
is a proper measure of distance.
We will show in the next section that this distance relates to the inverse problem by intuitively counting how many error bars away are two distributions away one from the other, given a fixed experiment length T . A related concept is the one of volume, which can be used to quantify the number distributions that cannot be distinguished one from the other on the basis of an experiment of finite length. The relevance of this measure will be elucidated in section 5.1.3.

Sanov theorem and distinguishable distributions
The metric introduced in section 5.1.1 can be justified by providing an intuitive interpretation in terms of distinguishable distribution, a concept which we will present starting from a simple consistency requirement. Suppose to be given a datasetŝ of length T generated by an underlying (unknown) distribution. Then, given an operator set φ it is possible to construct the empirical averagesφ and to infer the maximum likelihood estimate of the couplings g = g (φ) describing the data, and to use them to generate a different datasetŝ of the same length asŝ. The maximum likelihood estimator g = g (φ ) ofŝ will, in general, be different from g. Thus, distributions labeled by g and g cannot be distinguished on the basis of a dataset of length T , as sketched in figure 5.1. What one expects is that by increasing T , the model g gets closer and closer to g. This idea can be rigorously formulated by means of Sanov theorem (presented in section 2.2.4), which allows to prove the following corollary. Corollary 5.1. Consider a statistical model (φ, g) associated with a probability density p. Then, given a set of empirical averagesφ generated by p and a maximal likelihood estimator g , the probability that the maximum likelihood estimator g takes a value close to g on the dataset associated withφ is given by where p is defined by the statistical model (φ, g ) and δg = [−δ, δ] M .
The proof of this corollary is presented in appendix (E.3). What it implies is that the Kullback-Leibler divergence controls the probability that after the resampling procedure explained above one ends in a model very different from the starting one.
As expected, such probability is exponentially small in T . We will informally rewrite above corollary in the form implying a choice of δ enforcing δg very close to 0. This will allow us to characterize the concept of indistinguishable distribution.
Definition 5.5. Consider two models g and g within the same family of operators φ. Then, given a dataset of length T and empirical averagesφ sampled by the model (φ, g) and an accuracy > 0, we say that g and g are indistinguishable if the maximum likelihood estimator g satisfies Given corollary (5.11), it is easy to prove (appendix E.3) that for large T the distinguishability of two distributions is determined by the generalized susceptibility, as stated in the next proposition.
Proposition 5.2. Given two models (φ, g) and (φ, g ), in the limit of large T they are indistinguishable if (5.14) Remark 5.1. Although the notion of indistinguishability inherits asymmetry in g and g from the Kullback-Leibler divergence D KL (p||p ), above proposition shows that for large T the definition symmetrizes.
Remark 5.2. This proposition clarifies the role of the Fisher metric (5.1): it shows that the distance among two close-by distributions is proportional to the log-probability that the maximum likelihood estimator of a statistical model (φ, g) takes value g = g .
From this perspective, it is non-trivial to notice that this result is invariant after reparametrization of the probability densities.
This last property identifies an approximatively elliptical region of indistinguishability in the space M(φ) around each statistical model (φ, g), whose volume V T, (g) can be easily calculated in the large T limit, and is given by 2π T Finally, the volume V T, (g) allows to define the concept of density of models, and to link it to the metrics described in section 5.1.1.
Then for any fixed T and > 0 we define the density of states ρ T, (g) as (5.20)

Complexity measures and criticality
One of the most relevant problems in the field of statistical learning is the one of choosing the most appropriate model in order to fit an empirical datasetŝ generated by an unknown distribution. In particular it is well-known that models containing a large number of parameters typically lead to large values for the likelihood function P T (ŝ|g), while parsimonious models tend to produce worst in-sample values. Conversely, parsimonious models tend to generalize better, while complex models tend to fit noisy components of data leading to a poor out-of-sample performance. Using a prior function P 0 (φ, g) which keeps into account the complexity of the model itself is a practical strategy which can be used to find an optimal compromise between faithfulness to the data and generalizability of the model. Popular priors used to achieve those goals are: We will show in the following that the above results of information geometry allow to construct a measure of complexity which generalizes the BIC stated above, retaining the main feature of being completely invariant under reparametrization of the model [61,13]. In order to do this, we consider the prior: where the term N (φ) is the volume of M(φ) defined in (5.10).
Proposition 5.3. Consider the probability for an unknown datasetŝ of length T to belong to a given class of statistical models φ. Under the prior (5.25) this is given by (5. 26) In the limit T → ∞, this quantity concentrates according to: where g is the maximum likelihood estimator of g.
The proof of this result is completely analogous to the one shown in appendix A.6, and is obtained through a saddle-point expansion of the likelihood function P T (ŝ|g).
This result implies that the score assigned to the model φ converges to (up to an Our interest lies in the fact that, assuming that on the basis of dimensional analysis the complexity measure log N (φ) scales like where is a characteristic length scale, when high-dimensional models are considered, the scaling of the complexity might be anomalous, in the sense that can scale in the limit N → ∞ as a power of N . This argument additionally suggests that models φ containing critical points should be penalized by the prior (5.25), which assigns low scores to complex models. Intuitively, it has to become very costly to describe critical points even if the number of parameters of the model M is not large.
More specifically, if one assumes the scaling ∼ N α , then it is where one has the scaling H 0 (φ, g |ŝ) ∼ T . Then one can intuitively expect that a fixed scaling of T, N and M is required in order for the inverse problem to be meaningful (i.e., the left term side of (5.30) to dominate the score). Hence, when dealing with high-dimensional inference, avoiding overfitting requires not only to study how M scales with N , but also to consider that the geometric properties of the model themselves can play a role through the logarithmic correction in the last term of (5.30).

Examples
The independent spin case Consider the independent spin model described in section 2.2.5. By using equation (2.45) it is possible to find that Hence, the number of distinguishable independent spin models which can be described in an experiment of final length T with accuracy is (5.32) so that for example, just N T, ≈ 5 distinguishable models can be described by means of T = 100 observations of N = 1 spin with an accuracy of e − = 1%, while for T = 1000 and e − = 10% one gets N T, ≈ 23. The finiteness of N T, also implies that infinite regions of the h ∈ R N space belong to the same distinguishable distribution.
This can easily be checked, and one can see for example that for N = 1 the condition implies that for T = 100 and e − = 1% all models with h smaller than h min ≈ −1.16 (or h larger than h max = −h min ) belong to the same region of indistinguishability.

Fully connected ferromagnet
Let's consider the fully connected ferromagnet described in section 3.3. In that case the calculation of detχ is non-trivial, and requires an analysis of the finite N correc-

Inference of a non-equilibrium model
Many recent works in the field of neurobiology focus on neuronal ensembles which are described by means of strings of binary variables encoding the activity pattern of a set of N ∼ 10 1 or N ∼ 10 2 neurons [72,76,24]. Such compact description of the fundamental units of those system has been argued to be meaningful, triggering the expectation that techniques such as the ones described in chapter 2 might be applied on empirical data in order to extract relevant information about the interaction patterns of networks of real neurons. As a result of those expectations, striking features of neural ensembles started to emerge from the solution of the inverse problems applied to experimental data [82,58,78]. These findings posed a challenging question, whose answer has yet to be fully clarified in order to assess their validity, namely: how much of those emerging features depends on the inference procedure which has been applied, and how much is intrinsically associated with structural properties of the system? The implications of the answer go well beyond the field of neurobiology, and apply more generally to the field of statistical learning. In this section we want to provide a partial answer to this point, and show that procedures similar to the ones used to study such neural networks may generate spurious features in the inferred models, as well as genuine ones. We address in particular the issue of criticality, which we identify from the point of view of statistical mechanics with the presence of long-range correlations in a system as a result of strong collective interactions among its constituents. We apply those ideas to two datasets whose nature is similar to the one considered in [72,76,24]: a set of simulated realizations of a Hawkes point-process [38,37] and a dataset describing transactions in a financial market.

The Hawkes process
We will introduce the Hawkes point-process as a null-model to describe a system consisting of N interacting units which are able produce events in time and crossinfluence each other in absence of remarkable collective behaviors (i.e., the emergence of long-range correlations in time or space). The study of the discretized version of this model will allow an analysis of the genuine and the spurious features of the inferred model under the procedure described in chapter 2.

Definition and basic properties
We will briefly remind the notion of point process, which we will use to construct the Hawkes process, addressing the reader to [14,19] for a more detailed description.
Intuitively, X i (τ ) counts the number of events of type i falling in the interval [0, τ ], being the probability of the occurrence of an event in (τ, τ + dτ ] equal to λ i (τ )dτ , and being the one associated with the outcome of two events of order smaller than dτ . A well-known example is provided by the Poisson process, which is a pointprocesses specified by a constant, deterministic value for the intensity λ i (τ ) = µ i .
Definition 5.8. We will call a Hawkes point-process the stationary, N -variate pointprocess X(τ ) = (X 1 (τ ), . . . X N (τ )) defined by a stochastic intensity vector λ(τ ) = (λ 1 (τ ), . . . , λ N (τ )) of the form such thatK(τ ) is a positive matrix kernel satisfying We focus on the properties of this model in the stationary regime, which is guaranteed to exist for the choice of the spectral radius of the kernelK(τ ) that we specified through (5.42). Despite the fact that both averages and two point correlations of X(τ ) can be analytically computed for a large class of functionsK(τ ) [38,37], in the following discussion we will just require the knowledge of the average intensity Proposition 5.4. Given a stationary Hawkes point-process, the average intensity vector λ is given by as one can easily see by taking the expectation value of equation (5.39) and imposing the stationarity condition λ(τ ) = λ.
We employ the notation E [. . .] to indicate an average taken in the stationary state of the model, andδ denotes the identity matrix in dimension N .
We want to highlight some of the features of the Hawkes process which differentiate it from statistical models such as the ones described in section 2.1.
• Dynamics: The Hawkes process describes a stochastic process characterized by the dynamics (5.39), while a statistical model (φ, g) of the form (2.1) describes a stationary probability density. This implies that any information concerning the directionality in time (e.g., causality) of the interactions is lost when passing to a description in terms of i.i.d. binary strings. 2 • Non-stationarity: For any non-stationary generalization of the Hawkes process in which the kernel changes in time (i.e., it is of the formK(τ, τ )), or the exogenous intensity is a function µ(τ ), it is likely that inferring a stationary model may lead to errors in the interpretation of the results. In particular, what is described as an interacting, stationary system in the language of the inferred model p(s) ∝ exp J i<j s i s j + h i s i may correspond to a non-interacting, non-stationary real system [84].
• Criticality: The divergence of the mean intensity E [λ(τ )] doesn't indicate criticality of the statistical model describing the stationary state of the Hawkes process. In particular, the divergence of λ(τ ) is not linked to collective effects, as it is present even for finite N , while a proper phase transition in the statistical mechanics sense can arise just in the large N limit.
These considerations also apply when considering the binary encoding of the stochastic process describing spiking neurons, or more generally when considering any pointprocess which is binned and discretized in order to perform an inference procedure such as the one described in chapter 2.

The fully-connected Hawkes process
We introduce here the notion of fully-connected Hawkes process, which we will relate to the fully-connected pairwise model in the following part of the discussion.
Definition 5.9. Consider an N -dimensional Hawkes point-process, whose intensity vector λ(τ ) is defined by which corresponds to the choice of an exponentially decaying influence kernel K ij (τ − τ ) = α ij e −β(τ −τ ) θ(τ − τ ). Letα be a matrix of the form 45) and the vector µ to be equal to µ i = µ for each i. Then such process will be called a fully-connected Hawkes process.
For a fully-connected Hawkes process, it is easy to see by employing formula (5. 43) that while the stationarity condition (5.42) reduces to α < β.

Binning and discretization
In order to establish a connection between a spin system and a Hawkes process, we consider a discretization in time and a binarization of the signal dealt according to the following procedure. In order to shorten the notation, we will often write b i (X, δτ ) and s

Trades in a financial market
Financial markets are complex systems in which a large number of individuals interacts by buying and selling contracts at variable prices according to an unknown, dynamically varying set of criteria (e.g., their specific needs, their past experience, their future expectations). In this sense, markets can be seen as intermediary entities implicitly defined by a set of trading rules which mediate the interactions of individuals. Those rules should be such that efficient allocation of resources is achieved, so that price of traded goods reflects correct information about their fundamental value [33]. Evidence that this is not always the case has dramatically emerged in recent times [16,46,83]. Part of the responsibility has been attributed to the instability of the microscopic mechanism by means of which financial markets process information, producing prices and providing liquidity for investors [17]. Hence, it makes sense to characterize empirically how such mechanism operates, and to identify its weaknesses, its sources of inefficiencies and potential causes of its instability. With this ideas in mind, we want to characterize from the empirical point of view a part of the complex process leading to price formation.

Types of market data
In modern financial markets the action of the participants is constantly recorded, and most of the events taking place during its activity are electronically stored. In some cases, part of this data is available for investigation. In particular some main categories of datasets describing market activity which can be identified and classified on the basis of the timescale they are associated with. The most detailed level of description (timescales ranging from tens of milliseconds to the second) is achieved when informations about single market events triggered by individual agents are available [48,59,29]. A more coarse-grained description of the market is obtained by focusing on the price process and its variations. More precisely, it is possible to define an instantaneous price for any contract, and to keep track of all its variations (tick-bytick data) throughout the duration of the market activity. Data describing all events changing the price (called either trade or quote events) are necessary to achieve this level of description (being the typical time resolution required the one of the second) [28]. Finally, data corresponding to market behavior at lower frequencies are often publicly available, and involve, beyond the daily opening and closing price, the volume traded and the highest and lowest daily price for all traded goods (e.g., they can be found in [5]). In this discussion, we focus on data describing trade events, which belong to the intermediate regime in which the price process is monitored with the resolution of around one second. Any of those trade events corresponds to the transfer of a contract from a seller to a buyer at a given price, for a given quantity (volume) of a good.

Cross-correlation of trade events
It has been observed in empirical data across several markets that trade events of single securities are not independent one from another, rather they influence each other leading to interesting clustering phenomena. Moreover, by considering multiple securities traded in the same market venue, it is possible to check that even event times associated with the to trade of different instruments are strongly correlated among each other. Then, one can be interested in answering the following question: do correlations in trading times arise from correlated exogenous phenomena driving market activity, or do they form due to an endogenous contagion process spreading across the market? While the former scenario would correspond to a picture in which market activity reflects fundamental exogenous information, the latter would be associated with the scenario of a (potentially unstable) market which self-interacts without necessarily assimilating external information. Those scenarios can in principle coexist, although it is not easy to construct a quantitative, empirically measurable notion dis-tinguishing the two regimes [44]. It should also be added that part of the explanation for long-range correlation in trading times has been identified in the mechanism of order-splitting: the finite amount of liquidity available in the market forces traders to split large orders (meta-orders, care-orders or hidden-orders) in smaller lots which are traded incrementally, leading to long-range correlation of trading times (from hours to days, sometimes even up to weeks). Indeed a relevant role could also be played by collective interactions across different securities, which could lead to correlated order flow. This possibility is empirically inquired in section (5.3.2), where we apply the techniques described in chapter 2 to this type of financial system, and try to understand the results on the basis of what we presented in the early part of this chapter.

Applications
With these ideas in mind, we consider two sets of realizations of a point-processes X(τ ).  • Correlations drop to zero as the window δτ is made smaller (a phenomenon known in the field of finance as Epps effect [32]

Features of the inferred models
Extracting the couplings of a fully connected model from the values of magnetization and correlation described above leads to the results depicted in figure 5.6, where we consider both the disordered case g = (h,Ĵ ) and the two-parameter model g = (h, J).
We stress in the following some of the main features.   (3.11)) and TAP equations (equations (3.14) and (3.15)) -which produced consistent results. We stress some features of the results we obtained:

Pairwise fully-connected model for NYSE trade events
• The ratio h/J changes according to δt, so that it is not possible to interpret h as measuring exogenous driving factors and J as a genuine interaction. Moreover, as explained in section 3.3), this inference procedure may mix interactions with external fields due to the approximate symmetry g → g + δJ(−m, 1). What is possible to say is that the Hawkes process which best describes this (h, J) curve is defined by parameters µ ≈ 0.011 s −1 , α ≈ 0.022 s −1 , β ≈ 0.03 s −1 , so that the exogenous intensity µ corresponds to approximatively one fourth of the average intensity E[λ(τ )].
• Results describing the Hawkes process allow to understand the proximity of the inferred parameters (h, J) to the critical point as related to the divergence of the average intensity E[λ(τ )], rather than arising from a collective effect.
• As in the previous case, the inferred model doesn't lie on the critical line h = 0.
This is due to the fact (section 3.3) that correlations are of the order of 1/N , so that a description in terms of fully-connected ferromagnet leads to a nondegenerate description of the data.
Remark 5.4. A procedure which has been proposed to estimate the distance of an inferred model from the critical point consists in rescaling of all the couplings by a common factor β (i.e., performing a shift g → βg ) which is interpreted as a fictitious inverse temperature. Studying how the elements of the susceptibility matrix χ vary with respect to β should allow to identify criticality in the inferred model by the presence of peaks close to β = 1 in specific components of the matrix. We perform this procedure with our data and plot the results in figure 5.11, finding that: • This procedure is not isotropic, in the sense that the shift g → βg implicitly indicates that the direction (1, . . . , 1) should be the preferred one in order to evaluate distances in the coupling space.
• This type of measure does not describe the distance of the inferred model from the critical point in term of distinguishable models (equivalently, this measure of distance is not invariant under reparametrization of the statistical model (φ, g)).
These points can lead to problems when model condensation is present: very different models may be described by slightly shifting β. Moreover an inferred model may lie close to the critical point due not only due to model condensation, but also due to the choice of a stable inference procedure, so that it is likely thatχ attains large value in the point g , and that by moving from that point one can expect fluctuations to strongly decrease. A better measure of distance would be provided by considering geodesics in the coupling space under the Fisher metricsχ, as shown in section 5.1.2. We remind that properties (2.36) and (2.37) allows to informally identify this measure as counting how many error bars is one away from the critical point. This approach doesn't specify any privileged direction in the coupling space (as the geodesic distance is associated with whatever path in the coupling space is minimizing the number of such error bars), nor varies according to the reparametrization (as in that case also error bars are reparametrized accordingly). As an example, one finds that the distance d T, (p crit , p )

Chapter 6 Conclusion
In this work we have presented a general approach to the field of statistical learning, in which the problem of estimating parameters describing a complex system is seen as an inverse problem in the field of statistical mechanics. This perspective has been proven to be especially relevant in order to study extended systems, which in this language are associated to physical systems in the thermodynamic limit. This regime is wellknown in physics, and several techniques (mean-field approximations) are known to solve the inverse problem in this framework with high accuracy. Interesting collective features emerge in this regime: models can condensate leading to regions of the space of parameters which are able to describe anomalously well very diverse datasets, and null-modes can develop leading to degenerate representations of a dataset. All of these factors have to be kept into account when studying inverse problem for empirical datasets, in order to disentangle the genuine features of a system from the spurious ones depending on the inference procedure which is applied. We have also shown that complete representations of the inverse problem lead to the exact solution of several systems (complete systems, one-dimensional systems, tree-like interaction networks), and allow a general understanding of the locality and stability features of the inverse problems, which are easier and more resilient to noise than the direct CHAPTER 6. CONCLUSION ones. The notion of regularizer has also been discussed and its use has been clarified by specific, solvable examples, in order to show general features of non-parametric inference. We find that a symmetry property characterizes the regularizers, and a tradeoff between computational complexity and relevance of the inference procedure has to be sought on the basis of such symmetry in order to perform model selection.
Finally, we have shown how differential geometry can be used to understand the consistency of the inverse problem, and how the special features associated with the large N limit have a clear geometric interpretation in terms of distance and volume. In this language, criticality of the inferred model is related to the strong divergence of the number of datasets which can be described through a small shift of the inferred parameters. Finally, we have presented the application of these ideas to two datasets, a synthetic one describing a self-excitatory point process and an empirical one describing transactions in a financial market. We used those datasets in order to illustrate our ideas by separating genuine features of the inferred model and spurious ones, finding that the dataset describing financial transactions can be welldescribed by a fully-connected ferromagnet in which interactions play a prominent role with respect to external driving factors.
according to the specific goal of the inference problem which one is trying to solve (e.g., which observables are considered relevant for a particular application).
Proof. The proof of this result amounts to solve the constrained optimization problem Proof. First one can show that which allows to proof that for any vector x, the quadratic form x Tχ x is greater or equal than zero. In fact one has that µ,ν>0 Additionally, if the operators φ µ,s are minimal in the sense defined in section 2.1 above expression has to be strictly larger than zero for x = 0. In fact if x Tχ x = 0, then it must hold for each state s that which by minimality of φ implies that x µ = 0 for each µ.

A.3 Small deviations of the empirical averages
We want to show that given a statistical model (φ, g), equations (2.23) and (2.24) hold for the empirical averagesφ.
Proof. For the averages, it is sufficient to show that due to the factorization property of P T (ŝ|g) one has while for the covariances one can write By noting that due to independence all terms with t = t vanish from previous expression, one recovers equation (2.24).

A.4 Sanov theorem
We want to prove Sanov theorem (2.35), which states that given a probability distri-

A.5 Cramér-Rao bound
Cramér-Rao bound states that given a statistical model (φ, g) with F (g) strictly convex and an unbiased estimator of g denoted as g , the covariance matrix of g under the measure . . . T is bound according to equation (2.38).
Proof. First, it is necessary to prove that to prove that, after defining V µ = ∂ log P T (ŝ|g) ∂gµ , and using equation (2.15) one has where we also used equation (2.23) (i.e., φ µ = φ µ T ). Then, it is possible to compute the covariance and exploit Cauchy-Schwartz inequality, which implies that for any pair of vectors x, y it holds that

Equation (A.17) fixes the value of the left-hand side term of equation (A.19) to be
x T y, while the right-hand side can be expanded into where we used that V V T T = Tχ due to equation (2.24). Finally, by choosing the arbitrary vector x to be x =χ −1 y/T (χ is invertible due to strict concavity of F (g)), it holds for any y that which proves the thesis (2.38).

A.6 Convergence of the inferred couplings
Given a set of empirical frequenciesp, we want to prove that for a generic set of models described by an operator set φ the mean and the covariances of couplings g defining a probability distribution p, weighted by the measure provided by the posterior P T (g|p) are given in the limit of large T by equations (2.36) and (2.37).
Proof. To calculate them, we first notice that, by defining it is possible to write so that the calculation of the generating function log Z(φ) allows to find the required momenta of g. In the limit of large T it is possible to perform a saddle-point estimation of the function Z(φ) around the minimum of the convex function D KL (p||p) (or, equivalently, the maximum of the concave free energy F (g)), which requires the expansion of the Kullback-Leibler divergence. This procedure yields where -as shown in section 2.2 -the maximum likelihood estimator g can be defined as the minimizer of the Kullback Leibler divergence. The differentiation of log Z(φ) finally leads to Proof. This can be shown by noting that by using Stirling formula and approximating the sum with an integral one can write has two minima m + and m − , whose contribution can be kept into account through The contribution of F 0 (h, J) to the ensemble averages is while the one to the susceptibility matrixχ is given bŷ where χ s.p. = ∂m s.p. /∂h. Its eigenvalues are given by N χ s.p. (0, 1 + m 2 s.p. ).

The role of Gaussian fluctuations
The term F f luct (h, J) allows to compute the eigenvalue decomposition for the matrix χ, whose smallest eigenvalue receives a contribution which grows in N , and is related to the Gaussian integral (B.1). It results .

(B.11)
The contribution of F f luct (h, J) to the solution of the direct problem is The matrixb(h, J, m) (whose explicit form is not particularly illuminating) can be obtained by deriving above averages with respect to h and J.

Determinant of the generalized susceptibility
The term √ detχ is shown in chapter 5 to be relevant in order to count the number of where, after easy combinatorics, we used the fact that the correlation c (t) measured in the observation number t depends just upon the total magnetization m (t) . Finding a solution to this problem is easy due to convexity of t (m (t) ) 2 . In particular by taking the limit T → ∞ a solution can be found for any m, while then the minimum and the maximum value of c are given respectively by m 2 −1/N 1−1/N and 1. Interestingly, the same result can be obtained with more simplicity by exploiting the necessary condition Var[ i s i ] ≥ 0. Notice also that for large N the connected correlation coefficient c − m 2 is bound from below by (m 2 − 1)/N : equivalently no large system, subject to whatever type of interaction, can be globally anti-correlated.
3. The Hessian defines the descent direction v = −χ −1 (g)∇H(g). Algorithms exploiting this property usually go under the name of Newton's methods.
These properties are simples consequences of differentiability and convexity of H(g), and allow to solve the problem the problem of its minimization. The first property implies that given a g such that ∇H(g) = 0, g is a global minimum of H(g). If this equation can be explicitly solved, the minimum can be found. Indeed if, as it often is the case, the condition ∇H(g) = 0 is non-analytically solvable, it is possible to exploit properties 2. and 3. in order to build iterative algorithms which decrease the target function H(g) at each step. In particular, iterative algorithms exploiting property 2. are expected to achieve linear convergence to the minimum, while more sophisticated algorithms (Newton methods) constructed by using the Hessian can achieve quadratic convergence. More efficient schemes (quasi-Newton methods) such as the L-BFGS approximation [49,23] exploit an approximation for the Hessian in order to save memory and computational power by exploiting successive updates of the gradient. We present in the following an example of a simple algorithm which can be used to minimize a convex differentiable H(g), which we use mainly as a proof of principle for the solvability of this type of problem. Secondly, the efficiency of the Boltzmann learning algorithm presented in section 3.1 is rooted in the gradient descent method.

C.1.1 Gradient descent algorithm
Given a convex, differentiable H(g) and starting point g (0) , we consider a sequence {g (k) } K k=1 built according to the iterative scheme where we introduced the schedule { k } K k=1 . Suppose that each of the k is chosen in order to satisfy the (Armijo) condition for a given 0 < β < 1, by considering the initial value k = 1 and iterating the map k ← k /τ for τ < 1 until (C.4) is satisfied. Then it holds that either min g∈R M H(g) = −∞ or lim k→∞ ||∇H(g)|| 2 = 0, that is, if a minimum exists, the sequence {g (k) } K k=1 can approximate it with arbitrary precision.
Remark C.1. Searching the optimal k is usually called a line search, and the procedure that we introduce to find it is guaranteed to find an k satisfying (C.4) if the maximum eigenvalue ofχ is bounded by a given χ max . In particular the convexity of H(g) and a straightforward application of Taylor theorem allow to prove that any k in the interval 0 ≤ k ≤ 2(1 − β) χ max (C.5) satisfies the Armijo condition (C.4).
Such limit can be either −∞ (in which case H(g) has no minimum) or can be finite.
The finiteness of the limit implies that H(g (∞) ) − H(g (0) ) = lim K→∞ K k=0 k ||∇H(g (k) )|| 2 (C.7) which leads to lim K→∞ K ||∇H(g (K) )|| 2 = 0 . (C.8) Notice that the rate of convergence in K of this algorithm can be rather slow, which is the reason why more sophisticated algorithms are commonly used to perform this task (see [21]).

C.2 Non-differentiable target
If a convex function is not differentiable in all of its domain the solution of the minimization problem is technically more complicated, but it is still possible to take advantage of the convexity property in order to build efficient minimization algorithms (see [20,21]). Consider a convex function H(g) : R M → R. Then one can define a sub-gradient as any global underestimator of H(g), namely v ∈ R M is a subgradient of H(g) in g if for any g one has H(g) ≥ H(g ) + v T (g − g ) . (C.9) The set of all the sub-gradients of H(g) in g is called the sub-differential of H(g), and is denoted with∇H(g). One can show that •∇H(g) is non-empty if H(g) is locally convex and bounded around g.
•∇H(g) is closed and convex.
This properties characterize the sub-differential as a notion generalizing the ordinary differential, which is suitable to solve problems involving non-differentiable functions.
In particular the properties shown above for differentiable functions generalize to: 1. If 0 ∈∇H(g) then g is a global minimum of H(g).
2. The direction v = − ∇ H(g) is not in general a descent direction.
This implies that in order to minimize a non-differentiable function it is still possible to find the points whose sub-differential is equal to zero, but that a naive sub-gradients descent similar to (C.3) is not guaranteed to find a solution.
An example: the absolute value which is minimum for The notion of sub-gradient also allow us to generalize the gradient descent algorithm to non-differentiable functions, as shown in the following.
Notice that in this case the sequence k is not optimized on-line, rather it is fixed at the beginning of the algorithm. This is because the sub-gradient doesn't specify necessarily a descent direction, hence the sub-gradient descent may increase the function H(g), requiring the values g best and H(g best ) to be stored at each iteration step.

Appendix D
Complete families D.1 Rate of convergence for the complete inverse problem Consider a statistical model p in which all states have strictly positive probability (i.e. it exists a p min = 0 such that ∀s p min ≤ p s ). We want to show how within inference scheme (4.14) the inferred couplings concentrate around their actual values at fixed N in the limit T → ∞.
Proof. The expression for g µ is:

D.2 Factorization property for tree-like models
In this section we prove a fundamental property of statistical models whose the interaction structure is loop-less, which we call trees and rigorously define as follows. 1 Definition D.1. Consider a statistical model (φ, g) of the form (2.1), with g µ = 0 for all g µ ∈ g. Then the set φ is called a tree if it is not possible to find a cycle connecting any set of vertices, i.e., it doesn't exist a closed path {i 1 , . . . , i L−1 , i L = i 1 } ∈ V L such that for each couple {i n , i n+1 } there exist an operator φ in,i n+1 ∈ φ depending on both s in and s i n+1 , with φ in,i n+1 = φ im,i m+1 for all n = m ∈ {1, . . . , L − 1}.
For trees we will show along the lines of [54] that the following factorization property holds. Then it is possible to reabsorb the ψ(s i ) factor inside a new operator obtained by the following change on a generic φ ρ ∈ ∂i\φ µ : φ ρ (s ∂φρ ) → φ ρ (s ∂φρ ) = φ ρ (s ∂φρ ) + log ψ(s i ) g ρ . (D.10) The statistical model describing the reduced problem for Γ spins can thus be described by using M operators, so that it is possible to use the inductive hypothesis to show where only the term −S(q) depends on the distribution q. Then the maximum entropy principle (appendix A.1) states that the density q ∈ M(φ, g ) maximizing S(q) is the statistical model described (φ, g ), whose associated probability density has been called p . Then one has finally, by using property 2. of M(φ, g ) and the continuity property 3. one has that for δ sufficiently small, Prob[p ∈ M (φ, g )] is arbitrarily close to Prob[g (φ) − g ∈ δg], which together with (E.11) proves the thesis (5.11).

E.3 Expansion of the Kullback-Leibler divergence
We want to prove that, given a pair of statistical models (φ, g) and (φ, g ) and an accuracy parameter , for large T they are indistinguishable if condition (5.14) holds.

E.4 Volume of indistinguishability
Given the space M(φ) identified by the minimal operator set φ, we want to show that the volume of the space of indistinguishable distributions around a point g is given by equation (5.15), where T is the length of the dataset and > 0 is the accuracy parameter.

E.5 Estimation of the empirical observables for an
Hawkes point process Consider a fully connected Hawkes process defined as in (5.39), characterized by exogenous intensity µ and kernel parameters parameters α and β. We will show that the qualitative features of the fully connected pairwise model associated through the binning functions (5.47) and (5.48) can be obtained by using an approximate scheme.
More precisely, given a realization of a fully connected Hawkes point-process X and a bin size δτ we will calculate the quantities i (X, δτ ) s  where δX c (t δτ ) = X c (δτ (t + 1)) − X c (δτ t), while c ∈ {i, j, i + j} and \c ∈ {V \{i}, V \{j}, V \{i + j}} refer to the channels which one needs to take into account to calculate magnetizations and correlations. Above probability can be computed by taking into account that: • The convolution of a set of Hawkes processes is a Hawkes process. λ \c (τ ) = µ \c + e −βτ (λ \c − µ \c ) + K k=1 α \c,\c e −β(τ −τ k ) θ(τ − τ k ) (E. 28) and perform explicitly the integral. This is very hard to do analytically, so that we consider an approximate scheme in which it is possible to obtain a qualitatively correct result for the averages, motivated by the fact that in both the cases that we consider (c = i and c = i + j) we have that α \c,c , α \c,c α c,c , α c,\c and µ \c µ c . This regime justifies the approximation in which the trajectory λ \c (τ ) is described by the is the average in the stationary state in which channel c is conditioned in order not to produce events. Finally ψ(δτ ) is a generic function such that φ(0) = 1 and φ(∞) = 0.
Then, one can insert this approximation into equation (E. 26), supposing that the number of events K is deterministic and concentrated around its average number.