In an artificial neural network working as associative memory, one deals with a network of N neurons of which each one has state s_i (i = 1…N) that can be “active” (s_i = 1) or “quiescent” (s_i = −1). The configuration of the whole network is given by the vector s¯≡{s1,s2,..sN} and its temporal evolution follows the parallel non-linear dynamics: (1)si(t+Δt)=E[si(t)]≐sign[∑j=1NJijsj(t)],  where J = {J_ij} is the connection matrix. We set external inputs to be equal to 0. We assume a symmetric bimodal distribution for the synaptic polarities in the wiring matrix J, so 50% of the connections are excitatory and 50% inhibitory. After a transient time related to the finite value of N, the network reaches a fixed point, s_i = E[s_i], or a limit cycle of length L, s_i = E^(L)[s_i].

Previous work has studied the cycle length and transient time distribution as a function of the properties of J (Gutfreundt et al., 1988; Sompolinsky et al., 1988; Derrida, 1989; Bastolla et al., 1997). In order to work as an associative memory, the matrix J must be tailored in such a way that one or more patterns of neurons are fixed points of the dynamics in Equation (1), i.e., they are the “memory elements” stored in the network. To store one pattern ξ̄, the connection matrix is simply the dyadic form given by J_ij = ξ_iξ_j ¹ , while to store a generic number P of patterns ξ̄μ (μ = 1…P) one follows the storage prescription of Cooper (1973) and Cooper et al. (1979), who exploited an old idea which goes back to Hebb (1949) and Eccles (1953) and which states that the change in synaptic transmission is proportional to the product of the signals of pre and post-synaptic neurons. The process for which each matrix element is appropriately determined is called learning. Specifically, the “Hebbian” rule results in the following expression for the connectivity matrix, (2)Jij=1P∑μ=1Pξiμξjμ.  The set of vectors ξ̄mu is known as “training set.” In this case, it is not guaranteed that each ξ̄μ is a fixed point. In other words, ξ̄μ is stable in probabilistic sense. Further, the probability for ξ̄μ to be a fixed point depends on the values of P and N. This dependence has been first studied by Hopfield (1982); Hopfield et al. (1983); Hopfield (1984) who concluded that the retrieval of the memory stored in the Hebbian matrices is guaranteed up to a P-value which is a critical fraction on the number of network nodes N of the order of 10–20%. Above this value, the associative memory quickly degrades. Following these studies, Amit et al. (1985a,b), who noticed the similarity between the Hopfield model for the associative memory and the spin glasses, developed a statistical theory for the determination of the critical P/N ratio, that turned out to be ≈ 0.14, in good agreement with the previous Hopfield estimation. Above P=0.14N the number of errors is so large that the network based on the Hebbian matrix is no longer capable to work as an associative memory. All these studies assumed a modified form of Equation (2): the diagonal elements of J are forced to be zero.

In order to demonstrate the validity of our analytic results (see Section 3), we perform numerical simulations by evolving the network model as described in Equation (1). We design the default network by fixing the NxN recurrent connections as given in Equation (2), by randomly assigning the value ±1 to ξiμ and retaining the diagonal elements. So, the N(N−1) connections are 50% excitatory and 50% inhibitory and the N neurons can form self-connections. We then run simulations by varying the size of the network, N = 50, .., 200 and the number of stored memories in Equation (2), P = 1, …, 2000. Finally, for each pair of N and P, we perform 1000 different random realizations.

All P patterns introduced in Equation (2) are given as input to the network and their dynamics is followed until the network reaches the equilibrium state. The initial patterns are chosen among those that were stored in the adjacency matrix and that have been randomly chosen in the designing of the network. Evolved patterns were recorded at each time step and compared with the initial one. Then, if the evolved pattern is different from the initial state, we calculated the temporal evolution of the distortion (number of wrong bits) and determined the probability that one of the bits was wrong, the probability that the whole vector was exactly recovered, and the number of memory patterns that could not be recovered, as a function of N and P. Basically, to calculate the storage capacity, it is sufficient to determine all these quantities by using the distortion between the stimulus (the stored memory) and its first evolved pattern.

In order to investigate the maximum storage memory of our model, we calculate the one-step dynamical evolution. We give as input a vector of the training set and we calculate a single step of the dynamical evolution according to Equation (1). Then, we compare the output with the input. We aim to look whether or not a vector, ξ̄μ, belonging to the training set, is truly a fixed point. If ξ̄μ is a fixed point, the output coincides with the input, and the recovery has been successful. If ξ̄μ is not a fixed point, the two vectors differ for at least a single bit. We now derive an analytical expression for the probability that the recovery of a stored pattern was not successful. The first step is to find the probability p_B that -given the matrix J of Equation (2)- a single element of the vector (a “bit”) was wrong, i.e., the probability that E[ξiμ]≠ξiμ. Basically, we need to evolve a vector ξ̄mu (from the training set) for one step and count how many bits of its time evolution are different from the bits of ξ̄mu itself. Obviously, if ξ̄mu actually is a fixed point, this distance vanishes. On the contrary, ξ̄mu is not a fixed point, the network has made a recovery error. Thus, p_B (or better, p_V, as we see in the next paragraph) measures “how many” training set vectors are not fixed points. The argument of the sign function in Equation (1) is Aiμ=∑j=1N∑ν=1Pξiνξjνξjμ, this contains NP terms among which there are N+P−1 terms (those with j=i ² and those with ν=μ) where two out of the three ξ of the product are equals to each other ξiν and the third is ξiμ. Thus Aiμ=(N+P−1)ξiμ+Tiμ, with Tiμ=∑j≠iN∑ν≠μPξiνξjνξjμ. The first term is the “coherent” one, its sign is identical to ξiμ, and it will -if dominant- guarantee that ξ̄μ is a fixed point of the dynamics. The second term Tiμ , on the contrary, is “noise” and its presence can either reinforce or weaken the stability of ξiμ as fixed point. Specifically, if |Tiμ| >(N+P−1) and sign(Tiμ)≠ξiμ, then the i-th bit of the vector ξ̄μ will turn out to be wrong. The quantity Tiμ is the sum of (N−1)(P−1) statistically independent terms, each one being +1 or −1. Therefore, for large enough P and N, its distribution N(T) can be approximated by a gaussian with zero mean and standard deviation (N-1)(P-1): (3)N(T)=e−T2/(2(N−1)(P−1))2π(N​−​1)(P​−​1).  It is now straightforward to determine the probability that |Tiμ| >(N+P−1) and sign(Tiμ)≠ξiμ, thus that one of the bits of E[ξiμ] was wrong, as pB=∫N+P-1∞dTN(T). In conclusion: (4)pB=12[1−erf(N​+​P​−​12(N​−​1)(P​−​1))].  It is worth to note that this expression is symmetric under the exchange of P with N, and that for large P and N, with P/N = 1, it tends to (1−erf(2))/2≈0.02275 which corresponds to the maximum of probability in a wrong recovery of a single bit (see Figures 1, 2).

The second step is the determination of the probability p_V that one of the P vectors encoded into the connection matrix (the training set) turns out not be a fixed point. If only a single bit of the vector is wrong, the whole vector is considered “wrong.” Since there are N bits that can be wrong, the probability p_V will be much higher than p_B. The calculation is straightforward, in order not to be wrong, all the bits of the vector ξ̄μ must be right, thus p_V = 1−(1-pB)N, therefore: (5)pV=1−[12+12erf(N​+​P​−​12(N​−​1)(P​−​1))]N.  Finally, the number, N_V, of memory vectors that are not recovered, i.e., that are not true fixed points of the dynamics is given by Pp_V, that is: (6)NV=[1−[12+12erf(N​+​P​−​12(N​−​1)(P​−​1))]N]P

Equations (4), (5) and, in particular, Equation (6) represent the main result of this work. Before showing their validity, via a comparison with numerical simulations, and discussing their relevance in the framework of artificial neural networks, it is important to present the asymptotical approximation for N_V. The argument of the error function, for either P≫N or P≪N, is large, and can be expanded as erf(x) ≈ 1−exp(-x2)/(xπ). Furthermore, as p_B is exponentially small with N (or P) for large N (P), we use, (1-pB)N≈(1-NpB). Thus, for large N or large P: (7)pV≈N3/2P1/2e−(N+P)22NP2π(N+P) (8)NV≈N3/2P3/2e−(N+P)22NP2π(N+P) We note that, while in the exact expression for N_V (Equation 6) the P↔N exchange symmetry is lost, in the approximate form the symmetry is recovered.

For sake of comparison with the previous literature, it is also useful to express the main results as a function of α≐P/N. Equations (4) (for large N) and (8) read: (9)pB≈12[1−erf(1+α2α)]  (10)NV≈NP12πα1+αe−(1+α)22α. While p_B only depends on α, N_V clearly is an extensive observable, being proportional to P and N. Furthermore, both expressions keep their symmetry with respect to the exchange of P and N, thus to the exchange of α with 1/α. The last observation anticipates that there must exists a region at large α-values where the same features are observed as at small values of α.

To check the predictions of our network model, we have simulated the Model (1) and studied the dynamics for several values of N and P, in the range of few hundred, see Section 2.3 for details. In the numerical analysis, the P memory vectors have been randomly chosen and used to construct the connection matrix J. Next, we tested whether or not the stored memories were fixed points of the dynamics. The values of p_B, p_V and N_V were calculated by averaging over (up to) 1000 different random realizations of ξ̄μ. The results of the numerical simulations are reported (dots) in Figure 1, together with the analytical Expressions (4)–(6) (lines). The three panels refer to the three quantities p_B (Figure 1A), p_V (Figure 1B), and N_V (Figure 1C) as a function of P for the selected values of N, as reported in the legend. From Figure 1, we observe that on increasing P, at fixed N, both the single bit probability error, the probability of recovery error (P_V), and the number of wrong recoveries N_V, after a first fast increase, reach a maximum (equal to 0.02275 for p_B, close to one for p_V, and larger than N for N_V) then start to decrease, tending to zero for very large P-values.

To better emphasize this behavior, the same quantities are reported (analytic results only) as a function of α in Figure 2 (linear scale) and in Figure 3 (log scale) for selected N. The dotted lines in panels C of both figures represent N_V = P, i.e., indicate the case of “totally wrong recovery.” Due to the already observed α↔1/α symmetry, the asymptotic curve in Figure 3A appears with a left-right symmetry around α = 1. From Figures 2, 3, we can clearly identify two regions of high recovery efficiency. The low α region, already studied many years ago by Hopfield (1982); Hopfield et al. (1983); Hopfield (1984) and Amit et al. (1985a,b), shows the existence of a quick transition toward “loss of memory recover” on increasing α around α ≈ 0.14. The second region at large α-values is not yet explored.

Although the value α = 1 (P = N) represents traditionally a sort of limit in the computation of the storable memories in a RNN, there is no reason why not to store more than N memory elements in a network of N neurons, that by construction allows 2^N possible patterns. Indeed, the number of fixed points in a (random) symmetric matrix is known to be, for fully connected symmetric matrices as in our case, exponentially large with N (Tanaka and Edwards, 1980). Specifically, the number of fixed points P_o is equal to P_o = exp(γN), with γ ≈ 0.2. P_o, much larger than N, can be considered a natural limit for P.

The recovery efficiency increases for large P. In fact, the coherent term in the argument of the sign function increases linearly with P and the noise increases as P^1/2. For large P, the relative weight of the noise decreases as P^−1/2, this allows to store a large number of memories in a relatively small neural network.

For practical purposes, as for example in the design of an artificial neural network with high efficiency (large storage capacity) and effectiveness (low recovery error rates), it is important to study (Equation 6, and its approximation in Equation 8) and, in particular, to find the conditions for which the network shows “perfect recovery.” Let's define perfect recovery as the state where the number of retrieval errors N_V is smaller than one.

In Figure 4 we show the contour plot of the (decimal) logarithm of N_V, from Equations (6) and (8), in the P-N range [0–100]. The full lines are the loci of the points where log₁₀(N_V) equals 0, 0.4, 0.8, 1.2, and 1.6, as indicated on the right side of the figure. The dashed lines are the same level lines for the (logarithm of the) approximate form of N_V reported in Equation (8). As can be observed, for N_V ≈ 1, the approximation (Equation 8) for N_V is highly accurate, indicating that this approximation can be safely applied to find the “perfect recovery” condition.

In the P-N plane the existence of two regions (small and large α) where the perfect recovery (N_V = 1, red lines) takes place can be easily observed and the result is symmetric under the exchange of P and N. In the already explored small α region, we also show (full blue line) the P = 0.14N condition. Similar to the high α region, it is important to find a simple relation between N and P identifying the N_V = 1 condition. We aim, therefore, to obtain a function P(N) which returns, at given N, the P-value such that N_V = 1. We write the prefactor NP in Equation (11) as αN² and exploit the α≫1 limit, so to obtain NV≈N2α1/2 exp(-α/2)/2π. The equation N²α^1/2exp(-α/2)/2π=1 can be squared, αexp(−α) = 2π/N⁴, and solved with respect to α, to give α=-W-1(-2π/N4), where W₋₁(x) is the second real branch of the Lambert function (Olver et al., 2010). In conclusion, the “perfect recovery condition” is satisfied -for each N-value- if we chose to store a number of memories larger than P(N) given by: (11)P(N)=−NW−1(−2π/N4).  For practical purposes, for large enough N, we can use the small-argument expansion of the Lambert function −W₋₁(−x) ≈ −ln(x)+ln(−ln(x)) (Corless et al., 1996), to have: (12)P(N)=N[ln(N42π)+ln(ln(N42π))]. The results for P(N) are shown in Figure 5 as a function of N in the range 1–1000. The black line represents the exact, numerical, solution to N_V = 1, with N_V in Equation (6), the blue line is the expression for P(N) in Equation (11), while the red line is those in Equation (12).

It is important to note that the presence of a decrease of the retrieval error probability at high P, or α, values is due to the presence of non-zero diagonal elements in the J matrix that creates a coherent term of weight P. Indeed, repeating the rationale leading to Equation (4) with the assumption that J_ii = 0, would give rise to the same (Equations 4–6) but with the numerator of the argument of the error functions equal to N − 1 instead of to N + P − 1. This is shown graphically in Figure 6 where we compare for N = 50, both theoretically (full line) and numerically (full dots), the quantities p_B, p_V, and N_V as a function of P in the two cases: diagonal elements in Equation (2) (black) and diagonal elements forced to vanish (orange).

The stabilization of the fixed points ξ̄μ in the high storage region arises from the presence of the non-zero diagonal elements. Asymptotically, on increasing P, the diagonal elements growth coherently and the J matrix tends to become the unit matrix. However, the dynamics (see Equation 1) dictated by the matrix J does not tend to the dynamics dictated by the unit matrix. In the latter case, indeed, all the 2^N state vectors should become fixed points and the network should loose on important feature: the capability to distinguish between the stored memories (the vectors ξ̄μ, for μ = 1…P) and the spurious fixed points, all the vectors ζ̄ not belonging to the set ξ̄μ but such that E[ζ̄] = ζ̄. To study this property, we have calculated the probability that a (randomly chosen) vector ζ̄ (different from all the ξ̄μ used to build the J matrix) was recognized as a “memory” from the network dynamics. To be consistent with the previous notation (where we called p_B and p_V the probability of errors, not that of correct retrieval of the memory states) we define p̄B (p̄V) as the probability of correctly not retrieving a vector not belonging to the training set. More specifically, the quantity p̄V is the probability that one dynamical step after presenting a vector ζ not belonging to the training set to the network, the output a vector is different from ζ.” More specifically, the quantity p̄V is the probability that presenting a vector ζ not belonging to the training set to the network, after one dynamical step we found as output a vector different from ζ. Similarly for p̄B. It turns out that ³: (13)p¯B=12[1−erf(P2(N​−​1)(P​−​1))] (14)p¯V=1−[12+12erf(P2(N​−​1)(P​−​1))]N. In Figure 7 we report the comparison of the P dependence of p_B and p̄B (Figure 7A) and that of p_V and p̄V (Figure 7B). As usual, full lines are the theoretical results, while the full dots are the outcome of the numerical simulation. Black data are for the “memory states,” while the green ones are for the “spurious state.” As can be seen, the spurious state becomes more and more “present” in the set of memories stored by the network as P increases. It seems however that also at high P-values the retrieval of the memory states is reasonably good and that of the spurious states reasonably bad.

To be quantitative on this point, we rewrite Equation (14) in its large N limit: (15)p¯V≈N3/2P−1/2e − P2N2π  and compare it with Equation (7). In particular, is interesting to calculate the ratio, ρ, between the probability of wrong retrieval of a spurious state and that of a memory state: ρ = p̄V/pV. From Equations (7) and (15) it turns out: (16)ρ=(N+PP)  e(N+P)22NPe−P2N. This quantity only depends on α: (17)ρ=(1+αα)  e(1+12α) and ρ has a finite high α limit: (18)limα→∞ρ=e. In other words, although the number of spurious attractors tends to increase for P ≫ N, the vectors encoded into the system through the connection matrix are retrieved with an efficiency almost three times better than for the spurious states.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.