Parallel Learning by Multitasking Neural Networks
A modern challenge of Artificial Intelligence is learning multiple patterns
at once (i.e.parallel learning). While this can not be accomplished by standard
Hebbian associative neural networks, in this paper we show how the Multitasking
Hebbian Network (a variation on theme of the Hopfield model working on sparse
data-sets) is naturally able to perform this complex task. We focus on systems
processing in parallel a finite (up to logarithmic growth in the size of the
network) amount of patterns, mirroring the low-storage level of standard
associative neural networks at work with pattern recognition. For mild dilution
in the patterns, the network handles them hierarchically, distributing the
amplitudes of their signals as power-laws w.r.t. their information content
(hierarchical regime), while, for strong dilution, all the signals pertaining
to all the patterns are raised with the same strength (parallel regime).
Further, confined to the low-storage setting (i.e., far from the spin glass
limit), the presence of a teacher neither alters the multitasking performances
nor changes the thresholds for learning: the latter are the same whatever the
training protocol is supervised or unsupervised. Results obtained through
statistical mechanics, signal-to-noise technique and Monte Carlo simulations
are overall in perfect agreement and carry interesting insights on multiple
learning at once: for instance, whenever the cost-function of the model is
minimized in parallel on several patterns (in its description via Statistical
Mechanics), the same happens to the standard sum-squared error Loss function
(typically used in Machine Learning).
Elena Agliari, Andrea Alessandrelli, Adriano Barra, Federico Ricci-Tersenghi
# Parallel Learning by Multitasking Neural Networks

###### Abstract

A modern challenge of Artificial Intelligence is learning multiple patterns at once (i.e. _parallel learning_). While this can not be accomplished by standard Hebbian associative neural networks, in this paper we show how the Multitasking Hebbian Network (a variation on theme of the Hopfield model working on sparse data-sets) is naturally able to perform this complex task. We focus on systems processing in parallel a finite (up to logarithmic growth in the size of the network) amount of patterns, mirroring the low-storage level of standard associative neural networks at work with pattern recognition. For mild dilution in the patterns, the network handles them hierarchically, distributing the amplitudes of their signals as power-laws w.r.t. their information content (hierarchical regime), while, for strong dilution, all the signals pertaining to all the patterns are raised with the same strength (parallel regime).

Further, confined to the low-storage setting (i.e., far from the spin glass limit), the presence of a teacher neither alters the multitasking performances nor changes the thresholds for learning: the latter are the same whatever the training protocol is supervised or unsupervised. Results obtained through statistical mechanics, signal-to-noise technique and Monte Carlo simulations are overall in perfect agreement and carry interesting insights on _multiple learning at once_: for instance, whenever the cost-function of the model is minimized _in parallel on several patterns_ (in its description via Statistical Mechanics), the same happens to the standard sum-squared error Loss function (typically used in Machine Learning).
###### Contents

* 1 Introduction
* 2 Parallel learning in multitasking Hebbian neural networks
	* 2.1 A preliminary glance at the emergent parallel retrieval capabilities
	* 2.2 From parallel storing to parallel learning
* 3 Parallel Learning: the picture by statistical mechanics
	* 3.1 Study of the Cost function and its related Statistical Pressure
		* 3.1.1 Low-entropy data-sets: the Big-Data limit
		* 3.1.2 Ergodicity breaking: the critical phase transition
	* 3.2 Stability analysis via standard Hessian: the phase diagram
		* 3.2.1 Ergodic state: \(\bar{\mathbf{n}}=\bar{n}_{d,\rho,\beta}(0,\ldots,0)\)
		* 3.2.2 Pure state: \(\bar{\mathbf{n}}=\bar{n}_{d,\rho,\beta}(1,0,\ldots,0)\)
		* 3.2.3 Parallel state: \(\bar{\mathbf{n}}=\bar{n}_{d,\rho,\beta}(1,\ldots,1)\)
		* 3.2.4 Hierarchical state: \(\bar{\mathbf{n}}=\bar{n}_{d,\rho,\beta}((1-d),d(1-d),d^{2}(1-d),...)\)
	* 3.3 From the Cost function to the Loss function
* 4 Conclusions
* A A more general sampling scenario
* B On the data-set entropy \(\rho\)
* B.1 I: multitasking Hebbian network equipped with not-affecting-dilution noise
* B.2 II: multitasking Hebbian network equipped with not-preserving-dilution noise
* C Stability analysis: an alternative approach
* C.1 Stability analysis via signal-to-noise technique
* C.2 Evaluation of momenta of the effective post-synaptic potential
* D Explicit Calculations and Figures for the cases \(K=2\) and \(K=3\)
* D.1 \(K=2\)
* D.2 \(K=3\)
* E Proofs
* E.1 Proof of Theorem 1
* E.2 Proof of Proposition 1
Introduction

Typically, Artificial Intelligence has to deal with several inputs occurring at the same time: for instance, think about automatic driving, where it has to distinguish and react to different objects (e.g., pedestrians, traffic lights, riders, crosswalks) that may appear simultaneously. Likewise, when a biological neural network learns, it is rare that it has to deal with one single input per time1: for instance, while trained in the school to learn any single letter, we are also learning about the composition of our alphabets. In this perspective, when stating that neural networks operate _in parallel_, some caution in potential ambiguity should be paid. To fix ideas, let us focus on the Hopfield model [34], the _harmonic oscillator_ of associative neural networks accomplishing pattern recognition [8; 25]: its neurons indeed operate synergistically in parallel but with the purpose of retrieving one single pattern per time, not several simultaneously [8; 10; 36]. A parallel processing where multiple patterns are simultaneously retrieved cannot be accessible to the standard Hopfield networks as long as each pattern is fully informative, namely its vectorial binary representation is devoid of blank entries. On the other hand, when a fraction of entries can be blank [14] multiple-pattern retrieval is potentially achievable by the network. Intuitively, this can be explained by noticing that the overall number of neurons making up the networks - and thus available for information processing - equals the length of the binary vectors codifying the patterns to be retrieved, hence, as long as these vectors contain information in all their entries, there is no free room for dealing with multiple patterns. Conversely, the multitasking neural networks, introduced in [2], are able to overcome this limitation and have been shown to succeed in retrieving multiple patterns simultaneously, just by leveraging the presence of lacunae in the patterns stored by the network. The emerging pattern-recognition properties have been extensively investigated at medium storage (i.e., on random graphs above the percolation threshold) [23], at high storage (i.e., on random graphs below the percolation threshold) [24] as well as on scale-free [44] and hierarchical [3] topologies.

Footnote 1: It is enough to note that, should serial learning take place rather than parallel learning, Pavlov’s Classical Conditioning would not be possible [14].

However, while the study of the parallel retrieval capabilities of these multi-tasking networks is nowadays over, the comprehension of their parallel learning capabilities just started and it is the main focus of the present paper.

In these regards it is important to stress that the Hebbian prescription has been recently revised to turn it from a storing rule (built on a set of already definite patterns, as in the original Amit-Gutfreund-Sompolinksy (AGS) theory) into a genuine learning rule (where unknown patterns have to be inferred by experiencing solely a sample of their corrupted copies), see e.g., [5; 13; 27]2.

Footnote 2: While Statistical Learning theories appeared in the Literature a long time ago, see e.g. [1; 29; 43] for the original works and [6; 20; 26; 41] for updated references, however the statistical mechanics of Hebbian learning was not deepened in these studies.

In this work we merge these extensions of the bare AGS theory and use definite patterns (equipped with blank entries) to generate a sparse data-set of corrupted examples, that is the solely information experienced by the network: we aim to highlight the role of lacunae density and of the data-set size and quality on the network performance, in particular deepening the way the network learns simultaneously the patterns hidden behind the supplied examples. In this investigation we focus on the low-storage scenario (where the number of definite patterns grows sub-linearly with the volume of the network) addressing both the _supervised_ and the _unsupervised_ setting.

The paper is structured as follows: the main text has three Sections. Beyond this Introduction provided in Section 1, in Section 2 we revise the multi-tasking associative network; once briefly sum
marized its parallel retrieval capabilities (Sec. 2.1), we introduce a simple data-set the network has to cope with in order to move from the simpler storing of patterns to their learning from examples (Sec. 2.2). Next, in Section 3 we provide an exhaustive statistical mechanical picture of the network's emergent information processing capabilities by taking advantage of Guerra's interpolation techniques: in particular, focusing on the Cost function (Sec. 3.1), we face the _big-data_ limit (Sec. 3.1.1) and we deepen the nature of the phase transition the network undergoes as ergodicity breaking spontaneously takes place (Sec. 3.1.2). Sec. 3.2 is entirely dedicated to provide phase diagrams (namely plots in the space of the control parameters where different regions depict different global computational capabilities). Further, before reaching conclusions and outlooks as reported in Sec. 4, in Sec. 3.3 we show how the network's Cost function (typically used in Statistical Mechanics) can be sharply related to standard Loss functions (typically used in Machine Learning) to appreciate how parallel learning effectively lowers several Loss functions at once.

In the Appendices we fix a number of subtleties: in Appendix A we provide a more general setting for the sparse data-sets considered in this research3, while in Appendix B we inspect the relative entropies of these data-sets and finally in Appendix C we provide a revised version of the Signal-to-Noise technique (that allows to evaluate computational shortcuts beyond providing an alternative route to obtain the phase diagrams). Appendices D and E give details on calculations, plots and proofs of the main theorems.

Footnote 3: In the main text we face the simplest kind of pattern’s dilution, namely we just force to be blank the same fraction of their entries whose position is preserved in the generation of the data-sets (hence whenever the pattern has a zero, in all the examples it gives rise to, the zero will be kept), while in the appendix we relax this assumption (and blank entries can move along the examples yet preserving their amount). As in the thermodynamic limit the theory is robust w.r.t. these structural details we present as a main theme the simplest setting and in the appendix A the more cumbersome one.

## 2 Parallel learning in multitasking Hebbian neural networks

### A preliminary glance at the emergent parallel retrieval capabilities

Hereafter, for the sake of completeness, we briefly review the retrieval properties of the multitasking Hebbian network in the low-storage regime, while we refer to [2; 4] for an extensive treatment.

**Definition 1**.: _Given \(N\) Ising neurons \(\sigma_{i}=\pm 1\) (\(i=1,...,N\)), and \(K\) random patterns \(\mathbf{\xi}^{\mu}\) (\(\mu=1,...,K\)), each of length \(N\), whose entries are i.i.d. from_

\[\mathbb{P}(\xi_{i}^{\mu})=\frac{(1-d)}{2}\delta_{\xi_{i}^{\mu},-1}+\frac{(1-d) }{2}\delta_{\xi_{i}^{\mu},+1}+d\delta_{\xi_{i}^{\mu},0}, \tag{1}\]

_where \(\delta_{i,j}\) is the Kronecker delta and \(d\in[0,1]\), the Hamiltonian (or cost function) of the system reads as_

\[\mathcal{H}_{N}(\mathbf{\sigma}|\mathbf{\xi}):=-\frac{1}{2N}\sum_{\begin{subarray}{c}i,j\\ i\neq j\end{subarray}}^{N,N}\left(\sum_{\mu=1}^{K}\xi_{i}^{\mu}\xi_{j}^{\mu} \right)\sigma_{i}\sigma_{j}. \tag{2}\]

The parameter \(d\) tunes the "dilution" in pattern entries: if \(d=0\) the standard Rademacher setting of AGS theory is recovered, while for \(d=1\) no information is retained in these patterns: otherwise stated, these vectors display, on average, a fraction \(d\) of blank entries.
**Definition 2**.: _In order to assess the network retrieval performance we introduce the \(K\) Mattis magnetizations_

\[m_{\mu}:=\frac{1}{N}\sum_{i}^{N}\xi_{i}^{\mu}\sigma_{i},\ \mu=1,...,K, \tag{3}\]

_which quantify the overlap between the generic neural configuration \(\mathbf{\sigma}\) and the \(\mu^{th}\) pattern._

Note that the cost function (2) can be recast as a quadratic form in \(m_{\mu}\), namely

\[\mathcal{H}_{N}(\mathbf{\sigma}|\mathbf{\xi})=-\frac{N}{2}\sum_{\mu}m_{\mu}^{2}+\frac{ K}{2}, \tag{4}\]

where the term \(K/2\) in the r.h.s. stems from diagonal terms (\(i=j\)) in the sum at the r.h.s. of eq. 2 and in the low-load scenario (i.e., \(K\) grows sub-linearly with \(N\)) can be neglected in the thermodynamic limit (\(N\to\infty\)).

As we are going to explain, the dilution ruled by \(d\) is pivotal for the network in order to perform parallel processing. It is instructive to first consider a toy model handling just \(K=2\) patterns: let us assume, for simplicity, that the first pattern \(\mathbf{\xi}^{1}\) contains information (i.e., no blank entries) solely in the first half of its entries and the second pattern \(\mathbf{\xi}^{2}\) contains information solely in the second half of its entries, that is

\[\mathbf{\xi}^{1}=\underbrace{(\xi_{1}^{1},...,\xi_{N/2}^{1},0,...,0)}_{\in\{-1,+ 1\}^{\frac{N}{2}}},\quad\mathbf{\xi}^{2}=(\underbrace{0,...,0}_{\in\{0\}^{\frac{N }{2}}},\underbrace{\xi_{N/2+1}^{1},...,\xi_{N}^{1}}_{\in\{-1,+1\}^{\frac{N}{2} }}) \tag{5}\]

Unlike the standard Hopfield reference (\(d=0\)), where the retrieval of one pattern employs all the resources and there is no chance to retrieve any other pattern, not even partially (i.e., as \(m_{1}\to 1\) then \(m_{2}\approx 0\) because patterns are orthogonal for large \(N\) values in the standard random setting), here nor \(m_{1}\) neither \(m_{2}\) can reach the value \(1\) and therefore the complete retrieval of one of the two still leaves resources for the retrieval of the other. In this particular case, the minimization of the cost function \(\mathcal{H}_{N}(\mathbf{\sigma}|\mathbf{\xi})=-\frac{N}{2}\left(m_{1}^{2}+m_{2}^{2}\right)\) is optimal when _both_ the magnetizations are equal to one-half, that is when they both saturate their upper bound. In general, for arbitrary dilution level \(d\), the minimization of the cost function requires the network to be in one of the following regimes

* _hierarchical scenario_: for values of dilution not too high (i.e., \(d<d_{c}\), _vide infra_), one of the two patterns is fully retrieved (say \(m_{1}\approx 1-d\)) and the other is retrieved to the largest extent given the available resources, these being constituted by, approximately, the \(Nd\) neurons corresponding to the blank entries in \(\mathbf{\xi}^{1}\) (thus, \(m_{2}\approx d(1-d)\)), and so on if further patterns are considered.
* _parallel scenario_: for large values of dilution (i.e., above a critical threshold \(d_{c}\)), the magnetizations related to all the patterns raise and the signals they convey share the same amplitude.

In general, in this type of neural network, the _pure state ansatz4_\(\mathbf{m}=(1,0,0,...,0)\), that is \(\sigma_{i}=\xi_{i}^{1}\) for \(i=1,...,N\), barely works and parallel retrieval is often favored. In fact, for \(K\geq 2\), at relatively low values of pattern dilution \(d_{1}\) and in the zero-noise limit \(\beta\to\infty\), one can prove the validity of the so-called _hierarchical ansatz_[2] as we briefly discuss: one pattern, say \(\mathbf{\xi}^{1}\), is perfectly retrieved and displays a Mattis magnetization \(m^{1}\approx(1-d)\); a fraction \(d\) of neurons is not involved and is therefore available for further retrieval, with any remaining pattern, say \(\mathbf{\xi}^{2}\), which yields \(m_{2}\sim(1-d)d\); proceeding iteratively, one finds \(m_{\ell}=d^{\ell-1}(1-d)\) for \(\ell=1,...,\hat{K}\) and the overall number \(\hat{K}\) of patterns
simultaneously retrieved corresponds to the employment of all the resources. Specifically, \(\hat{K}\) can be estimated by setting \(\sum_{\ell=0}^{\hat{K}-1}(1-d)d^{\ell}=1\), with the cut-off at finite \(N\) as \((1-d)d^{\hat{K}-1}\geq N^{-1}\), due to discreteness: for any fixed and finite \(d\), this implies \(K\lesssim\log N\), which can be thought of as a "parallel low-storage" regime of neural networks. It is worth stressing that, in the above mentioned regime of low dilution, the configuration leading to \(m_{\ell}=d^{\ell-1}(1-d)\) for \(\ell=1,...,\hat{K}\) is the one which minimizes the cost function. The hierarchical retrieval state \(\mathbf{m}=(1-d)\left(1,d,d^{2},d^{3},\cdots\right)\) can also be specified in terms of neural configuration as [2]

\[\sigma_{i}^{*}=\xi_{i}^{1}+\sum_{\nu=2}^{\hat{K}}\xi_{i}^{\nu}\prod_{\rho=1}^{ \nu-1}\delta_{\xi_{i}^{\rho},0}\,. \tag{6}\]

This organization is stable until a critical dilution level \(d_{c}\) is reached where \(m_{1}\sim\sum_{k>1}m_{k}\)[2], beyond that level the network undergoes a rearrangement and a new organization called _parallel ansatz_ supplants the previous one. Indeed for high values of dilution (i.e \(d\to 1\)) it is immediate to check that the ratio among the various intensities of all the magnetizations stabilizes to the value one, i.e. \((m_{k}/m_{k-1})\sim d^{k-1}(1-d)/d^{k-2}(1-d)\to 1\), hence in this regime all the magnetizations are raised with the same strength and the network is operationally set in a fully parallel retrieval mode: the parallel retrieval state simply reads \(\mathbf{m}=(\bar{m})\left(1,1,1,1,\cdots\right)\). This picture is confirmed by the plots shown in Fig. 1 and obtained by solving the self-consistency equations for the Mattis magnetizations related to the multitasking Hebbian network equipped with \(K=2\) patterns that read as [2]

\[m_{1} = d(1-d)\tanh(\beta m_{1})+\frac{(1-d)^{2}}{2}\left\{\tanh[\beta(m _{1}+m_{2})]+\tanh[\beta(m_{1}-m_{2})]\right\}, \tag{7}\] \[m_{2} = d(1-d)\tanh(\beta m_{1})+\frac{(1-d)^{2}}{2}\left\{\tanh[\beta(m _{1}+m_{2})]-\tanh[\beta(m_{1}-m_{2})]\right\} \tag{8}\]

where \(\beta\in\mathbb{R}^{+}\) denotes the level of noise.

We remark that these hierarchical or parallel organizations of the retrieval, beyond emerging naturally

Figure 1: Numerical solutions of the two self-consistent equations (7) and (8) obtained for \(K=2\), see [2], as a function of \(d\) and for different choices of \(\beta\): in the \(d\to 0\) limit the Hopfield serial retrieval is recovered (one magnetization with intensity one and the other locked at zero), for \(d\to 1\) the network ends up in the parallel regime (where all the magnetizations acquire the same value), while for intermediate values of dilution the hierarchical ordering prevails (both the magnetizations are raised, but their amplitude is different).
within the equilibrium description provided by Statistical Mechanics, are actually the real stationary states of the dynamics of these networks at work with diluted patterns as shown in Figure 2.

### From parallel storing to parallel learning

In this section we revise the multitasking Hebbian network [2; 4] in such a way that it can undergo a _learning_ process instead of a simple _storing_ of patterns. In fact, in the typical learning setting, the set of definite patterns, hereafter promoted to play as "archetypes", to be reconstructed by the network is not available, rather, the network is exposed to examples, namely noisy versions of these archetypes.

As long as enough examples are provided to the network, this is expected to correctly form its own representation of the archetypes such that, in further expositions to a new example related to a certain archetype, it will be able to retrieve it and, since then, suitably generalize it.

This generalized Hebbian kernel has recently been introduced to encode unsupervised [5] and supervised [13] learning processes and, in the present paper, these learning rules are modified in order to deal with diluted patterns.

First, let us define the data-set these networks have to cope with: the archetypes are randomly drawn from the distribution (1). Each archetype \(\mathbf{\xi}^{\mu}\) is then used to generate a set of \(M_{\mu}\) perturbed

Figure 2: We report two examples of Monte Carlo dynamics until thermalization within the hierarchical (upper plots, dilution level \(d=0.2\)) and parallel (lower plots, dilution level \(d=0.8\)) scenarios respectively. These plots confirm that the picture provided by statistical mechanics is actually dynamically reached by the network. We initialize the network sharply in a pattern as a Cauchy condition (represented as the dotted blue Dirac delta peaked at the pattern in the second columns) and, in the first column, we show the stationary values of the various Mattis magnetizations pertaining to different patterns, while in the second column we report their histograms achieved by sampling 1000 independent Monte Carlo simulations: starting from a sequential retrieval regime, the network ends up in a multiple retrieval mode, hierarchical vs parallel depending on the level of dilution in the patterns.
versions, denoted as \(\mathbf{\eta}^{\mu,a}\) with \(a=1,...,M_{\mu}\) and \(\mathbf{\eta}^{\mu,a}\in\{-1,0,+1\}^{N}\). Thus, the overall set of examples to be supplied to the network is given by \(\mathbf{\eta}=\{\mathbf{\eta}^{\mu,a}\}_{\mu=1,...,K}^{a=1,...,M_{\mu}}\). Of course, different ways to sample examples are conceivable: for instance, one can require that the position of blank entries appearing in \(\mathbf{\xi}^{\mu}\) is preserved over all the examples \(\{\mathbf{\eta}^{\mu,a}\}_{a=1,...,M_{\mu}}\), or one can require that only the number of blank entries \(\sum_{i=1}^{N}\delta_{\xi_{i}^{\mu},0}\) is preserved (either strictly or in the average). Here we face the first case because it requires a simpler notation, but we refer to Appendix A for a more general treatment.

**Definition 3**.: _The entries of each examples are depicted following_

\[\mathbb{P}(\eta_{i}^{\mu,a}|\xi_{i}^{\mu})=\frac{1+r_{\mu}}{2}\delta_{\eta_{i }^{\mu,a},\xi_{i}^{\mu}}+\frac{1-r_{\mu}}{2}\delta_{\eta_{i}^{\mu,a},-\xi_{i}^ {\mu}}, \tag{9}\]

_for \(i=1,\ldots,N\) and \(\mu=1,\ldots,K\). Notice that \(r_{\mu}\) tunes the data-set quality: as \(r_{\mu}\to 1\) examples belonging to the \(\mu\)-th set collapse on the archetype \(\mathbf{\xi}^{\mu}\), while as \(r_{\mu}\to 0\) examples turn out to be uncorrelated with the related archetype \(\mathbf{\xi}^{\mu}\)._

As we will show in the next sections, the behavior of the system depends on the parameters \(M_{\mu}\) and \(r_{\mu}\) only through the combination \(\frac{1-r_{\mu}^{2}}{M_{\mu}r_{\mu}^{2}}\), therefore, as long as the ratio \(\frac{1-r_{\mu}^{2}}{M_{\mu}r_{\mu}^{2}}\) is \(\mu\)-independent, the theory shall not be affected by the specific choice of the archetype. Thus, for the sake of simplicity, hereafter we will consider \(r\) and \(M\) independent of \(\mu\) and we will pose \(\rho:=\frac{1-r^{2}}{Mr^{2}}\). Remarkably, \(\rho\) plays as an information-content control parameter [13]: to see this, let us focus on the \(\mu\)-th pattern and \(i\)-th digit, whose related block is \(\mathbf{\eta}_{i}^{\mu}=(\eta_{i}^{\mu},\eta_{i}^{\mu,2},\ldots,\eta_{i}^{\mu,M})\), the error probability for any single entry is \(\mathcal{P}(\xi_{i}^{\mu}\neq 0)\mathcal{P}(\eta_{i}^{\mu,a}\neq\xi_{i}^{\mu})=(1-d)( 1-r_{\mu})/2\) and, by applying the majority rule on the block, we get \(\mathcal{P}(\xi_{i}^{\mu}\neq 0)\mathcal{P}(\mathrm{sign}(\sum\limits_{a}\eta_{i}^{ \mu,a})\xi_{i}^{\mu}=-1)\underset{M\gg 1}{\approx}\frac{(1-d)}{2}\left[1-\mathrm{ erf}\left(1/\sqrt{2\rho}\right)\right]\) thus, by computing the conditional entropy \(H_{d}(\xi_{i}^{\mu}|\mathbf{\eta}_{i}^{\mu})\), that quantifies the amount of information needed to describe the original message \(\xi_{i}^{\mu}\) given the related block \(\mathbf{\eta}_{i}^{\mu}\), we get

\[H_{d}(\xi_{i}^{\mu}|\mathbf{\eta}_{i}^{\mu}) = -\left[\frac{1+d}{2}+\frac{1-d}{2}\mathrm{erf}\left(\frac{1}{ \sqrt{2\rho}}\right)\right]\ \log\left[\frac{1+d}{2}+\frac{1-d}{2}\mathrm{erf}\left(\frac{1}{ \sqrt{2\rho}}\right)\right]\] \[-\left[\frac{1-d}{2}-\frac{1-d}{2}\mathrm{erf}\left(\frac{1}{ \sqrt{2\rho}}\right)\right]\ \log\left[\frac{1-d}{2}-\frac{1-d}{2}\mathrm{erf}\left(\frac{1}{ \sqrt{2\rho}}\right)\right]\]

which is monotonically increasing with \(\rho\). Therefore, with a slight abuse of language, in the following \(\rho\) shall be referred to as _data-set entropy_.

The available information is allocated directly in the synaptic coupling among neurons (as in the standard Hebbian storing), as specified by the following supervised and unsupervised generalization of the multitasking Hebbian network:

**Definition 4**.: _Given \(N\) binary neurons \(\sigma_{i}=\pm 1\), with \(i\in(1,...,N)\), the cost function (or Hamiltonian) of the multitasking Hebbian neural network in the supervised regime is_

\[\mathcal{H}_{N,K,d,M,r}^{(sup)}(\mathbf{\sigma}|\mathbf{\eta})=-\frac{1}{2N}\frac{1}{ (1-d)(1+\rho)}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\left(\frac{1}{Mr}\sum_{a=1}^{ M}\eta_{i}^{\mu,a}\right)\left(\frac{1}{Mr}\sum_{b=1}^{M}\eta_{j}^{\mu,b}\right) \sigma_{i}\sigma_{j}. \tag{11}\]

**Definition 5**.: _Given \(N\) binary neurons \(\sigma_{i}=\pm 1\), with \(i\in(1,...,N)\), the cost function (or Hamiltonian) of the multitasking Hebbian neural network in the unsupervised regime is_

\[\mathcal{H}_{N,K,d,M,r}^{(unsup)}(\mathbf{\sigma}|\mathbf{\eta})=-\frac{1}{2N}\frac{1}{ (1-d)(1+\rho)}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\left(\frac{1}{Mr^{2}}\sum_{a =1}^{M}\eta_{i}^{\mu,a}\eta_{j}^{\mu,a}\right)\sigma_{i}\sigma_{j}. \tag{12}\]
**Remark 1**.: _The factor \((1-d)(1+\rho)\) appearing in (2.11) corresponds to \(\mathbb{E}_{\xi},\mathbb{E}_{(\eta|\xi)}\left[\sum\limits_{a}\eta_{i}^{\mu,a}/( Mr)\right]^{2}\) and it plays as a normalization factor. A similar factor is also inserted in (2.12)._

**Remark 2**.: _By direct comparison between (2.11) and (2.12), the role of the "teacher" in the supervised setting is evident: in the unsupervised scenario, the network has to handle all the available examples regardless of their archetype label, while in the supervised counterpart a teacher has previously grouped examples belonging to the same archetype together (whence the double sum on \(a=(1,...,M)\) and on \(b=(1,...,M)\) appearing in eq. (2.11), that is missing in eq. (2.12))._

We investigate the model within a canonical framework: we introduce the Boltzmann-Gibbs measure

\[\mathcal{P}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\sigma}|\mathbf{\eta}):=\frac{ \exp[-\beta\mathcal{H}^{(sup,unsup)}_{N,K,d,M,r}(\mathbf{\sigma}|\mathbf{\eta})]}{ \mathcal{Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta})}, \tag{2.13}\]

where

\[\mathcal{Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta}):=\sum_{\mathbf{\sigma}} \exp\left[-\beta\mathcal{H}^{(sup,unsup)}_{N,K,d,M,r}(\mathbf{\sigma}|\mathbf{\eta})\right] \tag{2.14}\]

is the normalization factor, also referred to as partition function, and the parameter \(\beta\in\mathbb{R}^{+}\), rules the broadness of the distribution in such a way that for \(\beta\to 0\) (infinite noise limit) all the \(2^{N}\) neural configurations are equally likely, while for \(\beta\to\infty\) the distribution is delta-peaked at the configurations corresponding to the minima of the Cost function.

The average performed over the Boltzmann-Gibbs measure is denoted as

\[\omega^{(sup,unsup)}_{N,K,\beta,d,M,r}[\cdot]:=\sum_{\mathbf{\sigma}}^{2^{N}}\, \cdot\,\mathcal{P}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\sigma}|\mathbf{\eta}). \tag{2.15}\]

Beyond this average, we shall also take the so-called _quenched_ average, that is the average over the realizations of archetypes and examples, namely over the distributions (2.1) and (2.9), and this is denoted as

\[\mathbb{E}[\cdot]=\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}[\cdot]. \tag{2.16}\]

**Definition 6**.: _The quenched statistical pressure of the network at finite network size \(N\) reads as_

\[\mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}=\frac{1}{N}\mathbb{E}\log\mathcal{ Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta}). \tag{2.17}\]

_In the thermodynamic limit we pose_

\[\mathcal{A}^{(sup,unsup)}_{K,\beta,d,M,r}=\lim_{N\to\infty}\mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}. \tag{2.18}\]

_We recall that the statistical pressure equals the free energy times \(-\beta\) (hence they convey the same information content)._

**Definition 7**.: _The network capabilities can be quantified by introducing the following order parameters, for \(\mu=1,\ldots,K\),_

\[m_{\mu} :=\frac{1}{N}\sum_{i=1}^{N}\xi_{i}^{\mu}\sigma_{i},\] \[n_{\mu,a} :=\frac{1}{(1+\rho)r}\frac{1}{N}\sum_{i=1}^{N}\eta_{i}^{\mu,a} \sigma_{i}, \tag{2.19}\] \[n_{\mu} :=\frac{1}{M}\sum_{a=1}^{M}n_{\mu,a}=\frac{1}{(1+\rho)r}\frac{1}{ NM}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu,a}\sigma_{i},\]
We stress that, beyond the fairly standard \(K\) Mattis magnetizations \(m_{\mu}\), which assess the alignment of the neural configuration \(\mathbf{\sigma}\) with the archetype \(\mathbf{\xi}^{\mu}\), we need to introduce also \(K\) empirical Mattis magnetizations \(n_{\mu}\), which compare the alignment of the neural configuration with the average of the examples labelled with \(\mu\), as well as \(K\times M\) single-example Mattis magnetizations \(n_{\mu,a}\), which measure the proximity between the neural configuration and a specific example. An intuitive way to see the suitability of the \(n_{\mu}\)'s and of the \(n_{\mu,a}\)'s is by noticing that the cost functions \(\mathcal{H}^{(sup)}\) and \(\mathcal{H}^{(unsup)}\) can be written as a quadratic form in, respectively, \(n_{\mu}\) and \(n_{\mu,a}\); on the other hand, the \(m_{\mu}\)'s do not appear therein explicitly as the archetypes are unknowns in principle.

Finally, notice that no spin-glass order parameters is needed here (since we are working in the low-storage regime [8; 25]).

## 3 Parallel Learning: the picture by statistical mechanics

### Study of the Cost function and its related Statistical Pressure

To inspect the emergent capabilities of these networks, we need to estimate the order parameters introduced in Equations (19) and analyze their behavior versus the control parameters \(K,\beta,d,M,r\). To this task we need an explicit expression of the statistical pressure in terms of these order parameters so to extremize the former over the latter. In this Section we carry on this investigation in the thermodynamic limit and in the low storage scenario by relying upon Guerra's interpolating techniques (see e.g., [30; 17; 31; 32]): the underlying idea is to introduce an interpolating statistical pressure whose extrema are the original model (which is the target of our investigation but we can be not able to address it directly) and a simple one (which is usually a one-body model that we can solve exactly). We then start by evaluating the solution of the latter and next we propagate the obtained solution back to the original model by the fundamental theorem of calculus, integrating on the interpolating variable. Usually, in this last passage, one assumes replica symmetry, namely that the order-parameter fluctuations are negligible in the thermodynamic limit as this makes the integral propagating the solution analytical. In the low-load scenario replica symmetry holds exactly, making the following calculation rigorous. In fact, as long as \(K/N\to 0\) while \(N\to\infty\), the order parameters self-average around their means [19; 46], that will be denoted by a bar, that is

\[\lim_{N\to\infty}\mathcal{P}_{N,K,\beta,d,M,r}(m_{\mu}) = \delta\left(m_{\mu}-\bar{m}_{\mu}\right),\quad\forall\mu\in(1,...,K), \tag{22}\] \[\lim_{N\to\infty}\mathcal{P}_{N,K,\beta,d,M,r}(n_{\mu}) = \delta\left(n_{\mu}-\bar{n}_{\mu}\right),\quad\forall\mu\in(1,...,K), \tag{23}\]

where \(\mathcal{P}_{N,K,\beta,d,M,r}\) denotes the Boltzmann-Gibbs probability distribution for the observables considered. We anticipate that the mean values of these distributions are independent of the training (either supervised or unsupervised) underlying the Hebbian kernel.

Before proceeding, we slightly revise the partition functions (14) by inserting an extra term in their exponents because it allows to apply the functional generator technique to evaluate the Mattis magnetizations. This implies the following modification, respectively in the supervised and unsupervised settings, of the partition function

**Definition 8**.: _Given the interpolating parameter \(t\in[0,1]\), the auxiliary field \(J\) and the constants \(\{\psi_{\mu}\}_{\mu=1,...,K}\in\mathbb{R}\) to be set a posteriori, Guerra's interpolating partition function for the supervised
_and unsupervised multitasking Hebbian networks is given, respectively, by_

\[\mathcal{Z}^{(sup)}_{N,K,\beta,d,M,r}(\mathbf{\eta};J,t)=\sum_{\{\mathbf{ \sigma}\}}\int\,d\mu(z_{\mu})\exp\Bigg{[}J\sum_{\mu,i}\xi_{i}^{\mu}\sigma_{i}+ \frac{t\beta N(1+\rho)}{2(1-d)}\sum_{\mu}n_{\mu}^{2}(\mathbf{\sigma})+(1-t)\frac{N }{2}\sum_{\mu}\psi_{\mu}\,n_{\mu}(\mathbf{\sigma})\Bigg{]}. \tag{10}\] \[\mathcal{Z}^{(unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta};J,t)=\sum_{\{ \mathbf{\sigma}\}}\int\,d\mu(z_{\mu})\exp\Bigg{[}J\sum_{\mu,i}\xi_{i}^{\mu}\sigma_{ i}+\frac{t\beta N(1+\rho)}{2(1-d)M}\sum_{\mu=1}^{K}\sum_{a=1}^{M}n_{\mu,a}^{2}( \mathbf{\sigma})+(1-t)N\sum_{\mu,a}\psi_{\mu}\,n_{\mu,a}(\mathbf{\sigma})\Bigg{]}. \tag{11}\]

More precisely, we added the term \(J\sum_{\mu}\sum_{i}\xi_{i}^{\mu}\sigma_{i}\) that allows us to to "generate" the expectation of the Mattis magnetization \(m_{\mu}\) by evaluating the derivative w.r.t. \(J\) of the quenched statistical pressure at \(J=0\). This operation is not necessary for _Hebbian storing_, where the Mattis magnetization is a natural order parameter (the Hopfield Hamiltonian can be written as a quadratic form in \(m_{\mu}\), as standard in AGS theory [8]), while for _Hebbian learning_ (whose cost function can be written as a quadratic form in \(n_{\mu}\), not in \(m_{\mu}\) as the network does not experience directly the archetypes) we need such a term for otherwise the expectation of the Mattis magnetization would not be accessible. This operation gets redundant in the \(M\to\infty\) limit, where \(m_{\mu}\) and \(n_{\mu}\) become proportional by a standard Central Limit Theorem (CLT) argument (see also Sec. 3.1.1 and [13]). Clearly, \(\mathcal{Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta})=\lim_{J\to 0} \mathcal{Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta};J)\) and these generalized interpolating partition functions, provided in eq.s (10) and (11) respectively, recover the original models when \(t=1\), while they return a simple one-body model at \(t=0\).

Conversely, the role of the \(\psi_{\mu}\)'s is instead that of mimicking, as close as possible, the true post-synaptic field perceived by the neurons.

These partition functions can be used to define a generalized measure and a generalized Boltzmann-Gibbs average that we indicate by \(\omega_{t}^{(sup,unsup)}[\cdot]\). Of course, when \(t=1\) the standard Boltzmann-Gibbs measure and related averages are recovered.

Analogously, we can also introduce a generalized interpolating quenched statistical pressures as

**Definition 9**.: _The interpolating statistical pressure for the multitasking Hebbian neural network is introduced as_

\[\mathcal{A}^{(sup,unsup)}_{N,K\beta,d,M,r}(J,t)\coloneqq\frac{1}{N}\mathbb{E }\left[\ln\mathcal{Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(\mathbf{\eta};J,t)\right], \tag{12}\]

_and, in the thermodynamic limit,_

\[\mathcal{A}^{(sup,unsup)}_{K,\beta,d,M,r}(J,t)\coloneqq\lim_{N\to\infty} \mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}(J,t). \tag{13}\]

_Obviously, by setting \(t=1\) in the interpolating pressures we recover the original ones, namely \(\mathcal{A}^{(sup,unsup)}_{K,\beta,d,M,r}(J)=\mathcal{A}^{(sup,unsup)}_{K, \beta,d,M,r}(J,t=1)\), which we finally evaluate at \(J=0\)._

We are now ready to state the next

**Theorem 1**.: _In the thermodynamic limit (\(N\to\infty\)) and in the low-storage regime (\(K/N\to 0\)), the quenched statistical pressure of the multitasking Hebbian network - trained under supervised or unsupervised learning - reads as_

\[\mathcal{A}^{(sup,unsup)}_{K,\beta,d,M,r}(J) = \mathbb{E}\left\{\ln\left[2\cosh\left(J\sum_{\mu=1}^{K}\xi^{\mu} +\frac{\beta}{1-d}\sum_{\mu=1}^{K}\bar{n}_{\mu}\tilde{\eta}^{\mu}\right)\, \right]\right\}-\frac{\beta}{1-d}(1+\rho)\sum_{\mu=1}^{K}\bar{n}_{\mu}^{2}. \tag{14}\]
_where \(\mathbb{E}=\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\), \(\hat{\eta}^{\mu}=\frac{1}{Mr}\sum_{a=1}^{M}\eta_{i}^{\mu,a}\), and the values \(\bar{n}_{\mu}\) must fulfill the following self-consistent equations_

\[\bar{n}_{\mu}=\frac{1}{(1+\rho)}\mathbb{E}\left\{\tanh\left[\frac{\beta}{(1-d) }\sum_{\nu=1}^{K}\bar{n}_{\nu}\hat{\eta}^{\nu}\right]\hat{\eta}^{\mu}\right\},\quad\forall\mu\in(1,...,K), \tag{3.8}\]

_as these values of the order parameters are extremal for the statistical pressure \(\mathcal{A}_{K,\beta,d,M,r}^{(sup,unsup)}(J=0)\)._

**Corollary 1**.: _By considering the auxiliary field \(J\) coupled to \(m_{\mu}\) and recalling that \(\lim_{N\to\infty}m_{\mu}=\bar{m}_{\mu}\), we can write down a self-consistent equation also for the Mattis magnetization as \(\bar{m}_{\mu}=\partial_{J}\mathcal{A}_{K,\beta,d,M,r}^{(sup,unsup)}(J)_{|J=0}\), thus we have_

\[\bar{m}_{\mu}=\mathbb{E}\left\{\tanh\left[\frac{\beta}{(1-d)}\sum_{\nu=1}^{K }\bar{n}_{\nu}\hat{\eta}^{\nu}\right]\xi^{\mu}\right\},\quad\forall\mu\in(1,...,K). \tag{3.9}\]

For the proof of Proposition 1 and of Corollary 1 we refer to Appendix E.1.

We highlighted that the expressions of the quenched statistical pressure for a network trained with or without the supervision of a teacher do actually coincide: intuitively, this happens because we are considering only a few archetypes (i.e. we work at low load), consequently, the minima of the cost function are well separated and there is only a negligible role of the teacher in shaping the landscape to avoid overlaps in their basins of attractions. Clearly, this is expected to be no longer true in the high load setting and, indeed, it is proven not to hold for non-diluted patterns, where supervised and unsupervised protocols give rise to different outcomes [5, 13]. By a mathematical perspective the fact that, whatever the learning procedure, the expression of the quenched statistical pressure is always the same, is a consequence of standard concentration of measure arguments [17, 45] as, in the \(N\to\infty\) limit, beyond eq. (3.2), it is also \(\mathcal{P}(n_{\mu,a})\to\delta(n_{\mu,a}-\bar{n}_{\mu})\).

The self-consistent equations (3.9) have been solved numerically for several values of parameters and results for K=2 and for \(K=3\) are shown in Fig. 3 (where also the values of the cost function is reported) and Fig. 4 respectively. We also checked the validity of these results by comparing them with the outcomes of Monte Carlo simulations, finding an excellent asymptotic agreement; further, in the large \(M\) limit, the magnetizations eventually converge to the values predicted by the theory developed in the storing framework, see eq. (2.6). Therefore, in both the scenarios, the hierarchical or parallel organization of the magnetization's amplitudes are recovered: beyond the numerical evidence just mentioned, in Appendix D an analytical proof is provided.

#### 3.1.1 Low-entropy data-sets: the Big-Data limit

As discussed in Sec. 2.2, the parameter \(\rho\) quantifies the amount of information needed to describe the original message \(\mathbf{\xi}^{\mu}\) given the set of related examples \(\{\mathbf{\eta}^{\mu,a}\}_{a=1,...,M}\). In this section we focus on the case \(\rho\ll 1\) that corresponds to a highly-informative data-set; we recall that in the limit \(\rho\to 0\) we get a data-set where either the items (\(r\to 1\)) or their empirical average (\(M\to\infty\), \(r\) finite) coincide with the archetypes, in such a way that the theory collapses to the standard Hopfield reference.

As explained in Appendix E.2, we start from the self-consistent equations (3.8)-(3.9) and we exploit the Central Limit Theorem to write \(\hat{\eta}^{\mu}\sim\xi^{\mu}\left(1+\lambda_{\mu}\sqrt{\rho}\right)\), where \(\lambda_{\mu}\sim\mathcal{N}(0,1)\). In this way we reach the simpler expressions given by the next

**Proposition 1**.: _In the low-entropy data-set scenario, preserving the low storage and thermodynamic limit assumptions, the two sets of order parameters of the theory, \(\bar{m}_{\mu}\) and \(\bar{n}_{\mu}\) become related by the
following equations_

\[\bar{n}_{\mu} = \frac{\bar{m}_{\mu}}{(1+\rho)}+\beta^{{}^{\prime}}\frac{\rho\,\bar{ n}_{\mu}}{(1+\rho)}\mathbb{E}_{\xi,Z}\left\{\left[1-\tanh^{2}\left(g(\beta,Z,\bar{ \boldsymbol{n}})\right)\right](\xi^{\mu})^{2}\right\}, \tag{3.10}\] \[\bar{m}_{\mu} = \mathbb{E}_{\xi,Z}\left\{\tanh\left[g(\beta,\boldsymbol{\xi},Z, \bar{\boldsymbol{n}})\right]\xi^{\mu}\right\}, \tag{3.11}\]

_where_

\[g(\beta,\boldsymbol{\xi},Z,\bar{\boldsymbol{n}})=\beta^{{}^{\prime}}\sum_{ \nu=1}^{K}\bar{n}_{\nu}\xi^{\nu}+\beta^{{}^{\prime}}\,Z\sqrt{\rho\,\sum_{\nu= 1}^{K}\bar{n}_{\nu}^{2}\left(\xi^{\nu}\right)^{2}} \tag{3.12}\]

_and \(Z\sim\mathcal{N}(0,1)\) is a standard Gaussian variable. Furthermore, to lighten the notation and assuming \(d\neq 1\) with no loss of generality, we posed_

\[\beta^{{}^{\prime}}=\frac{\beta}{1-d}. \tag{3.13}\]

The regime \(\rho\ll 1\), beyond being an interesting one (e.g., it can be seen as a _big data_\(M\to\infty\) limit of the theory), offers a crucial advantage because of the above emerging proportionality relation between \(\bar{n}\) and \(\bar{m}\) (see eq. 3.10). In fact, the model is supplied only with examples - upon which the \(n_{\mu}\)'s are defined - while it is not aware of archetypes - upon which the \(m_{\mu}\)'s are defined - yet we can

Figure 3: Snapshots of cost function (upper plots) -where we use the label \(E\) for _energy_- and magnetizations (lower plots) for data-sets generated by \(K=2\) archetypes and at different entropies, in the noiseless limit \(\beta\to\infty\). Starting by \(\rho=0.0\) we see that the hierarchical regime (black lines) dominates at relatively mild dilution values (i.e., the energy pertaining to this configuration is lower w.r.t. the parallel regime), while for \(d\to 1\) the hierarchical ordering naturally collapse to the parallel regime (red lines), where all the magnetizations acquire the same values. Further note how, by increasing the entropy in the data-set (e.g. for \(\rho=0.1\) and \(\rho=0.4\)), the domain of validity of the parallel regime enlarges (much as increasing \(\beta\) in the network, see Fig. 1). The vertical blue lines mark the transitions between these two regimes as captured by Statistical Mechanics: it corresponds to switching from the white to the green regions of the phase diagrams of Fig. 6.
use this relation to recast the self-consistent equation for \(\bar{n}\) into a self-consistent equation for \(\bar{m}\) such that its numerical solution in the space of the control parameters allows us to get the phase diagram of such a neural network more straightforwardly.

Further, we can find out explicitly the thresholds for learning, namely the minimal amount of examples (given the level of noise \(r\), the amount of archetype to handle \(K\), etc.) that guarantee that the network can safely infer the archetype from the supplied data-set. To obtain these thresholds we have to deepen the ground state structure of the network, that is, we now handle the Eqs. (20)-(21) to compute their zero fast-noise limit (\(\beta\to\infty\)). As detailed in the Appendix E.2 (see Corollary 3), by taking the limit \(\beta\to\infty\) in eqs. (20)-(21) we get

\[\bar{m}_{\mu}\,=\,\mathbb{E}_{\xi}\left\{\mathrm{erf}\left[\left(\sum_{\nu=1 }^{K}\bar{m}_{\nu}\xi^{\nu}\right)\left(2\rho\sum_{\nu=1}^{K}\bar{m}_{\nu}^{2} \left(\xi^{\nu}\right)^{2}\right)^{-1/2}\right]\xi^{\mu}\right\}\,. \tag{22}\]

Once reached a relatively simple expression for \(\bar{m}_{\mu}\), we can further manipulate it and try to get information about the existence of a lower-bound value for \(M\), denoted with \(M_{\otimes}\), which ensures that the network has been supplied with sufficient information to learn and retrieve the archetypes.

Figure 4: _Behaviour of the Mattis magnetizations as more and more examples are supplied to the network._ Monte Carlo numerical checks (colored dots, \(N=6000\)) for a diluted network with \(r=0.1\) and \(K=3\) are in plain agreement with the theory: solutions of the self-consistent equation for the Mattis magnetizations reported in the Corollary 1 are shown as solid lines. As dilution increases, the network behavior departs from a Hopfield-like retrieval (\(d=0.1\)) where just the blue magnetization is raised (serial pattern recognition) to the hierarchical regime (\(d=0.25\) and \(d=0.55\)) where multiple patterns are simultaneously retrieved with different amplitudes, while for higher values of dilution the network naturally evolves toward the parallel regime (\(d=0.75\)) where all the magnetizations are raised and with the same strength. Note also the asymptotic agreement with the dotted lines, whose values are those predicted by the multitasking Hebbian storage [2].
Setting \(\beta\to\infty\), we expect that the magnetizations fulfill the hierarchical organization, namely \((\bar{m}_{1},\bar{m}_{2},\ldots)=(1-d)(1,d,\ldots)\) and (3.14) becomes

\[\bar{m}_{\mu}\sim\frac{1-d}{2}\mathbb{E}_{\xi^{\nu\neq\mu}}\left\{\mathrm{erf} \left[\frac{d^{\mu}+\sum\limits_{\nu\neq\mu}^{K}d^{\nu}\xi^{\nu}}{\sqrt{2\rho} \sqrt{d^{2\mu}+\sum\limits_{\nu\neq\mu}^{K}d^{2\nu}\left(\xi^{\nu}\right)^{2}} }\right]+\mathrm{erf}\left[\frac{d^{\mu}-\sum\limits_{\nu\neq\mu}^{K}d^{\nu} \xi^{\nu}}{\sqrt{2\rho}\sqrt{d^{2\mu}+\sum\limits_{\nu\neq\mu}^{K}d^{2\nu} \left(\xi^{\nu}\right)^{2}}}\right]\right\}\,, \tag{3.15}\]

where we highlighted that the expectation is over all the archetypes but the \(\mu\)-th one under inspection.

Next, we introduce a confidence interval, ruled by \(\Theta\), and we require that

\[\bar{m}_{\mu}>(1-d)d^{\mu-1}\mathrm{erf}\left[\Theta\right]. \tag{3.16}\]

In order to quantify the critical number of examples \(M_{\otimes}^{\mu}\) needed for a successful learning of the archetype \(\mu\) we can exploit the relation

\[\mathbb{E}_{\xi^{\nu\neq\mu}}\Big{\{}\mathrm{erf}\Big{[}f(\xi)\Big{]}\Big{\}} \geq\min_{\xi^{\nu\neq\mu}}\!\Big{\{}\mathrm{erf}\Big{[}f(\xi)\Big{]}\Big{\}}\,, \tag{3.17}\]

where in our case

\[\begin{split}\min_{\xi^{\nu\neq\mu}}\!\Big{\{}\mathrm{erf} \Big{[}f(\xi)\Big{]}\Big{\}}&=\min_{\xi^{\nu\neq\mu}}\left\{ \mathrm{erf}\left[\frac{d^{\mu}+\sum\limits_{\nu\neq\mu}^{K}d^{\nu}\xi^{\nu}}{ \sqrt{2\rho}\sqrt{d^{2\mu}+\sum\limits_{\nu\neq\mu}^{K}d^{2\nu}\left(\xi^{\nu }\right)^{2}}}\right]+\mathrm{erf}\left[\frac{d^{\mu}-\sum\limits_{\nu\neq\mu}^{ K}d^{\nu}\xi^{\nu}}{\sqrt{2\rho}\sqrt{d^{2\mu}+\sum\limits_{\nu\neq\mu}^{K}d^{2\nu} \left(\xi^{\nu}\right)^{2}}}\right]\right\}\\ &=\,2\,\mathrm{erf}\left[\left(d^{\mu}-\sum\limits_{\nu\neq\mu}^{ K}d^{\nu}\right)\left(2\rho\sum\limits_{\nu=1}^{K}d^{2\nu}\right)^{-1/2}\right]. \end{split} \tag{3.18}\]

Thus, using the previous relation in (3.16), the following inequality must hold

\[\mathrm{erf}\left[\left(d^{\mu}-\sum\limits_{\nu\neq\mu}^{K}d^{\nu}\right) \left(2\rho\sum\limits_{\nu=1}^{K}d^{2\nu}\right)^{-1/2}\right]=\mathrm{erf} \left[\sqrt{\frac{1+d}{2\rho(1-d)}}\frac{2d^{\mu-1}-1-2d^{\mu}+d^{K}}{\sqrt{1- d^{2K}}}\right]>d^{\mu-1}\mathrm{erf}\left[\Theta\right] \tag{3.19}\]

and we can write the next

**Proposition 2**.: _In the noiseless limit \(\beta\to\infty\), the critical threshold for learning \(M_{\otimes}\) (in the number of required examples) depends on the data-set noise \(r\), the dilution \(d\), the amount of archetypes to handle \(K\) (and of course on the amplitude of the chosen confidence interval \(\Theta\)) and reads as_

\[M_{\otimes}^{\mu}(\Theta,r,d,K)>2\left(\mathrm{erf}^{-1}\left[d^{\mu-1} \mathrm{erf}\left[\Theta\right]\right]\right)^{2}\left(\frac{1-r^{2}}{r^{2}} \right)\frac{(1-d)(1-d^{2K})}{(1+d)(2d^{\mu-1}-1-2d^{\mu}+d^{K})^{2}} \tag{3.20}\]

_and in the plots (see Figure 5) we use \(\Theta=1/\sqrt{2}\) as this choice corresponds to the fairly standard condition \(\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}[\xi_{i}^{1}h_{i}^{(1)}(\boldsymbol{ \xi}^{1})]>\sqrt{\mathrm{Var}[\xi_{i}^{1}h_{i}^{(1)}(\boldsymbol{\xi}^{1})]}\) when \(\mu=1\)._

To quantify these thresholds for learning, in Fig. 5 we report the required number of examples to learn the first archetype (out of \(K=2,3,4,50\) as shown in the various panels) as a function of the dilution of the network.
#### 3.1.2 Ergodicity breaking: the critical phase transition

The main interest in the statistical mechanical approach to neural networks lies in inspecting their emerging capabilities, that typically appear once ergodicity gets broken: as a consequence, finding the boundaries of validity of ergodicity is a classical starting point to deepen these aspects.

To this task, hereafter we provide a systematic fluctuation analysis of the order parameters: the underlying idea is to check when, starting from the high noise limit (\(\beta\to 0\), where everything is uncorrelated and simple Probability arguments apply straightforwardly), these fluctuations diverge, as that defines the onset of ergodicity breaking as stated in the next

**Theorem 2**.: _The ergodic region, in the space of the control parameters \((\beta,d,\rho)\) is confined to the half-plane defined by the critical line_

\[\beta_{c}=\frac{1}{1-d}, \tag{3.21}\]

_whatever the entropy of the data-set \(\rho\)._

Proof.: The idea of the proof is the same we used so far, namely Guerra interpolation but on the rescaled fluctuations rather that directly on the statistical pressure.

Figure 5: We plot the logarithm of the critical number of example (required to raise the first magnetization) \(M_{\otimes}^{1}\) at different loads \(K=2,3,4,5\) and as a function of the dilution of the networks, for different noise values of the data-set (as shown in the legend). Note the divergent behavior of \(M_{\otimes}^{1}\) when approaching the critical dilution level \(d_{c}(K)=d_{1}\), as predicted by the parallel Hebbian storage limit [2; 4]: this is the crossover between the two multi-tasking regimes, hierarchical vs parallel, hence, solely at the value of dilution \(d_{1}\), there is no sharp behavior to infer and, correctly, the network can not accomplish learning. This shines by looking at (3.20) where the critical amount of examples to correctly infer the archetype is reported: its denominator reduces to \(1-2d+d^{K}\) and, for \(\mu=1\), it becomes zero when \(d\to d_{1}\).
The rescaled fluctuations \(\tilde{n}_{\nu}^{2}\) of the magnetizations are defined as

\[\tilde{n}_{\nu}=\sqrt{N}(n_{\nu}-\bar{n}_{\nu}). \tag{3.22}\]

We remind that the interpolating framework we are using, for \(t\in(0,1)\), is defined via

\[Z(t)=\sum_{\{\sigma\}}\exp\left[\frac{\beta}{2}tN(1+\rho)\sum_{\mu=1}^{K}n_{ \mu}^{2}+N(1-t)\beta(1+\rho)N_{\mu}n_{\mu}\right], \tag{3.23}\]

and it is a trivial exercise to show that, for any smooth function \(F(\sigma)\) the following relation holds:

\[\frac{d\langle F\rangle}{dt}=\frac{\beta}{2}(1+\rho)\left(\langle F\sum_{\nu} \tilde{n}_{\nu}^{2}\rangle-\langle F\rangle\langle\sum_{\nu}\tilde{n}_{\nu}^{ 2}\rangle\right), \tag{3.24}\]

such that by choosing \(F=\tilde{n}_{\mu}^{2}\) we can write

\[\begin{split}\frac{d\langle\tilde{n}_{\mu}^{2}\rangle}{dt}& =\frac{\beta}{2}(1+\rho)\left(\langle\tilde{n}_{\mu}^{2}\sum_{\nu} \tilde{n}_{\nu}^{2}\rangle-\langle\tilde{n}_{\mu}^{2}\rangle\langle\sum_{\nu }\tilde{n}_{\nu}^{2}\rangle\right)\\ &=\frac{\beta}{2}(1+\rho)\left(\langle\tilde{n}_{\mu}^{4}\rangle+ \langle\bar{n}_{\mu}^{2}\sum_{\nu\neq\mu}\tilde{n}_{\nu}^{2}\rangle-\langle \tilde{n}_{\mu}^{2}\rangle^{2}-\langle\tilde{n}_{\mu}^{2}\rangle\langle\sum_{ \nu\neq\mu}\tilde{n}_{\nu}^{2}\rangle\right)\\ &=\beta(1+\rho)\langle\tilde{n}_{\mu}^{2}\rangle^{2}\end{split} \tag{3.25}\]

thus we have

\[\langle\tilde{n}_{\mu}^{2}\rangle_{t}=\frac{\langle\tilde{n}_{\mu}^{2}\rangle _{t=0}}{1-t\beta(1+\rho)\langle\tilde{n}_{\mu}^{2}\rangle_{t=0}} \tag{3.26}\]

where the Cauchy condition \(\langle\tilde{n}_{\mu}^{2}\rangle_{t=0}\) reads

\[\begin{split}\langle\tilde{n}_{\mu}^{2}\rangle_{t=0}& =\,N\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\frac{\sum_{\{\sigma\}} \left(\frac{1}{N^{2}(1+\rho)^{2}}\sum_{i,j}\tilde{\eta}_{i}^{\mu}\tilde{\eta}_ {j}^{\mu}\sigma_{i}\sigma_{j}+\tilde{n}_{\mu}^{2}-2\frac{1}{N(1+\rho)}\sum_{i }\hat{\eta}_{i}^{\mu}\sigma_{i}\tilde{n}_{\mu}\right)\exp\left[\beta\sum_{\nu }\tilde{n}_{\nu}\sum_{i}\tilde{\eta}_{i}^{\nu}\sigma_{i}\right]}{\sum_{\{ \sigma\}}\exp\left[\beta\sum_{\nu}\tilde{n}_{\nu}\sum_{i}\tilde{\eta}_{i}^{ \nu}\sigma_{i}\right]}\\ &=\,\frac{1-d}{(1+\rho)}-N_{\mu}^{2}.\end{split} \tag{3.27}\]

Evaluating \(\langle\tilde{n}_{\mu}^{2}\rangle_{t}\) for \(t=1\), that is when the interpolation scheme collapses to the Statistical Mechanics, we finally get

\[\langle\tilde{n}_{\mu}^{2}\rangle_{t=1}=\frac{1-d-(1+\rho)N_{\mu}^{2}}{[1- \beta\left(1-d-(1+\rho)N_{\mu}^{2}\right)]} \tag{3.28}\]

namely the rescaled fluctuations are described by a meromorphic function whose pole is

\[\beta=\frac{1}{\left(1-d-(1+\rho)N_{\mu}^{2}\right)}\xrightarrow{N_{\mu}=0} \beta_{{}_{C}}=\frac{1}{1-d}, \tag{3.29}\]

that is the critical line reported in the statement of the theorem.
### Stability analysis via standard Hessian: the phase diagram

The set of solutions for the self-consistent equations for the order parameters (3.10) describes as candidate solutions a plethora of states whose stability must be investigated to understand which solution is preferred as the control parameters are made to vary: this procedure results in picturing the phase diagrams of the network, namely plots in the space of the control parameters where different regions pertain to different macroscopic computational capabilities.

Remembering that \(A_{K,\beta,d,M,r}(\bar{\mathbf{n}})=-\beta f_{K,\beta,d,M,r}(\bar{\mathbf{n}})\) (where \(f_{K,\beta,d,M,r}(\bar{\mathbf{n}})\) is the free energy of the model), in order to evaluate the stability of these solutions, we need to check the sign of the second derivatives of the free energy. More precisely, we need to build up the Hessian, a matrix \(\mathbf{A}\) whose elements are

\[\frac{\partial^{2}f(\bar{\mathbf{n}})}{\partial n^{\mu}\partial n^{\nu}}=A^{\mu\nu}\,. \tag{3.30}\]

Then, we evaluate and diagonalize \(\mathbf{A}\) at a point \(\tilde{\mathbf{n}}\), representing a particular solution of the self-consistency equation (3.10): the numerical results are reported in the phase diagrams provided in Fig.6.

We find straightforwardly

\[A^{\mu\nu}=(1+\rho)\left[[1-\beta(1-d)]+\rho\beta\mathbb{E}\left\{\mathcal{T} _{K\beta,\rho}^{2}(\bar{\mathbf{n}},z)(\xi^{\mu})^{2}\right\}\right]\delta^{\mu \nu}+Q^{\mu\nu} \tag{3.31}\]

where we set \(\mathcal{T}_{K\beta,\rho}(\bar{\mathbf{n}},z)=\tanh\left(\beta\sum_{\lambda=1}^{K} \bar{n}_{\lambda}\xi^{\lambda}+z\beta\sqrt{\rho\sum_{\lambda=1}^{K}(\bar{n}_{ \lambda}\xi^{\lambda})^{2}}\right)\) and

\[\begin{array}{rl}Q^{\mu\nu}&=\,\beta\mathbb{E}\left\{\left[\mathcal{T}_{K \beta,\rho}^{2}(\bar{\mathbf{n}},z)\right]\xi^{\mu}\xi^{\nu}\right\}(1-\delta^{\mu \nu})+2\rho\beta^{2}\mathbb{E}\left\{\left[\mathcal{T}_{K\beta,\rho}(\bar{\mathbf{ n}},z)\right]\left[1-\mathcal{T}_{K\beta,\rho}^{2}(\bar{\mathbf{n}},z)\right]\left[ \bar{n}_{\nu}\xi^{\nu}+\bar{n}_{\mu}\xi^{\mu}\right]\xi^{\mu}\xi^{\nu}\right\} \\ &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \
where \({\cal T}=\tanh\left[\beta\bar{n}\xi^{\mu}(1+z\sqrt{\rho})\right]\). It is easy to check that \(A\) becomes diagonal, with

\[A^{\mu\mu} = (1+\rho)\Big{[}1-\beta(1-d)+\beta(1-d)\,\mathbb{E}\left\{{\cal T}^{ 2}\right\}\Big{]}\] \[+4\beta^{2}\rho\bar{n}(1-d)\mathbb{E}\left\{{\cal T}\Big{[}1-{ \cal T}^{2}\Big{]}\right\}+2\beta^{3}\rho^{2}\bar{n}^{2}(1-d)\mathbb{E}\left\{ \Big{[}1-3{\cal T}^{2}\Big{]}\Big{[}1-{\cal T}^{2}\Big{]}\right\}\,,\] \[A^{\nu\nu\neq\mu} = (1+\rho)\Big{[}1-\beta(1-d)+\beta(1-d)^{2}\,\mathbb{E}\left\{{ \cal T}^{2}\right\}\Big{]}\,.\]

Notice that these eigenvalues do not depend on \(K\) since \({\cal T}\) does not depend on \(K\). Requiring the positivity for all the eigenvalues, we get the region in the plane \((d,\beta^{-1})\), where the pure state is stable: this correspond the blue region in the phase diagrams reported in Fig. 6.

We stress that these pure state solutions, namely the standard Hopfield-type ones, in the ground state (\(\beta^{-1}\to 0\)) are never stable whenever \(d\neq 0\) as the multi-tasking setting prevails. Solely at positive

Figure 6: Phase diagram in the dilution-noise (\(d-\beta^{-1}\)) plane for different values of \(K\) and \(\rho\). We highlight that different regions –marked with different colors– represent different operational behavior of the network: in yellow the ergodic solution, in light-blue the pure state solution (that is, solely one magnetization different from zero), in white the hierarchical regime (that is, several magnetizations differ from zero and they all assume different values) and in light-green the parallel regime (several magnetization differ from zero but their amplitude is the same for all).
values of \(\beta\), this single-pattern retrieval state is possible as the role of the noise is to destabilize the weakest magnetization of the hierarchical displacement, _vide infra_).

#### 3.2.3 Parallel state: \(\bar{\boldsymbol{n}}=\bar{n}_{d,\rho,\beta}(1,\ldots,1)\)

In this case the structure of the solution has the form of a symmetric mixture state corresponding to the unique self consistency equation for all \(\mu=1,\ldots K\), namely

\[\bar{n}\,=\,\frac{\mathbb{E}_{\xi,Z}\left\{\tanh\left[g(\beta,\boldsymbol{\xi},Z,\bar{n})\right]\xi^{\mu}\right\}}{(1+\rho)}+\beta\frac{\rho\,\bar{n}}{(1+ \rho)}\mathbb{E}_{\xi,Z}\left\{\left[1-\tanh^{2}\left(g(\beta,\boldsymbol{\xi} Z,\bar{n})\right)\right](\xi^{\mu})^{2}\right\}, \tag{3.37}\]

where

\[g(\beta,\boldsymbol{\xi},Z,\bar{n})=\beta\bar{n}\left[\sum_{\lambda=1}^{K}\xi^ {\lambda}+\beta\,Z\sqrt{\rho\,\sum_{\lambda=1}^{K}\left(\xi^{\lambda}\right)^ {2}}\right]. \tag{3.38}\]

In this case, the diagonal terms of \(\boldsymbol{A}\) are

\[\begin{array}{rcl}a=A^{\mu\mu}\,=&\Big{[}1-\beta(1-d)+\beta\,\mathbb{E}\left \{\left[\mathcal{T}^{2}\right](\xi^{\mu})^{2}\right\}\Big{]}(1+\rho)\\ \\ &&+4\beta^{2}\rho\bar{n}\mathbb{E}\left\{\mathcal{T}\Big{[}1-\mathcal{T}^{2} \Big{]}\xi^{\mu}\right\}+2\beta^{3}\rho^{2}\bar{n}^{2}\mathbb{E}\left\{\Big{[} 1-3\mathcal{T}^{2}\Big{]}\Big{[}1-\mathcal{T}^{2}\Big{]}(\xi^{\mu})^{2} \right\}\,,\end{array} \tag{3.39}\]

instead the off-diagonal ones are

\[b=A^{\mu\nu\neq\mu}\,=\,\beta\mathbb{E}\left\{\left[\mathcal{T}^{2}\right] \xi^{\mu}\xi^{\nu}\right\}+2\rho^{2}\beta^{3}\bar{n}^{2}\mathbb{E}\left\{\left[ 1-3\mathcal{T}^{2}\right]\left[1-\mathcal{T}^{2}\right](\xi^{\mu}\xi^{\nu})^{ 2}\right\}\,. \tag{3.40}\]

In general the following relationship holds

\[\boldsymbol{A}=\begin{pmatrix}a&b&\cdots&b&b\\ b&a&\cdots&b&b\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ b&b&\cdots&a&b\\ b&b&\cdots&b&a\end{pmatrix} \tag{3.41}\]

This matrix has always only two kinds of eigenvalues, namely \(a-b\) and \(a+(K-1)b\), thus, for the stability of the parallel state, after computing (3.39) and (3.40), we have only to check for which point in the \((d-\beta^{-1})\) plane both \(a-b\) and \(a+(K-1)b\) are positive. The region, in the phase diagrams of Fig. 6, where the parallel regime is stable is depicted in green.

2.4 Hierarchical state: \(\bar{\boldsymbol{n}}=\bar{n}_{d,\rho,\beta}((1-d),d(1-d),d^{2}(1-d),...)\)

In this case the structure of the solution has the hierarchical form \(\bar{\boldsymbol{n}}=\bar{n}_{d,\rho,\beta}((1-d),d(1-d),d^{2}(1-d),...)\) and the region left untreated so far in the phase diagram, namely the white region in the plots of Fig. 6, is the room left to such hierarchical regime.

### From the Cost function to the Loss function

We finally comment on the persistence -in the present approach- of quantifiers related to the evaluation of the pattern recognition capabilities of neural networks, i.e. the Mattis magnetization, also as quantifiers of a good learning process. The standard Cost functions used in Statistical Mechanics of neural networks (e.g., the Hamiltonians) can be related one-to-one to standard Loss functions used
in Machine Learning (i.e. the squared sum error functions), namely, once introduced the two Loss functions \(L_{\mu}^{+}:=(1/2N)||\xi^{\mu}-\sigma||^{2}=1-m_{\mu}\) and \(L_{\mu}^{-}=(1/2N)||\xi^{\mu}+\sigma||^{2}=1+m_{\mu}\)5, it is immediate to show that

Footnote 5: Note that in the last passage we naturally highlighted the presence of the Mattis magnetization in these Loss functions.

\[H(\mathbf{\sigma}|\mathbf{\xi})=\frac{-1}{2N}\sum_{i,j}^{N,N}\sum_{\mu}^{K}\xi_{i}^{\mu }\xi_{j}^{\mu}\sigma_{i}\sigma_{j}\equiv-N\sum_{\mu}^{K}\left(1-L_{\mu}^{+} \cdot L_{\mu}^{-}\right),\]

thus minimizing the former implies minimizing the latter such that, if we are extremizing w.r.t. the neurons we are performing machine retrieval (i.e. pattern recognition), while if we extremize w.r.t. the weights we perform machine leaning: indeed, at least in this setting, learning and retrieval are two faces of the same coin (clearly the task here, from a machine learning perspective, is rather simple as the network is just asked to correctly classify the examples and possibly generalize).

In Fig. 7 we inspect what happens to these Loss functions -pertaining to the various archetypes-as the Cost function gets minimized: we see that, at difference with the standard Hopfield model (where solely one Loss function per time diminishes its value), in this parallel learning setting several Loss functions (related to different archetypes) get simultaneously lowered, as expected by a parallel learning machine.

## 4 Conclusions

Since the AGS milestones on Hebbian learning dated 1985 [9; 10], namely the first comprehensive statistical mechanical theory of the Hopfield model for pattern recognition and associative memory, attractor neural networks have experienced an unprecedented growth and the bulk of techniques developed for spin glasses in these four decades (e.g. replica trick, cavity method, message passage,

Figure 7: Left: Parallel minimization of several (mean square-error) Loss functions \(L_{\pm}=||\xi^{\mu}\pm\sigma||^{2}\) (each pertaining to a different archetype) as the noise in the data-set \(r\) is varied. Here: \(M=25\), \(N=10000\). The horizontal gray dashed lines are the saturation level of the Loss functions, namely \(1-\frac{d}{2}-(1-d)d^{\mu-1}\). We get \(r_{\otimes}\) (the vertical black line) by the inversion of (3.20). Right: Parallel minimization of several (mean square-error) Loss function \(L_{\pm}=||\xi^{\mu}\pm\sigma||^{2}\) (each pertaining to a different archetype) as the data-set size \(M\) is varied: as M grows the simultaneous minimization of more than one Loss functions takes place, at difference with learning via standard Hebbian mechanisms where one Loss function -dedicate to a single archetype- is minimized per time. Orange and blue lines pertain to Loss functions of other patterns that, at these levels of dilution and noise, can not be minimized at once with the previous ones.
interpolation) acts now as a prosperous cornucopia for explaining the emergent information processing capabilities that these networks show as their control parameters are made to vary.

In these regards, it is important to stress how nowadays it is mandatory to optimize AI protocols (as machine learning for complex structured data-sets is still prohibitively expensive in terms of energy consumption [38]) and, en route toward a Sustainable AI (SAI), statistical mechanics may still pave a main theoretical strand: in particular, we highlight how the knowledge of the _phase diagrams_ related to a given neural architecture (that is the ultimate output product of the statistical mechanical approach) allows to set "a-priori" the machine in the optimal working regime for a given task thus unveiling a pivotal role of such a methodology even for a conscious usage of AI (e.g. it is useless to force a standard Hopfield model beyond its critical storage capacity)6.

Footnote 6: Further, still searching for optimization of resources, such a theoretical approach can also be used to inspect key computational shortcuts (as, e.g. early stopping criteria faced in Appendix C.1 or the role of the size of the mini-batch used for training [40] or of the flat minima emerging after training [15].

Focusing on Hebbian learning, however, while the original AGS theory remains a solid pillar and a paradigmatic reference in the field, several extensions are required to keep it up to date to deal with modern challenges: the first generalization we need is to move from a setting where the machine stores already defined patterns (as in the standard Hopfield model) toward a more realistic learning procedure where these patterns are unknown and have to be inferred from examples: the Hebbian storage rule of AGS theory quite naturally generalizes toward both supervised and an unsupervised learning prescriptions [5, 13]. This enlarges the space of the control parameters from \(\alpha,\beta\) (or \(K\), \(N\), \(\beta\)) of the standard Hopfield model toward \(\alpha,\beta,\rho\) (or \(K\), \(N\), \(\beta\), \(M\), \(r\)) as we now deal also with a data-set where we have \(M\) examples of mean quality r for each pattern (archetype) or, equivalently, we speak of a data-set produced at given entropy \(\rho\).

Once this is accomplished, the second strong limitation of the original AGS theory that must be relaxed is that patterns share the same length and, in particular, this equals the size of the network (namely in the standard Hopfield model there are N neurons to handle patterns, whose length is exactly N for all of them): a more general scenario is provided by dealing with patterns that contain different amounts of information, that is patterns are diluted. Retrieval capabilities of the Hebbian setting at work with diluted patterns have been extensively investigated in the last decade [2, 3, 23, 24, 28, 35, 42, 44, 47] and it has been understood how, dealing with patterns containing sparse entries, the network automatically is able to handle several of them in parallel (a key property of neural networks that is not captured by standard AGS theory). However the study of the parallel learning of diluted patterns was not addressed in the Literature and in this paper we face this problem, confining this first study to the low storage regime, that is when the number of patterns scales at most logarithmically with the size of the network. Note that this further enlarges the space of the control parameters by introducing the dilution \(d\): we have several control parameters because the network information processing capabilities are enriched w.r.t. the bare Hopfield reference7.

Footnote 7: However we clarify how it could be inappropriate to speak about structural differences among the standard Hopfield model and the present multitasking counterpart: ultimately these huge behavioral differences are just consequences of the different nature of the the data-sets provided to the same Hebbian network during training.

We have shown here that if we supply to the network a data-set equipped with dilution, namely a sparse data-set whose patterns contain -on average- a fraction \(d\) of blank entries (whose value is 0) and, thus, a fraction \((1-d)\) of informative entries (whole values can be \(\pm 1\)), then the network spontaneously undergoes parallel learning and behaves as a multitasking associative memory able to learn, store and retrieve multiple patterns in parallel. Further, focusing on neurons, the Hamiltonian of the model plays as the Cost function for neural dynamics, however, moving the att
have shown how the latter is one-to-one related to the standard (mean square error) Loss function in Machine Learning and this resulted crucial to prove that, by experiencing a diluted data-set, the network lowers in parallel several Loss functions (one for each pattern that is learning from the experienced examples).

For mild values of dilution, the most favourite displacement of the Mattis magnetizations is a _hierarchical ordering_, namely the intensities of these signals scale as power laws w.r.t. their information content \(m_{K}\sim d^{K}\cdot(1-d)\), while at high values of dilution a _parallel ordering_, where all these amplitudes collapse to the same value, prevails: the phase diagrams of these networks properly capture these different working regions.

Remarkably, confined to the low storage regime (where glassy phenomena can be neglected), the presence (or its lacking) of a teacher does not alter the above scenario and the threshold for a secure learning, namely the minimal required amount of examples (given the constraints, that is the noise in the data-set r, the amount of different archetypes \(K\) to cope with, etc.) \(M_{\otimes}\) that guarantees that the network is able to infer the archetype and thus generalize, is the same for supervised and unsupervised protocols and its value has been explicitly calculated: this is another key point toward sustainable AI. Clearly there is still a long way to go before a full statistical mechanical theory of extensive parallel processing will be ready, yet this paper acts as a first step in this direction and we plan to report further steps in a near future.
A more general sampling scenario

The way in which we add noise over the archetypes to generate the data-set in the main text (see eq. (9)) is a rather peculiar one as, in each example, it preserves the number but also the positions of lacunae already present in the related archetype. This implies that the noise can not affect the amplitudes of the original signal, i.e. \(\sum_{i}(\eta^{\mu,a}_{i})^{2}=\sum_{i}(\xi^{\mu}_{i})^{2}\) holds for any \(a\) and \(\mu\), while we do expect that with more general kinds of noise the dilution, this property is not preserved sharply.

Here we consider the case where the number of blank entries present in \(\mathbf{\xi}^{\mu}\) is preserved on average in the related sample \(\{\eta^{\mu,a}\}_{a=1,\dots,M}\) but lacunae can move along the examples: this more realistic kind of noise gives rise to cumbersome calculations (still analytically treatable) but should not affect heavily the capabilities of learning, storing and retrieving of these networks (as we now prove).

Specifically here we define the new kind of examples \(\tilde{\eta}^{\mu,a}_{i}\) (that we can identify from the previous ones \(\eta^{\mu,a}_{i}\) by labeling them with a tilde) in the following way

**Definition 10**.: _Given \(K\) random patterns \(\mathbf{\xi}^{\mu}\) (\(\mu=1,...,K\)), each of length \(N\), whose entries are i.i.d. from_

\[\mathbb{P}(\xi^{\mu}_{i})=\frac{(1-d)}{2}\delta_{\xi^{\mu}_{i},-1}+\frac{(1-d )}{2}\delta_{\xi^{\mu}_{i},+1}+d\delta_{\xi^{\mu}_{i},0}, \tag{10}\]

_we use these archetypes to generate \(M\times K\) different examples \(\{\tilde{\eta}^{\mu,a}_{i}\}^{a=1,\dots,M}\) whose entries are depicted following_

\[\begin{split}&\mathbb{P}(\tilde{\eta}^{\mu,a}_{i}|\xi^{\mu}_{i}= \pm 1)=A_{\pm}(r,s)\delta_{\tilde{\eta}^{\mu,a}_{i},\xi^{\mu}_{i}}+B_{\pm}(r,s) \delta_{\tilde{\eta}^{\mu,a}_{i},-\xi^{\mu}_{i}}+C_{\pm}(r,s)\delta_{\tilde{ \eta}^{\mu,a}_{i},0}\\ &\mathbb{P}(\tilde{\eta}^{\mu,a}_{i}|\xi^{\mu}_{i}=0)=A_{0}(r,s) \delta_{\tilde{\eta}^{\mu,a}_{i},\xi^{\mu}_{i}}+B_{0}(r,s)\delta_{\tilde{\eta }^{\mu,a}_{i},+1}+C_{0}(r,s)\delta_{\tilde{\eta}^{\mu,a}_{i},-1}\end{split} \tag{11}\]

_for \(i=1,\dots,N\) and \(\mu=1,\dots,K\), where we pose_

\[\begin{split}& A_{\pm}(r,s)=\frac{1+r}{2}\left[1-\frac{d}{1-d}(1 -s)\right]+\frac{d(1-s)(1-r)}{4(1-d)}\,,\qquad A_{0}(r,s)=\frac{1+s}{2}\,,\\ & B_{\pm}(r,s)=\frac{1-r}{2}\left[1-\frac{d}{1-d}(1-s)\right]+ \frac{d(1-s)(1+r)}{4(1-d)}\,,\qquad B_{0}(r,s)=\frac{1-s}{4}\,,\\ & C_{\pm}(r,s)=\frac{d}{2(1-d)}(1-s)\,,\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad C_{0}(r,s)=\frac{1-s}{4}\,,\end{split} \tag{12}\]

_with \(r,s\in[0;1]\) (whose meaning we specify soon,_ vide infra_)._

Equation (11) codes for the new noise, the values of the coefficients presented in (12) have been chosen in order that all the examples contain on average the same fraction \(d\) of null entries as the original archetypes. To see this it is enough to check that the following relation holds for each \(a=1,\dots,M\), \(i=1,\dots,N\) and \(\mu=1,\dots,K\)

\[\mathbb{P}(\tilde{\eta}^{\mu,a}_{i}=0)=\sum_{x\in\{-1,0,1\}}\mathbb{P}(\tilde{ \eta}^{\mu,a}_{i}=0|\xi^{\mu}_{i}=x)\mathbb{P}(\xi^{\mu}_{i}=x)=C_{\pm}(r,s)(1 -d)+A_{0}(r,s)d=d\,. \tag{13}\]

Defined the data-set, the cost function follows straightforwardly in Hebbian settings as

**Definition 11**.: _Once introduced \(N\) Ising neurons \(\sigma_{i}=\pm 1\) (\(i=1,...,N\)) and the data-set considered in the definition above, the Cost function of the multitasking Hebbian network equipped with not-preserving-dilution noise reads as_

\[\mathcal{H}^{(sup)}_{N,K,M,r,s,d}(\mathbf{\sigma}|\tilde{\mathbf{\eta}})=-\frac{1}{N} \frac{1}{(1-d)(1+\tilde{\rho})}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\left(\frac{1 }{\tilde{r}M}\sum_{a=1}^{M}\tilde{\eta}^{\mu,a}_{i}\right)\left(\frac{1}{ \tilde{r}M}\sum_{b=1}^{M}\tilde{\eta}^{\mu,b}_{j}\right)\sigma_{i}\sigma_{j}, \tag{14}\]
_where_

\[\tilde{r}=\frac{r}{(1-d)}\left[1-\frac{d}{2}(5-3s)\right]\] (A.6)

_and \(\tilde{\rho}\) is the generalization of the data-set entropy, defined as:_

\[\tilde{\rho}=\frac{1-\tilde{r}^{2}}{M\tilde{r}^{2}}\,.\] (A.7)

**Definition 12**.: _The suitably re-normalized example's magnetizations \(n_{\mu}\) read as_

\[n_{\mu}:=\frac{1}{(1+\tilde{\rho})}\frac{1}{N}\sum_{i=1}^{N}\left(\frac{1}{ \tilde{r}M}\sum_{a=1}^{M}\tilde{\eta}_{i}^{\mu,a}\right)\sigma_{i}\,.\] (A.8)

En route toward the statistical pressure, still preserving Guerra's interpolation as the underlying technique, we give the next

**Definition 13**.: _Once introduced the noise \(\beta\in\mathbb{R}^{+}\), an interpolating parameter \(t\in(0,1)\), the \(K+1\) auxiliary fields \(J\) and \(\psi_{\mu}\) (\(\mu\in(1,...,K)\)), the interpolating partition function related to the model defined by the Cost function (A.5) reads as_

\[\mathcal{Z}^{(sup)}_{\beta,N,K,M,r,s,d}(\boldsymbol{\xi},\boldsymbol{\tilde{ \eta}};J,t)=\sum_{\{\boldsymbol{\sigma}\}}\exp\Bigg{[}\ J\sum_{\mu,i=1}^{K,N} \xi_{i}^{\mu}\sigma_{i}+t\beta N\frac{(1+\tilde{\rho})}{2(1-d)}\sum_{\mu=1}^{ K}n_{\mu}^{2}(\boldsymbol{\sigma})+(1-t)N\sum_{\mu=1}^{K}\psi_{\mu}\,n_{\mu}( \boldsymbol{\sigma})\Bigg{]}.\] (A.9)

_and the interpolating statistical pressure \(\mathcal{A}_{\beta,K,M,r,s,d}=\lim_{N\to\infty}A_{\beta,N,K,M,r,s,d}\) induced by the partition function (A.9) reads as_

\[A_{\beta,N,K,M,r,s,d}(J,t)=\frac{1}{N}\mathbb{E}\Big{[}\ln\mathcal{Z}^{(sup) }_{\beta,N,K,M,r,s,d}(\boldsymbol{\xi},\boldsymbol{\tilde{\eta}};J,t)\Big{]}\] (A.10)

_where \(\mathbb{E}=\mathbb{E}_{\xi}\mathbb{E}_{(\tilde{\eta}|\xi)}\)._

**Remark 3**.: _Of course, as in the model studied in the main text still with Guerra's interpolation technique, we aim to find an explicit expression (in terms of the control and order parameters of the theory) of the interpolating statistical pressure evaluated at \(t=1\) and \(J=0\)._

We thus perform the computations following the same steps of the previous investigation: the \(t\) derivative of interpolating pressure is given by

\[\frac{d\mathcal{A}_{\beta,K,M,r,s,d}(J,t)}{dt}=\,\frac{\beta}{2(1-d)}(1+\tilde {\rho})\sum_{\mu=1}^{K}\langle n_{\mu}^{2}\rangle_{t}-\sum_{\mu=1}^{K}\psi_{ \mu}\langle n_{\mu}\rangle_{t}.\] (A.11)

fixing

\[\psi_{\mu}=\frac{\beta}{1-d}(1+\tilde{\rho})\bar{n}_{\mu}\] (A.12)

and computing the one-body term

\[\begin{split}\mathcal{A}_{\beta,K,M,r,s,d}(J,t=0)&= \mathbb{E}\ln\,\left[2\cosh\left(\sum_{\mu=1}^{K}\psi_{\mu}\frac{1}{(1+\tilde {\rho})}\frac{1}{\tilde{r}M}\sum_{a=1}^{M}\tilde{\eta}_{i}^{\mu,a}+J\sum_{\mu= 1}^{K}\xi^{\mu}\right)\right]\\ &=\mathbb{E}\ln\,\left\{2\cosh\left[\frac{\beta}{1-d}\sum_{\mu=1 }^{K}\bar{n}_{\mu}\left(\frac{1}{\tilde{r}M}\sum_{a=1}^{M}\tilde{\eta}_{i}^{ \mu,a}\right)+J\sum_{\mu=1}^{K}\xi^{\mu}\right]\right\}.\end{split}\] (A.13)

We get the final expression as \(N\to\infty\) such that we can state the next
**Theorem 3**.: _In the thermodynamic limit \((N\to\infty)\) and in the low load regime \((K/N\to 0)\), the quenched statistical pressure of the multitasking Hebbian network equipped with not-preserving-dilution noise, whatever the presence of a teacher, reads as_

\[\mathcal{A}_{\beta,K,M,r,s,d}(J)\,=\,\mathbb{E}\left\{\ln\left[2\cosh\left( \beta^{{}^{\prime}}\sum_{\mu=1}^{K}\bar{n}_{\mu}\tilde{\eta}^{\mu}+J\sum_{\mu= 1}^{K}\xi^{\mu}\right)\,\right]\right\}-\frac{\beta^{{}^{\prime}}}{2}(1+\tilde {\rho})\sum_{\mu=1}^{K}\bar{n}_{\mu}^{2}. \tag{110}\]

_where \(\beta^{{}^{\prime}}=\beta/(1-d)\), \(\mathbb{E}=\mathbb{E}_{\xi}\mathbb{E}_{(\tilde{\eta}|\xi)}\) and \(\tilde{\eta}^{\mu}=\frac{1}{\tilde{r}M}\sum_{a=1}^{M}\tilde{\eta}_{i}^{\mu,a}\)and the values \(\bar{n}_{\mu}\) must fulfill the following self-consistent equations_

\[\bar{n}_{\mu}=\frac{1}{(1+\tilde{\rho})}\mathbb{E}\left\{\left[\tanh\left( \beta^{{}^{\prime}}\sum_{\nu=1}^{K}\bar{n}_{\nu}\tilde{\eta}^{\nu}\right) \,\right]\tilde{\eta}^{\mu}\right\}\quad\text{for}\;\;\mu=1,\ldots,K\,, \tag{111}\]

_that extremize the statistical pressure \(\mathcal{A}_{\beta,K,M,r,s,d}(J=0)\) w.r.t. them._

Furthermore, the simplest path to obtain a self-consistent equation also for the Mattis magnetization \(m_{\mu}\) is by considering the auxiliary field \(J\) coupled to \(m_{\mu}\) namely \(\bar{m}_{\mu}=\nabla_{J}\mathcal{A}_{\beta,K,M,r,s,d}(J)|_{J=0}\) to get

\[\bar{m}_{\mu}=\mathbb{E}\left\{\tanh\left[\beta^{{}^{\prime}}\sum_{\nu=1}^{K }\bar{n}_{\nu}\tilde{\eta}^{\nu}\right]\xi^{\mu}\right\}\quad\text{for}\;\;\mu =1,\ldots,K\,. \tag{112}\]

We do not plot these new self-consistency equations as, in the large \(M\) limit, there are no differences w.r.t. those obtained in the main text.

## Appendix B On the data-set entropy \(\rho\)

In this appendix, focusing on a single generic bit, we deepen the relation between the conditional entropy \(H(\xi_{i}^{\mu}|\mathbf{\eta}_{i}^{\mu})\) of a given pixel \(i\) regarding archetype \(\mu\) and the information provided by the data-set regarding such a pixel, namely the block \(\left(\eta_{i}^{\mu,1},\eta_{i}^{\mu,2},\ldots,\eta_{i}^{\mu,M}\right)\) to justify why we called \(\rho\) the data-set entropy in the main text. As calculations are slightly different among the two analyzed models (the one preserving dilution position provided in the main text and the generalized one given in the previous appendix) we repeat them model by model for the sake of transparency.

### I: multitasking Hebbian network equipped with not-affecting-dilution noise

Let us focus on the \(\mu\)-th pattern and the \(i\)-th digit, whose related block is

\[\eta_{i}^{\mu}=\left(\eta_{i}^{\mu,1},\eta_{i}^{\mu,2},\ldots,\eta_{i}^{\mu,M }\right); \tag{113}\]

the error probability for any single entry is

\[\mathbb{P}(\xi_{i}^{\mu}\neq 0)\mathbb{P}(\eta_{i}^{\mu,a}\neq\xi_{i}^{\mu})=(1-d )(1-r)/2 \tag{114}\]

and, by applying the majority rule on the block, it is reduced to

\[\mathbb{P}(\xi_{i}^{\mu}\neq 0)\mathbb{P}\left(\text{sign}\Big{(}\sum_{a}\eta_ {i}^{\mu,a}\xi_{i}^{\mu}\Big{)}=-1\right)\underset{M\gg 1}{\approx}\frac{(1-d)}{2} \left[1-\text{erf}\left(\frac{1}{\sqrt{2\rho}}\right)\right]. \tag{115}\]
Thus

\[H_{d,r,M}(\boldsymbol{\xi}^{\mu}|\boldsymbol{\eta}^{\mu})=-\left[x(d,r,M)\log_{2}x( d,r,M)+y(d,r,M)\log_{2}y(d,r,M)\right]\] (B.4)

where

\[x(d,r,M)=\frac{(1-d)}{2}\left[1-\mathrm{erf}\left(\frac{1}{\sqrt{2\rho}}\right) \right]\;,\;\;y(d,r,M)=1-x(d,r,M)\,.\] (B.5)

### II: multitasking Hebbian network equipped with not-preserving-dilution noise

Let us focus on the \(\mu\)-th pattern and the \(i\)-th digit, whose related block is

\[\tilde{\eta}_{i}^{\mu}=\left(\tilde{\eta}_{i}^{\mu,1},\tilde{\eta}_{i}^{\mu,2 },\ldots,\tilde{\eta}_{i}^{\mu,M}\right);\] (B.6)

the error probability for any single entry is

\[\mathbb{P}(\xi_{i}^{\mu}\neq 0)\mathbb{P}(\tilde{\eta}_{i}^{\mu,a}\xi_{i}^{ \mu}\neq+1|\xi_{i}^{\mu}\neq 0)+\mathbb{P}(\xi_{i}^{\mu}=0)\mathbb{P}( \tilde{\eta}_{i}^{\mu,a}\neq 0|\xi_{i}^{\mu}=0)=d(1-s)\,.\] (B.7)

By applying the majority rule on the block, it is reduced to

\[\begin{split}&\mathbb{P}(\xi_{i}^{\mu}\neq 0)\left[1-\mathbb{P} \Big{(}\mathrm{sign}(\hat{\eta}_{i}^{\mu}\xi_{i}^{\mu})=+1\Big{|}\xi_{i}^{\mu }\neq 0\Big{)}\right]+\mathbb{P}(\xi_{i}^{\mu}=0)\mathbb{P}\Big{(}\mathrm{ sign}[|\hat{\eta}_{i}^{\mu}|]=+1\Big{|}\xi_{i}^{\mu}=0\Big{)}\\ &\underset{M\gg 1}{\approx}\frac{(1-d)}{2}\left\{1-\mathrm{erf} \left[\left(2\tilde{\rho}-\frac{d(1-s)}{(1-d)M\tilde{r}^{2}}\right)^{-1/2} \right]\right\}+\frac{d}{2}\;\left\{1-\mathrm{erf}\left[\left(\frac{1-s}{M \tilde{r}^{2}}\right)^{-1/2}\right]\right\}\,.\end{split}\] (B.8)

Figure 8: Comparison of the numerical solution of the self consistency equations related to the Mattis magnetization in the two models: upper panel is due to the first model (reported in the main text), lower panel reports on the second model (deepened here). Beyond a different transient at small \(M\) the two models behave essentially in the same way.
Thus

\[H_{d,r,s,M}(\xi_{i}^{\mu}|\tilde{\mathbf{\eta}}_{i}^{\mu})=-\left[x(d,r,s,M)\log_{2}x(d, r,s,M)+y(d,r,s,M)\log_{2}y(d,r,s,M)\right] \tag{114}\]

where

\[x(d,r,s,M)=\frac{(1-d)}{2}\left\{1-\text{erf}\left[\left(2\tilde{\rho}-\frac{d (1-s)}{(1-d)M\tilde{r}^{2}}\right)^{-1/2}\right]\right\}+\frac{d}{2}\ \left\{1-\text{erf}\left[\left(\frac{1-s}{M\tilde{r}^{2}}\right)^{-1/2} \right]\right\}\]

\[y(d,r,s,M)=1-x(d,\tilde{\rho}) \tag{115}\]

Whatever the model, the conditional entropies \(H_{d,r,M}(\xi_{i}^{\mu}|\mathbf{\eta}_{i}^{\mu})\) and \(H_{d,r,s,M}(\xi_{i}^{\mu}|\tilde{\mathbf{\eta}}_{i}^{\mu})\) are a monotonic increasing functions of \(\rho\) and \(\tilde{\rho}\) respectively, hence the reason to calling \(\rho\) and \(\tilde{\rho}\) the entropy of the data-set.

## Appendix C Stability analysis: an alternative approach

### Stability analysis via signal-to-noise technique

The standard signal-to-noise technique [8] is a powerful method to investigate the stability of a given neural configuration in the noiseless limit \(\beta\to\infty\): by requiring that each neuron is aligned to its field (the post-synaptic potential that it is experiencing, i.e. \(h_{i}\sigma_{i}\geq 0\ \ \forall i\in(1,...,N)\)) this analysis allows to correctly classify which solution (stemming from the self-consistent equations for the order parameters) is preferred as the control parameters are made to vary and thus it can play as an alternative route w.r.t the standard study of the Hessian of the statistical pressure reported in the main text (see Sec. 3.2).

In particular, recently, a revised version of the signal-to-noise technique has been developed [11; 12] and in this new formulation it is possible to obtain the self-consistency equations for the order parameters explicitly so we can compare directly outcomes from signal-to-noise to outcomes from statistical mechanics. By comparison of the two routes that lead to the same picture, that is statistical mechanics and revised signal-to-noise technique, we can better comprehend the working criteria of these neural networks.

We suppose that the network is in the hierarchical configuration prescribed by eq. (6), that we denote as \(\mathbf{\sigma}=\mathbf{\sigma}^{*}\), and we must evaluate the local field \(h_{i}(\mathbf{\sigma}^{*})\) acting on the generic neuron \(\sigma_{i}\) in this configuration to check that \(h_{i}(\mathbf{\sigma}^{*})\sigma_{i}^{*}>0\) is satisfied for any \(i=1,\ldots,N\): should this be the case, the configuration would be stable, vice versa unstable.

Focusing on the supervised setting with no loss of generality (as we already discussed that the teacher essentially plays no role in the low storage regime) and selecting (arbitrarily) the hierarchical ordering as a test case to be studied, we start by re-writing the Hamiltonian (11) as

\[-\mathcal{H}_{N,K,M,r}(\mathbf{\sigma}|\mathbf{\eta})=\sum_{i=1}^{N}h_{i}(\mathbf{\sigma}) \sigma_{i}\,, \tag{116}\]

where the local fields \(h_{i}\) appear explicitly and are given by

\[h_{i}(\mathbf{\sigma})= \frac{1}{2N\,r^{2}M^{2}(1-d)(1+\rho)}\sum_{\mu=1}^{K}\sum_{j\neq i }^{N}\sum_{a,b}^{M_{\mu},M_{\mu}}\eta_{i}^{\mu,a}\eta_{j}^{\mu,b}\sigma_{j}\,. \tag{117}\]
The updating rule for the neural dynamics reads as

\[\sigma_{i}^{(n+1)}=\sigma_{i}^{(n)}\text{sign}\left(\tanh\left[\beta\sigma_{i}^{(n) }h_{i}^{(n)}(\boldsymbol{\sigma}^{(n)})\right]+\Gamma_{i}\right)\ \ \text{with}\ \ \Gamma_{i}\sim\mathcal{U}[-1;+1]\,, \tag{104}\]

that, in the zero fast-noise limit \(\beta\to+\infty\), reduces to

\[\sigma_{i}^{(n+1)}=\sigma_{i}^{(n)}\text{sign}\left(\sigma_{i}^{(n)}h_{i}^{(n) }(\boldsymbol{\sigma}^{(n)})\right). \tag{105}\]

To inspect the stability of the hierarchical parallel configuration, we initialize the network in such configuration, i.e., \(\boldsymbol{\sigma}^{(1)}=\boldsymbol{\sigma}^{*}\), then, following Hinton's prescription [21, 48]8 the one-step n-iteration \(\boldsymbol{\sigma}^{(2)}\) leads to an expression of the magnetization that reads as

Footnote 8: The _early stopping prescription_ given by Hinton and coworkers became soon very popular, yet it has been criticized in some circumstances (in particular where glassy features are expected to be strong and may keep the network out of equilibrium for very long times, see e.g. [7, 18, 39]): we stress that, in the present section, we are assuming the network has already reached the equilibrium further, confined to the low storage inspection, spin glass bottlenecks in thermalization should not be prohibitive.

\[m_{\mu}^{(2)}=\frac{1}{N}\sum_{i=1}^{N}\xi_{i}^{\mu}\sigma_{i}^{(2)}=\frac{1}{ N}\sum_{i=1}^{N}\xi_{i}^{\mu}\left[\sigma_{i}^{*}\text{sign}\left(\sigma_{i}^{*}\,h_ {i}^{(1)}(\boldsymbol{\sigma}^{*})\right)\right]; \tag{106}\]

Next, using the explicit expression of the hierarchical parallel configuration (6), we get

\[m_{1}^{(2)}=\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i}^{1}\right)^{2 }\text{sign}\left(\sigma_{i}^{*}\,h_{i}^{(1)}(\boldsymbol{\sigma}^{*})\right); \tag{107}\] \[m_{\mu>1}^{(2)}=\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i}^{\mu} \right)^{2}\prod_{\rho=1}^{\mu-1}\delta\left(\xi_{i}^{\rho}\right)\text{sign} \left(\sigma_{i}^{*}\,h_{i}^{(1)}(\boldsymbol{\sigma}^{*})\right);\]

by applying the central limit theorem to estimate the sums appearing in the definition of \(h_{i}^{(1)}\) for \(i=1,\ldots,N\), we are able to split, mimicking the standard signal-to-noise technique, a signal contribution (\(\kappa_{1,\mu}^{(1)}\)) and a noise contribution (\(\kappa_{2,\mu}^{(1)}\)) as presented in the following

\[\sigma_{i}^{*}h_{i}^{(1)}(\boldsymbol{\sigma}^{*})\sim\kappa_{1,\mu}^{(1)}+z_ {i}\sqrt{\kappa_{2,\mu}^{(1)}}\ \ \ \ \text{with}\ \ z_{i}\sim\mathcal{N}(0,1) \tag{108}\]

where

\[\kappa_{1,\mu}^{(1)}\coloneqq\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\left[ \sigma_{i}^{*}\,h_{i}^{(1)}(\boldsymbol{\sigma}^{*})\right]\ \ \ \ \kappa_{2,\mu}^{(1)}\coloneqq\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\left[ \sigma_{i}^{*}\,h_{i}^{(1)}(\boldsymbol{\sigma}^{*})\right]^{2} \tag{109}\]

Thus, Eq. (106) becomes

\[m_{\mu}^{(2)} = \left[\frac{1}{N}\sum_{i=1}^{N}\left(\xi_{i}^{\mu}\right)^{2} \prod_{\rho=1}^{\mu-1}\delta\left(\xi_{i}^{\rho}\right)\text{sign}\left( \kappa_{1,\mu}^{(1)}+z_{i}\sqrt{\kappa_{2,\mu}^{(1)}}\right)\right]\ \ \ \ \text{ with}\ \ \mu=1,2,\ldots,K\,. \tag{110}\]

For large values of \(N\), the arithmetic mean coincides with the theoretical expectation, thus

\[\frac{1}{N}\sum_{i=1}^{N}g(\xi_{i},z_{i})\ \xrightarrow[N\to\infty]{} \mathbb{E}_{\xi,z}[g(\xi,z)]=\mathbb{E}_{\xi}\left[\int\frac{dz}{\sqrt{2\pi}} e^{-\frac{z^{2}}{2}}g(\xi,z)\right]\,. \tag{111}\]
therefore, we can rewrite Eq. (104) as

\[m_{\mu}^{(2)}\,=\,(1-d)d^{\mu-1}\mathrm{erf}\left[\frac{\kappa_{1,\mu}^{(1)}}{ \sqrt{2\Big{(}\kappa_{2,\mu}^{(1)}-\Big{(}\kappa_{1,\mu}^{(1)}\Big{)}^{2}\, \Big{)}}}\right]\quad\text{ with }\;\mu=1,2,\ldots,K\,. \tag{105}\]

While we carry out the computations of \(\kappa_{1,\mu}^{(1)}\) and \(\kappa_{2,\mu}^{(1)}\) in Appendix C.2, here we report only their values, which are

\[\kappa_{1,\mu}^{(1)}=\frac{1}{2}\frac{1}{(1+\rho)}d^{\mu-1}\,,\quad\kappa_{2, \mu}^{(1)}=\frac{1}{4}\frac{1}{1+\rho}d^{\mu-1}\left[\frac{1+d+d^{2}-d^{2K-2 \mu+2}}{1+d}\right]. \tag{106}\]

So we get

\[m_{\mu}^{(2)}(d,K,\rho)\,=\,(1-d)d^{\mu-1}\mathrm{erf}\left[\frac{1}{\sqrt{2}} \frac{\sqrt{(1+d)d^{\mu-1}}}{\sqrt{(1+\rho)(1+d+d^{2}-d^{2K-2\mu+2})-d^{\mu}-d ^{\mu-1}}}\right] \tag{107}\]

Figure 9: Signal to Noise numerical inspection of the Mattis magnetizations for a diluted network with \(r=0.1\) and \(K=3\) in the hierarchical regime (at levels of pattern’s dilution \(d<d_{c}\) as reported in their titles): we highlight the agreement, as the saturation level is reached, among Signal to Noise analysis (orange dots) and the value of the magnetization of the first pattern found by the mechanical statistical approach (reported as solid red line). The dashed lines represent the Hebbian storing prescriptions at which the values of the magnetizations converge. The vertical black line depicts the critical amount of examples \(M_{\otimes}\) that must be experienced by the network to properly depict the archetypes: note that this value is systematically above, in M, that the point where all the bifurcations happened, hence all the magnetizations stabilized on their hierarchical displacements.
As shown in Fig. 9, as a critical amount of perceived examples is collected, this expression is in very good agreement with the estimate stemming from the numerical solution of the self-consistent equations and indeed we can finally state the last

**Theorem 4**.: _In the zero fast-noise limit (\(\beta\to+\infty\)), if the neural configuration_

\[\tilde{\mathbf{\sigma}}=\tilde{\mathbf{\sigma}}(\mathbf{\xi}) \tag{102}\]

_is a fixed point of the dynamics described by the sequential spin upload rule_

\[\sigma_{i}^{(n+1)}=\sigma_{i}^{(n)}\mathrm{sign}\left[\beta\sigma_{i}^{(n)}h_{i }^{(n)}(\mathbf{\sigma}^{(n)})\right] \tag{103}\]

_where_

\[h_{i}^{(n)}(\mathbf{\sigma})= \frac{1}{N\,M^{2}(1+\rho)r^{2}}\sum\limits_{\mu=1}^{K}\sum\limits _{j\neq i}^{N}\sum\limits_{a,b}^{M_{\mu},M_{\mu}}\eta_{i}^{\mu,a}\eta_{j}^{\mu,b }\sigma_{j}^{(n)}\,, \tag{104}\]

_then the order parameters \(n_{\mu}(\mathbf{\sigma})=[NM(1+\rho)r]^{-1}\sum\limits_{i}^{N}\sum\limits_{a}^{M} \eta_{i}^{\mu,a}\sigma_{i}\) must satisfy the following self equations_

\[n_{\mu} = \frac{1}{(1+\rho)}\mathbb{E}_{\mathbf{\xi}}\mathbb{E}_{(\mathbf{\eta}|\bm {\xi})}\Bigg{\{}\hat{\eta}^{\mu}\tilde{\sigma}(\mathbf{\xi})\,\mathrm{sign}\left[ \sum\limits_{\nu=1}^{K}n_{\nu}\hat{\eta}^{\nu}\tilde{\mathbf{\sigma}}(\mathbf{\xi}) \right]\Bigg{\}}\,, \tag{105}\]

_where we set \(\hat{\eta}^{\mu}=(Mr)^{-1}\sum\limits_{a}^{M}\eta^{\mu,a}\)._

**Remark 4**.: _The empirical evidence that, via early stopping criterion, we still obtain the correct solution proves a posteriori the validity of Hinton's recipe in the present setting and it tacitly candidate statistical mechanics as a reference also to inspect computational shortcuts._

Proof.: The local fields \(h_{i}\) can be rewrite using the definition of \(n_{\mu}\) as

\[h_{i}^{(n)}(\mathbf{\sigma})=\sum\limits_{\mu=1}^{K}n_{\mu}^{(n)}(\mathbf{\sigma}) \hat{\eta}_{i}^{\mu}\,, \tag{106}\]

in this way the upload rule can be recast as

\[\sigma_{i}^{(n+1)}=\sigma_{i}^{(n)}\mathrm{sign}\left[\sum\limits_{\mu=1}^{K} n_{\mu}^{(n)}(\mathbf{\sigma})\hat{\eta}_{i}^{\mu}\sigma_{i}^{(n)}\right]. \tag{107}\]

Computing the value of the \(n_{\mu}\)-order parameters at the \((n+1)\) step of uploading process we get

\[n_{\mu}^{(n+1)}(\mathbf{\sigma}) = \frac{1}{N(1+\rho)}\sum\limits_{i}^{N}\hat{\eta}_{i}^{\mu} \sigma_{i}^{(n+1)}\] \[= \frac{1}{N(1+\rho)}\sum\limits_{i}^{N}\hat{\eta}_{i}^{\mu}\sigma_ {i}^{(n)}\mathrm{sign}\left[\,\sum\limits_{\nu=1}^{K}n_{\nu}^{(n)}(\mathbf{\sigma} )\hat{\eta}_{i}^{\nu}\sigma_{i}^{(n)}\right]\,.\]
If \(\tilde{\mathbf{\sigma}}(\mathbf{\xi})\) is a fixed point of our dynamics, we must have \(\tilde{\mathbf{\sigma}}^{(n+1)}\equiv\tilde{\mathbf{\sigma}}^{(n)}\) and \(n^{(n+1)}(\mathbf{\sigma})\equiv n^{(n)}(\mathbf{\sigma})\), thus (114) becomes

\[n_{\mu}(\mathbf{\sigma})\,=\,\frac{1}{N(1+\rho)}\sum_{i}^{N}\hat{\eta}_{i}^{\mu} \tilde{\sigma}_{i}(\mathbf{\xi})\operatorname{sign}\left[\sum_{\nu=1}^{K}n_{\nu}( \mathbf{\sigma})\hat{\eta}_{i}^{\nu}\tilde{\sigma}_{i}(\mathbf{\xi})\right]. \tag{115}\]

For large value of \(N\), the arithmetic mean coincides with the theoretical expectation, thus

\[\frac{1}{N}\sum_{i=1}^{N}g(\eta_{i})\xrightarrow[N\to\infty]{}\mathbb{E}_{\eta }\Big{[}g(\eta)\Big{]} \tag{116}\]

therefore, (115) reads as

\[n_{\mu}\,=\,\frac{1}{(1+\rho)}\mathbb{E}_{\mathbf{\xi}}\mathbb{E}_{(\mathbf{\eta}|\bm {\xi})}\left\{\hat{\eta}^{\mu}\tilde{\sigma}(\mathbf{\xi})\operatorname{sign}\left[ \,\sum_{\nu=1}^{K}n_{\nu}\hat{\eta}^{\nu}\tilde{\sigma}(\mathbf{\xi})\right]\right\}\,. \tag{117}\]

where we used \(\mathbb{E}_{\mathbf{\eta}}=\mathbb{E}_{\mathbf{\xi}}\mathbb{E}_{(\mathbf{\eta}|\mathbf{\xi})}\). 

**Corollary 2**.: _Under the hypothesis of the previous theorem, if the neural configuration coincides with the parallel configuration_

\[\tilde{\mathbf{\sigma}}(\mathbf{\xi})=\mathbf{\sigma}^{*}=\mathbf{\xi}^{1}+\sum_{\nu=2}^{K} \mathbf{\xi}^{\nu}\prod_{\rho=1}^{\nu-1}\delta\left(\mathbf{\xi}^{\rho}\right) \tag{118}\]

_then the order parameters \(n_{\mu}(\mathbf{\sigma})\) must satisfy the following self equation_

\[n_{\mu}\,=\,\frac{1}{(1+\rho)}\mathbb{E}_{\mathbf{\xi}}\mathbb{E}_{(\mathbf{\eta}|\bm {\xi})}\Bigg{\{}\hat{\eta}^{\mu}\left(\xi^{1}+\sum_{\lambda=2}^{K}\xi^{\lambda} \prod_{\rho=1}^{\lambda-1}\delta\left(\xi^{\rho}\right)\right)\,\operatorname {sign}\left[\sum_{\nu=1}^{K}n_{\nu}\hat{\eta}^{\nu}\left(\xi^{1}+\sum_{\lambda =2}^{K}\xi^{\lambda}\prod_{\rho=1}^{\lambda-1}\delta\left(\xi^{\rho}\right) \right)\right]\Bigg{\}}\,. \tag{119}\]

Proof.: We only have to replace in (109) the explicit form of \(\mathbf{\sigma}^{*}\) and we get the proof. 

### Evaluation of momenta of the effective post-synaptic potential

In this section we want to describe the computation of first and second momenta \(\kappa_{1,\mu}^{(1)}\) and \(\kappa_{2,\mu}^{(1)}\) in Sec. C.1, we will present only the case of \(\mu=1\).

Let us start from \(\kappa_{1,\mu}^{(1)}\):

\[\kappa_{1,\mu}^{(1)}\Big{|}_{\xi_{i}^{1}=\pm 1}\coloneqq\mathbb{E}_{\xi} \mathbb{E}_{(\eta|\xi)}\left[\sigma_{i}^{*}h_{i}^{(1)}(\mathbf{\sigma}^{*})\Big{|} _{\xi_{i}^{1}=\pm 1}\right]=\frac{1}{2}\sum_{j\neq i}^{N}\frac{1}{r_{1}M_{1}^{2}(1+ \rho)}\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\left[\sum_{a,b}^{M_{1},M_{1}} \eta_{j}^{1}\left(\xi_{j}^{1}+\sum_{\nu=2}^{K}\xi_{j}^{\nu}\prod_{\rho=1}^{ \nu-1}\delta_{\xi_{j}^{\rho},0}\right)\right];\]

since \(\mathbb{E}_{\xi}[\xi_{i}^{\mu}]=0\) the only non-zero terms are the ones with \(\mu=1\):

\[\kappa_{1,\mu}^{(1)}\Big{|}_{\xi_{i}^{1}=\pm 1} \coloneqq\,\frac{1}{2}\sum_{j\neq i}^{N}\frac{1}{r_{1}M_{1}^{2}(1+ \rho)}\mathbb{E}_{\xi}\left[\sum_{a,b}^{M_{1},M_{1}}r_{1}\xi_{j}^{1}\left(\xi_ {j}^{1}+\sum_{\nu=2}^{K}\xi_{j}^{\nu}\prod_{\rho=1}^{\nu-1}\delta_{\xi_{j}^{ \rho},0}\right)\right] \tag{120}\] \[=\,\frac{1}{2NM_{1}^{2}r_{1}(1+\rho)}\sum_{j\neq 1}^{N}\sum_{a,b}^{M_{ 1},M_{1}}r_{1}=\frac{1}{2(1+\rho)}\]
where we used \(\mathbb{E}_{(\eta|\xi)}[\eta_{i}^{\mu,a}]=r\xi_{i}^{\mu}\). Moving on, we start the computation of \(\kappa_{2,\mu}^{(1)}\), due to \(\mathbb{E}_{\xi}[\xi_{i_{1}}^{\mu}\xi_{i_{1}}^{\nu}]=\delta^{\mu\nu}\), the only non-zero terms are:

\[\kappa_{2,\mu}^{(1)}\Big{|}_{\xi_{i}^{1}=\pm 1} \coloneqq \mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\left[\{\sigma_{i}^{*}h_{ i}^{(1)}(\mathbf{\sigma}^{*})^{2}\Big{|}_{\xi_{i}^{1}=\pm 1}\right]\] \[\frac{1}{4N^{2}(1-d)^{2}}\sum_{\mu=1}^{K}\mathbb{E}_{\xi} \mathbb{E}_{(\eta|\xi)}\frac{1}{M_{\mu}^{4}r_{\mu}^{4}(1+\rho)^{2}}\sum_{k,j \neq i}^{N,N}\left(\sum_{\alpha_{1},a_{2},b_{1},b_{2}}^{M_{\mu}}\eta_{i}^{\mu,a _{1}}\eta_{i}^{\mu,a_{2}}\eta_{j}^{\mu,b_{1}}\eta_{k}^{\mu,b_{2}}\right)\] \[\left(\xi_{j}^{1}+\sum_{\nu_{1}=2}^{K}\xi_{j}^{\nu_{1}}\prod_{ \rho_{1}=1}^{\nu_{1}-1}\delta\left(\xi_{j}^{\rho_{1}}\right)\right)\left(\xi_ {k}^{1}+\sum_{\nu_{2}=2}^{K}\xi_{\nu_{2}}^{\nu_{2}-1}\prod_{\rho_{2}=1}^{\nu_{ 2}-1}\delta\left(\xi_{k}^{\rho_{2}}\right)\right)=A_{\mu=1}+B_{\mu>1}\]

namely we will analyze separately the case for \(\mu=1\) (\(A_{\mu=1}\)) and \(\mu>1\) (\(B_{\mu>1}\)).

\[A_{\mu=1} = \frac{1}{4N^{2}(1-d)^{2}}\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)} \frac{1}{M_{1}^{4}r_{1}^{4}(1+\rho)^{2}}\sum_{k,j\neq i}^{N,N}\left(\sum_{a_{1 },a_{2},b_{1},b_{2}}^{M_{1}}\eta_{i}^{1,a_{2}}\eta_{j}^{1,b_{1}}\eta_{k}^{1,b_{ 2}}\right)\left(\xi_{j}^{1}\xi_{k}^{1}\right) \tag{101}\] \[= \frac{1}{4N^{2}(1-d)^{2}}\frac{1}{M_{1}^{4}r_{1}^{4}(1+\rho)^{2} }\mathbb{E}_{(\eta|\xi=\pm 1)}\left[\sum_{a_{1},a_{2}}^{M_{1}}\eta_{i}^{1,a_{1}} \eta_{i}^{1,a_{2}}\right]\sum_{k,j\neq i}^{N,N}\mathbb{E}_{\xi}\mathbb{E}_{( \eta|\xi)}\left[\sum_{b_{1},b_{2}}^{M_{1}}\eta_{j}^{1,b_{1}}\eta_{k}^{1,b_{2} }\right]\left(\xi_{j}^{1}\xi_{k}^{1}\right)\] \[= \frac{1}{4}\frac{1}{(1+\rho)}\,;\]

\[B_{\mu>1} \coloneqq \frac{1}{4N^{2}(1-d)^{2}}\sum_{\mu>1}^{K}\mathbb{E}_{\xi} \mathbb{E}_{(\eta|\xi)}\frac{1}{M_{\mu}^{4}r_{\mu}^{4}(1+\rho)^{2}}\sum_{k,j \neq i}^{N,N}\left(\sum_{\alpha_{1},a_{2},b_{1},b_{2}}^{M_{\mu}}\eta_{i}^{\mu, a_{1}}\eta_{i}^{\mu,a_{2}}\eta_{j}^{\mu,b_{1}}\eta_{k}^{\mu,b_{2}}\right)\] \[\left(\prod_{\rho_{1}=1}^{\mu-1}\delta\left(\xi_{j}^{\rho_{1}} \right)\right)\left(\prod_{\rho_{2}=1}^{\mu-1}\delta\left(\xi_{k}^{\rho_{2}} \right)\right)\] \[= \frac{(1-d)^{-2}}{4N^{2}}\sum_{\mu=2}^{K}\frac{1}{M_{\mu}^{2}(1+ \rho)}\sum_{k,j\neq i}^{N,N}\sum_{b_{1},b_{2}}^{M_{\mu}}\mathbb{E}_{\xi} \left[\xi_{j}^{\mu}\xi_{k}^{\mu}\left(\prod_{\rho_{1}=1}^{\mu-1}\delta\left( \xi_{j}^{\rho_{1}}\right)\right)\left(\prod_{\rho_{2}=1}^{\mu-1}\delta\left( \xi_{k}^{\rho_{2}}\right)\right)\right]\] \[= \frac{1-d}{4(1+\rho)}\sum_{\mu=2}^{K}d^{2(\mu-1)}=\frac{1-d}{4(1+ \rho)}\frac{d^{2}-d^{2K}}{1-d^{2}}=\frac{1}{4(1+\rho)}\frac{d^{2}-d^{2K}}{1+d}\,.\]

Putting together Eq. (101) and Eq. (100) we reach the expression of \(\kappa_{2,\mu}^{(1)}\Big{|}_{\xi_{i}^{1}=\pm 1}\).

## Appendix D Explicit Calculations and Figures for the cases \(K=2\) and \(K=3\)

In this appendix, we collect the explicit expression of the self-consistent equations in (3.10) and (3.11) (focusing only on the cases of \(K=2\) and \(K=3\)) and some Figures obtained from their numerical solution.
### \(K=2\)

Fixing \(K=2\) and explicitly perform the mean with respect to \(\xi\), (3.10) and (3.11) read as

\[\begin{array}{rcl}\bar{n}_{1}&=&\frac{\bar{m}_{1}}{(1+\rho)}+\frac{\beta^{{}^{ \prime}}(1-d)\rho\,\bar{n}_{1}}{(1+\rho)}\bigg{[}1-d\,\mathcal{S}_{2}(\bar{n}_ {1},0)-\frac{(1-d)}{2}\mathcal{S}_{2}(\bar{n}_{1},-\bar{n}_{2})-\frac{(1-d)}{2 }\mathcal{S}_{2}(\bar{n}_{1},\bar{n}_{2})\bigg{]}\\ \\ \bar{n}_{2}&=&\frac{\bar{m}_{2}}{(1+\rho)}+\frac{\beta^{{}^{\prime}}(1-d)\rho \,\bar{n}_{2}}{(1+\rho)}\bigg{[}1-d\,\mathcal{S}_{2}(0,\bar{n}_{2})-\frac{(1-d) }{2}\mathcal{S}_{2}(\bar{n}_{1},-\bar{n}_{2})-\frac{(1-d)}{2}\mathcal{S}_{2}( \bar{n}_{1},\bar{n}_{2})\bigg{]}\\ \\ \bar{m}_{1}&=&\frac{(1-d)^{2}}{2}\bigg{[}\mathcal{T}_{2}(\bar{n}_{1},\bar{n}_ {2})+\mathcal{T}_{2}(\bar{n}_{1},-\bar{n}_{2})\bigg{]}+d(1-d)\mathcal{T}_{2}( \bar{n}_{1},0)\\ \\ \bar{m}_{1}&=&\frac{(1-d)^{2}}{2}\bigg{[}\mathcal{T}_{2}(\bar{n}_{1},\bar{n}_ {2})-\mathcal{T}_{2}(\bar{n}_{1},-\bar{n}_{2})\bigg{]}+d(1-d)\mathcal{T}_{2}(0,\bar{n}_{2})\end{array}\] (D.1)

where we used

\[\begin{array}{rcl}\mathcal{T}_{2}(x,y)&=&\mathbb{E}_{\lambda}\tanh\left[ \beta^{{}^{\prime}}\left(x+y+\lambda\sqrt{\rho\Big{(}x^{2}+y^{2}\Big{)}} \right)\right],\\ \\ \mathcal{S}_{2}(x,y)&=&\mathbb{E}_{\lambda}\tanh^{2}\left[\beta^{{}^{\prime}} \left(x+y+\lambda\sqrt{\rho\Big{(}x^{2}+y^{2}\Big{)}}\right)\right].\end{array}\] (D.2)

Solving numerically this set of equations we construct the plots presented in Fig.10.

### \(K=3\)

Moving on the case of \(K=3\) by following the same steps of the previous subsection, we get

\[\bar{n}_{1} = \frac{\bar{m}_{1}}{(1+\rho)}+\frac{\beta^{{}^{\prime}}(1-d)\rho\, \bar{n}_{1}}{(1+\rho)}\bigg{\{}1-d\,\frac{(1-d)}{2}\left[\mathcal{S}_{3}(\bar{ n}_{1},\bar{n}_{2},0)+\mathcal{S}_{3}(\bar{n}_{1},0,\bar{n}_{3})+\mathcal{S}_{3}( \bar{n}_{1},-\bar{n}_{2},0)+\mathcal{S}_{3}(\bar{n}_{1},0,-\bar{n}_{3})\right]\] \[-d^{2}\mathcal{S}_{3}(\bar{n}_{1},0,0)-\frac{(1-d)^{2}}{4}\left[ \mathcal{S}_{3}(\bar{n}_{1},\bar{n}_{2},\bar{n}_{3})+\mathcal{S}_{3}(\bar{n}_ {1},\bar{n}_{2},-\bar{n}_{3})+\mathcal{S}_{3}(\bar{n}_{1},-\bar{n}_{2},\bar{n} _{3})+\mathcal{S}_{3}(\bar{n}_{1},-\bar{n}_{2},-\bar{n}_{3})\right]\bigg{\}}\,,\] \[\bar{m}_{1} = \frac{(1-d)^{3}}{4}\bigg{[}\mathcal{T}_{3}(\bar{n}_{1},\bar{n}_{ 2},\bar{n}_{3})+\mathcal{T}_{3}(\bar{n}_{1},\bar{n}_{2},-\bar{n}_{3})+ \mathcal{T}_{3}(\bar{n}_{1},-\bar{n}_{2},\bar{n}_{3})+\mathcal{T}_{3}(\bar{n}_ {1},-\bar{n}_{2},-\bar{n}_{3})\bigg{]}\] \[+d\frac{(1-d)^{2}}{2}\bigg{[}\mathcal{T}_{3}(\bar{n}_{1},\bar{n}_ {2},0)+\mathcal{T}_{3}(\bar{n}_{1},0,\bar{n}_{3})+\mathcal{T}_{3}(\bar{n}_{1},- \bar{n}_{2})+\mathcal{T}_{3}(\bar{n}_{1},0,-\bar{n}_{3})\bigg{]}+d^{2}(1-d) \mathcal{T}_{3}(\bar{n}_{1},0,0)\,,\]

where we used

\[\begin{array}{rcl}\mathcal{T}_{3}(x,y,z)&=&\mathbb{E}_{\lambda}\tanh\left[ \beta^{{}^{\prime}}\left(x+y+z+\lambda\sqrt{\rho\Big{(}x^{2}+y^{2}+z^{2}\Big{)} }\right)\right]\,,\\ \\ \mathcal{S}_{3}(x,y,z)&=&\mathbb{E}_{\lambda}\tanh^{2}\left[\beta^{{}^{\prime}} \left(x+y+z+\lambda\sqrt{\rho\Big{(}x^{2}+y^{2}+z^{2}\Big{)}}\right)\right]\,. \end{array}\] (D.4)

In order to lighten the presentation, we reported only the expression of \(\bar{m}_{1}\) and \(\bar{n}_{1}\), the related expressions of \(\bar{m}_{2}(\bar{m}_{3})\) and \(\bar{n}_{2}(\bar{n}_{3})\) can be obtained by making the simple substitutions \(\bar{m}_{1}\longleftrightarrow\bar{m}_{2}(\bar{m}_{3})\), and \(\bar{n}_{1}\longleftrightarrow\bar{n}_{2}(\bar{n}_{3})\) in (D.3). The numerical solution of the previous set of equations is depicted in Fig.11.
## Appendix E Proofs

### Proof of Theorem 1

In this subsection we show the proof Proposition 1. In order to prove the aforementioned proposition, we put in front of it the following

**Lemma 1**.: _The \(t\) derivative of interpolating pressure is given by_

\[\frac{d\mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}}{dt} = \frac{\beta}{2(1-d)}(1+\rho)\sum_{\mu=1}^{K}\mathbb{E}\,\omega_{ t}[n_{\mu}^{2}]-\sum_{\mu=1}^{K}\psi_{\mu}\mathbb{E}\,\omega_{t}[n_{\mu}]. \tag{100}\]

Since the computation is lengthy but not cumbersome we decided to omit it.

**Proposition 3**.: _In the low load regime, in the thermodynamic limit the distribution of the generic order parameter \(X\) is centres at its expectation value \(\bar{X}\) with vanishing fluctuations. Thus, being

Figure 10: Numerical resolution of the system of equations (100) for \(K=2\): we plot the behaviour of the magnetization \(\tilde{m}\) versus degree of dilution \(d\) for fixed \(r=0.2\) and different value of \(\beta\) (from right to left \(\beta=1000,6.66,3.33\)) and \(\rho\) (from top to bottom \(\rho=0.8,0.2,0.0\)). We stress that for \(\rho=0.0\) we recover the standard diluted model presented in Fig.1.
\(\Delta X=X-\bar{X}\), in the thermodynamic limit, the following relation holds_

\[\mathbb{E}\,\omega_{t}[(\Delta X)^{2}]\xrightarrow[N\to+\infty]{}0\,. \tag{100}\]

**Remark 5**.: _We stress that afterwards we use the relations_

\[\mathbb{E}\,\omega_{t}[(n_{\mu}-\bar{n}_{\mu})^{2}]=\mathbb{E}\,\omega_{t}[n_{ \mu}^{2}]-2\,\bar{n}_{\mu}\mathbb{E}\,\omega_{t}[n_{\mu}]+\bar{n}_{\mu}^{2}\,. \tag{101}\]

_which are computed with brute force with Newton's Binomial._

_Now, using these relations, if we fix the constants as_

\[\psi_{\mu}=\frac{\beta}{1-d}(1+\rho)\bar{n}_{\mu} \tag{102}\]

_in the thermodynamic limit, due to Proposition 3, the expression of derivative w.r.t. \(t\) becomes_

\[\frac{d\mathcal{A}^{(sup,unsup)}_{K,\beta,d,M,r}}{dt} = -\frac{\beta}{2(1-d)}(1+\rho)\sum_{\mu=1}^{K}\bar{n}_{\mu}^{2}. \tag{103}\]

Proof.: Let us start from finite size \(N\) expression. We apply the Fundamental Theorem of Calculus:

\[\mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}=\mathcal{A}^{(sup,unsup)}_{N,K, \beta,d,M,r}(t=1)=\mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}(t=0)+\int\limits _{0}^{1}\left.\partial_{s}\mathcal{A}^{(sup,unsup)}_{N,K,\beta,d,M,r}(s) \right|_{s=t}dt. \tag{104}\]

Figure 11: Numerical solution of the system of equations (101) for \(K=3\): we plot the behavior of the magnetization \(\bar{m}\) versus the degree of dilution \(d\) for fixed \(r=0.2\) and different value of \(\beta\) (from left to right \(\beta=1000,6.66,3.33\)) and \(\rho\) (from top to bottom \(\rho=0.8,0.2,0.0\)).
We have already computed the derivative w.r.t. \(t\) in Eq. (102). It only remains to calculate the one-body term:

\[\mathcal{Z}^{(sup,unsup)}_{N,K,\beta,d,M,r}(t=0)=\sum_{\{\sigma\}}\exp\Bigg{[} \sum_{i=1}^{N}\left(\sum_{\mu=1}^{K}\frac{\psi_{\mu}}{2(1+\rho)}\hat{\eta}^{\mu }+J\xi^{\mu}\right)\sigma_{i}\Bigg{]}. \tag{103}\]

Using the definition of quenched statistical pressure (100) we have

\[\mathcal{A}^{(sup,unsup)}_{K,\beta,d,M,r}(J,t=0)=\ln\,\left[2\cosh \left(\sum_{\mu=1}^{K}\frac{\psi_{\mu}}{2(1+\rho)}\hat{\eta}^{\mu}+J\xi^{\mu} \right)\right] \tag{104}\] \[=\mathbb{E}\left\{\ln 2\cosh\left[\frac{\beta}{1-d}\sum_{\mu=1}^{K }\bar{n}_{\mu}\hat{\eta}^{\mu}+J\xi^{\mu}\right]\right\}\]

where \(\mathbb{E}=\mathbb{E}_{\xi}\mathbb{E}_{(\eta|\xi)}\). Finally, putting inside (100) (104) and (102), we reach the thesis. 

### Proof of Proposition 1

In this subsection we show the proof Proposition 1.

Proof.: For large data-sets, using the Central Limit Theorem we have

\[\hat{\eta}^{\mu}\sim\xi^{\mu}\left(1+\sqrt{\rho}\;Z_{\mu}\right)\,. \tag{105}\]

where \(Z_{\mu}\) is a standard Gaussian variable \(Z_{\mu}\sim\mathcal{N}(0,1)\). Replacing Eq. (105) in the self-consistency equation for \(\bar{n}\), namely Eq. (101), and applying Stein's lemma9

Footnote 9: This lemma, also known as Wick’s theorem, applies to standard Gaussian variables, say \(J\sim\mathcal{N}(0,1)\), and states that, for a generic function \(f(J)\) for which the two expectations \(\mathbb{E}\left(Jf(J)\right)\) and \(\mathbb{E}\left(\partial_{J}f(J)\right)\) both exist, then

\[\mathbb{E}\left(Jf(J)\right)=\mathbb{E}\left(\frac{\partial f(J)}{\partial J} \right)\,. \tag{106}\]

 in order to recover the expression for \(\bar{m}_{\mu}\), we get the large data-set equation for \(\bar{n}_{\mu}\), i.e. Eq. (104).

We will use the relation

\[\mathbb{E}_{\lambda_{\mu}}\left[F\left(a+\sum_{\mu=1}^{K}b_{\mu}\lambda_{\mu} \right)\right]=\mathbb{E}_{x}\,\left[F\left(a+Z\sqrt{\sum_{\mu=1}^{K}b_{\mu} ^{2}}\right)\right]\,, \tag{107}\]

where \(\lambda_{\mu}\) and \(Z\) are i.i.d. Gaussian variables. Doing so, we obtain

\[g(\beta,\boldsymbol{\xi},Z,\bar{n})=\beta^{{}^{\prime}}\sum_{\nu=1}^{K}\bar{n }_{\nu}\xi^{\nu}+\beta^{{}^{\prime}}\sqrt{\rho}\sum_{\nu=1}^{K}Z_{\nu}\bar{n} _{\nu}^{2}\left(\xi^{\nu}\right)^{2}=\beta^{{}^{\prime}}\left(\sum_{\nu=1}^{K }\bar{n}_{\nu}\xi^{\nu}+Z\sum_{\nu=1}^{K}\sqrt{\rho\,\bar{n}_{\nu}^{2}\left( \xi^{\nu}\right)^{2}}\right)\,, \tag{108}\]

thus we reach the thesis. 

**Corollary 3**.: _The self consistency equations in the large data-set assumption and null-temperature limit are_

\[\bar{m}_{\mu}\,=\,\mathbb{E}_{\xi}\left\{\operatorname{erf}\left[\left(\sum_{ \nu=1}^{K}\bar{m}_{\nu}\xi^{\nu}\right)\left(2\rho\sum_{\nu=1}^{K}\bar{m}_{\nu }^{2}\left(\xi^{\nu}\right)^{2}\right)^{-1/2}\right]\xi^{\mu}\right\}. \tag{109}\]
Proof.: In order to lighten the notation we rename

\[C=\tanh^{2}\left[g(\beta,\boldsymbol{\xi},Z,\bar{\boldsymbol{n}})\right]\,. \tag{101}\]

We start by assuming finite the limit

\[\lim_{\beta^{{}^{\prime}}\to\infty}\beta^{{}^{\prime}}(1-C)=D\in\mathbb{R} \tag{102}\]

and we stress that as \(\beta^{{}^{\prime}}\to\infty\) we have \(C\to 1\). As a consequence, the following reparametrization is found to be useful,

\[C=1-\frac{\delta C}{\beta^{{}^{\prime}}}\quad\text{as}\quad\beta^{{}^{\prime} }\to\infty. \tag{103}\]

Therefore, as \(\beta^{{}^{\prime}}\to\infty\), it yields

\[\bar{n}_{\mu}=\frac{\bar{m}_{\mu}}{1+\rho-\rho\,\delta C\,(1-d)} \tag{104}\]

to reach this result, we have also used the relation

\[\mathbb{E}_{z}\text{sign}[A+Bz]=\text{erf}\left[\frac{A}{\sqrt{2}B}\right]\, \tag{105}\]

where \(z\) is a Gaussian variable \(\mathcal{N}(0,1)\) and the truncated expression \(\bar{n}_{\mu}=\bar{m}_{\mu}/(1+\rho)\) for the first equation in (104). 

This research has been supported by Ministero degli Affari Esteri e della Cooperazione Internazionale (MAECI) via the BULBUL grant (Italy-Israel), CUP Project n. F85F21006230001, and has received financial support from the Simons Foundation (grant No. 454949, G. Parisi) and ICSC - Italian Research Center on High Performance Computing, Big Data and Quantum Computing, funded by European Union - NextGenerationEU.

Further, this work has been partly supported by The Alan Turing Institute through the Theory and Methods Challenge Fortnights event _Physics-informed machine learning_, which took place on 16-27 January 2023 at The Alan Turing Institute headquarters.

E.A. acknowledges financial support from Sapienza University of Rome (RM120172B8066CB0).

E.A., A.A. and A.B. acknowledges GNFM-INdAM (Gruppo Nazionale per la Fisica Mamematica, Istituto Nazionale d'Alta Matematica), A.A. further acknowledges UniSalento for financial support via PhD-AI and AB further acknowledges the PRIN-2022 project _Statistical Mechanics of Learning Machines: from algorithmic and information-theoretical limits to new biologically inspired paradigms_.
