# Supervised Hebbian learning: toward eXplainable AI
Supplementary Material

##  I Introduction

In this Supplementary Material we discuss in details the analytical and
computational techniques underlying the results presented in the main text.
In particular, Secs. II and III are dedicated to the analysis of the Hopfield
network (HN) equipped with the Hebbian learning rule (bio-inspired computing):
in the former we revise the signal-to-noise approach, while in the latter we
detail the statistical mechanics approach. Next, Secs. IV and V deal with
Restricted Boltzmann machines (RBMs, artificial computing): in the former we
prove that the scaling laws for a successful learning by the HN hold also for
the RBM, while in the latter we discuss why and how these two (apparently
distant) learning schemes can be related, tackling their equivalence from
Jayne’s maximum-entropy perspective [1].
Before starting, it is useful to recall the basic definitions concerning the
random dataset and the HN cost functions implementing the unsupervised and the
supervised Hebbian learning.

######  Definition 1.

The _archetype_ dataset is made of $K$ binary vectors of length $N$, denoted
as $\boldsymbol{\xi}^{\mu}=(\xi_{1}^{\mu},...,\xi_{N}^{\mu})$, for
$\mu=1,...,K$, whose entries are i.i.d. Rademacher random variables

$\mathcal{P}(\xi_{i}^{\mu}=+1)=\mathcal{P}(\xi_{i}^{\mu}=-1)=\frac{1}{2},$ (1)

for any $i=1,...,N$ and $\mu=1,...,K$.
The dataset available to the machine, denoted as
$\mathcal{S}:=\\{\boldsymbol{\eta}^{\mu a}\\}_{\mu=1,...,K}^{a=1,...,M}$, is
made of $M$ corrupted _examples_ of each archetype: $\boldsymbol{\eta}^{\mu
a}$ is the $a$-th example of the $\mu$-th archetype and its entries are
defined as

$\eta_{i}^{\mu a}=\chi_{i}^{\mu a}\xi_{i}^{\mu}$ (2)

where $\chi_{i}^{\mu a}$ is a symmetric Bernoullian random variable drawn as

$\mathcal{P}(\chi_{i}^{\mu a}=\pm 1)=\frac{1\pm r}{2},r\in[0,1]$ (3)

for any $i=1,...,N$, $\mu=1,...,K$, and $a=1,...,M$.

######  Definition 2.

The averaging with respect to $\boldsymbol{\xi}^{\mu}$ and
$\boldsymbol{\chi}^{\mu a}$ are

$\displaystyle\mathbb{E}_{\xi^{\mu}}f(\boldsymbol{\xi}^{\mu})$
$\displaystyle:=$
$\displaystyle\int\prod_{i=1}^{N}\left\\{\frac{d\xi_{i}^{\mu}}{2}\left[\delta(\xi_{i}^{\mu}-1)+\delta(\xi_{i}^{\mu}+1)\right]\right\\}f(\boldsymbol{\xi}^{\mu}),$
(4) $\displaystyle\mathbb{E}_{\chi^{\mu}}f(\boldsymbol{\chi}^{\mu})$
$\displaystyle:=$
$\displaystyle\int\prod_{i,a=1}^{N,M}\left\\{\frac{d\chi_{i}^{\mu
a}}{2}\left[(1+r)\delta(\chi_{i}^{\mu a}-1)+(1-r)\delta(\chi_{i}^{\mu
a}+1)\right]\right\\}f(\boldsymbol{\chi}^{\mu}),$ (5)
$\displaystyle\mathbb{E}_{\xi}$ $\displaystyle:=$
$\displaystyle\prod_{\mu=1}^{K}\mathbb{E}_{\xi^{\mu}},\qquad\mathbb{E}_{\chi}:=\prod_{\mu=1}^{K}\mathbb{E}_{\chi^{\mu}},$
(6)

where $f$ is a generic function.

######  Definition 3.

The Hamiltonian of the _unsupervised_ Hopfield model is defined as

$\mathcal{H}^{\textrm{\tiny{HN},unsup}}(\boldsymbol{\sigma}|\mathcal{S}):=-\frac{N}{2M\mathcal{R}}\sum_{\mu=1}^{K}\sum_{a=1}^{M}\left(\frac{1}{N}\sum_{i=1}^{N}\eta_{i}^{\mu
a}\sigma_{i}\right)^{2},$ (7)

while the Hamiltonian of the _supervised_ Hopfield model is defined as

$\mathcal{H}^{\textrm{\tiny{HN},sup}}(\boldsymbol{\sigma}|\mathcal{S}):=-\frac{N}{2\mathcal{R}}\sum_{\mu=1}^{K}\left(\frac{1}{NM}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu
a}\sigma_{i}\right)^{2}$ (8)

where $\sigma_{i}\in\\{-1,+1\\}$ for $i=1,\cdots,N$ is a binary neuron (i.e.,
an Ising spin) and $\mathcal{R}:=r^{2}+\frac{1-r^{2}}{M}$ is a normalization
factor.

######  Remark 1.

In the physical model under investigation self-interaction terms are excluded,
that is, a neuron $\sigma_{i}$ interacts with any other neuron $\sigma_{k}$
with $k\neq i$. This should be accounted for in (7) and (8) by inserting a
corrective contribution that neutralizes diagonal terms, however, since this
contribution is constant, here it is neglected in order to retain the notation
simple.

##  II Signal-to-noise Analysis

By the signal-to-noise technique we find the conditions for the HN to
successfully generalize from examples, namely to retrieve one of the
archetypes although it was never presented to them, in the noiseless
$\beta\to\infty$ limit. The idea is to assume that the network is in a
retrieval state of an archetype, say $\boldsymbol{\xi}^{1}$, evaluate the
internal fields $h_{i}(\boldsymbol{\xi}^{1})$ acting on each neuron in such a
configuration, and check that the constraints for stability
$h_{i}(\boldsymbol{\xi}^{1})\sigma_{i}>0$ are all satisfied. The analysis is
led in the high-storage regime, where the load $\alpha:=\lim_{N\to\infty}K/N$
is finite. Clearly, we expect an interplay between $\alpha$, $M$ and $r$ as,
by raising $\alpha$ we need a large and good-quality dataset (i.e., large $M$
and $r$) to disentangle examples, and, for a given load, by reducing $r$ we
need a larger $M$ to retain enough information in the dataset. In fact, as
proved in the following subsections II.1 and II.2, respectively, the optimal
tradeoff in the unsupervised regime reads as

$\alpha\left(1+\frac{1-r^{4}}{Mr^{4}}\right)+\frac{1-r^{2}}{Mr^{2}}\lesssim
1,$ (9)

while in the supervised regime it reads as

$\alpha\left(1+\frac{1-r^{2}}{Mr^{2}}\right)^{2}+\frac{1-r^{2}}{Mr^{2}}\lesssim
1.$ (10)

Before proceeding, we anticipate that these scalings are in excellent
agreement with Monte Carlo (MC) simulations (see Fig. S1) and obviously in
full accordance with the statistical mechanical predictions reported in the
Sec. III. Moreover, in the present treatment, the standard signal-to-noise
approach is recast into the one-step Hinton’s prescription for fast MC
sampling [2]: the reward in this reformulation of the signal-to-noise
technique is that it is predictive also for these numerical algorithms that,
in turn, we implement in the simulations presented in Sec. IV, see Eq. (118).

Figure S1: Comparison between signal-to-noise predictions (solid lines) and MC
simulations (symbols) for the unsupervised (dark color, $\circ$) and the
supervised (bright color, $\square$) Hebbian learning in the noiseless
$\beta\to\infty$ limit. More precisely, the theoretical estimates are obtained
by eq. (9) and eq. (10), while the numerical estimates are obtained by fixing
a certain number $M$ of examples and determining the minimal value of the
dataset quality $r$ such that the mean overlap between the neural
configuration $\boldsymbol{\sigma}$ and a chosen archetype, say
$\boldsymbol{\xi}^{1}$, i.e., the magnetization $\langle m\rangle$, is
(approximately) unitary, the operation is then repeated for several values of
$M$ ranging exponentially from $2$ to $2^{12}$. Note that the theoretical
estimates perfectly interpolate the numerical outcomes.

###  II.1 Signal vs Noise Analysis: Unsupervised Hebbian Learning

The unsupervised case has been previously addressed by signal-to-noise
technique in [3] whose results are here recovered and refined, further, the
previous investigation is now reformulated in analogy with Hinton’s one-step
MC recipe [2].

######  Proposition 1.

In the unsupervised setting, the condition determining if the HN successfully
retrieves one of the archetypes is

$\frac{1}{\sqrt{2\left[\alpha\left(1+\frac{1-r^{4}}{r^{4}M}\right)+\frac{1-r^{2}}{Mr^{2}}\right]}}>\theta$
(11)

where $\theta\in\mathbb{R}^{+}$ is a tolerance level, $\textrm{erf}(\theta)$
representing a lower-bound for the overlap between the neural configuration
and the retrieved archetype (we omit plots in this case as the results are
qualitatively the same as those already reported in [3]).

###### Proof.

First we recall the unnormalized log-density of the model:

$-\beta\mathcal{H}^{\textrm{\tiny{HN},
unsup}}(\boldsymbol{\sigma}|\boldsymbol{\eta})=\frac{\beta
N}{2M}\sum_{\mu,a=1}^{K,M}\left(\frac{1}{N}\sum_{i=1}^{N}\eta_{i}^{\mu
a}\sigma_{i}\right)^{2}=\frac{\beta}{2NM}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\sum_{a=1}^{M}\eta_{i}^{\mu
a}\eta_{j}^{\mu
a}\sigma_{i}\sigma_{j}=\frac{1}{2}\sum_{k=1}^{N}h_{k}(\boldsymbol{\sigma})\sigma_{k}$
(12)

where

$h_{k}(\boldsymbol{\sigma}):=\frac{\beta}{NM}\sum_{\mu=1}^{K}\sum_{i\neq
k}^{N}\sum_{a=1}^{M}\eta_{i}^{\mu a}\eta_{k}^{\mu a}\sigma_{i}$ (13)

is the internal field acting on $\sigma_{k}$.
The discrete time MC dynamics is given by the following update rule

$\sigma_{k}^{(n+1)}=\sigma_{k}^{(n)}\mathrm{sign}\left\\{\tanh\left[h_{k}^{(n)}(\boldsymbol{\sigma}^{(n)})\sigma_{k}^{(n)}\right]+\zeta_{k}^{(n)}\right\\},\quad\zeta_{k}^{(n)}\sim\mathcal{U}(-1,+1).$
(14)

where $n\in\mathbb{N}$, $k=1,\dots,N$ and $\mathcal{U}(-1,1)$ is the uniform
probability density in the interval $(-1,+1)$.
By performing the zero fast-noise limit $\beta\to\infty$, equation (14)
becomes

$\sigma_{k}^{(n+1)}=\sigma_{k}^{(n)}\mathrm{sign}\left(h_{k}^{(n)}(\boldsymbol{\sigma}^{(n)})\sigma_{k}^{(n)}\right).$
(15)

Since we want to study the stability of the archetype retrieval, we set as the
initial configuration for the MC dynamics $\sigma_{k}^{(1)}=\xi_{k}^{1}$ for
$k=1,...,N$. The one-step MC approximation for the magnetization is then

$m_{1}^{(2)}=\frac{1}{N}\sum_{k=1}^{N}\xi_{k}^{1}\sigma_{k}^{(2)}=\frac{1}{N}\sum_{k=1}^{N}\mathrm{sign}\left(h_{k}^{(1)}(\boldsymbol{\xi}^{1})\xi_{k}^{1}\right)$
(16)

and, if $N\gg 1$, by the central limit theorem, equation (16) can be
approximated as

$m_{1}^{(2)}\underset{N\gg
1}{\approx}\int\frac{dz}{\sqrt{2\pi}}\exp\left(-\frac{z^{2}}{2}\right)\mathrm{sign}\left(\mu_{1}+z\sqrt{\mu_{2}-\mu_{1}^{2}}\right)=\mathrm{erf}\left(\frac{\mu_{1}}{\sqrt{2(\mu_{2}-\mu_{1}^{2})}}\right),$
(17)

where we posed

$\displaystyle\mu_{1}$ $\displaystyle:=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}^{(n)}(\boldsymbol{\xi}^{1})\xi_{k}^{1}\right],$
(18) $\displaystyle\mu_{2}$ $\displaystyle:=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}^{(n)}(\boldsymbol{\xi}^{1})^{2}\right].$
(19)

First, we calculate $\mu_{1}$:

$\mu_{1}=\frac{\beta}{NM}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu=1}^{K}\sum_{i\neq
k}^{N}\sum_{a=1}^{M}\xi_{k}^{\mu}\chi_{k}^{\mu a}\xi_{i}^{\mu}\chi_{i}^{\mu
a}\xi_{i}^{1}\xi_{k}^{1}\right]=\frac{\beta}{NM}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu>1}^{K}\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu
a}\chi_{i}^{a\mu}+\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1a}\right].$ (20)

Exploiting the fact that, by construction, archetypes are orthogonal in the
average, that is $\mathbb{E}[\boldsymbol{\xi}^{\mu}\boldsymbol{\xi}^{\nu}]=0$
for any $\mu\neq\nu$,

$\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu>1}^{K}\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu a}\right]=0,$
(21)

and, recalling (3),

$\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1a}\right]=(N-1)Mr^{2},$ (22)

we get

$\mu_{1}=\frac{\beta}{N}(N-1)r^{2}\underset{N\gg 1}{\approx}\mu_{1}=\beta
r^{2}.$ (23)

Then, we calculate $\mu_{2}$:

$\displaystyle\mu_{2}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}(\boldsymbol{\boldsymbol{\xi}^{1}})^{2}\right]=\frac{\beta^{2}}{N^{2}M^{2}}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\left(\sum_{\mu>1}^{K}\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu
a}+\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1a}\right)^{2}\right]=$ (24) $\displaystyle=$
$\displaystyle\frac{\beta^{2}}{\mathcal{R}^{2}N^{2}M^{4}}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\left(\sum_{\mu>1}^{K}\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu
a}\right)^{2}+\left(\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1a}\right)^{2}\right].$ (25)

For the sake of simplicity we pose
$A:=\left(\sum_{\mu>1}^{K}\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu a}\right)^{2}$
and $B:=\left(\sum_{a=1}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1a}\right)^{2}$such that

$\mu_{2}=\frac{\beta^{2}}{N^{2}M^{2}}\left(\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[A\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[B\right]\right).$
(26)

After minimal manipulations we can write these terms as

$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[A\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\nu>1}\sum_{a,b}\sum_{i\neq
k}\chi_{k}^{\mu a}\chi_{i}^{\mu a}\chi_{k}^{\mu b}\chi_{i}^{\mu b}\right],$
(27) $\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[B\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1a}\chi_{k}^{1b}\chi_{i}^{1b}\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b}^{M}\sum_{i,j\neq
k}^{N}\lambda_{ij}\chi_{k}^{1a}\chi_{i}^{1a}\chi_{k}^{1b}\chi_{j}^{1b}\right],$
(28)

where $\lambda_{ij}:=1-\delta_{ij}$ . By merging the results in (27) and (28)
we get

$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[A\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[B\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\nu=1}\sum_{a,b}\sum_{i\neq
k}\chi_{k}^{\mu a}\chi_{i}^{\mu a}\chi_{k}^{\mu b}\chi_{i}^{\mu
b}\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b}^{M}\sum_{i,j\neq
k}^{N}\lambda_{ij}\chi_{k}^{1a}\chi_{i}^{1a}\chi_{k}^{1b}\chi_{j}^{1b}\right]=$
(29) $\displaystyle=$ $\displaystyle
K(N-1)\left[M+r^{4}M(M-1)\right]+(N-1)(N-2)\left[r^{2}M+r^{4}M(M-1)\right]$
(30) $\displaystyle\underset{N\gg 1}{\approx}$ $\displaystyle
KNr^{4}M^{2}\left(1+\frac{1-r^{4}}{r^{4}M}\right)+N^{2}r^{4}M^{2}\left(1+\frac{1-r^{2}}{r^{2}M}\right)$
(31)

and, by direct substitution of (31) into (26), we get

$\mu_{2}\underset{N\gg
1}{\approx}\beta^{2}r^{4}\left[\alpha\left(1+\frac{1-r^{4}}{r^{4}M}\right)+1+\frac{1-r^{2}}{r^{2}M}\right].$
(32)

Finally, by plugging equations (23) and (32) into the expression of the
magnetization (17), we get

$m_{1}^{(2)}\underset{N\gg
1}{\approx}\mathrm{erf}\left(\frac{1}{\sqrt{2\left[\alpha\left(1+\frac{1-r^{4}}{r^{4}M}\right)+\frac{1-r^{2}}{Mr^{2}}\right]}}\right)$
(33)

and, by requiring that this one-step MC magnetization is larger than
$\textrm{erf}(\theta)$ we recover Eq. (11). ∎

###  II.2 Signal vs Noise Analysis: Supervised Hebbian learning

We now repeat analogous passages for the supervised setting that lead to

######  Proposition 2.

In the supervised setting, the condition determining if the HN successfully
retrieves one of the archetypes is

$\frac{1}{\sqrt{2\alpha\left(1+\rho\right)^{2}+2\rho}}>\theta$ (34)

where $\theta\in\mathbb{R}^{+}$ is a tolerance level, $\textrm{erf}(\theta)$
representing a lower-bound for the overlap between the neural configuration
and the retrieved archetype (for instance, the case $\theta=1/\sqrt{2}$ is
shown in Fig. S2).

Figure S2: Signal-to-Noise analysis. Left panel: Contour plot for the one-
step-MC magnetization $m^{(2)}$ defined in (40) and evaluated by MC
simulations, run at different values of $\alpha$ and $\rho$; the solid line
corresponds to the curve $\alpha\left(1+\rho\right)^{2}+\rho=1$ obtained by
eq. (34) by setting as tolerance level $\theta=\frac{1}{\sqrt{2}}$. Right
panel: Contour plot for the one-step-MC magnetization $m^{(2)}$ defined in
(40) and evaluated by MC simulations, run at different values of $r$ and $M$
(notice the log-scale); the solid lines are again obtained by eq. (34) but now
plotted as a function of $(r,\log_{2}M)$ and we see that, as the quality $|r|$
reaches zero, more and more examples are needed to ensure the retrieval of the
archetype. The colormap on the right is shared by the two panels.

###### Proof.

First we recall the unnormalized log-density of the model:

$-\beta\mathcal{H}^{\textrm{\tiny{HN},
sup}}(\boldsymbol{\sigma}|\boldsymbol{\eta})=\frac{\beta
N}{2\mathcal{R}}\sum_{\mu=1}^{K}\left(\frac{1}{NM}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu
a}\sigma_{i}\right)^{2}=\frac{\beta}{2\mathcal{R}NM^{2}}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N}\sum_{a,b=1}^{M}\eta_{i}^{\mu
a}\eta_{j}^{\mu
b}\sigma_{i}\sigma_{j}=\frac{1}{2}\sum_{k=1}^{N}h_{k}(\boldsymbol{\sigma})\sigma_{k}$
(35)

where

$h_{k}(\boldsymbol{\sigma}):=\frac{\beta}{\mathcal{R}NM^{2}}\sum_{\mu=1}^{K}\sum_{i\neq
k}^{N}\sum_{a,b=1}^{M}\eta_{i}^{\mu a}\eta_{k}^{\mu b}\sigma_{i}$ (36)

is the internal field acting on $\sigma_{k}$.
The discrete-time MC dynamics is given by the following update rule

$\sigma_{k}^{(n+1)}=\sigma_{k}^{(n)}\mathrm{sign}\left\\{\tanh\left[h_{k}^{(n)}(\boldsymbol{\sigma}^{(n)})\sigma_{k}^{(n)}\right]+\zeta_{k}^{(n)}\right\\},\quad\zeta_{k}^{(n)}\sim\mathcal{U}(-1,+1),$
(37)

where $n\in\mathbb{N}$ and $k=1,\dots,N$. By performing the low temperature
limit $\beta\to\infty$ equation (37) becomes

$\sigma_{k}^{(n+1)}=\sigma_{k}^{(n)}\mathrm{sign}\left(h_{k}^{(n)}(\boldsymbol{\sigma}^{(n)})\sigma_{k}^{(n)}\right).$
(38)

We want to study the stability of the archetype retrieval thus we set as an
initial configuration for the MC dynamics $\sigma_{k}^{(1)}=\xi_{k}^{1}$; the
one-step MC approximation for the magnetization is then

$m_{1}^{(2)}=\frac{1}{N}\sum_{k=1}^{N}\xi_{k}^{1}\sigma_{k}^{(2)}=\frac{1}{N}\sum_{k=1}^{N}\mathrm{sign}\left(h_{k}^{(1)}(\boldsymbol{\xi}^{1})\xi_{k}^{1}\right).$
(39)

If $N\gg 1$, by the central limit theorem, equation (39) can be approximated
as follows

$m_{1}^{(2)}\underset{N\gg
1}{\approx}\int\frac{dz}{\sqrt{2\pi}}\exp\left(-\frac{z^{2}}{2}\right)\mathrm{sign}\left(\mu_{1}+z\sqrt{\mu_{2}-\mu_{1}^{2}}\right)=\mathrm{erf}\left(\frac{\mu_{1}}{\sqrt{2(\mu_{2}-\mu_{1}^{2})}}\right)$
(40)

where we posed

$\displaystyle\mu_{1}$ $\displaystyle:=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}^{(n)}(\boldsymbol{\xi}^{1})\xi_{k}^{1}\right],$
(41) $\displaystyle\mu_{2}$ $\displaystyle:=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}^{(n)}(\boldsymbol{\xi}^{1})^{2}\right],$
(42)

with $\langle\cdot\rangle$ the average over the noise affecting the patterns.
First, we calculate $\mu_{1}$:

$\mu_{1}=\frac{\beta}{\mathcal{R}NM^{2}}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu=1}^{K}\sum_{i\neq
k}^{N}\sum_{a,b=1}^{M}\xi_{k}^{\mu}\chi_{k}^{\mu a}\xi_{i}^{\mu}\chi_{i}^{\mu
b}\xi_{i}^{1}\xi_{k}^{1}\right]=\frac{\beta}{\mathcal{R}NM^{2}}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu>1}^{K}\sum_{a,b=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu
b}+\sum_{a,b}^{M}\sum_{i\neq k}^{N}\chi_{k}^{1a}\chi_{i}^{1b}\right].$ (43)

Exploiting the fact that, by construction, archetypes are orthogonal in the
average, that is $\mathbb{E}_{\xi}[\xi^{\mu}_{i}\xi^{\mu}_{k}]=0$ for any
$i\neq k$,

$\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu>1}^{K}\sum_{a,b=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu b}\right]=0,$
(44)

and recalling (3)

$\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b}\sum_{i\neq
k}\chi_{k}^{1a}\chi_{i}^{1b}\right]=(N-1)M^{2}r^{2},$ (45)

we get

$\mu_{1}=\frac{\beta}{(1+\rho)N}(N-1)\underset{N\gg
1}{\approx}\mu_{1}=\frac{\beta}{(1+\rho)}.$ (46)

Then, we calculate $\mu_{2}$:

$\displaystyle\mu_{2}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}(\boldsymbol{\boldsymbol{\xi}}^{1})^{2}\right]=\frac{\beta^{2}}{\mathcal{R}^{2}N^{2}M^{4}}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\left(\sum_{\mu>1}^{K}\sum_{a,b=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu
b}+\sum_{a,b}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1b}\right)^{2}\right]=$ (47) $\displaystyle=$
$\displaystyle\frac{\beta^{2}}{\mathcal{R}^{2}N^{2}M^{4}}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\left(\sum_{\mu>1}^{K}\sum_{a,b=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu
b}\right)^{2}+\left(\sum_{a,b}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1b}\right)^{2}\right].$ (48)

For the sake of simplicity we pose
$A:=\left(\sum_{\mu>1}^{K}\sum_{a,b=1}^{M}\sum_{i\neq
k}^{N}\xi_{k}^{\mu}\xi_{i}^{\mu}\chi_{k}^{\mu a}\chi_{i}^{\mu b}\right)^{2}$
and $B:=\left(\sum_{a,b}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{1a}\chi_{i}^{1b}\right)^{2}$such that

$\mu_{2}=\frac{\beta^{2}}{\mathcal{R}^{2}N^{2}M^{4}}\left(\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[A\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[B\right]\right).$
(49)

After minimal manipulations we can write these terms as

$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[A\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\nu>1}\sum_{a,b,c,d}\sum_{i\neq
k}\chi_{k}^{\mu a}\chi_{i}^{\mu b}\chi_{k}^{\mu c}\chi_{i}^{\mu d}\right],$
(50) $\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[B\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b,c,d}^{M}\sum_{i\neq
k}^{N}\chi_{k}^{a1}\chi_{i}^{1b}\chi_{k}^{1c}\chi_{i}^{1d}\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b,c,d}^{M}\sum_{i,j\neq
k}^{N}\lambda_{ij}\chi_{k}^{1a}\chi_{i}^{1b}\chi_{k}^{1c}\chi_{j}^{1d}\right],$
(51)

where $\lambda_{ij}:=1-\delta_{ij}$ . By merging the results in (50) and (51)
we get

$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[A\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[B\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{\mu=1}\sum_{a,b,c,d,}\sum_{i\neq
k}\chi_{k}^{\mu a}\chi_{i}^{\mu b}\chi_{k}^{\mu c}\chi_{i}^{\mu
d}\right]+\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\sum_{a,b,c,d}^{M}\sum_{i,j\neq
k}^{N}\lambda_{ij}\chi_{k}^{1a}\chi_{i}^{1b}\chi_{k}^{1c}\chi_{j}^{1d}\right]=$
$\displaystyle=$ $\displaystyle
K(N-1)M^{4}r^{4}\left(1+\rho\right)^{2}+(N-1)(N-2)M^{4}r^{4}\left(1+\rho\right)$
$\displaystyle\underset{N\gg 1}{\approx}$ $\displaystyle
KNM^{4}r^{4}\left(1+\rho\right)^{2}+N^{2}M^{4}r^{4}\left(1+\rho\right).$ (53)

By direct substitution of (53) into (49) we obtain

$\mu_{2}\underset{N\gg
1}{\approx}\beta^{2}\left(\alpha+\frac{1}{1+\rho}\right).$ (54)

Finally, by plugging equations (46) and (54) into the expression of the
magnetization (40) we get

$m_{1}^{(2)}\underset{N\gg
1}{\approx}\mathrm{erf}\left(\frac{1}{\sqrt{2\alpha\left(1+\rho\right)^{2}+2\rho}}\right),$
(55)

and, by requiring that this one-step MC magnetization is larger than
$\textrm{erf}(\theta)$ we recover Eq. (34). ∎

######  Remark 2.

Recalling Eq. (40), one can see that setting $\theta=\frac{1}{\sqrt{2}}$ in
(34) corresponds to the condition
$\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[h_{k}^{(n)}(\boldsymbol{\xi}^{1})\xi_{k}^{1}\right]>\sqrt{\textrm{Var}(h_{k}^{(n)}(\boldsymbol{\xi}^{1})\xi_{k}^{1})}$,
that is the standard condition used in signal-to-noise analysis; this is also
depicted in Fig. S2

##  III Statistical mechanics of Supervised Hebbian Learning

The route to pave in this approach is the standard one in glassy statistical-
mechanics [4]: at first we have to outline the control parameters of the
system (i.e., $\alpha,\beta,\rho$), introduce the partition function
$\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S})$ and quenched free-
energy $f_{\alpha,\beta,\rho}$ related to the cost function (8), as well as
the order parameters necessary to describe the macroscopic behaviour of the
system. We must then achieve an explicit expression of the quenched free-
energy in terms of these parameters: this will be accomplished in the replica
symmetric (RS) regime (namely under the assumption that the order parameters
do not fluctuate in the thermodynamic limit) and, finally, we must extremize
the free energy w.r.t. the order parameters to force the thermodynamic
requirements of minimum energy and maximum entropy. This extremization returns
a set of self-consistent equations for the evolution of the order parameters
in the space of the control parameters whose inspection allows us to draw a
phase diagram for the network and thus to obtain an exhaustive
characterization of its emergent information-processing skills.
We will address this investigation in a mathematical rigorous way, by relying
on the generalized Guerra’s interpolation scheme [5, 6, 7], rather than
pseudo-heuristic tools like replica-trick: the basic idea of this approach is
to introduce a suitable functional that interpolates between the original
model and a series of simpler (and solvable) models and then to solve for
these simpler models and propagate the solution back to the other extremum of
the interpolation, namely the model under study.
The statistical mechanical treatment of the unsupervised HN has been already
addressed in [3] hence here we will just focus on the supervised case and we
can drop the superscript “sup” without ambiguity.

######  Definition 4.

The control parameters that tune the performance of the HN implementing a
supervised Hebbian learning (8) are

  * •

the network _storage_ $\alpha:=\lim_{N\to\infty}\frac{K}{N}$, that is, the
ratio between the number $K$ of archetypes that we aim to retrieve and the
number $N$ of neurons employed for this task,

  * •

the inverse _noise_ $\beta:=T^{-1}$ tuning the stochasticity in the network
dynamics,

  * •

the dataset _entropy_ 111Strictly speaking, $\rho$ is not an entropy, yet here
we allow ourselves for this slight abuse of language because, as discussed in
the main text, the conditional entropy
$H(\xi_{i}^{\mu}|\boldsymbol{\eta}_{i}^{\mu})$ is a monotonically increasing
function of $\rho$. $\rho:=\frac{1-r^{2}}{Mr^{2}}$, that is the amount of
information needed to describe the archetype set
$\\{\boldsymbol{\xi}^{\mu}\\}_{\mu=1,...,K}$ given the sample $\mathcal{S}$.

######  Definition 5.

The partition function of the supervised HN (8) is defined as

$\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S}):=\sum_{\\{\sigma\\}}^{2^{N}}\exp\left[-\beta\mathcal{H}^{\textrm{\tiny{HN}}}(\boldsymbol{\sigma}|\mathcal{S})\right]=\sum_{\\{\sigma\\}}^{2^{N}}\exp\left[\frac{\beta
N}{2\mathcal{R}}\sum_{\mu=1}^{K}\left(\frac{1}{NM}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu
a}\sigma_{i}\right)^{2}\right]$ (56)

and the related quenched free-energy is defined as

$f_{\alpha,\beta,\rho}:=-\frac{1}{\beta
N}\mathbb{E}_{\xi}\mathbb{E}_{\chi}\left[\log\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S})\right]$
(57)

where the average $\mathbb{E}$ is specified in Definition (2).

######  Definition 6.

The order parameters required to describe the performance of the HN
implementing a supervised Hebbian learning (8) are

  * •

the _archetype magnetization_
$m:=\frac{1}{N}\sum_{i=1}^{N}\xi_{i}^{1}\sigma_{i}$,

  * •

the _mean example magnetization_
$n:=\frac{r}{\mathcal{R}}\frac{1}{NM}\sum_{i,a=1}^{N,M}\eta_{i}^{1a}\sigma_{i}$,

  * •

the two-replica _overlap_ for the Ising neurons
$q_{ab}:=\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}^{(a)}\sigma_{i}^{(b)},$

  * •

the two-replica _overlap_ for the Gaussian neurons (vide infra)
$p_{ab}:=\frac{1}{K-1}\sum_{\mu=2}^{K}z_{\mu}^{(a)}z_{\mu}^{(b)}$.

Notice that $m$ and $n$ are referred to the first archetype without loss of
generality.
Further, being $x$ the generic order parameter, we denote with $\langle
x\rangle$ its expectation under the Boltzmann-Gibbs measure
$[\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S})]^{-1}\exp\left[-\beta\mathcal{H}^{\textrm{\tiny{HN}}}(\boldsymbol{\sigma}|\mathcal{S})\right]$.

######  Proposition 3.

The integral representation of the partition function of the HN implementing
supervised Hebbian learning (8) reads as

$\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\boldsymbol{\lambda},\boldsymbol{\eta}^{1})=\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\left[-\sum_{\mu=2}^{K}\frac{z_{\mu}^{2}}{2}+\sqrt{\frac{\beta}{N}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\lambda_{i}^{\mu}z_{\mu}\sigma_{i}+\frac{\beta
N}{2}{n}^{2}(1+\rho)\right],$ (58)

where $\boldsymbol{\lambda}=\\{\lambda_{i}^{\mu}\\}_{i=1,...,N}^{\mu=1,...,K}$
represent i.i.d standard-Gaussian random fields.

###### Proof.

First, we split the Hamiltonian
$\mathcal{H}^{\textrm{\tiny{HN}}}(\boldsymbol{\sigma}|\mathcal{S})$ into two
contributions: one containing the terms related to examples with label $\mu=1$
(playing as a “signal”) and the other containing the terms related to examples
with labels $\mu\neq 1$ (playing as slow noise). Then, exploiting the relation
$\int dz\exp(-z^{2}/2+Az)=\sqrt{2\pi}\exp(A^{2}/2)$, we linearise the noise
terms inside (56) obtaining

$\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S})=\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\left[-\sum_{\mu=2}^{K}\frac{z_{\mu}^{2}}{2}+\sqrt{\frac{\beta}{N\mathcal{R}}}\frac{1}{M}\sum_{\mu=2}^{K}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu
a}\sigma_{i}z_{\mu}+\frac{\beta
N}{2\mathcal{R}}\left(\frac{1}{NM}\sum_{i,a=1}^{N,M}\eta_{i}^{1a}\sigma_{i}\right)^{2}\right],$
(59)

and, recalling Definition (6), we rewrite the signal contribution, namely

$\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S})=\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\left[-\sum_{\mu=2}^{K}\frac{z_{\mu}^{2}}{2}+\sqrt{\frac{\beta}{N\mathcal{R}}}\frac{1}{M}\sum_{\mu=2}^{K}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu
a}\sigma_{i}z_{\mu}+\frac{\beta N}{2}(1+\rho){n}^{2}\right].$ (60)

We now handle the noise contribution trying to make it analytically more
treatable: from (60) we extract the noise contribution, namely

$\frac{1}{M}\sqrt{\frac{\beta}{N\mathcal{R}}}\sum_{\mu=2}^{K}\sum_{i,a=1}^{N,M}\eta_{i}^{\mu
a}z_{\mu}\sigma_{i}=\sqrt{\frac{\beta}{N\mathcal{R}}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\left(\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu
a}\right)z_{\mu}\sigma_{i}$ (61)

and we see that the random field acting on the pairs $\sigma_{i}z_{\mu}$ is
$\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu a}$, hence, in the large dataset
scenario $M\gg 1$, this random field can be replaced by a Gaussian-distributed
random field with the same mean and variance that turn out to be,
respectively,

$\mathbb{E}_{\chi}\mathbb{E}_{\xi}\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu
a}=\frac{1}{M}\sum_{a=1}^{M}\mathbb{E}_{\xi}\xi_{i}^{\mu}\,\mathbb{E}_{\chi}\chi_{i}^{\mu
a}=0$ (62)
$\begin{split}\mathbb{E}_{\chi}\mathbb{E}_{\xi}\left(\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu
a}\right)^{2}=&\mathbb{E}_{\chi}\mathbb{E}_{\xi}\left(\frac{1}{M}\sum_{a=1}^{M}\chi_{i}^{\mu
a}\xi_{i}^{\mu}\right)^{2}=\\\
=&\mathbb{E}_{\chi}\left(\frac{1}{M}\sum_{a=1}^{M}\chi_{i}^{\mu
a}\right)^{2}=\\\
=&\frac{1}{M^{2}}\mathbb{E}_{\chi}\left[\sum_{a=1}^{M}\left(\chi_{i}^{\mu
a}\right)^{2}+\sum_{a\neq b=1}^{M,M}\chi_{i}^{\mu a}\chi_{i}^{\mu
b}\right]=\\\
=&\frac{1}{M^{2}}\left[M+M(M-1)r^{2}\right]=r^{2}+\frac{1-r^{2}}{M}=\mathcal{R}\end{split}$
(63)

where $\mathbb{E}_{\xi}$ and $\mathbb{E}_{\chi}$ stand for averages over the
distribution of the random tensors $\xi$ and $\chi$ whose probability density
was specified in Definition (3).

Therefore, the noise term in the free energy can be replaced by

$\sqrt{\frac{\beta}{N}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\lambda_{i}^{\mu}z_{\mu}\sigma_{i}\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{with}\leavevmode\nobreak\ \leavevmode\nobreak\
\lambda_{i}^{\mu}\sim\mathcal{N}(0,1)$ (64)

recovering eq. (58). ∎

######  Remark 3.

The partition function $\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta}(\mathcal{S})$
given in (60) can be seen as the partition function of a two-species spin-
glass model, where one species is made of binary spins
$\boldsymbol{\sigma}\in\\{-1,+1\\}^{N}$ and the other species is made of
Gaussian spins $\boldsymbol{z}\sim\mathcal{N}(0,1)^{K-1}$; spins of different
nature interact pairwisely by a random coupling given by
$\sqrt{\frac{\beta}{N}}\lambda_{i}^{\mu}$, binary spins interact pairwisely by
a Hebbian-like coupling
$\frac{\beta}{2\mathcal{R}}\bar{\eta}_{i}^{1}\bar{\eta}_{j}^{1}$; interactions
between real spins are absent.

######  Remark 4.

The previous proposition is consistent with the universality of the quenched
noise proved for the Sherrington-Kirkpatrick spin-glasses [8] and extended to
bipartite spin-glasses and neural networks [9, 10]. Remarkably, it allows us
to substitute the digital i.i.d. entries that contribute to the slow noise
with Gaussian random variables and this, in turn, allows us to apply the Wick-
Isserlis theorem to analytically treat the expression of the quenched free-
energy.

As a result of the previous proposition we can recast the quenched free-energy
(57) as

$f_{\alpha,\beta,\rho}:=-\frac{1}{\beta
N}\mathbb{E}\log\left\\{\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\left[-\sum_{\mu=2}^{K}\frac{z_{\mu}^{2}}{2}+\sqrt{\frac{\beta}{N}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\lambda_{i}^{\mu}z_{\mu}\sigma_{i}+\frac{\beta
N}{2}{n}^{2}(1+\rho)\right]\right\\}$ (65)

where the operator
$\mathbb{E}:=\mathbb{E}_{\chi^{1}}\mathbb{E}_{\xi^{1}}\mathbb{E}_{\lambda}$ is
defined by

$\displaystyle\mathbb{E}_{\lambda}f(\boldsymbol{\lambda})$ $\displaystyle:=$
$\displaystyle\int\prod_{a=1}^{M}\prod_{i=1}^{N}\left\\{\frac{d\lambda_{i}^{\mu}}{\sqrt{2\pi}}\exp\left[-\frac{\left(\lambda_{i}^{\mu}\right)^{2}}{2}\right]\right\\}f(\boldsymbol{\lambda}),$
(66)

where $f$ is a generic function and by Definition (2).

We are now ready to introduce the Guerra functional that will guide us toward
the solution of the model. In a nutshell, the idea is that the Guerra
functional interpolates between two extrema: one corresponds to the original
model (and the functional coincides with the model free-energy), the other
corresponds to a solvable model (and the functional coincides with the free-
energy of a one-body model); therefore, one can solve for the Guerra
functional in this second extremum and then propagate the solution back to the
original model via the fundamental theorem of calculus. We start this journey
by giving the following

######  Definition 7.

The Guerra interpolating functional $G_{\alpha,\beta,\rho}(t,J_{m})$ is
defined as:

$\begin{split}G_{\alpha,\beta,\rho}(t,J_{m})=-\frac{1}{\beta
N}\mathbb{E}\log\Big{\\{}\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\Big{[}-\frac{1-(1-t)\beta(1-\langle
q\rangle)}{2}\sum_{\mu=2}^{K}z_{\mu}^{2}+\sqrt{t}\sqrt{\frac{\beta}{N}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\lambda_{i}^{\mu}z_{\mu}\sigma_{i}+\\\
+t\frac{\beta N}{2}{n}^{2}(1+\rho)+\beta\langle
n\rangle(1+\rho)\left(1-t\right)N{n}+\sqrt{\alpha\beta\langle
p\rangle\left(1-t\right)}\sum_{i=1}^{N}\theta_{i}\sigma_{i}+\sqrt{\beta\langle
q\rangle\left(1-t\right)}\sum_{\mu=2}^{K}\psi_{\mu}z_{\mu}-J_{m}\beta
Nm\Big{]}\Big{\\}}\end{split}$ (67)

where $t\in[0,1]$ is the interpolation parameter, $J_{m}$ is an auxiliary
field coupled to the Mattis magnetization of the archetype, the operator
$\mathbb{E}$ is defined as

$\mathbb{E}:=\mathbb{E}_{\chi^{1}}\mathbb{E}_{\xi^{1}}\mathbb{E}_{\lambda}\mathbb{E}_{\theta}\mathbb{E}_{\psi}$
(68)

with

$\displaystyle\mathbb{E_{\theta}}f(\boldsymbol{\theta})$ $\displaystyle:=$
$\displaystyle\int\prod_{\mu=2}^{K}\left\\{\frac{d\theta_{\mu}}{\sqrt{2\pi}}\exp\left[-\theta_{\mu}^{2}/2\right]\right\\}f(\boldsymbol{\theta}),$
(69) $\displaystyle\mathbb{E_{\psi}}f(\boldsymbol{\psi})$ $\displaystyle:=$
$\displaystyle\int\prod_{i=1}^{N}\left\\{\frac{d\psi_{i}}{\sqrt{2\pi}}\exp\left[-\psi_{i}^{2}/2\right]\right\\}f(\boldsymbol{\psi}),$
(70)

and the brackets $\langle\cdot\rangle$ represent the Boltzmann-Gibbs measure
(as further specified in Definition (9)).

######  Remark 5.

By setting $t=1$ and $J_{m}=0$ in Guerra’s functional, we recover immediately
the correct expression for the free energy of the model (65):
$G_{\alpha,\beta,\rho}(t=1,J_{m}=0)=f_{\alpha,\beta,\rho}$. On the other hand,
by setting $t=0$ we end up with a one-body system that is exactly solvable.
The fictitious field yielded by $J_{m}$ will allow us to determine the
expectation of the magnetization by deriving the free energy, namely, $\langle
m\rangle=\left.\frac{dG_{\alpha,\beta,\rho}(t=1,J_{m})}{dJ_{m}}\right|_{J_{m}=0}$.

######  Definition 8.

The functional (67) yields an interpolating Boltzmann-Gibbs measure under
which the expectation of the generic function
$f(\boldsymbol{\sigma},\boldsymbol{z})$ reads as

$\left\langle
f(\boldsymbol{\sigma},\boldsymbol{z})\right\rangle_{t}:=\frac{\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)f(\boldsymbol{\sigma},\boldsymbol{z})\cdot\exp\left[\mathcal{B}(\boldsymbol{\sigma},\boldsymbol{z})\right]}{\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\left[\mathcal{B}(\boldsymbol{\sigma},\boldsymbol{z})\right]}$
(71)

where $\mathcal{B}(\boldsymbol{\sigma},\boldsymbol{z})$ is the generalized
Boltzmann-Gibbs weight given by

$\displaystyle\mathcal{B}(\boldsymbol{\sigma},\boldsymbol{z}):=$
$\displaystyle-\frac{1-(1-t)\beta(1-\langle
q\rangle)}{2}\sum_{\mu=2}^{K}z_{\mu}^{2}+\sqrt{t}\sqrt{\frac{\beta}{N}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\lambda_{i}^{\mu}z_{\mu}\sigma_{i}+$
(72) $\displaystyle+t\frac{\beta N}{2}{n}^{2}(1+\rho)+\beta\langle
n\rangle(1+\rho)\left(1-t\right)N{n}+\sqrt{\alpha\beta\langle
p\rangle\left(1-t\right)}\sum_{i=1}^{N}\theta_{i}\sigma_{i}+\sqrt{\beta\langle
q\rangle\left(1-t\right)}\sum_{\mu=2}^{K}\psi_{\mu}z_{\mu}-J_{m}\beta Nm.$

######  Definition 9.

Theoretical setting (Replica Symmetric Ansatz): In the replica-symmetric
framework we assume that, in the thermodynamic limit $N\to\infty$, the
following variance and covariance

$\displaystyle\Delta\left[n_{1}^{2}\right]$ $\displaystyle:=$
$\displaystyle\mathbb{E}\left\langle\left({n}-\langle
n\rangle\right)^{2}\right\rangle_{t}=\mathbb{E}\left\langle{n}^{2}\right\rangle_{t}+2\langle
n\rangle\mathbb{E}\left\langle n\right\rangle_{t}+\langle n\rangle^{2}$ (73)
$\displaystyle\Delta\left[q_{12}p_{12}\right]$ $\displaystyle:=$
$\displaystyle\mathbb{E}\left\langle\left(q_{12}-\langle
q\rangle\right)\left(p_{12}-\langle
p\rangle\right)\right\rangle_{t}=\mathbb{E}\left\langle
q_{12}p_{12}\right\rangle_{t}-\langle q\rangle\mathbb{E}\left\langle
p_{12}\right\rangle_{t}-\langle p\rangle\mathbb{E}\left\langle
q_{12}\right\rangle_{t}+\langle p\rangle\langle q\rangle$ (74)

go to $0$, namely $\lim_{N\to\infty}\Delta\left[n_{1}^{2}\right]=0$ and
$\lim_{N\to\infty}\Delta\left[q_{12}p_{12}\right]=0$. Which implies

$\displaystyle\lim_{N\to\infty}\langle n\rangle_{t}$ $\displaystyle=$
$\displaystyle\langle n\rangle,$ (75) $\displaystyle\lim_{N\to\infty}\langle
q_{12}\rangle_{t}$ $\displaystyle=$ $\displaystyle\langle q\rangle,$ (76)
$\displaystyle\lim_{N\to\infty}\langle p_{12}\rangle_{t}$ $\displaystyle=$
$\displaystyle\langle p\rangle.$ (77)

As we want to calculate the free energy (i.e., Guerra’s functional at $t=1$)
by evaluating $G_{\alpha,\beta,\rho}(t=0,J_{m})$ and then propagating it to
$t=1$, a technical aspect we need to preliminary address is the
$t-$derivative, or the streaming, of Guerra’s interpolating functional and
this is achieved in the next

######  Lemma 1.

The streaming equation for Guerra’s functional
$G_{\alpha,\beta,\rho}(t,J_{m})$ is

$\frac{d}{dt}G_{\alpha,\beta,\rho}(t,J_{m})=\frac{\langle
n\rangle^{2}}{2}(1+\rho)+\frac{\alpha}{2}\langle p\rangle(1-\langle
q\rangle)+\Delta\left[q_{12}p_{12}\right]-\frac{1}{2}\Delta\left[{n}^{2}\right].$
(78)

###### Proof.

By a direct evaluation

$\displaystyle\frac{d}{dt}G_{\alpha,\beta,\rho}(t,J_{m})$ $\displaystyle=$
$\displaystyle\frac{1}{\beta N}\mathbb{E}\Biggr{(}\frac{\beta\left(1-\langle
q\rangle\right)}{2}\sum_{\mu=2}^{K}\langle
z_{\mu}^{2}\rangle_{t}-\frac{1}{2\sqrt{t}}\sqrt{\frac{\beta}{N}}\sum_{\mu=2}^{K}\sum_{i=1}^{N}\lambda_{i}^{\mu}\langle
z_{\mu}\sigma_{i}\rangle_{t}+$ $\displaystyle-$ $\displaystyle\frac{\beta
N(1+\rho)}{2}\langle{n}^{2}\rangle_{t}+\beta\langle
n\rangle(1+\rho)N\Big{\langle}n\Big{\rangle}_{t}+\frac{\sqrt{\beta\alpha\langle
p\rangle}}{2\sqrt{1-t}}\sum_{i=1}^{N}\theta_{i}\langle\sigma_{i}\rangle_{t}+\frac{\sqrt{\beta\langle
q\rangle}}{2\sqrt{1-t}}\sum_{\mu=2}^{K}\psi_{\mu}\langle
z_{\mu}\rangle_{t}\Biggr{)},$

Now we can introduce the order parameters in the above expression, obtaining

$\displaystyle\frac{d}{dt}G_{\alpha,\beta,\rho}(t,J_{m})$ $\displaystyle=$
$\displaystyle\mathbb{E}\Biggr{(}\frac{\beta\alpha\left(1-\langle
q\rangle\right)}{2}\left\langle
p_{11}\right\rangle_{t}-\frac{\alpha\beta}{2}\left(\left\langle
p_{11}\right\rangle_{t}-\left\langle p_{12}q_{12}\right\rangle_{t}\right)+$
$\displaystyle-$
$\displaystyle\frac{1+\rho}{2}\left\langle{n}^{2}\right\rangle_{t}+(1+\rho)\langle
n\rangle\left\langle n\right\rangle_{t}+\frac{\beta\alpha\langle
p\rangle}{2}\left(1-\left\langle
q_{12}\right\rangle_{t}\right)+\frac{\beta\alpha\langle
q\rangle}{2}\sum_{\mu=2}^{K}\left(\left\langle
p_{11}\right\rangle_{t}-\left\langle p_{12}\right\rangle_{t}\right)\Biggr{)}.$

By expressing in terms of mean values, variances and covariances all the
correlation functions contained in the above expression and forcing all the
$\Delta$’s to go to $0$ in the thermodynamic limit (i.e., by assuming the RS
ansatz), we obtain (78). ∎

We are now ready to state the main

######  Theorem 1.

The replica symmetric free-energy of the HN implementing the supervised
Hebbian learning (8), in the infinite volume limit $N\to\infty$ and large
dataset scenario $M\gg 1$, can be expressed in terms of control and order
parameters as

$\displaystyle f_{\alpha,\beta,\rho}=$
$\displaystyle\frac{\alpha}{2\beta}\log\left[1-\beta\left(1-\langle
q\rangle\right)\right]-\frac{\alpha}{2}\frac{\langle
q\rangle}{1-\beta\left(1-\langle q\rangle\right)}+\frac{\langle
n\rangle^{2}}{2}(1+\rho)+\frac{\alpha}{2}\langle p\rangle(1-\langle q\rangle)$
(81) $\displaystyle-\frac{1}{\beta}\log
2-\frac{1}{\beta}\mathbb{E}_{\chi}\mathbb{E}_{\theta}\log\cosh\left.\Big{(}\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}-J_{m}\beta+\sqrt{\alpha\beta\langle
p\rangle}\theta\Big{)}\right|_{J_{m}=0}.$

###### Proof.

Via the fundamental theorem of calculus we write

$f_{\alpha,\beta,\rho}(J_{m})=G_{\alpha,\beta,\rho}(t=1,J_{m})=G_{\alpha,\beta,\rho}(t=0,J_{m})+\int_{0}^{1}dt^{\prime}\frac{dG_{\alpha,\beta,\rho}}{dt}(t^{\prime},J_{m})$
(82)

so that we are left with the evaluation of two terms, the Cauchy datum
$G_{\alpha,\beta,\rho}(t=0,J_{m})$, that can be faced by direct calculation,
and the integral of the streaming of Guerra’s functional, that was provided in
Lemma 1. It will be in this last computation that the RS assumption will allow
us to obtain an explicit expression for the free energy in terms of the order
parameters and thus solve for the model, as shown hereafter.
We start by evaluating $G_{\alpha,\beta,\rho}(t=0,J_{m})$:

$\displaystyle G_{\alpha,\beta,\rho}(t=0,J_{m})$ $\displaystyle=$
$\displaystyle-\frac{1}{\beta
N}\mathbb{E}\log\Big{\\{}\sum_{\\{\sigma\\}}\int\prod_{\mu=2}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\Big{[}-\frac{1-\beta\left(1-\langle
q\rangle\right)}{2}\sum_{\mu=2}^{K}z_{\mu}^{2}+\sqrt{\beta\langle
q\rangle}\sum_{\mu=2}^{K}\psi_{\mu}z_{\mu}$ (83) $\displaystyle+$
$\displaystyle\frac{\beta\langle
n\rangle}{rM}\sum_{i,a=1}^{N,M}\xi_{i}^{1}\chi_{i}^{1a}\sigma_{i}+\sqrt{\alpha\beta\langle
p\rangle}\sum_{i=1}^{N}\theta_{i}\sigma_{i}-J_{m}\beta\sum_{i=1}^{N}\xi_{i}^{1}\sigma_{i}\Big{]}\Big{\\}},$

as this is a trivial one-body problem (i.e., its probability structure is
completely factorized), we immediately can write

$\displaystyle G_{\alpha,\beta,\rho}(t=0,J_{m})$ $\displaystyle=$
$\displaystyle-\frac{1}{\beta
N}\sum_{\mu=2}^{K}\mathbb{E}\log\Big{\\{}\int\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)\exp\Big{[}-\frac{1-\beta\left(1-\langle
q\rangle\right)}{2}z_{\mu}^{2}+\sqrt{\beta\langle
q\rangle}\psi_{\mu}z_{\mu}\Big{]}\Big{\\}}$ (84) $\displaystyle-$
$\displaystyle\frac{1}{\beta
N}\sum_{i=1}^{N}\mathbb{E}\log\Big{\\{}\sum_{\sigma_{i}=\pm
1}\exp\Big{[}\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\xi_{i}^{1}\chi_{i}^{1a}\sigma_{i}+\sqrt{\alpha\beta\langle
p\rangle}\theta_{i}\sigma_{i}-J_{m}\beta\xi_{i}^{1}\sigma_{i}\Big{]}\Big{\\}},$

hence, by direct evaluation of the Gaussian integral in $z$ and sum in
$\sigma$, we reach

$\displaystyle G_{\alpha,\beta,\rho}(t=0,J_{m})$ $\displaystyle=$
$\displaystyle-\frac{\alpha}{\beta}\mathbb{E}_{\psi}\log\Big{\\{}\frac{1}{\sqrt{1-\beta\left(1-\langle
q\rangle\right)}}\exp\left[\frac{\beta\langle
q\rangle\psi^{2}}{2\left[1-\beta\left(1-\langle
q\rangle\right)\right]}\right]\Big{\\}}$ (85) $\displaystyle-$
$\displaystyle\frac{1}{\beta}\log
2-\frac{1}{\beta}\mathbb{E}_{\chi}\mathbb{E}_{\theta}\log\Big{\\{}\cosh\Big{[}\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}-J_{m}\beta+\sqrt{\alpha\beta\langle
p\rangle}\theta\Big{]}\Big{\\}}$

where the first average can be carried out directly returning

$\displaystyle G_{\alpha,\beta,\rho}(t=0,J_{m})$ $\displaystyle=$
$\displaystyle\frac{\alpha}{2\beta}\log\left[1-\beta\left(1-\langle
q\rangle\right)\right]-\frac{\alpha}{2\beta}\frac{\beta\langle
q\rangle}{1-\beta\left(1-\langle q\rangle\right)}$ (86) $\displaystyle-$
$\displaystyle\frac{1}{\beta}\log
2-\frac{1}{\beta}\mathbb{E}_{\chi}\mathbb{E}_{\theta}\log\cosh\left(\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}-J_{m}\beta+\sqrt{\alpha\beta\langle
p\rangle}\theta\right).$

Now, by using Eq. (78) into the initial scheme (82) we reach

$\displaystyle f_{\alpha,\beta,\rho}(J_{m})=$
$\displaystyle\frac{\alpha}{2\beta}\log\left[1-\beta\left(1-\langle
q\rangle\right)\right]-\frac{\alpha}{2\beta}\frac{\beta\langle
q\rangle}{1-\beta\left(1-\langle q\rangle\right)}+\frac{\langle
n\rangle^{2}}{2}(1+\rho)+\frac{\alpha}{2}\langle p\rangle(1-\langle
q\rangle)+\Delta\left[q_{12}p_{12}\right]-\frac{1}{2}\Delta\left[n_{1}^{2}\right]$
(87) $\displaystyle-\frac{1}{\beta}\log
2-\frac{1}{\beta}\mathbb{E}_{\chi}\mathbb{E}_{\theta}\log\cosh\left(\beta\langle
n\rangle\frac{1}{M}\sum_{a=1}^{M}\chi^{a}-J_{m}\beta+\sqrt{\alpha\beta\langle
p\rangle}\theta\right),$

such that, finally, sending $N\to\infty$, under the RS assumption, at
$J_{m}=0$ we obtain Eq. (81). ∎

######  Corollary 1.

In the infinite volume limit $N\to\infty$ and large dataset scenario $M\gg 1$,
the replica-symmetric self-consistent equations for the evolution of the order
parameters of the HN (8) implementing the supervised Hebbian learning in the
space of the tunable parameters are

$\displaystyle\langle p\rangle$ $\displaystyle=\frac{\beta\langle
q\rangle}{\left[1-\beta\left(1-\langle q\rangle\right)\right]^{2}},$ (88)
$\displaystyle\langle n\rangle$
$\displaystyle=\frac{1}{1+\rho}\mathbb{E}_{\chi}\mathbb{E}_{\theta}\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}\tanh\left(\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}+\sqrt{\alpha\beta\langle
p\rangle}\theta\right),$ $\displaystyle\langle q\rangle$
$\displaystyle=\mathbb{E}_{\chi}\mathbb{E}_{\theta}\tanh^{2}\left(\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}+\sqrt{\alpha\beta\langle
p\rangle}\theta\right),$ $\displaystyle\langle m\rangle$
$\displaystyle=\mathbb{E}_{\chi}\mathbb{E}_{\theta}\tanh\left(\beta\langle
n\rangle\frac{1}{Mr}\sum_{a=1}^{M}\chi^{a}+\sqrt{\alpha\beta\langle
p\rangle}\theta\right).$

###### Proof.

Let us consider the explicit expression for the quenched free-energy given in
(81) and extremize w.r.t. the order parameters:

$\frac{\partial f_{\alpha,\beta,\rho}(J_{m})}{\partial\langle
n_{1}\rangle}=\frac{\partial f_{\alpha,\beta,\rho}(J_{m})}{\partial\langle
q\rangle}=\frac{\partial f_{\alpha,\beta,\rho}(J_{m})}{\partial\langle
p\rangle}=0,$ (89)

we also pose, by construction,

$\frac{\partial f_{\alpha,\beta,\rho}(J_{m})}{\partial J_{m}}=\langle
m\rangle.$ (90)

With some algebra and by sending $J_{m}\to 0$, we obtain eqs. (88). ∎

######  Proposition 4.

In the large dataset limit $M\to\infty$, the replica-symmetric self-consistent
equations (88) for the HN (8) implementing the supervised Hebbian learning can
be expressed as

$\displaystyle\langle n\rangle$ $\displaystyle=\frac{1}{1+\rho}\frac{\langle
m\rangle}{1-\frac{\rho}{1+\rho}\beta(1-\langle q\rangle)},$ (91)
$\displaystyle\langle q\rangle$
$\displaystyle=\mathbb{E}_{z}\tanh^{2}\left[\beta\langle
n\rangle+z\beta\sqrt{\langle n\rangle^{2}\rho+\frac{\alpha\langle
q\rangle}{\left[1-\beta\left(1-\langle q\rangle\right)\right]^{2}}}\right],$
$\displaystyle\langle m\rangle$
$\displaystyle=\mathbb{E}_{z}\tanh\left[\beta\langle
n\rangle+z\beta\sqrt{\langle n\rangle^{2}\rho+\frac{\alpha\langle
q\rangle}{\left[1-\beta\left(1-\langle q\rangle\right)\right]^{2}}}\right].$

###### Proof.

Note that, for $M\to\infty$, we can write

$\frac{1}{M}\sum_{a=1}^{M}\chi^{a}\sim
r+\lambda\sqrt{\frac{1-r^{2}}{M}}\leavevmode\nobreak\
\mathrm{with}\leavevmode\nobreak\ \lambda\sim\mathcal{N}(0,1).$ (92)

At first let us tackle $\langle n_{1}\rangle$ as follows

$\displaystyle(1+\rho)\langle n\rangle$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\lambda}\mathbb{E}_{\theta}\left(1+\sqrt{\rho}\lambda\right)\tanh\left(\beta\langle
n\rangle\left(1+\sqrt{\rho}\lambda\right)+\sqrt{\alpha\beta\langle
p\rangle}\theta\right)=$ (93) $\displaystyle=$ $\displaystyle\langle
m\rangle+\sqrt{\rho}\mathbb{E}_{\lambda}\mathbb{E}_{\theta}\partial_{\lambda}\tanh\left[\beta\langle
n\rangle\left(1+\sqrt{\rho}\lambda\right)+\sqrt{\alpha\beta\langle
p\rangle}\theta\right]=$ $\displaystyle=$ $\displaystyle\langle
m\rangle+\beta\langle n\rangle\rho(1-\langle q\rangle)$

then we move on to $\langle m\rangle$

$\displaystyle\langle m\rangle$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\lambda}\mathbb{E}_{\theta}\tanh\left[\beta\langle
n\rangle\left(1+\sqrt{\rho}\lambda\right)+\sqrt{\alpha\beta\langle
p\rangle}\theta\right]=$ (94) $\displaystyle=$
$\displaystyle\mathbb{E}_{\lambda}\mathbb{E}_{\theta}\tanh\left[\beta\langle
n\rangle+\sqrt{\left(\beta\langle
n\rangle\right)^{2}\rho}\lambda+\sqrt{\alpha\beta\langle
p\rangle}\theta\right]=$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\tanh\left[\beta\langle
n\rangle+z\sqrt{\left(\beta\langle n\rangle\right)^{2}\rho+\alpha\beta\langle
p\rangle}\right],$

and, analogously for $\langle q\rangle$, we obtain (91). ∎

As a last manipulation we work out explicitly the self-consistencies in the
zero fast-noise limit $\beta\to\infty$ as this allows a characterization of
the ground state of the network too; results are provided in the next

######  Proposition 5.

The zero-temperature $\beta\to\infty$, large dataset $M\to\infty$ and large
network $N\to\infty$ limit of the self-consistencies for the order parameters
of the Hopfield model (8) implementing the supervised Hebbian learning read as

$\displaystyle\langle m\rangle$ $\displaystyle=\langle
n\rangle\left[1+\rho(1-\Delta)\right],$ (95) $\displaystyle G$
$\displaystyle=\sqrt{2\langle
n\rangle^{2}\rho+2\frac{\alpha}{\left(1-\Delta\right)^{2}}},$
$\displaystyle\Delta$
$\displaystyle=\frac{1}{G}\partial\mathrm{erf}\left(\frac{\langle
n\rangle}{G}\right),$ $\displaystyle\langle n\rangle$
$\displaystyle=\frac{1}{1+\rho(1-\Delta)}\mathrm{erf}\left(\frac{\langle
n\rangle}{G}\right),$

where

$\displaystyle\mathrm{erf}(x):=\int_{0}^{x}\frac{2}{\sqrt{\pi}}\exp(-t^{2})\,dt,$
(96) $\displaystyle\mathrm{\partial
erf}(x):=\frac{2}{\sqrt{\pi}}\exp(-x^{2}).$

###### Proof.

Following the strategy already used by AGS (e.g., see [11]) we pose
$\Delta:=\beta(1-\langle q\rangle)$ and we use this definition in Eq. (91)
obtaining

$\displaystyle\langle n\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{1+\rho}\frac{\langle
m\rangle}{1-\frac{\rho}{1+\rho}\Delta},$ (97) $\displaystyle\Delta$
$\displaystyle=$
$\displaystyle\beta\left\\{1-\mathbb{E}_{z}\tanh^{2}\left[\beta\langle
n\rangle+z\beta\sqrt{\langle
n\rangle^{2}\rho+\frac{\alpha\left(1-\frac{\Delta}{\beta}\right)}{\left[1-\Delta\right]^{2}}}\right]\right\\},$
(98) $\displaystyle\langle m\rangle$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\tanh\left[\beta\langle
n\rangle+z\beta\sqrt{\langle
n\rangle^{2}\rho+\frac{\alpha\left(1-\frac{\Delta}{\beta}\right)}{\left[1-\Delta\right]^{2}}}\right].$
(99)

Then, we add a field $\tilde{\beta}x$ inside the hyperbolic tangents, as this
allows to rewrite $\Delta$ as $\frac{\partial m(x)}{\partial x}$ and thus

$\displaystyle\langle n\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{1+\rho}\frac{\langle
m\rangle}{1-\frac{\rho}{1+\rho}\Delta},$ (100) $\displaystyle\Delta$
$\displaystyle=$ $\displaystyle\frac{\partial m(x)}{\partial x},$ (101)
$\displaystyle\langle m\rangle(x)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{z}\tanh\left[\beta x+\beta\langle
n\rangle+z\beta\sqrt{\langle
n\rangle^{2}\rho+\frac{\alpha\left(1-\frac{\Delta}{\beta}\right)}{\left[1-\Delta\right]^{2}}}\right],$
(102)

such that by taking the limit $\tilde{\beta}\to\infty$ and afterwards $x\to 0$
we get (95). ∎

The numerical solution of the self-consistent equations (88) were used to draw
the phase diagram presented in the main text, while here we explicitly show
the numerical solution for the magnetization in Fig. S3 $a$ and results
obtained by MC simulations in Fig. S3 $b$-$c$, and discuss their consistency.
In particular, in Fig. S3 $b$-$c$ we compare the outcomes for the estimates of
the magnetization $\langle m\rangle$ and its derivative $\partial\langle
m\rangle/\partial\rho$ w.r.t. $\rho$ that plays as a susceptibility of system
performance versus the dataset entropy. Notice that $\langle m\rangle$
exhibits a flex at a point $\rho_{c}$ that corresponds to the peak in the
susceptibility and this point gets closer and closer to the transition point
derived analytically. Further, by increasing the network size $N$, peaks in
the susceptibility get sharper. Therefore, these sizes can nicely spot the
existence and the location of the phase transition, on the other hand, the
first-order nature of the transition is well evidenced by the numerical
solution of self-consistencies, as shown in Fig. S3 $a$.

Figure S3: Theoretical solution and finite-size scaling by MC simulations.
Panel $a$: numerical solution of the self-consistent equation (91) for the
archetype magnetization versus $\sqrt{\rho}$. Notice that, as expected for
first-order transition, the order parameter exhibits a first-order singularity
and this happens at a value $\rho_{c}(\alpha,\beta)$ that is the one
highlighted by vertical lines in panels $b$ and $c$. Panel $b$: expected
magnetization for the archetype versus the entropy of the dataset resulting
from MC simulations. Panel $c$: Susceptibility w.r.t. $\rho$ versus the
dataset entropy resulting from MC simulations. For different choices of
$\alpha$, corresponding to curves of different colors as explained by the
legend, we varied the network size $N$ and the network load $K$ from $N\times
K=3000$ (solid line) to $N\times K=14000$ (dotted line) and to $N\times
K=43000$ (dashed line); vertical lines represent the transition points
$\rho_{c}(\alpha,\beta)$ as predicted theoretically. Note that, in panel $b$,
the curves exhibit a flex approximately corresponding to the intersection with
the vertical lines and in panel $c$ the curves peak approximately at the
intersection with the vertical lines, and these matches get sharper as the
networks size is increased. Simulations and numerical solutions are obtained
for $\beta=10$.  Figure S4: Schematic representation of the models
investigated in this work. In the upper part of this figure we show the HN
implementing the supervised Hebbian learning discussed in Secs. II and III.
Notice that the coupling between each neuron pair can be written as the scalar
product between the empirical averages of the blocks
$\boldsymbol{\eta}_{i}^{\mu}=(\eta_{i}^{\mu 1},...,\eta_{i}^{\mu M})$. In the
lower part of the figure we show a RBM which is trained according to the
grandmother-cell protocol: while the examples pertaining to the first
archetype are presented and the visible neurons are clamped as
$\boldsymbol{\sigma}=\boldsymbol{\eta}^{1a}$ for $a=1,...,M$, the first hidden
neuron is clamped as active, while the remaining hidden neurons are set
quiescent. The same operation is repeated for the other group pf examples.

##  IV Boltzmann machine Training

In this section we briefly address machine learning by one of its classical
architecture that is the RBM – in a statistical-mechanics jargon these are
nothing but bipartite spin-glasses [12] – in order to derive its supervised
learning rules and compare their learning skills with those shown by the HN
described in the previous sections, see Fig. S4. For the sake of simplicity
here, as standard in machine learning investigations, we fix $\beta=1$, with
no loss of generality as the fast noise can always be re-introduced in the
network by rescaling the couplings.

######  Definition 10.

We consider a RBM built of two-layers made of, respectively, $N$ binary
neurons $\sigma_{i}$, $i\in(1,...,N)$ and $K$ real-valued neurons $z_{\mu}$,
$\mu\in(1,...,K)$ equipped with a standard Gaussian prior and whose
Hamiltonian reads as

$\mathcal{H}^{\textrm{\tiny{RBM}}}(\boldsymbol{\sigma},\boldsymbol{z}|\boldsymbol{W})=-\frac{1}{\sqrt{N}}\sum_{i,\mu=1}^{N,K}W_{i\mu}\sigma_{i}z_{\mu}$
(103)

where $\boldsymbol{W}\in\mathbb{R}^{N\times K}$ is the matrix of weights among
the two layers (and there are no interactions within the same layer whence the
restriction).

######  Definition 11.

The partition function of the RBM (103) is

$\mathcal{Z}^{\textrm{\tiny{RBM}}}(\boldsymbol{W})=\sum_{\\{\sigma\\}}\int\prod_{\mu=1}^{K}dz_{\mu}\exp\left(\sum_{i,\mu=1}^{N,K}W_{i\mu}\sigma_{i}z_{\mu}\right)e^{-\left(\frac{1}{2}\sum_{\mu=1}^{K}z_{\mu}^{2}\right)},$
(104)

where the factor $\propto
e^{-\left(\frac{1}{2}\sum_{\mu=1}^{K}z_{\mu}^{2}\right)}$ is the Gaussian
prior for the $z$ neurons. The joint probability density related to the
partition function (104) assumes the following form

$\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z}|\boldsymbol{W})=\frac{1}{\mathcal{Z}^{\textrm{\tiny{RBM}}}(\boldsymbol{W})}\exp\left(\sum_{i,\mu=1}^{N,K}W_{i\mu}\sigma_{i}z_{\mu}-\sum_{\mu=1}^{K}\frac{z_{\mu}^{2}}{2}\right).$
(105)

######  Remark 6.

As a learning criterion we minimize the Kullback-Leibler cross-entropy
$D_{KL}(\mathcal{Q}\|\mathcal{P})$, defined as

$D_{KL}(\mathcal{Q}\|\mathcal{P})=-\sum_{\\{\sigma,z\\}}\mathcal{Q}(\sigma,z)\log\frac{\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z}|\boldsymbol{W})}{\mathcal{Q}(\boldsymbol{\sigma},\boldsymbol{z})}$
(106)

where $\mathcal{Q}(\boldsymbol{\sigma},\boldsymbol{z})$ is the empirical
probability density of the experimental dataset $\mathcal{S}$. The
minimization takes place by properly tuning the weights $\boldsymbol{W}$ in
order for $\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z}|\boldsymbol{W})$ to
be as “close” as possible to
$\mathcal{Q}(\boldsymbol{\sigma},\boldsymbol{z})$. This can be accomplished by
the gradient descent method that yield the following iterative rule for
$\boldsymbol{W}$

$W_{i\mu}^{(n+1)}=W_{i\mu}^{(n)}-\epsilon\frac{dD_{KL}(\mathcal{Q}\|\mathcal{P})}{dW_{i\mu}},$
(107)

where $\epsilon>0$ is the learning rate.

######  Definition 12.

Given an observable $O(\boldsymbol{\sigma},\boldsymbol{z})$ depending on the
variables $\\{\boldsymbol{\sigma},\boldsymbol{z}\\}$, we indicate with
$\langle\cdot\rangle_{\textrm{free}}$ the standard average with respect to
probability density (105)

$\langle
O\rangle_{\textrm{free}}:=\sum_{\\{\sigma,z\\}}\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z}|\boldsymbol{W})O(\boldsymbol{\sigma},\boldsymbol{z}),$
(108)

and with $\langle\cdot\rangle_{\textrm{clamped}}$ the average with respect to
the empirical probability density
$\mathcal{Q}(\boldsymbol{\sigma},\boldsymbol{z})$

$\langle
O\rangle_{\textrm{clamped}}:=\sum_{\\{\sigma,z\\}}\mathcal{Q}(\boldsymbol{\sigma},\boldsymbol{z})O(\boldsymbol{\sigma},\boldsymbol{z}).$
(109)

######  Remark 7.

By using the fact that

$\frac{dD_{KL}(\mathcal{Q}\|\mathcal{P})}{dW_{i\mu}}=\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{free}}-\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{clamped}},$ (110)

iteration (107) can be written as

$W_{i\mu}^{(n+1)}=W_{i\mu}^{(n)}+\epsilon\left(\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{clamped}}-\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{free}}\right).$ (111)

In order to perform the free average, we apply alternating Gibbs sampling [13]
and draw $\boldsymbol{\sigma}$ via the following probability density

$\mathcal{P}(\boldsymbol{\sigma}|\boldsymbol{z},\boldsymbol{W})=\frac{\mathcal{P}(\sigma,z|W)}{\sum_{\\{\sigma\\}}\mathcal{P}(\sigma,z|W)}=\prod_{i=1}^{N}\frac{1}{2}\left[1+\tanh\left(\sum_{\mu=1}^{K}z_{\mu}W_{i\mu}\sigma_{i}\right)\right],$
(112)

and analogously $\boldsymbol{z}$ via

$\mathcal{P}(\boldsymbol{z}|\boldsymbol{\sigma},\boldsymbol{W})=\frac{\mathcal{P}(\sigma,z|W)}{\sum_{\\{z\\}}\mathcal{P}(\sigma,z|W)}=\prod_{\mu=1}^{K}\frac{1}{2}\left[1+\tanh\left(\sum_{i=1}^{N}z_{\mu}W_{i\mu}\sigma_{i}\right)\right].$
(113)

Computationally it is convenient to realize the sampling according to the
simplified rules:

$\displaystyle\sigma_{i}^{(t+1)}$ $\displaystyle=$
$\displaystyle\mathrm{sign}\left[\tanh\left(\sum_{\mu=1}^{K}W_{i\mu}z_{\mu}^{(t)}\right)+\zeta_{i}^{t+1}\right]\quad\zeta_{i}^{t+1}\sim\mathcal{U}(-1,+1),$
(114) $\displaystyle z_{\mu}^{(t+1)}$ $\displaystyle=$
$\displaystyle\mathrm{sign}\left[\tanh\left(\sum_{i=1}^{N}W_{i\mu}\sigma_{i}^{(t+1)}\right)+\zeta_{\mu}^{t+1}\right]\quad\zeta_{\mu}^{t+1}\sim\mathcal{U}(-1,+1).$
(115)

######  Definition 13.

Computational setting (Grandmother Cell ansatz): All computational experiments
with the RBM (103) have been carried out by generating a synthetic dataset of
Rademacher archetypes $\\{\boldsymbol{\xi}^{\mu}\\}_{\mu=1,...,K}$ (see eq.
1), whence a set of corrupted examples $\mathcal{S}:=\\{\boldsymbol{\eta}^{\mu
a}\\}_{\mu=1,...,K}^{a=1,...,M}$ as

$\eta_{i\mu}^{a}=\xi_{i}^{\mu}\chi_{i}^{\mu a}\quad\chi_{i}^{\mu
a}\sim\frac{1+r}{2}\delta(\chi_{i}^{\mu
a}-1)+\frac{1-r}{2}\delta(\chi_{i}^{\mu a}+1)$ (116)

where $i=1,\cdots,N\quad\mu=1,\cdots,K,\quad a=1,\cdots,M$ and $r$ tunes the
quality of the dataset (see eqs. 2-3). We envisage the following supervised-
learning protocol to clamp the $\sigma$’s and the $z$’s variables: for every
new example presented to the networks, say $\boldsymbol{\eta}^{\nu,a}$, we set
$\boldsymbol{\sigma}=\boldsymbol{\eta}^{\nu,a}$ and
$\boldsymbol{z}=\boldsymbol{z}^{(\nu)}$, where $\boldsymbol{z}^{(\nu)}$ has
entries $z^{(\nu)}_{\mu}=\delta_{\mu\nu}$.

In this way the clamped average appearing in the learning rule (107) reads as

$\left\langle\sigma_{i}z_{\mu}\right\rangle_{\textrm{clamped}}=\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu
a}=\bar{\eta}_{i}^{\mu},$ (117)

while the free average appearing in the learning rule (107) can be estimated
by MC simulations or similar computational routes (here we followed the
“Persistent Contrastive Divergence” criterion, see e.g., [14]). In particular,
at the beginning of the simulation (in the first epoch) the set of variables
$(\boldsymbol{\sigma},\boldsymbol{z})$ is sampled randomly (i.e.,
$\sigma_{i}\in\\{-1,+1\\}$ is sampled by a symmetric Bernoulli distribution
with parameter $1/2$ and $z_{\mu}\in\mathbb{R}$ is sampled by a standard
Gaussian distribution, i.i.d. for any $i=1,...,N$ and $\mu=1,...,K$), then
they are updated and averaged at each epoch according to the rules (114)-(115)
for 10 steps, finally the learning rate $\epsilon$ are chosen according to the
Robbins-Monro scheme [15]: $\epsilon_{n}=\frac{50}{100+n}$ that is
$\sum_{n=1}^{\infty}\epsilon_{n}=\infty,\sum_{n=1}^{\infty}\epsilon_{n}^{2}<\infty$.
Finally, the magnetization of the archetype $\langle m\rangle$ is cheaply
evaluated via Hinton’s one-step MC approximation

$\langle
m\rangle=\frac{1}{K}\sum_{\mu=1}^{K}\left|\frac{1}{N}\sum_{i=1}^{N}\xi_{i}^{\mu}\tanh\left(W_{i\mu}\right)\right|.$
(118)

Even if this is only an approximation of the magnetization, it has very low
variance and computational cost and still captures the essential behavior of
the model. Results are presented in figure $2$ in the main text, where we
showed that, remarkably, the values of $\alpha,\beta,\rho$ that ensure a
successful information-processing are the same for both biological learning
(i.e., by the HN via Hebb’s rule) and artificial learning (i.e., by the RBM
via contrastive divergence).

As evidenced in [16, 17, 18] for structureless datasets, by training the RBM
following the grandmother-cell scheme one ends up with weight distributions
that are sharply peaked and, in particular, the expected value for $W_{i\mu}$
corresponds to $\bar{\eta}_{i}^{\mu}$, therefore, in this case, the setting
$W_{i\mu}=\bar{\eta}_{i}^{\mu}$ corresponds to a trained machine that is able
to reconstruct archetypes and to generate new examples. Notably, in the
structured case, the setting $W_{i\mu}=\bar{\eta}_{i}^{\mu}$ does not
correspond to an accomplished training but still it provides an effective pre-
training as highlighted in [18, 19].

We conclude this section with a last remark on the duality between HNs and
RBMs [20, 16, 21]. In fact, we understood that, as long as we identify the
weights in the latter with the patterns in the former, the RBM results to have
learnt the same patterns that the HN can retrieve and vice versa, also, the
two models share the same phase diagrams, hence where the HN is able to
retrieve, the RBM is able to correctly infer the archetype and perform
generalizations. However, one could argue that the HN must deal with binary
patterns (being these archetypes or examples) while, at the end of the RBM
training, we have no guarantee that the entries of the weights aree digital,
i.e. $W_{i\mu}\neq\pm 1$; this is not a real problem as we have shown in the
past that the duality is robust w.r.t. this kind of perturbation, namely the
HN keeps retrieving even if its entries are not sharply $\pm 1$ [22, 23].

##  V Maximum Entropy Approach for the supervised Hebbian learning

In this section we aim to reach the expression for the probability
distribution of the trained RBM, that is
$\mathcal{P}^{\textrm{\tiny{RBM}}}(\boldsymbol{\sigma},\boldsymbol{z}|\boldsymbol{W}=\bar{\boldsymbol{\eta}})$,
which, in turns, coincides with the probability distribution stemming from the
linearized partition function
$\mathcal{Z}^{\textrm{\tiny{HN}}}_{\beta=1}(\mathcal{S})$ found for the
supervised HN in Proposition 3, from another perspective. Let us recall the
learning rule for the RBM,

$W_{i\mu}^{(n+1)}=W_{i\mu}^{(n)}+\epsilon\left(\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{clamped}}-\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{free}}\right),$ (119)

which makes use of the correlation function $\langle
z_{\mu}\sigma_{i}\rangle_{\textrm{clamped}}$, thus this is the only
information that the RBM requires from a dataset. In this section we will show
how it is possible to retrieve both the RBM and the HN by looking for the
least structured density function that reproduces the clamped correlation
function. This construction requires the minimal number of constraints in
order to recover the key models of this work. Let us spell out the crucial
assumptions of this derivation:

  * •

the density function depends on the set of binary variables
$\boldsymbol{\sigma}\in\\{-1,+1\\}^{N}$ and on the set of real variables
$\boldsymbol{z}\in\mathbb{R}^{K}$,

  * •

we require that the expected value of the $L_{2}-$norm of the $\boldsymbol{z}$
variables is bounded;

  * •

we require that the density function correctly reproduces the correlations
between the $\boldsymbol{z}^{\prime}s$ and the $\boldsymbol{\sigma}$’s
variables.

######  Definition 14.

We introduce the following two-point empirical correlation functions

$\displaystyle C_{z^{2}}$ $\displaystyle=$
$\displaystyle\sum_{\mu=1}^{K}\overline{z_{\mu}^{2}},$ (120) $\displaystyle
C_{\sigma z}^{i,\mu}$ $\displaystyle=$
$\displaystyle\overline{\sigma_{i}z_{\mu}}.$ (121)

where the bar denotes the empirical average evaluated over the examples making
up the dataset $\mathcal{S}$.

######  Definition 15.

Given a probability density function
$\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z})$, the expectation of an
observable $O(\boldsymbol{\sigma},\boldsymbol{z})$ is defined as

$\langle
O\rangle_{\mathcal{P}}=\sum_{\sigma}\int\prod_{\mu=1}^{K}\left(\frac{dz_{\mu}}{\sqrt{2\pi}}\right)O(\boldsymbol{\sigma},\boldsymbol{z})\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z}),$
(122)

and, to lighten the notation, we introduce the trace operator Tr

$\text{Tr}[O(\boldsymbol{\sigma},\boldsymbol{z})]:=\sum_{\boldsymbol{\sigma}}\int\left[\prod_{\mu=1}^{K}\frac{dz_{\mu}}{\sqrt{2\pi}}\right]O(\boldsymbol{\sigma},\boldsymbol{z})$
(123)

such that

$\langle
O\rangle_{\mathcal{P}}=\text{Tr}[O(\boldsymbol{\sigma},\boldsymbol{z})\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z})].$
(124)

######  Theorem 2.

The least structured probability distribution reproducing the sets of
correlation functions $C_{z^{2}}$ and $C_{\sigma z}^{i,\mu}$ is

$\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z})=\frac{1}{\mathcal{Z}}\exp\left[-\frac{\lambda_{1}}{2}\sum_{\mu=1}^{K}z_{\mu}^{2}+\sum_{i,\mu=1}^{N,K}\Lambda_{i\mu}\sigma_{i}z_{\mu}\right]$
(125)

where $\mathcal{Z}$ is the normalizing factor, the parameters
$\lambda_{1},\\{\Lambda_{i\mu}\\}_{i=1,...,N}^{\mu=1,...,K}$ are Lagrangian
multipliers to be set in order constrain the correlation functions as
prescribed in Definition (14) and therefore their value is dataset dependent;
this probability distribution corresponds to the density function of a RBM
and, by marginalising w.r.t. the $z_{\mu}$ variables, we recover the density
function of a HN

$\mathcal{P}(\boldsymbol{\sigma})=\frac{1}{\mathcal{Z}}\exp\left[\frac{1}{2\lambda_{1}}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\Lambda_{i\mu}\sigma_{i}\Lambda_{j\mu}\sigma_{j}\right].$
(126)

###### Proof.

We introduce the following Lagrangian

$\begin{split}&S[\mathcal{P}]=-\text{Tr}(\mathcal{P}\log\mathcal{P})+\lambda_{0}[\text{Tr}(\mathcal{P})-1]+\frac{\lambda_{1}}{2}\left[\text{Tr}\left(\sum_{\mu=1}^{K}z_{\mu}^{2}\mathcal{P}\right)-C_{z^{2}}\right]+\\\
&+\sum_{i,\mu=1}^{N,K}\Lambda_{i\mu}\left[\text{Tr}(\sigma_{i}z_{\mu}\mathcal{P})-C_{\sigma
z}^{i,\mu}\right]\end{split}$ (127)

that accounts for maximum entropy and the above-mentioned constraints. By
extremizing $S[\mathcal{P}]$ with respect to $\mathcal{P}$ and
$\lambda_{0,1}$, $\Lambda_{i\mu}$, playing as the Lagrangian multipliers, we
obtain

$\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z})=\exp\left[1-\lambda_{0}-\frac{\lambda_{1}}{2}\sum_{\mu=1}^{K}z_{\mu}^{2}+\sum_{i,\mu=1}^{N,K}\Lambda_{i,\mu}\sigma_{i}z_{\mu}\right],$
(128)

and, posing $\mathcal{Z}=\exp(\lambda_{0}-1)$, we recover (125). By
integrating (128) with respect to $z_{\mu}$ we get the marginal, that can be
recast in the probability distribution for the HN

$\mathcal{P}(\boldsymbol{\sigma})\propto\exp\left[\frac{1}{2\lambda_{1}}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\Lambda_{i,\mu}\sigma_{i}\Lambda_{j,\mu}\sigma_{j}\right].$
(129)

∎

######  Remark 8.

The dataset-dependent expression for the Lagrangian multipliers $\lambda_{1}$,
$\Lambda_{i\mu}$ that reproduce the specific correlation functions of our
trained RBM model are

$\lambda_{1}=1,$ (130)
$\Lambda_{i\mu}=\sqrt{\frac{\beta}{N\mathcal{R}}}\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu
a}.$ (131)

In fact, with this choice, we get

$\mathcal{P}(\boldsymbol{\sigma},\boldsymbol{z}|\mathcal{S})\propto\exp\left[-\sum_{\mu=1}^{K}\frac{z_{\mu}^{2}}{2}+\sqrt{\frac{\beta}{N\mathcal{R}}}\frac{1}{M}\sum_{\mu=1}^{K}\sum_{i=1}^{N}\left(\sum_{a=1}^{M}\eta_{i}^{\mu
a}\right)\sigma_{i}z_{\mu}\right].$ (132)

and, by marginalizing over the $\boldsymbol{z}$’s variable, we recover the
probability distribution for the supervised HN defined in Eq. (8)

$\mathcal{P}(\boldsymbol{\sigma}|\mathcal{S})\propto\exp\left[\frac{1}{2}\frac{\beta}{N\mathcal{R}}\sum_{\mu=1}^{K}\sum_{i,j=1}^{N,N}\left(\frac{1}{M}\sum_{a=1}^{M}\eta_{i}^{\mu
a}\right)\left(\frac{1}{M}\sum_{a=1}^{M}\eta_{j}^{\mu
a}\right)\sigma_{i}\sigma_{j}\right].$ (133)

## References

  * [1] E. Schneidman, Elad, et al., Weak pairwise correlations imply strongly correlated network states in a neural population, Nature 440.7087:1007-1012, (2006).
  * [2] G.E. Hinton, Training products of experts by minimizing contrastive divergence, Neur. Comp. 14(8), 1771, (2002).
  * [3] E. Agliari, et al., The emergence of a concept in shallow neural networks, Neural Networks in press (available at arXiv: 2109.00454), (2022).
  * [4] M. Mezard, G. Parisi, M.A. Virasoro, Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, World Scientific Publishing Company, (1987).
  * [5] F. Guerra, Broken replica symmetry bounds in the mean field spin glass model, Comm. Math. Phys. 233(1):1-12, (2003).
  * [6] F. Guerra, et al., The replica symmetric approximation of the analogical neural network, J. Stat. Phys. 140(4):784-796, (2010).
  * [7] E.Agliari, F. Alemanno, A. Barra, A. Fachechi, Generalized Guerra’s interpolation schemes for dense associative neural networks, Neur. Netw. 128 254-267 (2020).
  * [8] P. Carmona, Y. Hu, Universality in Sherrington-Kirkpatrick spin glass model, Ann. Inst. H. Poincare (B) Prob. $\&$ Stat. 42, 2, (2006).
  * [9] G. Genovese, Universality in bipartite mean field spin glasses, J. Math. Phys. 53(12), 123304, (2012).
  * [10] E. Agliari, A. Barra, C. Longo, D. Tantati, Neural Networks Retrieving Boolean Patterns in a Sea of Gaussian Ones, J. Stat. Phys. 168,1085-1104 (2017).
  * [11] A.C.C. Coolen, R. Kühn, P. Sollich, Theory of neural information processing systems, Oxford Press (2005).
  * [12] A. Barra, G. Genovese, F. Guerra, Equilibrium statistical mechanics of bipartite spin systems, J. Phys. A 44(24):245002, (2011).
  * [13] C. Roussel, S. Cocco, R. Monasson, Barriers and dynamical paths in alternating Gibbs sampling of restricted Boltzmann machines, Phys. Rev. E 104(3), 034109 (2021) .
  * [14] T. Tieleman, Training Restricted Boltzmann Machines using Approximations to the Likelihood Gradient, Proceedings of the 25th International Conference on Machine Learning, Helsinki, Finland, 2008
  * [15] L. Cai, High-dimensional exploratory item factor analysis by a Metropolis-Hasting Robbins-Monro algorithm, Psychometrika 75.1:33-57, (2010).
  * [16] F.E. Leonelli, et al., On the effective initialisation for restricted Boltzmann machines via duality with Hopfield model, Neural Networks 143, 314-326, (2021).
  * [17] E. Agliari, F.E. Leonelli, C. Marullo, Storing, learning and retrieving biased patterns, Appl. Math. and Comp. 415 126716 (2022)
  * [18] A. Fachechi, E. Agliari, F. Alemanno, A. Barra, Dreaming Boltzmann machines outperform standard ones, submitted IEEE Trans. Neur. Netw. (2022)
  * [19] M. Smart, A. Zilman, On the mapping between Hopfield networks and Restricted Boltzmann Machines, International Conference on Learning Representations (2021).
  * [20] A. Barra, A. Bernacchia, E. Santucci, P. Contucci, On the equivalence among Hopfield neural networks and restricted Boltzmann machines, Neur. Net. 34:1-9, (2012).
  * [21] C. Marullo, E. Agliari, Boltzmann Machines as Generalized Hopfield Networks: A Review of Recent Results and Outlooks, Entropy 23, 34 (2021).
  * [22] A. Barra, G. Genovese, P. Sollich, D. Tantari, Phase transitions in Restricted Boltzmann Machines with generic priors, Phys. Rev. E 96(4), 042156, (2017).
  * [23] A. Barra, G. Genovese, P. Sollich, D. Tantari, Phase diagram of restricted Boltzmann machines and generalized Hopfield networks with arbitrary priors, Phys. Rev. E 97(2), 022310, (2018).