The pixels of a recurrent weight matrix almost refuse to tell a story. A low-rank term is spread across all \(N^2\) couplings, so even as its strength grows, any one entry can remain smaller than the random texture around it. The animation below follows the controlled experiment in this post: the random realization and structured direction stay fixed while \(\rho\) rises. Coarse-graining makes the planted geometry visible on the left, the static spectrum lets its eigenvalue escape the bulk in the middle, and the driven diagnostics on the right show where that clean matrix story begins to fray.

What I wanted to keep from Rainer's movie was its visual grammar rather than its data: a matrix, its spectrum, and one low-dimensional trace sharing the same clock. Here that clock is \(\rho\). No network is being trained, no neurons are reclustered between frames, and the random bulk and structured vectors do not change. The low-rank signal first becomes visible because sorting and averaging let it add coherently while the random component cancels; the spectrum is always computed before that display-only compression.

The difficulty is hidden in the phrase the network's spectrum. A driven nonlinear recurrent network offers several spectra, and they answer different questions. The eigenvalues of its connectivity describe a static operator. The eigenvalues of an averaged Jacobian describe a trajectory after its changing gain mask has been compressed into one mean. A scalar outlier compresses that averaged spectrum again. Stability, meanwhile, belongs to the ordered sequence of Jacobians encountered along the trajectory. Low-rank structure can survive all the way through the average without the stability problem collapsing to that structure.

The rest of the post follows this diagram from left to right. The first step is clean finite-rank random-matrix theory. The second produces two cheap scalar approximations, one structured and one global. The full averaged spectrum then recovers much of what those scalars discard. The final time-ordered step explains why even that stronger agreement remains a diagnostic rather than a stability theorem.

I. When Structure Leaves the Bulk

Consider a continuous-time recurrent rate network with currents \(h_i(t)\), \[ \tau\frac{dh_i}{dt} =-h_i+\sum_{j=1}^N J_{ij}\phi(h_j)+bI(t), \] where \(\phi(x)=\max(0,x)\) and \(I(t)\) is a common Ornstein–Uhlenbeck drive. The reported experiments set independent additive noise to zero and use \(\tau=1\). Their connectivity is \[ J=gW-bqq^\top+\rho uv^\top, \qquad q=\frac{1}{\sqrt N}\mathbf 1, \] with \(W_{ij}\sim\mathcal N(0,1/N)\). Each term gives the matrix a different geometry. The random matrix \(gW\) produces a circular eigenvalue bulk of radius approximately \(g\); the mean term contributes a real balance outlier near \(-b\); and the last term inserts a structured direction with strength \(\rho\), right vector \(u\), left vector \(v\), and alignment \[ \alpha=v^\top u. \]

The balance outlier is a useful precedent. In a tightly balanced network, recurrent inhibition cancels most of the common input in the population mean, leaving a fast mean mode whose characteristic timescale shrinks with \(b\). An eigenvalue near \(-b\) then has a fairly direct dynamical interpretation. It is tempting to expect the added low-rank outlier to provide the same kind of handle, but that expectation is justified first for the matrix, not yet for the driven network.

For the static connectivity, the calculation is clean. Write a finite-rank perturbation as \(S=UMV^\top\). The matrix determinant lemma gives \[ \det\!\left(zI-(gW+S)\right) =\det(zI-gW) \det\!\left(I-MV^\top(zI-gW)^{-1}U\right). \] Outside the random bulk, the projected resolvent has leading term \(z^{-1}V^\top U\), so a detached rank-one candidate lies at \[ z_{\mathrm{spike}}=\rho\alpha. \] The sign matters: a positive aligned spike produces an outlier on the positive real axis, while the balance mode lies at \(-b\) because its signed amplitude is negative. The formula also comes with a geometric boundary. When \(|\rho\alpha|\leq g\), the candidate lies in the Ginibre disk and should not be described as a detached outlier.

At this level, the ghost behaves as hoped. A high-dimensional random matrix accepts a low-rank perturbation, and one eigenvalue carries the perturbation's signed overlap out of the bulk. The next step is harder because the network does not evolve by repeatedly applying \(J\).

II. When the Network Moves, the Matrix Changes

Along a nonlinear trajectory, a small perturbation \(\delta h(t)\) evolves according to \[ \tau\frac{d\,\delta h}{dt}=A(t)\delta h, \qquad A(t)=-I+JD(t), \] where \[ D(t)=\operatorname{diag}\!\left(\phi'(h_1(t)),\ldots,\phi'(h_N(t))\right). \] For ReLU, \(D(t)\) is a changing binary gate. Active neurons pass their columns of \(J\) into the linearization; inactive neurons remove theirs. The structured mode in \(J\) is therefore filtered differently at different times, and an eigenvector of the connectivity need not remain an eigenvector of the instantaneous Jacobian.

One natural compression is to replace the moving gate by its trajectory average, \[ \bar D=\langle D(t)\rangle_t, \qquad A_{\mathrm{avg}}=-I+J\bar D. \] This matrix remembers which neurons were active often and forgets the order in which they were active. That distinction is easy to overlook because \(A_{\mathrm{avg}}\) has an ordinary spectrum and admits the same outlier machinery as \(J\). It is nevertheless a different operator, with a bulk shifted by the leak and rescaled by heterogeneous gains.

The actual perturbation after time \(T\) is instead \[ \delta h(T)=\mathcal T\exp\!\left( \frac{1}{\tau}\int_0^T A(t)\,dt \right)\delta h(0), \] where \(\mathcal T\) preserves temporal order. Its asymptotic growth rate is the top Lyapunov exponent, \[ \lambda_1 =\lim_{T\to\infty}\frac{1}{T} \log\left\|\mathcal T\exp\!\left( \frac{1}{\tau}\int_0^T A(t)\,dt \right)\right\|. \] A positive value means nearby trajectories separate exponentially; a negative value means perturbations decay. This is the stability verdict for the driven system. Neither an eigenvalue of \(J\) nor an eigenvalue of \(A_{\mathrm{avg}}\) inherits that verdict automatically.

III. A Cheap Outlier for an Averaged World

Although \(A_{\mathrm{avg}}\) cannot settle stability, its low-rank outlier is still a useful intermediate object. Separate the unstructured base from the rank-one term, \[ B=-I+(gW-bqq^\top)\bar D, \qquad A_{\mathrm{avg}}=B+\rho u(v^\top\bar D). \] Applying the determinant identity again gives an exact scalar equation for an outlier \(z\), \[ 1=\rho v^\top\bar D\,(zI-B)^{-1}u. \] The object has already changed. The static overlap \(\alpha=v^\top u\) has become the trajectory-dependent contraction \(v^\top\bar D u\), while the random balanced base enters through the resolvent \((zI-B)^{-1}\). The trajectory has filtered the structure, but the filtered term remains rank one.

If we now replace \(B\) by the leak \(-I\), the resolvent becomes \((z+1)^{-1}I\), and the equation collapses to the cheap scalar proxy \[ \hat z_{\mathrm{scalar}}=-1+\rho v^\top\bar D u. \] It costs only an inner product once a trajectory has supplied \(\bar D\), and it preserves how the mean gain mask aligns with the structured vectors. That alignment sounds informative, but it needs a baseline. The still cheaper alternative is \[ \hat z_{\mathrm{mean}}=-1+\rho\alpha\bar p, \qquad \bar p=\frac{1}{N}\operatorname{tr}\bar D, \] which retains only the population's active fraction. For a delocalized aligned mode that is weakly correlated with per-neuron gain, these two quantities can be almost indistinguishable.

IV. The Scalar Shortcut Breaks First

The corrected one-dimensional sweep follows the aligned mode through every level of that ladder. We fix \(N=500\), \(g=1.6\), \(b=5\), and \(\alpha=1\), then vary \(\rho\) from \(0\) to \(4\) across eight paired seeds. Pairing matters here: each seed reuses the same random bulk, structured vectors, initial condition, and drive increments across \(\rho\), so changes along a curve are not blurred by a fresh network at every point. The static spike detaches only after \(\rho>g\), and the overlap-tracked mode of \(A_{\mathrm{avg}}\) remains weakly identified until roughly \(\rho=2.4\). Before then, it is a bulk-tracking heuristic rather than a separated eigenmode.

Across all 88 paired samples, the structured scalar moves with the corrected finite-time Lyapunov estimate, with correlation \(0.79\). The agreement is real, but it is not distinctive. The global mean-gain baseline has lower RMSE, \(0.22\) rather than \(0.26\), and fewer positive-sign disagreements, about \(3\%\) rather than \(9\%\). In this aligned experiment, retaining the directional contraction \(v^\top\bar D u\) does not improve on knowing the active fraction alone.

The full averaged operator retains considerably more. Its spectral abscissa \[ s(A_{\mathrm{avg}})=\max_i\operatorname{Re}\lambda_i(A_{\mathrm{avg}}) \] has correlation \(0.96\) with the finite-time Lyapunov estimate, RMSE \(0.08\), and one sign disagreement across the sweep. At \(\rho=4\), the seed-mean scalar has moved to \(0.19\), the right edge of the averaged spectrum sits near \(0.00\), and the finite-time Lyapunov estimate remains at \(-0.12\). The paired bias \(\lambda_1-\hat z_{\mathrm{scalar}}\) is \(-0.31\), with a 95% interval from \(-0.47\) to \(-0.15\). Even the smaller residual \(\lambda_1-s(A_{\mathrm{avg}})\) remains negative there. Because the averaged spectral estimate is itself close to zero with non-negligible uncertainty, this comparison establishes a calibration gap, not the location of a boundary.

A longer-window check preserves the mean calibration gap. At \(g=1.6\) and \(\rho=4\), extending the measurement to \(60\tau\) gives \(\lambda_1=-0.146\pm0.071\) SEM, compared with a scalar proxy near \(0.20\) and a tracked averaged-Jacobian outlier near \(-0.04\). Halving the integration step gives the same ordering. Two of the eight seeds remain positive, so this is a persistent bias in the seed mean rather than a claim that every realization is stable.

This changes the story of the experiment. The structured mode survives the nonlinear gain average, yet compressing the whole averaged spectrum to that mode discards information that the spectral abscissa still carries. The first important loss does not wait for temporal order to be removed; it begins when one outlier is asked to stand in for the rest of \(A_{\mathrm{avg}}\).

We can describe that loss more carefully by separating two residuals, \[ \Delta_{\mathrm{scalar}} =\operatorname{Re}z_{\mathrm{avg}}-\hat z_{\mathrm{scalar}}, \qquad \Delta_{\mathrm{dyn}} =\lambda_1-\operatorname{Re}z_{\mathrm{avg}}. \] The scalarization gap records what is lost when the base resolvent is replaced by the leak. The dynamical gap is everything left after the averaged-Jacobian outlier has been computed: temporal ordering, gain fluctuations, products of non-normal matrices, finite-time and finite-\(N\) effects, and any remaining tracking error. At \(\rho=4\), their mean magnitudes are approximately \(0.19\) and \(0.12\), and neither dominates uniformly across the sweep. The second residual is descriptive; it is not a causal measurement of temporal ordering alone.

A corrected \(9\times9\) sweep over \((g,\rho)\) preserves the same ordering without producing a phase boundary. The mean finite-time Lyapunov estimate remains negative in all 81 sampled cells, while the cheap summaries become positive in part of the low-\(g\), high-\(\rho\) corner. Across 405 paired seed–grid-point samples, the structured scalar has correlation \(0.70\) with \(\lambda_1\) and a positive-sign disagreement rate near \(7\%\); the mean-gain baseline is again slightly better. That broader sweep is useful corroboration, but turning it into another large figure would suggest a measured transition that the data do not contain.

V. What Averaging Forgets

The strong agreement between \(s(A_{\mathrm{avg}})\) and \(\lambda_1\) on the one-dimensional slice says that averaging is informative in this regime. It does not make the two objects identical, and the mathematical reason appears even in a trajectory with only two gain masks. Suppose \(D_1\) acts for one short interval and \(D_2\) for the next, producing Jacobians \(A_1\) and \(A_2\). The perturbation evolves as \[ e^{\Delta t A_2/\tau}e^{\Delta t A_1/\tau}. \] The Baker–Campbell–Hausdorff expansion rewrites this product as \[ \exp\!\left( \frac{\Delta t}{\tau}(A_1+A_2) +\frac{\Delta t^2}{2\tau^2}[A_2,A_1] +\cdots \right), \] where \([A_2,A_1]=A_2A_1-A_1A_2\). Averaging keeps the first term and discards the commutators. If the Jacobians commute, the missing terms vanish. In a random recurrent network with changing ReLU gates they generally do not, so the route taken through the active subspaces can survive even after their mean has been fixed.

Non-normality adds a related distinction between eventual decay and temporary growth. The spectral abscissa \(s(A)\) tracks the rightmost eigenvalue, while the numerical abscissa \[ \omega(A)=\lambda_{\max}\!\left(\frac{A+A^\top}{2}\right) \] asks whether some direction initially grows in Euclidean norm. A matrix can have \(s(A)<0\) and \(\omega(A)>0\), so perturbations may grow transiently even when every eigenmode decays under that fixed operator. A moving gain mask compounds this geometry by changing the directions available for amplification. The mere fact that the full spectral edge lies to the right of a tracked outlier does not diagnose non-normality; that claim would require the numerical-abscissa comparison.

The diagnostic ladder is now easier to read. The structured scalar preserves one contraction \(v^\top\bar D u\). The full averaged Jacobian preserves the entire mean gain profile and, on this slice, most of the Lyapunov trend. Both have forgotten when those gates acted. Low-rank structure survives averaging; stability does not collapse to that structure, and any agreement with the averaged spectrum remains empirical rather than automatic.

VI. What the Ghost Can Tell Us

The scalar proxy has a narrower role after these corrections, but it still has one. During a large parameter sweep or a training run, computing a full averaged spectrum or a long Lyapunov estimate at every checkpoint can be prohibitive. The contractions \(v^\top\bar D u\) and \(\alpha\bar p\) are cheap once a representative trajectory has been collected. In this aligned model the global baseline is at least as informative, so the structured contraction should be treated as one candidate summary rather than the privileged coordinate of stability.

What emerges is a diagnostic policy rather than a new theorem. Track the static outlier to verify that the intended low-rank structure exists. Track the scalar and mean-gain contractions to see how the trajectory's average gain filters it. Compute the spectrum of \(A_{\mathrm{avg}}\) when those scalar approximations need calibration. Use a time-ordered Lyapunov estimate when the sign of stability matters. The cheap objects decide where to look; the last object decides what the dynamics do.

The opening parameter sweep suggests how this policy might be used during training. If a learned mode begins to separate from the bulk, it could trigger targeted trajectory and Lyapunov diagnostics, although the controlled experiment here does not establish that an outlier predicts a learning transition. The resemblance to the staged emergence of modes in deep linear networks is suggestive, but for now it remains an analogy rather than a mechanism.

A full low-rank non-stationary dynamical mean-field theory would try to close this gap instead of averaging it away. It would retain low-dimensional overlaps for the structured modes, couple them self-consistently to the time-dependent Gaussian bulk, and derive a reduced time-ordered stability problem. That theory remains unfinished. For now, the outlier is best understood as a ghost in the machine in a literal mathematical sense: it makes hidden structure visible, while the averaged spectrum and the tangent flow show that it does not act alone.

This project began in a graduate course on dynamical systems, neural networks, and theoretical neuroscience. I am grateful to Rainer Engelken for the conversations that shaped it, including the training movie whose visual grammar inspired the opening animation.

References

1. R. Engelken and S. Goedeke, “A non-stationary dynamic mean field theory for strongly interacting neuronal networks,” NeurIPS 2022.
2. F. Mastrogiuseppe and S. Ostojic, “Linking connectivity, dynamics, and computations in low-rank recurrent neural networks,” Neuron 2018.
3. T. Tao, “Outliers in the spectrum of iid matrices with bounded rank perturbations,” Probability Theory and Related Fields 2013.
4. G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, “Lyapunov characteristic exponents for smooth dynamical systems,” Meccanica 1980.
5. A. M. Saxe, J. L. McClelland, and S. Ganguli, “Exact solutions to the nonlinear dynamics of learning in deep linear neural networks,” arXiv:1312.6120.