Lyapunov spectrum of random neural networks
David G. Clark
cond-mat.dis-nn, q-bio.NC
2026-10-09
Derives the full Lyapunov spectrum of random recurrent networks at N→∞; theory matches N=4096 simulations and proves chaos is extensive. GPT-6 produced the initial derivation in 100 minutes.
In 1988, Sompolinsky, Crisanti and Sompolinsky introduced the random recurrent neural network: N neurons coupled through an unstructured, asymmetric Gaussian matrix, each low-pass filtering its input and applying a tanh nonlinearity. Past a critical coupling g = 1 the network enters a chaotic phase, which has served for decades as the standard model of spontaneous cortical activity and as the substrate for trained RNNs. Their paper left one problem open: evaluate the whole distribution of Lyapunov exponents, the rates at which infinitesimal perturbations grow or shrink along different directions. In the 38 years since, the maximum exponent yielded to analysis (a Schrödinger-type ground-state energy in continuous time), while the full spectrum existed only as numerical QR computations on large networks. This paper computes the full spectrum in the N → ∞ limit.
The spectrum matters because the largest exponent only says whether the dynamics are chaotic. The full set gives the attractor dimension (Kaplan–Yorke) and the entropy rate (sum of positive exponents), both invariant under smooth coordinate changes, so they describe the dynamics rather than the coordinates you picked.
Three moves.
First, an identity exact at finite N. Counting exponents below a threshold s is the same as counting decaying directions after every exponent is shifted down by s. Instead of evolving tangent vectors from an initial condition, the paper drives the shifted tangent dynamics with a source on a doubly infinite time window and selects the minimum-norm solution. That solution grows along the unstable subspace before the source and decays along the stable subspace after it; it is the only bounded solution, and its one-step response matrix is exactly the projector onto the stable subspace. The trace of a projector is the dimension of its range, so the normalized trace gives the cumulative distribution F(s).
Second, the minimum-norm solution is written as the zero-regularization limit of a regularized least-squares problem, whose optimality condition forms a forward–backward system: an N-dimensional field running forward in time, another running backward, coupled by the regularizer η, which also guarantees uniqueness. The construction is a cousin of Hermitization from non-Hermitian random-matrix theory.
Third, the cavity method (add one neuron to a large reservoir and solve for its effective environment self-consistently) is applied jointly to the network and the forward–backward system. At large N the added neuron sees a Gaussian cavity field with statistics fixed by dynamical mean-field theory, plus two self-consistent kernels; its gain trajectory, the slope of the nonlinearity along its activity, drives the two fields. The single-site problem is iterated numerically, η is lowered to zero, and F(s) is read off.
One load-bearing assumption: the limit N → ∞ commutes with η → 0. The paper states plainly that this exchange is unproven, supported by agreement with simulations and by precedent from random-matrix theory.
The theory first reproduces two known answers, extending both from the endpoints δ → 0 and δ = 1 to arbitrary time step δ: the circular-law spectrum at the trivial fixed point, and the maximum exponent (Schrödinger ground-state formula in continuous time, log g·C^d(0) at δ = 1), now unified as one nonlinear eigenvalue problem at finite δ.
Against simulations of N = 4096 networks using the standard QR method, the single-site theory matches the full spectrum closely for g = 3 and g = 5 across time steps δ from 0.05 to 0.5.
| Comparison | Before | This paper |
| Spectrum at fixed point | circular law, only δ→0 and δ=1 | extended to all δ |
| Maximum exponent | two separate limit formulas | one finite-δ eigenvalue problem |
| Full spectrum | numerical only | analytic large-N theory, matches N=4096 |
Two consequences follow directly. The chaos is extensive: attractor dimension per neuron and entropy rate per neuron converge to self-averaging values, and the number of positive exponents scales with N, turning Sompolinsky's mean-field argument and Engelken's numerics into a derivation. The attractor dimension rises with g, saturates, and at δ = 0.5 peaks near g = 10 before declining; the entropy rate increases monotonically over the range computed.
For theoretical neuroscience, this closes a marquee open problem. Combined with earlier participation-ratio calculations, every quantity in the phase diagram now comes from theory where large simulations used to be the only route.
For machine learning, the connection is shorter than it looks. The order-to-chaos transition at initialization in deep networks, pushing Lyapunov exponents toward zero to learn long-time dependencies in RNNs, and transiently chaotic reasoning dynamics in looped transformers all currently rely on single exponents or numerical spectra. An analytic spectrum is a sharper instrument for these analyses.
For research practice, the paper is an unusually transparent specimen of frontier models doing research-grade mathematics. The author had assumed the problem was analytically intractable; the target was proposed by colleague Litwin-Kumar. GPT-6 Astra produced the initial derivation working independently for 100 minutes inside a session with prior related work. The author then worked with GPT-6 and Claude Opus 5.5 to strip the original's heavy machinery (Deninger's theorem, a Fuglede–Kadison log-determinant) down to the minimum-norm argument, and Claude Opus 5.5 wrote all code, generated all figures, and checked the derivation numerically. The AI methodology section reproduces the prompt verbatim.
Stated by the author: the central exchange of limits has no mathematical justification, and a rigorous treatment is, in the author's words, "beyond the expertise of the author"; formalization in a proof assistant is suggested as future work.
Scope is narrow: i.i.d. Gaussian couplings, tanh nonlinearity, rate neurons. Multiple populations, low-rank couplings, correlated reciprocal couplings, and spiking networks are listed as extensions, not results. The final step is a numerical iteration of the single-site problem, so this is an analytic reduction rather than a closed form. Validation covers g = 3 and g = 5 with δ from 0.05 to 0.5; extremes of coupling and time step are not shown. And the AI-produced derivation is verified the way this field always verifies: recovery of known limits plus agreement with simulation, with no independent third-party audit.