Learning and organization of memory for evolving patterns

A network $J$ (with elements $J_{ij}$) is presented with $N$ random (orthogonal) patterns $\ket{\sigma^{\alpha}}$ (with $\alpha=1,\dots N$), with entries $\sigma^{\alpha}_{i}\in\{-1,1\}$, reflecting the $N$ pattern classes. For a network with $C$ compartments (with $1\leq C\leq N$), we initialize each sub-network $J^{s}$ at time $t_{0}$ as $J^{s}_{i,j}(t_{0})=\frac{1}{N/C}\sum_{\alpha\in\mathcal{A}_{s}}\sigma^{\alpha}_{i}\sigma^{\alpha}_{j}$ and $J^{s}_{ii}(t_{0})=0$; here, $\mathcal{A}_{s}$ is a set of $N/C$ randomly chosen (without replacement) patterns initially assigned to the compartment (sub-network) $s$. We then let the network undergo an initial learning process. At each step an arbitrary pattern $\sigma^{\nu}$ is presented to the network and a sub-network $J^{s}$ is chosen for an update with a probability

$\displaystyle P_{s}=\frac{\exp\left[{-\beta_{\mathrm{S}}E\left(J^{s}(t),\sigma^{\nu}(t)\right)}\right]}{\sum_{r=1}^{C}\exp\left[{-\beta_{\mathrm{S}}E\left(J^{s}(t),\sigma^{\nu}(t)\right)}\right]},$

(S5)

where the energy is defined as,

$$E\left(J^{s}(t),\sigma^{\nu}(t)\right)=\frac{-1}{2L}\sum_{i,j}J^{s}_{i,j}(t)\sigma^{\nu}_{i}(t)\sigma^{\nu}_{j}(t)\equiv\frac{-1}{2}\braket{\sigma^{\nu}(t)}{J^{s}(t)}{\sigma^{\nu}(t)}$$

(S6)

and $\beta_{\mathrm{S}}$ is the inverse temperature associated with choosing the right compartment. We then update the selected sub-network $J^{s}$, using the Hebbian update rule,

$$J^{s}_{i,j}(t+1)=\begin{cases}(1-\lambda)\,J^{s}_{i,j}(t)+\lambda\,\sigma^{\nu}_{i}\sigma^{\nu}_{j},&\text{if }i\neq j;\\ 0,&\text{otherwise.}\end{cases}$$

(S7)

For dynamic patterns, the presented patterns undergo evolution with mutation rate $\mu$, which reflects the average number of spin flips in a given pattern per network’s update event (Fig. 1).

Our goal is to study the memory retrieval problem in a network that has reached its steady state. The state of a network $J(t_{n})$ at the time step $n$ can be traced back to the initial state $J(t_{0})$ as,

$$J(t_{n})=(1-\lambda)^{n}J(t_{0})+\lambda\sum_{i=1}^{n}(1-\lambda)^{n-i}\ket{\sigma(t_{i})}\bra{\sigma(t_{i})}$$

(S8)

The contribution of the initial state $J(t_{0})$ to the state of the network at time $t_{n}$ decays as $(1-\lambda)^{n}$ (eq. S8). Therefore, we choose the number of steps to reach the steady state as $n_{\rm stat.}=\max\left[10N,2C\,\operatorname{ceil}\left(\frac{\log 10^{-5}}{\log(1-\lambda)}\right)\right]$. This criteria ensures that $(1-\lambda)^{n_{\rm stat.}}\leq 10^{5}$ and the memory of the initial state $J(t_{0})$ is removed from the network $J(t)$. We will then use this updated network to collect the statistics for memory retrieval. To report a specific quantity from the network (e.g., the energy), we pool the $n_{\rm stat.}$ samples collected from each of the 50 realizations.

Once the trained network approaches a stationary state, we collect the statistics of the stored memory.

To find a memory attractor $\sigma^{\alpha}_{\text{att}}$ for a given pattern $\sigma^{\alpha}$ we use a Metropolis algorithm in the energy landscape $E(J^{s},\sigma^{\alpha})$ (eq. S6). To do so, we make spin-flips in a presented pattern $\sigma^{\alpha}\rightarrow\tilde{\sigma}^{\alpha}$ and accept a spin-flip with probability

$\displaystyle P(\sigma^{\alpha}\rightarrow\tilde{\sigma}^{\alpha})=\min\left(1,e^{-\beta_{\mathrm{H}}\Delta E}\right),$

(S9)

where $\Delta E=E(J^{s},\tilde{\sigma}^{\alpha})-E(J^{s},\sigma^{\alpha})$ and $\beta_{\mathrm{H}}$ is the inverse (Hopfield) temperature for pattern retrieval in the network (see Fig. 1). We repeat this step for $2\times 10^{6}$ steps, which is sufficient to find a minimum of the landscape (see Fig. S3).

For systems with more than one compartment $C$, we first choose a compartment according to eq. S5, and then perform the Metropolis algorithm within the associated compartment.

After finding the energy minima, we update the systems for $n^{\prime}_{\rm stat.}=\max[2\cdot 10^{3},n_{\rm stat.}]$ steps. At each step we present patterns as described above and collect the statistics of the recognition energy $E(J^{s}(t),\sigma^{\alpha}(t))$ between a presented pattern $\sigma^{\alpha}$ and the memory compartment $J^{s}(t)$, assigned according to eq. S5. These measurements are used to describe the energy statistics (Figs. 2,S2) of the patterns and the accuracy of pattern-compartment association (Fig. 4B,D). After the $n^{\prime}_{\rm stat.}$ steps, we again use the Metropolis algorithm to find the memory attractors associated with the presented patterns. We repeat this analysis for $50$ independent realizations of the initializing pattern classes $\{\sigma^{\alpha}(t_{0})\}$, for each set of parameters $\{L,N,C,\lambda,\mu,\beta_{\mathrm{S}},\beta_{\mathrm{H}}\}$.

When calculating the mean performance ${\mathcal{Q}}$ of a strategy (see Figs. 2,3,4,S7), we set the overlap between attractor and pattern $q^{\alpha}=|\braket{\sigma^{\alpha}_{\text{att}}}{\sigma^{\alpha}}|$ equal to zero when patterns are not recognized ($q^{\alpha}<0.8$). As a result, systems can only achieve a non-zero performance if they recognize some of the patterns. This choice eliminates the finite size effect of a random overlap $\sim 1/\sqrt{L}$ between an attractor and a pattern (see Fig. S3). This correction is especially important when comparing systems with different sub-network sizes ($L_{\mathrm{c}}\equiv L/C$) in the $\beta_{\mathrm{H}}<1$ regime (Figs. 4,S7), where random overlaps for small $L_{\mathrm{c}}$ could otherwise result in a larger mean performance compared to larger systems that correctly reconstruct a fraction of the patterns.

We use the mutual information ${\text{MI}}(\Sigma,\mathcal{C})$ between the set of pattern classes $\Sigma\equiv\{\sigma^{\alpha}\}$ and the set of all compartments $\mathcal{C}$ to quantify the accuracy in associating a presented pattern with the correct compartment,

$$\begin{split}{\text{MI}}(\Sigma,\mathcal{C})&=H(\mathcal{C})-H(\mathcal{C}|\Sigma)\\ &=-\sum_{c\in\mathcal{C}}P(c)\log P(c)-\left[-\sum_{\sigma^{\alpha}\in\Sigma}P(\sigma^{\alpha})\sum_{c\in\mathcal{C}}P(c|\sigma^{\alpha})\log P(c|\sigma^{\alpha})\right].\end{split}$$

(S10)

Here $H(\mathcal{C})$ and $H(\mathcal{C}|\Sigma)$ are the entropy of the compartments, and the conditional entropy of the compartments given the presented patterns, respectively. If chosen randomly, the entropy associated with choosing a compartment is $H_{\text{ random}}(\mathcal{C})=\log C$. The mutual information (eq. S10) measures the reduction in the entropy of compartments due to the association between the patterns and the compartments, measured by the conditional entropy $H(\mathcal{C}|\Sigma)$. Figure 4B,D shows the mutual information ${\text{MI}}(\Sigma,\mathcal{C})$ scaled by its upper bound $H(\mathcal{C})$, in order to compare networks with a different number of compartments.

The optimal learning rate is determined by maximizing the network’s performance ${\mathcal{Q}}(\lambda)$ (eq. 2) against evolving patterns with a specified mutation rate:

$$\lambda^{*}=\operatorname*{argmax}_{\lambda}{{\mathcal{Q}}}(\lambda)$$

(S11)

We can numerically estimate the optimal learning rate as defined eq. S11; see Figs. 2,3. However, the exact analytical evaluation of the optimal learning rate is difficult and we use an approximate approach and find the learning rate that minimizes the expected energy of the patterns in the stationary state $\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]$, assuming that patterns are shown to the network at a fixed order. Here, the subscripts explicitly indicate the learning rate of the network $\lambda$, and the evolutionary overlap of the pattern $\rho$. To evaluate an analytical approximation for the energy, we first evaluate the state of the network $J(t)$ at time step $t$, given all the prior encounters of the networks with patterns shown at a fixed order.

	$\displaystyle\frac{1}{L}\left(J(t)+\mathds{1}\right)$	$\displaystyle=\lambda\sum_{j=1}^{\infty}(1-\lambda)^{(j-1)}\ket{\sigma(t-j)}\bra{\sigma(t-j)}$		(S12)
		$\displaystyle=\lambda\underbrace{\sum_{i=1}^{\infty}(1-\lambda)^{(i-1)N}\underbrace{\sum_{\alpha=1}^{N}(1-\lambda)^{\alpha-1}\ket{\sigma^{\alpha}(t-\alpha-(i-1)N)}\bra{\sigma^{\alpha}(t-\alpha-(i-1)N)}}_{\text{sum over {\it N} pattern classes}}}_{\text{sum over time (generations, }i)}$		(S13)
		$\displaystyle=\lambda\underbrace{\sum_{\alpha=1}^{N}\underbrace{\sum_{i=0}^{\infty}(1-\lambda)^{(\alpha-1)+iN}\ket{\sigma^{\alpha}(t-\alpha-iN)}\bra{\sigma^{\alpha}(t-\alpha-iN)}}_{\text{sum over time}}}_{\text{sum over patterns}}.$		(S14)

Here, we referred to the (normalized) pattern vector from the class $\alpha$ presented to the network at time step $t$ by $\ket{\sigma^{\alpha}(t)}\equiv\frac{1}{\sqrt{L}}\sigma^{\alpha}(t)$. Without a loss of generality, we assumed that the last pattern presented to the network at time step $t-1$ is from the first pattern class $\ket{\sigma^{1}(t-1)}$, which enabled us to split the sum in eq. S12 into two separate summations over pattern classes and $N$ time-steps generations (eq. S13). Adding the identity matrix $\mathds{1}$ on the left-hand side of eq. S12 assures that the diagonal elements vanish, as defined in eq. S7.

The mean energy of the patterns, which in our setup is the energy of the pattern from the $N^{th}$ class at time $t$, follows

$$\begin{split}\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]&=\mathbb{E}\left[{-\frac{1}{2}\braket{\sigma^{N}(t)}{J(t)}{\sigma^{N}(t)}}\right]\\ &=\mathbb{E}\left[{-\frac{L-1}{2}\lambda\sum_{\alpha=1}^{N}\sum_{i=0}^{\infty}(1-\lambda)^{(\alpha-1)+iN}\braket{\sigma^{N}(t)}{\sigma^{\alpha}(t-\alpha-iN)}\braket{\sigma^{\alpha}(t-\alpha-iN)}{\sigma^{N}(t)}}\right].\end{split}$$

(S15)

Since the pattern families are orthogonal to each other, we can express the overlap between patterns at different times as $\braket{\sigma^{\alpha}(t_{1})}{\sigma^{\nu}(t_{2})}=\delta_{\alpha,\nu}(1-2\mu)^{|t_{2}-t_{1}|}\equiv\delta_{\alpha,\nu}\rho^{|t_{2}-t_{1}|}$, and simplify the energy function in eq. LABEL:eq:energypattern,

$$\begin{split}\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]&=-\frac{L-1}{2}\lambda\sum_{i=0}^{\infty}(1-\lambda)^{(N-1)+iN}\rho^{2(N+iN)}\\ &=-\frac{L-1}{2}\lambda(1-\lambda)^{(N-1)}\rho^{2N}\sum_{i=0}^{\infty}\left((1-\lambda)^{N}\rho^{2N}\right)^{i}\\ &=-\frac{L-1}{2}\,\lambda\,\frac{(1-\lambda)^{(N-1)}\rho^{2N}}{1-(1-\lambda)^{N}\rho^{2N}}.\end{split}$$

(S16)

Since accurate pattern retrieval depends on the depth of the energy valley for the associative memory, we will use the expected energy of the patterns as a proxy for the performance of the network. We can find the approximate optimal learning rate that minimized the expected energy by setting $\partial\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]/{\partial\lambda}=0$, which results in

$$\begin{split}(1-2\mu)^{2N}=(1-\lambda^{*})^{-N}(1-N\lambda^{*})&\Longrightarrow 1-4N\mu+\mathcal{O}(\mu^{2})=1+\frac{1}{2}(N-N^{2})(\lambda^{*})^{2}+\mathcal{O}(\lambda^{3});\\ &\Longrightarrow\lambda^{*}(\mu)\simeq\sqrt{8\mu/(N-1)}.\end{split}$$

(S17)

where we used the fact that both the mutation rate $\mu$ and the learning rate $\lambda$ are small, and therefore, expanded eq. S17 up to the leading orders in these parameters.

In addition, eq. S17 establishes an upper bound for the learning rate: $\lambda<\frac{1}{N}$. Therefore, our expansion in mutation rate (eq. S17) is only valid for $8\mu<\frac{1}{N}$, or equivalently for $\mu_{\text{eff}}=N\mu<12.5\%$; the rates used in our analyses lie far below these upper bounds.

The memory attractors associated with a given pattern class can lag behind and only reflect the older patterns presented prior to the most recent encounter of the network with the specified class. As a result, the upper bound for the performance of a network ${\mathcal{Q}}_{{}_{\text{lag}}}=\rho^{g_{{}_{\text{lag}}}N}\approx 1-2g_{{}_{\text{lag}}}\mu_{{\text{eff}}}$ is determined by both the evolutionary divergence of patterns between two encounters $\mu_{\text{eff}}$ and number of generations $g_{{}_{\text{lag}}}$ by which the stored memory lags behind; we measure $g_{{}_{\text{lag}}}$ in units of generations; one generation is defined as the average time between a network’s encounter with the same pattern class i.e., $N$. We characterize this lag $g_{{{}_{\text{lag}}}}$ by identifying the past pattern (at time $t-g_{{}_{\text{lag}}}N$) that has the maximum overlap with the network’s energy landscape at given time $t$:

$$g_{{{}_{\text{lag}}}}=\operatorname*{argmax}_{g\geq 0}\mathbb{E}\left[\braket{\sigma(t-g\,N)}{J(t)}{\sigma(t-g\,N)}\right]\equiv\operatorname*{argmin}_{g\geq 0}\mathbb{E}\left[{E_{{}_{\text{lag}}}(g)}\right]$$

(S18)

where we introduced the expected lagged energy $\mathbb{E}\left[{E_{{}_{\text{lag}}}(g)}\right]$. Here, the vector $\ket{\sigma(t)}$ refers to the pattern $\sigma$ presented to the network at time $t$, which can be from any of the pattern classes. Because of the translational symmetry in time in the stationary state, the lagged energy can also be expressed in terms of the overlap between a pattern at time $t$ and the network at a later time $t+gN$. We evaluate the lagged energy by substituting the expression for the network’s state $J(t+gN)$ from eq. S12, which entails,

	$\displaystyle\frac{2}{L-1}\mathbb{E}\left[{E_{\text{lag}}(g)}\right]$	$\displaystyle=-\frac{1}{L-1}\mathbb{E}\left[{\braket{\sigma(t)}{J(t+g\,N)}{\sigma(t)}}\right]$		(S19)
		$\displaystyle=-\mathbb{E}\left[{\frac{1}{L-1}(1-\lambda)^{Ng}\braket{\sigma(t)}{J(t)}{\sigma(t)}+\lambda\sum_{j=0}^{gN-1}(1-\lambda)^{gN-1-j}\braket{\sigma(t)}{\sigma(t+j)}^{2}}\right]$		(S20)
		$\displaystyle=\frac{2}{L-1}(1-\lambda)^{Ng}\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]-\lambda\sum_{i=0}^{g-1}\sum_{\alpha=0}^{N-1}(1-\lambda)^{gN-1-Ni-\alpha}\braket{\sigma^{N}(t)}{\sigma^{N-\alpha}(t+Ni+\alpha)}^{2}$		(S21)
		$\displaystyle=\frac{2}{L-1}(1-\lambda)^{Ng}\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]-\lambda\sum_{i=0}^{g-1}(1-\lambda)^{gN-1-Ni}\rho^{2Ni}$		(S22)
		$\displaystyle=-\lambda\left(\frac{(1-\lambda)^{N(g+1)}\rho^{2N}}{1-(1-\lambda)^{N}\rho^{2N}}+\frac{(1-\lambda)^{N(g+1)-1}-(1-\lambda)^{N-1}\rho^{2Ng}}{(1-\lambda)^{N}-\rho^{2N}}\right).$		(S23)

Here, we used the expression of the network’s matrix $J$ in eq. S8 to arrive at eq. S20, and then followed the procedure introduced in eq. S13 to arrived at the double-summation in eq. S21. We then used the equation for pattern overlap $\braket{\sigma^{\alpha}(t_{1})}{\sigma^{\nu}(t_{2})}=\delta_{\alpha,\nu}\rho^{|t_{2}-t_{1}|}$ to reduce the sum in eq. S22 and arrived at the result in eq. S23 by evaluating the geometric sum and substituting the empirical average for the energy $\mathbb{E}\left[{E_{{}_{\lambda,\rho}}(J,\sigma)}\right]$ from eq. LABEL:eq:energypattern_final.

We probe this lagged memory by looking at the performance ${\mathcal{Q}}$ for patterns that are correctly associated with their memory attractors (i.e., those with $\braket{\sigma_{\text{att}}}{\sigma}>0.8$). As shown in Fig. S1, for a broad parameter regime, the mean performance for these correctly associated patterns agrees well with the theoretical expectation ${\mathcal{Q}}_{{}_{\text{lag}}}=\rho^{g_{{}_{\text{lag}}}N}$, which is lower than the naive expectation ${\mathcal{Q}}_{0}$.

As shown in Fig. 1, large learning rates in networks with memory for evolving patterns result in the emergence of narrow connecting paths between the minima of the energy landscape. We refer to these narrow connecting paths as mountain passes. In pattern retrieval, the Monte Carlo search can drive a pattern out of one energy minimum into another minimum and potentially lead to pattern misclassification.

We use two features of the energy landscape to probe the emergence of the mountain passes.

First, we show that if a pattern is misclassified, it has fallen into a memory attractor associated with another pattern class and not a spuriously made energy minima. To do so, we compare the overlap of the attractor with the original pattern $|\braket{\sigma^{\alpha}_{\text{att}}}{\sigma^{\alpha}}|$ (i.e., the reconstruction performance of the patterns) with the maximal overlap of the attractor with all other patterns $\max_{\nu\neq\alpha}|\braket{\sigma^{\alpha}_{\text{att}}}{\sigma^{\nu}}|$. Indeed, as shown in Fig. S3A for evolving patterns, the memory attractors associated with $99.4\%$ of the originally stored patterns have either a large overlap with the correct pattern or with one of the other previously presented pattern. 71.3% of the patterns are correctly classified (stable patterns in sector I in Fig. S3A), whereas 28.1% of them are associated with a secondary energy minima after equilibration (unstable patterns in sector II in Fig. S3A). A very small fraction of patterns ($<1\%$) fall into local minima given by the linear combinations of the presented patterns (sector IV in Fig. S3A). These minima are well-known in the classical Hopfield model [44, 45]. Moreover, we see that equilibration of a random pattern (i.e., a pattern orthogonal to all the presented classes) in the energy landscape leads to memory attractors for one of the originally presented pattern classes. The majority of these random patterns lie in sector II of Fig. S3A), i.e., they have a small overlap with the network since they are orthogonal to the originally presented pattern classes, and they fall into one of the existing memory attractors after equilibration.

Second, we characterize the possible paths for a pattern to move from one valley to another during equilibration, using Monte Carlo algorithm with the Metropolis acceptance probability,

$$\rho(\sigma\to\sigma^{\prime})=\text{min}\left(1,e^{-\beta(E(J,\sigma^{\prime})-E(J,\sigma))}\right)$$

(S24)

We estimate the number of beneficial spin-flips (i.e., open paths) that decrease the energy of a pattern at the start of equilibration (Fig. S3B). The average number of open paths is smaller for stable patterns compared to the unstable patterns that are be miscalssified during retrieval (Fig. S3B). However, the distributions for the number of open paths largely overlap for stable and unstable patterns. Therefore, the local energy landscape of stable and unstable patterns are quite similar and it is difficult to discriminate between them solely based on the local gradients in the landscape. Fig. S4A shows that the average number of beneficial spin-flips grows with the mutation rate of the patterns but this number is comparable for stable and unstable patterns. Moreover, the unstable stored patterns (blue) have far fewer open paths available to them during equilibration compared to random patterns (red) that are presented to the network for the first time (Figs. S3B & S4A). Notably, on average half of the spin-flips reduce the energy of for random patterns, irrespective of the mutation rate. This indicates that even though previously presented pattern classes are statistically distinct from random patterns, they can still become unstable, especially in networks which are presented with evolving patterns.

It should be noted that the evolution of the patterns only indirectly contribute to the misclassification of memory, as it necessitates a larger learning rate for the networks to optimally operate, which in turn results in the emergence of mountain passes. To clearly demonstrate this effect, Figs. S3C,D, and S4D shows the misclassification behavior for a network trained to store memory for static pattern, while using a larger learning rate that is optimized for evolving patterns. Indeed, we see that pattern misclassification in this case is consistent with the existence of mountain passes in the network’s energy landscape.

We use spectral decomposition of the energy landscape to characterize the relative positioning and the stability of patterns in the landscape. As shown in Figs. S3, S4, destabilization of patterns due to equilibration over mountain passes occurs in networks with high learning rates, even for static patterns. Therefore, we focus on how the learning rate impacts the spectral decomposition of the energy landscape in networks presented with static patterns. This simplification will enable us to analytically probe the structure of the energy landscape, which we will compare with numerical results for evolving patterns.

We can represent the network $J$ (of size $L\times L$) that store a memory of $N$ static patterns with $N$ non-trivial eigenvectors $\ket{\Phi^{i}}$ with corresponding eigenvalues $\Gamma_{i}$, and $N-L$ degenerate eigenvectors, $\ket{\Psi^{k}}$ with corresponding trivial eigenvalues $\gamma_{k}=\gamma=-1$:

$\displaystyle J=\sum_{i=1}^{N}\Gamma_{i}\ket{\Phi^{i}}\bra{\Phi^{i}}+\sum_{k=1}^{L-N}\gamma_{k}\ket{\Psi^{k}}\bra{\Psi^{k}}.$

(S25)

The non-trivial eigenvectors span the space of the presented patterns, for which the recognition energy can be expressed by

$\displaystyle E(J,\sigma^{\alpha})=-\frac{1}{2}\sum_{i=1}^{N}\Gamma_{i}\braket{{\sigma}^{\alpha}}{\Phi^{i}}\braket{\Phi^{i}}{{\sigma}^{\alpha}}.$

(S26)

An arbitrary configuration $\chi$ in general can have components orthogonal to the $N$ eigenvectors $\ket{\Phi^{i}}$, as it points to a vertex of the hypercube, and should be expressed in terms of all the eigenvectors $\{\Phi^{1},\dots,\Phi^{N},\Psi^{1},\dots,\Psi^{L-N}\}$:

$\displaystyle E(J,\chi)=-\frac{1}{2}\left(\underbrace{\sum_{i=1}^{N}\Gamma_{i}\braket{{\chi}}{\Phi^{i}}\braket{\Phi^{i}}{{\chi}}}_{\text{stored patterns}}+\underbrace{\sum_{k=1}^{L-N}\gamma\braket{{\chi}}{\Psi^{k}}\braket{\Psi^{k}}{{\chi}}}_{\text{trivial space}}\right).$

(S27)

Any spin-flip in a pattern (e.g., during equilibration) can be understood as a rotation in the eigenspace of the network (eq. S27). As a first step in characterizing these rotations we remind ourselves of the identity

$\displaystyle\ket{{\chi}}=\sum_{i=1}^{N}\braket{\Phi^{i}}{{\chi}}\ket{\Phi^{i}}+\sum_{k=1}^{L-N}\braket{\Psi^{k}}{{\chi}}\ket{\Psi^{k}},$

(S28)

with the normalization condition

$\displaystyle\sum_{i=1}^{N}\left(\braket{\Phi^{i}}{{\chi}}\right)^{2}+\sum_{k=1}^{L-N}\left(\braket{\Psi^{k}}{{\chi}}\right)^{2}=1.$

(S29)

In addition, since the diagonal elements of the network are set to $J_{ii}=0$ (eq. S7), the eigenvalues should sum to zero, or alternatively,

$$\sum_{i=1}^{N}\Gamma_{i}=-\sum_{k=1}^{L-N}\gamma_{k}=L-N.$$

(S30)

To asses the stability of a pattern ${\sigma}^{\nu}$, we compare its recognition energy $E(J,{\sigma}^{\nu})$ with the energy of the rotated pattern after a spin-flip $E(J,\tilde{\sigma}^{\nu})$. To do so, we first consider a simple scenario, where we assume that the pattern $\sigma^{\nu}$ has a large overlap with one dominant non-trivial eigenvector $\Phi^{A}$ (i.e., $\braket{\sigma^{\nu}}{\Phi^{A}}^{2}=m^{2}\approx 1$). The other components of the pattern can be expressed in terms of the remaining $N-1$ non-trivial eigenvectors with mean squared overlap $\frac{1-m^{2}}{N-1}$. The expansion of the recognition energy for the presented pattern is restricted to the $N$ non-trivial directions (eq. S27), resulting in

$\displaystyle E(J,{\sigma}^{\nu})=-\frac{1}{2}\left(m^{2}\Gamma_{A}+\sum_{i\neq A}\frac{1-m^{2}}{N-1}\Gamma_{i}\right)=-\frac{1}{2}\left(m^{2}\Gamma_{A}+(1-m^{2})\tilde{\Gamma}\right),$

(S31)

where $\tilde{\Gamma}=\frac{1}{N-1}\sum_{i\neq A}\Gamma_{i}=\frac{N\bar{\Gamma}-\Gamma_{A}}{N-1}$ is the mean eigenvalue for the non-dominant directions.

A spin-flip ($\ket{\sigma^{\nu}}\rightarrow\ket{\tilde{\sigma}^{\nu}}$ ) can rotate the pattern out of the dominant direction $\Phi^{A}$ and reduce the squared overlap by $\epsilon^{2}$. The rotated pattern $\ket{\tilde{\sigma}^{\nu}}$ in general lies in the $L$-dimensional space and is not restricted to the $N$-dimensional (non-trivial) subspace. We first take a mean-field approach in describing the rotation of the pattern after a spin-flip. Because of the normalization condition (eq. S29), the loss in the overlap with the dominant direction should result in an average increase in the overlap with the other $L-1$ eigenvectors by $\frac{\epsilon^{2}}{L-1}$. The energy of the rotated pattern after a spin-flip $E(J,\tilde{\sigma}^{\nu})$ can be expressed in terms of all the $L$ eigenvectors (eq. S27),

	$\displaystyle E(J,\tilde{\sigma}^{\nu})$	$\displaystyle=-\frac{1}{2}\left[(m^{2}-\epsilon^{2})\Gamma_{A}+\sum_{i\neq A}\left(\frac{1-m^{2}}{N-1}+\frac{\epsilon^{2}}{L-1}\right)\Gamma_{i}+\sum_{k}\frac{\epsilon^{2}}{L-1}\gamma_{k}\right]$
		$\displaystyle=E(J,{\sigma}^{\nu})+\frac{\epsilon^{2}}{2}\left[\Gamma_{A}-\frac{1}{L-1}\left(\sum_{i\neq A}\Gamma_{i}+\sum_{k}\gamma_{k}\right)\right].$		(S32)
		$\displaystyle=E(J,{\sigma}^{\nu})+\frac{\epsilon^{2}}{2}\Gamma_{A}\left(1+\frac{1}{L-1}\right).$		(S33)

where in eq. S33 we used the fact that the eigenvalues should sum up to zero. On average, a spin-flip $\ket{{\sigma}^{\nu}}\rightarrow\ket{\tilde{\sigma}^{\nu}}$ increases the recognition energy by $E(J,\tilde{\sigma}^{\nu})-E(J,{\sigma}^{\nu})=\frac{\epsilon^{2}}{2}\Gamma_{A}\left[1+\mathcal{O}\left(L^{-1}\right)\right]$. This is consistent with the results shown in Figs. S3B,D and Figs. S4A,D, which indicate that the majority of the spin-flips keep a pattern in the original energy minimum and only a few of the spin-flips drive a pattern out of the original attractor.

In the analysis above, we assumed that the reduction in overlap with the dominant eigenvector $\epsilon^{2}$ is absorbed equally by all the other eigenvectors (i.e., the mean-field approach). In this case, the change in energy is equally distributed across the positive and the negative eigenvalues ($\Gamma$’s and $\gamma$’s in eq. S32), resulting in an overall increase in the energy due to the reduced overlap with the dominant direction $\ket{\Phi^{A}}$. The destabilizing spin-flips are associated with atypical changes that rotate a pattern onto a secondary non-trivial direction $\ket{\Phi^{B}}$ (with positive eigenvalue $\Gamma_{B}$), as a result of which the total energy could be reduced. To better characterize the conditions under which patterns become unstable, we will introduce a perturbation to the mean-field approach used in eq. S33. We will assume that a spin-flip results in a rotation with a dominant component along a secondary non-trivial direction $\ket{\Phi^{B}}$. Specifically, we will assume the reduced overlap $\epsilon^{2}$ between the original pattern $\ket{{\sigma}^{\nu}}$ and the dominant direction $\ket{\Phi^{A}}$ is distributed in an imbalanced fashion between the other eigenvectors, with a fraction $p$ projected onto a new (non-trivial) direction $\ket{\Phi^{B}}$, while all the other $L-2$ directions span the remaining $(1-p)\epsilon^{2}$. In this case, the energy of the rotated pattern is given by

	$\displaystyle E(J,\tilde{\sigma}^{\nu})$	$\displaystyle=-\frac{1}{2}\left[(m^{2}-\epsilon^{2})\Gamma_{A}+\left(\frac{1-m^{2}}{N-1}+p\epsilon^{2}\right)\Gamma_{B}+\sum_{i\neq A,B}\left(\frac{1-m^{2}}{N-1}+\frac{(1-p)\epsilon^{2}}{L-2}\right)\Gamma_{i}+\sum_{k}\frac{(1-p)\epsilon^{2}}{L-2}\gamma_{k}\right]$
		$\displaystyle=E(J,{\sigma}^{\nu})+\frac{\epsilon^{2}}{2}\left[\Gamma_{A}-p\,\Gamma_{B}+\mathcal{O}\left(L^{-1}\right)\right].$		(S34)

Therefore, a spin-flip is beneficial if $\Gamma_{A}<p\,\Gamma_{B}$. To further concretize this condition, we will estimate the typical loss $\epsilon^{2}$ and gain $p\epsilon^{2}$ in the squared overlap between the pattern and its dominating directions due to rotation by a single spin-flip.

Let us consider a rotation $\ket{\sigma^{\nu}}\to\ket{\tilde{\sigma}^{\nu}}$ by a flip in the $n^{th}$ spin of the original pattern $\ket{\sigma^{\nu}}$. This spin flip reduces the original overlap of the pattern $m=\braket{\sigma^{\nu}}{\Phi^{A}}$ with the dominant direction $\ket{\Phi^{A}}$ by the amount $\frac{2}{\sqrt{L}}\Phi_{n}^{A}$, where $\Phi_{n}^{A}$ is the $n^{th}$ entry of the eigenvector $\ket{\Phi^{A}}$. Since the original overlap is large (i.e., $m\simeq 1$), all entries of the dominant eigenvector are approximately $\Phi_{i}^{A}\simeq{1}/{\sqrt{L}},\,\forall i$, resulting in a reduced overlap of the rotated pattern $\braket{\sigma^{\nu}}{\Phi^{A}}=m-\frac{2}{L}$. Therefore, the loss in the squared overlap $\epsilon^{2}$ by a spin flip is given by

$\displaystyle\epsilon^{2}=\braket{{\sigma}^{\nu}}{\Phi^{j}}^{2}-\braket{\tilde{\sigma}^{\nu}}{\Phi^{j}}^{2}=m^{2}-\left(m^{2}-4\frac{m}{L}+4\frac{1}{L^{2}}\right)=4\frac{m}{L}+\mathcal{O}(\frac{1}{L^{2}}).$

(S35)