Paper deep dive
Multiplicity of Stable Attractors in Disordered Neural Models
Raffaele Marino, Roberto Livi, Antonio Politi
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/2/2026, 12:56:14 PM
Summary
This paper investigates the multiplicity of stable fixed points (SFPs) in a model of neural ordinary differential equations (nODEs) with disordered random coupling matrices. Using large-deviation statistics and a perturbative method, the authors estimate the number of surviving SFPs as coupling strength increases. They find that for moderate coupling, the number of SFPs decreases via saddle-node bifurcations, but an exponentially large number survive. The study compares asymmetric (Ginibre ensemble) and symmetric coupling matrices, finding similar qualitative behavior in SFP multiplicity despite different dynamical regimes (chaos vs. gradient flow).
Entities (7)
Relation Signals (6)
Neural Ordinary Differential Equations → hasfeature → Stable Fixed Points
confidence 95% · The model exhibits an exponentially large number of equilibrium points (stable fixed points).
Large-Deviation Statistics → usedfor → Estimating SFP Multiplicity
confidence 94% · We show how large-deviation statistics allows one to obtain reliable estimates of the multiplicity of stable fixed-points
Ginibre Ensemble → characterizes → Asymmetric Coupling
confidence 92% · The entries A_ij follow a standard Gaussian distribution... so that G = 1/sqrt(N) A is an element of the real Ginibre ensemble.
Coupling strength → causes → Saddle-Node Bifurcation
confidence 90% · The depletion of SFPs is due to a sequence of tangent (also called saddle-node or fold) bifurcations... upon increasing g... one SFP approaches one of the neighbouring saddles... and eventually the two points mutually annihilate
Symmetric Coupling → leadsto → Gradient Dynamics
confidence 88% · Symmetric matrices represent an important class of models, characterized by a strictly gradient dynamics
Asymmetric Coupling → allows → Chaos
confidence 87% · the coupling matrix is non-symmetric, thus allowing for a chaotic evolution in the limit of large coupling values.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We show how large-deviation statistics allows one to obtain reliable estimates of the multiplicity of stable fixed-points in a model of neural ordinary differential equations previously employed in computational tasks. The result is obtained by developing a suitable perturbative method in the amplitude of the disorder. It turns out that for not-too-large coupling strengths there are no qualitative differences between the symmetric case, when the dynamics is a purely gradient evolution, and the asymmetric case, when limit cycles and chaos can, in principle, arise. The selection of this specific model is dictated by pedagogical reasons, but we are confident that the approach can be extended to other many-degree-of-freedom dynamical models characterized by different classes of random coupling matrices.
Tags
Links
- Source: https://arxiv.org/abs/2607.22047v1
- Canonical: https://arxiv.org/abs/2607.22047v1
Trouble viewing inline? Open PDF directly →
Full Text
32,987 characters extracted from source content.
Expand or collapse full text
Multiplicity of Stable Attractors in Disordered Neural Models Raffaele Marino raffaele.marino@butterflydecisions.com Butterfly Decisions srl, Via dei Principati 74 - 84122 Salerno, Italy Roberto Livi roberto.livi@unifi.it University of Florence, Department of Physics and Astronomy, Via G. Sansone 1 - 50019 Sesto Fiorentino (FI), Italy Istituto dei Sistemi Complessi, CNR, Via Madonna del Piano 10 – 50019 Sesto Fiorentino (FI), Italy Antonio Politi a.politi@abdn.ac.uk Istituto dei Sistemi Complessi, CNR, Via Madonna del Piano 10 – 50019 Sesto Fiorentino (FI), Italy Institute of Pure and Applied Mathematics, Department of Physics, University of Aberdeen, Aberdeen AB24 3UE, United Kingdom (July 24, 2026) Abstract We show how large-deviation statistics allows one to obtain reliable estimates of the multiplicity of stable fixed-points in a model of neural ordinary differential equations previously employed in computational tasks. The result is obtained by developing a suitable perturbative method in the amplitude of the disorder. It turns out that for not-too-large coupling strengths there are no qualitative differences between the symmetric case, when the dynamics is a purely gradient evolution, and the asymmetric case, when limit cycles and chaos can, in principle, arise. The selection of this specific model is dictated by pedagogical reasons, but we are confident that the approach can be extended to other many-degree-of-freedom dynamical models characterized by different classes of random coupling matrices. Stable Attractors, Disordered Neural Networks, Order-Statistics, Large-Deviation Introduction - Models of random neural networks defined in terms of first-order ordinary differential equations have been recently recognized as effective candidates for solving various computational problems. For instance, in Wainrib and Touboul (2013) a continuous version of the Sherrington-Kirkpatrick model has been investigated as a testbed of the relations between topological and dynamical complexity. There, by tuning the variance σ of the i.i.d. Gaussian entries of the random coupling matrix among the continuous variables (neurons), one finds a phase transition from a single equilibrium state to a chaotic one. The main result of that study is that for any finite number N of neurons in the network the model exhibits an exponentially large number of equilibrium points. Following quite different motivations, another class of these models, known as Coherent Ising Machines (CIM’s), has been studied as nonconventional architectures for finding approximate solutions of large-scale combinatorial optimization problems (e.g. Ghimenti et al. (2026); Syed and Berloff (2023); Syed et al. (2026)). In CIM’s the continuous neural variables are subject to a local double-well potential and they are coupled via a symmetric matrix, whose random entries are i.i.d. Gaussian variables. Due to its gradient-dynamic structure, CIM’s are typically employed for identifying global energy minima by making use of an annealing process, starting from a single trivial ground state and evolving towards stable fixed points. Another model, inspired by Recurrent Neural Networks Kim et al. (2019); Rungratsameetaweemana et al. (2025) and much similar to CIM’s, is the set of neural Ordinary Differential Equations (nODEs), recently analyzed in Marino et al. (2024). It again deals with continuous neural variables in a double-well local potential. However, the coupling matrix is non-symmetric, thus allowing for a chaotic evolution in the limit of large coupling values. Such a model has been employed in a region of parameters space containing spontaneous or planted stable fixed-points. This dynamical system has been found able to perform standard learning tasks, making use of implicit feed-forward modules Marino et al. (2024). The procedure exploits the presence of large sets of stable attractors (typically, fixed-points), employed as targets of a learnable dynamics, where the Euler numerical integration outlines the recurrent architecture of deep-learning algorithms. Moreover, it has been shown that effective training strategies can be applied to enforce the access to planted attractors, representative of classification patterns Gagliani et al. (2026); Chicchi et al. (2025); Marino et al. (2025). All of these models were investigated having in mind specific computational tasks, but no systematic effort has been devoted to understanding the underlying dynamical structure, common to this wide class of models of nODEs. In this Letter we provide a preliminary account of such aspects, by focusing our attention on the dynamical and statistical complexity of the model studied in Marino et al. (2024). More precisely, we describe how the exponentially large number, 2N2^N (N being the system size), of stable fixed points present for zero coupling, g=0g=0, decreases, when g is switched on. An exponentially large number of stable fixed-points survive, at least up to some finite value of g∼(1)g O(1). In this range of g values we construct a large-deviation functional, by exploiting the extreme-value statistics, associated with the fixed-points depletion mechanism. All of this is achieved thanks to a perturbative criterion. To our knowledge, this is a fully novel strategy, which discloses yet unexplored perspectives in relating many-degree-of-freedom dynamical systems with the statistics of disordered models, such as spin-glasses. Final remarks about the dynamics in the large coupling limit envisage future studies about the dynamical phases of the model, characterized by the presence of stable limit cycles and chaotic attractors. The model - We study a set of N particles/neurons xix_i, which satisfy the following dynamical equations dxidt=−xi(xi2−1)+gN∑j=1NAijxj. dx_idt=-x_i(x_i^2-1)+ g N\, _j=1^N\,A_ijx_j\,. (1) Each particle is subject to a local double-well potential, whose minima are located at ±1± 1 and separated by a maximum in 0, namely V(xi)=14(xi2−1)2V(x_i)= 14(x_i^2-1)^2. Additionally, the particles mutually interact via an N×N× N random matrix A, while g denotes the coupling strength. The entries AijA_ij follow a standard Gaussian distribution (0,1)N(0,1), so that =1NG= 1 N\,A is an element of the real Ginibre ensemble. Given the symmetry of the nODEs, dynamics (1) is invariant under the change of sign, i.e. xi→−xi\x_i\→\-x_i\. For g=0g=0 the dynamics reduces to a Cartesian product of overdamped processes. Given a generic initial condition, each xix_i moves to its closest minimum and the overall dynamical system collapses onto one of the 2N2^N stable fixed points (SFPs), x→∗(k)(0)\ x^*(k)(0)\ with 1≤k≤2N1≤ k≤ 2^N, corresponding to all sequences of ±1± 1. As soon as the coupling g is switched on, mutual interactions arise and it is natural to expect deviations of the coordinates of the SFPs from the initial ±1± 1 values. Results - To build intuition on the qualitative organization of the dynamics, we first analyze the system for small sizes N. In this setting, the attractor landscape can be examined systematically as a function of the coupling parameter g. Since we expect the number of SFPs to depend on the realization of the matrix, the data reported in Fig. 1 are obtained by averaging over different matrices. Figure 1: Fraction of PN(g)P_N(g) of SFPs as a function of g for several sizes. Solid curves are simulation estimates averaged over different matrices (1010 for N≤20N≤ 20, 22 for N=25N=25); shaded bands indicate the uncertainty on the mean (error bar on the average). Dashed curves are the perturbative predictions obtained by integrating the minima density. There we report the average fraction PN(g)P_N(g) of SFPs as a function of g, for several values of N. For small g, the curves remain close to 11 suggesting that the number of SFPs is constant. Upon increasing g, a significant drop (notice the logarithmic vertical scale) is observed, and the decrease becomes progressively steeper as N grows. Beyond this drop, there are much fewer SFPs, and for g>1g>1 limit cycles as well as chaotic attractors may appear (data not reported). So long as SFPs persist, we have observed that they retain the same sign code: namely, the vector s→(k)(g)=sign(x→∗(k)(g))∈±1N s^(k)(g)\;=\;sign\! ( x^*(k)(g) )∈\± 1\^N (2) coincides with the sign pattern present at g=0g=0. In this sense, the surviving fixed points can be interpreted as descendants of the SFPs at g=0g=0, even though the values of their components change with g. Figure 2: For small g, the stable fixed point (blue, lower branch) is a descendant of the g=0g=0 sign pattern with xi∗(k)≈−1x_i^*(k)≈-1, and the component drifts continuously as g increases. The unstable equilibrium branch (red) approaches the stable one and annihilates with it near g≃0.4g 0.4. In green (orange) is presented the stable (unstable) path computed using our perturbative method. The depletion of SFPs is due to a sequence of tangent (also called saddle-node or fold) bifurcations. In fact, for g=0g=0 there exists a much larger number of unstable fixed points (UFPs) y→∗(k)(0)\ y^*(k)(0)\, usually saddles of different order, corresponding to the (3N−2N)(3^N-2^N) length-N sequences of 0s and ±1± 1s, containing at least one 0. For each SFP there exist N neighbouring UFP’s, obtained by setting one of its components equal to 0. Analogously to the SFPs, saddles can be continued analytically starting from g=0g=0. Upon increasing g, it may happen that one SFP x→∗(k)(g) x^*(k)(g) approaches one of the neighbouring saddles y→∗(l)(g) y^*(l)(g) and eventually the two points mutually annihilate for some SFP-dependent critical value gcg_c. This scenario is clearly illustrated in Fig. 2, where one component of an SFP is reported as a function of g. The blue line corresponds to the SFP; beyond gc≈0.4g_c≈ 0.4 an initial condition in the vicinity of the no-longer existing SFP converges towards another SFP, which inherits the basin of attraction of the disappeared SFP. The disappearance of the SFP could be equally monitored by following any component of the SFP, since the bifurcation is a general phenomenon, which implies that all components of the two colliding fixed points simultaneously collapse with one another. We conjecture that the bifurcation is driven by a single component, the most sensitive to the coupling. The bifurcation mechanism is well captured by a perturbative approach. Given an SFP at g=0g=0 with components xj∗(0)∈±1Nx^*_j(0)∈\± 1\^N, we introduce the induced field (see Eq.(2) ) Si=1N∑j=1NAijsj(k)S_i\;=\; 1 N _j=1^NA_ij\,s^(k)_j, and thereby approximate the dynamics of the i-th component as x˙i=xi(1−xi2)+gSi x_i=x_i(1-x_i^2)+g\,S_i. In this approximation, the evolution still factorizes into the independent relaxation of the various components within quartic potentials. For a given g, an SFP exists so long as its components correspond to a local minimum. Since a positive (negative) SiS_i tends to destabilize a negative (positive) si(k)s^(k)_i, it is convenient to introduce ui=Sisi(k)u_i=S_is^(k)_i (for the sake of simplicity we drop the index k). Let us now denote with i the component of the SFP characterized by the most negative uiu_i and be u its value. Given u, it is readily seen that the critical value where such a minimum disappears is gc=−233/2u,g_c\;=\;- 23^3/2\,u, (3) The accuracy of this approach can be appreciated in Fig. 2, where the j-th component (therein i=1i=1) of an SFP determined via the perturbative approach (see the green line) can be compared with the actual exact solution. The two curves are very close to one another up to the critical point, in spite of gcg_c being not too small. The advantage of this approximate approach is that one can infer the range of existence of the various SFPs directly from the knowledge of the matrix A, without the need of performing simulations for different coupling strengths. In fact, if we denote with ρN(u) _N(u) the normalized empirical density of u values (averaged over disorder realizations), we can express the fraction PN(g)P_N(g) of SFPs for a given value of g, as PN(g)=∫uc(g)+∞ρN(u)uP_N(g)= _u_c(g)^+∞ _N(u)\,du where uc(g)=−2/(33/2g)u_c(g)=-2/(3^3/2g) is obtained by inverting Eq. (3). Figure 3: Main panel: the symbols denote numerical estimates of ρN(u) _N(u) for N=20N=20. Data have been obtained by averaging over 1000 realizations of the matrix A. Full circles, crosses, and squares correspond to the Ginibre ensemble, a uniform distribution of entries, and to the symmetric case, respectively. The solid curves identify the corresponding theoretical predictions. Inset: the triangles correspond to h(u)h(u) estimated as lnρN/N+ln2 _N/N+ 2 for the Ginibre ensemble. The black and red solid curves correspond to the theoretical predictions for asymmetric and symmetric matrices, respectively. The resulting fractions PNP_N correspond to the dashed line in Fig. 1: they agree pretty well with the direct numerical data over several decades, confirming the validity of the perturbative scheme. Since ρN(u) _N(u) contains the relevant information to characterize the stability of the various SFPs, we now focus on its structure. The numerical estimates for N=20N=20 are reported in Fig. 3 (see full circles). Notice that positive u values are irrelevant in the context of the stability analysis, since they correspond to SFPs which become more, rather than less, stable, under the action of the mutual coupling. Moreover, given the inverse proportionality between u and g, the very negative u values identify the most fragile SFPs: those which disappear for very small coupling strengths. Given the matrix A, u is the minimum value within a set of N variables SiS_i, each one obtained by summing N independent elements, multiplied by de facto random independent signs (we sum over all SFPs and hence all sign patterns). Therefore, we can invoke order-statistics identities for generic distributions (see I.1). In particular, by denoting with φ , the probability density function (PDF) of the SiS_i elements and with Φ the corresponding cumulative distribution function (CDF), it is known Ross (2020) that ρN(u)=Nφ(u)(1−Φ(u))N−1. _N(u)=N\, (u)\, (1- (u) )^N-1. (4) In the present case, ϕ(u)φ(u) is the unit-variance Gaussian. There is only a doubt on the validity of general theorems. Since the minimum is taken after multiplying the variables SiS_i by the pattern of signs employed in the definition of SiS_i itself, we cannot exclude subtle correlations sneak in. Hence, we have directly tested the validity of Eq. (4). The comparison can be appreciated in Fig. 3: the theoretical prediction is indistinguishable from the direct numerical results. The same is true for other values of N (data not shown). It is natural to invoke the large-deviation theory and conjecture that lnρN _N is asymptotically proportional to N; in mathematical terms, ψ(u)=limN→∞(lnρN(u))/Nψ(u)= _N→∞( _N(u))/N. This is indeed correct, and from Eq. (4), we see that ψ(u)=ln(1−Φ(u))ψ(u)= (1- (u)). In fact, it is more instructive to refer to the actual (average) number EN(u)=2NPN(u)E_N(u)=2^NP_N(u) of expected SFPs, rather than to their fraction. Their growth rate is h(u)=ln2+ψ(u)=ln(2(1−Φ(u)))h(u)= 2+ψ(u)= (2(1- (u))). The resulting distribution is shown in the inset of Fig. 3 (see the black solid curve). The triangles are obtained directly from (lnρN)/N( _N)/N (see the triangles): the initial rising part is a clearcut finite-size effect. Figure 4: Comparison between the perturbative theoretical prediction h(g)=ln2+ln[1−Φ(uc(g))]h(g)= 2+ [1- (u_c(g))] (black solid line) and the rate extracted from direct numerical simulations in the large-N limit (red squares with error bars). Finally, we have plotted the multiplicity index h as a function of the coupling strength g (via Eq. (3)) to allow for a comparison with direct numerical simulations. The solid curve in Fig. 4 corresponds to the perturbative estimate of h(g)h(g), which is strictly larger than 0 for arbitrarily large g values. This is a consequence of the implicit assumption that the only way SFPs disappear is via tangent bifurcations, which let anyhow N SFPs survive. This is not true, as we have spotted various homoclinic and heteroclinic bifurcations, which lead to either the disappearance or the destabilization of some SFPs (eventually to the onset of limit cycles and chaos). Whether or not these additional mechanisms can lead to a vanishing or even possibly negative h(g)h(g) for a finite coupling strength is hard to say. Direct numerical estimates of h(g)h(g), obtained via a fit in the range N∈[10,25]N∈[10,25], are in fact consistently smaller than the theoretical prediction (see the red squares) and the difference becomes substantial for g larger than 0.5. However, it is appropriate to notice that numerical data are unavoidably affected by finite-size corrections, quite difficult to quantify. Different ensembles - If the matrix entries do not follow a Gaussian distribution, the central-limit theorem anyhow implies that, in the large N limit, u should still be distributed in a Gaussian way, though with likely deviations in the tails. The relevance of such deviations, can be appreciated in Fig. 3 (see the blue crosses), where we report data obtained by assuming a uniform distribution in the interval [−3,3][- 3, 3] (still unit variance) and N=20N=20. Appreciable differences can be seen only for u smaller than −4-4 (i.e. g<0.1g<0.1). Symmetric matrices represent an important class of models, characterized by a strictly gradient dynamics, which implies that the evolution can only converge onto an SFP. Since the symmetry Aij=AjiA_ij=A_ji introduces strong correlations among the matrix entries, we do not expect our theory to reproduce exactly the multiplicity of SFPs. This is confirmed by the refined perturbative arguments developed in I.2. However, as shown in Fig. 3), the numerical data obtained for N=20N=20 (see the red squares) do not differ significantly from those for the Ginibre ensemble. However, an important qualitative discrepancy can be appreciated in the inset, where one sees that the exponential growth rate h(u)h(u) remains strictly finite for u→0u→ 0 (i.e. g→∞g→∞). Since in the large-g limit, there is no guarantee that the perturbative approach provides sufficiently accurate results, we have performed direct simulations of the full model for increasing N and g=5g=5. The data reported at the end of I.2 reveal a clear exponential growth, while no evidence can be found in the asymmetric case, where periodic cycles and chaotic attractors eat a large fraction of the phase space. Conclusions and perspectives - In this Letter we have focused our attention on the survival mechanism of stable fixed points in model (1), showing that a suitable perturbative approach allows for reliable estimates, making use of large-deviation theory. We have also shown that such a perturbative approach is effective for different classes of random coupling matrices, including symmetric ones. Preliminary results (to appear in a forthcoming publication) indicate that the perturbative method is effective also when the random coupling matrix contains planted attractors or has been passed through a training procedure for solving specific computational tasks. In these cases, the matrix entries cannot be assumed anymore to be i.i.d. variables and the standard order-statistics (see I.1) does not apply. In fact, the probability distributions ϕ(x)φ(x) obtained from ensembles of the above mentioned matrices typically acquire long-tails, testifying at their intrinsic non-Gaussian nature and to the presence of correlations among the matrix elements. However, we can argue that the perturbative method still applies to sum of variables and this is why it works pretty well, although providing less accurate estimates of the survival probability of SFP’s. In the limit of large coupling strengths, preliminary studies show the appearance of limit cycles and even chaotic attractors (see also I.3). Their relation with the survival of SFPs and with the presence of unstable fixed points is quite an intricate problem. In order to tackle it successfully, statistical and dynamical concepts and tools have to be suitably combined. For instance, in a more elaborated model than (1) it has been found that in the ferromagnetic and paramagnetic chaotic phases there is no direct correspondence between the presence of an exponentially large number of unstable fixed points and the attractor manifold of these phases Fournier et al. (2026). References L. Chicchi, D. Fanelli, D. Febbe, L. Buffoni, F. Di Patti, L. Giambagli, and R. Marino (2025) Deterministic versus stochastic dynamical classifiers: opposing random adversarial attacks with noise. 6 (3), p. 035054. External Links: Document, Link Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. S. J. Fournier, A. Pacco, V. Ros, and P. Urbani (2026) Nonreciprocal interactions and high-dimensional chaos: comparing dynamics and statistics of equilibria in a solvable class of models. 113, p. 044139. External Links: Document, Link Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. S. Gagliani, F. Giuseppe Pacifico, L. Chicchi, D. Fanelli, D. Febbe, L. Buffoni, and R. Marino (2026) Train stochastic non linear coupled odes to classify and generate. 7 (2), p. 025020. External Links: Document, Link Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. F. Ghimenti, A. Sriram, A. Yamamura, H. Mabuchi, and S. Ganguli (2026) Geometry and dynamics of annealed optimization in the coherent ising machine with hidden and planted solutions. Phys. Rev. E 113, p. 054123. External Links: Document, Link Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. R. Kim, Y. Li, and T. J. Sejnowski (2019) Simple framework for constructing functional spiking recurrent neural networks. Proceedings of the National Academy of Sciences 116 (45), p. 22811–22820. External Links: Document, Link, https://w.pnas.org/doi/pdf/10.1073/pnas.1905926116 Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. R. Marino, L. Buffoni, L. Chicchi, L. Giambagli, and D. Fanelli (2024) Stable attractors for neural networks classification via ordinary differential equations (sa-node). Machine Learning: Science and TechnologyMachine Learning: Science and TechnologyPhys. Rev. EMachine Learning: Science and TechnologyNeural Computation 5 (3), p. 035087. External Links: Document, Link Cited by: Multiplicity of Stable Attractors in Disordered Neural Models, Multiplicity of Stable Attractors in Disordered Neural Models. R. Marino, L. Buffoni, L. Chicchi, F. D. Patti, D. Febbe, L. Giambagli, and D. Fanelli (2025) Learning in wilson-cowan model for metapopulation. 37 (4), p. 701–741. External Links: ISSN 0899-7667, Document, Link, https://direct.mit.edu/neco/article-pdf/37/4/701/2506385/neco_a_01744.pdf Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. S. M. Ross (2020) A first course in probability. 10th edition, Pearson Benelux. External Links: ISBN 978-1292269207 Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. N. Rungratsameetaweemana, R. Kim, T. Chotibut, and T. J. Sejnowski (2025) Random noise promotes slow heterogeneous synaptic dynamics important for robust working memory computation. Proceedings of the National Academy of Sciences 122 (3), p. e2316745122. External Links: Document, Link, https://w.pnas.org/doi/pdf/10.1073/pnas.2316745122 Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. M. Syed and N. G. Berloff (2023) Physics-enhanced bifurcation optimisers: all you need is a canonical complex network. IEEE Journal of Selected Topics in Quantum Electronics 29 (2: Optical Computing), p. 1–6. External Links: Document Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. M. Syed, R. Z. Wang, and N. G. Berloff (2026) Soft vector spins with dimensional annealing for combinatorial optimization. arXiv preprint arXiv:2604.01003. Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. G. Wainrib and J. Touboul (2013) Topological and dynamical complexity of random neural networks. Phys. Rev. Lett. 110, p. 118101. External Links: Document, Link Cited by: Multiplicity of Stable Attractors in Disordered Neural Models. I END MATTER I.1 Mean and variance of the minima distribution Here, we summarize analytical formulae about the distribution of the minima u≡min1≤i≤NSiu≡ _1≤ i≤ NS_i (5) that enters the perturbative scheme. The key point is that, for a stable fixed point at g=0g=0, the SiS_i are i.i.d. Gaussian variables. As a consequence, u is the minimum of N i.i.d. standard normal variables, whose distribution, mean, and variance are controlled by classical order-statistics and extreme-value theory. I.1.1 From Ginibre to i.i.d. Gaussians Let ∈ℝN×NG ^N× N be a real Ginibre matrix with Aij∼i.i.d.(0,1)A_ij .i.d. N(0,1), i.e., =/NG=A/ N. For a stable fixed point at g=0g=0 we set s→∈±1N s∈\± 1\^N and define Si≡∑j=1NGijsj=1N∑j=1NAijsj,u≡min1≤i≤NSi.S_i\;≡\; _j=1^NG_ijs_j\;=\; 1 N _j=1^NA_ijs_j, u\;≡\; _1≤ i≤ NS_i. (6) For each fixed i, SiS_i is a linear combination of independent Gaussians, hence Gaussian. Moreover, [Si]=0,Var(Si)=1N∑j=1NVar(Aij)=δii,E[S_i]=0, (S_i)= 1N _j=1^NVar(A_ij)= _i, (7) so that Si∼(0,1).S_i (0,1). (8) Since different rows of A are independent, Sii=1N\S_i\_i=1^N are independent as well. Therefore, for any stable fixed point s→ s, the random variable u is the minimum of N i.i.d. standard normal variables. I.1.2 Exact finite-N distribution of u Let Φ and φ denote the CDF and PDF of the standard normal. By independence, the probability of the minimum in a sample of N i.i.d. random variables from a standard normal is (u<x)=1−(1−Φ(x))N.P(u<x)=1- (1- (x) )^N. (9) Hence, the CDF and PDF of u are FN(u) F_N(u) =(u≤x)=1−(1−Φ(u))N, =P(u≤ x)=1- (1- (u) )^N, (10) ρN(u) _N(u) =FN′(u)=Nφ(x)(1−Φ(x))N−1. =F_N (u)=N\, (x)\, (1- (x) )^N-1. (11) These are standard order-statistics identities. I.2 Symmetric Matrices Here, we sketch the extension to the case of symmetric matrices of the perturbative criterion discussed in the manuscript for asymmetric ones. Let A be an N×N× N symmetric Gaussian matrix with off-diagonal entries Aij=Aji∼(0,1)A_ij=A_ji (0,1) (i≠ji≠ j) and diagonal variance d=Var(Aii)d=Var(A_i). The two conventions of interest are d=0(Aii=0)d=0\ (A_i=0), and d=1(Aii∼(0,1))d=1\ (A_i (0,1)). For a stable fixed point at g=0g=0 we set s→∈±1N s∈\± 1\^N. Let s=diag(s1,…,sN)D_s=diag(s_1,…,s_N) and ~=ss A=D_s\,A\,D_s, i.e. A~ij=siAijsj A_ij=s_iA_ijs_j. The map ↦~A A leaves the symmetric Gaussian ensemble invariant (gauge invariance), and we have ui=siSi=1N∑jA~ij.u_i=s_iS_i= 1 N _j A_ij. (12) Hence, the joint law of ui\u_i\ does not depend on the pattern s→ s, which also justifies averaging over all 2N2^N patterns. In this case, the only source of correlation is the symmetry constraint A~ik=A~ki A_ik= A_ki. For i≠ki≠ k, the double sum Cov(ui,uk)=1N∑j,l[A~ijA~kl]Cov(u_i,u_k)= 1N _j,lE[ A_ij A_kl] receives a single nonvanishing contribution (j=k,l=ij=k,\ l=i), so that Cov(ui,uk)=Cij=1N(i≠k),Var(ui)=Cii=(N−1)⋅1+dN=1+d−1N. split&Cov(u_i,u_k)=C_ij= 1N (i≠ k), \\ &Var(u_i)=C_i= (N-1)· 1+dN=1+ d-1N. split (13) In a more compact form we can write Cij=σN2δij+1N,σN2=1+d−2N,C_ij= _N^2\, _ij+ 1N, _N^2=1+ d-2N, (14) which admits the exact representation ui=σNvi+z/Nu_i= _Nv_i+z/ N with v1,…,vN,zv_1,…,v_N,z i.i.d. (0,1)N(0,1). Since the common shift z/Nz/ N factors out of the minimum mN=minivim_N= _iv_i, we have u=dσNmN+zN,z⟂mN.u d= _N\,m_N+ z N, z m_N. (15) By convolving Eq. (15), we obtain the probability density function (PDF) ρNsym(u)=∫−∞+∞zN2πe−Nz2/21σNρN(u−zσN). _N^sym(u)= _-∞^+∞\!dz\, N2π\,e^-Nz^2/2\, 1 _N\, _N\! ( u-z _N ).\\ (16) The corresponding cumulative density function (CDF) reads P(u≤x)=1−∫−∞+∞zN2πe−Nz2/2(1−Φ(x−zσN))N.P(u≤ x)=1- _-∞^+∞\!dz\, N2π\,e^-Nz^2/2\, (1- \! ( x-z _N ) )^N. (17) Since the dominant common-mode fluctuation is z=cNz=c N, with c=O(1)c=O(1), the diagonal convention drops out and the growth rate of the expected number of surviving SFPs, h(u)=limN→∞1Nln(2NρN(u))h(u)= _N→∞ 1N \! (2^N _N(u) ), becomes hsym(u)=ln2+maxc≥0[ln(1−Φ(u−c))−c22],h_sym(u)= 2+ _c≥ 0 [ (1- (u-c) .)- c^22 ], (18) which coincides in the limit u→−∞u→-∞ with the result ln[2(1−Φ(u))] [2(1- (u))] obtained for asymmetric matrices. On the other hand, at variance with the asymmetric case, hsym(u)h_sym(u) does not vanish at u=0u=0 (see the inset of Fig. 3 in the manuscript): hsym(0)=ln2+lnΦ(c∗)−c∗22=0.1992…,φ(c∗)=c∗Φ(c∗)⇒c∗≃0.5061. split&h_sym(0)= 2+ (c^*)- c^*22=0.1992…,\\ & (c^*)=c^* (c^*)\ \ c^* 0.5061. split (19) Note that these analytic formulae hold in the limit N→∞N→∞. In order to appreciate the reliability of the perturbative criterion, in Fig. 5 we report an example of the exponential growth of the number of SFPs of model (1) with symmetric coupling matrices for g=5g=5, i.e. a coupling value quite far from the perturbative regime. Taking into account that this numerical estimate of the growth rate, 0.159, is unavoidably affected by hard-to-quantify finite size effects, it is remarkable its closeness to the theoretical prediction 0.1992. Figure 5: Number of stable fixed points (SFP) as a function of the system size N, shown on a linear–log scale. Blue markers denote the mean SFP averaged over disorder realizations, with error bars. The red line is a least-squares fit that confirms an exponential growth SFP∼e0.159NSFP e^0.159N. I.3 Large-g limit In order to analyze dynamics (1) in the limit of large values of g it is worth rescaling xi→gxix_i→ gx_i and t→gt→ gt, thus obtaining dxidt=−xi(xi2−1g)+1N∑j=1NAijxj,i=1,2,⋯,N dx_idt=-x_i(x_i^2- 1g)+ 1 N\, _j=1^N\,A_ijx_j,\,i=1,2,·s,N\, (20) In the limit of g→+∞g→+∞, while inducing also the asymptotic limit t→+∞t→+∞, this set of equations simplifies to dxidt=−xi3+1N∑j=1NAijxj,i=1,2,⋯,N dx_idt=-x_i^3+ 1 N\, _j=1^N\,A_ijx_j,\,i=1,2,·s,N\, (21) If the coupling matrix AijA_ij would be diagonal with real eigenvalues λii=1N\ _i\_i=1^N these equations represent a set of N Duffing oscillators, with stable fixed point in 0 if λi<0 _i<0 or in ±λi± _i if λi>0 _i>0. In the case of asymmetric coupling matrices we have found numerically that dynamics (21) yields an asymptotic evolution converging to a low-dimensional chaotic attractor (data not shown).