On Maximum Entropy and Inference

Maximum Entropy is a powerful concept that entails a sharp separation between relevant and irrelevant variables. It is typically invoked in inference, once an assumption is made on what the relevant variables are, in order to estimate a model from data, that affords predictions on all other (dependent) variables. Conversely, maximum entropy can be invoked to retrieve the relevant variables (sufficient statistics) directly from the data, once a model is identified by Bayesian model selection. We explore this approach in the case of spin models with interactions of arbitrary order, and we discuss how relevant interactions can be inferred. In this perspective, the dimensionality of the inference problem is not set by the number of parameters in the model, but by the frequency distribution of the data. We illustrate the method showing its ability to recover the correct model in a few prototype cases and discuss its application on a real dataset.


Introduction
Statistical mechanics stems from classical (or quantum) mechanics. The latter prescribes which are the relevant quantities (i.e., the conserved ones). The former brings this further, and it predicts that the probability to observe a system in a microscopic state s, in thermal equilibrium, is given by: where H[s] is the energy of configuration s. The inverse temperature β is the only relevant parameter that needs to be adjusted, so that the ensemble average H matches the observed energy U. It has been argued [1] that the recipe that leads from H[s] to the distribution P(s) is maximum entropy: among all distributions that satisfy H = U, the one maximizing the entropy S = − ∑ s P(s) log P(s) should be chosen. Information theory clarifies that the distribution Equation (1) is the one that assumes nothing else but H = U, or equivalently, that all other observables can be predicted from the knowledge of H[s]. This idea carries through more generally to inference problems: given a dataset of N observationŝ s = {s (1) , . . . , s (N) } of a system, one may invoke maximum entropy to infer the underlying distribution P(s) that reproduces the empirical averages of a set M of observables φ µ (s) (µ ∈ M). This leads to Equation (1) with: − βH[s] = ∑ µ∈M g µ φ µ (s) (2) that result from entropy maximization and are also known to coincide with maximum likelihood estimation (see [2][3][4]). For example, in the case of spin variables s ∈ {±1} n , the distribution that reproduces empirical averages s i and correlations s i s j is the pairwise model: which in the case J ij ≡ J ∀i, j and h i ≡ h ∀i is the celebrated Ising model. The literature on inference of Ising models, stemming from the original paper on Boltzmann learning [5] to early applications to neural data [6] has grown considerably (see [7] for a recent review), to the point that some suggested [8] that a purely data-based statistical mechanics is possible.
Research has mostly focused on the estimate of the parameters g = {h i , J ij }, which itself is a computationally challenging issue when n 1 [7], or in recovering sparse models, distinguishing true interactions (J ij = 0) from spurious ones (J ij = 0; see, e.g., [9]). Little has been done to go beyond pairwise interactions (yet, see [10][11][12][13]). This is partly because pairwise interactions offer a convenient graphical representation of statistical dependences; partly because -th order interactions require ∼n parameters and the available data hardly ever allow one to go beyond = 2 [14]. Yet, strictly speaking, there may be no reason to believe that interactions among variables are only pairwise. The choice of form (4) for the Hamiltonian represents an assumption on the intrinsic laws of motion, which reflects an a priori belief of the observer on the system. Conversely, one would like to have inference schemes that certify that pairwise interactions are really the relevant ones, i.e., those that need to be included in H in order to reproduce correlations s i 1 s i 2 · · · s i of arbitrary order [15].
We contrast a view of inference as parameter estimation of a preassigned (pairwise) model, where maximum entropy serves merely an ancillary purpose, with the one where the ultimate goal of statistical inference is precisely to identify the minimal set M of sufficient statistics that, for a given datasetŝ, accurately reproduces all empirical averages. In this latter perspective, maximum entropy plays a key role in that it affords a sharp distinction between relevant variables (φ µ (s), µ ∈ M) (which are the sufficient statistics) and irrelevant ones, i.e., all other operators that are not a linear combination of the relevant ones, but whose values can be predicted through theirs. To some extent, understanding amounts precisely to distinguishing the relevant variables from the "dependent" ones: those whose values can be predicted.
Bayesian model selection provides a general recipe for identifying the best model M; yet, as we shall see, the procedure is computationally unfeasible for spin models with interactions of arbitrary order, even for moderate dimensions (n = 5). Our strategy will then be to perform model selection within the class of mixture models, where it is straightforward [16], and then, to project the result on spin models. The most likely models in this setting are those that enforce a symmetry among configurations that occur with the same frequency in the dataset. These symmetries entail a decomposition of the log-likelihood with a flavor that is similar to principal component analysis (yet of a different kind than that described in [17]). This directly predicts the sufficient statistics ψ λ (s) that need to be considered in maximum entropy inference. Interestingly, we find that the number of sufficient statistics depends on the frequency distribution of observations in the data. This implies that the dimensionality of the inference problem is not determined by the number of parameters in the model, but rather by the richness of the data.
The resulting model features interactions of arbitrary order, in general, but is able to recover sparse models in simple cases. An application to real data shows that the proposed approach is able to spot the prevalence of two-body interactions, while suggesting that some specific higher order terms may also be important.

Spin Models with Interactions of Arbitrary Order
Consider a system of n spin variables s i = ±1, a state of which is defined by a configuration s = (s 1 , . . . , s n ). The number of all possible states is 2 n . A generic model is written as: where φ µ (s) = ∏ i∈µ s i is the product of all spins involved in the corresponding interaction g µ , the sum on µ runs on a subset of operators φ µ and Z ensures normalization. We follow the same notation as in [18]: There are 2 n possible such operators, which can be indexed by an integer µ = 0, . . . , 2 n − 1 whose binary representation indicates those spins that occur in the operator φ µ . Therefore, µ = 0 corresponds to the constant operator φ 0 (s) = 1, and for µ = 11, the notation i ∈ µ is equivalent to i ∈ {1, 2, 4}, i.e., φ 11 (s) = s 1 s 2 s 4 . Given a datasetŝ of N observations of s, assuming them to be i.i.d. draws from P(s|g, M), the parameters g µ are determined by solving Equation (3). Bayesian inference maintains that different models should be compared on the basis of the posterior P{M|ŝ} that can be computed integrating the likelihood over the parameters g (we refer to [18] for a discussion of Bayesian model selection within this setup). This can already be a daunting task if n is large. Note that each operator φ µ (µ > 0) can either be present or not in M; this implies that the number of possible models is 2 2 n −1 . Finding the most likely model is impossible, in practice, even for moderately large n.

Bayesian Model Selection on Mixture Models
Let us consider mixture models in which the probability of state s is of the form: where 1 A (s) is the indicator function 1 A (s) = 1 if s ∈ A, and 1 A (s) = 0 otherwise. The prior (and posterior) distributions of take a Dirichlet form; see Appendix B. In other words, model Q assigns the same probability ρ j to all configurations in the same set Q j . Formally, Q = {Q i } is a partition of the set of configurations s, i.e., a collection of subsets such that {±1} n = j Q j and Q j Q j = ∅ ∀j = j . The model's parameters ρ j are subject to the normalization constraint ∑ j |Q j |ρ j = 1, where |Q| stands for the number of elements within Q. We denote by q = |Q| the number of subsets in Q. The number of independent parameters in model Q is then q − 1.
The number of possible models Q is the number of partitions of a set of 2 n elements, which is the Bell number B 2 n . This grows even faster than the number of spin models M. Yet, Bayesian model selection can be easily carried out, as shown in [16], assuming Dirichlet's prior. In brief, the most likely partition Q * depends on the assumed prior, but it is such that if two states s and s are observed a similar number of times k s k s , then the most likely model places them in the same set Q q and assigns them the same probability ρ q . In other words, considering the frequency partition: that groups in the same subset all states s that are observed the same number k s of times, the optimal partition Q * is always a coarse graining of K, likely to merge together subsets corresponding to similar empirical frequencies.
We refer the interested reader to Appendix B and [16] for more details, as well as for a heuristic for finding the Q * model.

Mapping Mixture Models into Spin Models
Model Q allows for a representation in terms of the variables g µ , thanks to the relation: which is of the same nature of the one discussed in [11] and whose proof is deferred to Appendix A. The index in g µ Q indicates that the coupling refers to model Q and merely corresponds to a change of variables → g; we shall drop it in what follows, if it causes no confusion.
In Bayesian inference, should be considered as a random variable, whose posterior distribution for a given datasetŝ can be derived (see [16] and Appendix B). Then, Equation (8) implies that also g is a random variable, whose distribution can also be derived from that of .
Notice, however, that Equation (8) spans only a q − 1-dimensional manifold in the 2 n − 1-dimensional space g, because there are only q − 1 independent variables . This fact is made more evident by the following argument: Let v be a 2 n − 1 component vector such that: Then, we find that: In other words, the linear combination of the random variables g µ with coefficients v µ that satisfy Equation (9) is not random at all. There are (generically) 2 n − 1 − q vectors v that satisfy Equation (9) each of which imposes a linear constraint of the form of Equation (10) on the possible values of g.
In addition, there are q orthogonal directions u λ that can be derived from the singular value decomposition of χ µ j : This in turn implies that model Q can be written in the exponential form (see Appendix D for details): where: The exponential form of Equation (12) identifies the variables ψ λ (s) with the sufficient statistics of the model. The maximum likelihood parametersĝ λ can be determined using the knowledge of empirical averages of ψ λ (s) alone, solving the equations ψ λ = ψ λ for all λ = 1, . . . , q. The resulting distribution is the maximum entropy distribution that reproduces the empirical averages of ψ λ (s). In this precise sense, ψ λ (s) are the relevant variables. Notice that, the variables ψ λ (s) are themselves an orthonormal set: In particular, if we focus on the K partition of the set of states, the one assigning the same probability ρ k to all states s that are observed k times, we find that P(s|ĝ, Q) = k s /N exactly reproduces the empirical distribution. This is a consequence of the fact that the variablesĝ λ that maximize the likelihood must correspond to the maximum likelihood estimatesρ k = k/N, via Equation (8). This implies that the maximum entropy distribution Equation (12) reproduces not only the empirical averages ψ λ (s), but also that of the operators φ µ (s) for all µ. A direct application of Equation (8) shows that the maximum entropy parameters are given by the formula: Similarly, the maximum likelihood parametersĝ λ are given by: Notice that, when the set Q 0 = {s : k s = 0} of states that are not observed is not empty, all couplingsĝ µ with χ µ 0 = 0 diverge. Similarly, allĝ λ with w λ,0 = 0 also diverge. We shall discuss later how to regularize these divergences that are expected to occur in the under-sampling regime (i.e., when N ≤ 2 n ).
It has to be noted that, of the q parameters ρ q , only q − 1 are independent. Indeed, we find that one of the q singular values Λ λ in Equation (11) is practically zero. It is interesting to inspect the covariance matrix C µ,ν = E[δg µ δg ν ] of the deviations δg µ = g µ − E[g µ ] from the expected values computed on the posterior distribution. We find (see Appendices C and D) that C µ,ν has eigenvalues Λ 2 λ along the eigenvectors u λ and zero eigenvalues along the directions v. The components λ with the largest singular value Λ λ are those with the largest statistical error, so one would be tempted to consider them as "sloppy" directions, as in [19]. Yet, by Equation (17), the value ofĝ λ itself is proportional to Λ λ , so the relative fluctuations are independent of Λ λ . Indeed "sloppy" modes appear in models that overfit the data, whereas in our case, model selection on mixtures ensures that the model Q * does not overfit. This is why relative errors on the parameters g λ are of comparable magnitude. Actually, variables ψ λ that correspond to the largest eigenvalues Λ λ are the most relevant ones, since they identify the directions along which the maximum likelihood distribution Equation (12) tilts most away from the unconstrained maximal entropy distribution P 0 (s) = 1/2 n . A further hint in this direction is that Equation (17) implies that variables ψ λ (s) with the largest Λ λ are those whose variation across states s typically correlates mostly with the variation of log k s in the sample.
Notice that the procedure outlined above produces a model that is sparse in the variables, i.e., it depends only on q − 1 parameters, where q is, in the case of the K partition, the number of different values that k s takes in the sample. Yet, it is not sparse in the g µ variables. Many of the results that we have derived carry through with obvious modifications if the sums over µ are restricted to a subset M of putatively relevant interactions. Alternatively, the results discussed above can be the starting point for the approximate scheme to find sparse models in the spin representation.

Illustrative Examples
In the following, we present simple examples clarifying the effects of the procedure outlined above.

Recovering the Generating Hamiltonian from Symmetries: Two and Four Spins
As a simple example, consider a system of two spins. The most general Hamiltonian that should be considered in the inference procedure is: Imagine the data are generated from the Hamiltonian: and let us assume that the number of samples is large enough, so that the optimal partition Q * groups configurations of aligned spins Q = = {s : s 1 = s 2 } distinguishing them from the configuration of unaligned ones Q = = {s : s 1 = −s 2 }. Following the strategy explained in Section 2.1, we observe that χ µ = = χ µ = = 0 for both µ = 1 and 2. Therefore, Equation (8) implies g 1 = g 2 = 0. Therefore, the Q * model only allows for g 3 to be nonzero. In this simple case, symmetries induced by the Q * model (i.e., (s 1 , s 2 ) → (−s 1 , −s 2 )) directly produce a sparse model where all interactions that are not consistent with them are set to zero.
Consider now a four-spin system. Suppose that the generating Hamiltonian is that of a pairwise fully-connected model as in Figure 3 (left), with the same couplings g 3 = g 5 = g 6 = g 9 = g 1 0 = g 1 2 = J. With enough data, we can expect that the optimal model is based on the partition Q * that distinguishes three sets of configurations: depending on the absolute value of the total magnetization. The Q * model assigns the same probability ρ j to configurations s in the same set Q j . Along similar lines to those in the previous example, it can be shown that any interaction of order one is put to zero (g 1 = g 2 = g 4 = g 8 = 0), as well as any interaction of order three (g 7 = g 11 = g 13 = g 14 = 0), because the corresponding interactions are not invariant under the symmetry s → −s that leaves Q * invariant. The interactions of order two will on the other hand correctly be nonzero and take on the same value g 3 = g 5 = g 6 = g 9 = g 1 0 = g 1 2.
The value of the four-body interaction is: This, in general, is different from zero. Indeed, a model with two-and four-body interactions shares the same partition Q * in Equation (20). Therefore, unlike the example of two spins, symmetries of the Q * model do not allow one to recover uniquely the generative model (Figure 3, left). Rather, the inferred model has a fourth order interaction (Figure 3, right) that cannot be excluded on the basis of symmetries alone. Note that there are 2 2 4 −1 = 32,768 possible models of four spins. In this case, symmetries allow us to reduce the set of possible models to just two.

Exchangeable Spin Models
Consider models where P(s) is invariant under any permutation π of the spins, i.e., P(s 1 , . . . , s n ) = P(s π 1 , . . . , s π n ). For these models, P(s) only depends on the total magnetization ∑ i s i . For example, the fully-connected Ising model: belongs to this class. It is natural to consider the partition where Q q contain all configurations with q spins s i = −1 and n − q spins s j = +1 (q = 0, 1, . . . , n). Therefore, when computing χ µ q , one has to consider |Q q | = ( n q ) configurations. If µ involves m spins, then ( m j )( n−m q−j ) of them will involve j spins s i = −1, and the operator φ µ (s) takes the value (−1) j on these configurations. Therefore, χ µ q only depends on the number m = |µ| of spins involved and: This implies that the coefficients g µ of terms that involve m spins must all be equal. Indeed, for any two operators µ = µ with |µ| = |µ |: Therefore, the proposed scheme is able, in this case, to reduce the dimensionality of the inference problem dramatically, to models where interactions g µ only depend on the number m = |µ| of spins involved in φ µ .
Note also that any non-null vector v µ such that ∑ µ v µ = 0 and v µ = 0 if |µ| = m > 0 satisfies: The vectors u λ corresponding to the non-zero singular values ofχ need to be orthogonal to each of these vectors, so they need to be constant for all µ that involve the same number of spins. In other words, u µ λ = u λ (|µ|) only depend on the number |µ| of spins involved. A suitable choice of a set of n independent eigenvectors in this space is given by u λ (m) = aδ λ,m that correspond to vectors that are constant within the sectors of µ with |µ| = λ and are zero outside. In such a case, the sufficient statistics for models of this type are: as it should indeed be. We note in passing that terms of this form have been used in [20].
Inference can also be carried out directly. We first observe that the g λ are defined up to a constant. This allows us to fix one of them arbitrarily, so we will take g 0 = 0. If K λ is the number of observed configurations with λ spins s i = −1, then the equation φ λ = φ λ (for λ > 0) reads: n λ e g λ so that, after some algebra, From this, one can go back to the couplings of operators g µ using: Figure 4 illustrates this procedure for the case of the mean field (pairwise) Ising model Equation (22). As this shows, the procedure outlined above identifies the right model when the number of samples is large enough. If N is not large enough, large deviations from theoretical results start arising in couplings of highest order, especially if β is large.

The Deep Under-Sampling Limit
The case where the number N of sampled configurations is so small that some of the configurations are never observed deserves some comments. As we have seen, taking the frequency partition K, where Q k = {s : k s = k}, if Q 0 = ∅, then divergences can manifest in those couplings where χ µ 0 = 0.
It is instructive to consider the deep under-sampling regime where the number N of visited configurations is so small that configurations are observed at most once in the sample. This occurs when N 2 n . In this case, the most likely partitions are (i) the one where all states have the same probability Q 0 and (ii) the one where states observed once have probability ρ/N and states not yet observed have probability (1 − ρ)/(2 n − N), i.e., Q 1 = {Q 0 , Q 1 } with Q k = {s : k s = k}, k = 0, 1. Following [16], it is easy to see that that generically, the probability of model Q 0 overweights the one of model Q 1 , because P{ŝ|Q 1 } P{ŝ|Q 0 }. Under Q 0 , it is easy to see that χ µ 0 = 0 for all µ > 0. This, in turn, implies that g µ = 0 exactly for all µ > 0. We reach the conclusion that no interaction can be inferred in this case [21].
Taking instead the partition Q 1 , a straightforward calculation shows that Equation (8) leads to g µ = aφ µ . Here, a should be fixed in order to solve Equation (3). It is not hard to see that this leads to a → ∞. This is necessary in order to recover empirical averages, which are computed assuming that unobserved states s ∈ Q 0 have zero probability.
This example suggests that the divergences that occur when assuming the K partition, because of unobserved states (k s = 0) can be removed by considering partitions where unobserved states are clamped together with states that are observed once.
These singularities arise because, when all the singular values are considered, the maximum entropy distribution exactly reproduces the empirical distribution. This suggests that a further method to remove these singularities is to consider only the first singular values (those with largest Λ λ ) and to neglect the others, i.e., to set g λ = 0 for all other λ's. It is easy to see that this solves the problem in the case of the deep under-sampling regime considered above. There, only one singular value exists, and when this is neglected, one derives the resultĝ µ = 0, ∀µ again. In order to illustrate this procedure in a more general setting, we turn to the specific case of the U.S. Supreme Court data [8].

A Real-World Example
We have applied the inference scheme to the data of [8] that refer to the decisions of the U.S. Supreme Court on 895 cases. The U.S. Supreme Court is composed of nine judges, each of whom casts a vote against (s i = −1) or in favor (s i = +1) of a given case. Therefore, this is a n = 9 spin system for which we have N = 895 observations. The work in [8] has fitted this dataset with a fully-connected pairwise spin model. We refer to [8] for details on the dataset and on the analysis. The question we wish to address here is whether the statistical dependence between judges of the U.S. Supreme Court can really be described as a pairwise interaction, which hints at the direct influence of one judge on another one, or whether higher order interactions are also present.
In order to address this issue, we also studied a dataset n = 9 spins generated form a pairwise interacting model, Equation (22), from which we generated N = 895 independent samples. The value of β = 2.28 was chosen so as to match the average value of two-body interactions fitted in the true dataset. This allows us to test the ability of our method to recover the correct model when no assumption on the model is made.
As discussed above, the procedure discussed in the previous section yields estimatesĝ µ that allow us to recover empirical averages of all the operators. These, for a finite sample size N, are likely to be affected by considerable noise that is expected to render the estimatedĝ µ extremely unstable. In particular, since the sample contains unobserved states, i.e., states with k s = 0, we expect some of the parameters g µ to diverge or, with a finite numerical precision, to attain large values. Therefore, we also performed inference considering only the components with largest Λ λ in the singular value decomposition. Table 1 reports the values of the estimated parametersĝ λ obtained for the U.S. Supreme Court considering only the top = 2 to 7 singular values, and it compares them to those obtained when all singular values are considered. We observe that when enough singular values are considered, the estimated couplings converge to stable values, which are very different from those obtained when all 18 singular values are considered. This signals that the instability due to unobserved states can be cured by neglecting small singular values Λ λ 1.  This is confirmed by Figure 5, which shows that estimates ofĝ µ are much more stable when few singular values are considered (top right panel). The top left panel, which refer to synthetic data generated from Equation (22), confirms this conclusion. The estimatesĝ µ are significantly larger for a two-body interaction than for higher order and one-body interactions, as expected. Yet, when all singular values are considered, the estimated values of a two-body interaction fluctuate around values that are much larger than the theoretical one (β/n 0.2533 . . .) and the ones estimated from fewer singular values. In order to test the performance of the inferred couplings, we measure for each operator µ the change: (24) in log-likelihood when g µ is set to zero. If ∆ µ is positive or is small and negative, the coupling g µ can be set to zero without affecting much the ability of the model to describe the data. A large and negative ∆ µ instead signals a relevant interaction g µ . Clearly, ∆ µ ≤ 0 for all µ whenĝ µ is computed using all the q components. This is because in that case, the log-likelihood reaches the maximal value it can possibly achieve. When not all singular values are used, ∆ µ can also attain positive values. Figure 5 confirms our conclusions that inference using all the components is unstable. Indeed for the synthetic data, the loss in likelihood is spread out on operators of all orders, when all singular values are considered. When few singular values are considered, instead, the loss in likelihood is heavily concentrated on two body terms ( Figure 5, bottom left). Pairwise interactions stick out prominently because ∆ µ < 0 for all two-body operators µ. Still, we see that some of the higher order interactions, with even order, also generate significant likelihood losses.
With this insight, we can now turn to the U.S. Supreme Court data, focusing on inference with few singular values. Pairwise interactions stick out as having both sizableĝ µ ( Figure 5, top right) and significant likelihood loss ( Figure 5, bottom right). Indeed, the top interactions (those with minimal ∆ µ ) are prevalently pairwise ones. Figure 6 shows the hypergraph obtained by considering the top 15 interactions [22], which are two-or four-body terms (see the caption for details). Comparing this with synthetic data, where we find that the top 19 interactions are all pairwise, we conjecture that four-body interactions may not be spurious. The resulting network clearly reflects the orientation of individual judges across an ideological spectrum going from liberal to conservative positions (as defined in [8]). Interestingly, while two-body interactions describe a network polarized across this spectrum with two clear groups, four-body terms appear to mediate the interactions between the two groups. The prevalence of two-body interactions suggests that direct interaction between the judges is clearly important, yet higher order interactions seem to play a relevant role in shaping their collective behavior. Judges are represented as nodes with labels referring to the initials (as in [8]). Two-body interactions are represented by (red) links of a width that increases with |∆ µ |, whereas four-body interactions as (green) shapes joining the four nodes. The shade of the nodes represents the ideological orientation, as reported in [8], from liberal (black) to conservative (white).
As in the analysis in [8], single-body terms, representing a priori biases of individual judges, are not very relevant [23].

Conclusions
The present work represents a first step towards a Bayesian model selection procedure for spin models with interactions of arbitrary order. Rather than tackling the problem directly, which would imply comparing an astronomical number of models even for moderate n, we show that model selection can be performed first on mixture models, and then, the result can be projected in the space of spin models. This approach spots symmetries between states that occur with a similar frequency, which impose constraints between the parameters g µ . As we have seen, in simple cases, these symmetries are enough to recover the correct sparse model, imposing that g µ = 0 for all those interactions φ µ that are not consistent with the symmetries. These symmetries allow us to derive a set of sufficient statistics ψ λ (s) (the relevant variables) whose empirical values allow one to derive the maximum likelihood parametersĝ λ . The number q of sufficient statistics is given by the number of sets in the optimal partition Q * of states. Therefore, the dimensionality of the inference problem is not related to the number of different interaction terms φ µ (s) (or equivalently, of parameters g µ ), but it is rather controlled by the number of different frequencies that are observed in the data. As the number N of samples increases, q increases and so does the dimensionality of the inference problem, until one reaches the well-sampled regime (N 2 n ) when all states s are well resolved in frequency. It has been observed [11] that the family of probability distributions of the form (5) is endowed with a hierarchical structure that implies that high-order and low-order interactions are entangled in a nontrivial way. For example, we observe a non-trivial dependence between two-and four-body interactions. On the other hand, [18] shows that the structure of interdependence between operators in a model is not simply related to the order of the interactions and is invariant with respect to gauge transformations that do not conserve the order of operators. This, combined with the fact that our approach does not introduce any explicit bias to favor an interaction of any particular order, suggests that the approach generates a genuine prediction on the relevance of interactions of a particular order (e.g., pairwise). Yet, it would be interesting to explore these issues further, combining the quasi-orthogonal decomposition introduced in [11] with our approach.
It is interesting to contrast our approach with the growing literature on sloppy models (see, e.g., [19]). Transtrum et al. [19] have observed that inference of a given model is often plagued by overfitting that causes large errors in particular combinations of the estimated parameters.
Our approach is markedly different in that we stem right from the beginning from Bayesian model selection, and hence, we rule out overfitting from the outset. Our decomposition in singular values identifies those directions in the space of parameters that allow one to match the empirical distribution while preserving the symmetries between configurations observed with a similar frequency.
The approach discussed in this paper is only feasible when the number of variables n is small. Yet, the generalization to a case where the set M of interactions is only a subset of the possible interactions is straightforward. This entails setting to zero all couplings g µ relative to interactions µ ∈ M. Devising decimation schemes for doing this in a systematic manner, as well as combining our approach with regularization schemes (e.g., LASSO) to recover sparse models comprise a promising avenue of research for exploring the space of models.

Appendix A. The Completeness Relation (8)
We notice that the set of operators satisfies the following orthogonality relations: Then, taking the logarithm of Equation (5), multiplying by φ ν (s)/2 n with ν > 0 and summing over s, one finds that: Combining the above identity with the expression with Equation (6) finally yields Equation (8).

Appendix B. The Posterior Distribution of
Following [16], we assume a Dirichlet prior for the parameter vector , i.e., This is a natural choice due to the fact that it is a conjugate prior for the parameters of the multinomial distribution [24]. This means that the posterior distribution has the same functional form as the prior, and the a parameters can be interpreted as pseudocounts. In other words, the posterior probability of ρ is still a Dirichlet distribution, i.e., where k s is the number of times state s i was observed in the sample. We remind that the choice a = 0.5 corresponds to the least informative (Jeffrey's) prior [25], whereas with a = 0, the expected values of the parameters q coincide with the maximal likelihood estimates. The likelihood of the sampleŝ under model Q is given by: where: is the number of sample points in partition Q j and m j = |Q j | is the number of states in partition Q j . The work in [16] shows that the partition that maximizes the likelihood is the one where states with similar frequencies are grouped together, which is a coarse-graining of the K partition.
(A15) Therefore, the covariance matrix among the elements of ρ is composed of a diagonal part plus a part proportional to the identity matrix.