Disordered Dynamics in High Dimensions: Connections to Random Matrices and Machine Learning

TL;DR

This work presents a pedagogicalDMFT framework for analyzing high‑dimensional, disordered dynamical systems with random matrices, focusing on learning and generalization processes in machine learning. By deriving and solving DMFT equations via cavity and path‑integral methods, it recovers central spectral laws (e.g., Wigner semicircle and Marchenko–Pastur) and characterizes loss dynamics through two‑time correlation and response functions, revealing nontrivial phenomena like non‑monotone losses and multistage learning. It applies the formalism to linear regression, kernel methods, random feature models, online SGD, and evolving networks, highlighting how bias/variance decompositions and non‑Hermitian dynamics arise from random projections and structured covariates. The framework exposes how spectral content, finite‑sample effects, and non‑normality shape training dynamics, interpolation thresholds, and generalization, offering a principled route to quantify learning curves in high‑dimensional random data. Overall, the DMFT approach provides a unifying lens to connect random matrix theory with practical learning dynamics, with potential for broad applicability beyond the linear/Gaussian settings considered here.

Abstract

We provide an overview of high dimensional dynamical systems driven by random matrices, focusing on applications to simple models of learning and generalization in machine learning theory. Using both cavity method arguments and path integrals, we review how the behavior of a coupled infinite dimensional system can be characterized as a stochastic process for each single site of the system. We provide a pedagogical treatment of dynamical mean field theory (DMFT), a framework that can be flexibly applied to these settings. The DMFT single site stochastic process is fully characterized by a set of (two-time) correlation and response functions. For linear time-invariant systems, we illustrate connections between random matrix resolvents and the DMFT response. We demonstrate applications of these ideas to machine learning models such as gradient flow, stochastic gradient descent on random feature models and deep linear networks in the feature learning regime trained on random data. We demonstrate how bias and variance decompositions (analysis of ensembling/bagging etc) can be computed by averaging over subsets of the DMFT noise variables. From our formalism we also investigate how linear systems driven with random non-Hermitian matrices (such as random feature models) can exhibit non-monotonic loss curves with training time, while Hermitian matrices with the matching spectra do not, highlighting a different mechanism for non-monotonicity than small eigenvalues causing instability to label noise. Lastly, we provide asymptotic descriptions of the training and test loss dynamics for randomly initialized deep linear neural networks trained in the feature learning regime with high-dimensional random data. In this case, the time translation invariance structure is lost and the hidden layer weights are characterized as spiked random matrices.

Figures (14)

• Figure 1: Cavity derivation of the marginal dynamics for a single site of the system as $N \to \infty$. Adding a new site to the system comes with $N$ reciprocal couplings $\bm a_0 \in \mathbb{R}^N$ to the original variables which are now perturbed $\bm h(t) \to \tilde{\bm h}(t)$. In the large system size limit $N \to \infty$ the system can be viewed as a single-variable stochastic process driven by a colored noise process with a delayed feedback through response function $R(t,t')$.
• Figure 2: The DMFT response function for the random Wigner matrix encodes the semicircle eigenvalue law. (a) The spectral distribution $\rho(\lambda)$ for a randomly sampled $N = 8000$ Wigner matrix (red) is compared to the asymptotic theoretical density $\rho(\lambda) = \frac{1}{2\pi} \sqrt{4-\lambda^2}$ (dashed black lines). (b) The response function $R(\tau)$ as a function of the time lag $\tau$ for the dynamics $\frac{d}{dt} \bm h(t) = - \bm M \bm h(t) - z \bm h(t)$ where $\bm M$ is a Wigner matrix and $z=2$. For this system, the relaxation rate at $z=2$ is powerlaw with $R(\tau) = \frac{1}{2\pi} \int e^{-( \lambda + 2)\tau}\sqrt{4-\lambda^2 } \propto \tau^{-3/2}$ for large $\tau$.
• Figure 3: The cavity method for the linear regression problem can proceed in two steps. First, a computation of the marginals for the weight discrepancy $h_0 = w_0 - w_0^*$ when a new $N+1$st feature is added requires considering feedback through the perturbed training errors $\tilde{\Delta}_\mu(t)$. Second, the training error made on an added $P+1$st data point $\Delta_0(t)$ requires considering feedback through the perturbed weight discrepancies $\tilde{\bm h}(t)$. Under the joint limit $N,P \to \infty$ with $P = \alpha N$, the loss
• Figure 4: The dynamics of linear regression with random dataset of side $P = \alpha N$ are governed by the DMFT response function $R(\tau)$ which encodes the spectral properties of a Wishart matrix. Experiments with $N=1000$ are shown in solid lines while the DMFT is plotted in black dashed lines. (a) For $\alpha < 1$, the response function $R(\tau)$ saturates to $1-\alpha$ at large time lag $\tau \to \infty$. Alternatively if $\alpha > 1$, the response function at large $\tau$ relaxes exponentially with timescale set by the minimum eigenvalue $R(\tau) \sim \exp\left(- \left[1 - \alpha^{-1/2} \right]^2 \tau \right)$. (b) From the Fourier transform of the response $\mathcal{R}(\omega)$, we can recover the eigenvalue density $\rho(\lambda) = \frac{1}{\pi} \ \text{Im} \mathcal{R}( i\lambda -\epsilon)$ which we plot without the singularity at $\lambda= 0$. (c) The dynamics of gradient flow accurately describe the effect of subsampled data in the proportional regime $P/ N = \alpha$.
• Figure 5: The effective densities $\rho_k(z)$ and corresponding response functions $\mathcal{H}_k(\tau)$ for structured kernel regression with eigenvalues $\lambda_k = k^{-b}$ with $b=1.25$. We plot the first $8$ eigemodes $k \in [10]$ for $P \in \{64, 512\}$. (a) The effective densities $\rho_k(z)$ exhibit larger spread for smaller $P$. As $P \to \infty$ these converge to Dirac masses $\rho_k(z) \to \delta(z-\lambda_k)$. (b) The response functions $\mathcal{H}_k(\tau)= \int dz e^{-z\tau} \rho_k(z)$ as a function of the timelag $\tau$. For small $\tau$, these functions scale as $\mathcal{H}_k(\tau) \sim e^{-\lambda_k \tau}$, while for large $\tau$ they relax to a constant.