Version:

###

Entropy-based Characterization
of Modeling Constraints

Orestis Loukas111E-mail<EMAIL_ADDRESS>Ho Ryun
Chung222E-mail<EMAIL_ADDRESS>

Institute for Medical Bioinformatics and Biostatistics
Philipps-Universität Marburg
Hans-Meerwein-Straße 6, 35032 Germany

### Abstract

In most data-scientific approaches, the principle of Maximum Entropy (maxent)
is used to a posteriori justify some parametric model which has been already
chosen based on experience, prior knowledge or computational simplicity. In a
perpendicular formulation to conventional model building, we start from the
linear system of phenomenological constraints and asymptotically derive the
distribution over all viable distributions that satisfy the provided set of
constraints. The maxent distribution plays a special role, as it is the most
typical among all phenomenologically viable distributions representing a good
expansion point for large-$N$ techniques. This enables us to consistently
formulate hypothesis testing in a fully-data driven manner. The appropriate
parametric model that is supported by the data can be always deduced at the
end of model selection. In the maxent framework, we recover major scores and
selection procedures used in multiple applications and assess their ability to
capture associations in the data-generating process and identify the most
generalizable model. This data-driven counterpart of standard model selection
demonstrates the unifying prospective of the deductive logic advocated by
maxent principle, while potentially shedding new insights to the inverse
problem.

##  1 Introduction

Given phenomenological information in form of marginals, averages or higher
moments of features modeling seeks to determine some class of statistical
models (i.e. probability distributions) that best captures trends or laws in
the underlying system that generated the observed data. Since phenomenological
data is by definition a physical invariant of the problem at hand, it is most
naturally expressed by the invariant form of the defining set of equations
which constrain the space of model distributions to reproduce the observed
moments. As long as constraints remain of linear form, at least, this poses a
well-defined problem eliminating ambiguities related to reparametrization
symmetries in phenomenological information.

For a long time, it has been known that the distribution which satisfies the
problem-defining constraints while maximizing the entropy functional is in a
sense the least-biased distribution [1, 2] that incorporates the provided
information and nothing more. Aside from this conceptually appealing
interpretation, the entropy-based approach to modeling offers many benefits of
practical significance. Given certain information from phenomenology entropy
maximization uniquely selects one model distribution as the most
representative among all phenomenologically viable distributions offering a
clear deductive reasoning to modeling. At the same time, the maxent logic
combinatorially assigns probabilities to any other distribution in the most
intuitive way, namely according to its “distance” from the provided set of
observations. As this point appears to be of fundamental importance in the
maxent logic, we start our more formal parameter-agnostic exploration by
investigating the distribution over phenomenologically viable distributions.

Since entropy measures information flowing into our model, the maxent
principle applied on the linear system of phenomenological constraints turns
also classification tasks into clear-cut deductive procedures. In this work,
we explain how to unambiguously define model architectures and directly count
the physical degrees of freedom defining model dimensionality in order to
ultimately perform selection of significant constraints that can be
distinguished from sampling noise. Furthermore, we examine how large-order
expansions of information-theoretic measures and benchmark scores can be
meaningfully evaluated, when performed around the distribution of maximum
entropy.

More specifically, we rigorously show that expansions in the fluctuations of
probabilities around their maxent estimate at large sample size $N$ result in
leading Gaussian kernels with terms that are quadratic in the fluctuating
degrees of freedom. Such formal observation justifies writing discrete
distributions over distributions into multivariate normal distributions in the
fluctuations around the maxent saddle point. Asymptotically, this suggests the
integral representation of expectation sums over viable distributions which
can be in turn evaluated using standard large-$N$ localization techniques.
Eventually, we arrive at hypothesis testing in the maxent world and suggest
significance level(s) compatible with the large-$N$ expansion. In addition, we
review standard modeling scores outlining their meaning in the present
framework.

Under the unifying umbrella of maxent reasoning, we test accuracy of
(hyper)graph reconstruction and generalizability of models (i.e. maxent
distributions satisfying sets of constraints) trained on synthetic data from
an effective Ising magnet. Taking $N$ large verifies the theoretically
anticipated estimates of error bounds. Most crucially, it reveals the
universality of many measures which is expected to manifest at large order:
Information-theoretic “distances” between distributions asymptotically depend
only on the number of physical degrees of freedom and not on further details
of the phenomenological problem at hand. In particular, most model selection
procedures seem to merge into a universal rule that accurately detects the
asymptotic structure.

All in all, the constraint optimization of entropy gives a robust definition
of the notion of a model assigning to any distinct set of consistent
constraints its maxent distribution. At an operational level, the simplicity
and flexibility of directly working with physical degrees of freedom, i.e.
probabilities, could help efficiently tackle inverse problems in data science
for which parameters (lacking microscopic justification) often lead to
distractingly cumbersome inference. In addition, it enables the large-$N$
estimation of important information-theoretic metrics, such as the expected
entropy, training and test errors. Under a unifying framework, the maxent
principle seems to naturally provide the tools to consistently perform
application-adaptive selection of the most suitable model reflecting
underlying mechanisms in a fully data-driven fashion.

##  2 Phenomenological constraints

Any phenomenological problem is related to some probability distribution over
physical realizations. In this work, we concretely investigate statistical
systems of categorical features each assuming a distinct number of states. All
admissible realizations of these features comprise the space of microstates
$\mathcal{A}$. A distribution $\mathfrak{p}$ over $\mathcal{A}$ can be thought
of as an ${\left|\mathcal{A}\right|}$-dimensional vector living on the
probability simplex $\mathcal{P}$ which assigns a probability to each
microstate. Modeling seeks to determine some optimal distribution(s) in
$\mathcal{P}$ that satisfies $D\geq 1$ phenomenological constraints (one is
the normalization condition in $\mathcal{P}$) which we denote by a vector
$\widehat{\mathfrak{m}}$.

The most obvious and easy to handle type of constraints is of linear form,

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}R_{a\alpha}\,\mathfrak{p}^{\alpha}\overset{!}{=}\widehat{m}_{a}\quad\text{for}\quad
a=1,\ldots D~{},$ (1)

based on an architecture matrix $\mathbf{R}$ with each row representing a
constraint and columns corresponding to the microstates in $\mathcal{A}$. In
the following, we take the architecture matrix to be of full row rank
$D=\text{rank}\mathbf{\mathbf{R}}$ directly specifying the number of
independent constraints. In realistic applications however, a redundant set of
linear constraints summarized by a coefficient matrix $\mathbf{C}$ acting on
$\mathfrak{p}$ is often provided, instead. Whenever the phenomenological input
concretely consists of a set of marginal relative frequencies, $\mathbf{C}$
becomes a binary matrix dictating the relative frequencies of microstates to
be summed over to obtain marginal relative frequencies in the empirical set
(i.e. each row of the coefficient matrix describes a marginal constraint sum).
To obtain a robust model definition it is instructive to transform the
coefficient matrix $\mathbf{C}$ to its reduced row-echelon form. Dropping any
all-zero row corresponding to redundant constraints, one then obtains the
architecture matrix $\mathbf{R}$ which uniquely characterizes the linear
system at hand.

Noting in particular that there are by construction $D$ columns in
$\mathbf{R}$ with all-zeros but one (the leading) 1, it can be easily seen
that normalization of $\mathfrak{p}\in\mathcal{P}$ implies

$\displaystyle\sum_{a=1}^{D}R_{a\alpha}=1\quad\forall\,\alpha=1,\ldots,{\left|\mathcal{A}\right|}\quad\text{and}\quad\sum_{a=1}^{D}\widehat{m}_{a}=1~{}.$
(2)

The solution set of linear system Eq. (1), which is unambiguously defined by
architecture matrix $\mathbf{R}$ and inhomogeneous vector
$\widehat{\mathfrak{m}}$, induces an equivalence class in
$\mathbb{R}^{\left|\mathcal{A}\right|}$. To see this first note that the
empirical distribution $\mathfrak{f}\in\mathcal{P}$ which summarizes the
observed relative frequencies of all microstates in $\mathcal{A}$ always
satisfies the coupled equations, thus being member of the equivalence class
which we denote by $[\mathfrak{f}]_{\mathbf{R}}$ in the following. By
construction, any distribution $\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}$
shares the same statistics $\widehat{\mathfrak{m}}$ as the empirical
distribution. Whenever only normalization constraint is specified as an all-
one row in $\mathbf{R}$, the equivalence class covers the full simplex,
$[\mathfrak{f}]_{(1,\dots,1)}=\mathcal{P}$, whereas in the other extreme case
where the architecture matrix becomes the identity, all frequencies are
learned by heart, so that $[\mathfrak{f}]_{\text{\small
1}\text{1}}=\\{\mathfrak{f}\\}$.

In a general setup, there might exist microstates whose probability is fixed
in the problem statement to some (non)-finite empirical value. If some
empirical marginal is known to be zero by logical deduction, for example, all
the associated probabilities appearing in this marginal sum are necessarily
set to zero due to the non-negativity on the simplex $\mathcal{P}$. In the
following, we concentrate on all constraints and associated microstates that
are not a priori determined by some empirical or logical argument removing any
probabilities that are learned by heart from empirical data. The latter can be
trivially reinstated at the end of model selection without modifying
conclusions. In other words, we understand $\mathcal{A}$ to only include those
microstates with $\mathfrak{p}_{\alpha}>0$ that retain some probabilistic
interest for model building.

##  3 Entropy concentration

Mathematically, the notion of an optimal solution amounts to selecting a
target functional to be extremized over the equivalence class associated to
the system of phenomenological constraints. Choosing [3, 4, 5] as a measure of
uncertainty Shannon’s differential entropy,

$\displaystyle
H[\mathfrak{p}]=-\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\mathfrak{p}_{\alpha}\log\mathfrak{p}_{\alpha}\quad\text{with}\quad\mathfrak{p}\in\mathcal{P}~{},$
(3)

we are naturally led to the principle of maximum entropy (maxent) under linear
constraints Eq. (1). In this way, we ensure that only the information
specified by the linear system in Eq. (1) enters into the model distribution
and nothing more. If no constraints were specified besides the dimensionality
of $\mathcal{A}$, then the distribution of maximum entropy would equal the
uniform distribution $\mathfrak{u}_{\alpha}={\left|\mathcal{A}\right|}^{-1}$
as a manifestation [6] of having the same degree of ignorance about any two
microstates in $\mathcal{A}$. In the other extreme scenario that all
${\left|\mathcal{A}\right|}$ empirical frequencies were constrained, Eq. (1)
would unambiguously fix the maxent distribution to the empirical distribution
$\mathfrak{f}$, as previously discussed. Most interestingly, we are concerned
with all those cases were $\mathbf{R}$ is column-rank deficit thus leading to
a non-trivial equivalence class $[\mathfrak{f}]_{\mathbf{R}}$ from which the
maxent principle selects a representative distribution that we henceforth
denote by $\hat{\mathfrak{p}}$. As established in [7] the maxent distribution
always exists and is unique for any consistent set of linear constraints.

Suppose we are now given a dataset with $N$ records. From elementary
combinatorics, we know that the multinomial coefficient (r.h.s. gives the
large-$N$ expansion)

$\displaystyle
W[N\mathfrak{p}]=\frac{N!}{\prod_{\alpha=1}^{\left|\mathcal{A}\right|}(N\mathfrak{p}^{\alpha})!}=\exp\left\\{NH[\mathfrak{p}]-\frac{\left|\mathcal{A}\right|-1}{2}\log
N+\mathcal{O}\left(1\right)\right\\}$ (4)

counts the possible ways to separate those $N$ individuals into
${\left|\mathcal{A}\right|}$ teams such that the counts $N\mathfrak{p}$ are
reproduced. As we have no way to distinguish among realizations of the
$N\mathfrak{p}$ counts, we should assign equal probabilities to them. Under
this thermodynamic assumption, member distributions in the equivalence class
formally distribute according to

$\displaystyle
p\,(\mathfrak{p}\;|\;N\mathfrak{f},\mathbf{R})=\left(\sum_{\mathfrak{p}^{\prime}\in[\mathfrak{f}]_{\mathbf{R}}}\exp\left\\{NH[\mathfrak{p}^{\prime}]\right\\}\right)^{-1}\exp\left\\{NH[\mathfrak{p}]\right\\}~{},$
(5)

at sample sizes $N\gg 1$. Incorporating model-dependent effects of
$\mathcal{O}\left(1\right)$ in the large-$N$ expansion of Eq. (4) would result
in $\mathcal{O}\left(1/N\right)$ effects in any relevant expectation over the
distribution of distributions in the equivalence class. At sufficiently large
$N$ we may thus neglect subleading corrections in Stirling’s formula used to
expand factorials and directly start from Eq. (5).

In particular, the ratio $W[N\hat{\mathfrak{p}}]/W[N\mathfrak{p}]$ for any
distribution $\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}$ measures the
relative abundance of the maxent distribution over another member in the
equivalence class. The large-$N$ expansion of Eq. (4) leads to

$\displaystyle\log\frac{W[N\hat{\mathfrak{p}}]}{W[N\mathfrak{p}]}=N\left(H[\hat{\mathfrak{p}}]-H[\mathfrak{p}]\right)+\mathcal{O}\left(1\right)~{}.$
(6)

Taking the $N\rightarrow\infty$ limit, we recognize [8] that the maxent
distribution can be realized in more ways compared to any other distribution
in $[\mathfrak{f}]_{\mathbf{R}}$, as the entropy difference will be by
definition positive. Consequently, $\hat{\mathfrak{p}}$ is the distribution we
should overwhelmingly anticipate given the empirical constraints, the latest
asymptotically, since any other distribution would be exponentially
suppressed, justifying our choice of target functional. Intuitively, the
dominating role of $\hat{\mathfrak{p}}$ over any other distribution in the
equivalence class shows that the maxent distribution is in a sense typical of
this set.

In order to prove the main and subsequent theorems in this paper we shall need
a fairly general

######  Lemma 1.

Let $\hat{\mathfrak{p}}$ be the maxent distribution in equivalence class
$[\mathfrak{f}]_{\mathbf{R}}$ induced by architecture $\mathbf{R}$ on data
described by empirical distribution $\mathfrak{f}\in\mathcal{P}$. For any
distribution $\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}$ with linearized
fluctuations

$\displaystyle\mathfrak{p}_{\alpha}=\hat{\mathfrak{p}}_{\alpha}+\sqrt{\frac{\hat{\mathfrak{p}}_{\alpha}}{N}}\,\boldsymbol{\pi}_{\alpha}$
(7)

around the maxent distribution in the equivalence class, the large-$N$
expansion of the entropy difference

$\displaystyle
2N\left(H[\hat{\mathfrak{p}}]-H[\mathfrak{p}]\right)=\sum_{\alpha=1}^{|\mathcal{A}|}\boldsymbol{\pi}_{\alpha}^{2}+\mathcal{O}\left(N^{-1/2}\right)$
(8)

results into a leading term which is quadratic in the fluctuations
$\boldsymbol{\pi}\in\mathbb{R}^{\left|\mathcal{A}\right|}$.

#### Proof

As long as111The case of $\hat{\mathfrak{p}}_{\alpha}=0$ exhibits no
probabilistic fluctuations and does not contribute to our formulas. Formally,
this is manifested by taking $0\log 0=0$ in any information-theoretic measure.
$\mathfrak{p}_{\alpha}\neq 0$, parametrizing fluctuations according to Eq. (7)
does not degenerate. Generically, normalization on the simplex $\mathcal{P}$
means

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\boldsymbol{\pi}_{\alpha}\sqrt{\hat{\mathfrak{p}}_{\alpha}}=0~{}.$
(9)

Furthermore, invoking a well-known fact [9, 10] from information theory about
maxent distributions $\hat{\mathfrak{p}}$,

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\left(\mathfrak{p}_{\alpha}-\hat{\mathfrak{p}}_{\alpha}\right)\log\hat{\mathfrak{p}}_{\alpha}=0\quad\forall\,\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}~{},$
(10)

we recognize that

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\boldsymbol{\pi}_{\alpha}\sqrt{\hat{\mathfrak{p}}_{\alpha}}\log\hat{\mathfrak{p}}_{\alpha}=0~{},$
(11)

whenever the maxent distribution in $[\mathfrak{f}]_{\mathbf{R}}$ is chosen as
expansion point. Combining Eq. (9) with Eq. (11) in the Taylor expansion of
the entropy difference results into Eq. (8). The factor of $\sqrt{N}$ in
parametrization of Eq. (7) is chosen [11] so that the quadratic term in Eq.
(8) scales as $\mathcal{O}\left(1\right)$ in the large-$N$ expansion ensuring
the proper asymptotic limit by suppressing higher-order corrections.

We are now in a position to state

######  Theorem 1.

Asymptotically, the scaled difference
$2N(H[\hat{\mathfrak{p}}]-H[\mathfrak{p}])$ of the entropy of any distribution
$\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}$ from the maximum entropy in the
equivalence class induced by $\mathbf{R}$ on data described by $\mathfrak{f}$
is distributed as chi-squared with
$\left|\ker\mathbf{R}\right|={\left|\mathcal{A}\right|}-D$ degrees of freedom,
in short as $\chi_{{\left|\mathcal{A}\right|}-D}^{2}$. For any positive
$\delta$ the likelihood to find a fraction of distributions in
$[\mathfrak{f}]_{\mathbf{R}}$ having entropies in the range

$\displaystyle H[\hat{\mathfrak{p}}]-\delta\leq H[\mathfrak{p}]\leq
H[\hat{\mathfrak{p}}]$

is accordingly given by the cumulative distribution
$F_{{\left|\mathcal{A}\right|}-D}(2N\delta)$ of
$\chi_{{\left|\mathcal{A}\right|}-D}^{2}$ distribution.

#### Proof

At sufficiently large $N$, we are interested to evaluate the fraction of
distributions in the equivalence class having entropies in a range $\delta>0$
below maxent $H[\hat{\mathfrak{p}}]$, i.e. in the region

$\displaystyle\mathcal{R}=\left\\{\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}\,\big{|}\,H[\mathfrak{p}]\geq
H[\hat{\mathfrak{p}}]-\delta\right\\}~{}.$ (12)

Formally using the distribution $p$ over distributions introduced in Eq. (5)
this fraction is written as

$\displaystyle F_{k}(2N\delta)=\sum_{\mathfrak{p}\in\mathcal{P}}\text{\small
1}\text{1}_{\mathcal{R}}\,p\,(\mathfrak{p}\;|\;N\mathfrak{f},\mathbf{R})=\left(\sum_{\mathfrak{p}^{\prime}\in[\mathfrak{f}]_{\mathbf{R}}}\exp\left\\{-N(H[\hat{\mathfrak{p}}]-H[\mathfrak{p}^{\prime}])\right\\}\right)^{-1}\sum_{\mathfrak{p}\in\mathcal{R}}\exp\left\\{-N(H[\hat{\mathfrak{p}}]-H[\mathfrak{p}])\right\\}~{},$
(13)

($\text{\small 1}\text{1}_{\mathcal{R}}$ denotes the indicator function),
which is clearly governed by the entropy difference
$H[\hat{\mathfrak{p}}]-H[\mathfrak{p}]$ of each member distribution from
maxent. For the purpose of evaluating $p$ as an expansion in powers of $1/N$,
we set

$\displaystyle\boldsymbol{\pi}_{\alpha}=\sqrt{\frac{N}{\hat{\mathfrak{p}}_{\alpha}}}\left[\mathfrak{p}_{\alpha}-\hat{\mathfrak{p}}_{\alpha}\right]\quad\forall\,\alpha=1,\ldots,{\left|\mathcal{A}\right|}$
(14)

such that $\boldsymbol{\pi}_{\alpha}=\mathcal{O}\left(\sqrt{N}\right)$.
Clearly, Lemma 1 applies for $\boldsymbol{\pi}$, so that the governing entropy
difference results in a quadratic term as in Eq. (8). On account of linear
system Eq. (1), fluctuations around any reference distribution in the
equivalence class satisfy homogeneous system

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\left(R_{a\alpha}\sqrt{\hat{\mathfrak{p}}_{\alpha}}\right)\boldsymbol{\pi}_{\alpha}=0~{}~{}\forall\,a=1,\ldots,D~{}.$
(15)

Due to these homogeneous constraints any
$\boldsymbol{\pi}\in\mathbb{R}^{|\mathcal{A}|}$ lives in the
$({\left|\mathcal{A}\right|}-D)$-dimensional cone spanned by the kernel of
architecture matrix $\mathbf{R}$ after column-wise rescaling by
$\sqrt{\hat{\mathfrak{p}}_{\alpha}}>0$. At an operational level, the quadratic
sum in Eq. (8) expressed in any valid coordinate system parameterizing the
kernel of the rescaled architecture matrix can be diagonalized with unity-det
Jacobian. To this end, we can write the fluctuations in Eq. (7) around maxent
solution as

$\displaystyle\boldsymbol{\pi}_{\alpha}=\sum_{i=1}^{{\left|\mathcal{A}\right|}-D}x_{i}X_{i\alpha}~{},$
(16)

where $\mathbf{X}$ orthonormally encodes the kernel space, i.e. for
$a=1,\dots,D$ and $i,j=1,\dots,{\left|\mathcal{A}\right|}-D$ it is

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\left(R_{a\alpha}\sqrt{\hat{\mathfrak{p}}_{\alpha}}\right)X_{i\alpha}=0\quad\text{and}\quad\sum_{\alpha=1}^{\left|\mathcal{A}\right|}X_{i\alpha}X_{j\alpha}=\delta_{ij}~{}.$
(17)

Hence, the leading quadratic term in the entropy difference becomes

$\displaystyle\sum_{\alpha=1}^{|\mathcal{A}|}\boldsymbol{\pi}_{\alpha}^{2}=\sum_{i=1}^{{\left|\mathcal{A}\right|}-D}x_{i}^{2}\,\sim\,\chi^{2}_{{\left|\mathcal{A}\right|}-D}~{},$
(18)

readily following a chi-square distribution with
${\left|\mathcal{A}\right|}-D$ degrees of freedom.

Exploiting the radial symmetry of the large-$N$ expansion we subsequently
introduce spherical coordinates such that any
$\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}$ living on the
$({\left|\mathcal{A}\right|}-D-1)$ – sphere with squared radius
$\chi^{2}=\sum_{i=1}^{{\left|\mathcal{A}\right|}-D}x_{i}^{2}$ has leading
entropy

$\displaystyle H[\mathfrak{p}]=H[\hat{\mathfrak{p}}]-(2N)^{-1}\chi^{2}~{}.$
(19)

Notice that the uniqueness of the optimal solution unambiguously fixes the
origin of the chosen spherical coordinate system to correspond to the maxent
distribution in the given equivalence class. Conveniently, the spherical
parametrization allows us to write the integral representation of Eq. (13)
after performing the angular integration,

$\displaystyle
F_{{\left|\mathcal{A}\right|}-D}(2N\delta)=\frac{1}{2^{\frac{{\left|\mathcal{A}\right|}-D}{2}-1}\,\Gamma\left(\frac{{\left|\mathcal{A}\right|}-D}{2}\right)}\int_{0}^{\sqrt{2N\delta}}\text{d}\chi\,\chi^{{\left|\mathcal{A}\right|}-D-1}\,e^{-\chi^{2}/2}~{},$
(20)

as the cumulative distribution of chi squared. This concludes the proof.

Based on radial symmetry, Theorem 1 provides [12] a neat geometrization of the
distribution over distributions Eq. (5) in the equivalence class. Around the
maxent distribution $\hat{\mathfrak{p}}$ in a
$\left|\ker\mathbf{R}\right|$-dimensional space, any spherical shell of radius
$\sqrt{2N\Delta H}$ traces over all those distributions
$\mathfrak{p}\in[\mathfrak{f}]_{\mathbf{R}}$ whose entropy differs from
$H[\hat{\mathfrak{p}}]$ by the same amount $\Delta H$. Clearly, the interval
within which the entropy of a specified portion of distributions is to be
found becomes narrower with larger $N$, showing that the maxent distribution
has an entropy more and more representative of the whole equivalence class.
Hence, the allowed by $\mathbf{R}$ entropies concentrate around the maxent in
a geometric sense. A posteriori, this formally justifies our choice of entropy
(which in the present setting rather emerges at large-$N$ than being
postulated) as a target functional for selecting a model distribution.

###  3.1 Calculating the optimal solution

The preceding discussion about the concentration of entropies around the
maxent and the application of chi-squared statistics presupposes that the
distribution of maximum entropy has already been determined. Of course, one
can parametrize the solution space of linear system Eq. (1) and try to
maximize analytically the entropy. For small dimensionality of
$\ker\mathbf{R}$, the solution could have a neat closed form. In other setups,
there might exist some perturbative method to proceed.

In a reverse approach, we directly incorporate the optimization objective into
the phenomenological system to obtain a non-linear system of coupled
equations, which can be iteratively solved via Newton-Raphson method. Given
architecture matrix $\mathbf{R}$ and empirical inhomogeneous vector
$\widehat{\mathfrak{m}}$, we start from the uniform distribution
$\mathfrak{p}^{(0)}=\mathfrak{u}$. Upon using the exponential form [9] of
maxent distribution in absence of zero-probability microstates,

$\displaystyle\hat{\mathfrak{p}}_{\alpha}=\exp\left(\sum_{a=1}^{D}\hat{\theta}_{a}R_{a\alpha}\right)~{},$
(21)

for some vector $\hat{\boldsymbol{\theta}}\in\mathbb{R}^{D}$ of Lagrange
multipliers, the update rule for Newton estimate of the maxent solution to
system Eq. (1) at the $n$-th step reads

$\displaystyle\mathfrak{p}_{\alpha}^{(n+1)}=\mathfrak{p}_{\alpha}^{(n)}\exp\left\\{-\sum_{a,b=1}^{D}R_{a\alpha}((\mathbf{J}^{(n)})^{-1})_{ab}\left(m_{b}^{(n)}-\widehat{m}_{b}\right)\right\\}\quad\text{for
}\alpha=1,\ldots,{\left|\mathcal{A}\right|}$ (22)

in terms of the previously estimated empirical moment

$\displaystyle
m_{a}^{(n)}=\sum_{\alpha=1}^{\left|\mathcal{A}\right|}R_{a\alpha}\,\mathfrak{p}^{(n)}_{\alpha}\quad,\quad
a=1,\ldots,D$ (23)

and Jacobian matrix

$\displaystyle
J^{(n)}_{ab}=\sum_{\alpha=1}^{\left|\mathcal{A}\right|}R_{a\alpha}\,\mathfrak{p}^{(n)}_{\alpha}\,R_{b\alpha}\quad,\quad
a,b=1,\ldots,D~{}.$ (24)

The latter matrix is invertible, as it can be easily seen by noting that
vanishing of its quadratic form for any vector $\mathfrak{u}\in\mathbb{R}^{D}$
implies

$\displaystyle
0\overset{!}{=}\sum_{a,b=1}^{D}u_{a}J^{(n)}_{ab}u_{b}=\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\mathfrak{p}^{(n)}_{\alpha}\left(\sum_{a=1}^{D}u_{a}R_{a\alpha}\right)^{2}\quad\Rightarrow\quad\sum_{a=1}^{D}u_{a}R_{a\alpha}=0\quad\Rightarrow\quad
u^{a}=0\quad\forall a=1,\ldots D~{},$ (25)

since all $\mathfrak{p}^{(n)}_{\alpha}>0$ (having excluded structural zeros
whose Newton estimate becomes trivially zero) while the rows of $\mathbf{R}$
are per assumption linearly independent. As customary, Newton method converges
fast to the unique maxent distribution $\hat{\mathfrak{p}}$, at the cost of a
$D\times D$ matrix inversion at each iteration step.

There exist many variations of the main iterative scheme. Suppose that instead
of $\mathbf{R}$ a binary matrix $\mathbf{C}$ which describes a possibly
redundant set of $M$ marginal constraints is used. Cycling through each
marginal constraint $a=1,\ldots M$ to iteratively update probabilities via
one-dimensional Newton method (as the multidimensional Newton method would be
plagued by the singularity of redundancies) gives the rule

$\displaystyle\mathfrak{p}^{(nM+a)}_{\alpha}=\mathfrak{p}^{(nM+a-1)}_{\alpha}\exp\left\\{C_{a\alpha}\left(\frac{\widehat{m}_{a}}{m^{(nD+a-1)}_{a}}-1\right)\right\\}\approx\mathfrak{p}^{(nM+a)}_{\alpha}\left(\frac{\widehat{m}_{a}}{m^{(nM+a-1)}_{a}}\right)^{C_{a\alpha}}~{},$
(26)

which is determined by the ratio of empirical marginals to the running
estimate of model marginals

$\displaystyle
m^{(nM+a-1)}_{a}=\sum_{\alpha=1}^{\left|\mathcal{A}\right|}C_{a\alpha}\mathfrak{p}^{(nM+a-1)}_{\alpha}~{}.$
(27)

One can show [9, 13] that the Taylor-expanded r.h.s. of Eq. (26) always
converges to the maxent solution. This approach leads to the algorithm of
iterative proportional fitting [14].

##  4 Constraint selection

In this section, we exploit the logic advocated in Section 3 as our starting
point in order to motivate a fully data-driven selection of phenomenological
constraints. In most cases, we are given summary statistics of an empirical
distribution, i.e. a set of phenomenological constraints summarized by
$\widehat{\mathbf{C}}$ whose invariant form $\widehat{\mathbf{R}}$ induces an
equivalence class $[\mathfrak{f}]_{\widehat{\mathbf{R}}}\subseteq\mathcal{P}$,
as outlined in Section 2. This robustly defines the notion of a model in a
purely data-driven fashion. In addition, a chi-square distribution
asymptotically assigns probabilities to all member distributions in the
equivalence class based on their entropy relative to maxent. With increasing
sample size $N$, an increasing number of distributions has entropies
concentrating around maxent, making $\hat{\mathfrak{p}}$ typical of
$[\mathfrak{f}]_{\widehat{\mathbf{R}}}$ and hence a good candidate of a model
distribution.

###  4.1 Hypothesis testing

#### The hyper MaxEnt.

If the given architecture matrix $\widehat{\mathbf{R}}$ defined all our
knowledge about the system under investigation, we first ask how likely it
would be to observe the empirical distribution $\mathfrak{f}$ whose entropy
differs from the maxent by

$\displaystyle\hat{\delta}=H[\hat{\mathfrak{p}}]-H[\mathfrak{f}]~{}.$ (28)

On account of Theorem 1, the asymptotic likelihood to observe the empirical
distribution – if the $D$ constraints in $\widehat{\mathbf{R}}$ are all we
know/trust from the data – is naturally given as the $p$-value under
cumulative $\chi^{2}_{{\left|\mathcal{A}\right|}-D}$ distribution to record a
value $2N\hat{\delta}$ or more extreme. If the empirical $p$-value lies above
a pre-specified threshold, we could argue in favor of the maxent distribution
$\hat{\mathfrak{p}}$ as capturing all essential information, any discrepancy
to $\mathfrak{f}$ attributed to mere statistical noise. In the opposite case
of a $p$-value below the threshold, we would have to deduce that the provided
constraints were not sufficient to capture mechanisms of the data-generating
process, $\hat{\mathfrak{p}}$ systematically failing to adequately
characterize $\mathfrak{f}$.

In principle, full access to empirical distribution $\mathfrak{f}$ makes a
classification of all possible – given the microstate space $\mathcal{A}$ –
architectures $\widehat{\mathbf{R}}$ and corresponding generalized moments
$\widehat{\mathfrak{m}}$ amenable. Any set of constraints which fails to
capture the systematics of the data-generating process, so that the observed
data $\mathfrak{f}$ becomes very atypical in the induced equivalence class, is
immediately vetoed. In the spirit of Occam’s razor, we select the simplest
model with the highest entropy which still makes the empirical dataset likely
(empirical $p$-value above the threshold) as our starting point. We shall call
the model incorporating all significant constraints – i.e. constraints whose
removal leads to significantly low $p$-values – and nothing more, the hyper-
maxent.

#### Likelihood ratio test.

For those constraints whose removal does not bring any significant result, as
far as the empirical entropy difference
$H[\hat{\mathfrak{p}}]-H[\mathfrak{f}]$ is concerned, we could proceed by
relative comparison of constraints according to

######  Theorem 2.

Given is a base architecture $\widehat{\mathbf{R}}$ and a more complex
$\widehat{\mathbf{R}}^{\prime}$ which implies the former. Asymptotically, the
scaled difference $2N(H[\hat{\mathfrak{p}}]-H[\hat{\mathfrak{p}}^{\prime}])$
of maxent in
$[\mathfrak{f}]_{\widehat{\mathbf{R}}^{\prime}}\subset[\mathfrak{f}]_{\widehat{\mathbf{R}}}$
from maxent in $[\mathfrak{f}]_{\widehat{\mathbf{R}}}$ is distributed as chi-
squared with $\widehat{D}^{\prime}-\widehat{D}$ degrees of freedom.

Evidently, Theorem 2 reduces to Theorem 1 using as most complex model
architecture $\widehat{\mathbf{R}}^{\prime}=\text{\small 1}\text{1}$ with
$\widehat{D}^{\prime}={\left|\mathcal{A}\right|}$ that learns by heart.

#### Proof

This theorem describes a frequently encountered task in model building, to
namely quantify the suitability of a simpler set of constraints
$\widehat{\mathbf{R}}$ over a more stringent $\widehat{\mathbf{R}}^{\prime}$
which implies the former. At the level of architecture matrices, such nested
models are related via a linear map

$\displaystyle\widehat{R}_{a\alpha}=\sum_{a^{\prime}=1}^{\widehat{D}^{\prime}}T_{aa^{\prime}}\widehat{R}^{\prime}_{a^{\prime}\alpha}\quad,\quad
a=1,\dots,\widehat{D}$ (29)

for some transformation matrix $\mathbf{T}$ of
$\text{rank}\mathbf{T}=\widehat{D}$ (otherwise $\widehat{\mathbf{R}}$ would be
row-rank deficit). In our framework, both models are unambiguously
characterized by maxent distributions $\hat{\mathfrak{p}}$ and
$\hat{\mathfrak{p}}^{\prime}$, respectively. Clearly, the form Eq. (21) of
maxent distribution $\hat{\mathfrak{p}}^{\prime}$ automatically accommodates
the reference distribution of the simpler model, so that there always exists
vector $\hat{\boldsymbol{\theta}}$ with
$\hat{\mathfrak{p}}^{\prime}=\hat{\mathfrak{p}}$ given empirical data.

In a designer setup, different datasets that respect the reference constraints
in $\widehat{\mathbf{R}}$ are produced. By definition, a maxentist trusting
$\widehat{\mathbf{R}}$ would always conclude out of all those datasets that
$\hat{\mathfrak{p}}$ is the appropriate distribution. On the other hand, a
modeler trusting $\widehat{\mathbf{R}}^{\prime}$ would compute an a priori
different maxent distribution $\hat{\mathfrak{p}}^{\prime}$ depending on the
received dataset, but respecting the reference constraints from
$\widehat{\mathbf{R}}$. Writing $\hat{\mathfrak{p}}^{\prime}$ as in Eq. (7)
this translates into fluctuations $\boldsymbol{\pi}$ around the reference
maxent distribution $\hat{\mathfrak{p}}$ satisfying projection conditions

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\left(\widehat{R}_{a\alpha}\sqrt{\hat{\mathfrak{p}}_{\alpha}}\right)\boldsymbol{\pi}_{\alpha}=0\quad\forall\,a=1,\dots,\widehat{D}~{}.$
(30)

We are interested in all those cases that the latter modeler would select a
different distribution whereas the first would insist on $\hat{\mathfrak{p}}$.
To properly quantify this discrepancy we need to filter out redundancies
whenever different datasets lead to the same model distributions, as they
differ in constraints that none of the two modelers can see so that both
insist on their respective maxent models. In other words, fluctuations from
$\hat{\mathfrak{p}}^{\prime}$ around $\hat{\mathfrak{p}}$ must be
perpendicular to the kernel of $\widehat{\mathbf{R}}^{\prime}$. Hence, its
linearized fluctuations must look like

$\displaystyle\boldsymbol{\pi}_{\alpha}=\sum_{a^{\prime}=1}^{\widehat{D}^{\prime}}x_{a^{\prime}}\widehat{R}^{\prime}_{a^{\prime}\alpha}~{}.$
(31)

Reference projection Eq. (30) then constrains the vector
$x\in\mathbb{R}^{\widehat{D}^{\prime}}$ via

$\displaystyle\sum_{a^{\prime}=1}^{\widehat{D}^{\prime}}T_{aa^{\prime}}x_{a^{\prime}}=0\quad\forall\,a=1,\dots,\widehat{D}~{},$
(32)

upon using Eq. (29) Hence, $x$ lies in the kernel of $\mathbf{T}$ which is of
dimension $\widehat{D}^{\prime}-\widehat{D}$, so that in total

$\displaystyle\boldsymbol{\pi}_{\alpha}=\sum_{i=1}^{\widehat{D}^{\prime}-\widehat{D}}y_{i}\sum_{a^{\prime}=1}^{\widehat{D}^{\prime}}Y_{ia^{\prime}}\widehat{R}^{\prime}_{a^{\prime}\alpha}~{},$
(33)

where $\mathbf{Y}$ orthonormally spans the kernel of $\mathbf{T}$. Analogously
to the proof of Theorem 1, this form of fluctuations with
$\boldsymbol{\pi}\cdot\boldsymbol{\pi}=y\cdot y$ directly leads to the claimed
chi-squared distribution by applying Lemma 1.

In the large-$N$ expansion, the multinomial distribution for any two
distributions $\mathfrak{f},\mathfrak{q}\in\mathcal{P}$,

$\displaystyle\log\text{mult}(N\mathfrak{f};\mathfrak{q})=-ND_{\textsc{kl}}(\mathfrak{f}\;\|\;\mathfrak{q})-\frac{\left|\mathcal{A}\right|-1}{2}\log
N+\mathcal{O}\left(1\right)~{},$ (34)

is governed222The same applies for the multivariate hypergeometric
distribution, as long as the population size and population counts are much
larger than the size and frequencies of samples from this population. by the
kl divergence of distribution $\mathfrak{q}$ from empirical $\mathfrak{f}$. In
particular, the log-likelihood ratio to leading order in $N$ upon using Eq.
(10) equals

$\displaystyle\log\frac{\text{mult}(N\mathfrak{f};\hat{\mathfrak{p}}^{\prime})}{\text{mult}(N\mathfrak{f};\hat{\mathfrak{p}})}=N\left(H[\hat{\mathfrak{p}}]-H[\hat{\mathfrak{p}}^{\prime}]\right)+\mathcal{O}\left(1\right)~{},$
(35)

i.e. the observed entropy difference between maxent distributions of the more
complex to simpler model. Hence, we recover the likelihood ratio test [15]
(see also Lagrange multiplier test [16]) in the maxent framework, where
$\hat{\mathfrak{p}}$ is preferred, whenever it is sufficiently descriptive of
all those distributions satisfying the additional constraints in
$\widehat{\mathbf{R}}^{\prime}$ and nothing more. Conversely, we can formulate

######  Corollary 2.1.

For any two sets of constraints described by architectures
$\widehat{\mathbf{R}}$ and $\widehat{\mathbf{R}}^{\prime}$ where the former is
implied by the latter, the likelihood ratio of an observed dataset of sample
size $N$ is governed by the entropy difference
$\hat{\delta}=H[\hat{\mathfrak{p}}]-H[\hat{\mathfrak{p}}^{\prime}]$ between
the maxents in the induced equivalence classes. The maxent distribution
$\hat{\mathfrak{p}}\in[\mathfrak{f}]_{\widehat{\mathbf{R}}}$ is rejected as a
good model distribution in favor of
$\hat{\mathfrak{p}}^{\prime}\in[\mathfrak{f}]_{\widehat{\mathbf{R}}^{\prime}}$
at $\alpha$ significance level, whenever
$1-F_{\widehat{D}^{\prime}-\widehat{D}}(2N\hat{\delta})<\alpha$.

###  4.2 Information criteria

#### Maximum likelihood and BIC score.

In the equivalence class induced by architecture matrix $\widehat{\mathbf{R}}$
on counts $N\mathfrak{f}$, the probability density $\rho$ of observing entropy
difference $H[\hat{\mathfrak{p}}]-H[\mathfrak{p}]\equiv r^{2}/2$ is given
according to

$\displaystyle\rho\,(\mathfrak{p}\;|\;N\mathfrak{f},\widehat{\mathbf{R}})\,r^{{\left|\mathcal{A}\right|}-\widehat{D}-1}\text{d}r\,\text{d}\Omega=N^{\frac{{\left|\mathcal{A}\right|}-\widehat{D}}{2}}e^{-Nr^{2}/2}\,r^{{\left|\mathcal{A}\right|}-\widehat{D}-1}\text{d}r\,\text{d}\Omega~{},$
(36)

after redefining the entropic radius $\chi=\sqrt{N}r$ in Eq. (19) to extract
the large-order scaling which remains constant throughout the equivalence
class. Taking the logarithm of Eq. (36) evaluated at the empirical entropy
difference,

$\displaystyle\log\rho\,(\mathfrak{f}\;|\;N\mathfrak{f},\widehat{\mathbf{R}})=-N\left(H[\hat{\mathfrak{p}}]-H[\mathfrak{f}]\right)+\frac{{\left|\mathcal{A}\right|}-\widehat{D}}{2}\log
N=\mbox{$\frac{1}{2}$}{\left|\mathcal{A}\right|}\log
N+NH[\mathfrak{f}]-\mbox{$\frac{1}{2}$}\widehat{\textsc{bic}}~{},$ (37)

we recover from the model-dependent part the Bayesian Information Criterion
(in short bic score), first introduced in [17]. In the spirit of mle, the set
of constraints must be selected which make the data most likely to be
observed. One thus needs to choose the architecture $\widehat{\mathbf{R}}$
which assigns highest weight to the observed dataset $\mathfrak{f}$ or in
other words whose penalized empirical entropy difference maximizes the density
$\rho$ in the induced equivalence class
$[\mathfrak{f}]_{\widehat{\mathbf{R}}}$, conversely minimizing the bic score.

#### The expectation of entropy and AIC score.

Another way to characterize an equivalence class and assess the quality of
modeling is via expectations over member distributions
$\mathfrak{p}\in[\mathfrak{f}]_{\widehat{\mathbf{R}}}$. For instance, the
expected entropy in $[\mathfrak{f}]_{\widehat{\mathbf{R}}}$ equals

$\displaystyle N\langle H[\mathfrak{p}]\rangle=NH[\hat{\mathfrak{p}}]-\langle
N(H[\hat{\mathfrak{p}}]-H[\mathfrak{p}])\rangle=NH[\hat{\mathfrak{p}}]-\frac{{\left|\mathcal{A}\right|}-\widehat{D}}{2}~{},$
(38)

using the entropic radius from Eq. (19) and taking the expectation over chi-
squared density. Simplicity (lower dimension $\widehat{D}$) seems to dictate
the expected entropy to lie further away from maxent, cf. chi-squared
probability density spreading out with increasing degrees of freedom. On the
other hand, goodness of fit translates into the expected entropy approaching
the empirical one. Combining both heuristic objectives we require the expected
entropy to be closer to empirical $H[\mathfrak{f}]$ than to maxent. In other
words, we could look for the minimum of the difference of entropy differences

$\displaystyle N\left(\langle
H[\mathfrak{p}]\rangle-H[\mathfrak{f}]\right)-N\left(H[\hat{\mathfrak{p}}]-\langle
H[\mathfrak{p}]\rangle\right)=ND_{\textsc{kl}}(\mathfrak{f}\;\|\;\hat{\mathfrak{p}})-({\left|\mathcal{A}\right|}-\widehat{D})=\mbox{$\frac{1}{2}$}\widehat{\textsc{aic}}-NH[\mathfrak{f}]-{\left|\mathcal{A}\right|}~{}.$
(39)

In the spirit of Eq. (6), minimization of aic score (Akaike information
criterion [18]) can be understood as selecting the set of constraints whose
maxent-to-expected-entropy ratio of combinatorial possibilities is larger than
the expected-to-empirical-entropy ratio.

##  5 Training and test metrics

For the purposes of this section, we suppose that some form of “truth”
asymptotically exists. Let $\mathfrak{q}\in\mathcal{P}$ describe the
asymptotic model as some maxent distribution of architecture $\mathbf{R}$. To
assess the quality of modeling at a given sample size $N$ we first define the
expected training error of an architecture $\widehat{\mathbf{R}}$ by

$\displaystyle\varepsilon_{\text{train}}=\sum_{\mathfrak{f}\in\mathcal{P}}\text{mult}(N\mathfrak{f};\mathfrak{q})\,ND_{\textsc{kl}}(\mathfrak{q}\;\|\;\hat{\mathfrak{p}})~{},$
(40)

where $\hat{\mathfrak{p}}\in[\mathfrak{f}]_{\widehat{\mathbf{R}}}$ is
understood to be the maxent distribution in the equivalence class induced by
$\widehat{\mathbf{R}}$ on samples $\mathfrak{f}$ drawn from $\mathfrak{q}$.
This error is governed by the expected Kullback-Leibler (kl) divergence [19]
of model from sample distribution for samples drawn from the “true”
distribution $\mathfrak{q}$. Intuitively, the training error quantifies the
average proximity to $\mathfrak{q}$ achieved by models $\hat{\mathfrak{p}}$
that were trained on a sample of size $N$. It depends on $N$, the training
architecture and the unknown asymptotic model. To evaluate Eq. (40) we
consider two distinct cases.

If model architecture $\widehat{\mathbf{R}}$ implies all asymptotic
constraints $\mathbf{R}$, then there must exist parameter vector
$\hat{\boldsymbol{\theta}}\in\mathbb{R}^{\hat{D}}$ such that

$\displaystyle\mathfrak{q}_{\alpha}=\exp\left(\sum_{a=1}^{\widehat{D}}\hat{\theta}_{a}\widehat{R}_{a\alpha}\right)~{}.$
(41)

As data varies, linearized fluctuations of model distribution
$\hat{\mathfrak{p}}$ around $\mathfrak{q}$,

$\displaystyle\hat{\mathfrak{p}}_{\alpha}=\mathfrak{q}_{\alpha}\left(1+\frac{1}{\sqrt{N}}\sum_{a=1}^{\widehat{D}}y_{a}\widehat{R}_{a\alpha}\right)$
(42)

(column-wise multiplication by $\mathfrak{q}_{\alpha}\neq 0$ is a rank-
preserving operation) induce a parametrization of some generic distribution
$\mathfrak{f}$ around the asymptotic distribution $\mathfrak{q}$,

$\displaystyle\mathfrak{f}_{\alpha}=\hat{\mathfrak{p}}_{\alpha}+\frac{1}{\sqrt{N}}\sum_{i=1}^{{\left|\mathcal{A}\right|}-\widehat{D}}x_{i}\widehat{X}_{i\alpha}=\mathfrak{q}_{\alpha}+\frac{1}{\sqrt{N}}\left[\sum_{a=1}^{\widehat{D}}y_{a}\left(\widehat{R}_{a\alpha}\mathfrak{q}_{\alpha}\right)+\sum_{i=1}^{{\left|\mathcal{A}\right|}-\widehat{D}}x_{i}\widehat{X}_{i\alpha}\right]~{},$
(43)

where $\widehat{\mathbf{X}}$ orthogonally encodes the kernel of
$\widehat{\mathbf{R}}$, similarly to Eq. (17). Note in particular that
normalization condition on the simplex $\mathcal{P}$ realized by $\bar{a}$-th
row in the coefficient matrix with $\widehat{C}_{\bar{a}\alpha}=1$
automatically ensures

$\displaystyle\sum_{\alpha=1}^{\left|\mathcal{A}\right|}X_{i\alpha}=\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\widehat{C}_{\bar{a}\alpha}\widehat{X}_{i\alpha}=\sum_{a=1}^{\widehat{D}}(\mathbf{T}^{-1})_{\bar{a}a}\underbrace{\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\widehat{R}_{a\alpha}\widehat{X}_{i\alpha}}_{=0}=0\quad\forall\,i=1,\ldots{\left|\mathcal{A}\right|}-\widehat{D}~{},$
(44)

where $\mathbf{T}$ represents the Gaussian-Jordan transformation matrix from
$\widehat{\mathbf{C}}$ to $\widehat{\mathbf{R}}$. Enforcing thus
$\mathfrak{f}\in\mathcal{P}$ requires that the components of
$y\in\mathbb{R}^{\widehat{D}}$ are linearly dependent,

$\displaystyle\sum_{a=1}^{\widehat{D}}y_{a}\underbrace{\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\widehat{R}_{a\alpha}\mathfrak{q}_{\alpha}}_{\text{``true''
$\widehat{\mathfrak{m}}$}}=0$ (45)

by vanishing projection on the moments of $\widehat{\mathbf{R}}$ induced on
the “true” distribution. Hence, there remain $\widehat{D}-1$ directions along
which model distribution $\hat{\mathfrak{p}}$ can vary.

At large $N$, expansion Eq. (34) states that multinomial distribution is
governed by the “distance” of asymptotic $\mathfrak{q}$ to empirical
distribution $\mathfrak{f}$. Upon using Eq. (43) the information-theoretic
approximation to multinomial distribution gives a leading quadratic (in the
fluctuations) term,

$\displaystyle-\log\text{mult}(N\mathfrak{f};\mathfrak{q})=\mbox{$\frac{1}{2}$}\sum_{a,b=1}^{\widehat{D}}y_{a}y_{b}\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\widehat{R}_{a\alpha}\mathfrak{q}_{\alpha}\widehat{R}_{b\alpha}+\mbox{$\frac{1}{2}$}\sum_{i,j=1}^{{\left|\mathcal{A}\right|}-\widehat{D}}x_{i}x_{j}\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\widehat{X}_{i\alpha}\frac{1}{\mathfrak{q}_{\alpha}}\widehat{X}_{j\alpha}+\mathcal{O}\left(\log
N\right)~{}.$ (46)

As argued in Eq. (25), precision matrices of the above form are invertible, as
long as $\mathbf{R}$ and $\mathbf{X}$ are of full row-rank. Any linear term in
the fluctuations vanishes identically due to normalization condition Eq. (44)
or drops after imposing Eq. (45), whereas mixed terms in $x$ and $y$ vanish
due to the normality of $\widehat{\mathbf{X}}$ space to
$\widehat{\mathbf{R}}$. This large-$N$ expansion of multinomial distribution
shows that it asymptotically goes into a multivariate normal distribution,
thus enabling us to nicely write the integral representation of Eq. (40),
which takes the form of an $({\left|\mathcal{A}\right|}-1)$-dimensional
integral in $x$ and $y$. In particular, the Gaussian kernel associated to Eq.
(46) can be always diagonalized by an orthogonal transformation in the
fluctuations. At the same time, the integrand becomes

$\displaystyle
ND_{\textsc{kl}}(\mathfrak{q}\;\|\;\hat{\mathfrak{p}})=\mbox{$\frac{1}{2}$}\sum_{a,b=1}^{\widehat{D}}y_{a}y_{b}\sum_{\alpha=1}^{\left|\mathcal{A}\right|}\widehat{R}_{a\alpha}\mathfrak{q}_{\alpha}\widehat{R}_{b\alpha}+\mathcal{O}\left(N^{-1/2}\right)~{},$
(47)

rendering straight-forward the computation of

$\displaystyle\varepsilon_{\text{train}}=\frac{\widehat{D}-1}{2}+\mathcal{O}\left(N^{-1/2}\right)~{}.$
(48)

Whenever our model fully captures the asymptotic architecture, the expected
error converges as $N\rightarrow\infty$. It is optimal precisely when
$\widehat{\mathbf{R}}=\mathbf{R}$.

On the other hand, if $\widehat{\mathbf{R}}$ misses any of the constraints
from asymptotic architecture $\mathbf{R}$, the large-$N$ localization of the
integral representation of Eq. (40) at the saddle point would give

$\displaystyle\varepsilon_{\text{train}}=ND_{\textsc{kl}}(\mathfrak{q}\;\|\;\hat{\mathfrak{p}})+\mathcal{O}\left(\sqrt{N}\right)~{},$
(49)

where $\hat{\mathfrak{p}}$ is the model distribution induced by
$\widehat{\mathbf{R}}$ on the asymptotic counts $N\mathfrak{q}$. Since the kl
divergence would be finite – model architecture lacking some constraints – the
training error diverges linearly in $N$ in this case.

Similarly, one can estimate the expected test (or generalization) error of a
model architecture $\widehat{\mathbf{R}}$,

$\displaystyle\varepsilon_{\text{test}}=\sum_{\mathfrak{f},\mathfrak{g}\in\mathcal{P}}\text{mult}(N\mathfrak{g};\mathfrak{q})\text{mult}(N\mathfrak{f};\mathfrak{q})ND_{\textsc{kl}}(\mathfrak{g}\;\|\;\hat{\mathfrak{p}})~{},$
(50)

where $\hat{\mathfrak{p}}\in[\mathfrak{f}]_{\widehat{\mathbf{R}}}$ represents
the model distribution, as before. Parameterizing fluctuations around the
asymptotic distribution $\mathfrak{q}$ along the lines of Eq. (42) and Eq.
(43), we obtain

$\displaystyle\varepsilon_{\text{test}}=\frac{{\left|\mathcal{A}\right|}+\widehat{D}-2}{2}+\mathcal{O}\left(N^{-1/2}\right)~{},$
(51)

if test model $\widehat{\mathbf{R}}$ implies the asymptotic architecture
$\mathbf{R}$ and an estimate that diverges as $N\rightarrow\infty$, whenever
any of the asymptotic constraints is missed by $\widehat{\mathbf{R}}$.

##  6 Model selection in the inverse Ising problem

A classical paradigm used both in benchmark tests and to illustrate various
properties of magnetic materials is the Ising model. Suppose that a group of
Ising spins $\sigma_{i}\in\\{0,1\\}$ in lattice-gas gauge labeled by site
index $i=1,\ldots L$, which could have been e.g. extracted from a much bigger
Ising crystal, had energy (Hamiltonian)

$\displaystyle
E(\vec{\sigma})=\sum_{i=1}^{L}h_{i}\sigma_{i}+\mbox{$\frac{1}{2}$}\sum_{i,j=1}^{L}J_{ij}\sigma_{i}\sigma_{j}+\mbox{$\frac{1}{6}$}\sum_{i,j,k=1}^{L}T_{ijk}\sigma_{i}\sigma_{j}\sigma_{k}~{}.$
(52)

As we marginalize over lattice regions, there is nothing that prevents us from
writing effective associations of higher orders among the remaining spins, but
for concreteness we consider up to cubic interactions. Furthermore, we take
$L=5$ to speed up evaluation of multiple instances tested in parallel. For our
benchmark tests below, we shall randomly generate different values for bias
vector $\mathbf{h}$, symmetric coupling matrix $\mathbf{J}$ and cubic tensor
$\mathbf{T}$ (excluding self-interactions which are presently meaningless)
based on the hypergraph

$\displaystyle\mathcal{G}_{\text{Ising}}=\\{\\{1,2,3\\},\\{1,2,4\\},\\{3,5\\},\\{4,5\\}\\}~{}.$
(53)

For simple Ising systems, any $n$-point feature association implies all lower
terms. For example, a trilinear association $\\{i,j,k\\}$ automatically
induces all pairwise associations among spins $\\{i,j\\}$, $\\{i,k\\}$ and
$\\{j,k\\}$ as well as the one-site biases for $i$, $j$ and $k$. Hence, the
Ising system has $5$ non-vanishing biases, $7$ couplings and 2 cubic
interactions. The Boltzmann distribution

$\displaystyle\mathfrak{q}(\vec{\sigma})=\left(\sum_{\sigma_{1}^{\prime}=0,1}\cdots\sum_{\sigma^{\prime}_{L}=0,1}e^{E(\vec{\sigma}^{\prime})}\right)^{-1}e^{E(\vec{\sigma})}$
(54)

of statistical system Eq. (52) defines an asymptotic model with333In the
parametric formulation of thermodynamics, the normalization condition on
$\mathcal{P}$ is realized by the partition function. $D=15$ in a spin space of
${\left|\mathcal{A}\right|}=2^{5}=32$ microstates.

From this Ising system in the inverse problem, a modeler obtains samples at
different sizes $N$. Knowing only that this is a system of potentially
interacting spins leaves $7\,580$ candidate hypergraphs with distinct
architectures $\widehat{\mathbf{R}}$ reflecting different interaction
structures. For each architecture, ipf or any Newton-based method can quickly
deduce the corresponding maxent distribution on a given training sample, as
outlined in Section 3.1. In such setting, applications are usually concerned
with two intimately related aspects: (i) the reconstruction of the original
hypergraph $\mathcal{G}_{\text{Ising}}$ from a sample of size $N$ and (ii) the
generalization error of a trained model on different samples drawn from the
same system.

To answer those questions by better understanding the various information-
theoretic metrics discussed in Section 4, we have independently generated 10
000 samples from $1\,000$ different realizations of the Ising system Eq. (52)
with interaction structure Eq. (53) at various sample sizes. By a realization,
we simply mean a set of fields and couplings parametrizing the Ising
Hamiltonian, which were independently drawn from a normal distribution of zero
mean and $\mathcal{O}\left(1\right)$ variance. On each training dataset of
size $N$ the maxent distribution associated to all $7\,580$ architectures is
deduced. Subsequently, the architecture $\widehat{\mathbf{R}}$ is recorded
each time which exhibits (a) lowest bic, (b) lowest aic, (c) lowest
$\widehat{D}$ and empirical $p$-value above threshold (hyper-maxent), (d)
lowest $\widehat{D}$ and typical for all more complex models implying
$\widehat{\mathbf{R}}$ (hyper-maxent + lrt).

In order to determine a threshold for (c) and (d) below which a $p$-value is
considered to be significant (see Section 4.1) we exploit two facts. According
to Eq. (48) and Eq. (51), expected errors tend to decrease for simpler models,
hence the cutoff should increase for higher-dimensional kernels Eq. (17) for
which we anticipate on average smaller atypical fluctuations. On the other
hand, larger datasets tend to eliminate sampling noise sharpening the
distribution of $p$-values. The empirical $p$-value associated to any
architecture that misses some asymptotic constraint goes faster to zero, while
any slightly atypical dataset sampled from $\mathfrak{q}$ could lead to
smaller $p$-values under the asymptotic model. The acceptance threshold should
thus scale with inverse sampling size to ensure a proper asymptotic limit.
Combining these observations we conclude that the fraction

$\displaystyle\alpha_{\text{empirical}}=\frac{{\left|\mathcal{A}\right|}-\widehat{D}}{N}$
(55)

has a scaling which seems reasonable for a cut-off in the maxent hypothesis
testing. In a similar spirit, we use in the lrt part of (d) as significance
level

$\displaystyle\alpha_{\textsc{lrt}}=\frac{2{\left|\mathcal{A}\right|}-\widehat{D}-\widehat{D}^{\prime}}{N}$
(56)

to decide about the $p$-value when comparing two nested model architectures
$\widehat{\mathbf{R}}$ and $\widehat{\mathbf{R}}^{\prime}$.

Figure 1: Left: Average “distance” from the asymptotic distribution associated
to $\mathcal{G}_{\text{Ising}}$ of model distribution that is selected by the
various scores as a function of sample size $N$. Below: Average reconstruction
accuracy determined by the number of correctly identified hypergraphs
$\mathcal{G}_{\text{Ising}}$ in the generated samples. Right: Average
generalization error (and its dispersion) of model distribution on unseen data
sampled from the asymptotic system as a function of $N$. Figure 2:
$N$-dependent average rates of true/false positive associations against
asymptotic graph $\mathcal{G}_{\text{Ising}}$ characterizing the various
selection procedures.

In Figure 1, we plot the training error Eq. (40) averaged over all sampled
datasets of various sizes $N=10^{2},10^{3},\ldots,10^{11}$. To get a better
idea on the training quality the average graph reconstruction that was
achieved at each sample size is comparatively plotted. This shows how often
model architecture $\widehat{\mathbf{R}}$ selected by the various metrics
precisely captures the underlying Ising model associated to graph Eq. (53),
but nothing more. On the right plot, the test error Eq. (50) on $1\,000$
samples is given again averaged over the $10^{7}$ training datasets at each
sample size. To obtain a feeling about the dispersion of the generalization
error we alongside plot the difference of the test error of the other scores
from the consistently best-scoring lrt model within each sample.

We observe that hyper-maxent and lrt with the acceptance threshold introduced
in Eq. (55) and Eq. (56) as well as the bic score are all scoring metrics that
are consistent with large-$N$ expansion(s). As $N$ increases, they clearly
exhibit the correct asymptotic behavior eventually achieving $100\%$
reconstruction accuracy. At the same time, the generalization error saturates
the bound in Eq. (50) for $\widehat{D}=D$. It is important to stress that
modifying the significance level, does not change the asymptotic behavior, as
long as the scaling goes (at least) as $1/N$. This demonstrates the
universality of large-$N$ expansion which renders any compatible selection
algorithms indistinguishable as $N$ tends to infinity. On the other hand, aic
consistently appears sub-optimal failing to recover the underlying model.
Instead, the reasoning of maximizing the expected entropy Eq. (38) seems to
capture models that are generically more complex than necessary and are thus
plagued by false positive (fp) associations.

To get a handle on true/false positive (tp/fp) rates, we plot444A Python
script to randomly generate Ising systems of a given architecture, sample from
them and benchmark the presented model-selection procedures can be found under
github.com/imbbLab/ModelSelection. the percentages of correctly identified
interactions, which are achieved by the various selection procedures, averaged
over all sampled datasets at each size $N$. As a base reference, we use the
probability that a randomly selected constraint out of the 31 possible is a
true (false) positive, i.e. falls (does not fall) within the set of the 14
asymptotic constraints. This roc-curve inspired plot verifies our previous
finding on the large-$N$ compatibility of all three selection procedures
except for the aic score.

However, it is noteworthy that aic manages to quickly reach a higher tp rate
at the expense of keeping on average a finite fp rate, which in the present
setting remains lower than the base rate (dotted gray line). By modifying
acceptance threshold Eq. (56), one can adjust the fp tolerance of lrt
selection in the small-$N$ regime. Most interestingly, we recognize that bic
and the hyper-maxent consistently have (almost) vanishing fp rates, even at
sample sizes of order $100$. In particular, this gives a strong hint in favor
of the $\mathcal{O}\left(1/N\right)$ scaling of acceptance threshold Eq. (55).
To further support the suggested scaling for the threshold, we mention the
observation that the architecture $\widehat{\mathbf{R}}=\mathbf{R}$
corresponding to the asymptotic graph Eq. (53) always has an empirical
$p$-value Eq. (28) above $\alpha_{\text{empirical}}$ at any sample size $N$.

## References

  * [1] C. E. Shannon “A mathematical theory of communication” The Bell System Technical Journal 27 (1948) no. 3, 379–423.
  * [2] E. T. Jaynes Probability theory: The logic of science.  Cambridge university press 2003.
  * [3] J. Shore and R. Johnson “Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy” IEEE Transactions on Information Theory 26 (1980) no. 1, 26–37.
  * [4] J. B. Paris and A. Vencovská “A note on the inevitability of maximum entropy” International Journal of Approximate Reasoning 4 (1990) no. 3, 183–223.
  * [5] I. Csiszár “Maxent, mathematics, and information theory” in Maximum entropy and Bayesian methods pp. 35–50.  Springer 1996.
  * [6] J. M. Keynes “Chapter iv: The principle of indifference” A treatise on probability 4 (1921) 41–64.
  * [7] O. Loukas and H. R. Chung “Categorical distributions of maximum entropy under marginal constraints” arXiv preprint arXiv:2204.03406 (2022).
  * [8] E. T. Jaynes “Prior probabilities” IEEE Transactions on systems science and cybernetics 4 (1968) no. 3, 227–241.
  * [9] I. Csiszár “I-divergence geometry of probability distributions and minimization problems” The annals of probability (1975) 146–158.
  * [10] S.-i. Amari and H. Nagaoka Methods of information geometry vol. 191.  American Mathematical Soc. 2000.
  * [11] C. Geyer and G. Meeden “Asymptotics for constrained dirichlet distributions” Bayesian Analysis 8 (2013) no. 1, 89–110.
  * [12] E. T. Jaynes Concentration of Distributions at Entropy Maxima (1979) pp. 315–336.  Springer Netherlands Dordrecht 1989.
  * [13] C. T. Ireland and S. Kullback “Contingency tables with given marginals” Biometrika 55 (1968) no. 1, 179–188.
  * [14] Y. M. Bishop, S. E. Fienberg, and P. W. Holland Discrete multivariate analysis: theory and practice.  Springer Science & Business Media 2007.
  * [15] S. Silvey Statistical Inference.  Chapman & Hall/CRC Monographs on Statistics & Applied Probability. Taylor & Francis 1975.
  * [16] A. Buse “The likelihood ratio, wald, and lagrange multiplier tests: An expository note” The American Statistician 36 (1982) no. 3a, 153–157 [arXiv:https://doi.org/10.1080/00031305.1982.10482817].
  * [17] G. Schwarz “Estimating the Dimension of a Model” The Annals of Statistics 6 (1978) no. 2, 461 – 464.
  * [18] H. Akaike “A new look at the statistical model identification” IEEE transactions on automatic control 19 (1974) no. 6, 716–723.
  * [19] S. Kullback Information theory and statistics.  Courier Corporation 1997.