# Learning in associative networks through Pavlovian dynamics

1
Learning in associative networks through
Pavlovian dynamics

Miriam Aquaro1,2, Daniele Lotito2,3, Chiara Marullo2,4
1Dipartimento di Matematica, Sapienza Università di Roma, Roma, Italy
2GNFM-INdAM, Gruppo Nazionale di Fisica Matematica, Istituto Nazionale di Alta
Matematica, Italy
3Dipartimento di Informatica, Università di Pisa, Pisa, Italy
4National Research Council (CNR), Institute for High Performance Computing and
Networking (ICAR), Napoli, Italy

Keywords: Pavlov Conditioning, Hebbian Learning, Statistical Mechanics

Abstract

Hebbian learning theory is rooted in Pavlov’s Classical Conditioning. In
addition, while mathematical models of the former have been proposed and
studied in the past decades, especially in spin glass theory, only recently it
has been numerically shown that the Pavlovian neural and synaptic dynamics
mechanisms give rise to synaptic weights that correspond to the Hebbian
learning rule.
In this paper, we show that we can analytically prove this result. Moreover,
we obtain the Pavlovian dynamics directly with statistical mechanics
considerations, so that we can model the whole process within the same
framework. Then, we test and study the model in different numerical settings.
Finally, drawing from evidence on pure memory reinforcement during sleep
stages, we show how the proposed model can simulate neural networks that
undergo sleep-associated memory consolidation processes, thereby proving the
compatibility of Pavlovian learning.

##  1 Introduction

Hebbian learning theory has been extensively studied in various fields, and
the mathematical implementation of Hebb’s rule is crucial in many spin glass
models, the most famous being Hopfield model [Amit et al., 1985]. Hebb’s rule
is usually stated as “neurons that fire together wire together” [Shatz, 1992]
or, more formally, “any two cells or systems of cells that are repeatedly
active at the same time will tend to become ‘associated’ so that activity in
one facilitates activity in the other” [Cooper et al., 2013]. However, what
Hebb wrote in The Organization of Behaviour [Hebb, 1949] is “when an axon of
cell A is near enough to excite a cell B and repeatedly or persistently takes
part in firing it, some growth process or metabolic change takes place in one
or both cells such that A’s efficiency, as one of the cells firing B, is
increased”. The most known statement differs from the original one in the
sense that the former is based on mere contiguity, while the latter includes
temporal precedence, causality and contingency [Keysers and Gazzola, 2014].
The mathematical expression of Hebb’s rule in the Hopfield model is so
straightforward that it could be considered derived from either statement.
This is evident in the Hebbian prescription for the coupling matrix
$J_{ij}\propto\sum_{\mu=1}^{K}\xi_{i}^{\mu}\xi_{j}^{\mu}$, where the matrix
element $J_{ij}$ models the association between neuron $i$ and neuron $j$, and
this element is proportional to the number of patterns
$\\{\boldsymbol{\xi}^{\mu}\\}_{\mu=1,2,...,K}$ fed to the network whose $i$
and $j$ entries assume the same value.

However, when modeling neurophysiological phenomena such as Pavlov’s classical
conditioning, the simplification “neurons that fire together wire together”
falls short[Keysers and Gazzola, 2014]. Nonetheless, a rigorous mathematical
model of Pavlov’s conditioning that accounts for temporal precedence,
causality, and contingency can be constructed. Recent work in [Agliari et al.,
2023] demonstrates that systems of Ising neurons, influenced by an external
field, provide a meaningful mathematical representation of traditional
Pavlovian conditioning, and numerical experiments show that the related
learning process leads to the Hebbian prescription for storing network
memories. In particular, the network dynamically learns memories and stores
them in a coupling matrix whose function is coordinating the neural activity
and which is equal to the coupling matrix of the Hopfield model, deeply
studied in statistical mechanics.
It is thus natural to wonder if we can obtain the same Pavlovian dynamics
introduced in [Agliari et al., 2023] entirely within the statistical mechanics
framework. Here, with basic modeling assumptions, we write the partition
function for a network of neurons and synapses that allows to setup of the ODE
system that describes the coupled evolution of both neurons and synapses. This
leads to a mathematical model of Pavlov’s Classical Conditioning using only
spin glass theory tools, which we will outline briefly in Sec. 2. The
resulting Pavlovian dynamics is the same found in [Agliari et al., 2023]. We
then analytically prove its convergence to the Hebbian kernel, a result that
was only numerically supported in the aforementioned work. Subsequently, we
analyze settings that more closely mirror what happens in nature, in which the
stimuli presented to the neurons are not chosen from the uniform distribution
across possible stimuli.
Once having established that this Pavlovian model for associative memory
implies Hebbian learning, we argue that this model is suited for describing
biological networks. When modelling biological learning networks, we need to
check the robustness of the resulting model with respect to the fundamental
features of biological learning. One of these features is undoubtedly sleep,
which biological evidence has shown promotes the reactivation of neural
circuits involved in learning, leading to the strengthening of synaptic
connections associated with newly acquired information [Crick and Mitchison,
1983]. This process, often referred to as memory replay or neuronal replay,
occurs predominantly during specific stages of sleep, such as slow wave sleep
(SWS) and rapid eye movement sleep (REM). In particular, these two sleep
phases serve a critical function in filtering and prioritizing memories, which
helps prevent cognitive overload and ensures efficient storage of relevant
information. This selective retention mechanism involves the strengthening of
important memories while weakening or discarding irrelevant or spurious ones,
a process essential for optimal cognitive functioning. This led to the
development of computational models that aimed to mimic these biological
learning dynamics [Hassabis et al., 2017].

In the associative neural networks literature this function is accomplished by
the dreaming mechanism, a well-established and biologically plausible way to
mitigate crosstalk effects between memories [Fachechi et al., 2019a, Agliari
et al., 2019, Aquaro et al., 2022]. Thus, in the spirit of describing
biological networks within our framework, we have to demand that the model can
describe neural networks that sleep, we show that this is possible. The
resulting network can handle larger information content and more effectively
than in the Hebbian case.

The paper is structured as follows: Sec. 2 introduces the spin glass modeling
of Pavlov’s Classical Conditioning, which enables the formulation of a
multiscale system of ordinary differential equations representing neural and
synaptic Pavlovian dynamics. Sec. 3 discusses the convergence of the dynamics
to the Hebbian prescription for the coupling matrix of associative neural
networks and provides a discrete version of this Pavlovian dynamics. Sec. 4
generalizes the dynamics of the model, showing that, using the same
mathematical tools, it is possible to derive the so-called dreaming
prescription for the coupling matrix. Sec. 5 is dedicated to numerical checks
and experiments analyzing the convergence of the coupling matrix under
different scenarios where the stimuli are not uniformly selected. Finally,
Sec. 6 is left for conclusions and outlooks.

##  2 A statistical mechanics view on neural and synaptic Pavlovian dynamics

In this section, we introduce the statistical mechanics formulation of
Pavlov’s classical conditioning, a more formal derivation is present in
[Agliari et al., 2023], and it is discussed in this paper in Appendix A.
We consider a system composed of $N$ binary neurons whose configurations are
described by $\boldsymbol{\sigma}\in\Sigma_{N}:=\\{-1,1\\}^{N}$. The system’s
thermodynamics is characterized by the following 2-point Hamiltonian:

$H_{N}(\boldsymbol{\sigma},\boldsymbol{J})=-\frac{1}{2}\sum_{i\neq
j=1}^{N,N}J_{ij}\sigma_{i}\sigma_{j}-u\sum_{i=1}^{N}h_{i}\sigma_{i},$ (1)

where $\boldsymbol{J}\in\mathbb{R}^{N\times N}$ represents the symmetric
coupling matrix describing pairwise neuron interactions,
$h_{i}\in\\{-1,0,1\\}$ is the bias or, in more biological terms, the firing
threshold acting on each neuron and $u$ is a global amplification factor that
modulates the strength of this bias. Let us now extend the configuration space
to $\Sigma_{N}\times\mathcal{J}_{N}$, where $\mathcal{J}_{N}$ represent the
space of the $N\times N$ symmetric matrices with null diagonal elements and
equip this space with a Boltzmann-Gibbs measure

$\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J})=\frac{1}{Z_{N,\beta}}\exp\big{[}-\beta
H_{N}(\boldsymbol{\sigma},\boldsymbol{J})\big{]},$ (2)

where

$Z_{N,\beta}=\sum_{\boldsymbol{J}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma}\in\Sigma_{N}}\exp\big{[}-\beta
H_{N}(\boldsymbol{\sigma},\boldsymbol{J})\big{]}$ (3)

is the associated partition function. The parameter $\beta$ governs the
broadness of the distribution, dictating the noise level affecting the neuron
dynamics. As $\beta\to 0$, the distribution becomes flat, and each neuron is
randomly oriented. Conversely, as $\beta\to\infty$, the distribution sharply
peaks at the configurations that minimize the Hamiltonian. In the classical
statistical mechanics framework applied to neural networks, such as in pattern
recognition applications [Amit, 1989, Coolen et al., 2005], synapses are
considered as _quenched_ variables, meaning that they remain fixed.
Consequently, the summation over $\boldsymbol{J}\in\mathcal{J}_{N}$ will not
appear in the definition of the partition function (3). However, in this work,
we aim to explore the coupled dynamical evolution of neurons
$\boldsymbol{\sigma}$ and synapses $\boldsymbol{J}$. This scenario requires
that we do not treat the couplings as quenched variables. Instead, the
partition function should now include the summation over synapses
configurations, and this corresponds to the so-called _annealed expression_
[Marc et al., 1987].

We now consider an interaction matrix space characterized by the Rademacher
probability measure, i.e.

$P(J_{ij}=\pm 1)=1/2,\quad\text{for}\ i,j=1,\dots,N,\quad i\neq j.$ (4)

In other words, we are choosing a coupling matrix whose entries are i.i.d.
random variables drawn from a Rademacher distribution. To simplify
computations, we add in the partition function from (3) the two source terms
$\boldsymbol{t}$ and $\boldsymbol{T}$ for neurons and weights respectively

${Z}_{N}(\beta,\boldsymbol{t},\boldsymbol{T})=\sum_{\boldsymbol{\sigma}}\sum_{\boldsymbol{J}}\exp\Big{(}\beta\sum_{i<j}^{N}J_{ij}\sigma_{i}\sigma_{j}+\beta
u\sum_{i=1}^{N}h_{i}\sigma_{i}+\sum_{i=1}^{N}t_{i}\sigma_{i}+\sum_{i<j}^{N}T_{ij}J_{ij}\Big{)},$
(5)

where the sum over the entries of the coupling matrix is explicitly given by

$\sum_{\boldsymbol{J}}\equiv\prod_{i<j}\sum_{J_{ij}=\pm 1}.$

Summing over all $J_{ij}$ variables in eq. (5) results in

${Z}_{N}(\beta,\boldsymbol{t},\boldsymbol{T})=\sum_{\boldsymbol{\sigma}}\exp\Big{(}\sum_{i=1}^{N}\left(\beta
uh_{i}+t_{i}\right)\sigma_{i}\Big{)}\prod_{i<j}2\cosh\left(\beta\sigma_{i}\sigma_{j}+T_{ij}\right).$
(6)

By deriving with respect to the source $T_{ij}$ and then setting both the
source terms to zero, we find

$\langle
J_{ij}\rangle=\frac{\partial\log{{Z}_{N}(\beta,\boldsymbol{t},\boldsymbol{T})}}{\partial
T_{ij}}\Big{|}_{\boldsymbol{T}=0,\boldsymbol{t}=0}=\langle\sigma_{i}\sigma_{j}\rangle\tanh\beta\approx\langle\sigma_{i}\rangle\langle\sigma_{j}\rangle\tanh\beta,$
(7)

where we used the fact that $\tanh(\sigma x)=\sigma\tanh x$ for $\sigma=\pm 1$
and the mean-field assumption under which each neuron is subject to a net
effective external field, rather than to the interactions with other spins.
Similarly, the spin activity expectation $\langle\sigma_{i}\rangle$ yields:

$\langle\sigma_{i}\rangle=\frac{\partial\log{{Z}_{N}(\beta,\boldsymbol{t},\boldsymbol{T})}}{\partial
t_{i}}\Big{|}_{\boldsymbol{T}=0,\boldsymbol{t}=0}\approx\tanh\Big{(}\beta\sum_{j\neq
i}^{N}\langle J_{ij}\rangle\langle\sigma_{j}\rangle+\beta uh_{i}\Big{)},$ (8)

where the mean-field assumption is applied again in the last equality.
Given the equivalence of all spin indices, the results from Eqs. (7) and (8)
remain consistent for all $i$ and $j$, with $i\neq j$.

In order to lighten the notation we set $\langle J_{ij}\rangle=J_{ij}$ and
$\langle\sigma_{i}\rangle=\sigma_{i}$, and we directly write the dynamics of
the system as follows

$\displaystyle\dot{{\sigma_{i}}}$
$\displaystyle=-\frac{1}{\tau}\sigma_{i}+\frac{1}{\tau}\tanh\Big{(}\beta\sum_{j\neq
i}^{N}\ J_{ij}\sigma_{j}\beta uh_{i}\Big{)},$ (9)
$\displaystyle\dot{{J_{ij}}}$
$\displaystyle=-\frac{1}{\tau^{\prime}}J_{ij}+\frac{1}{\tau^{\prime}}\sigma_{i}\sigma_{j}\tanh\beta,$
(10)

with $\tau$ and $\tau^{\prime}$ being the neural and synaptic time-scales. We
assume that $\tau^{\prime}\gg\tau$, as synaptic dynamics is characterized by a
longer time-scale in nature [Maass and Zador, 1997].

##  3 Pavlovian dynamics implies Hebbian learning

In this section, we will focus on the evolution of the synaptic coupling
matrix derived in the previous section, showing that the network learns the
associations between patterns and converges to the Hebbian prescription. We
will prove this convergence both in the case of the continuous dynamics
described by the Eq.(10) (see Sec. 3.1) and for its discrete approximation
(see Sec. 3.2 and Appendix B). Finally, in Sec. 3.3 we will present the
setting in which the patterns are extracted from a given probability
distribution. In this case, the coupling matrix probability distribution
converges to the distribution of the asymptotic states, whose mean and
variance are analytically computed.

###  3.1 Convergence in the continuous case

Let us start by analysing the dynamics described by Eq. (10). Assuming that
the synaptic dynamics leads to a steady state, without loss of generality, we
can set $\dot{{J_{ij}}}=0$ thus obtaining
$J_{ij}=\tanh(\beta)\sigma_{i}\sigma_{j}$. We can also enforce two singular-
bit stimuli, $\xi_{i}$ and $\xi_{j}$, onto the neurons $\sigma_{i}$ and
$\sigma_{j}$, which means that we can set $\sigma_{i}=\xi_{i}$ and
$\sigma_{j}=\xi_{j}$. We can also directly set
$\boldsymbol{\sigma}=\boldsymbol{\xi}$; we call this procedure presenting a
pattern to the network. We observe that in this case the synaptic matrix
elements

$J_{ij}=\tanh(\beta)\,\xi_{i}\xi_{j}\,$ (11)

resembles the Hebbian prescription. In Eq. (11), the parameter $\beta$
controls the coupling strength. When $\beta\to 0$, the noise dominates the
stimuli, preventing any learning. Otherwise, if $\beta\to+\infty$ we have that
$\tanh(\beta)\approx 1$ and (11) corresponds to the Hebbian prescription

$J_{ij}=\frac{1}{K}\sum_{\mu=1}^{K}\xi^{\mu}_{i}\xi^{\mu}_{j},$ (12)

in case of a single pattern $K=1$.
We also notice that, in this simple case, we have full knowledge of the
synaptic matrix dynamics. In fact, if we fix
$\boldsymbol{\sigma}=\boldsymbol{\xi}^{1}$, and denote $J(t):=J_{ij}(t),\
J_{0}=J(t=0)\ c:=\xi_{i}^{1}\xi_{j}^{1}\tanh\beta$, Eq. (10) reduces to the
differential equation

$\frac{dJ}{dt}=-\frac{1}{\tau^{\prime}}J(t)+\frac{c}{\tau^{\prime}}.$ (13)

We identify this equation as a first-order linear ordinary differential
equation. The general solution to such an equation combines the homogeneous
solution and a particular solution. We observe that the solution of the
homogeneous equation $\frac{dJ}{dt}+\frac{1}{\tau^{\prime}}J(t)=0$ is
$J(t)=Ae^{-\frac{1}{\tau^{\prime}}t}$, where $A$ is a constant determined by
initial conditions. At the same time, the solution of the particular equation
$-\frac{1}{\tau^{\prime}}J+\frac{c}{\tau^{\prime}}=0$ is the constant value
$c$, we get $J(t)=Ae^{-\frac{1}{\tau^{\prime}}t}+c$, where $A=J_{0}-c$. The
solution to the differential equation (13) can then be written as

$J(t)=(J_{0}-c)e^{-\frac{1}{\tau^{\prime}}t}+c.$ (14)

Thus, $J(t)$ approaches exponentially the steady-state value $c$, starting
from $J_{0}$, with a rate determined by $\tau^{\prime}$.
In the case of multiple patterns (i.e. $K>1$) it is sufficient to replace $c$
with its temporal average over the synaptic scale $\tau^{\prime}$. We thus
have

$c=c(t)\equiv\frac{1}{\tau^{\prime}}\int_{t-\tau^{\prime}/2}^{t+\tau^{\prime}/2}\text{d}\tau\,\sigma_{i}(\tau)\sigma_{j}(\tau)\tanh\beta.$
(15)

Since the neural time scale is well separated from the synaptic one, assuming
that, regardless of $t$ and $\tau^{\prime}$, in the temporal window neurons
are aligned to each pattern $\boldsymbol{\xi}^{\mu}$ an amount of time that is
proportional to the probabilities $\\{p_{\mu}\\}_{\mu=1,\ldots,K}$, we can
simply set in (14)
$c=\sum_{\mu=1}^{K}p_{\mu}\xi^{\mu}_{i}\xi^{\mu}_{j}\tanh{\beta}$. Then, if
$\beta\gg 1$ and $p_{\mu}=\frac{1}{K}$ for each pattern, we recover the
Hebbian coupling prescription (12).

###  3.2 Convergence in the discrete case

Building upon the convergence results established in the preceding section,
the objective now is to move towards a discrete formulation of Eqs. (9) and
(10), showing that, under different settings, the convergence to the Hebbian
prescription (12) continues to yield.

Let us start by writing the discrete version of Eqs. (9) and (10) as in
[Agliari et al., 2023]:

$\displaystyle\sigma_{i}^{(n+1)}$
$\displaystyle=\sigma_{i}^{(n)}\left(1-\frac{\delta
t}{\tau}\right)+\frac{\delta t}{\tau}\tanh\left(\beta\sum_{j\neq
i}^{N}J_{ij}^{(n)}\sigma_{j}^{(n)}+\beta uh_{i}^{(n)}\right),$ (16)
$\displaystyle J_{ij}^{(n+1)}$ $\displaystyle=J_{ij}^{(n)}\left(1-\frac{\delta
t}{\tau^{\prime}}\right)+\frac{\delta
t}{\tau^{\prime}}\sigma_{i}^{(n)}\sigma_{j}^{(n)}\tanh\beta,$ (17)

where $n$ denotes the discrete time step, $\delta t$ represents the
differential time unit responsible for system evolution and $h_{i}^{(n)}$
corresponds to an external stimulus that varies over time.
Since we are mainly interested in studying the convergence of the coupling,
and we discussed that $\tau^{\prime}$ is substantially larger than $\tau$, we
choose to set $\delta t=\tau$, that corresponds to allow the spins to align
instantly with the stimulus. In this way we obtain

$\displaystyle\sigma_{i}^{(n+1)}$ $\displaystyle=$
$\displaystyle\tanh\big{(}\beta\sum_{j\neq
i}^{N}J_{ij}^{(n)}\sigma_{j}^{(n)}+\beta uh_{i}^{(n)}\big{)},$ (18)
$\displaystyle J_{ij}^{(n+1)}$ $\displaystyle=$ $\displaystyle
J_{ij}^{(n)}\left(1-\frac{\tau}{\tau^{\prime}}\right)+\frac{\tau}{\tau^{\prime}}\sigma_{i}^{(n)}\sigma_{j}^{(n)}\tanh\beta.$
(19)

We remark that the derivation presented in the previous section is not unique.
For the reader that seeks a more formal derivation, we have sketched in
Appendix A the master equation approach presented in [Agliari et al., 2023].

It has been shown numerically that if one selects external stimuli
$\\{\boldsymbol{\xi}^{\mu}\\}_{\mu=1}^{K}$ by sampling uniformly the index
$\mu$ from the index set $\\{1,\dots,K\\}$, the dynamics governed by
(18)-(19), yields a coupling matrix $\boldsymbol{J}$ that, for long times,
tends to converge to Hebb’s prescription [Agliari et al., 2023].

We now show that the analytical results obtained in Sec. 3.1 hold also in the
discrete case. Moreover, we compute the amplitude of the fluctuation around
the steady state in the general case of non-uniformly presented patterns.
Let us focus on (19), we first analyze the situation where
$\frac{\tau}{\tau^{\prime}}<1$ and the neurons are clamped to certain fixed
values111Previous work only report analytical estimation in this case, see
[Agliari et al., 2023] Appendix A.2.
$\boldsymbol{\sigma}=\boldsymbol{\xi}^{1}$, then we generalize to the case
where neuron values change. Similarly to what we have done in the continuous
case let us define, for the sake of convenience, the following auxiliary
parameters

$\epsilon:=\frac{\tau}{\tau^{\prime}},\quad
c:=\xi^{1}_{i}\xi^{1}_{j}\tanh\beta$

so that Eq. (19) can be written as

$J_{ij}^{(n+1)}=(1-\epsilon)J_{ij}^{(n)}+\epsilon c\,.$ (20)

Observing that no stochasticity is involved in the dynamics, for each $n$, we
can replace the value of $J^{(n)}_{ij}$ with that from the previous iterations
$n-1$. Proceed recursively, we obtain the following

$\displaystyle J_{ij}^{(n+1)}$
$\displaystyle=(1-\epsilon)[(1-\epsilon)J_{ij}^{(n-1)}+\epsilon c]+\epsilon c$
$\displaystyle\ \ \vdots$ $\displaystyle=(1-\epsilon)^{n}J_{ij}^{(0)}+\epsilon
c\sum_{k=0}^{n}(1-\epsilon)^{k}.$

Let us now observe that, since $J_{ij}^{(0)}$ is finite and $|1-\epsilon|<1$
the first term on the right-hand side goes to zero as $n$ approaches infinity.
The second term, instead, being a geometric series with a ratio
$1-\epsilon<1$, will converge to $\frac{1}{\epsilon}$.
Bringing everything together and substituting the explicit value of $c$, we
get

$\lim_{n\to\infty}J_{ij}^{(n+1)}=\xi^{1}_{i}\xi^{1}_{j}\tanh\beta,$ (21)

which corresponds to the single pattern Hebbian prescription. Let us remark
that, the hyperbolic tangent in Eq.(21) acts as a global attenuating factor
that accounts for the thermal noise; for sufficiently low temperatures we have
$\tanh\beta\approx 1$.

If we have two patterns $K=2$, then the coupling matrix converges to the mixed
Hebbian kernel

$J_{ij}:=\frac{1}{2}\left(\xi_{i}^{1}\xi_{j}^{1}+\xi_{i}^{2}\xi_{j}^{2}\right)\tanh\beta\,.$
(22)

In Appendix B, we show how this result can be easily achieved by clamping
neurons to pattern $\boldsymbol{\xi}^{1}$ during even iteration steps and to
pattern $\boldsymbol{\xi}^{2}$ during odd steps. We also discussed how the
same logic can be extended to the case with an arbitrary number of patterns
cyclically presented to the network.

###  3.3 Characterization of the steady state

Suppose stimuli are not presented in periodic cycles. In that case, we can
only make statements about the probability distribution of the coupling matrix
elements as the number of iteration steps goes to infinity. Furthermore, in
the real world, neural systems do not receive uniformly distributed external
signals. In fact, certain concepts or stimuli will appear more frequently than
others. Hence, it’s relevant to ask if similar results are obtainable with a
heterogeneous presentation of external stimuli.
To deal with this more realistic setting, now each pattern, denoted as
$\boldsymbol{\xi}^{\mu}$ with $\mu\in\\{1,2,\ldots,K\\}$, has a certain
probability $p_{\mu}$ of being presented at any iteration step. In principle,
this probability can vary with the iteration step index. Thus, the update rule
for the coupling matrix is

$J_{ij}^{(n+1)}=(1-\epsilon)J_{ij}^{(n)}+\epsilon c^{\mu(n)},$ (23)

where $c^{\mu(n)}:=\xi_{i}^{\mu(n)}\xi_{j}^{\mu(n)}\tanh\beta$ corresponds to
the pattern presented at step $n$, and $\mu(n)$ is now a random variable
representing the index of the pattern chosen according to the probabilities
$\\{p_{\mu}(n)\\}_{\mu=1,\ldots,K}$. Obviously, it always holds
$\sum_{\mu}p_{\mu}(n)=1$, regardless of $n$.
To analyze the expected behavior of the coupling matrix over time, we compute
the expected value $\mathbb{E}[J_{ij}^{(n+1)}]$,

$\mathbb{E}[J_{ij}^{(n+1)}]=(1-\epsilon)\mathbb{E}[J_{ij}^{(n)}]+\epsilon\mathbb{E}[c^{\mu(n)}].$
(24)

From this equation, we see that the expected contribution of a randomly chosen
pattern at step $n$ is
$\mathbb{E}[c^{\mu(n)}]=\sum_{\mu=1}^{K}p_{\mu}(n)c^{\mu(n)}$.
Moreover, we are interested in the case where the pattern distribution does
not vary with time so that $c^{\mu(n)}\equiv c_{\mu}$, $p_{\mu}(n)\equiv
p_{\mu}$, for every $\mu$ and $n$, and
$\mathbb{E}[c_{\mu}]=\sum_{\mu=1}^{K}p_{\mu}c^{\mu}$. Then, we can simplify
the equation to:

$\mathbb{E}[J_{ij}^{(n+1)}]=(1-\epsilon)\mathbb{E}[J_{ij}^{(n)}]+\epsilon\sum_{\mu=1}^{K}p_{\mu}c^{\mu}.$
(25)

Thus the expected value of the coupling matrix elements evolves towards the
weighted average of the contributions from all $K$ patterns, with the weights
given by their respective probabilities. This can be shown with recursive
substitutions, as before. Instead, here we observe that, as $n$ approaches
infinity, in the steady state
$\mathbb{E}[J_{ij}^{(n+1)}]=\mathbb{E}[J_{ij}^{(n)}]=\mathbb{E}[J_{ij}^{(\infty)}]$,
we easily obtain:

$\mathbb{E}[J_{ij}^{(\infty)}]=\frac{\epsilon}{1-(1-\epsilon)}\sum_{\mu=1}^{K}p_{\mu}c^{\mu}=\sum_{\mu=1}^{K}p_{\mu}c^{\mu}.$
(26)

We underline that the expected value of $c^{\mu}$ is equal to the expected
value of the coupling matrix element $J_{ij}^{(\infty)}$222To shorten the
notation we have not displayed the neuron indexes in $c^{\mu}\equiv
c^{\mu}_{ij}$. , and in the low-temperature regime they coincide with the
generalized Hebbian prescription

$J_{ij}=\sum_{\mu=1}^{K}p_{\mu}\xi^{\mu}_{i}\xi^{\mu}_{j}.$ (27)

To compute the variance, we start with the definition:

$\text{Var}[J_{ij}^{(n+1)}]=\mathbb{E}[(J_{ij}^{(n+1)}-\mathbb{E}[J_{ij}^{(n+1)}])^{2}].$
(28)

Substituting the update rule for $J_{ij}^{(n+1)}$ into the variance formula
and expanding, we get:

$\displaystyle\text{Var}[J_{ij}^{(n+1)}]$
$\displaystyle=\mathbb{E}\left[\left((1-\epsilon)J_{ij}^{(n)}+\epsilon
c^{\mu(n)}-\mathbb{E}[(1-\epsilon)J_{ij}^{(n)}+\epsilon
c^{\mu(n)}]\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left((1-\epsilon)(J_{ij}^{(n)}-\mathbb{E}[J_{ij}^{(n)}])+\epsilon(c^{\mu(n)}-\mathbb{E}[c^{\mu(n)}])\right)^{2}\right].$

We expand this expression and apply the linearity of expectation,

$\displaystyle\text{Var}[J_{ij}^{(n+1)}]$
$\displaystyle=(1-\epsilon)^{2}\text{Var}[J_{ij}^{(n)}]+\epsilon^{2}\text{Var}[c^{\mu(n)}]$
$\displaystyle\quad+2(1-\epsilon)\epsilon\mathbb{E}[(J_{ij}^{(n)}-\mathbb{E}[J_{ij}^{(n)}])(c^{\mu(n)}-\mathbb{E}[c^{\mu(n)}])].$

Given that $J_{ij}^{(n)}$ and $c^{\mu(n)}$ are independent (assuming the
pattern at step $n$ is independent of the coupling matrix’s previous state),
the last term vanishes.

$\text{Var}[J_{ij}^{(n+1)}]=(1-\epsilon)^{2}\text{Var}[J_{ij}^{(n)}]+\epsilon^{2}\text{Var}[c^{\mu(n)}].$
(29)

The variance of $c^{\mu(n)}$ is equal to the expected value of the square of
$c^{\mu(n)}$ minus the square of the expected value of $c^{\mu(n)}$.

$\text{Var}[c^{\mu(n)}]=\sum_{\mu=1}^{K}p_{\mu}\left(c^{\mu(n)}\right)^{2}-\left(\sum_{\mu=1}^{K}p_{\mu}c^{\mu(n)}\right)^{2}.$
(30)

Under the assumption that the pattern distribution does not vary with time, as
$n$ approaches infinity, since the system reaches a steady state, the variance
of the coupling matrix elements stabilize, and we obtain

$\text{Var}[J_{ij}^{(\infty)}]=\frac{\epsilon^{2}}{1-(1-\epsilon)^{2}}\left[\sum_{\mu=1}^{K}p_{\mu}\left(c^{\mu}\right)^{2}-\left(\sum_{\mu=1}^{K}p_{\mu}c^{\mu}\right)^{2}\right].$
(31)

Equivalently, we can leave implicit the variance of $c^{\mu}$ and simplify the
initial factor.

$\text{Var}[J_{ij}^{(\infty)}]=\frac{\epsilon}{2-\epsilon}\text{Var}[c^{\mu}]=\left(\frac{2}{2-\epsilon}-1\right)\text{Var}[c^{\mu}]$
(32)

When the ratio between the neuron and the synaptic time scales goes to zero,
the variance of the steady-state coupling matrix element goes to zero. This is
because $\epsilon$ regulates the system response for any single pattern
presentation, the lower $\epsilon$ the lower the variation of
$J_{ij}^{(\infty)}$ during an updating step. Thus, in the low $\epsilon$
regime, to make $J_{ij}^{(\infty)}$ significantly different from its mean, the
network should be exposed to an atypical pattern sequence, e.g. the index
random variable $\mu$ is uniform and we draw the same index multiple times.
Actually, in this limit, we can Taylor expand Eq. (32) to get

$\text{Var}[J_{ij}^{(\infty)}]\approx\left(1+\frac{\epsilon}{2}+O(\epsilon^{2})\right)\text{Var}[c^{\mu}]-\text{Var}[c^{\mu}]=\frac{\epsilon}{2}\text{Var}[c^{\mu}]+O(\epsilon^{2})\text{Var}[c^{\mu}].$
(33)

Neglecting the terms of $O(\epsilon^{2})$, we first observe that the amplitude
of the fluctuation is proportional to
$\sqrt{\epsilon}=\sqrt{\frac{\tau}{\tau^{\prime}}}$, as discussed in [Agliari
et al., 2023] under stronger assumptions. Then, we consider the limit
$\beta\rightarrow\infty$ to write

$\text{Var}[J_{ij}^{(\infty)}]\approx\frac{1}{2}\frac{\tau}{\tau^{\prime}}\left[1-\left(\sum_{\mu=1}^{K}p_{\mu}\xi^{\mu}_{i}\xi^{\mu}_{j}\right)^{2}\right],$
(34)

where we used the fact that in the low-temperature limit $(c^{\mu})^{2}=1$. In
Eq. 34, we recognize the Hebbian matrix element 27. In particular, the square
root of (34), i.e. $\Delta_{ij}:=\sqrt{\text{Var}[J_{ij}^{(\infty)}]}$ is the
expectation of the displacement between $J_{ij}^{(\infty)}$ and the Hebbian
prescription given in Eq. 27. Under stronger assumptions, the authors of
[Agliari et al., 2023] obtain the following formula for $\Delta_{ij}$ in the
single pattern case

$\Delta_{ij}=\sqrt{\frac{1}{2}\frac{\tau}{\tau^{\prime}}\tanh^{2}\beta\left[1-\left(\xi^{\mu}_{i}\xi^{\mu}_{j}\right)^{2}\right]+\left[\xi^{\mu}_{i}\xi^{\mu}_{j}(1-\tanh\beta)\right]^{2}}.$
(35)

We have just proven that the multiple patterns generalization

$\Delta_{ij}=\sqrt{\frac{1}{2}\frac{\tau}{\tau^{\prime}}\tanh^{2}\beta\left[1-\left(\sum_{\mu=1}^{K}p_{\mu}\xi^{\mu}_{i}\xi^{\mu}_{j}\right)^{2}\right]+\left[\sum_{\mu=1}^{K}p_{\mu}\xi^{\mu}_{i}\xi^{\mu}_{j}(1-\tanh\beta)\right]^{2}}$
(36)

holds, at least in the limit $\beta\rightarrow\infty$, since it is equivalent
to (34).
Finally, we can obtain either from Eq. 34, or from the low-temperature limit
of Eq. 35 the normalized Frobenius norm333Let us recall that, given a $n\times
m$ matrix $\boldsymbol{A}$, the Frobenius norm is defined as
$\lVert\boldsymbol{A}\rVert:=\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{m}A_{ij}^{2}}$,
the adjective normalized comes from the fact that we have divided the value of
the Frobenius norm by $\sqrt{m\times n}$.

$\lVert\boldsymbol{\Delta}\lVert_{F}=\frac{1}{N}\sqrt{\frac{1}{2}\frac{\tau}{\tau^{\prime}}\sum_{i,j=1}^{N}\left[1-\left(\sum_{\mu=1}^{K}p_{\mu}\xi^{\mu}_{i}\xi^{\mu}_{j}\right)^{2}\right]}.$
(37)

##  4 Beyond Hebb’s paradigm: the dreaming prescription

Until now, we have been working in the mean-field limit, the same
approximation used on the standard Hopfield network in Amit, Gutfreund and
Sompolinsky (AGS) theory [Amit, 1989]. Indeed, it is well-known that the
Hebbian prescription (resulting in the celebrated Hopfield model [Hopfield,
1982]) is not the optimal construction for the coupling matrix
$\boldsymbol{J}$. In particular, from AGS theory [Amit, 1989, Coolen et al.,
2005] we see that the network can handle at most with $K\sim 0.14N$ patterns,
while – from Gardner’s theory – we know that spin-glass models with pairwise
interactions can manage a number of patterns $K\sim 2N$ (or $K\sim N$ for
symmetric coupling matrices) [Gardner and Derrida, 1988]. In particular,
Hopfield model fails to retrieve stored information when the spurious states,
corresponding to mixtures of memories, get exponentially more abundant than
pure memories. Various approaches have been used to solve this problem, we
cite unlearning algorithms [Kohonen, 1984, Personnaz et al., 1985, Kanter and
Sompolinsky, 1987a, Hemmen et al., 1990, Plakhov and Semenov, 1992, Fachechi
et al., 2019b], which are based on the idea of increasing the energy of
spurious state configurations to make them less stable. In particular, the
unlearning algorithm of Plakhov and Semenov [Plakhov and Semenov, 1992]
assumes that the synaptic matrix is updated by progressively removing the
correlations between the total fields acting on the spins, rather than the
correlations between spins directly. Inspired by this algorithm, we generalize
the modelling discussed in previous sections to make it suitable for framing
more performing, and biologically plausible, learning mechanisms. To this
extent, we aim to adapt the dynamics (9)-(10) so that the coupling matrix, by
learning directly from the total fields acting on the spins, converges to the
so-called Kanter and Sompolinsky projector [Kanter and Sompolinsky, 1987b]
which is given by

$J_{ij}^{(Dream)}:=\frac{1}{K}\sum_{\mu=1}^{K}\sum_{\nu=1}^{K}\xi_{i}^{\mu}{(C^{-1})_{\mu\nu}}\xi_{j}^{\nu}$
(38)

where

$C_{\mu\nu}:=\frac{1}{N}\sum_{i=1}^{N}\xi_{i}^{\mu}\xi_{i}^{\nu},$ (39)

is the pattern correlation matrix. This network can handle a number of
patterns $K\sim N$, thus saturating Gardner’s theoretical bound.
The generalisation of Kanter and Sompolinsky’s projector is the Reinforcement
& Removal algorithm [Fachechi et al., 2019b], also called the dreaming
prescription; hence the label Dream in Eq. (38). If we want the synaptic
dynamics given in Eq. (10) to converge to Eq. (38), we need to modify the way
in which external stimuli act on the system. To illustrate this, let us notice
that

$J_{ij}^{(Dream)}=\frac{1}{K}\boldsymbol{\xi}_{i}^{T}{\boldsymbol{C}}^{-1}\boldsymbol{\xi}_{j}=\frac{1}{K}\left({\boldsymbol{C}^{-1/2}}\boldsymbol{\xi}_{i}\right)^{T}{\boldsymbol{C}^{-1/2}}\boldsymbol{\xi}_{j}=\frac{1}{K}\tilde{\boldsymbol{\xi}}_{i}^{T}\tilde{\boldsymbol{\xi}}_{j},$
(40)

where $\boldsymbol{\xi}_{i}=(\xi_{i}^{1},\dots,\xi_{i}^{K})$, and
$\tilde{\boldsymbol{\xi}}_{i}:=\boldsymbol{C}^{-1/2}\boldsymbol{\xi}_{i}$.
This holds because the pattern correlation matrix is positive-definitive, thus
the (principal) inverse square-root is well-defined. Now we see that, in order
to recover the kernel (38), we need to apply stimuli which are aligned to the
patterns $\tilde{\boldsymbol{\xi}}$. This modification of the stimuli cannot
be directly learnt by the synaptic matrix $J_{ij}$ and it is not enough to
recover the dreaming prescription, this is because here the spins align
towards the sign of the field but they are not exposed to the pattern
correlation matrix. This is to say that we need to modify the dynamical
equations to allow the neurons to perceive the change in the stimuli from
$\boldsymbol{\xi}\rightarrow\tilde{\boldsymbol{\xi}}$, and we choose to modify
the dynamics of the coupling matrix, as learning, as well as the forgetting of
spurious memories, happens in the synapses time-scale.
We take inspiration by Plakhov’s algorithm [Plakhov and Semenov, 1992] and let
the synaptic matrix learn directly from the field as follows

$\displaystyle\sigma_{i}^{(n+1)}$ $\displaystyle=$
$\displaystyle\tanh\left(\beta\sum_{i\neq
j}^{N}J_{ij}^{(n)}\sigma_{j}^{(n)}+\beta
u\tilde{h}_{i}^{(n)}\right):=\tanh(\beta H_{i}^{(n)}),$ (41) $\displaystyle
J_{ij}^{(n+1)}$ $\displaystyle=$ $\displaystyle
J_{ij}^{(n)}\left(1-\frac{\tau}{\tau^{\prime}}\right)+\frac{\tau}{\tau^{\prime}}\frac{H_{j}^{(n)}}{V^{(n)}}\frac{H_{i}^{(n)}}{V^{(n)}}\tanh\beta.$
(42)

Here, $\tilde{h}_{i}$ is the external stimuli obtained by rotating the
patterns according to the inverse square root of the correlation matrix, while
$H_{i}:=\sum_{i\neq j}^{N}J_{ij}^{(n)}\sigma_{j}^{(n)}+u\tilde{h}_{i}^{(n)}$
is the total field acting on the $i$-th spin. The normalisation factor
$V^{(n)}$ is fixed as follows
$V^{(n)}=\frac{\sqrt{N}}{\left\|\boldsymbol{H}^{(n)}\right\|}$ in such a way
that the rescaled total field has norm
$\left\|\frac{\boldsymbol{H}^{(n)}}{V^{(n)}}\right\|=\sqrt{N}$ that is the
norm of a vector in an $N$-dimensional space whose every entry is either $+1$
or $-1$, and the neural configuration is one such vector.
Similar to the previous approach, one can once again prove the convergence of
the dynamics to

$J_{ij}^{(\infty)}=\sum_{\mu=1}^{K}p_{\mu}\tilde{\xi}^{\mu}_{i}\tilde{\xi}^{\mu}_{j}\tanh{\beta},$

provided we are in our typical learning setting. We remind that in this
setting the neural time scale is well separated from the synaptic one, and we
assume that in the temporal window, regardless of $t$ and $\tau^{\prime}$,
neurons are aligned to each pattern $\boldsymbol{\xi}^{\mu}$ an amount of time
that is proportional to the probabilities $\\{p_{\mu}\\}_{\mu=1,\ldots,K}$. In
particular, if $\beta\gg 1$ and $p_{\mu}=\frac{1}{K}$ for each pattern, we
recover the Hebbian kernel (40). We illustrate this point in our numerical
experiments, presented in Sec. 5.

##  5 Results and numerical evidences

We have analytically characterized the model’s behaviour in the case where the
presented stimuli are non-uniformly distributed, with a certain pattern being
more recurrent, and we also have modified the dynamical equations to converge
to the dreaming recipe. In this section, we illustrate these results with
numerical experiments. In particular, we will show the convergence of the
synaptic matrix $\boldsymbol{J}^{(n)}$ to different Hebbian kernels
$\boldsymbol{J}^{(final)}$, whose explicit expression depends on the scenario
under consideration. Moreover, we will check the robustness of the convergence
by studying the fluctuation amplitude around the expected convergence value
(36), that we recast here as

$\|\boldsymbol{\Delta}\|_{F}=\sqrt{\frac{1}{2}\frac{\tau}{\tau^{\prime}}(1-\|\boldsymbol{J}^{(final)}\|^{2}_{F})+(1-\tanh\beta)^{2}\|\boldsymbol{J}^{(final)}\|_{F}^{2}}\,.$
(43)

As a first experiment, we partition the set of $K$ patterns in two families
$\\{\boldsymbol{\xi}^{1,\mu}\\}_{\mu=1}^{K_{1}}$ and
$\\{\boldsymbol{\xi}^{2,\mu}\\}_{\mu=1}^{K_{2}}$, where $K=K_{1}+K_{2}$. We
proceed by randomly selecting a pattern from the first family with a
probability of $p$, and from the second family with a probability of $1-p$.
Once the class is determined, the pattern is uniformly chosen from within that
class and presented to the network. The pattern is extracted uniformly once
the class is chosen. Since this procedure is semantic (i.e. it only involves
index $\mu$) and the Hebbian schemes do not allow temporal correlation between
the patterns (i.e. there is no correlation between patterns with different
indices), it is quite natural to expect that the whole training procedure
converges toward a matrix $J^{(final)}$ that takes the form of mixed Hebb
interaction matrix

$J_{ij}^{(mixed)}:=pJ_{ij}^{(1)}+(1-p)J_{ij}^{(2)},$ (44)

where

$J^{(a)}_{ij}:=\frac{1}{K_{a}}\sum_{\mu=1}^{K_{a}}\xi^{a,\mu}_{i}\xi^{a,\mu}_{j},$

for $a\in\\{1,2\\}$.
In order to study the convergence of the coupling matrix to the prescription
(44), we first analyze the evolution of the synaptic matrix
$\boldsymbol{J}^{(n)}$ by computing, for each iteration, the normalized
Frobenius norm of the difference between the synaptic matrix obtained from the
numerical solution of equations (18)-(19) and the corresponding mixed Hebbian
kernel $\mathbf{J}^{(mixed)}$. The convergence of the coupling matrix is
reported in the left panel of Fig.1 where learning performances of the
dynamics (18)-(19) is shown for different choices of the ratio
$\tau/\tau^{\prime}$. More precisely, the plot is obtained by selecting the
pattern family according to a Bernoulli experiment with a fixed parameter $p$,
subsequently presenting to the network a uniformly selected pattern from the
chosen family. As expected, the different choice of the ratio
$\tau/\tau^{\prime}$ affects both the convergence rate of the algorithm and
the amplitude of the fluctuations around the mean values of the normalized
Frobenius distance Eq. 37 as $t\to\infty$. It is then clear that the closer
the timescales are, the less accurate the retrieval will be. Only when the
timescales are completely separated, the normalized Frobenius distance
converges to zero (see Eq. (43)). The convergence of the coupling matrix to
the prescription (44) is reported in Fig. 1 (left panel) for various values of
the ratio $\tau/\tau^{\prime}$, and for families with the same number of
patterns, i.e. $K_{1}=K_{2}=K/2$. In Fig. 1 (right panel) we report the
retrieval capabilities of the resulting model. To test the retrieval
capabilities of the network, for each pattern, we sample $M$ noisy version of
it by randomly flipping a percentage of bits equal to $(1-r)/2$, with
$r\in(0,1]$; in this context, $r$ can be interpreted as the quality of the
noisy pattern generated. These examples are used as initial condition for the
dynamics. We let the spins evolve according to Eq. (18) for $5$ iterations,
with no stimuli (i.e $h_{i}=0$) and with the synaptic matrix fixed and equal
to $\boldsymbol{J}=\boldsymbol{J}^{mixed}$.The final neural configuration
$\boldsymbol{\sigma}^{(\infty)}$ depends on the initial configuration which,
in turns, depends on both the original pattern $\boldsymbol{\xi}^{a,\mu}$ and
on the realization of the noise injection process that we just describes and
that creates $M$ noisy patterns for each $\mu\in\\{1,\dots,K_{a}\\}$ that we
can label with the new subscript $\bar{a}$. Thus,
$\boldsymbol{\sigma}^{(\infty)}\equiv\boldsymbol{\sigma}^{(\infty)}_{\bar{a}}(\boldsymbol{\xi}^{a,\mu})$
with $\bar{a}\in\\{1,\dots,M\\},a\in\\{1,2\\},\mu\in\\{1,\dots,K_{a}\\}$. As a
performance metric we use Mattis magnetization [Camilli et al., 2022] which is
defined as follows

$m(\boldsymbol{\sigma}^{1},\boldsymbol{\sigma}^{2}):=\frac{\sum_{i=1}^{N}\sigma^{1}_{i}\sigma^{2}_{i}}{N}$
(45)

where $\boldsymbol{\sigma}^{1},\boldsymbol{\sigma}^{2}\in\Sigma_{N}$. We
observe that $m\in[-1,1]$ and if the two configurations are equal then their
Mattis magnetization is equal to $1$ on the contrary if they are orthogonal
then the Mattis magnetization is zero. We evaluate the Mattis magnetization
between each of the $M$ final neural configurations and the corresponding
pattern and then we take the average over the index $\bar{a}$ to get
$m_{\mu,a}=\frac{1}{M}\sum_{\bar{a}=1}^{M}m(\boldsymbol{\sigma}^{(\infty)}_{\bar{a}}(\boldsymbol{\xi}^{a,\mu}),\boldsymbol{\xi}^{a,\mu}))$.
We repeat this algorithm for every pattern, and we take the average also over
$\mu$ index inside each family to get
$\bar{m}_{a}=\frac{1}{K_{i}}\sum_{\mu=1}^{K_{i}}m_{\mu,a}$, for $a=1,2$. The
Mattis magnetizations for the two families are reported as a function of $r$
and for various values of the probability $p$. For $p=1/2$, we recover a
learning kernel which is proportional to the Hebbian kernel of the Hopfield
model, and the magnetizations follow the same behaviour. For $p\neq 1$, we
have a separation between the two families instead. In particular, the higher
$p$, the higher the average Mattis magnetization of the first family of
patterns, meaning that the corresponding patterns are associated with more
stable configurations for the energy function. These results can be easily
extended to the case of linear combinations of $\Omega$ Hebbian kernels with
individual probabilities $p_{a}$:

$J_{ij}^{(final)}=\sum_{a=1}^{\Omega}p_{a}J^{(a)}_{ij},$ (46)

with $\Omega$ being the number of families, and clearly
$p_{1}+\dots+p_{\Omega}=1$.

Figure 1: Results related to the dynamics for binomially extracted patterns.
Left panel: Learning performances of the dynamics given by (18) and (19),
quantified by the normalized Frobenius distance of the synaptic matrix w.r.t.
the mixed Hebb’s kernel (44) (i.e.
$\boldsymbol{J}^{(final)}=\boldsymbol{J}^{(mixed)}$) for
$\tau/\tau^{\prime}=10^{-2},10^{-3},10^{-4}$. The asymptotes correspond to the
values of the normalized Frobenius norm of $\boldsymbol{\Delta}$ given in Eq.
(37) and adapted to the specific case in which the Hebbian kernel is
$\boldsymbol{J}^{mixed}$. Right panel: Retrieval capabilities as a function of
the quality parameter $r$ for $p=0.5,0.6,1$ and $M=100$. The solid curves
refer to the Mattis magnetizations of patterns belonging to the first family,
while dashed ones to the second family. The network parameters are $N=150$,
$u=150$, $\beta=100$, $\alpha=0.05$, for both the panels, with $\alpha$ being
the storage capacity $\alpha=K/N$.

For our second experiment, we explore the setting, where, instead of families
of patterns chosen according to a discrete probability distribution, we have
$K$ patterns that we present to the network with different probabilities
$p_{\mu}$, with $\mu=1,\,...,\,K$. We chose the pattern, that acts as the
external stimulus, based on a scale-free distribution. This choice is
particularly interesting since many real-world distributions adhere to power
laws. In other words, the index $\mu$ is chosen based on the following
probability assignment:

$p_{\mu}=\frac{c_{K,\gamma}}{\mu^{\gamma}},\quad\mu=1,\dots,K,$ (47)

where $c_{K,\gamma}$ is a normalization constant depending on the number of
stored patterns and on the coefficient $\gamma$. This scenario corresponds to
the special case where we have a single pattern per family and $K$ total
families. So by using the formula (46) with
$K_{a}=1,\Omega=K,p_{a}=\frac{c_{K,\gamma}}{a^{\gamma}}$ we get

$J^{(power)}_{ij}=c_{K,\gamma}\sum_{\mu=1}^{K}\xi^{\mu}_{i}\frac{1}{\mu^{\gamma}}\xi^{\mu}_{j}.$
(48)

As in the previous case, we study the convergence of the synaptic matrix
towards the power-law Hebbian kernel by computing at each iteration the
normalized Frobenius norm of the difference between $\boldsymbol{J}$ and
$\boldsymbol{J}^{(power)}$ (see Fig. 2, top panel.

Figure 2: Top panel: Convergence of the synaptic matrix towards the power-law
Hebbian kernel (i.e. $\boldsymbol{J}^{(final)}=\boldsymbol{J}^{(power)}$) for
the ratio $\tau/\tau^{\prime}=10^{-2},10^{-3},10^{-4}$. The curves are
obtained by averaging the results of $1500$ different realizations of the
retrieval process, and the horizontal dashed lines correspond to the expected
normalized Frobenius distance, Eq. 37. Bottom left panel: Mattis magnetization
curves as a function of the exponent $\gamma$ for $q=0.08,0.25,\dots,0.92$.
Bottom right panel: Mattis magnetization as a function of the pattern index
for $\gamma=0,0.5,1,2$. The network parameters are $N=150$, $\alpha=0.05$,
$\beta=100$ for both panels. The network parameters are $N=150$, $u=150$,
$\beta=100$, $\alpha=0.05$, for all the panels, with $\alpha$ being the
storage capacity $\alpha=K/N$.

Then, we analyze the retrieval capabilities of the resulting model as a
function of both the scaling exponent $\gamma$ (Fig. 2 left panel) and the
index $\mu$ related to patterns (Fig. 2 right panel).
We expect that, for a given $\gamma$, there will be some patterns that the
network retrieves correctly and others that are not retrieved. This is because
the retrieval depends on the probability given by Eq. (47): given a pattern
$\boldsymbol{\xi}^{\mu}$, the higher the label $\mu$ the more rarely it will
be presented to the network and the more difficult it will be for the network
to retrieve it. To test our claims we select patterns at random, we sample the
initial neural configurations with quality $r$ w.r.t. the patterns (just as we
did before), we run the neural dynamics given in Eq. (18) with no stimuli and
with the synaptic matrix fixed and equal to Eq. (48) and we plot the
distribution of the resulting magnetizations
$\bar{m}=m(\boldsymbol{\sigma}^{(\infty)}(\boldsymbol{\xi}^{\mu}),\boldsymbol{\xi}^{\mu})$.
We see a peak at $\bar{m}=1$, that corresponds to the patterns that are
correctly retrieved, then we locate the remaining probability mass at
$\bar{m}\gtrsim 0$ because, for the rest of the patterns, the retrieval fails.
The results of this experiment for various $\gamma$ are shown in Fig. 2,
bottom left panel. There, we plot the quantiles $q$ of the distribution of
magnetization values for increasing values of $\gamma$, each of them
corresponds to a different simulation. The so obtained curves start from
$\bar{m}=1$, at $\gamma=0$, because in this case, the uniform distribution is
recovered and the parameters of the network correspond to the retrieval phase
of the Hopfield model with load $\alpha\ll\alpha_{c}$, where
$\alpha_{c}\approx 0.14$ is the critical load of the Hopfield model. Then,
$\bar{m}$ will eventually drop because some patterns are not retrieved. This
drop occurs for lower $\gamma s$ for the lower quantiles because the lower
$\gamma$ is, the more probability mass will be close to $\bar{m}=1$ 444For all
of these experiments, we also varied the network parameters (i.e. the network
size $N$, the storage capacity $\alpha$, the thermal noise $\beta$) both for
the binomially and power-law extracted patterns, and we found that the
training performances are analogous to those of the pure Hebbian case..

Figure 3: Results for the modified synaptic dynamics that leads to the
dreaming prescription. Left panel: normalized Frobenius distance, Eq. 37, of
the synaptic matrix w.r.t. the dreaming kernel for
$\tau/\tau^{\prime}=10^{-2},10^{-3},10^{-4}$. Right Panel: Mean Mattis
magnetizations as a function of the quality parameter $r$ for
$\alpha=0.03,\,0.07,\ldots,\,0.6$. The other network parameters are $N=150$,
$u=150$, and $\beta=100$ for both the panels.

Moving to a numerical analysis of the dreaming prescription, in Fig. 3, top
panel, we plot the distance between the coupling matrix and the dreaming
prescription (40) as a function of the evolution time, for different choices
of the ratio $\tau/\tau^{\prime}$, finding that the algorithm converges, as
expected. We observe that the norm of the difference between the dreaming and
the numerical coupling matrices is given by (37). We also observe that the
dynamics follows (14), as it happened in all the precedent numerical
scenarios. In Fig. 3, bottom left panel, we present the results concerning the
system’s retrieval capabilities after training. This scenario is achieved by
keeping the resultant coupling matrix fixed, and the results are plotted as a
function of the quality $r$ for different values of the storage capacity
$\alpha=K/N$. Evidently, the model can manage a larger information content
compared to the Hebbian case reported in [Agliari et al., 2023]. Additionally,
at a sufficiently low storage capacity, the system handles low-quality
examples more effectively than the standard Hopfield model. It is noteworthy
that when comparing Fig. 3 in [Agliari et al., 2023] (representing the Hebbian
case) with Fig. 3 (bottom-left panel, representing the dreaming kernel), the
latter model exhibits greater stability against perturbations in the initial
conditions for retrieval dynamics, especially at higher values of the quality
parameter $r$.
This implies that the model, when equipped with the dreaming kernel, features
broader attraction basins related to the stored patterns. As a result, the
retrieval quality increases, especially at very high storage capacities where
the Hebbian counterpart fails to achieve effective retrieval. This result is
in perfect agreement with the analysis for the attraction basins performed in
[Fachechi et al., 2019b].
Finally, we emphasize that when the load of the network is below the higher
critical storage capacity, the results are analogous to the pure Hebbian case.
Whilst, above this threshold, the network properly converges to the dreaming
kernel.

##  6 Conclusions

In this paper, we examined the modeling of Pavlov’s Conditioning theory using
a statistical mechanics approach.
We have analytically proven a Pavlovian model for associative memory implies
the Hebbian learning prescription, this statement is the main claim of
[Agliari et al., 2023], and up to this point, was only characterized
numerically. We illustrated the analytical results in the setting in which the
stimuli are divided into two families and the network sees one family of
stimuli more frequently than the other. In such scenarios, the coupling matrix
converges to a weighted element-wise sum of the coupling matrices of the two
families, with the weights determined by the probabilities of selecting each
family.
Finally, in the spirit of describing biological networks within our framework,
we have shown that the model can describe neural networks that sleep.
Pavlovian learning rules are indeed compatible with the dreaming mechanism
used in the associative neural network literature. In this case, the dynamics
of the synapses allow our model’s coupling matrix to converge towards the so-
called dreaming kernel. The network not only handles more patterns compared to
its purely Hebbian counterpart but also exhibits enhanced retrieval quality
and resilience against initial condition perturbations, especially at higher
storage capacities. We have discussed this setting, emphasizing that all
analytical considerations regarding convergence remain valid in this case. In
conclusion, this paper provides a comprehensive analytical treatment of the
learning process in Pavlovian learning machines.

## Acknowledgments

The authors acknowledge Adriano Barra, Elena Agliari, and Alberto Fachechi for
the precious discussions about this work.
D.L. acknowledges INdAM-GNFM and C.N.R. (National Research Council).

##  Appendix A The master equation derivation of the coupled dynamics

A more formal approach to derive the evolution equations (18)-(19) for the
expectation values of $\sigma_{i}$ and $J_{ij}$ over time makes use of the
master equation related to their stochastic process [Agliari et al., 2023]. We
sketch this formulation in this appendix and remand to the paper for the full
formal discussion.
We start by introducing
$w_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J}\to\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime})$,
the transition rate between two system states
$(\boldsymbol{\sigma},\boldsymbol{J})$ and
$(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime})$, and let
$\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)$ denote the
probability of the system being in the configuration
$(\boldsymbol{\sigma},\boldsymbol{J})$ at time $t$.
The master equation governing the time evolution of the probability
$\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)$ is given by:

$\displaystyle\frac{d}{dt}\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)=$
$\displaystyle\sum_{\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}\sum_{\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\Bigg{[}\mathcal{P}_{N,\beta}(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime},t)w_{N,\beta}(\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime}\to\boldsymbol{J},\boldsymbol{\sigma})$
$\displaystyle-
P_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)w_{N,\beta}(\boldsymbol{J},\boldsymbol{\sigma}\to\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime})\Bigg{]}$
(49)

Reference [Agliari et al., 2023] proceeds assuming that neural activities and
synaptic weights are updated simultaneously to factorize the transition rate
$w(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime}\to\boldsymbol{\sigma},\boldsymbol{J})$
as

$w_{N,\beta}(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime}\to\boldsymbol{\sigma},\boldsymbol{J})=\prod_{l,m}\frac{1}{2}\left[1+\tanh\left(\beta
J_{lm}\sigma_{l}^{\prime}\sigma_{m}^{\prime}\right)\right]\prod_{i}\frac{1}{2}\left[1+\sigma_{i}\tanh\left(\beta\sum_{j}J_{ij}^{\prime}\sigma_{j}^{\prime}+\beta
h_{i}\right)\right].$ (50)

The transition rate is thus expressed as a product of single-variable
conditional probabilities.

Focusing on neural dynamics first, the mean neural activity is defined as

$\langle\sigma_{i}\rangle=\sum_{\boldsymbol{\sigma}\in\Sigma_{N}}\sum_{\boldsymbol{J}\in\mathcal{J}_{N}}\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)\,\sigma_{i},\leavevmode\nobreak\
\leavevmode\nobreak\ \textrm{for}\leavevmode\nobreak\ \leavevmode\nobreak\
i=1,...,N$ (51)

and set the typical timescale for its dynamics to be $\tau$.
Multiplying the master equation (49) by $\sigma_{i}$, and summing over all
possible neural and synaptic configurations, we arrive at the following
stochastic process:

$\displaystyle\tau\frac{d\langle\sigma_{i}\rangle}{dt}=$ $\displaystyle+$
$\displaystyle\sum_{\boldsymbol{J},\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}\sigma_{i}[\mathcal{P}_{N,\beta}(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime},t)w_{N,\beta}(\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime}\to\boldsymbol{J},\boldsymbol{\sigma})]$
$\displaystyle-$
$\displaystyle\sum_{\boldsymbol{J},\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}\sigma_{i}[\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)w_{N,\beta}(\boldsymbol{J},\boldsymbol{\sigma}\to\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime})].$

It is observed that, by definition,

$\sum_{\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}\sum_{\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}w_{N,\beta}(\boldsymbol{J},\boldsymbol{\sigma}\to\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime})=1.$
(53)

This identity can be used to simplify (A),

$\tau\frac{d\langle\sigma_{i}\rangle}{dt}=-\langle\sigma_{i}\rangle+\sum_{\boldsymbol{J},\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}\sigma_{i}\mathcal{P}_{N,\beta}(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime},t)w_{N,\beta}(\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime}\to\boldsymbol{J},\boldsymbol{\sigma})$
(54)

which, in turn, is equivalent to

$\tau\frac{d\langle\sigma_{i}\rangle}{dt}=-\langle\sigma_{i}\rangle+\left\langle\tanh\left(\beta\sum_{j}J_{ij}\sigma_{j}+\beta
h_{i}\right)\right\rangle.$ (55)

By applying the mean field approximation, equation (55) can be recast as

$\tau\frac{d\langle\sigma_{i}\rangle}{dt}=-\langle\sigma_{i}\rangle+\tanh\Big{(}\beta\sum_{j}\Big{\langle}J_{ij}\Big{\rangle}\langle\sigma_{j}\rangle+\beta
h_{i}\Big{)}.$ (56)

It is noted that, while at the individual neuron level, neural firing is
typically described by Langevin equations [Tuckwell, 1988], at this coarse-
grained treatment, the evolution of the mean firing rate is governed by an
ODE: the random terms are averaged out, but the noise $\beta$ is still a
central quantity. Indeed, for $\beta\to 0$ neither the post-synaptic
potentials nor the external stimuli can be perceived, and it is the term
$-\langle\sigma_{i}$ in (56) that ensures the mean firing rate to return to
zero.

Moving to the synaptic dynamics, its typical timescale can be assumed to be
$\tau^{\prime}>\tau$, and the dynamical evolution of the average of synaptic
weights, defined by

$\langle
J_{ij}\rangle=\sum_{\boldsymbol{\sigma}\in\Sigma_{N}}\sum_{\boldsymbol{J}\in\mathcal{J}_{N}}\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)J_{ij}$

is obtained by multiplying the master equation by $J_{{}_{ij}}$ and summing
over all possible neural configurations and synaptic matrices.

$\displaystyle\tau^{\prime}\frac{d\langle J_{ij}\rangle}{dt}=$
$\displaystyle+\sum_{\boldsymbol{J},\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}J_{ij}\left[\mathcal{P}_{N,\beta}(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime},t)w_{N,\beta}(\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime}\to\boldsymbol{J},\boldsymbol{\sigma})\right]$
$\displaystyle-\sum_{\boldsymbol{J},\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma},\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}J_{ij}\left[\mathcal{P}_{N,\beta}(\boldsymbol{\sigma},\boldsymbol{J},t)w_{N,\beta}(\boldsymbol{J},\boldsymbol{\sigma}\to\boldsymbol{J}^{\prime},\boldsymbol{\sigma}^{\prime})\right],$
(57)

which simplifies to

$\tau^{\prime}\frac{d\langle J_{ij}\rangle}{dt}=-\langle
J_{ij}\rangle+\sum_{\boldsymbol{J}^{\prime}\in\mathcal{J}_{N}}\sum_{\boldsymbol{\sigma}^{\prime}\in\Sigma_{N}}\mathcal{P}_{N,\beta}(\boldsymbol{\sigma}^{\prime},\boldsymbol{J}^{\prime},t)\tanh\left(\beta\right)\sigma_{i}^{\prime}\sigma_{j}^{\prime},$

which under the mean field approximation, leads to

$\tau^{\prime}\frac{d\langle J_{ij}\rangle}{dt}=-\langle
J_{ij}\rangle+\tanh\left(\beta\right)\langle\sigma_{i}\rangle\langle\sigma_{j}\rangle.$
(58)

Equations (9) and (58) constitute the system of coupled ODEs, and are
equivalent to (18) and (19).

##  Appendix B Analytical derivation of the Hebbian kernel from Pavlovian
dynamics

In Section 3.2 of the main text we proved the convergence of the discrete
dynamics to the Hebbian kernel in the simple case where the neurons are
clamped to a fixed pattern. Here we will analyze the more realistic case in
which we let the neurons change. In this appendix, we focus on the unrealistic
case where the patterns are presented cyclically to the network.
In the case of two patterns the neurons are clamped to $\boldsymbol{\xi}^{1}$
when the iteration step is even and to $\boldsymbol{\xi}^{2}$ when the
iteration step is odd.

To simplify the notation, let us denote by
$c_{1}:=\xi_{i}^{1}\xi_{j}^{1}\tanh\beta$,
$c_{2}:=\xi_{i}^{2}\xi_{j}^{2}\tanh\beta$ and
$\epsilon:=\frac{\tau}{\tau^{\prime}}$. The dynamics of the coupling matrix
can then be written as

$\begin{cases}J_{ij}^{(n+1)}=(1-\epsilon)J_{ij}^{(n)}+\epsilon
c_{1}\quad\text{if }n\equiv 0\pmod{2}\\\
J_{ij}^{(n+1)}=(1-\epsilon)J_{ij}^{(n)}+\epsilon c_{2}\quad\text{if }n\equiv
1\pmod{2}\,.\end{cases}$ (59)

Following the same approach used in Section 3.2, let us express
$J_{ij}^{(n+1)}$ in terms of its previous values while considering the
alternating influence of $c_{1}$ and $c_{2}$. To simplify, we consider the
sequence of updates as pairs, each pair covering a full cycle of the two
patterns. Starting from an odd iteration step,

$\displaystyle J_{ij}^{(n+2)}$
$\displaystyle=(1-\epsilon)[(1-\epsilon)J_{ij}^{(n)}+\epsilon c_{2}]+\epsilon
c_{1}$
$\displaystyle=(1-\epsilon)^{2}J_{ij}^{(n)}+\epsilon[(1-\epsilon)c_{2}+c_{1}].$

After $2n$ steps (or $n$ cycles), the coupling matrix evolves as:

$J_{ij}^{(2n)}=(1-\epsilon)^{2n}J_{ij}^{(0)}+\epsilon\left[c_{1}+(1-\epsilon)c_{2}\right]\sum_{k=0}^{n-1}(1-\epsilon)^{2k}.$
(60)

The sum in the last term is a geometric series with ratio $(1-\epsilon)^{2}$,
which converges, as $n$ approaches infinity to
$\frac{1}{2\epsilon-\epsilon^{2}}$. Therefore we get

$\lim_{n\to\infty}J_{ij}^{(2n)}=\frac{1}{2-\epsilon}\left[c_{1}+(1-\epsilon)c_{2}\right]$
(61)

that, when the timescales are well-separated (i.e. $\epsilon\rightarrow 0$),
is exactly mixed Hebbian kernel defined in (22)

Similarly, we compute the limiting coupling matrix, but starting, from an even
iteration step. By denoting with $J_{ij,\texttt{even}}^{(\infty)}$ the “even
update” convergence matrix and with $J_{ij,\texttt{odd}}^{(\infty)}$ the “odd
update” convergence matrix, we obtain

$\displaystyle\lim_{n\to\infty}J_{ij}^{(n)}$
$\displaystyle=\frac{1}{2}\left(J_{ij,\texttt{even}}^{(\infty)}+J_{ij,\texttt{odd}}^{(\infty)}\right)$
$\displaystyle=\frac{1}{2}\frac{1}{2-\epsilon}\left[c_{1}+(1-\epsilon)c_{2}+c_{2}+(1-\epsilon)c_{1}\right]$
$\displaystyle=\frac{1}{2}\frac{1}{2-\epsilon}\left[c_{1}\left(2-\epsilon\right)+c_{2}\left(2-\epsilon\right)\right]$
$\displaystyle=\frac{1}{2}\left(c_{1}+c_{2}\right).$

We are now ready to extend this conclusion to an arbitrary number of patterns.
Let’s consider a scenario where we have $K$ patterns, denoted by
$\\{\boldsymbol{\xi}^{\mu}\\}_{\mu=1,\dots K}$, and the system alternates
among these patterns in a sequential manner over time.

For each pattern $\mu$, let $c_{\mu}:=\xi_{i}^{\mu}\xi_{j}^{\mu}\tanh\beta$.
The dynamics of the coupling matrix when the system is exposed to these
patterns in sequence can be generalized as follows:

$J_{ij}^{(n+1)}=\begin{cases}(1-\epsilon)J_{ij}^{(n)}+\epsilon
c_{1}&\quad\text{if }n\equiv 0\pmod{K}\\\ (1-\epsilon)J_{ij}^{(n)}+\epsilon
c_{2}&\quad\text{if }n\equiv 1\pmod{K}\\\ \quad\vdots\\\
(1-\epsilon)J_{ij}^{(n)}+\epsilon c_{M}&\quad\text{if }n\equiv
K-1\pmod{K}.\end{cases}$ (62)

Following a similar approach as with two patterns, we consider the updates
over $K$ cycles to capture the influence of all patterns. After $Kn$ steps (or
$n$ complete cycles through all $K$ patterns), the coupling matrix evolves as:

$J_{ij}^{(Kn)}=(1-\epsilon)^{Kn}J_{ij}^{(0)}+\epsilon\left[\sum_{\mu=1}^{K}c_{\mu}\left(1-\epsilon\right)^{\mu-1}\right]\sum_{\eta=0}^{n-1}(1-\epsilon)^{K\eta}.$
(63)

The sum in the last term, $\sum_{\eta=0}^{n-1}(1-\epsilon)^{K\eta}$, is a
geometric series with ratio $(1-\epsilon)^{K}$ which converges, as $n$
approaches infinity, to $\frac{1}{1-(1-\epsilon)^{K}}$. Therefore we can write

$\lim_{n\to\infty}J_{ij}^{(Kn)}=\frac{\epsilon}{1-(1-\epsilon)^{K}}\sum_{\mu=1}^{K}c_{\mu}\left(1-\epsilon\right)^{\mu-1}=\frac{1}{K}\sum_{\mu=1}^{K}\xi_{i}^{\mu}\xi_{j}^{\mu}\tanh\beta,$
(64)

where in the last equality we used the fact that the first-order Taylor
expansion for $\frac{\epsilon}{1-(1-\epsilon)^{K}}$ simplifies to the
normalization factor $\frac{1}{K}$, and that neglecting terms with the factor
$\epsilon$ in the round brackets of the sum, we can further simplify this
expression in the limit as $\epsilon\rightarrow 0$, which corresponds to well-
separated timescales.

More rigorously, we can sum over all the possible $K$ cycle realizations, each
of which starts from a different pattern $\boldsymbol{\xi}^{\mu}$ and finishes
with the same pattern. We can compute

$\lim_{n\to\infty}J_{ij}^{(Kn)}=\frac{1}{K}\frac{\epsilon}{1-(1-\epsilon)^{K}}\sum_{\mu=1}^{K}c_{\mu}\left[\sum_{\eta=1}^{K}\left(1-\epsilon\right)^{\eta-1}\right]=\frac{1}{K}\frac{\epsilon}{1-(1-\epsilon)^{K}}\sum_{\mu=1}^{K}c_{\mu}\frac{1-(1-\epsilon)^{K}}{1-(1-\epsilon)}$
(65)

to reveal that Eq.(64) is an exact result.
Thus, we have demonstrated that, in the infinite time limit and for a system
that sequentially alternates among $M$ patterns, the coupling matrix converges
to the Hebbian prescription.

## References

  * [Agliari et al., 2019] Agliari, E., Alemanno, F., Barra, A., and Fachechi, A. (2019).  Dreaming neural networks: Rigorous results.  Journal of Statistical Mechanics: Theory and Experiment, 2019(8).
  * [Agliari et al., 2023] Agliari, E., Aquaro, M., Barra, A., Fachechi, A., and Marullo, C. (2023).  From Pavlov conditioning to Hebb Learning.  Neural Computation, 35(5):930–957.
  * [Amit et al., 1985] Amit, D., Gutfreund, H., and Sompolinsky, H. (1985).  Storing infinite numbers of patterns in a spin-glass model of neural networks.  Physical Review Letters, 14(55):123304.
  * [Amit, 1989] Amit, D. J. (1989).  Modeling brain function: The world of attractor neural networks.  Cambridge university press.
  * [Aquaro et al., 2022] Aquaro, M., Alemanno, F., Kanter, I., Durante, F., Agliari, E., and Barra, A. (2022).  Recurrent neural networks that generalize from examples and optimize by dreaming.  arXiv.
  * [Camilli et al., 2022] Camilli, F., Contucci, P., and Mingione, E. (2022).  An inference problem in a mismatched setting: a spin-glass model with mattis interaction.  SciPost Physics, 12(4):125.
  * [Coolen et al., 2005] Coolen, A. C. C., Kühn, R., and Sollich, P. (2005).  Theory of Neural Information Processing Systems.  Oxford University Press.
  * [Cooper et al., 2013] Cooper, R. P., Cook, R., Dickinson, A., and Heyes, C. M. (2013).  Associative (not Hebbian) learning and the mirror neuron system.  Neuroscience Letters, 540:28–36.
  * [Crick and Mitchison, 1983] Crick, F. and Mitchison, G. (1983).  The function of dream sleep.  Nature, 304(5922):111–114.
  * [Fachechi et al., 2019a] Fachechi, A., Agliari, E., and Barra, A. (2019a).  Dreaming neural networks: Forgetting spurious memories and reinforcing pure ones.  Neural Networks, 112:24–40.
  * [Fachechi et al., 2019b] Fachechi, A., Agliari, E., and Barra, A. (2019b).  Dreaming neural networks: Forgetting spurious memories and reinforcing pure ones.  Neural Networks, 112:24–40.
  * [Gardner and Derrida, 1988] Gardner, E. and Derrida, B. (1988).  Optimal storage properties of neural network models.  J. Phys. A, 21:271.
  * [Hassabis et al., 2017] Hassabis, D., Kumaran, D., Summerfield, C., and Botvinick, M. (2017).  Neuroscience-Inspired Artificial Intelligence.  Neuron, 95(2):245–258.
  * [Hebb, 1949] Hebb, D. (1949).  The organization of behavior.  Wiley Press.
  * [Hemmen et al., 1990] Hemmen, J. V., Ioffe, L., Kühn, R., and Vaas, M. (1990).  Increasing the efficiency of a neural network through unlearning.  Physica A, 163:386–392.
  * [Hopfield, 1982] Hopfield, J. J. (1982).  Neural networks and physical systems with emergent collective computational abilities.  Proceedings of the national academy of sciences, 79(8):2554–2558.
  * [Kanter and Sompolinsky, 1987a] Kanter, I. and Sompolinsky, H. (1987a).  Associative recall of memory without errors.  Phys. Rev. A, 35:380.
  * [Kanter and Sompolinsky, 1987b] Kanter, I. and Sompolinsky, H. (1987b).  Associative recall of memory without errors.  Physical Review A, 35(1):380.
  * [Keysers and Gazzola, 2014] Keysers, C. and Gazzola, V. (2014).  Hebbian learning and predictive mirror neurons for actions, sensations and emotions.  Philosophical Transactions of the Royal Society B: Biological Sciences, 369(1644).
  * [Kohonen, 1984] Kohonen, T. (1984).  Self organization and associative memory.  Springer-Verlag.
  * [Maass and Zador, 1997] Maass, W. and Zador, A. (1997).  Dynamic stochastic synapses as computational units.  Advances in neural information processing systems, 10.
  * [Marc et al., 1987] Marc, M., Giorgio, P., and Virasoro, M. A. (1987).  Spin glass theory and beyond.  World Scientific Publishing.
  * [Personnaz et al., 1985] Personnaz, L., Guyon, I., and Dreyfus, G. (1985).  Information storage and retrieval in spin-glass like neural networks.  Journal of Physics Letters, 46:L–359–365.
  * [Plakhov and Semenov, 1992] Plakhov, A. and Semenov, S. (1992).  The modified unlearning procedure for enhancing storage capacity in hopfield network.  IEEE Transactions, page 242.
  * [Shatz, 1992] Shatz, C. J. (1992).  The developing brain.  Scientific American, 267:60–67.
  * [Tuckwell, 1988] Tuckwell, H. C. (1988).  Introduction to theoretical neurobiology: linear cable theory and dendritic structure, volume 1.  Cambridge University Press.