Predictable Mean-Field Chaos in Random Recurrent Neural Networks

TL;DR: Take the canonical chaotic random RNN — positive Lyapunov exponent, trajectories that decorrelate in O(1) time — and ask: given the exact continuous past of a single neuron, and nothing else, can you predict its entire future? For a broad class of smooth activation functions (erf-like, but notably not tanh or ReLU), the answer is yes, exactly, with zero error. Chaos doesn't vanish; it gets converted into a resolution cost. The information needed to forecast is pushed into ever-higher time derivatives of the observed trajectory at a rate \alpha that the authors show is a different quantity from the Lyapunov exponent \lambda_L, with different critical scaling. Resolving p orthogonalized derivative-modes buys you a prediction horizon \sim \log p. This cleanly separates two things that are usually conflated — microscopic instability and irreducible unpredictability — and it turns the power spectrum of neural activity into a diagnostic for whether observed variability is genuine noise or hidden deterministic structure.

The standard story, and the loophole in it

The Sompolinsky–Crisanti–Sommers (SCS) model is the hydrogen atom of chaotic neural dynamics: N rate neurons, \dot{x}_i = -x_i + \sum_j J_{ij}\phi(x_j), with i.i.d. Gaussian couplings J_{ij} \sim \mathcal{N}(0, g^2/N). Above a critical gain g_c, the network is deterministically chaotic. In the N\to\infty limit, dynamical mean-field theory (DMFT) replaces the whole network with a single effective neuron driven by a self-consistent stationary Gaussian process: \dot{x} = -x + \eta(t), where \eta's statistics are fixed by the network's own autocorrelation \Delta(\tau) = \langle x(t)x(t+\tau)\rangle.

The universal interpretation of this reduction: once you've integrated out the other N-1 neurons and the connectivity, \eta(t) is noise. It's an ensemble description. One neuron's history can't tell you its future, because the future depends on everything you threw away.

The loophole is that "Gaussian process" is not synonymous with "noise-driven." Wiener, Kolmogorov, and Krein classified stationary Gaussian processes into two types. Regular processes continually generate innovations — new randomness that no amount of past observation can anticipate (Brownian motion, Ornstein–Uhlenbeck). Singular processes are random across realizations, but any single realization's continuous past determines its entire future exactly. The dividing line is the Paley–Wiener criterion on the power spectral density \tilde{\Delta}(\omega):

Intuitively: if the spectrum decays exponentially or faster at high frequencies, the log-integral diverges and the process is singular. An exponentially-decaying spectrum means the trajectory is an analytic function of time — and an analytic function is globally determined by its behavior on any interval, the same way analytic continuation works for complex functions. All the "future" information is sitting in the infinite tower of time derivatives at the present moment.

Crucially, this is logically independent of the Lyapunov exponent. \lambda_L > 0 measures how fast nearby but distinct microscopic network states diverge under perturbation. Singularity asks about the conditional uncertainty of one realized trajectory given its exact past — and if the past is exact, there is no perturbation to amplify. The paper's contribution is to show that the chaotic SCS mean-field process actually lands on the singular side, to trace exactly when and why, and to quantify what "exact past" costs at finite resolution.

The mechanism: activation smoothness → complex-time singularity → spectral tail

The chain of logic runs through complex analysis, and it's genuinely pretty. DMFT gives a closed equation for the autocorrelation: \Delta(\tau) behaves like a classical particle in an effective potential, \ddot{\Delta} = -\partial_\Delta V(\Delta), where V is built from Gaussian expectations of the antiderivative of \phi. The high-frequency tail of the spectrum \tilde{\Delta}(\omega) is controlled by how far \Delta(\tau) can be analytically continued into complex time — specifically, by the distance t_S from the real axis to the nearest singularity: \tilde{\Delta}(\omega) \sim e^{-t_S \omega}.

And where that singularity sits is dictated by the Fourier decay of the activation function:

• Gaussian-or-faster Fourier decay (e.g. \phi = \operatorname{erf}(x/\sqrt{2}\sigma)): the potential V(\Delta) has its nearest branch point at S = \Delta_0 + \sigma^2 — strictly outside the physically visited range |\Delta| \le \Delta_0. The particle takes a finite imaginary time t_S > 0 to reach it, the spectrum decays as e^{-t_S\omega}, Paley–Wiener diverges: perfectly predictable.
• tanh (exponential Fourier decay, from its poles at \pm i\pi/2): the singularity gets pinned exactly at the physical boundary S = \Delta_0. The result is a stretched-exponential spectral tail — smooth trajectories, infinitely differentiable, but not analytic. Paley–Wiener converges: innovations persist.
• ReLU / hard-tanh (finite smoothness): algebraic spectral tail, unpredictability fails outright.

Convolving any non-smooth activation with a Gaussian of width \sigma restores predictability — but as \sigma \to 0 the analytic strip collapses and the practical prediction scale diverges, which is a nice consistency check.

Organizing the determinism: the Krylov chain

Perfect predictability from an infinite derivative tower is a mathematician's victory, not an engineer's. Raw derivatives X^{(n)}(0) are a horribly conditioned basis. The paper's second move is to orthogonalize them via a Lanczos construction on the power spectrum: treat \tilde{\Delta}(\omega)\,d\omega as a spectral measure, build its orthonormal polynomials \phi_n(\omega), and note that "multiply by \omega" is tridiagonal in this basis. This yields a semi-infinite Jacobi matrix J with hopping coefficients \alpha_n — a 1-D chain.

The payoff is a state-space realization of the mean-field process. Define Krylov coordinates s_n(t) = \phi_n(i\partial_t)X(t) from the trajectory. Then |s(t)\rangle evolves under the deterministic linear equation \frac{d}{dt}|s\rangle = -iJ|s\rangle, the observed neuron is just the boundary coordinate X(t) = \operatorname{Re}\langle e_0|s(t)\rangle, and — for the Gaussian process — the initial coordinates s_n(0) are i.i.d. standard Gaussians. All the randomness has been relocated to the initial condition; the dynamics are noiseless. Prediction becomes: read off the Krylov state from the observed past, evolve linearly.

For exponential spectral tails, the Lanczos coefficients grow linearly, \alpha_n \sim \alpha n, with \alpha = \pi/(2t_S). This \alpha is the paper's central new quantity: the rate at which predictive information is pumped into higher temporal modes (equivalently, the growth rate of Krylov complexity — readers of the quantum operator-growth literature will recognize the machinery of Parker et al.).

Now compare \alpha to \lambda_L near the chaos transition. For smooth odd nonlinearities, the near-critical kernel is \Delta_0\operatorname{sech}(\tau/\tau_\infty), locking \alpha = 1/\tau_\infty. Since \tau_\infty \sim (g-g_c)^{-1}, we get \alpha \sim (g - g_c) — while the classic result is \lambda_L \sim (g-g_c)^2. Different exponents: these are parametrically distinct quantities, not two views of one "chaos rate."

Across all regularized nonlinearities and parameter ranges they tested, \alpha > \lambda_L — analytically near the transition, numerically away from it. Intriguingly, this sits inside the bound \lambda_L \le 2\alpha proved by Parker et al. for quantum Hamiltonian operator dynamics, even though this system is dissipative and non-Hamiltonian and the bound doesn't formally apply. The authors flag this as suggestive, not established.

The chaos tax: horizon \sim \log p

Here's where chaos re-enters, in its proper place. Suppose you can only resolve the first p Krylov modes (equivalently, roughly p orthogonalized derivatives of the past). The unresolved tail acts as a bath — a Mori–Zwanzig-style split, with the tail integrated out as a self-energy and approximated by an OU closure. Uncertainty then comes only from the unknown initial tail coordinates.

A WKB analysis of the chain gives the picture: an excitation launched at the boundary site e_0 propagates down the chain with an accelerating front, k_{\rm front}(t) \sim e^{2\alpha t}, because the hopping \alpha_n \sim \alpha n grows with depth. Your forecast stays accurate until the information front reaches the unresolved region at site p:

This is the honest reconciliation with chaos: to double your prediction horizon you need to square the number of resolved modes. Exponential resolution demand is exactly what a positive Lyapunov exponent should cost you — except the relevant exponent is \alpha, not \lambda_L.

Does it survive contact with a finite network?

The most important sanity check is Section S11 / Figure 3: they simulate an actual SCS network with N = 10^4 neurons, g = 1.5, \phi = \operatorname{erf}, observe one neuron over a past window of T_{\rm past} = 3 time units sampled at \delta t = 0.005 (M = 601 samples), estimate Krylov coordinates by regularized least squares (fitting K = 100 coordinates, rank-30 truncated SVD, propagating p = 20), and forecast forward — knowing nothing about the connectivity and observing none of the other 9,999 neurons. The forecasts track the true continuation well past the decorrelation time, and the empirical error matches the analytical error of knowing roughly the first 12 Krylov coordinates exactly. A matched OU process with the same correlation time is decisively worse, which is the key control: the win comes from the derivative structure of the past, not from long correlations.

That the estimation needed rank-truncated SVD on an ill-conditioned design matrix is itself informative: the higher Krylov coordinates are exponentially hard to read from noisy, discretely sampled data, which is the finite-precision shadow of the \log p law.

What this changes, and what to be skeptical of

The conceptual shift is real. "Chaotic" and "unpredictable from its own past" are routinely treated as synonyms, in the RNN literature and in interpretations of neural variability. This paper shows they factor apart, and gives each its own rate with its own critical exponent. It also hands you an operational diagnostic: measure the high-frequency tail of a unit's power spectrum. Exponential tail → the apparent variability is deterministic structure recoverable from that unit's own history; power-law or stretched-exponential tail → genuine innovations. For neuroscience, where "is trial-to-trial variability noise or unobserved deterministic dynamics?" is a live and consequential question, that's a genuinely usable idea — with the massive caveat that real neurons have intrinsic noise and inputs, and any additive white noise makes a process regular in the Wiener sense. Section S10 addresses this partially: with weak white noise the spectrum splits into deterministic and stochastic components, and the Krylov construction still bounds the reducible part of the uncertainty. But this is perturbative; how gracefully the horizon degrades with realistic noise levels is not quantified.

Other things to keep skeptical about. The perfect-predictability class excludes tanh and ReLU — the two activations everyone actually uses — and the fix (Gaussian smoothing) has a prediction scale that blows up as \alpha \sim \sigma^{-1} when smoothing vanishes, so the result is about a class of models, not a universal RNN property. The ordering \alpha > \lambda_L is proved only near criticality for smooth odd nonlinearities; elsewhere it's numerics on a few activation families. The finite-network validation is one parameter setting, and the estimator's effective resolution (~12 modes) times a logarithmic horizon means the practical forecasting gain over naive baselines, while clear, is bounded — this is a theory paper whose "practical" section is a proof of concept, not a forecasting tool. Finally, everything rides on the DMFT limit; finite-N corrections reintroduce interactions the effective Gaussian process ignores, and the paper doesn't chart how predictability degrades with N.

Where to spend your time

Read Section II for the Paley–Wiener argument and the activation-function trichotomy — that's the core "aha" and it's compact. Section III is where the new quantity \alpha is defined and the \alpha-vs-\lambda_L scaling separation is derived; it's the part most likely to seed follow-up work (does \lambda_L \le 2\alpha hold generally in dissipative mean-field systems?). If you care about applicability, Supplement S10 (noise) and S11 (finite-network forecasting protocol) are where the idealizations meet reality, and where the next paper's open problems live.