# Maximal Relevance and Optimal Learning Machines

O Duranthon¹, M Marsili², and R Xie³

¹Department of Physics, École Normale Supérieure, 24 rue Lhomond, 75005 Paris, France

²The Abdus Salam International Centre for Theoretical Physics, Strada Costiera 11, 34151 Trieste, Italy

³Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics, Central China Normal University (CCNU), Wuhan, China

## Abstract

We explore the hypothesis that learning machines extract representations of maximal relevance, where the relevance is defined as the entropy of the energy distribution of the internal representation. We show that the mutual information between the internal representation of a learning machine and the features that it extracts from the data is bounded from below by the relevance. This motivates our study of models with maximal relevance – that we call Optimal Learning Machines – as candidates of maximally informative representations. We analyse how the maximisation of the relevance is constrained both by the architecture of the model used and by the available data, in practical cases. We find that sub-extensive features that do not affect the thermodynamics of the model, may affect significantly learning performance, and that criticality enhances learning performance, but the existence of a critical point is not a necessary condition. On specific learning tasks, we find that *i*) the maximal values of the likelihood are achieved by models with maximal relevance, *ii*) internal representations approach the maximal relevance that can be achieved in a finite dataset and *iii*) learning is associated with a broadening of the spectrum of energy levels of the internal representation, in agreement with the maximum relevance hypothesis.

## 1 Introduction

Statistical mechanics models have been used in statistical learning since the pioneering works on associative memory [1] and on Boltzmann Learning Machines [2]. Many attempts to make sense of the spectacular performance of learning machines such as deep neural networks, have focused on understanding their statistical mechanics properties, as compared e.g. to spin glasses [3]. Several attempts have been made to understand the typical properties of learning machines [4, 5, 6, 7]. These properties emerge after the interaction parameters have been learned from highly structured datasets [8, 9, 10], and they depend on what measure of error or likelihood is used in supervised or unsupervised learning, respectively. This makes it hard to discuss the typical properties of learned models without reference to the particular dataset on which the model is trained, the objective function assumed, the algorithms or the architecture used. This suggests that, in order to understand learning, we should focus on the ingredients which makes it possible [11].
This paper takes a different perspective, and it explores the claim of Refs. [12, 13] that learning maximises the *relevance* of the representation that the machine extracts from the data. Let us briefly recall the rationale behind this conclusion, which is articulated in Ref. [13], to which we refer for more details. In order to define the relevance, let us consider a generic unsupervised learning task that maps the data to the internal states of a learning machine, in order to capture the structure of statistical dependencies of the data. We assume that the data *x* is generated from an unknown distribution $\mathcal{P}(\mathbf{x})$, that we shall call the *generative model*. A learning machine maps the data *x* to the internal states *s* of the hidden layer(s) of the learning machine. Training with data drawn from $\mathcal{P}(\mathbf{x})$ induces a distribution $p(\mathbf{s})$ on the hidden states. We define the energy $E_{\mathbf{s}} = -\log p(\mathbf{s})$ as the coding cost associated to a microscopic state $\mathbf{s}$. The average energy

$$ \langle E \rangle = - \sum_{\mathbf{s}} p(\mathbf{s}) \log p(\mathbf{s}) \equiv H[\mathbf{s}] \quad (1) $$

coincides with the entropy of $\mathbf{s}$, and it is a measure of the *resolution* of the representation. Ref. [13] argues that the entropy of the energy distribution

$$ H[E] = - \sum_{E} p(E) \log p(E), \quad p(E) = \sum_{\mathbf{s}} p(\mathbf{s}) \delta(E - E_{\mathbf{s}}) \quad (2) $$

can be taken as a quantitative measure of the information that the model contains on the structure of the data, or of its generative model. In order to motivate this conclusion, we observe that $H[E]$ is typically small for Gibbs distributions, which encode states of maximal ignorance, because they have a very narrow distribution of energies. The internal representations $p(\mathbf{s})$ of a learning machine should be as far as possible from a maximum entropy model. Hence it should feature a broad spectrum of energies $E_{\mathbf{s}}$, that corresponds to a large value of $H[E]$. In this representation, two states $\mathbf{s}$ and $\mathbf{s}'$ with very different energies $E_{\mathbf{s}}$ and $E_{\mathbf{s}'}$ necessarily code for significantly different patterns in the data, because they cannot occur typically, within the same maximum entropy model.

Following Cubero et al. [13], we shall henceforth call $H[E]$ *resolution* and $H[E]$ *relevance*. The above discussion leads to the hypothesis that machines that are trained on highly structured data should approach the limit where $H[E]$ is as large as possible. In order to explore this claim, we define the class of Optimal Learning Machines (OLM) as the ideal limit of those statistical models that maximise the relevance at a given resolution $H[\mathbf{s}]$

$$ \max_{\{E_{\mathbf{s}}\}:H[\mathbf{s}]=\bar{E}} H[E] \quad (3) $$

over the distribution of energy levels $E_{\mathbf{s}}$. The optimisation in Eq. (3) is constrained, in practical cases, by the architecture of the machine. When these constraints are neglected, machines that satisfy Eq. (3) have a broad distribution of energies with an exponential density of states [13]. This, in turn, is equivalent to statistical criticality [14, 13], which is a widely observed statistical regularity in efficient representations [7, 12, 15]. This corroborates the conjecture that the principle of maximal relevance applies to real learning machines.

The contributions of this paper is twofold: first we show that, in a generic learning task, the relevance lower bounds the mutual information between the internal representation and the features that the machine extracts from the data. This provides a rationale for why learning machines should approach the limit of maximal relevance, in order to extract the maximal amount of information from the data.
Secondly, we investigate how the optimisation in Eq. (3) is constrained both by the architecture of the machine used and by the amount of data on which it is trained. We start exploring the relevance of different architectures. Gaussian learning machines [16] offer an interesting counter-example, where the relevance $H[E]$ does not vary during training. This is consistent with the fact that the learned distribution $p(s)$ remains a Gaussian, irrespective of the generative model $\mathcal{G}$. Hence, Gaussian learning machines do not learn anything on the generative model $\mathcal{G}$, beyond its parameters.

Taking $H[E]$ as a quantitative measure of learning performance, we study a class of toy learning machines with discrete energy levels, in order to characterise the architectural features associated with a high value of $H[E]$. Our results suggest that features that confer superior learning performance may be associated to sub-extensive features that are not accessible to standard statistical mechanics approaches. A superior learning performance is not necessarily related to the existence of a critical point separating two different phases, but when a critical point exists, learning performance improves when the model is tuned to its critical point.

Finally, we present a series of numerical experiments, in order to understand how closely real learning machines approximate the ideal limit of OLM when the dataset is finite. Within the class of toy learning machines with discrete energy levels discussed earlier, we show that the models that have superior learning performance in two distinct (unsupervised learning) tasks are those with the highest values of $H[E]$. Then we move to Restricted Boltzmann Machines (RBM). There we argue that the number of samples on which the machine is trained limits the resolution of the energy scale. When the resolution and relevance are evaluated at the scale dictated by the available data, our numerical experiments show that RBM converge to states of maximal relevance, in agreement with [12]. For small systems, it is possible to compute the full spectrum of energy levels and to follow its evolution. We show that during learning the distribution of energy levels indeed broadens. We summarise and comment our results in the concluding section.

## 2 Relevance and hidden features

Let us assume that data points $\mathbf{x} = (x_1, ..., x_d)$ are generated as draws from an unknown distribution<fnref n="1" /> $\mathcal{P}(\mathbf{x})$ that is characterised by a rich structure of dependence between the components $x_a$. For example, the MNIST dataset [17] of hand written digits has $d = 28 \times 28 = 784$ pixels, for each of which, the grey scale is codified in an integer $x_a$ in the range [0, 255]. In spite of the fact that $\mathbf{x}$ belongs to a very high-dimensional space, the MNIST dataset spans a manifold of only $d_{int} \approx 13$ intrinsic dimensions [18]. This implies that $\mathcal{P}(\mathbf{x})$ should be a model of strongly interacting variables, and that there should be a set of coordinates $\mathcal{Q}(\mathbf{x}) = (\mathcal{Q}_1(\mathbf{x}), ..., \mathcal{Q}_{d_{int}}(\mathbf{x}))$ that describes the relevant variation in the dataset. We refer to $\mathcal{Q}$ as the *hidden features*. As for $\mathcal{P}$, the backslash indicates that $\mathcal{Q}$ is an unknown theoretical construct that represents statistically significant patterns in the dataset<fnref n="2" />.

We focus on stochastic models, such as Restricted Boltzmann Machines (RBM) or Deep Belief Networks (DBN), with a finite number of discrete internal states. The goal of learning is to find a

<footnote n="1">¹Here and in the rest of the paper, backslashed symbols (e.g. $\mathcal{P}$) refer to unknown entities.</footnote>

<footnote n="2">²In the case of the MNIST dataset, a draw from the generative model $\mathcal{P}$ is a theoretical abstraction for the process of hand writing a digit by a human and features are not necessarily related to “meaningful” patterns such as the values of digits. Deep learning machines also rely on features that look like noisy patterns, such as those used in adversarial examples. As shown by Ilyas et al. [19], machines that only use features which are reasonable for humans perform worse, in terms of generalisation, than those that use all features.</footnote>
statistical model $p(s)$ over a discrete variable $s$, and a mapping $p(x|s)$ from $s$ to $x$, such that the generating distribution

$$p(x) = \sum_{s} p(x|s)p(s) \qquad (4)$$

is as close as possible to the empirical one (in unsupervised learning). The objective function employed may differ, yet in all cases, training aims at approximating the unknown generative process $p(x) \approx \mathcal{P}(x)$ as closely as possible<fnref n="3" />. Here we focus on the outcome of this process, abstracting from algorithmic details on how the representation is learned. We remark that, unlike $\mathcal{P}$ which is a theoretical abstraction, $p(s)$ is a proper statistical model, as well as the induced distribution $p(x)$ in Eq. (4). We define the energy of state $s$ as

$$E_s = -\log p(s) \qquad (5)$$

which is the cost (in nats) for coding state $s$.

<figure><img src="image_1.png" /><figcaption>Figure 1: (Top) Sketch of the structure of statistical dependencies of the data in terms of hidden features $\mathcal{X}$. Points with similar features $\mathcal{X} \approx \bar{\mathcal{X}}$ are statistically identical. Statistically distinct points are distinguished by a significant variation of $\mathcal{X}$. (Bottom) Learning aims at representing this structure in the internal representation $s$ of a learning machine, in terms of extracted features $\phi$ that approximately reproduce the hidden features ($\phi \approx \mathcal{X}$).</figcaption></figure>

If the learning machine captures correctly the structure of the data $x$ with which it has been trained, then it must extract features $\phi$ that approximate well the hidden features $\mathcal{X}$. Hence there

<footnote n="3">³In supervised learning the representation aims at reproducing a functional relation $\underline{x}_{\text{out}} = f(\underline{x}_{\text{in}})$ between two parts of the data $x = (\underline{x}_{\text{in}}, \underline{x}_{\text{out}})$ as well as possible. This is similar to unsupervised learning with $\mathcal{P}(x) = \mathcal{P}(\underline{x}_{\text{in}})\delta(\underline{x}_{\text{out}} - f(\underline{x}_{\text{in}}))$.</footnote>
must be a function $\phi(s)$ defined on the states of the hidden layers – the extracted features – such that its projection

$$ \langle \phi | x \rangle \equiv \sum_s \phi(s) p(s | x) \qquad (6) $$

is close to the hidden features $\mathcal{N}(x)$. We call $\phi$ the extracted features, to contrast them with the hidden ones $\mathcal{N}$.

In an ideal situation, the distribution $p(s)$ should contain no further information than that contained in $\phi(s)$. In other words, $\phi$ should be the sufficient statistics for the distribution $p(s)$. This means that the distribution of $s$ on all states with $\phi(s) = \bar{\phi}$ must be a maximum entropy distribution, i.e. $p(s|\bar{\phi}) = A(\bar{\phi})\delta_{\phi(s),\bar{\phi}}$, where $A(\bar{\phi})$ is the inverse of the number of states $s$ with $\phi(s) = \bar{\phi}$. Put differently, states with the same features $\phi(s) = \phi(s')$ should be statistically indistinguishable, i.e. $p(s) = p(s')$ for all $s$ and $s'$ such that $\phi(s) = \phi(s')$. Therefore $E_s = E_{s'}$ for any two states $s$ and $s'$ such that $\phi(s) = \phi(s')$, which means that $E_s$ must be a function of $\phi(s)$. In other words, the statistical dependence of the variables can be represented as the Markov chain$^4$

$$ s \to \phi \to E. \qquad (7) $$

The data processing inequality [20] then implies that

$$ I(s, \phi) \ge I(s, E) = H[E], \qquad (8) $$

where $I(X, Y)$ is the mutual information between variables $X$ and $Y$ [20], and the last equality comes from the fact that $I(s, E) = H[E] - H[E|s]$ and $H[E|s] = 0$. Therefore, $H[E]$ provides a lower bound to the information $I(s, \phi)$ that the representation contains on the extracted features. Notice that $\phi$ is a function of $s$, and hence $I(s, \phi) = H[\phi]$.

Going back to the hidden features $\mathcal{N}$, we observe that the training process can be represented by the Markov chain

$$ \mathcal{N} \to x \to s \to \phi \qquad (9) $$

where $s$ is the representation that unfolds from training. If $\phi$ are sufficient statistics, which means that $x$ and $s$ are independent conditional on $\phi$, then the last step of the chain can be inverted, i.e. $\mathcal{N} \to x \to \phi \to s$. Then the data processing inequality implies that $I(s, \mathcal{N}) \le I(s, \phi) = H[\phi]$. Under the assumption that the variation $H[\phi]$ of the extracted features is entirely induced by the hidden features, it is reasonable to expect that this inequality is tight$^5$

The inequality (8) is the main result of this Section. It connects the mutual information $I(s, \phi)$ of the extracted features with the relevance $H[E]$. Apart from simple cases, such as the Gaussian learning machine discussed below, it is not clear how to identify the sufficient statistics of a distribution $p(s)$. So the extracted feature $\phi$ are hardly accessible in practice, let alone the hidden ones $\mathcal{N}$. The distribution $p(s)$ instead can be computed for learning machines trained on structured data, and hence $H[E]$ can be estimated, as we shall see in the following Sections. The inequality (8) suggests that statistical models with a high value of $H[E]$ are natural candidates for good learning machines, even though there is no guarantee that the inequality (8) is tight.

<footnote>$^4$The notation $z \to y \to w$ implies that, if $y$ is known, $w$ and $z$ are independent, i.e. that $I(z, w|y) = 0$ [20]. This holds trivially true if $w$ is a function of $y$, because then $w$ is not a random variable, when $y$ is known.</footnote>

<footnote>$^5$If the last step of the Markov chain (9) can be reversed, then $I(s, \mathcal{N}) = I(\phi, \mathcal{N}) = H[\phi] - H[\phi|\mathcal{N}] = I(s, \phi) - H[\phi|\mathcal{N}]$. Therefore $I(s, \mathcal{N}) \approx I(s, \phi)$ is equivalent to the statement that the extracted features are approximately a function of the hidden features, i.e. $H[\phi|\mathcal{N}] \approx 0$.</footnote>
The comparison between different learning machines on the basis of $H[E]$ makes sense only if the representations that they generate have the same resolution $H[s]$, i.e. the same coding cost. For this reason, we shall consider the relevance $H[E]$ as a measure of the learning capacity, at a given level of resolution $H[s]$. An absolute maximum value for the maximum of $H[E]$ for a given $H[s]$ can be obtained by studying the problem

$$ \max_{\{W(E)\}: H[s] = E} H[E] \qquad (10) $$

where $W(E)$ is the number of energy levels with $E_s = E$, which is constrained by the two conditions

$$ H[s] = -\sum_{s} p(s) \log p(s) = \sum_{E} W(E)e^{-E}E \qquad (11) $$

$$ N = \sum_{s} 1 = \sum_{E} W(E) \qquad (12) $$

where $N$ is the number of available states. This problem is studied in Appendix A, where we show that the solution features and exponential density of states $W(E) = W_0e^{(1+\nu)E}$, where $\nu$ is a constant that depends on $H[s]$. This is an ideal limit, because in practical cases the maximisation is additionally constrained by the architecture of the learning machine and by the available data.

Notice that $H[s|\phi]$ provides a measure of that part of $H[s]$ which is not informative on the features. By Eq. (8), this is upper bounded by $H[s|E]$. Maximising $H[E]$ at fixed resolution $H[s]$ implies minimising $H[s|E] \ge H[s|\phi]$, i.e. squeezing noise out of the representation. An alternative definition of OLM, are machines that extract the most compressed representations from data, at a minimal information content $H[E]$ on the features

$$ \min_{I(s, \phi) \ge H[E]} H[s]. \qquad (13) $$

We remark that $H[E]$ is an unambiguous measure of information content when $E_s$ takes values in a discrete set. When $E_s$ can take any value on the real axis (as e.g. in RBMs), it is necessary to define the relevance with respect to a precision $\Delta$ [20]. In order to do this, we define the probability

$$ p_{\Delta}(E) = \sum_{s:|E_s - E| \le \Delta/2} e^{-E_s} \qquad (14) $$

that a state $s$ has energy in the interval of width $\Delta$ around $E$. The relevance at precision $\Delta$ is defined as

$$ H_{\Delta}[E] = - \sum_{E} p_{\Delta}(E) \log p_{\Delta}(E) \qquad (15) $$

where the sum on $E$ is restricted to integer multiples of $\Delta$. For a real random variable $E$ we expect that $H_{\Delta}[E] \approx h[E] - \log \Delta$ for $\Delta \to 0$, where $h[E]$ is the differential entropy [20]. We shall take $H_{\Delta}[E]$ as a measure of relevance, by making sure that $\Delta$ is chosen in the range where the relation $H_{\Delta}[E] \approx h[E] - \log \Delta$ is satisfied.

## 3 Exploring different architectures

In this section, we explore classes of statistical models in order to understand how the relevance $H[E]$ depends on the architecture. Under the hypothesis that good learning machines are those with high values of $H[E]$, this sheds light on the architectural features that confer superior learning capacity to models.
## 3.1 The Gaussian learning machine and maximum entropy models

Consider a Gaussian learning machine with $m$ visible units $\boldsymbol{x} = (x_1, \dots, x_m)$ and $n$ hidden units $\boldsymbol{s} = (s_1, \dots, s_n)$ which are all continuous variables. The joint probability distribution is given by

$$p(\boldsymbol{x}, \boldsymbol{s}) = \frac{1}{Z(\boldsymbol{b}, \hat{\boldsymbol{w}})} e^{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{b})^2 - \frac{1}{2}\boldsymbol{s}^2 + \boldsymbol{x}'\hat{\boldsymbol{w}}\boldsymbol{s}} \quad (16)$$

where prime denote transpose, $\boldsymbol{b} = (b_1, \dots, b_m)$ is a vector of parameters, $\hat{\boldsymbol{w}}$ is an $n \times m$ matrix of couplings, and $Z(\boldsymbol{b}, \hat{\boldsymbol{w}})$ is a normalisation constant. The couplings $\boldsymbol{b}$ and $\hat{\boldsymbol{w}}$ are adjusted so as to learn data generated from $\mathcal{H}(\boldsymbol{x})$. The marginal distribution of $\boldsymbol{x}$ is given by a Gaussian

$$p(\boldsymbol{x}) = \frac{(2\pi)^{n/2}}{Z(\boldsymbol{b}, \hat{\boldsymbol{w}})} e^{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{b})^2 + \frac{1}{2}\boldsymbol{x}'\hat{\boldsymbol{w}}\boldsymbol{x}}. \quad (17)$$

As shown in [16], this machine learns the first two moments of the distribution of $\boldsymbol{x}$. The dynamics of the parameters $\boldsymbol{b}$ and $\hat{\boldsymbol{w}}$ under different training algorithms has been discussed in Ref. [16] in detail. In what follows, $\boldsymbol{b}$ and $\hat{\boldsymbol{w}}$ are the parameters of the trained machine, so they depend on the data. The distribution in the internal layer is$^6$

$$p(\boldsymbol{s}) = \frac{(2\pi)^{m/2}}{Z(\boldsymbol{b}, \hat{\boldsymbol{w}})} e^{-\frac{1}{2}\boldsymbol{s}'(1-\hat{\boldsymbol{w}}'\hat{\boldsymbol{s}})\boldsymbol{s} + \boldsymbol{b}'\hat{\boldsymbol{w}}\boldsymbol{s}} \quad (18)$$

so the energy is

$$E(\boldsymbol{s}) = \frac{1}{2}\boldsymbol{s}'(1-\hat{\boldsymbol{w}}'\hat{\boldsymbol{s}})\boldsymbol{s} - \boldsymbol{b}'\hat{\boldsymbol{w}}\boldsymbol{s} - \frac{m}{2}\log(2\pi) + \log Z(\boldsymbol{b}, \hat{\boldsymbol{w}}) \quad (19)$$

$$= E_0(\hat{\boldsymbol{w}}) + \frac{1}{2}(\boldsymbol{s} - \boldsymbol{s}_0)'(1 - \hat{\boldsymbol{w}}'\hat{\boldsymbol{s}})(\boldsymbol{s} - \boldsymbol{s}_0) \quad (20)$$

where

$$E_0(\hat{\boldsymbol{w}}) = \frac{n}{2}\log(2\pi) - \log\sqrt{\det(1 - \hat{\boldsymbol{w}}'\hat{\boldsymbol{w}})} \quad (21)$$

is a constant (the ground state energy) and $\boldsymbol{s}_0(\boldsymbol{b}, \hat{\boldsymbol{w}}) = (1 - \hat{\boldsymbol{w}}'\hat{\boldsymbol{w}})^{-1}\hat{\boldsymbol{w}}\boldsymbol{b}$ is a vector. Notice that the resolution is given by:

$$H[\boldsymbol{s}] = E_0(\hat{\boldsymbol{w}}) + \frac{n}{2}. \quad (22)$$

The distribution of energies is

$$p(E) = \frac{1}{\Gamma(n/2)} (E - E_0)^{n/2-1} e^{-(E-E_0)}. \quad (23)$$

Remarkably, $p(E)$ is independent of $\boldsymbol{b}$ and it depends on $\hat{\boldsymbol{w}}$, i.e. on the data, only through the constant $E_0$. Hence the differential entropy [20] of this distribution

$$h(E) = \frac{n}{2} + \log \Gamma\left(\frac{n}{2}\right) - \left(\frac{n}{2} - 1\right) \log \psi\left(\frac{n}{2}\right) \approx \frac{1}{2} \log(4\pi en) + O(1/n) \quad (24)$$

is independent of the parameters $\boldsymbol{b}, \hat{\boldsymbol{w}}$, i.e. it is independent of the data. The scaling $h(E) \approx \frac{1}{2} \log n$ for large $n$ is consistent with the fact that asymptotically the distribution of energies is a Gaussian with variance proportional to $n$.

<footnote>$^6$The largest singular value of $\hat{\boldsymbol{w}}$ can be shown to be less than one, so that this is well defined.</footnote>
Let us now comment these results in terms of the discussion in Section 2. In the Gaussian learning machine, the extracted features $\phi(s) = \{s_i, i = 1, \dots, n, s_i s_j, i \le j = 1, \dots, n\}$ coincide with the sufficient statistics. A Gaussian learning machine is a maximum entropy model of the form

$$p(s|\theta) = \frac{1}{Z(\theta)} e^{\theta \cdot \phi(s)}, \qquad (25)$$

where the parameters $\theta = (b, \hat{w})$ are adjusted so as to reproduce the means and covariances of the data. As a consequence, the generative model $p(x)$ is always Gaussian, irrespective of what $\lambda$ is. A Gaussian learning machine does not learn anything on the shape of $\lambda$, because it assumes its shape from the outset. It only learns the parameters $\theta$. The amount of information learned on the features $I(s, \phi)$ equals $H[s]$, because $s$ is a part of $\phi$. Therefore, training increases the resolution $H[s]$ but it leaves the relevance $H[E]$ unchanged. The inequality (8) $I(s, \phi) \ge H[E]$ remains true, but is not informative.

An explicit manifestation of the poor performance of Gaussian learning machines, as compared e.g. to RBMs, in reproducing a dataset of handwritten digits, is given in Fig. 2. These findings suggest that the Gaussian multivariate model has a special place in statistical inference and the conclusions drawn from it should be taken with care.

<figure><img src="image_2.png" /><figcaption>Figure 2: Comparison between sample digits of the (reduced) MNIST dataset (top row) and those generated by a Gaussian learning machine (left) or a RBM (right) with $n$ hidden nodes, for $n = 10, \dots, 70$ (from bottom to top). Details on the dataset and on the RBM training experiments is provided in Section 4.2.</figcaption></figure>

With respect to the picture in Fig. 1, the Gaussian learning machine can only learn one point$^7$ $\bar{\phi}$, and it cannot reproduce a wide variation of features, as in Fig. 1.

We expect the same to be generally true for exponential models of the type of Eq. (25). Indeed, for large $n$, the distribution of the energy $E_s$ is very narrow, because with very high probability

$$E_s \approx \langle E \rangle_{\theta} = H[s].$$

<footnote>$^7$Conversely, the expected value $\bar{\phi}$ over the whole dataset of $\phi(s)$ is sufficient to estimate the parameters of this model.</footnote>
Schwab et al. [21] have shown that if the parameters $\theta$ themselves vary, as e.g. in

$$p(s) = \int d\theta p(s|\theta)p(\theta). \qquad (26)$$

then the energy acquires a broad variation and it exhibits an exponential density of state$^8 \log W(E) \propto E$. This is a defining feature of representations $p(s)$ that maximise the relevance. In the present picture, the variation of the hyper-parameters $\theta$ is induced by the distribution $\chi(x)$ of the data $x$. In other words, the hidden features $\eta$ play the role of the “hidden variables” of Ref. [21].

## 3.2 Ising Learning Machines

In this Section we shall consider a learning machines whose hidden variables $s$ are described by a fully connected model of Ising spins. These are maximum entropy models. So they are not particularly efficient as learning machines. Yet they provide a well defined playground for investigating the properties that make some architectures more efficient (i.e. with a larger $H[E]$) than others. This also allows us to contrast the properties of models with large $H[E]$ with those of well known models, such as the mean field Ising ferromagnet and the $J_{i,j} = \pm J$ spin glass. These insights make it possible to test whether models with the maximal $H[E]$ perform better than others in this class, on a specific application to real data. This will be the subject of Section 4.1.

We consider learning machines with a distribution of internal states given by

$$p(s) = \frac{1}{Z} e^{-\mathcal{H}(s)}, \quad Z = \sum_{s} e^{-\mathcal{H}(s)}, \qquad (27)$$

where $s = (s_1, \dots, s_n)$, $s_i = \pm 1$ are $n$ spins. The Hamiltonian takes the form

$$\mathcal{H}(s) = - \sum_{i<j} J_{i,j} s_i s_j, \qquad (28)$$

where the matrix $\hat{J}$ of couplings has elements $J_{i,j} \in \{\pm J\}$. The joint distribution over the visible variables $x \in \{0, 1\}^m$ and the internal state $s$ is given by $p(x, s) = p(x|s)p(s)$, where

$$p(x|s) = \frac{1}{Qs} \exp \left[ \sum_{j=1}^{m} \left( b_j + \sum_{i=1}^{n} s_i w_{i,j} \right) x_j \right], \quad Qs = \prod_{i=1}^{m} \left[ 1 + e^{b_j + \sum_i s_i w_{i,j}} \right]. \qquad (29)$$

The parameters $b$, $\hat{w}$ and $J$ are learned from the data, e.g. by maximum likelihood maximisation. Here we concentrate on the distribution $p(s)$ of the internal variables $s$.

We define the energy as the coding cost $E_s = -\log p(s) = \mathcal{H}(s) + \log Z$, and address the optimisation problem

$$\mathcal{J} = \arg \max_{\hat{J}: H[\hat{J}] = \bar{E}} H[E], \qquad (30)$$

where the maximum is taken over all matrices with elements $J_{i,j} = \pm J$, and on $J$. $\mathcal{J}$ is the subset of such matrices that achieve a maximal value of $H[E]$.

<footnote>$^8$Ref. [13] shows that this result is a general consequence of the Asymptotic Equipartition Property in e.g. [20]). If all typical states have the same probability $p(s|\theta) \sim e^{-nH[s]}$, then their number should be $W \sim e^{nH[s]}$.</footnote>
In order to compute $H[E]$ as a function of $\hat{J}$, we define the distribution of the energy as

$$p(E) = \sum_{\mathbf{s}} p(\mathbf{s}) \delta(E - E_{\mathbf{s}}) = W(E)e^{-E} \qquad (31)$$

where $W(E)$ is the number of states $\mathbf{s}$ with energy $E_{\mathbf{s}} = E$. From this and Eq. (2), we compute the value of $H[E]$.

The rest of this section is devoted to describe the trade-off between $H[E]$ and $H[\mathbf{s}]$ for models in the class of Eq. (28). Since $H[\mathbf{s}] \le n \log 2$ within this class, we shall discuss the rescaled resolution

$$h_s = \frac{H[\mathbf{s}]}{n \log 2}. \qquad (32)$$

Analogously, we observe that $(E_{\mathbf{s}} - \log Z)/J$ takes values only on the integers $[-\frac{n(n-1)}{2}, \frac{n(n-1)}{2}]$. This justifies the introduction of the rescaled relevance

$$h_E = \frac{H[E]}{\log n} \qquad (33)$$

that takes values in the interval $[0, 2)$. The sub-extensive nature of $H[E]$ suggests that properties that contribute to the relevance of a model may not be accessible to saddle point analysis of the partition function, that focuses only on the leading extensive terms. The behaviour of $H[E]$ on the number $n$ of spins depends on the fact that the model (28) contains only one parameter $J^9$.

The information content of a model can be quantified in terms of the Stochastic Complexity in Minimum Description Length [22]. This in our case is given by

$$SC = \int \frac{dJ}{J} \sqrt{V_J[E]}, \qquad (34)$$

where $V_J[E]$ is the variance of the distribution $p(E)$, which is closely related to the specific heat.

A second measure of the information content of a model can be obtained by the area under the $h_E$ vs $h_s$ curve. This scores models according to their relevance across different levels of resolution. We shall call this measure Multi-Scale Relevance (MSR), following [23], that introduced a similar measure on single samples.

Both $h_s$ and $h_E$ only depend on $\hat{J}$ through the degeneracy $W(E)$ of energy levels. $W(E)$ is invariant under the permutation of the spins. Hence if $\hat{J} \in \mathcal{J}$ then also the matrix with elements $J_{\pi(i),\pi(j)}$ is a solution, where $\pi(i)$ is a permutation of the indices $i = 1, \dots, n$. Also, matrices $\hat{J}$ and $\hat{J}'$ that are related by a gauge symmetry $J'_{i,j} = \tau_i J_{i,j} \tau_j$, with $\tau_i = \pm 1$, have the same $W(E)$. Hence, if $\hat{J} \in \mathcal{J}$, then also its gauge transformed $\hat{J}'$ belongs to $\mathcal{J}$. We partially exploit these symmetries by fixing the gauge with the choice $J_{1,j} = 1$ for all $j \ne 1$. Henceforth we shall focus on the reduced set $\mathcal{J}$ of matrices with this choice of the gauge.

An interesting point that we shall discuss is how the set $\mathcal{J}$ of solutions $\hat{J}$ are organised, i.e. whether they correspond to isolated maxima or whether they form a connected set of nearby solutions in terms of local moves. We consider local moves that correspond to “flipping” one

<footnote n="9">⁹For example, in a model with $H(s) = -\sum_{i<j} J_{i,j} s_i s_j - \sum_i h_i s_i$, $J_{i,j} = \pm J$ and $h_i = \pm h$, the energy takes values on a set of point $E = Jn_J + hn_h$ with $n_J \in [-n(n-1)/2, n(n-1)/2]$ and $n_h \in [-n, n]$. So the number of possible values of $E$ grows at most as $n^3$, and we expect the $h_E \in [0, 3)$. This suggests that $H[E]$ grows in general with $\log n$, with a coefficient that is bounded by a constant that grows linearly with the number of parameters.</footnote>
coupling $J_{i,j} \rightarrow -J_{i,j}$, or that “swap” a pair of opposite couplings $J_{i,j} = -J_{k,l}$, changing their sign ($J_{i,j}, J_{k,l} \rightarrow -J_{i,j}, -J_{k,l}$). Two models $\hat{J}$ and $\hat{J}'$ are connected if there is a sequence of such local moves that changes $\hat{J}$ into $\hat{J}'$. In this way, the set $\mathcal{J}$ can be divided into connected components. This analysis may shed light on the accessibility of solutions of Eq. (30) to local learning rules during training, and it connects to the finding of wide flat minima [6] in real learning machines<fnref n="10" />.

### 3.2.1 Exact enumeration for small systems

Each model in Eq. (28) is defined by a choice of the sign of the $n(n-1)/2$ couplings $J_{i,j}$ and by the strength $J$ of the couplings. There are $2^{n(n-1)/2}$ possible ways of choosing the signs. Yet, this number can be reduced by fixing the gauge. We were able to compute $h_E$ as a function of $h_s$ for all models up to $n=9$ spins. This allows us to find those models that achieve a maximal $h_E$ for a given resolution $h_s$.

<figure><img src="image_3.png" /></figure>

Figure 3: (Left) $h_E$ vs $h_s$ for all models with $n=8$ spins. The curves highlighted correspond to the models shown on the right plot (only $J_{i,j} < 0$ are shown as red links). The dotted line corresponds to the ferromagnet. The inset reports the relation between the stochastic complexity (SC) and the Multi-Scale Relevance (MSR), which is the area under the $h_E$ vs $h_s$ curve. The points corresponding to the models highlighted in the main plot are shown in green. (Right) Degeneracy $W(E)$ of energy levels for the model $a$ (top left) and for the ferromagnet. The dotted line corresponds to an exponential $W(E) = W_0 e^{\nu E}$ with $\nu = 0.2$.

Representative results are shown in Fig. 3 for $n=8$. There are $M=219$ different curves in the left panel, each corresponds to a different degeneracy $W(E)$ of energy levels<fnref n="11" />. Fig. 3(left) shows, in particular, the six functions $h_E$ that achieve maximal relevance for some value of $h_s$. A representative choice of the corresponding matrix $\hat{J} \in \mathcal{J}$ is shown in the right part of Fig. 3. We remark that there are many models with the same $W(E)$, that cannot be reduced to one another by a permutation of the spins (for model $a$ there are 325 different models). Fig. 3(right) shows the

<footnote n="10">¹⁰The latter refers to the fact that the landscape of the function that is optimised during training – e.g. the error or (minus) the likelihood – is characterised by wide flat minima in the space of parameters. If optimally trained machines correspond to models $\hat{J}$ of maximal relevance, then we expect that the set $\mathcal{J}$ has a large connected component.</footnote>

<footnote n="11">¹¹For $n=3, 4, 5, 6, 7, 8$ and $9$, we find $M=2, 3, 7, 16, 54, 219$ and $1625$, respectively.</footnote>
<table><caption>Table 1: Structure of the set of $\hat{J}$ corresponding to the OLM in Fig. 3 (right). The second column yields the size of the set $J$ of matrices $\hat{J}$ that achieve maximal $h_E$ in a given interval (4th column). The third column lists the connected components in $J$ under rewiring of a single negative coupling. The format used is size × multiplicity + .... All models in the same class share the same (ferromagnetic) ground state.</caption><thead><tr><th>OLM</th><th>|J|</th><th>connected components</th><th>h<sub>s</sub></th></tr></thead><tbody><tr><td>f</td><td>80640</td><td>80640 × 1</td><td>[0.999, 1]</td></tr><tr><td>e</td><td>40320</td><td>35280 × 1 + 5040 × 1</td><td>[0.995, 0.999]</td></tr><tr><td>d</td><td>10080</td><td>10080 × 1</td><td>[0.984, 0.995]</td></tr><tr><td>a</td><td>3360</td><td>2940 × 1 + 12 × 35</td><td>[0.71, 0.984]</td></tr><tr><td>c</td><td>10080</td><td>7560 × 1 + 1260 × 1 + 12 × 105</td><td>[0.62, 0.71]</td></tr><tr><td>b</td><td>3360</td><td>2100 × 1 + 420 × 2 + 4 × 105</td><td>[0.22, 0.62]</td></tr></tbody></table>

degeneracy $W(E)$ for model a) (top left in Fig. 3 right). As a comparison, Fig. 3 also shows the curve $h_E(h_s)$ (left) and the degeneracy $W(E)$ for the ferromagnet ($J_{i,j} = J \forall i, j$).

The architectures that achieve maximal relevance, that we identify with OLM, have the following distinctive properties: i) the large majority of the couplings $J_{i,j}$ are positive, which means that OLM are close to ferromagnetic models. Indeed, the ground state of all OLM are the same as that of the ferromagnet. Yet, ii) the relevance for OLM achieves a much higher value than that of the ferromagnetic model (dotted line in Fig. 3 left). This is consistent with the fact that $W(E)$ for the ferromagnet is concentrated on few values of $E$, whereas for OLM it spreads over a larger number of values of $E$. For model a) $W(E)$ is consistent with an exponential behaviour, as predicted in [13]. Finally iii) the architectures of OLM appear to be rather inhomogeneous, with negative $J_{i,j}$ impinging on a small subset of nodes. Such a low degree of symmetry, may be the origin of the fact that the number $|J|$ of matrices $\hat{J}$ that share the same $W(E)$ of OLM is rather large (see Table 1).

It is interesting to analyse the properties of the set $J$ with respect to local dynamics, as discussed earlier. We find (see Table 1 for $n = 8$) that $J$ is dominated by a single largest component, which is reminiscent of the wide flat minima discussed in [6]. As $h_s$ decreases, we see a tendency of the set to shrink and to fragment in different connected components. Yet, the union $J_a \cup J_b \cup \dots \cup J_f$ of all OLM $\in [a, b, \dots, f]$, is composed of a single connected component.

Finally, the inset of Fig. 3 left) shows that models with a higher relevance, as measured by the MSR, typically have a larger Stochastic Complexity (SC).

### 3.2.2 Some particular architectures for large $n$

This section analyses specific architectures for which the curve $h_E(h_s)$ can be computed for moderately large values of $n$. The matrices $\hat{J}$ of the four models that we consider are schematically depicted in Fig. 4 (left), together with the corresponding estimate of their capacity $h_E$. We discuss here only the main results and refer to the Appendix B for more details.

**The mean field Ising ferromagnet and spin glasses** The mean field Ising ferromagnet (MFIFM) corresponds to $J_{i,j} = J$ for all $i, j$ (model a in the left panel of Fig. 4). For large $n$, this is characterised by a high temperature disordered phase for $J < J_c/n$ and by a low temperature ordered phase for $J > J_c/n$ where a non-zero magnetisation spontaneously appears. The
energy for the MFIFM can take at most $n+1$ different values, which gives an upper bound $h_E \le 1$.

<figure><img src="image_4.png" /><figcaption>Figure 4: (left) The coupling matrices $\hat{J}$ of the four models discussed in Section 3.2.2: a) the MFIFM, b) the spin glass, c) the star model and d) the nested model. Red squares represent negative interactions. The maximal estimated value of $h_E$ is also reported for each model. (Right) $h_E$ vs $h_s$ for the MFIFM for different values of $n$. Inset: maximal value of $h_E$ as a function of $1/\log n$.</figcaption></figure>

Fig. 4 (right) reports the curves $h_E(h_s)$ for several values of $n$. The curve exhibits a maximum in the neighbourhood of the critical point $J \approx J_c/n$. The decreasing part to the right of the maximum corresponds to the high temperature phase, whereas the one to the left of the maximum to the low temperature phase. As $n$ increases, the maximum shifts to values of $h_s$ closer and closer to one. At the same time, the maximum gets sharper and sharper and its value slowly decreases towards a finite limit (see inset).

The limiting form of the curve $h_E(h_s)$ can be computed in the limit $n \to \infty$. This shows that

$$ \lim_{n \to \infty} h_E = \begin{cases} 1/2 & J \neq J_c/n \\ 3/4 & J = J_c/n, \end{cases} \qquad (35) $$

which is consistent with the bound $h_E \le 1$. The analysis leading to this result (see Appendix B) confirms that $h_E$ depends on sub-leading contributions that arise from the integration of fluctuations around the saddle point value. At the critical point $J = J_c/n$, the resolution converges to $h_s \to 1$. Yet, depending on how the limit $h_s \to 1$ is taken, all values of $h_E \in [1/2, 3/4]$ can be achieved. Indeed, for $J = 0$ we have $h_s = 1$ and $h_E = 1/2$, exactly. Fig. 4 (right) is fully consistent with this result, although the convergence is very slow.

When $J_{i,j} = \pm J$ are chosen at random with equal probability (model b in the left panel of Fig. 4), the model coincides with a spin glass. We know that the relevant scale for $J$ is $1/\sqrt{n}$. The ground state energy is extensive so $E_S$ spans a number of energy levels that is at most of order $n^{3/2}$. Therefore we expect that $h_E \le 3/2$ for a spin glass. In fact, for $J = c/\sqrt{n}$ with $c$ finite, the distribution of the energy extends over a range of order $\Delta E \propto \sqrt{n}$. Since the level spacing equals $2J \sim 1/\sqrt{n}$, this suggests that the number of energy levels on which $p(E)$ extends grows linearly with $n$. Hence we conjecture that $h_E \to 1$ as $n \to \infty$ for a model with randomly chosen $J_{i,j}$. This conclusion is corroborated by numerical results for systems of up to $n = 28$ spins.
**The star model** Let us now consider a model that differs from a MFIFM by that fact that one spin has antiferromagnetic interactions with $\ell < n$ other spins (model c in the left panel of Fig. 4).

<figure><img src="image_5.png" /><figcaption>Figure 5: (Left) Degeneracy of energy levels for $n = 256$ and $\ell = 2$ (green) and 128 (violet). Energies are rescaled by $E_0 = Jn(n-1)/2$. For $\ell = 2$ there are 385 distinct energy levels, whereas for $\ell = 128$ there are 8511. The inset magnifies the data in the range $E \in [-0.03E_0, 0]$. (Right) curves $h_E(h_s)$ for $\ell = n/2$ (full lines) and different values of $n$. Inset: maximal (over $h_s$) value of $h_E$ as a function of $\ell$, for $n = 256$ and 257.</figcaption></figure>

The degeneracy $W(E)$ of energy levels is shown in Fig. 5(left) for $n = 256$, $\ell = 2$ and $\ell = 128$. As $\ell$ increases, the number of different values that $E$ takes also increases. For the MFIFM ($\ell = 0$) $E$ can take only $n/2 + 1$ different values (for even $n$), whereas for $\ell = n/2$ we find that the number of possible values of $E$ grows almost as $n^2$. At the same time, $W(E)$ acquires a rapid variation as a function of $E$, so that for large $n$, $W(E)$ becomes a space filling curve (see Fig. 5 left and inset). Interestingly, the thermodynamics of the model is dominated by the convex envelope of $\log W(E)$, which is the same as that of the MFIFM. Therefore the star model and the MFIFM have indistinguishable thermodynamic properties. This is consistent with the fact that the number of $J_{i,j} = -1$ is a fraction of order $1/n$ of the total number of interactions and that they impinge on one out of the $n$ spins. It is easy to check that, for $\ell \le n/2$, the ground state of the star model is the ferromagnetic state $s_i = s_j$ for all $i \ne j$.

The maximal value of $h_E^*(\ell) = \max_{h_s} h_E$ for even values of $n$ is achieved at $\ell = n/2$ (see inset of Fig. 5 right), whereas for odd values of $n$, $h_E^*$ reaches a lower maximum at $\ell \approx n/3$. The symmetry for $\ell \to n-\ell$ of these curves is a consequence of the gauge transformation discussed above. The curves $h_E(h_s)$ for $\ell = n/2$ and for different values of $n$ are shown in Fig. 5(right). We find that the highest relevance $H[E]$ for even $n$ and $\ell = n/2$, grows faster than $\log n$. An extrapolation of the maximal value of $h_E$ is consistent with $h_E \to 5/4$ as $n \to \infty$. At fixed values of $h_s$ instead, we observe a slower growth of $H[E]$ with $\log n$, that is consistent with $h_E \to 1$ for large $n$. Note that the degeneracy of each energy level $E_m$ of the MFIFM is spread over $n$ other energy levels, $\sqrt{n}$ of which contribute to $H[E]$ (see inset of 5 left). This suggests a relation $h_E^{\text{Star}} = h_E^{\text{MFIFM}} + 1/2$, in agreement with numerical results.

Finally we remark that the set $\mathcal{J}$ of star models is connected in a single component, under a local dynamics updates of $\hat{J}$.
**The nested model** In the next architecture that we consider, the couplings in the lower right corner of the matrix $\hat{J}$ are negative, and those on the diagonal or above are positive (model *d*) of Fig. 4 left), i.e. $J_{i,j} = \text{sign}(n+\ell-i-j+1/2)$. In these models, the set of negative $J_{i,j}$ is organised according to a nested network [24]. Their properties are discussed in Appendix B where we also derive a recursion relation for $W(E)$ that allows us to compute $h_E$ for moderately large values of $n$. Fig. 6 shows the results for different values of $n$ and for $\ell = 1$ (left) and $n/2$ (right).

<figure><img src="image_6.png" /><figcaption>Figure 6: Relevance $h_E$ as a function of $h_s$ for the recursive model with $\ell = 1$ (left) and $\ell = n/2$ (right), for different values of $n$ (see legend). The inset shows the degeneracy $\log W(E)$ as a function of $E$ for the corresponding models.</figcaption></figure>

Although we could not access values of $n$ as large as in the previous cases, Fig. 6 shows that this model has qualitatively distinct features. First, we found that $\log W(E)$ approaches a continuous function of $E$ and it has support on a number of energy levels that increases as $n^2$. For $\ell = 1$ we found that $W(E) = W(-E)$ is a symmetric function of $E$ (see inset of Fig. 6 left). Second, in both cases, we found that $h_E$ increases monotonically with $n$ and it approaches a limiting value already for relatively small values of $n$. The maximal limiting values of $h_E$ is higher than that achieved with the other architectures ($h_E \approx 1.57$). We note that the number of anti-ferromagnetic interactions in these models is proportional to $n^2$, and hence the thermodynamic properties of nested models in the $n \to \infty$ limit differ from those of the MFIFM (even though the ground state is the same for $\ell \le n/2$).

There are $\frac{n!}{(n-\ell)!}$ possible models depending on how the set of $n-\ell$ fully ferromagnetic spins are chosen and on the ranking of the $\ell$ spins. This set of models forms a single connected component under local updates of the $J_{i,j}$. This is consistent with the conjecture that models with a high learning capacity are easily accessible by local updates of the couplings.

# 4 Numerical experiments

This section addresses the issue of whether the maximum relevance principle applies to learning experiments on real datasets. On the basis of this principle, we expect that *i*) models that achieve a higher value of $H[E]$ should describe better complex datasets, *ii*) within the precision allowed by a finite dataset, training converges to models of maximal relevance and *iii*) that learning induces broad distributions of energy levels. We address these three issues in turn.
## 4.1 Complex data is best described by models with maximal relevance

A prediction of the maximum relevance principle is that models that better fit a dataset are those with the highest relevance $H[E]$, provided the data contains a rich statistical structure. In order to test this prediction, we compared the ability of Ising models with different matrices $\hat{J}$ to represent a dataset, maximising the likelihood at different levels of resolution.

Because of the exponential growth<fnref n="12" /> in the number of models with the number $n$ of spins, we limit our analysis to the case $n = 5$. There are seven different classes of models (see Table 2). Each class contain models with the same $W(E)$ and hence with the same curve $h_E - h_s$. The class that achieves maximal value of $H[E]$ in a wide range of $H[s]$, that we shall call class A, contains two models<fnref n="13" />. Class A also has the largest values of the SC and of the MSR. The other classes – that we call B, C, D, E, FM and AF – contain one (E, AF and FM) or two models (B, C and D). Class FM contains the ferromagnetic model and AF the anti-ferromagnetic one.

<table><caption>Table 2: The seven classes of models for $n = 5$ spins. The last column lists the values $(i, j)$ for which $J_{i,j} < 0$. Some of the classes (FM, AF and E) correspond to a single model, the others to two. The second and third column lists the values of the MSR and SC respectively.</caption><thead><tr><th>Model</th><th>MSR</th><th>SC</th><th>{(i, j) : J<sub>i,j</sub> &lt; 0}</th></tr></thead><tbody><tr><td>FM</td><td>4.176</td><td>2.768</td><td>-</td></tr><tr><td>A</td><td>4.771</td><td>2.919</td><td>{(1, 2)}, {(1, 2), (1, 3), (1, 4)}</td></tr><tr><td>B</td><td>3.037</td><td>2.595</td><td>{(1, 2), (1, 3)}, {(1, 2), (1, 3), (2, 3), (3, 4)}</td></tr><tr><td>C</td><td>1.510</td><td>2.154</td><td>{(1, 2), (1, 3), (2, 3)}, {(1, 2), (1, 3), (1, 4), (2, 3), (2, 4)}</td></tr><tr><td>D</td><td>3.656</td><td>2.861</td><td>{(1, 4), (2, 3)}, {(1, 2), (1, 3), (1, 4), (2, 3))}</td></tr><tr><td>E</td><td>1.329</td><td>2.022</td><td>{(1, 2), (2, 3), (3, 4)}</td></tr><tr><td>AF</td><td>0.310</td><td>1.356</td><td>all</td></tr></tbody></table>

We focus on two datasets of $m = 9$ binary variables $x \in \{0, 1\}^m$. One is that used in [25] on the voting patterns in the US Supreme Court. This reports the votes (in favour or against) of the $m = 9$ judges of the US Supreme Court on $N = 895$ cases. The second is obtained from the MNIST dataset [17], focusing on a patch of $3 \times 3$ pixels of the $N = 6 \cdot 10^4$ hand written digits<fnref n="14" />, as in [7]. The joint distribution is built by taking $p(s)$ as in Eqs. (27,28) and $p(x|s)$ as in Eq. (29). For every model and each value of $J$ we maximise the log-likelihood over the parameters $b$ and $\hat{w}$ by gradient ascent. Fig. 7 shows that over a wide range of resolutions, the maximum of the likelihood is attained for models in class A. In particular, this interval contains the point where the log-likelihood achieves its maximum.

<footnote n="12">¹²To give an idea of the combinatorial complexity of the class of models, for $n = 8$ there are 219 different classes, each with many distinct models that cannot be reduced one to the other by a gauge transformation or by a permutation of the spins. For example, there are 325 models with the same $H[E] - H[s]$ curve as model a in Fig. 3.</footnote>

<footnote n="13">¹³One with only one negative coupling and the other with three negative couplings incident on the same node, e.g. $J_{1,2} = J_{1,3} = J_{1,4} < 0$.</footnote>

<footnote n="14">¹⁴Each MNIST data point contains $28 \times 28$ pixels, with an integer value in the range [0, 255]. We first transform the dataset into coarse grained pixels of $2 \times 2$ original pixels, and binarise the result applying a threshold. Then we focus on the $3 \times 3$ patch with pixels coordinates $(i, j)$ ranging from 5 to 7.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 7: Log likelihood for the USSC (left) and MNIST9 data for models of $n = 5$ spins, as a function of $H[\mathbf{s}]$. Models are classified according to their $h_E - h_s$ curve (see Table 2).</figcaption></figure>

## 4.2 RBMs achieve the maximal relevance compatible with a dataset

In order to test whether learning machines converge to models of maximal relevance during training, we focus on the architecture of Restricted Boltzmann Machines (RBM). These have $n$ hidden nodes and $n_v$ visible nodes, where $n_v$ matches the dimensionality of the dataset. We study a reduced MNIST dataset with $n_v = 144$ and $N = 6 \cdot 10^4$ data points. Details on RBMs and on the dataset are provided in Appendix C.

In general, the $2^n$ energy levels are learned from a sample of $N \ll 2^n$ data points. The finiteness of the data limits the precision with which the energy spectrum can be estimated. In order to take this into account, we observe that, with a change of variables $f = e^{-E}$, the differential entropies of $f$ and $E$ stand in the relation

$$h[E] = - \int_{0}^{1} df p(f) \log \left[ p(f) \left| \frac{df}{dE} \right| \right] = h[f] + H[\mathbf{s}] \quad (36)$$

where $H[\mathbf{s}] = \langle E \rangle$ is the average energy. The same relation holds for the variables at any precision $\Delta$, i.e.

$$H_{\Delta}[E] = H_{\Delta}[f] + H[\mathbf{s}] \quad (37)$$

The right hand side of Eq. (37) can be estimated in a finite sample $\hat{s}$, using the empirical distribution $\hat{f}_{\mathbf{s}} = k_{\mathbf{s}}/N$, where $k_{\mathbf{s}}$ is the number of times that state $\mathbf{s}$ occurs in the sample. Specifically, $H_{\Delta}[f]$ can be estimated using the empirical distribution of frequencies $P\{f_{\mathbf{s}} = k/N\} = km_k/N$, where $m_k$ is the number of states $\mathbf{s}$ that occur $k$ times in the sample. Hence

$$H_{\Delta}[f] \approx - \sum_{k} \frac{km_k}{N} \log \frac{km_k}{N} \equiv \hat{H}[k]. \quad (38)$$

Likewise, the second term in Eq. (37) can be computed as

$$H[\mathbf{s}] \approx - \sum_{\mathbf{s}} \frac{k_{\mathbf{s}}}{N} \log \frac{k_{\mathbf{s}}}{N} \equiv \hat{H}[\mathbf{s}]. \quad (39)$$
The two quantities $\hat{H}[s]$ and $\hat{H}[k]$ have been introduced earlier as a measure of resolution and relevance, respectively, within a sample [26, 27, 13]. Taken together, Eqs. (37, 38) and (39), imply

$$H_{\Delta}[E] \approx \hat{H}[k] + \hat{H}[s]. \qquad (40)$$

We remark that the precision $\Delta$ is implicitly defined by the sample. The fact that both $\hat{H}[k]$ and $\hat{H}[s]$ are biased estimates of the entropy is immaterial for our purposes. We remind that our goal is that of establishing whether a learning machine approaches states of maximal relevance, for a given dataset. This can be done comparing the values computed from a sample $\hat{s}$ drawn from $p(s)$, with the theoretical maximal value of $\hat{H}[k] + \hat{H}[s]$ that can be achieved in a sample of $N$ points. This has been computed in Ref. [27].

<figure><img src="image_8.png" /><figcaption>Figure 8: Relevance $H_{\Delta}[E] \approx \hat{H}[k] + \hat{H}[s]$ as a function of $\hat{H}[s] \approx H[s]$, for different training sessions of RBMs' with $n = 15, 20, 25, 30, 35, 40, 50, 60$ and $70$ hidden nodes (from left to right). Each session was run using persistent contrastive divergence [28] for 300 epochs with mini batches of 10 samples. Sample digits drawn from $p(s)$ at the end of training for different $n$ are shown in Fig. 2.</figcaption></figure>

Fig. 8 reports the results of numerical experiments on the reduced MNIST dataset, where the number of hidden nodes of the RBM ranges from $n = 15$ to $70$. Each training session lasts 300 epochs with persistent contrastive divergence [28]. Sample digits generated by the RBM are displayed in Fig. 2 (right). As a reference, the curve of maximal $\hat{H}[k] + \hat{H}[s]$ [27] is shown as a full black line. For each value of $n$, we observe trajectories that converge very fast to points that are close to the theoretical curve of maximal $\hat{H}[k] + \hat{H}[s]$, as first observed in [12]. As learning proceeds,
the internal state of the machine evolves towards more and more compressed representations (i.e. $\hat{H}[s]$ decreases).

### 4.2.1 RBM develop broad distribution of energies

Here we analyse the evolution of the energy distribution to test whether the distribution of energy broadens with learning, as predicted by relevance maximisation. The setting is the same of the previous section. We focus of RBM for small values of $n$, for which the distribution of the energy levels $E_s$, and hence $H_Δ[E]$, can be computed exactly during training. This allows us to monitor the evolution of the system as well as to appreciate the constraints imposed by the finiteness of the dataset and by the architecture. In order to do that, we take as a reference an ideal learning machine, that corresponds to the unconstrained maximisation of $H_Δ[E]$, at a fixed $H[s]$. The way in which this is computed parallels Ref. [13] and is detailed in Appendix A. As discussed in [13], the distinguishing feature of an ideal learning machine is that the degeneracy $W(E)$ of energy levels follows an exponential behaviour $W(E) \sim W_0e^{\nu E}$.

RBMs with $n \le 24$ hidden nodes have been trained with ($K = 10$) contrastive divergence [29] on the data. We used mini-batches of $10^3$ samples for the results shown in Fig. 9. Although these are not ideal for learning<fnref n="15" />, they allow us to monitor the dynamics of the RBM during learning, in the plane spanned by $H_Δ[E]$ and $H[s]$. As a function of $Δ$, the resolution is expected to behave as $H_Δ[E] = h[E] - \log Δ$, where $h[E]$ is the differential entropy [20] of the probability density function of energy levels. Indeed we find that $H_Δ[E] + \log Δ$ is nearly constant in the range $0.1 < Δ < 1$.

<figure><img src="image_9.png" /></figure>

Figure 9: a) typical trajectory during learning in the plane ($H[s], H_Δ[E]$) for a RBM with $n = 24$ trained on the reduced MNIST data ($Δ = 0.5$). The ideal limit of maximal $H_Δ[E] = h[E] - \log Δ$ (see Appendix A) is shown as a dashed line. b,c,d) Evolution of the distribution of energy levels for an RBM with $n = 22$ hidden units, after $n_e = 1, 10, 20, 50, 100, 200$ and 300 epochs. Panel b) shows the density $W(E)$ of energy levels, panel c) reports the probability density $p(E) = W(E)e^{-E}$. Panel d) shows the number $N(E)$ of data points that map to internal states of the machine with energy $E_s \in [E - Δ/2, E + Δ/2]$.

At different points in time, we compute the full spectrum of energy levels and compute $H_Δ[E]$.

<footnote n="15">¹⁵We confirmed that, as argued in [29], best performance is obtained with mini-batches of ~10 samples. Yet the trajectories obtained when the mini-batch size is small are very noisy.</footnote>
and $H[s]$. Fig. 9 (left) suggests that learning occurs in cycles, reminiscent of those discussed in Ref. [30]: A compression phase where $H[s]$ decreases, followed by an expansion phase where the system comes closer to the ideal limit of a machine with maximal $H_Δ[E]$ (dashed line).

Fig. 9 (right) shows the evolution of the spectrum of energy levels. As the bottom panel (b) shows, the spectrum of energy levels expands more and more during training. At times ($n_e = 1, n_e = 10$) when the system is closer to the ideal limit (dashed line in Fig. 9 left), the degeneracy of energy levels follows approximately an exponential behaviour $W(E) \approx W_0e^{\nu E}$, whereas it departs from it in the successive compression phases. In the late stages of the dynamics, for a given input $x$, the distribution $p(s|x)$ of internal states $s$ is sharply peaked. This makes it possible to map the dataset into a set of (clamped) internals states of the machine, as done in [12], and to compute the corresponding distribution $N(E)$ of energies. As shown in the top right panel d) of Fig. 9, the evolution of this distribution during learning follows the one of $p(E)$ (panel c) and converges to a similar shape.

# 5 Conclusion

Maximal relevance has been proposed as a principle that characterises optimal learning machines [13]. A distinctive feature of this approach is that it allows us to discuss the properties of learning machines without any reference to the data they have been trained with, as long as this contains a sufficiently rich structure.

The first contribution of this paper is to provide a relation between the relevance and feature extraction in learning. We show that the relevance $H[E]$ provides a lower bound to the mutual information between the state of the machine and the features that a learning machine extracts (see Section 2). This does not exclude the existence of efficient learning machines with low values of $H[E]$. Rather it ensures that learning machines with high values of $H[E]$ provide at least that much information on the features.

The second contribution of this paper is to investigate architectures that achieve maximal relevance, and to confirm that, when the constraints imposed by the finiteness of the dataset and by the architecture are taken into account, the internal representations of well trained learning machines approach the limit of maximal relevance.

We show that properties that confer a superior learning capacity, as measured by $H[E]$, may be related to sub-extensive quantities in the statistical mechanics treatment. This suggests that the properties responsible for the spectacular performance of learning machines may not be accessible to approaches that focus only on leading thermodynamic terms. A proper understanding of these systems requires a detailed analysis of saddle point fluctuations. This result also tallies with the observation that complexity reveals itself in the sub-extensive contribution to the entropy, both in statistical models<fnref n="16" /> [31, 22] and in time series<fnref n="17" /> [32].

For Ising systems, we find that local update rules can explore a large part of the set of OLM, suggesting that the relevance $H[E]$ is characterised by wide flat maxima in the space of the parameters. It is tempting to relate this property to the wide flat minima discussed by Baldassi et

<footnote n="16">¹⁶Minimum Description Length, as well as Bayesian inference, identify the leading order term that should be used in model selection to penalise the likelihood of models for their complexity. In both cases this term is $\frac{d}{2} \log N$, where $d$ is the number of parameters of the model and $N$ is the number of data points. This term is sub-leading with respect to the log-likelihood, which is proportional to $N$.</footnote>

<footnote n="17">¹⁷The complexity of a time series is defined in [32] as the predictive information, which is the mutual information between the recent past, in a window of $T$ time points, and the future of the time series.</footnote>
al. [6].

Our results also contribute to the discussion on the criticality of learned models [33, 7, 34]. We confirm that models have a superior learning performance when poised at the critical point. Yet, it is not strictly necessary for an OLM to be poised at a critical point marking the transition between two different phases. Indeed, we find recursive models with a learning performance superior to that of models that feature a phase transition, in a wide range of the resolution scales.

We hope that these results will not only contribute to our understanding of learning machines, but that they may also pave the way to applications that may improve further the performance of learning machines.

## Acknowledgments

We benefitted from discussions with Jean Barbier, Riccardo Zecchina, Federico Ricci-Tersenghi and Yasser Roudi. We are grateful to the authors of Re. [25] for sharing their data with us. We thank an anonymous referee of MSML2020 for pointing out the example in Section 3.1.

# A The ideal limit of an Optimal Learning Machine

For a learning machine with many hidden units ($n \gg 1$), we can approximate the distribution of energies with a continuous probability density $p(E)$, that we assumed to be defined in an interval $E \in [0, \bar{E}]$. For $\Delta$ small enough, the behaviour of $H_\Delta[E]$ with $\Delta$ is given by

$$H_\Delta[E] \approx h[E] - \log \Delta \qquad (41)$$

where $h[E]$ is the differential entropy [20]

$$h[E] = - \int_{0}^{\bar{E}} dE p(E) \log p(E). \qquad (42)$$

For a fixed $\Delta$, the ideal limit of a learning machine with maximal $H_\Delta[E]$ is then obtained by finding the $p(E)$ that makes $h[E]$ maximal.

Let $W(E)$ be the density of energy levels, i.e. $W(E)dE$ is the number of states $s$ with $E_s \in [E, E + dE)$. The probability density function $p(E)$ is then given by $p(E) = W(E)e^{-E}$. When no other constraint is imposed, Optimal Learning Machines are defined by the solution of the optimisation problem

$$W^*(E) = \arg \max_{W(E)} h[E] \qquad (43)$$

$$\text{s.t.} \quad H[\mathbf{s}] = \langle E \rangle \qquad (44)$$

$$\int_{0}^{\bar{E}} W(E) dE = N \qquad (45)$$

$$\int_{0}^{\bar{E}} W(E)e^{-E} dE = 1, \qquad (46)$$

where

$$h[E] = - \int_{0}^{\bar{E}} dE p(E) \log p(E) = \langle E \rangle - \langle \log W(E) \rangle, \qquad (47)$$
and N is the total number of states s. The solution of this problem is found introducing Lagrange multipliers to enforce the constraints (44,46), and fixing $\bar{E}$ to satisfy Eq. (45). The solution reads

$$W(E) = \frac{\nu}{e^{\nu\bar{E}} - 1} e^{(1+\nu)\bar{E}}, \qquad (48)$$

$$H[s] = \langle E \rangle = \frac{\bar{E}}{1 - e^{-\nu\bar{E}}} - \frac{1}{\nu}, \qquad (49)$$

$$h[E] = \log \frac{e^{\nu\bar{E}} - 1}{\nu} - \frac{\nu\bar{E}}{1 - e^{-\nu\bar{E}}} + 1, \qquad (50)$$

$$N = \frac{\nu}{1 + \nu} \frac{e^{(1+\nu)\bar{E}} - 1}{e^{\nu\bar{E}} - 1}, \qquad (51)$$

where Eq. (51) needs to be solved to obtain $\bar{E}$ as a function of N and Eq. (49) yields $\nu$ as a function of $H[s]$. The limit $\nu \to 0$ is revealing to understand the scaling of $H[s]$ and $h[E]$ on N. We find

$$\bar{E} = \log(1 + N\bar{E}) \approx \log N + \log \log N + O\left(\frac{\log \log N}{\log N}\right) \qquad (\nu = 1) \qquad (52)$$

$$H[s] = \frac{1}{2}\bar{E} \approx \frac{1}{2}\log N \qquad (53)$$

$$h[E] = \log \bar{E} \approx \log \log N. \qquad (54)$$

At this point, the curve $h[E]$ vs $H[s]$ is close to its maximum. So Eq. (54) provides an upper bound on the possible value of $h[E]$. In the case of a Boltzmann machine with $n$ hidden binary units, there are $N = 2^n$ hidden states. Hence $H[s]$ grows linearly with $n$ and $H[E] \propto \log n$ grows with the logarithm of the number of hidden units $n$. This is the same scaling that we have found for discrete energy levels, with $\Delta = J$.

## B Spin models

**The MFIFM model** The energy for the MFIFM ranges over an interval $[-Jn(n-1)/2, Jn(n-1)/2]$ of order $n^2$, but energies $E_s$ can take at most $n+1$ different values

$$E_m = -\frac{J}{2}(m^2 - n), \quad m = \sum_{i=1}^{n} s_i = -n, -n+2, \dots, n-2, n. \qquad (55)$$

This implies that $h_E \le \log(n+1)/\log n$. Note also that there is a single matrix $\hat{J}$ which corresponds to the degeneracy of energy levels $W(E)$ of the MFIFM (apart from gauge transformations)..

The distribution of energy levels is given by $p(E) = W(E)e^{-E}/Z$, where $Z = \sum_E W(E)e^{-E}$. Hence the relevance is given by

$$H[E] = \log Z - \langle \log W(E) \rangle + \langle E \rangle. \qquad (56)$$

The degeneracy of energy level $E_m$ in Eq. (55) is

$$W(E_m) = 2 \left( \frac{n}{\frac{n+m}{2}} \right) \approx \sqrt{\frac{2}{\pi n(1 - \mu^2)}} e^{nS(\mu)}, \quad \mu = \frac{m}{n} \in [-1, 1]$$
where $S(\mu) = -\frac{1+\mu}{2} \log \frac{1+\mu}{2} - \frac{1-\mu}{2} \log \frac{1-\mu}{2}$. In the regime where $h_s \in [0, 1]$ is finite, $J$ is of order $1/n$ and, for $n \to \infty$, the partition function is dominated by the saddle point $\mu^* = \arg \min_{\mu} f(\mu)$ with $f(\mu) = Jn\mu^2/2 - S(\mu)$. For $Jn \neq 1$, the first term in Eq. (56) can be computed with integration over the gaussian fluctuations

$$Z \approx \sqrt{\frac{2n}{\pi}} \int_{-1}^{1} \frac{d\mu}{\sqrt{1-\mu^2}} e^{-nf(\mu)} \approx \sqrt{\frac{2}{\pi(1-\mu^{*2})f''(\mu^*)}} e^{-nf(\mu^*)}$$

Here, the $\sqrt{n}$ factor is canceled because $f(\mu) - f(\mu^*) \sim (\mu - \mu^*)^2$, and the change of variables $z = \sqrt{n}(\mu - \mu^*)$ generates a $1/\sqrt{n}$ term.

Hence $\log Z \approx -nf(\mu^*) + \text{const.}$ The extensive term in $\log Z$ is canceled by an analogous term in

$$\langle \log W(E) \rangle - \langle E \rangle \approx nf(\mu^*) + \log \sqrt{2\pi n(1-\mu^2)}$$

so, to leading order, $H[E] \approx \frac{1}{2} \log n + \text{const.}$, which implies $h_E \to 1/2$ as $n \to \infty$.

For $Jn = 1$, instead, $f(\mu) - f(\mu^*) \approx a(\mu - \mu^*)^4 + \dots$. This implies that in the calculation of $Z$ we need a change of variables $z = n^{1/4}(\mu - \mu^*)$ that yields $\log Z \approx \frac{1}{4} \log n + \text{const.}$ This additional term is responsible for the asymptotic behaviour $h_E \to 3/4$ for $Jn = 1$.

**The spin glass model** We computed numerically the value of $\max_J h_E$ for 100 realisations of models with random $J_{i,j} = \pm 1$ (with $P\{J_{i,j} = +1\} = 1/2$), for systems of up to $n = 28$ spins. A linear fit of $H[E]$ versus $\log n$ yields $h_E \approx 1.02 \pm 0.01$.

**The star model** In the start model we take $J_{i,j} = -1$ if $i = 1$ and $j \le \ell + 1$, and $J_{i,j} = +1$ otherwise. We divide the $n-1$ group of spins, except spin $s_1$, into two groups: the $\ell$ ones with anti-ferromagnetic interaction with $s_1$ and the remaining $n-\ell-1$. If $q$ and $k$ are the number of positive spins in the first and second group, respectively, the energy is given by

$$\epsilon(k, q, s_1) = J \left[ s_1(n-1-2\ell-2k+2q) - \frac{1}{2}(2k+2q-n+1)^2 + \frac{n-1}{2} \right]. \quad (57)$$

The degeneracy of energy levels is given by

$$W(E) = \sum_{k=1}^{n-1-\ell} \sum_{q=1}^{\ell} \sum_{s_1=\pm 1} \binom{\ell}{q} \binom{n-1-\ell}{k} \delta_{E,\epsilon(k,q,s_1)}. \quad (58)$$

Notice that for $\ell = 0$ one recovers the MFIFM. A gauge transformation $J_{i,j} \to \tau_i J_{i,j} \tau_j$, with $\tau_1 = -1$ and $\tau_i = +1$ for $i > 1$, maps a model with $\ell > n/2$ into a model with $\ell' = n - \ell \le n/2$. So it is sufficient to study the model for $\ell \le n/2$.

Each models $\mathcal{J}$ with the star structure can be “reached” from any other model, under a dynamics where a single $J_{i,j} \mapsto -J_{i,j}$ changes sign. If we consider a “rewiring” dynamics where the signs of a pair $(J_{i,j}, J_{i,k}) = (-1, +1)$ of links are swapped $(J_{i,j}, J_{i,k}) \mapsto (-J_{i,j}, -J_{i,k})$, each of the $\ell$ negative links can be rewired to any of the $n-\ell-1$ spins with all ferromagnetic interactions. This dynamics, spans the whole subset of models in $\mathcal{J}$ with the same hub node. These subsets are related by a permutation of the nodes.
**The nested model** The triangle model has

$$J_{i,j} = \begin{cases} -1 & \text{for } i > \ell, j > n + \ell - i \\ +1 & \text{otherwise} \end{cases}$$

where $\ell = 1, \dots, n-1$. For $\ell = 1$ this corresponds to a matrix $\hat{J}$ where all the elements below the diagonal $i+j=n+1$ are negative, as shown in Fig. 10(left). In these models, the set of negative $J_{i,j}$ is organised according to a nested network [24].

<figure><img src="image_10.png" /><figcaption>Figure 10: Left: Sketch of the structure of the matrix $\hat{J}$ for structure of recursive models for $\ell = 1$. Red indicates negative $J_{i,j}$. Right: sketch of the recursion that relates a model with $n$ spins to a model with $n+2$ spins.</figcaption></figure>

These matrices can be obtained through a recursive procedure (see Fig. 10 right): we start from $n=2$ with $J_{1,2} = +1$. At each iteration, from a system of $n$ spins, we add two spins, $s_0$ and $s_{n+1}$. The spin $s_0$ has $J_{0,j} = -1$ with all other spins $j \le n$. The spin $s_{n+1}$ has $J_{i,n+1} = +1$ with all spins, including $s_0$. This generates a matrix with $n+2$ spins with the desired structure$^{18}$.

It is possible to generalise the recursion relation to obtain matrices with different values of $\ell$. For example, the case $\ell = n/2$ can be obtained by separating the set of spins into two equal parts, $\mathcal{I}_{\le} = \{i \le n/2\}$ and $\mathcal{I}_{>} = \{i > n/2\}$. To each of the two sets we add two spins at each iteration, thereby obtaining a system with $n+4$ spins. Of the two spins added to $\mathcal{I}_{\le}$, one has $J_{0,j} = -1$ with all spins $j \in \mathcal{I}_{\le}$ and $J_{0,j} = +1$ for all $j \in \mathcal{I}_{>}$. The other three spins have ferromagnetic interactions with all other spins, as well as among themselves. It is easy to see by induction, that the ground state of this model is the ferromagnetic one $s_i = s_j$ for all $i, j$. This holds true for all $\ell \le n/2$.

Numerical iteration of the recursion relations allowed us to compute the degeneracy $W_n(E)$ and the curves $h_E(h_s)$ for different values of $n$ and $\ell$.

$$E' = E - J[s_0s_{n+1} - s_0M + s_{n+1}M] \quad (59)$$

and magnetisation $M' = M + s_0 + s_{n+1}$. Thus the energy degeneracy of a system of $n$ spins can be computed as $W_n(E) = \sum_M w_n(E, M)$, where $w_n$ satisfies the recursive equation

$$w_{n+2}(E', M') = \sum_{s_0, s_{n+1} = \pm 1} \sum_{E, M} w_n(E, M) \delta_{E'(E)} \delta_{M', M+s_0+s_{n+1}} \quad (60)$$

Here the shorthand $\delta_{E'(E)}$ reduces the sum to only those terms where Eq. (59) is satisfied.

<footnote>$^{18}$For each state at energy $E$ and magnetisation $M = \sum_{i=1}^n s_i$, and each choice of $s_0, s_{n+1}$, we obtain a state of the system with $n+2$ spins with energy</footnote>
For a given value of $\ell$, there are $\frac{n!}{(n-\ell)!}$ nested models depending on how the set of $n-\ell$ fully ferromagnetic spins are chosen and on the ranking of the $\ell$ remaining spins in the nested hierarchy. This set of models forms a single connected component under a dynamics where a negative link $J_{i,j} = -1$ is “rewired” to a different node $k$, i.e. $(J_{i,j}, J_{k,j}) = (-1, +1) \mapsto (+1, -1)$. In order to show this, let $\partial_i$ be the set of spins with $J_{i,j} = -1$. If we sort spins in such a way that $\delta_1 = (1, \dots, \ell)$, $\delta_2 = (1, \dots, \ell-1)$, $\dots$, $\delta_k = (1, \dots, \ell-k+1)$, $\dots$, then the rewiring of the link $(k, \ell-k+1)$ to $(k+1, \ell-k+1)$ interchanges the order of spins $k$ and $k+1$ in the hierarchy. The rewiring of the link $(1, \ell)$ to $(1, k)$ with $k > \ell$ changes the membership of spins in the network connected by $J_{i,j} = -1$. Repeated application of these two moves in the appropriate order makes it possible to reach any model in the set of nested models, from any other model.

## C RBMs and the reduced MNIST dataset

We focus on the architecture of Restricted Boltzmann Machines (RBM). These have $n_v$ visible units $\mathbf{x} = (x_1, \dots, x_{n_v})$, and $n$ hidden binary units $\mathbf{s} = (s_1, \dots, s_n)$ with $s_i = 0, 1$. The joint probability distribution is

$$p(\mathbf{x}, \mathbf{s}) = \frac{1}{Z} e^{\sum_i b_i s_i + \sum_j c_j x_j + \sum_{i,j} s_i w_{i,j} x_j} \quad (61)$$

where $Z = Z(\mathbf{b}, \mathbf{c}, \hat{\mathbf{w}})$ is the partition function. The parameters $\mathbf{b}, \mathbf{c}$ and $\hat{\mathbf{w}}$ are learned from the data, in such a way as to maximise the log-likelihood $p(\mathbf{x})$. The energy of hidden state $\mathbf{s}$ is defined as

$$E_{\mathbf{s}} = -\log p(\mathbf{s}) = \log Z - \sum_{i} b_{i} s_{i} - \sum_{j} \log \left[ 1 + e^{c_{j} + \sum_{i} s_{i} w_{i,j}} \right]. \quad (62)$$

Since $E_{\mathbf{s}}$ can take any real value, we measure the relevance with respect to a precision $\Delta$ of the energy spectrum, as in Eq. (15).

In what follows, we study a reduced MNIST dataset [17], where we first coarse grain the original 28 × 28 pixels data in cells of 2 × 2 pixels, that are transformed into binary values by applying a threshold<fnref n="19" />. Finally we discard the first and the last rows and column of cells, thereby leaving us with a dataset of $N = 60000$ images of 12 × 12 binary pixels, corresponding to $n_v = 144$. The reason for working with such a reduced dataset is that it allows us to compute the exact distribution of energies, for values of $n$ in a relevant range. Indeed, the evaluation of the energy requires the calculation of the partition function $Z$, which can be computed only for small systems. Sect. 4.2.1 presents results for RBMs with a number of hidden units such that the exact calculation of $p_{\Delta}(E)$ on all $2^n$ of states is within reach of our computational resources. For larger systems we sample $N \ll 2^n$ states $\mathbf{s}$ from $p(\mathbf{s})$ with MCMC method and estimate $H_{\Delta}(E)$ from the sample (see Sect. 4.2).

In order to test whether learning machines converge to models of maximal relevance during training, we focus on the architecture of Restricted Boltzmann Machines (RBM). These have $n_v$ visible units $\mathbf{x} = (x_1, \dots, x_{n_v})$, and $n$ hidden binary units $\mathbf{s} = (s_1, \dots, s_n)$ with $s_i = 0, 1$. The joint probability distribution is

$$p(\mathbf{x}, \mathbf{s}) = \frac{1}{Z} e^{\sum_i b_i s_i + \sum_j c_j x_j + \sum_{i,j} s_i w_{i,j} x_j} \quad (63)$$

<footnote n="19">¹⁹If the sum of the grey levels of the four pixels exceeds 400, the coarse grained pixel is assigned the value one, otherwise it is zero.</footnote>
where Z is the partition function, which depends on the parameters b, c and $\hat{w}$ that are learned from the data, in such a way as to maximise the log-likelihood $p(\boldsymbol{x})$. The energy of hidden state s is defined as

$$E_s = -\log p(s) = \log Z - \sum_i b_i s_i - \sum_j \log [1 + e^{c_j + \sum_i s_i w_{i,j}}]. \quad (64)$$

Since $E_s$ can take any real value, we measure the relevance with respect to a precision $\Delta$ of the energy spectrum, as in Eq. (15).

## References

[1] J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proc Natl Acad Sci U S A, vol. 79, no. 8, pp. 2554–8, 1982.

[2] D. H. Ackley, G. E. Hinton, and T. J. Sejnowski, “A learning algorithm for boltzmann machines*,” Cognitive Science, vol. 9, no. 1, pp. 147–169, 1985.

[3] M. Baity-Jesi, L. Sagun, M. Geiger, S. Spigler, G. B. Arous, C. Cammarota, Y. LeCun, M. Wyart, and G. Biroli, “Comparing dynamics: Deep neural networks versus glassy systems,” in Proceedings of the 35th International Conference on Machine Learning (J. Dy and A. Krause, eds.), vol. 80 of Proceedings of Machine Learning Research, (Stockholmsmassan, Stockholm Sweden), pp. 314–323, PMLR, 10–15 Jul 2018.

[4] J. Tubiana and R. Monasson, “Emergence of compositional representations in restricted boltzmann machines,” Physical review letters, vol. 118, no. 13, p. 138301, 2017.

[5] A. Decelle, G. Fissore, and C. Furtlehner, “Thermodynamics of restricted boltzmann machines and related learning dynamics,” Journal of Statistical Physics, vol. 172, pp. 1576–1608, Sep 2018.

[6] C. Baldassi, C. Borgs, J. T. Chayes, A. Ingrosso, C. Lucibello, L. Saglietti, and R. Zecchina, “Unreasonable effectiveness of learning neural networks: From accessible states and robust ensembles to basic algorithmic schemes,” Proceedings of the National Academy of Sciences, vol. 113, no. 48, pp. E7655–E7662, 2016.

[7] M. E. Rule, M. Sorbaro, and M. H. Hennig, “Optimal encoding in stochastic latent-variable Models,” ArXiv e-prints, p. arXiv:1802.10361, Feb. 2018.

[8] M. Mézard, “Mean-field message-passing equations in the hopfield model and its generalizations,” Phys. Rev. E, vol. 95, p. 022117, Feb 2017.

[9] S. Goldt, M. Mézard, F. Krzakala, and L. Zdeborová, “Modelling the influence of data structure on learning in neural networks,” arXiv preprint arXiv:1909.11500, 2019.

[10] P. Rotondo, M. C. Lagomarsino, and M. Gherardi, “Counting the learnable functions of geometrically structured data,” Phys. Rev. Research, vol. 2, p. 023169, May 2020.

[11] L. Zdeborová, “Understanding deep learning is also a job for physicists,” Nature Physics, pp. 1–3, 2020.
[12] J. Song, M. Marsili, and J. Jo, "Resolution and relevance trade-offs in deep learning," *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2018, p. 123406, dec 2018.

[13] R. J. Cubero, J. Jo, M. Marsili, Y. Roudi, and J. Song, "Statistical criticality arises in most informative representations," *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2019, p. 063402, jun 2019.

[14] T. Mora and W. Bialek, "Are biological systems poised at criticality?," *Journal of Statistical Physics*, vol. 144, no. 2, pp. 268-302, 2011.

[15] R. Cubero, M. Marsili, and Y. Roudi, "Minimum description length codes are critical," *Entropy*, vol. 20, p. 755, Oct 2018.

[16] R. Karakida, M. Okada, and S.-i. Amari, "Dynamical analysis of contrastive divergence learning: Restricted boltzmann machines with gaussian visible units," *Neural Networks*, vol. 79, pp. 78-87, 2016.

[17] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, "Gradient-based learning applied to document recognition," *Proceedings of the IEEE*, vol. 86, no. 11, pp. 2278-2324, 1998.

[18] A. Ansuini, A. Laio, J. H. Macke, and D. Zoccolan, "Intrinsic dimension of data representations in deep neural networks," in *Advances in Neural Information Processing Systems*, pp. 6111-6122, 2019.

[19] A. Ilyas, S. Santurkar, L. Engstrom, B. Tran, and A. Madry, "Adversarial examples are not bugs, they are features," *Advances in neural information processing systems*, vol. 32, pp. 125-136, 2019.

[20] T. M. Cover and J. A. Thomas, *Elements of information theory*. John Wiley & Sons, 2012.

[21] D. J. Schwab, I. Nemenman, and P. Mehta, "Zipf's law and criticality in multivariate data without fine-tuning," *Phys. Rev. Lett.*, vol. 113, p. 068102, Aug 2014.

[22] I. J. Myung, V. Balasubramanian, and M. A. Pitt, "Counting probability distributions: Differential geometry and model selection," *Proceedings of the National Academy of Sciences*, vol. 97, no. 21, pp. 11170-11175, 2000.

[23] R. J. Cubero, M. Marsili, and Y. Roudi, "Finding informative neurons in the brain using Multi-Scale Relevance," *ArXiv e-prints*, Feb. 2018.

[24] S. Jonhson, V. Dominguez-Garcia, and M. A. Muñoz, "Factors determining nestedness in complex networks," *PLoS One*, vol. 8, no. 9, p. e74025, 2013.

[25] E. D. Lee, C. P. Broedersz, and W. Bialek, "Statistical Mechanics of the US Supreme Court," *Journal of Statistical Physics*, vol. 160, pp. 275-301, July 2015.

[26] M. Marsili, I. Mastromatteo, and Y. Roudi, "On sampling and modeling complex systems," *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2013, no. 09, p. P09003, 2013.

[27] A. Haimovici and M. Marsili, "Criticality of mostly informative samples: a bayesian model selection approach," *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2015, no. 10, p. P10013, 2015.
[28] T. Tieleman, “Training restricted boltzmann machines using approximations to the likelihood gradient,” in *Proceedings of the 25th international conference on Machine learning*, pp. 1064–1071, 2008.

[29] G. E. Hinton, “A practical guide to training restricted boltzmann machines,” in *Neural networks: Tricks of the trade*, pp. 599–619, Springer, 2012.

[30] O. Shamir, S. Sabato, and N. Tishby, “Learning and generalization with the information bottleneck,” *Theoretical Computer Science*, vol. 411, no. 29-30, pp. 2696–2711, 2010.

[31] G. Schwarz, “Estimating the dimension of a model,” *Ann. Statist.*, vol. 6, pp. 461–464, 03 1978.

[32] W. Bialek, I. Nemenman, and N. Tishby, “Predictability, complexity, and learning,” *Neural Computation*, vol. 13, no. 11, pp. 2409–2463, 2001.

[33] I. Mastromatteo and M. Marsili, “On the criticality of inferred models,” *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2011, no. 10, p. P10012, 2011.

[34] M. A. Muñoz, “Colloquium: Criticality and dynamical scaling in living systems,” *Rev. Mod. Phys.*, vol. 90, p. 031001, Jul 2018.