How extensive is chaos in random nonlinear recurrent neural networks?

Why the whole spectrum matters, not just the top rate

A chaotic system separates nearby trajectories exponentially, and the rate of separation depends on the direction of the initial displacement. The largest such rate, the maximum Lyapunov exponent, governs the late-time separation of a generic displacement. The full set of rates, the Lyapunov spectrum, describes how an infinitesimal ball of initial conditions is stretched and squeezed along every direction, so it carries far more information than the single largest exponent.

The paper frames its target as an open problem posed nearly four decades earlier for the random network of Sompolinsky et al.: evaluating the whole distribution of Lyapunov exponents, not only the maximum. The network studied is the Euler discretization, with time step 0<δ≤1, of that model, x_i(n+1) = (1−δ)x_i(n) + δ Σ_j J_ij φ(x_j(n)), with couplings J_ij drawn i.i.d. from a Gaussian of variance g²/N and φ(x) = tanh(x). The coupling strength g and the step δ are the two parameters that set the dynamics.

Two derived quantities make the spectrum worth computing. The sum λ_1+…+λ_k is the growth rate of k-dimensional volumes in tangent space, which feeds the attractor dimension; the spectrum also gives access to the entropy rate. Because the spectrum is invariant under diffeomorphisms of phase space, these quantities are intrinsic properties of the dynamics rather than artifacts of a coordinate choice.

Counting exponents below a threshold as a response

The central object is the cumulative distribution F(s) = (1/N)·#{i : λ_i < s}, the fraction of Lyapunov exponents below a threshold s. The paper's first move is an exact identity at finite N that turns this counting problem into a response function of the tangent-space dynamics.

The trick is to shift every exponent by s. Multiplying the Jacobian M(n) at every step by e^{−δs} lowers every Lyapunov exponent by s, giving a shifted Jacobian M_s(n) = α_s I_N + β_s J D(n) with α_s = (1−δ)e^{−δs} and β_s = δ e^{−δs}, where D(n) holds the gains d_i(n) = φ′(x_i(n)). Under this shifted dynamics, directions with λ_i < s decay forward in time and directions with λ_i > s grow, so N·F(s) is exactly the number of directions that decay forward.

To count those directions, the shifted dynamics is driven by a source and written as a block-bidiagonal linear system K V = I_v over all integer times. That system imposes no initial condition, so its solutions form an affine family; the paper removes the freedom by selecting the solution of minimum norm. The minimum-norm solution is the one that decays in both time directions, and it is built explicitly by splitting the source at time k+1 into stable and unstable components using the Oseledets splitting of tangent space. The resulting one-step response of this solution then yields F(s).

A cavity calculation reduces the response to one site

The finite-N identity alone does not give a number, because the response still depends on the full network. The second step is a cavity method that determines this response at large N through a self-consistent single-site problem. One neuron, labeled 0, is treated as a cavity surrounded by the rest of the network, and its input ξ_0(n) becomes a zero-mean Gaussian process whose covariance is fixed by the self-consistency condition ⟨ξ_0(n)ξ_0(n′)⟩ = g²⟨φ(x_0(n))φ(x_0(n′))⟩.

Solving that condition fixes the statistics of the gain trajectory d_0(n) = φ′(x_0(n)), which then drives forward and backward fields for neuron 0. The paper writes these fields as a 2×2 block system and inverts it by blocks to obtain the single-site responses R^{vu}_{00}, R^{uv}_{00}, and R^{vv}_{00}. Two kernels, A(n,n′) and P(n,n′), are then determined self-consistently from averages over the single-site trajectory, with the gains appearing inside the average for A but outside it for P.

The cumulative distribution is read out as F(s) = lim_{η→0+} ⟨R^{vv}_{00}(k+1,k)⟩, where η is a regularization that is taken to zero. The threshold s enters only through α_s and γ_s, and in the stationary state the kernels depend only on the time lag n−n′. The paper states its main assumption plainly: the derivation assumes the large-N limit and the vanishing-regularization limit can be exchanged, and a more mathematical treatment would need to justify this exchange along with other steps of the cavity calculation.

The theory reproduces two known limits

The calculation is checked against cases where the answer is already known. For g<1 the fixed point x_i = 0 is stable, every gain equals φ′(0) = 1, and the Jacobian is constant, so the Lyapunov exponents follow from the eigenvalues μ_i of J through λ_i = (1/δ)log|1−δ+δμ_i|. Because the eigenvalues of J fill the disk |μ|<g by the circular law as N→∞, F(s) becomes the fraction of that disk lying inside a second disk |μ−c_δ|<r_s, with c_δ = 1−1/δ and r_s = e^{δs}/δ.

The single-site theory reproduces this geometric answer for every δ. The paper shows the agreement by an electrostatic interpretation: the readout integral equals the flux of a two-dimensional electric field through the circle |μ−c_δ|=r_s, and by Gauss's law that flux counts the charge enclosed, which is exactly the overlap area of the two disks. As δ→0 the second disk approaches the half-plane Re μ < 1+s, recovering the semicircle law; at δ=1 the exponents reduce to log|μ_i|.

The maximum exponent is recovered too. The upper edge of the spectrum in the single-site theory is a threshold s*, and the paper shows that the solution with vanishing kernels exists for s>s* but cannot be continued below it, so the theory gives λ_1 = s* for every δ. This recovers the earlier results for δ→0 and for δ=1. The paper reports that the resulting theory agrees closely with simulations, and that together with recent calculations of the participation-ratio dimension, all quantities in its summary figure can now be computed from theory.

Extensive chaos means the spectrum scales with network size

The headline claim is that the chaos in this network is extensive. In this setting that means the Lyapunov spectrum is a nontrivial distribution over an N-independent fraction of directions as N grows, so the count of exponents below any threshold s is N·F(s) with F(s) a well-defined function in the N→∞ limit. Chaos is not confined to a handful of directions; a finite fraction of the network's degrees of freedom carries exponents of each sign.

This is why the cumulative distribution, rather than the maximum exponent alone, is the informative object. A single λ_1 says how fast a generic perturbation separates, but it does not say how many directions expand, how many contract, or how the rates are distributed. The paper's F(s) answers those questions directly, and the derived attractor dimension and entropy rate inherit that directional detail.

The paper also notes what the spectrum buys conceptually: because it is diffeomorphism-invariant, quantities derived from it describe the dynamics itself rather than a particular coordinate representation. That is the sense in which the calculation gives access to intrinsic properties of the dynamics, not just to a numerical summary of one simulation.

Where the calculation extends and where it stops

The derivation is specific: a rate network with i.i.d. Gaussian asymmetric couplings, tanh nonlinearity, Euler time stepping, and a single ergodic invariant measure assumed for the dynamics. The paper states that the approach extends straightforwardly to rate networks with structure beyond i.i.d. couplings, naming several populations such as separate excitatory and inhibitory groups, couplings with a low-rank component, and correlations between reciprocal couplings J_ij and J_ji. It also flags spiking networks as an extension one could attempt, without claiming a result there.

Beyond recurrent networks, the paper suggests computing Lyapunov spectra for the reasoning dynamics of looped transformers and for the dynamics of weights during learning. These are stated as directions the method might reach, not as findings of this paper, and should be read that way.

The main formal limitation sits with the method rather than the conclusion: the exchange of the large-N and vanishing-regularization limits is assumed, not proven, and the author notes that a more mathematical treatment would be needed to justify it and the other cavity steps. That caveat qualifies the derivation's rigor, while the agreement with the trivial-fixed-point spectrum, the recovered maximum exponent, and the reported simulation agreement are the evidence supporting the extensive-chaos conclusion.

The spectrum, not the top exponent, is the deliverable

The paper's distinctive contribution is not a new value for λ_1 but a closed route to the entire distribution F(s) for this random network, obtained by combining an exact finite-N identity with a self-consistent single-site problem. The most decision-relevant comparison is between what the maximum exponent alone tells you and what the full spectrum adds: the maximum sets the late-time separation rate of a generic displacement, while F(s) determines how many directions expand or contract at each rate and therefore fixes the volume-growth sums behind the attractor dimension and entropy rate.

For a reader deciding whether this result applies to a given network, the operative variables are the coupling strength g, the time step δ, and the i.i.d. Gaussian asymmetric coupling assumption. The theory is anchored where those assumptions hold: it matches the circular-law spectrum at the trivial fixed point for every δ, recovers the known maximum exponent in the δ→0 and δ=1 limits, and is reported to agree closely with simulations. Structured couplings, spiking dynamics, and transformer or learning dynamics lie outside what this paper establishes, and the unproven exchange of limits is the formal gap to keep in view when weighing how far the extensive-chaos conclusion can be pushed.

Sources used in this answer