2602.01434

# ArXiv Paper 2602.01434v1 ## Phase Transitions for Feature Learning in Neural Networks Andrea Montanari ∗ Zihao Wang † February 3, 2026 ## Abstract According to a popular viewpoint, neural networks learn from data by first identifying effective low-dimensional representations, and subsequently fitting the best model in this lowdimensional space. A sequence of recent works provides a rigorous formalization of this phenomenon when learning multi-index models. In this setting, we are given n i.i.d. pairs of covariate vectors and responses ( x i , y i ), where the vectors x i ∈ R d are isotropic, and responses y i only depend on x i through a k -dimensional projection Θ T ∗ x i . Feature learning amounts to learning the latent space spanned by Θ ∗ . In this context, we study the gradient descent dynamics of two-layer neural networks under the proportional asymptotics n, d → ∞ , n/d → δ ∈ (0 , ∞ ), while the dimension of the latent space k and the number of hidden neurons m are kept fixed. Earlier work establishes that feature learning via polynomial-time algorithms is possible if δ > δ alg , for δ alg a threshold depending on the data distribution, and is impossible (within a certain class of efficient algorithms) below δ alg . Here we derive an analogous threshold δ NN for a class of two-layer networks. Our characterization of δ NN is sufficiently explicit to open the way to study the dependence of learning dynamics on the target function, loss, activation, network width, and initialization. The threshold δ NN is determined by the following training scenario. Gradient descent first visits points for which the gradient of the empirical risk is large and learns the directions spanned by these gradients. Then the gradient becomes smaller and the dynamics becomes dominated by negative directions of the Hessian. The threshold δ NN corresponds to a phase transition in the spectrum of the Hessian after the first phase of training, and provides a quantitative explanation for empirical phenomena such as grokking. ## Contents | 1 | Introduction | 3 | |-----|-----------------------------------|------------| | | 1.1 Main results: An informal | . . . . 6 | | | 1.2 Related literature . . | . . . . 10 | | | 1.3 Summary of notations | . . . . 12 | | 2 | Results 12 | Results 12 | | | 2.1 Problem setup and assumptions | . . . . 12 | | | 2.2 Preliminaries: Discrete | . . . . 13 | | | 2.3 Result I: Hard directions | . . . . 14 | | | 2.4 Result II: Bulk of the | . . . . 15 | | | 2.5 Main results: Outliers | . . . . 16 | ∗ Department of Statistics and Department of Mathematics, Stanford University † Department of Mathematics, Stanford University | 3 Consequences and discussion of main results | 3 Consequences and discussion of main results | 19 | |-------------------------------------------------|-----------------------------------------------------------------------------------------------|------| | 3.1 | The single-neuron case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 19 | | 3.2 | Grokking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 19 | | 3.3 | Interpretation of the outliers equation . . . . . . . . . . . . . . . . . . . . . . . . . . | 20 | | 3.4 | Suboptimality among polynomial-time algorithms . . . . . . . . . . . . . . . . . . . . | 21 | | 4 Examples and numerical illustrations | 4 Examples and numerical illustrations | 22 | | 4.1 | General setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 22 | | 4.2 | Results for GeLU and Quad activations . . . . . . . . . . . . . . . . . . . . . . . . . . | 23 | | 4.3 | Training dynamics and grokking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 25 | | 5 Outline of the proofs | 5 Outline of the proofs | 26 | | 5.1 | Step 1: Gaussian conditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 26 | | 5.2 | Step 2: Outlier matrix and concentration to deterministic limit . . . . . . . . . . . . | 27 | | 5.3 | Step 3: Eigenvector alignment via residue calculation . . . . . . . . . . . . . . . . . . | 29 | | 6 Discussions | and conclusions | 30 | | A Technical | tools | 35 | | A.1 | Rouch´ e's theorem for random functions . . . . . . . . . . . . . . . . . . . . . . . . . | 35 | | A.2 | From pointwise convergence to uniform convergence . . . . . . . . . . . . . . . . . . | 36 | | A.3 | Asymptotics of minimum eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . | 36 | | A.4 | Concentration inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 38 | | B Proofs | for Sections 2.3 and 2.4 | 38 | | B.1 | Proofs of Lemma 2.1 and Lemma 2.2: Discrete time DMFT . . . . . . . . . . . . . . | 38 | | B.2 | Proof of Lemma 2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 40 | | B.3 | Proof of Lemma 2.4: Bulk of the Hessian . . . . . . . . . . . . . . . . . . . . . . . . | 41 | | C Proofs | of Theorems 1 and 2 | 43 | | C.1 | Outliers via resolvent analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 43 | | C.2 | Resolvent concentration via leave-one-out . . . . . . . . . . . . . . . . . . . . . . . . | 47 | | C.3 | Convergence of the outlier matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 53 | | C.4 | Eigenvalue localization via determinant zeros . . . . . . . . . . . . . . . . . . . . . . | 57 | | C.5 | Weak recovery of hard directions from eigenvectors . . . . . . . . . . . . . . . . . . . | 59 | | C.6 | Proof of Theorem 1: Finite time spectral phase transition . . . . . . . . . . . . . . . | 62 | | C.7 | Proof of Theorem 2: Large time limit . . . . . . . . . . . . . . . . . . . . . . . . . . | 62 | | D Proofs | of supporting lemmas | 66 | | D.1 | Proof of Lemma C.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 66 | | D.2 | Proof of Lemma C.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 67 | | D.3 | Proof of Lemma C.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 68 | | D.4 | Proof of Lemma C.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 69 | | D.5 | Proof of Lemma C.15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 70 | | D.6 | Proof of Lemma C.19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . | 71 | | E Additional | experiments | 72 | | E.1 | Results for networks with moderately large width . . . . . . . . . . . . . . . . . . . . | 72 | | E.2 | Grokking across activation functions . . . . . . . . . . . . . . . . . . . . . . . . . . . | 74 | ## 1 Introduction Pressed to explain the success of modern deep learning systems, many practitioners might mention the ability of these systems to 'learn low-dimensional representation of the data,' during training. This is often contrasted to linear methods (e.g. kernel methods or neural networks under 'lazy training') which work in a fixed representation that is determined by a human's choice of the kernel or featurization map. An important line of recent work has provided a formalization of this intuition for the problem of learning multi-index models, see [BES + 22, DLS22, ABAM23, BMZ24, ZWFL25] for a few pointers to this literature. This is a classical supervised learning problem. We observe n i.i.d. samples ( x i , y i ), i ≤ n with d -dimensional covariates x i which we assume don't have low-dimensional structure: x i ∼ N ( 0 , I d ). The responses y i depend on a k -dimensional projection of the covariates: $$y _ { i } = h ( \Theta ^ { T } x _ { i } , e _ { i } )$$ Here Θ ∗ ∈ R d × k has orthonormal columns, noise ε i is independent of x i , and h : R k +1 → R is a link function. We will attempt to learn a predictive model using a two-layer neural network: $$f ( x ) := - \sum _ { j = 1 } ^ { m } a _ { j } g ( b _ { j } + b _ { j } ) , \theta : = [ \theta _ { 1 } \ldots \theta _ { m } ] \in y ^ { x + b _ { j } }$$ In order to focus on the feature-learning aspects of this problem and avoid inessential technicalities, we will only train the first-layer weights Θ and and regard ( a j , b j ) m j =1 as fixed. We also emphasize that the number of neurons m is in general unrelated to the latent dimension k . Throughout, we will regard k as constants independent of n, d →∞ . We will eventually let m →∞ after n, d →∞ . Given a loss function ℓ ( · , · ) : R × R → R , we are interested in the dynamics of gradient descent (GD) on the empirical risk $$R is k ( \theta ) := - 1 _ { n = 1 } ^ { n } \sum _ { i = 1 } ^ { n } l ( y _ { i } , f e ( x _ { i } ) ) .$$ The neural network (2) is a natural model for the multi-index target distribution. Indeed any distribution (1) (more precisely, the conditional expectation E [ y | x ] = h E ( Θ T ∗ x )) can be approximated by a network of the form (2), with ∥ a ∥ 1 /m ≤ R 0 (for R 0 a sufficiently large constant), resulting in generalization error bounded by CR 0 √ d/n [Bar96, Bac17]. Vice versa, any such neural network is well approximated by a k -index model with k = O ( R 0 ∨ 1). Since the generalization error is O ( R 0 √ d/n ), learning by empirical risk minimization succeeds with n = O ( d ) samples. Hence, the key question is whether GD (with respect to the empirical risk (3)) can achieve the same and whether it is optimal among algorithms with comparable complexity. Despite a copious literature on this model, no sharp results exist about this question. In order to explain the issue, consider for the moment the case of a single-index model, k = 1, and assume h ( z ) = h ( z, ε ) is a generic even 1 function of z . Letting n, d →∞ with δ n := n/d → δ , then there exists an estimator ̂ θ = ̂ θ ( y , X ) with |⟨ ̂ θ , θ 0 ⟩| / ∥ ̂ θ ∥ bounded away from zero (i.e. it achieves 'weak recovery') for any δ > δ IT , while such an estimator does not exist for δ < δ IT . The 'information theoretic' (or 'statistical') δ IT = δ IT ( h ) was characterized in [BKM + 19]. 1 If h is not even, then for any δ > 0 the latent direction θ 0 can be learned with non-trivial accuracy by linear regression, after preprocessing the y i 's. Figure 1: Empirical success probabilities and predicted thresholds for learning a single-index model with link function h ( z, ε ) = z 2 (noiseless phase retrieval). The model is a single-neuron network with GeLU activation, trained via gradient descent; we declare success when the correlation between the learned feature and the target direction exceeds 1 / 2. Orange curves correspond to optimal spectral initialization, blue curves to random initialization; within each color family, lighter to darker shades indicate increasing dimensions d = 1000 , 2000 , 3000 , 4000. Dots are success rates over 60 independent trials (error bars show the standard deviation), and solid lines are sigmoid fits. <details> <summary>Image 1 Details</summary>  ### Visual Description ## Line Graph: Success Rate vs. δ = n/d ### Overview The graph compares the success rate of two initialization methods ("Spectral Init." and "Random Init.") across varying values of δ = n/d. Success rate is measured against a threshold of 0.5, with performance trends visualized through error bars and threshold markers. ### Components/Axes - **X-axis**: δ = n/d (ranging from 0 to 10, labeled in increments of 2). - **Y-axis**: Success Rate (Threshold = 0.5), scaled from 0.0 to 1.0 in 0.2 increments. - **Legend**: - Red dashed lines: δ_alg (≈1.5) and δ_nn (≈6.5). - Orange squares: Spectral Init. (solid line with error bars). - Blue circles: Random Init. (solid line with error bars). ### Detailed Analysis 1. **Spectral Init.** (Orange): - Starts at δ=0 with a success rate near 0.0. - Sharp upward transition to 1.0 by δ≈1.5 (marked by red dashed line δ_alg). - Remains flat at 1.0 for δ > 1.5. - Error bars are small and consistent, indicating low variability. 2. **Random Init.** (Blue): - Begins at δ=0 with a success rate near 0.0. - Gradual increase, surpassing Spectral Init. around δ≈6.5 (marked by red dashed line δ_nn). - Plateaus near 1.0 by δ=8. - Error bars grow larger for δ > 6.5, showing increased variability. ### Key Observations - Spectral Init. achieves near-perfect success rapidly (δ≈1.5) but plateaus early. - Random Init. lags initially but outperforms Spectral Init. for δ > 6.5. - Thresholds δ_alg and δ_nn demarcate critical performance shifts for each method. - Error bars suggest Spectral Init. is more reliable at low δ, while Random Init. becomes less consistent at higher δ. ### Interpretation The graph demonstrates a trade-off between initialization methods: - **Spectral Init.** excels in low-δ regimes (δ < 6.5), likely due to structured initialization aligning with problem constraints. - **Random Init.** becomes competitive at higher δ (δ > 6.5), possibly due to stochastic exploration overcoming initialization biases. - The thresholds δ_alg and δ_nn may represent problem-specific breakpoints where algorithmic efficiency or data distribution characteristics shift. - Variability in Random Init. at high δ suggests sensitivity to hyperparameters or noise, whereas Spectral Init.'s consistency implies robustness in controlled settings. </details> On the other hand, there exists a threshold δ alg = δ alg ( h ) ∈ (0 , ∞ ] and a spectral algorithm 2 that achieves weak recovery for any δ > δ alg [MM19], and it is conjectured that no polynomial time algorithm achieves weak recovery for δ < δ alg . In general, the two algorithmic and information theoretic thresholds do not coincide, i.e. δ IT ( h ) ≤ δ alg ( h ) always, and it is quite easy to come up with functions h for which δ IT ( h ) < δ alg ( h ) strictly. Finally (still assuming k = 1 for now), it is possible to construct an empirical risk of the form (3), with f θ ( x ) = ⟨ θ , x ⟩ and a gradient-based algorithm that achieves weak recovery for any δ > δ alg . This construction essentially mimics power iteration for the ideal spectral method. It is also known that for generic multi-index models (with 'generative exponent' at most two), a standard two-layer neural network can achieve weak recovery from n = ˜ O ( d ) samples 3 [DTA + 24, ZWFL25]. These results leaves open two crucial questions: - Q1. Do standard neural networks of the form (2) achieve weak recovery for all δ above the algorithmic threshold δ alg ? If not, what is the gap? How does it depend on architecture (activation function, width, . . . ) and algorithm (loss function, initialization, stepsize, . . . )? - Q2. What is the learning mechanism implemented by gradient descent for the two-layer network? Is it an approximation of a spectral method? If so, which spectral method? As a preview of our results, Fig. 1 compares the success probabilities of two algorithms for the case k = 1 and h ( z, ε ) = z 2 (noiseless phase retrieval), whereby we perform GD with respect to 2 Existence of a spectral algorithm that matches the algorithmic threshold δ alg is an open problem for k > 1. 3 We use ˜ O to hide logarithmic factors. Figure 2: Grokking phenomenon for a GeLU neuron learning the noiseless phase retrieval problem ( d = 5000, δ = n/d = 17 . 5, single run). Here ρ ( t ) := |⟨ θ ( t ) , θ ∗ ⟩| / ∥ θ ( t ) ∥ denotes the cosine similarity between the learned and target parameters. Light-colored markers show data points; solid lines connect them for visual clarity. Initially, training loss decreases while test loss stays near its initial value and misalignment remains high (overfitting phase). Subsequently, the misalignment drops to zero-indicating successful learning of the signal direction-while both losses plateau at around 0 . 6 and the generalization error drops to zero (grokking phase). <details> <summary>Image 2 Details</summary>  ### Visual Description ## Line Graph: Model Risk and Probability Over Training Time ### Overview The graph depicts three metrics over training time (t) from 0 to 100: 1. **1 - ρ(t)** (blue line): Likely represents confidence or error rate reduction. 2. **Risk(θ(t))/Risk(θ(0))** (orange line): Training risk ratio relative to initial risk. 3. **Risk_test(θ(t))/Risk_test(θ(0))** (red line): Test risk ratio relative to initial test risk. ### Components/Axes - **X-axis**: "Training Time t" (0–100, linear scale). - **Y-axis**: Unitless ratio (0.0–1.0, increments of 0.2). - **Legend**: Top-right corner, with three entries: - Blue: `1 - ρ(t)` - Orange: `Risk(θ(t))/Risk(θ(0))` - Red: `Risk_test(θ(t))/Risk_test(θ(0))` ### Detailed Analysis 1. **Blue Line (1 - ρ(t))**: - Starts at 1.0 (t=0), drops sharply to ~0.2 by t=60, then plateaus near 0. - Suggests rapid improvement in a metric (e.g., error rate) during early training. 2. **Orange Line (Training Risk Ratio)**: - Begins at ~0.98 (t=0), declines gradually to ~0.6 by t=80, then stabilizes. - Indicates diminishing returns in training risk reduction after ~40 units of training. 3. **Red Line (Test Risk Ratio)**: - Starts at 1.0 (t=0), dips to ~0.95 by t=20, then rises to ~0.65 by t=80. - Shows initial improvement in test risk, followed by stabilization at a higher value than training risk. ### Key Observations - All lines start near 1.0, reflecting baseline risk/confidence. - **Blue line** declines fastest, implying rapid early gains in the modeled metric. - **Test risk (red)** plateaus higher than training risk (orange), suggesting potential overfitting. - **Training risk (orange)** shows smoother convergence compared to test risk. ### Interpretation - The sharp drop in `1 - ρ(t)` (blue) may represent a confidence threshold or error rate reduction, peaking early in training. - The divergence between training and test risk ratios (orange vs. red) highlights a gap in generalization, with test risk stabilizing at ~0.65 vs. training risk at ~0.6. This could indicate overfitting or dataset shift. - The plateauing behavior of all lines after ~60–80 training steps suggests diminishing returns, where further training yields minimal improvements. ## Critical Notes - **Legend Accuracy**: Colors match lines exactly (blue=blue, orange=orange, red=red). - **Spatial Grounding**: Legend is top-right; axes are labeled with clear increments. - **No Outliers**: All lines follow smooth, expected trends without anomalies. - **Language**: English only; no non-English text present. </details> the risk of Eq. (3), with m = 1, Huber loss ℓ and f θ ( x ) = GeLU ( ⟨ θ , x ⟩ ). The two sets of curves correspond to different initializations: ( i ) The optimal spectral initialization given in [MM19]; ( ii ) Random initialization. (See Section 4 for further details.) Success is declared when GD achieves |⟨ θ , θ ∗ ⟩| / ∥ θ ∥ ≥ 1 / 2. We also report the algorithmic threshold for weak recovery, which is δ alg = 1 / 2 in this case. We observe a factor-5 gap between the optimal initialization and the random initialization, and a factor-12 difference between the weak recovery threshold δ alg and the observed threshold for the one-neuron network. Our prediction for the latter is δ NN ≈ 6 . 0 and matches well the empirical results. The threshold at δ NN also predicts the emergence of grokking as demonstrated in Fig. 2 (same setting of Fig. 1, at δ = 17 . 5). 'Grokking' refers to the widely reported empirical phenomenon according to which generalization error (the gap between test and train error) first increases to reach a plateau and then (often abruptly) decreases, with a corresponding decrease in test error [PBE + 22]. This phenomenology is clarified by our theoretical analysis: 1. Weobserve grokking for δ > δ NN . Indeed, feature learning (which takes place only for δ > δ NN ) drives the abrupt decrease of generalization error in the late phase of training. 2. The training time required for the drop in generalization error increases as δ ↓ δ NN , (making the drop difficult to observe for δ close to δ NN ). This is because feature learning is driven by a spectral phase transition in the Hessian of the empirical risk, and the spectral gap vanishes as δ ↓ δ NN . 3. The maximum of generalization error (the generalization error at the plateau) decreases for δ ≫ δ NN , in agreement with uniform convergence arguments. We will provide a conceptually simple quantitative explanation of the threshold δ NN for neural networks, and its gap with respect to the optimal algorithmic threshold δ alg . Namely, we will study the spectrum of the Hessian ∇ 2 Risk ( Θ ( t )) along the GD trajectory Θ ( t ) and determine a threshold δ NN ( t ) for the emergence of a descent direction (an eigenvector with negative eigenvalue) aligned with the latent subspace. Our prediction for the feature learning phase transition is δ NN = lim t →∞ δ NN ( t ). Here, the limit t →∞ is taken after n, d →∞ . In other words, we are studying the Hessian along the GD trajectory on time horizons of order O (1). This picture provides a conceptual explanation of the gap between optimal algorithmic threshold δ alg and neural network δ NN , since the latter corresponds to a spectral phase transition for a 'sub-optimal' preprocessing of the data. We emphasize that, while our determination of the spectral threshold δ NN (i.e. the study of the Hessian along the GD trajectory for t = O (1)) is fully rigorous, its connection to the GD behavior on longer time scales is a hypothesis based on a heuristic argument. Numerical simulations support this hypothesis but its mathematical proof is an open problem. ## 1.1 Main results: An informal overview ## 1.1.1 Background We train the model (2) using full-batch GD with step size η > 0: $$\theta ( t + 1 ) - \theta ( t ) = - n \gamma e Risk ( \theta )$$ and initialization ( θ j ) j ≤ m i.i.d. ∼ Unif ( S d -1 ), where S d -1 denotes the unit sphere in d dimensions (a constant scaling of θ j can be absorbed in the activation function). Throughout, we will focus on the proportional asymptotics $$\frac { n , d } { 5 } \rightarrow \infty , n / d = \delta _ { n } \rightarrow \delta \in (0 , \infty ) .$$ We begin by decomposing the latent space span( Θ ∗ ) into a subspace that (as we will see) can be learned by O (1) gradient evaluations ('easy') and the complementary subspace ('hard'). ̸ For z ∼ N (0 , I k ) independent of ε ∼ N (0 , 1), let y := h ( z, ε ). In the case k = 1, there is an obvious difference between data distributions such that E [ T ( y ) z ] = 0 for some function T : R → R , and others, for which E [ T ( y ) z ] = 0 for all such T . In the former case ̂ θ ( y , X ) = n -1 ∑ n i =1 T ( y i ) x i achieves weak recovery for all δ > 0. In the latter, weak recovery is information theoretically impossible for δ small enough. For k > 1, the two situations can coexist for different subspaces of span( Θ ∗ ). For k ′ ≤ k , let O ( k, k ′ ) ⊆ R k × k ′ be the set of matrices with with orthonormal columns ( O ( k, k ′ ) is also known as Stiefel manifold). For any U ∈ O ( k, k ′ ) ⊆ R k × k ′ , define the projectors P U := UU T and P ⊥ U := I k -P U , and the linear subspace spanned by the columns of U , span( U ). Definition 1.1. Let the collection of hard orthonormal sets be $$u _ { H } := \{ U \in O ( k , k ^ { 2 } ) : E$$ where ξ ∼ N (0 , 1) is independent of ( x, z, ε ) . We then define the hard subspace as the span of the maximal hard orthonormal set $$\frac { 1 } { n } - \frac { 1 } { n + 1 } ( U _ { H } ) ^ { 2 } . U _ { H } \in arg max _ { V \in V _ { H } } \mu _ { H } = span ( U _ { H } ) .$$ With an abuse of terminology, we also refer to the linear subspace { Θ ∗ u : u ∈ U H } ⊆ R d as to the hard subspace. We denote by r = r := max dim( U ) the dimension of the hard subspace. We refer to the orthogonal complement U E := U ⊥ H as to the easy subspace . H U ∈ U H As an example, consider the case k = 2 and y = h ( z 1 , z 2 ), for a generic function h that is symmetric under permutations of its arguments h ( z 1 , z 2 ) = h ( z 2 , z 1 ). Then U H = span((+1 , -1)) and U E = span((+1 , +1)). It is immediate to show that the definition of U H is unique, because the union of hard subspaces is a hard subspace. In the terminology of [DLB25], U H is the maximal subspace with generative index greater or equal to two (denoted by k ( U ⊥ H ) in [DLB25]). We emphasize that the term 'hard' here is not to be interpreted in its meaning in computational complexity theory. With an abuse of terminology, we also refer to any orthonormal basis U ∗ H of U H as to the hard subspace: there will be no risk of confusion because all statements will be invariant under change of basis. Similarly, we denote by U ∗ E an (arbitrary) orthonormal basis of U E . In alternative to Eq. (6), we use the equivalent definition $$y if E \left[ P _ { U } z \right ] = 0 .$$ In words, along a hard subspace, the expectation of P U z carries no information even after conditioning on both the target and the orthogonal complement. As shown in [DTA + 24], the network (weakly) learns the easy subspace U E within O (1) iterations of gradient descent (for generic a i 's) as long as n/d → α > 0. Namely, there exists ε 0 > 0 and t = O (1) such that, for m enough, and under certain conditions on the m , the step size η , and the activation function σ : $$\sum _ { n = 1 } ^ { \infty } P ( s _ { \min } ( t ) \in U _ { E } ) \geq e _ { 0 }$$ where s min ( · ) denotes the minimum singular value. On the other hand, as stated in the next section, after any constant (in n, d ) number of gradient steps, the weights Θ ( t ) remain roughly orthogonal to the subspace of hard directions. A related result was proved in [DTA + 24] for the case in which all latent directions are hard (i.e. U H = R k ). Our main results are summarized in Section 1.1.3. We characterize a spectral phase transition in the Hessian of the empirical risk, along the GD trajectory. Namely, we analyze the Hessian ∇ 2 Risk ( Θ ( t )) after t = O (1) iterations, and show that, for δ above a certain threshold δ NN ( t ), it develops negative outlier eigenvalues, with the corresponding eigenvector aligned with the hard subspace. This phase transition drives feature learning of the hard directions. ## 1.1.2 Hard directions are not learned within O (1) steps Our first result establishes that hard directions are not learned within O (1) steps of gradient descent. In particular for any fixed t , Θ ( t ) is asymptotically orthogonal to U H . We use an exact characterization of gradient descent in the limit of n, d →∞ , for any constant number of iterations t , which goes under the name of dynamical mean field theory 4 (DMFT) [SZ81, SZ82, CMW20, CCM21]. For any fixed t , DMFT construct a sequence of random vectors Θ ∗ ∈ R k , Θ(0) , . . . , Θ( t ) ∈ R m such that, for any suitably regular test function ϕ : $$\sum _ { i = 1 } ^ { d } \limits _ { n , d \rightarrow \infty } \sum _ { j = 1 } ^ { d } \limits _ { i = 1 } ^ { d } \phi ( \sqrt { d \theta ( t ) } )$$ 4 As shown in [CMW20] this is in fact a consequence of the state evolution analysis of approximate message passing algorithms. where Θ ( t ) i is the i -th row of Θ ( t ) and same for Θ ∗ i . (Throughout the paper, we use boldface to denote pre-limit vectors or matrices and normal font to denote their asymptotic counterpart.) Using this characterization, we prove that, for any bounded continuous function T : R m ( t +1) → R (see Lemma 2.2): $$\rho _ { n , d } \rightarrow \infty$$ In words, any estimator ̂ Θ computed in constant number of gradient steps and possibly postprocessing remains uncorrelated with the hard directions Θ ∗ P U ∗ H . We note that the DMFT construction of the random vectors (Θ(0) , . . . , Θ( t ) , Θ ∗ ) (as well as the random vectors ( V (0) , . . . , V ( t ) , V ∗ , ε ) which characterize the limit of ( Θ (0) T x i , . . . , Θ ( t ) T x i , Θ T ∗ x i , ε i )) is explicit and can be implemented numerically, see Section 2.2. This allows us to study the Hessian at the current position Θ ( t ), which is our main goal, as we discuss next. ## 1.1.3 Hessian phase transitions after gradient descent Figure 3: Cartoon of feature learning via GD on the empirical risk landscape. Within O (1) steps, the network learns easy directions (blue arrow) to a saddle point. Escaping this point requires moving along hard directions (green arrow), governed by the Hessian's negative eigenvalues. <details> <summary>Image 3 Details</summary>  ### Visual Description ## 3D Surface Plot: Optimization Landscape with Critical Points ### Overview The image depicts a 3D surface plot representing an optimization landscape. A red dot labeled "∇ = 0" marks a saddle point at the center of the surface. Two directional arrows emanate from this point: a blue arrow labeled "easy" pointing toward a blue dot labeled "Θ(0)" (likely a local minimum), and a green arrow labeled "hard" pointing downward into a steeper region of the surface. The surface itself has a grid pattern to indicate curvature and gradient direction. --- ### Components/Axes - **Surface**: A 3D terrain-like plot with a grid overlay, representing the function's topography. - **Red Dot (Center)**: Labeled "∇ = 0", indicating a critical point (saddle point) where the gradient vanishes. - **Blue Arrow ("easy")**: Points from the saddle point to a blue dot labeled "Θ(0)", suggesting a direction of easier descent toward a local minimum. - **Green Arrow ("hard")**: Points downward from the saddle point into a steeper region of the surface, labeled "hard" to denote a more challenging descent path. - **Grid Pattern**: Implicitly encodes curvature and gradient magnitude across the surface. --- ### Detailed Analysis 1. **Saddle Point (∇ = 0)**: - Located at the center of the surface. - Represents a critical point where the gradient is zero but is neither a local minimum nor maximum. - The surface curves upward in one direction (toward the blue dot) and downward in another (toward the green arrow's path). 2. **Blue Arrow ("easy")**: - Color: Blue (matches legend label "easy"). - Direction: Points from the saddle point to the blue dot labeled "Θ(0)". - Implies a path of least resistance toward a local minimum. 3. **Green Arrow ("hard")**: - Color: Green (matches legend label "hard"). - Direction: Points downward from the saddle point into a steeper region. - Indicates a direction of higher gradient magnitude, requiring more computational effort to navigate. 4. **Surface Curvature**: - The grid pattern shows opposing curvatures around the saddle point: - Concave upward toward the blue dot (local minimum). - Concave downward toward the green arrow's path. --- ### Key Observations - The saddle point ("∇ = 0") is a critical point with zero gradient but non-extremal behavior. - The "easy" direction (blue) aligns with a local minimum ("Θ(0)"), suggesting a favorable descent path. - The "hard" direction (green) leads to a steeper region, representing a more challenging optimization trajectory. - The surface's grid confirms the saddle point's mixed curvature properties. --- ### Interpretation This diagram illustrates the challenges of optimization in non-convex landscapes: 1. **Saddle Points as Bottlenecks**: The red dot ("∇ = 0") represents a critical point where gradient-based methods stall, as the gradient vanishes but the point is not a local extremum. 2. **Directional Bias**: The "easy" (blue) and "hard" (green) labels highlight the asymmetry in optimization paths. Algorithms may favor the "easy" direction toward "Θ(0)" but risk getting trapped in suboptimal regions if the "hard" path is not explored. 3. **Curvature Dynamics**: The opposing curvatures around the saddle point emphasize the need for second-order methods (e.g., Hessian-aware optimization) to escape such regions effectively. The visualization underscores the importance of understanding landscape geometry in optimization, particularly for non-convex problems where saddle points and mixed curvature dominate. </details> In order to understand the dynamics beyond time t = O (1), we need to consider deviation from the DMFT trajectory. These deviations are of order 5 1 / √ d , but their dynamics depends on the structure of the Hessian at Θ ( t ). It will be convenient to study the rescaled Hessian $$H ( t ) = m \sqrt { 2 } Risk ( \theta ( t ) ) .$$ Our main result is a sharp characterization asymptotic of the smallest eigenvalues of H ( t ) and associated eigenvectors under the proportional asymptotics n, d →∞ , n/d → δ , in two limit cases m = 1 and m ≫ 1. Surprisingly, these cases can be treated in unified fashion from a mathematical perspective. Case m = 1 . In this case H ( t ) ∈ R d × d takes the form $$H ( t ) := - 1 _ { n } \sum _ { i = 1 } ^ { n } g ( y _ { i } , \theta )$$ 5 This can be proven using the finite sample analysis of AMP of [RV18] and the mapping between AMP and generalized first order methods of [CMW20]. $$g ( y , z ) := a ^ { 2 } l ( y , a o ( z )$$ Case m ≫ 1 . In this case H ( t ) ∈ R md × md has a block structure ( m 2 blocks of dimensions d × d ). Lemma 2.3 below shows that, for a convex loss ℓ ( y, · ), the smallest eigenvalues and corresponding eigenvectors of m ∇ 2 Risk ( Θ ) are well-approximated by a block diagonal matrix H diag ( t ), up to an O (1 /m ) correction. (Here we are considering large width m after the limit n, d →∞ ). The j -th block of H diag ( t ) is a j H j ( t ) where $$H _ { j } ( t ) := - \sum _ { i = 1 } ^ { n } g ( y _ { j } , \theta ( t ) , x _ { j } ; j )$$ $$, f e ( x ) ) ^ { \prime } ( \theta _ { j } x + b _ { j } ) .$$ Hence, characterizing the low-lying eigenvalues of H ( t ) amounts to analyzing a random matrix with similar structure in the two cases m = 1 and m ≫ 1. From the mathematical viewpoint, the challenge lies in the fact that the parameters Θ ( t ) depend in an intricate way on the random data { x i , y i } n i =1 through gradient descent updates. Our main results can be summarized as follows (we refer for concreteness to the case m ≫ 1, and the block H j ( t ), and assume throughout the proportional asymptotics (5)). - Asymptotics of the empirical spectral distribution (ESD). For any t ≥ 0, let µ j n,d ( t ) denote the empirical spectral distribution of block j . We then determine a distribution µ j ∞ ( t ) such that $$u _ { n , d } ( t ) \rightarrow u ^ { i } _ { 0 } ( t ) .$$ (We use ⇒ to denote weak convergence and the convergence is understood to hold in probability). Under certain regularity assumptions on the loss and activation functions, the support of µ j ∞ ( t ) is bounded below, and we denote by c j ( t ) its infimum. Thresholds for emergence of informative eigenvectors. We characterize thresholds δ ∗ j for each block j such that, for δ > δ ∗ j , the following happens: - There exists t 0 < ∞ such that, for all t ≥ t 0 , the spectrum of H j ( t ) contains a negative eigenvalue z n,d j ( t ) below the infimum of the support of µ j ∞ . Namely, for some fixed ∆ > 0 $$\frac { p - \lim _ { n \rightarrow \infty } z _ { j } ^ { n } ( t ) } { n , d \rightarrow \infty } = z _ { j } ^ { n } ( t ) \leq \min ( c , ( 1 , 0 ) - \Delta ( 1 8 )$$ - The corresponding eigenvector ξ n,d j ( t ) has a projection on the hard subspace with magnitude bounded away from zero. No informative eigenvector below spectral thresholds. Viceversa, for δ < δ ∗ j and for all t large enough, no eigenvector of H j ( t ) has projection onto the hard subspace which is asymptotically bounded away from zero. This leads us to define the overall feature learning threshold as δ NN := min j ∈ [ m ] δ ∗ j . In summary, when the sample size is large enough, the Hessian has a negative isolated eigenvalue whose corresponding eigenvector is correlated with the hard directions. When the sample size is small, the Hessian has no eigenvector correlated with the hard directions. Our results provide explicit formulas for these thresholds, enabling sharp and numerically computable predictions for this phenomenon. We note that the threshold δ ∗ j is dependent on the details of the learning procedure, including the loss function, learning rate, initialization, and activation function. Our results allow to study the dependence of feature learning with respect to these factors that can be designed to minimize the sample size for learning. Outline of the paper. The remainder of the paper is organized as follows. Section 1.2 reviews related literature, and Section 1.3 summarizes our notation. Section 2 presents our main results and Section 3 discusses some of their consequences. Section 4 illustrates our theoretical predictions through numerical experiments. Section 5 provides an outline of the proofs. Section 6 concludes with a discussion of implications and future directions. Some technical preliminaries are provided in Appendix A. Detailed proofs are provided in Appendices B and C, supporting technical lemmas are collected in Appendix D, and additional numerical experiments are presented in Appendix E. ## 1.2 Related literature Optimal feature learning. The literature on multi-index models is quite extensive. Here we limit ourselves to results relevant to the specific setting of our paper, referring to [BH25] for a recent review. In particular we assume k fixed, isotropic covariates and focus on the weak recovery problem, i.e., estimating the latent directions with non-vanishing correlation. We generically refer to weak recovery of the latent space as 'feature learning'. Information-theoretically, the sample size required to learn a k -index model is linear in the dimension d [BKM + 19, AMB + 19], as expected from counting degrees of freedom. A significant amount of work has been devoted to understanding whether this can be achieved by polynomialtime algorithms [BKM + 19, MM19, LL20, CM20, DPVLB24, DLB25, TDD + 25, KZM25, DDM + 25]. We will refer to this as 'efficient feature-learning'. The situation is particularly simple in the case k = 1. If E [ z | y ] is not identically zero (in the notation of Section 1.1.1), then weak recovery of θ ∗ is possible for any δ > 0. If on the other hand E [ z | y ] = 0 (this is the case if y only depends on z via | z | ), [MM19] showed that efficient feature learning is possible when δ crosses a certain threshold δ alg by using an optimal spectral method. It has been conjectured that no polynomial-time algorithm succeeds for δ < δ alg [MM19, BKM + 19]. While it is possible to construct examples (functions h ) for which δ alg = ∞ , this situation is nongeneric in the sense that it requires the function h to satisfy an additional equality constraints. For k > 1, the simplest case is the one in which all directions are easy, in our language U E = R k . Starting with [ABAM23], several results established that efficient feature learning is possible from n = O ( d ) samples for certain classes of distributions of this type. On the other hand, for several interesting problems U E ⊊ R k , e.g. when h ( z, ε ) is invariant under permutations of some of the entries of z . Recently, [DLB25] proved that sequential spectral method achieve efficient feature learning from linearly many samples even when hard directions are present, generically. More generally, they related the scaling of sample size to a property of the link function h , known as the 'generative exponent'. The generic situation (covering common examples) corresponds to generative exponent at most two. The results of [DLB25] do not yield a sharp threshold δ alg for polynomial-time feature learning. For certain classes of functions h ( z, ε ) [TDD + 25] put forward a conjectured value for the threshold δ alg , based on the analysis of an AMP algorithm 'with side information.' A proof of this conjecture remains an outstanding problem. Optimal spectral methods are developed in [KZM25, DDM + 25]. These works establish that weak recovery is feasible for δ > δ alg but do not provide evidence for hardness for δ < δ alg . Feature learning in neural networks. Empirically, gradient-based training of neural networks efficiently finds hidden representations aligned with low-dimensional structure in the data. A large number of theoretical works studies how and when neural networks can efficiently perform feature learning in single-index, multi-index, and other models with planted structure [AGJ21, DLS22, ABAM22, ABAM23, BES + 22, BBPV23, WNL24, BMZ24, LOSW24, DTA + 24, FWNL25, ADK + 24, MU25, ZWFL25, FW26]. For single-index models, [LOSW24] showed that neural networks can efficiently learn singleindex polynomial targets with ˜ O ( d ) samples through sample reuse, and [MU25] provided a precise picture of the dynamics of learning single-index models (with information exponent one, within O (1) time) using wide networks under proportional scaling, demonstrating that overfitting occurs at a slower timescale than feature learning. Early results [AGJ21, DLS22, ABAM22, ABAM23] established that neural networks can achieve feature learning in multi-index models, albeit with suboptimal sample complexity. For multi-index models, [DTA + 24] proved that n = O ( d ) samples are sufficient for efficient feature learning using neural networks and batch gradient descent when the generative exponent is one. In other words, gradient descent can match the ideal sample size scaling as predicted by the generative exponent. In concurrent work, [ZWFL25] showed that neural networks trained via layer-wise gradient descent can efficiently learn generative-exponent-two, non-generative-staircase multi-index models from ˜ O ( d ) samples. Despite such a substantial amount of work, sharp thresholds for feature learning in neural networks remain unknown. In addition, the precise interplay between activation functions, data structure, loss functions, and initialization in determining the feature learning behavior is not well understood. Spiked Models. Spiked random matrix models provide a mathematical framework for understanding when and how informative directions can be detected from noisy data. In a standard spiked model, the spectrum of the bulk follows the generalized Marchenko-Pastur law [MP67, BY93, BS10]. When low-rank signals are introduced, the seminal work of [BAP05] identified the Baik-Ben Arous-P´ ech´ e (BBP) transition, which precisely characterizes the critical signal-to-noise ratio at which outlier eigenvalues emerge from the bulk spectrum. Above this threshold, the top eigenvectors become correlated with the latent signal, enabling spectral estimators to recover the hidden directions. Subsequent work has refined this picture through local laws and improved deformation analyses [BGN11, KY13, KY17]. Spiked models naturally emerge when designing and analyzing spectral methods used as initialization strategies for nonconvex problems such as phase retrieval [CC17, MM19]. The structure of the Hessian under generalized linear models also takes the same form and was further studied (for m = k = 1) in [LM21]. However, this type of analysis only provides information on the Hessian at a point θ (0) independent of the data itself. A key technical challenge that we face is the one of analyzing the Hessian at a point Θ ( t ) that depend on the data X , y . The recent paper [AGHJ25] analyzes the Hessian after a certain number of stochastic gradient descent (SGD) steps. However they assume that the Hessian is computed on fresh data, independent of the data used for SGD, hence sidestepping the challenge of dependencies between Θ ( t ) and the data. Considering the setting in which the same data is queried multiple times (and in particular Hessian is computed on the same data used for GD) poses technical challenges but is crucial to capture several crucial learning phenomena 6 . 6 For instance, the generalization error vanishes if data is not re-used, and the required sample size can be significantly larger than with multiple-pass GD [DTA + 24]. ## 1.3 Summary of notations We use n, d, m, k to denote the sample size, input dimension, network width, and number of target features, respectively. We use c, C, c i , C i , . . . to denote constants that may change from line to line; these constants are independent of n, d, m but may depend on k and on the constants that appear in the regularity assumptions on h, σ, ℓ . We use boldface to denote pre-limit vectors and matrices, and normal font to denote their asymptotic limiting counterparts. We use ∥ · ∥ to denote the Euclidean norm when applied to vectors and the operator norm when applied to matrices. We denote by I p the p × p identity matrix. A function ϕ : R p → R is called pseudo-Lipschitz of order s (for s ≥ 1), denoted ϕ ∈ PL s , if there exists a constant L > 0 such that for all x, y ∈ R p , we have $$\left | \phi ( x ) - \phi ( y ) \right | \leq L ( 1$$ We use ≲ and ≳ to hide constants as above. We write o P (1) for a quantity that converges to zero in probability as n, d →∞ . A sequence X n = O P (1) means that { X n } is tight: for every ε > 0, there exists M < ∞ such that sup n P ( | X n | > M ) ≤ ε . We denote by W p ( µ, ν ) the Wasserstein p -distance; when p = 2, we drop the subscript. We use λ min ( · ) to denote the minimal eigenvalue of a real symmetric matrix and S d -1 to denote the unit sphere in R d . We use ⇒ to denote weak convergence and p - - → or p-lim n,d →∞ to indicate convergence in probability. We use H := { z ∈ C : Im( z ) > 0 } to denote the upper half plane. Given an infinite symmetric positive semidefinite matrix S = ( S ( t, s ) : t, s ≥ 0) and a sequence a = ( a ( t ) : t ∈ N ), we write Z ∼ GP ( a, S ) if Z is a Gaussian process with mean a and covariance S . ## 2 Results ## 2.1 Problem setup and assumptions This subsection introduces the notation for the gradient descent (GD) updates and states the regularity assumptions that will be used in our analysis. Throughout, we use ℓ ′ and ℓ ′′ to denote the derivatives of ℓ with respect to its second argument. To express the update of GD in a compact form, we introduce the following notation. Define F j ( Θ T x i , Θ T ∗ x i , ε i ) := ( η/m ) a j ℓ ′ ( y i , f Θ ( x i )) σ ′ ( θ T j x i + b j ) and X T = [ x 1 . . . x n ] ∈ R d × n . Then, the GD update reads $$\theta ( t + 1 ) - \theta ( t ) = - \frac { 1 } { n }$$ where F j ( · ) ∈ R n denotes the coordinate-wise application of F j ( · ). Wealso define F T := [ F 1 . . . F m ] . We now state our regularity assumptions. Assumption 1. We assume the following. - (a) The data { ( x i , y i ) } n i =1 are i.i.d. with x i ∼ N ( 0 , I d ) and y i = h ( Θ T ∗ x i , ε i ) , where ε i ∼ N (0 , 1) is independent of x i and Θ ∗ ∈ R d × k has orthonormal columns. - (b) h ( · ) is locally Lipschitz continuous. - (c) σ ( · ) ∈ C 2 ( R ) , ∥ σ ′ ∥ ∞ ⩽ C , ∥ σ ′′ ∥ ∞ ⩽ C . - (d) 0 ⩽ a j ⩽ C , | b j | ⩽ C . - (e) For each y ∈ R , ℓ ( y, · ) ∈ C 2 ( R ) . Further, ℓ ( y, y ) = 0 for all y and ∥ ℓ ′′ ∥ ∞ ⩽ C , ∥ ℓ ′ ∥ ∞ ⩽ C . - (f) The hard subspace has dimension at least one. Remark 2.1. By rotation invariance of the distribution of the covariates x i , the learning dynamics is unchanged if we replace Θ ∗ , Θ ( t ) by Q Θ ∗ , Q Θ ( t ), for Q an arbitrary orthogonal matrix. By choosing a uniformly random Q , there is no loss of generality in assuming the following. For each ϕ ∈ PL 2 , we have $$\sum _ { i = 1 } ^ { d } \phi ( \sqrt { d \theta _ { i } } ) \rightarrow E [ \phi ( \theta _ { i } ) ]$$ where Θ ∗ ∼ N (0 , I k ), and Θ ∗ i denotes the i -th row of Θ ∗ . We note that the assumptions ∥ ℓ ′ ∥ ∞ ⩽ C and ∥ ℓ ′′ ∥ ∞ ⩽ C can be relaxed to polynomial growth conditions | ℓ ′ ( y, z ) | ⩽ C (1 + | y | + | z | ) C and | ℓ ′′ ( y, z ) | ⩽ C (1 + | y | + | z | ) C if we additionally assume | y | ⩽ C and ∥ σ ∥ ∞ ⩽ C . This relaxation allows us to include the square loss. ## 2.2 Preliminaries: Discrete time DMFT To characterize the gradient descent dynamics in the high-dimensional limit, we employ the discrete time dynamical mean field theory (DMFT) framework, which has a rich history in physics and mathematics, and has been studied extensively in the literature. DMFT provides a low-dimensional characterization of training dynamics, whose complexity does not scale with n, d . We briefly review the main setup below, which is mostly a restatement of results from [CCM21]. The DMFT asymptotics are determined by a set of self-consistent systems of stochastic processes. We define the DMFT process via $$\sum _ { s = 0 } ^ { t - 1 } R _ { s } ( t , s ) = - \frac { 1 } { 6 } \sum _ { s = 0 } ^ { t - 1 } R _ { s } ( t , * ) + n Q ( t ) ,$$ Here, ( W ∗ , W ( t ); t ∈ N ) and ( Q ( t ); t ∈ N ) are Gaussian processes taking values in R m with covariances determined by kernels C θ and C ℓ : $$\begin{aligned} E \left[ W \cdot W ^ { T } \right ] = I _ { k } , E \left[ W \right ] & = C _ { a } ( t , s ) . \\ & = C _ { a } ( t , s ) . \end{aligned}$$ and $$\sum _ { t \in ( 2 , s ) } [ q ( t ) q ( s ) ] = - \frac { 1 } { g } c ( t , s ) , t , s$$ The kernels C θ , C ℓ , R θ , R ℓ are determined as the unique solution of the following self-consistency equations $$\begin{aligned} C_0( t , s ) = E \left[ e^{(t)e^{(s)}} \right] , \\ C_0( t *, s ) = E \left[ e^{(t)e^{*}} \right] , \\ R_0( t , s ) = \frac{1}{n} \left[ \int_{-\infty}^{\infty} aF(V(t)) dW \right] , \\ C_0( t , s ) = - 1_{T^2} E \left[ F(V(t), W*, e) F( \cdot | \cdot ) \right] , \\ R_0( t , s ) = - 1_{T^2} E \left[ \frac{aQ(s)}{aW(s)} \right] , \\ R_0( t *, s ) = - 1_{T^2} E \left[ \frac{aF(V(t), W*, e) aW(s)}{aW(s)} \right] . \end{aligned}$$ We also denote V ∗ := W ∗ , and initialize the DMFT process with Θ(0) ∼ N (0 , I m ) independent of Θ ∗ . When computing partial derivatives with respect to W ( s ), W ∗ above, it is understood that V ( t ) is a function of W ∗ , ( W ( s ) : s ≤ t ) via Eq. (22), and similarly for the partial derivative of Θ( t ) with respect to Q ( s ). It is also useful to recall at this point that F ( V, W, ε ) is of order η and the above quantities have well defined limits as η → 0 + and ηt = Θ(1), which is described in [CCM21]. For any fixed time horizon T , the coupled system (22)-(25) can be solved forward in time: the kernels { C θ , R θ , C ℓ , R ℓ } determine the Gaussian processes ( W,Q ), which in turn determine ( V ( t ) , Θ( t )) t ≤ T via (22); plugging these processes back into (25) yields the same kernels but for the next time step. Hence the DMFT process and kernels are uniquely determined by Eqs. (22)-(25) . ## 2.3 Result I: Hard directions cannot be learned in O (1) time This subsection analyzes the gradient descent dynamics over a finite number of steps. We establish that hard directions cannot be learned in O (1) time. The results are summarized in the following two lemmas (the first parts of these lemmas are restatements of a result from [CCM21]). Lemma 2.1. Let ϕ be a test function which is locally Lipschitz continuous and has at most quadratic growth. Under proportional scaling n/d := δ n → δ ∈ (0 , + ∞ ) and Assumption 1, we have $$\sum _ { i = 1 } ^ { n } \phi ( e ^ { i T } x _ { i } , \ldots , \phi ( t ) x _ { i } , \ldots , \phi ( V ( 0 ) , \ldots , \phi ( t ) x _ { i } )$$ where ( V (0) , . . . , V ( t ) , V ∗ , ε ) is the discrete-time DMFT process. Furthermore, the joint distribution satisfies $$, e , h ( V , \epsilon ), P U ^ { + } H ) P U ^ { - } H ] = 0$$ for every measurable ψ ( · ) such that the expectation makes sense. Lemma 2.2. Let ϕ be a test function which is locally Lipschitz continuous and has at most quadratic growth. Under proportional scaling n/d := δ n → δ ∈ (0 , + ∞ ) and Assumption 1, we have $$\sum _ { i = 1 } ^ { d } \phi ( \sqrt { d \theta ( t ) _ { i } } , \ldots , \sqrt { d \theta ( t ) _ { i + 1 } } )$$ where (Θ(0) , . . . , Θ( t ) , Θ ∗ ) is the discrete-time DMFT process. Furthermore, we have $$\left[ \begin{matrix} T ( \theta ( 0 ) , \ldots , \theta ( t ) ) P _ { u } ^ { H } \end{matrix} \right] = 0$$ for each t and for every measurable T : R m ( t +1) → R such that the expectation makes sense. Collectively, these lemmas establish that hard directions cannot be learned in constant time. Specifically, for each fixed time t , the iterates remain orthogonal to the hard directions. The proofs of Lemma 2.1 and Lemma 2.2 are deferred to Appendix B. Intuitively, the proof of the second part of Lemma 2.1 is based on the fact that we can represent V ( t ) as a deterministic function of h ( V ∗ , ε ) and another random vector independent of P U ∗ H V ∗ . We point out that related statements were proven in [DTA + 24], but only for the case in which all latent directions are hard (i.e. U H = R k ). We also note that the joint distributions of ( V (0) , . . . , V ( t ) , V ∗ , ε ) and (Θ(0) , . . . , Θ( t ) , Θ ∗ ) can be computed explicitly via the discrete time DMFT formulas for each finite t (see also in Appendix B). Importantly, the complexity of these formulas scales only with m and t , rather than with n and d . These evolution formulas also provide a precise description of how the easy directions are learned in constant time. Finally, we emphasize an important point. Even if all directions are hard, i.e., U H = R k , and hence GD does not develop positive correlation with the latent space Θ ∗ for any t = O (1), the GD dynamics remains highly non-trivial on the same time horizon. Indeed, under the proportional asymptotics n = O ( d ), the parameters Θ ( t ) evolve to overfit the training data to a certain extent. ## 2.4 Result II: Bulk of the Hessian This subsection characterizes the bulk spectrum of the Hessian, along the GD trajectory. As explained in Section 1.1, the motivation of this analysis is that we expect the learning of the hard subspace to be driven by a spectral phase transition in the Hessian after easy directions have been learned. We first recall the form of the Hessian of our empirical loss for parameter Θ ( t ). Under Assumption 1, we have $$\sum _ { i = 1 } ^ { n } ( c _ { i } ) ^ { 2 } f e ( x _ { i } )$$ where $$\sum _ { i = 1 } ^ { n } f e ( x _ { i } ) = - \frac { 1 } { m } \left[ a _ { 1 } o ^ { r } ( a _ { m } o ^ { r } ) \right]$$ $$\sum _ { i = 1 } ^ { n } f e ( x _ { i } ) = \frac { 1 } { m } \left[ a _ { 1 } o ^ { n } ( \theta T _ { i } x _ { i } + b _ { 1 } ) x _ { i } T _ { i } \cdots \right]$$ for each Θ . We then denote by H diag ( Θ ) := ( m/n ) ∑ n i =1 ℓ ′ ( y i , f Θ ( x i )) ∇ 2 f Θ ( x i ) the rescaled block diagonal part of the Hessian. Lemma 2.3. Under Assumption 1 and additionally ℓ ′′ ( y, z ) ≥ 0 for all y, z , there exists a constant C independent of m,d,n such that, with high probability, $$\lambda _ { \min } ( H _ { diag } ( \Theta ) ) \leq \lambda _ { m }$$ Furthermore, let ξ with ∥ ξ ∥ = 1 be a unit eigenvector associated with the minimum eigenvalue, i.e., ∇ 2 Risk ( Θ ) ξ = λ min ( ∇ 2 Risk ( Θ ) ) ξ . Then there exists a constant C independent of m,d,n such that, with high probability, $$\frac { ( t , H _ { diag } ( \Theta ) \xi ) ^ { C } } { m } \leq \min ( H _ { diag } ( \Theta ) )$$ This establishes that the block diagonal term H diag ( Θ ) plays a dominant role in determining the minimum eigenvalue. This motivates us to analyze the spectrum of the following matrix: $$H _ { j } ( t ) = - \sum _ { i = 1 } ^ { n } g ( e _ { x ; } y ; f _ { e _ { x ; } y ; } )$$ We note that a j H j ( t ) is precisely the j -th block of m H diag ( Θ ( t )). Lemma 2.4. Let µ j n,d ( t ) denote the empirical spectral distribution (ESD) measure of H j ( t ) . Under Assumption 1 and the proportional scaling n/d := δ n → δ ∈ (0 , + ∞ ) , we have µ j n,d ( t ) ⇒ µ j ∞ ( t ) in probability. The limiting measure µ j ∞ ( t ) is characterized by its Stieltjes transform α j t ( z ) for z ∈ H , which is the unique solution α j t ∈ H of $$z + \frac { 1 } { a ^ { d } t } = s ; E [ \frac { G _ { i } ^ { c } } { s + C _ { i } a ^ { d } t } ]$$ where G j t := g ( V ( t ) , h ( V ∗ , ε ); j ) . This lemma characterizes the bulk of the spectrum: µ j ∞ ( t ) follows a generalized MarchenkoPastur law. In order to provide some intuition of this result, note that that H j = X T D j X , for a diagonal matrix D j (with ( D j ) ii := g ( Θ ( t ) T x i , y i ; j ) /n ). If D j were independent of X , then Lemma 2.4 would follow from textbook theory [BS10]. The content of Lemma 2.4 is that the bulk of the spectrum remains the same despite the dependencies between Θ ( t ) T x i , y i and X (which arise in particular because Θ ( t ) depends on X ). The proof, deferred to Appendix B, relies on Gaussian conditioning to show that these dependencies can be accounted (for any t = o ( d )) by a low rank deformation of X . As a consequence (by Cauchy interlacing theorem), the spectrum of the Hessian is only modified by the emergence of outliers, which we analyze next. ## 2.5 Main results: Outliers of the Hessian and phase transitions This subsection characterizes the isolated eigenvalues (outliers) that detach from the bulk, as well as the alignment of the associated eigenvectors with the hard subspace. We establish a phase transition at a computable threshold δ ∗ j ( t ). Recall that µ j ∞ ( t ) is the limiting empirical spectral distribution of H j ( t ) with Stieltjes transform α j t ( z ), and let c j ( t ) be its left edge: $$( 3 7 )$$ This can also be computed via the following edge equation which follows e.g. from [SC95, Section 4] (here we use the convention 1 / 0 = + ∞ ): $$\frac { G _ { i } ^ { j } } { \delta } = \sup \{ - 1 + o \cdot E _ { i } \} , A _ { i } := \frac { \delta + G _ { i } a } { ( - inf _ { j } ) } .$$ Define Θ ∗ H := Θ ∗ U ∗ H ∈ R d × r which is a basis of the hard subspace. Weimpose one more regularity assumption. To state it, we denote O T t := [ Θ T ∗ Θ(0) T . . . Θ( t ) T ] and ˜ F T t := [ F ( V (0) , V ∗ , ε ) T . . . F ( V ( t ) , V ∗ , ε ) T ] , where we recall F T = [ F 1 , . . . , F m ] . Assumption 2. For each t , the matrices E [ O t O T t ] and E [ ˜ F t ˜ F T t ] are positive definite. This assumption requires that no two features are perfectly parallel during training, and no two gradients are perfectly parallel as well. This assumption holds generically and may be removed via a more refined analysis. In the next subsections, we will present our results about outlier eigenvalues/eigenvectors, first for a fixed time t , and then in the large-time limit as t →∞ . (The limit t →∞ is taken after n, d →∞ .) ## 2.5.1 Results at fixed time t For real z < min(0 , c j ( t )), consider the equation (recall that r H = dim( U H ), and U H ∈ R k × r H ) $$\det ( - z I _ { r H } + E \left[ \frac { s } { G _ { t } ^ { 2 } d _ { t } ( z ) } U _ { H } T V _ { * } ( u _ { H } ^ { T } V _ { * } ) ^ { T } \right] = 0$$ where we recall G j t = g ( V ( t ) , h ( V ∗ , ε ); j ). Define the threshold (as usual, δ ∗ j ( t ) = ∞ if the condition is never verified) $$f ^ { \prime } ( t ) := i n f \{ \delta \in ( 0 , \infty ) \}$$ Theorem 1. Under proportional scaling n/d := δ n → δ ∈ (0 , + ∞ ) , Assumption 1, and Assumption 2, we have the following: - (a) When δ > δ ∗ ( t ) , equation (39) has at least one real solution z ∗ < min(0 , c j ( t )) In this case, denote the distinct real solutions as z ∗ 1 , . . . , z ∗ p < min(0 , c j ( t )) with multiplicities k 1 , . . . , k p . Then, for each i = 1 , . . . , p , there exist eigenvalues z ( i ) n, 1 , . . . , z ( i ) n,k i ∈ R and (orthonormal) eigenvectors ξ ( i ) , . . . , ξ ( i ) such that - j . n, 1 n,k i $$\sum _ { q = 1 , \ldots , k _ { i } } | z _ { n , q } - z _ { i } | ^ { ( i ) } P _ { i } \rightarrow 0$$ Furthermore, for each i = 1 , . . . , p , let S i t ∈ R r × k i be a full-rank matrix whose columns are orthonormal and form a basis of $$\sum _ { i = 1 } ^ { n } \left( - z _ { i } ^ { * } I _ { r H } + E \right) \delta G _ { i } ^ { T }$$ We have the following convergence $$\sum _ { q = 1 } ^ { p - \lim _ { n , d \rightarrow \infty } | \theta _ { T } h _ { S _ { n } , q } | ^ { 2 } } | \theta _ { T } h _ { S _ { n } , q } | ^ { 2 }$$ $$\sigma ( G ^ { i } ) ^ { 2 } \alpha ^ { i } ( z _ { t } ) , \ldots T V , ( U ^ { i } T V , ) ^ { - E } \left[ \frac { \sigma ( G ^ { i } ) ^ { 2 } \alpha ^ { i } ( z _ { t } ) } { ( s + G ^ { i } \alpha ^ { i } ( z _ { t } )) ^ { 2 } } U ^ { i } T V , ( U ^ { i } T V , ) ^ { - E } \right]$$ Let ξ p be the unit norm eigenvector corresponding to the p -th smallest eigenvalue of H j ( t ) . Then, for every fixed p , p-lim Θ T ξ p = 0 . - (b) When δ < δ ∗ j ( t ) , equation (39) has no negative solution for z < c j ( t ) . n,d →∞ ∗ H ## 2.5.2 Results for large time t The next assumption states that the discrete-time DMFT process converges to a stationary point. Assumption 3. The sequence of random vectors { V ( t ) } t ≥ 0 defined by the DMFT process has a limit in L 2 , namely lim E [ V ( t ) V 2 ] = 0 for some random vector V . Equivalently, { V ( t ) } t ≥ 0 is Cauchy in L 2 , i.e. lim t,s →∞ E [ ∥ V ( t ) -V ( s ) ∥ 2 ] = 0 t →∞ ∥ -∞ ∥ ∞ . Note that the second form of this assumption can be checked using the DMFT equations, which give access to the quantity E [ ∥ V ( t ) -V ( s ) ∥ 2 ] . Define G j ∞ := g ( V ∞ , h ( V ∗ , ε ); j ) and the limiting measure µ j ∞ ( ∞ ) via its Stieltjes transform α j ∞ ( z ), z ∈ H as the unique solution α j ∞ ∈ H (uniqueness from Theorem 4.3 in [BS10]): $$z + \frac { 1 } { a ^ { d } _ { \infty } } = \delta ; E [ \frac { G ^ { z } _ { \infty } } { b ^ { s } + G ^ { z } _ { \infty } d ^ { s } ( z ) } ]$$ We define c j ( ∞ ) to be the left edge of the limiting measure µ j ∞ ( ∞ ): $$c _ { j } ( \infty ) := i n f s u p p r ^ { j } _ { \infty } ( \infty ) .$$ This can be computed via an analogous of Eq. (38). The following lemma establishes the convergence of the spectral measures. Lemma 2.5. Under Assumption 1 and 3, we have µ j ∞ ( t ) ⇒ µ j ∞ ( ∞ ) as t →∞ and c j ( t ) → c j ( ∞ ) . For z < min(0 , c j ( ∞ )), consider the following equation: $$\frac { \delta G _ { o } } { G _ { o o } d _ { o o } ( z ) } U ^ { H V } ( u ^ { H V } ) ^ { T } = 0$$ Define the threshold (as usual, δ ∗ j ( ∞ ) = ∞ if the condition is never verified) $$f ^ { \prime } ( \infty ) = i n f \{ f \in ( 0 ,$$ Theorem 2. Under proportional scaling n/d := δ n → δ ∈ (0 , + ∞ ) and Assumption 1, 2, and 3, we have the following: - (a) When δ > δ ∗ j ( ∞ ) , equation (47) has at least one real solution z ∗ < c j ( ∞ ) . In this case, denote the distinct real solutions as z ∗ 1 , . . . , z ∗ p < c j ( ∞ ) with multiplicities k 1 , . . . , k p . Then, for each t and each i = 1 , . . . , p , there exist eigenvalues z ( i ) n, 1 , . . . , z ( i ) n,k i ∈ R and (orthonormal) eigenvectors ξ ( i ) n, 1 , . . . , ξ ( i ) n,k i such that $$\sum _ { i = 0 } ^ { \infty } ( t ) g _ { n , q } = z _ { n , q } g _ { n , q }$$ Furthermore, for each i = 1 , . . . , p , let S i ∞ ∈ R r × k i be a full-rank matrix whose columns are orthonormal and form a basis of $$\sum _ { i = 1 } ^ { n } \left( - z _ { i } I _ { r _ { H } } + E \right) \left[ \frac { s G _ { o _ { i } } } { G _ { o _ { o } } a _ { o _ { i } } ( z _ { i } ) } U ^ { T _ { H } } V _ { ( u * T _ { V } ) } + G _ { o _ { o } } a _ { o _ { o } } ( z _ { i } ) \right]$$ We have the following convergence $$\sum _ { n = 1 } ^ { \infty } p _ { n } q _ { n } = \sum _ { q = 1 } ^ { \infty } | e ^ { i q } T _ { H q } |$$ $$E \left[ \frac { \delta ( G_{\infty } ) ^ { 2 } \alpha ^ { 2 } \delta ( z _ { i } ) } { ( s + G_{\infty } \alpha \infty ( z _ { i } ) ) ^ { 2 } } U ^ { * } H ^ { * } V _ { \ast } ( U ^ { * } H ^ { * } V _ { \ast } ) \right]$$ - (b) When δ < δ ∗ j ( ∞ ) , equation (47) has no negative solution for z < c j ( ∞ ) . Let ξ p be the unit norm eigenvector corresponding to the p -th smallest eigenvalue of H j ( t ) . Then, there exists a finite t ∗ j such that, for every t > t ∗ j , we have p-lim n,d →∞ Θ T ∗ H ξ p = 0 for every fixed p . We defer all proofs to Appendix C. Remark 2.2. Assumptions 1.( c ), ( e ) are essential for Theorem 1 and Theorem 2 in that they ensure that ℓ ′ ( y i , f Θ ( t ) ( x i )) and σ ′′ ( · ) are uniformly lower bounded and therefore H j = X T D j X where ( D j ) ii ≥ -C/n . Without such bounds, the bulk spectrum of H j ( t ) can have unbounded support on the left, preventing any meaningful outlier/bulk decomposition and making the sharp phase transition in Theorems 1 and 2 impossible. ## 3 Consequences and discussion of main results ## 3.1 The single-neuron case Our general theory applies to the m = 1 (single-neuron) case, which is already interesting. Consider the case in which the target is a single-index model where the only latent direction is hard. Further assume a = 1 with bias b ∈ R . The Hessian with respect to θ ( t ) is $$\sum _ { i = 1 } ^ { n } Risk ( \theta ( t ) ) = - \frac { 1 } { n } \sum _ { i = 1 } ^ { n } ( c _ { i } ( x _ { i } )$$ Equivalently, the diagonal entries become ̂ g = ℓ ′′ σ ′ 2 + ℓ ′ σ ′′ . Define $$( 54 ) , G _ { t } := g ( V ( t ), h ( V *, e ) ).$$ Since m = 1, we can drop the subscript j and all results in the previous section hold with G j t replaced by ̂ G t ; we denote the corresponding Stieltjes transform by ̂ α t ( z ). In particular, the outlier equation (39) at fixed time reads $$\frac { \hat { G } _ { t } } { \tilde { U } ( \hat { z } ) } = 0 ,$$ with the large-time analogue obtained by replacing ( ̂ G t , ̂ α t ) by ( ̂ G ∞ , ̂ α ∞ ). The eigenvector-alignment and threshold statements follow by the same substitutions: replace G j t (resp. G j ∞ ) by ̂ G t (resp. ̂ G ∞ ) and α j t (resp. α j ∞ ) by ̂ α t (resp. ̂ α ∞ ) in (39)-(48) and in the corresponding alignment formulas. ## 3.2 Grokking As discussed in the introduction, our results suggest the emergence of grokking for δ > δ NN . (See Figure 2). 1. During the first training stage ( t = O (1)), the network cannot learn the hard directions, and therefore the parameters overfit slightly to the data whenever δ is finite (even when there is no easy directions). 2. During the second training stage ( t ≫ 1), as long as δ > δ NN , the emergence of negative directions of the Hessian drives the learning of hard features. Consequently the generalization error decreases. This phenomenon becomes less pronounced for δ ≫ δ NN because the generalization error decreases by uniform convergence arguments. We elaborate on this phenomenon further in the numerical experiments of Section 4. ## 3.3 Interpretation of the outliers equation Useful intuition into the outlier equations (39) and (47) is gained by recognizing that they coincide with the outlier equations of a simpler random matrix model. We will describe this interpretation for t fixed (Theorem 1); the analogous interpretation for t = ∞ (Theorem 2) is immediate. In order to motivate the simplified random matrix model, note that, by the DMFT formulas, there exists a measurable map φ t and an auxiliary vector ζ independent of U ∗ T H V ∗ such that $$V ( t ) = \phi _ { t } ( h ( V _ { x } , e ), s ).$$ Therefore, writing G j t = g ( V ( t ) , h ( V ∗ , ε ); j ), there exists a map Γ j t (obtained by composing g with φ t ) such that, jointly with U ∗ T H V ∗ , $$( V _ { \cdot } , ε ) , ( S ) , U ^ { H } V ^ { \cdot } , \{ \mu U ^ { H } T V ^ { \cdot } .$$ Motivated by these remarks, we define the random matrix model $$H ^ { j } ( t ) := - \sum _ { i = 1 } ^ { n } r ^ { j } _ { i , d } ( v _ { i } , e _ { i } , G _ { i } ) i < n L a v$$ Note that this differs from the actual distribution of H j ( t ) in two ways: ( i ) We replaced the diagonal terms g ( Θ ( t ) T x i , h ( Θ T ∗ x i , ε i ); j ) (which can be written in the form Γ j n,t ( h ( Θ T ∗ H x i , v n,i , ε n,i ) , ζ n,i ) for some Γ j n,t and ( v n,i , ε n,i , ζ n,i )) by its asymptotic counterpart Γ j t ( h ( Θ T ∗ H x i , v i , ε i ) , ζ i ) ; ( ii ) We drop dependencies between ( v n,i , ε n,i , ζ n,i ) and x i . By construction, Γ j t ( h ( V ∗ , ε ) , ζ ) d = G j t , so the bulk Stieltjes transform and edge equation of H sim j ( t ) obey the same generalized Marchenko-Pastur law as in Lemma 2.4. Moreover, the standard resolvent computation for outliers in this multi-index model yields the outlier condition $$\det ( - z I + E \left[ \frac { s } { \alpha ^ { 2 } + r ^ { 2 } } \right] ) = 0$$ which coincides with (39) upon identifying Γ j t d = G j t . Hence, the critical aspect ratio δ ∗ j ( t ) defined from (39) is exactly the spectral (BBP-type) threshold for detecting the hard subspace in the multiindex spiked covariance model H sim j ( t ) (with the additional requirement z < 0) (see [KZM25]). In words, we claim that the informative outlier eigenvalue of the Hessian block H j ( t ) are the same as the ones of the simplified model H sim j ( t ) which retains only the dependencies between Θ T ∗ H x i and x i . We emphasize that this does not mean that full random matrices H j ( t ) and H sim j ( t ) are equal in law, or even that they have the same outliers asymptotically. In particular, H sim j ( t ) can have at most r H outliers corresponding to the hard directions, whereas H j ( t ) might have additional outliers whose eigenvectors are not aligned with the hard directions. ## 3.4 Suboptimality among polynomial-time algorithms From previous discussion, we see that feature learning via gradient descent can be interpreted as a two-stage polynomial-time algorithm: 1. Preprocessing stage: Run O (1) steps of gradient descent, which learns the easy directions but cannot recover the hard subspace (Lemma 2.1). 2. Feature learning stage: Learn the hard directions via the spectral algorithm induced by the Hessian at the current iterate, as characterized in Section 2.5. We note that the first stage also induces the data preprocessing for the spectral feature learning algorithm in the second stage, as seen in the simplified model (58). This preprocessing is in general sub-optimal. This can be more clearly illustrated in the case when all latent directions are hard ( U H = R k ). In this setting, there are no easy directions to learn, and the simplified model (58) (with the same learning threshold) reduces to $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } R ^ { i , j } ( y _ { i } )$$ Here ζ i is independent random noise that does not depend on the data ( x i , y i ), and without loss of generality we can assume it Unif ([0 , 1]). The function Γ j t depends on the network architecture σ , a j , b j ), the loss function ( ℓ ), the training hyperparameters ( η ). A general spectral method uses instead an arbitrary preprocessing function Ψ [MM19, KZM25]: $$H _ { y , s } = - \sum _ { i = 1 } ^ { n } \varphi ( y _ { i } , G _ { i } )$$ The corresponding spectral threshold δ sp (Ψ) is typically minimized at a unique Ψ opt (up to measure zero modifications) and hence in general we have the strict inequality: $$\delta ^ { \ast } ( t ) = \delta _ { sp } ( T ^ { j } ) > \delta _ { opt } : = i n f \delta _ { sp } ( \Psi$$ ̸ and in fact, inf t ≥ 0 δ ∗ j ( t ) > δ opt . As a simple illustration, consider a single-neuron network ( m = 1) learning a single-index model ( k = 1), with a = 1 and b = 0. Suppose the activation function σ is locally quadratic near zero (i.e., σ (0) = σ ′ (0) = 0 and σ ′′ (0) = 0), and the weight vector is initialized with ∥ θ (0) ∥ = d -C for some large constant C > 0. In this regime, the gradient remains negligible for any O (1) steps, and the dynamics is dominated by the Hessian. Consequently, gradient descent effectively performs a power method on the Hessian at zero: $$\sum _ { i = 1 } ^ { n } \ell ( y _ { i } ; 0 ) x _ { i } z _ { i }$$ The feature learning threshold δ ∗ j ( t ) coincides with the spectral threshold δ ∗ sp ( ℓ ′ ( · , 0)) for the preprocessing function y ↦→ ℓ ′ ( y, 0). This is in general not optimal (consider for instance h ( z, ε ) = z 2 ), demonstrating that inequality (62) can be strict. ## 4 Examples and numerical illustrations This section instantiates the theory on the following data distribution with k = 1 ('real phase retrieval') $$y = h ( z , e ) = z ^ { 2 } ,$$ In this example, the information-theoretic threshold and the algorithmic threshold are known to be δ IT = δ alg = 1 / 2 [MM19]. We derive explicit thresholds for this task and compare them with numerical experiments with various activation functions. As we will see, the predicted thresholds precisely match empirical phase transitions and explain the observed grokking behavior when δ > δ NN . We extend the numerical analysis to multi-neuron networks in Appendix E. ## 4.1 General setting In this setting, the hard subspace coincides with the entire latent space, i.e. r = k = 1. Consequently, Lemmas 2.1 and 2.2 imply that no correlation with θ ∗ can be learned in finitely many gradient steps. We use the standard Huber loss with parameter M > 0, which ensures bounded diagonal terms in the Hessian and thus a finite left edge of the bulk spectrum: $$y ( y , z ) = \{ \frac { 1 } { M } ( x - y ) ^ { 2 } , \frac { 1 } { M } ( z - y ) ^ { 2 } M ^ { 2 } , \frac { 1 } { z } -$$ For m = 1 with ( a, b ) = (1 , 0), the empirical risk and Hessian read $$\begin{aligned} Risk ( \theta ) = - 1 _ { n } \sum _ { i = 1 } ^ { n } g ( \theta ^ { T } x ; i , \sigma ^ { 2 } _ { i } ) \\ &= - 1 _ { n } \sum _ { i = 1 } ^ { n } g ( \theta ^ { T } x ; i , \sigma ^ { 2 } _ { i } ) \end{aligned}$$ where $$f ( u ) ^ { 2 } + l ^ { " } ( y , \sigma ( u ) ) d ^ { " } ( u ) .$$ Using the discrete-time DMFT random variables ( V ( t ) , V ∗ ) from Lemma 2.1 (here y = h ( V ∗ , ε ) = V 2 ∗ ), we define $$\sum _ { i = 1 } ^ { n } g ( V ( t ) ; V ^ { 2 } _ { i } ) , G _ { \infty } := g ( V$$ Formulas for the asymptotic spectrum. The limiting spectral distribution of the Hessian follows a generalized Marchenko-Pastur law with Stieltjes transform ̂ α t ( z ) determined by (cf. Lemma 2.4): $$z + \frac { 1 } { a _ { i } ( z ) } = s \cdot E [ \frac { G _ { t } } { \delta + G _ { t } a _ { i } ( z ) } ]$$ we denote the left edge of the support of the asymptotic ESD by c ( t ). In the large-time limit, replace ( ̂ G t , ̂ α t , c ( t )) by ( ̂ G ∞ , ̂ α ∞ , c ( ∞ )) as in (45)-(46). Since r = 1, there is only one hard direction. The outlier equation (39) specializes to the scalar equation $$\frac { \delta G _ { t } V ^ { 2 } } { \delta G _ { t } \Delta t ( z ) } = 0 , z < \min ( c ( t ), + G _ { t } \Delta t ( z ) )$$ Figure 4: Phase transitions for GeLU with single-neuron learner. Left: Success rate with sigmoid fits; each data point shows mean ± standard deviation across 60 independent runs. Right: Final correlation with median and 30th/70th percentile error bars (smooth splines fitted). Both metrics exhibit sharp phase transitions matching the predicted threshold δ ∗ ≈ 6 . 0 (dash-dot line). <details> <summary>Image 4 Details</summary>  ### Visual Description ## Line Graphs: Success Rate vs. δ = n/d and Final Correlation vs. δ = n/d ### Overview The image contains two side-by-side line graphs. The left graph plots "Success Rate (Threshold = 0.5)" against δ = n/d (x-axis), while the right graph plots "Final Correlation" against the same x-axis. Both graphs include multiple data series differentiated by line color/shade and error bars. A vertical dashed line at δ = 6 is labeled as the "Predicted threshold" in both graphs. --- ### Components/Axes #### Left Graph (Success Rate) - **X-axis**: δ = n/d (ranges from 3 to 9, labeled in increments of 1). - **Y-axis**: Success Rate (Threshold = 0.5) (ranges from 0.0 to 1.0, labeled in increments of 0.2). - **Legend**: Located at the bottom-left corner. Entries: - `d = 1000` (light blue line with shaded area). - `d = 2000` (medium blue line with shaded area). - `d = 3000` (dark blue line with shaded area). - **Dashed Line**: Vertical line at δ = 6 (black, labeled "Predicted threshold"). #### Right Graph (Final Correlation) - **X-axis**: δ = n/d (same scale as left graph). - **Y-axis**: Final Correlation (ranges from 0.0 to 1.0, labeled in increments of 0.2). - **Legend**: Located at the bottom-right corner. Entries: - `d = 1000` (light yellow line with shaded area). - `d = 2000` (orange line with shaded area). - `d = 3000` (dark orange line with shaded area). - `d = 4000` (red line with shaded area). - **Dashed Line**: Vertical line at δ = 6 (black, labeled "Predicted threshold"). --- ### Detailed Analysis #### Left Graph (Success Rate) - **Trend**: All lines show a sigmoidal increase. Success rate remains near 0 for δ < 5, then rises sharply after δ = 6. - **Data Points**: - **d = 1000**: Crosses threshold (0.5) at δ ≈ 6.5. Final value ≈ 0.95 at δ = 9. - **d = 2000**: Crosses threshold at δ ≈ 6.3. Final value ≈ 0.98 at δ = 9. - **d = 3000**: Crosses threshold at δ ≈ 6.2. Final value ≈ 0.99 at δ = 9. - **Error Bars**: Smaller for higher d values (e.g., d = 3000 has tighter error margins than d = 1000). #### Right Graph (Final Correlation) - **Trend**: All lines plateau near 1.0 after δ = 6. Correlation remains low (≈0.1–0.2) for δ < 5. - **Data Points**: - **d = 1000**: Crosses threshold at δ ≈ 6.3. Final value ≈ 0.98 at δ = 9. - **d = 2000**: Crosses threshold at δ ≈ 6.1. Final value ≈ 0.99 at δ = 9. - **d = 3000**: Crosses threshold at δ ≈ 6.0. Final value ≈ 1.0 at δ = 9. - **d = 4000**: Crosses threshold at δ ≈ 5.9. Final value ≈ 1.0 at δ = 9. - **Error Bars**: Larger for smaller d values (e.g., d = 1000 has wider error margins than d = 4000). --- ### Key Observations 1. **Threshold Consistency**: Both graphs show the predicted threshold at δ = 6, but actual crossing points vary slightly by d value. 2. **d Value Impact**: Larger d values (e.g., 4000) achieve higher success rates and correlations at lower δ values compared to smaller d values (e.g., 1000). 3. **Convergence**: Lines for higher d values cluster closer to the predicted threshold line, suggesting diminishing returns beyond δ = 6. --- ### Interpretation - **Success Rate**: The threshold δ = 6 marks the point where success rate exceeds 0.5. Larger d values (more data) achieve this threshold earlier, indicating improved performance with increased data size. - **Final Correlation**: Similarly, higher d values reach near-perfect correlation (≈1.0) sooner, suggesting that δ = 6 is a critical inflection point for model reliability. - **Error Margins**: Smaller d values exhibit greater variability (wider error bars), reinforcing the need for larger datasets to stabilize results. - **Practical Implication**: The predicted threshold δ = 6 may serve as a guideline for balancing computational efficiency (lower δ) and performance (higher success/correlation). However, real-world application should account for d-dependent variability. </details> and the finite-time weak-recovery threshold is $$\delta ^ { \ast } ( t ) = i n f \{ s \in ( 0 , \infty ) : z < m i n ( t ) , 0 \}$$ The obvious analogous of these formulas hold for t = ∞ . Computing the threshold. Since ̂ α ′ t ( z ) > 0 for z < c ( t ) by standard properties of Stieltjes transforms, differentiating yields ∂ ∂z S t ( z ; δ ) < 0 on ( -∞ , min( c ( t ) , 0)). Hence, for each fixed δ , S t ( z ; δ ) is strictly decreasing in z , and a root exists for some z < min( c ( t ) , 0) if and only if S t ( z † ( t ); δ ) < 0 at the boundary value $$z ( t ) := \min ( c ( t ), 0 ) , \overline { a } _ { t } = \lim _ { t \rightarrow + \infty } z ( t )$$ Therefore, the finite-time threshold δ ∗ ( t ) solves the boundary equation $$E \left[ \frac { s G _ { i } V ^ { 2 } } { s + G _ { t } a _ { i } } \right ] = z ^ { t } ( t )$$ for δ ∈ (0 , ∞ ]. If (73) has no positive solution, then δ ∗ ( t ) = + ∞ . The large-time threshold δ ∗ ( ∞ ) is computed analogously by replacing ( ̂ G t , ̂ α t , c ( t )) with ( ̂ G ∞ , ̂ α ∞ , c ( ∞ )). ## 4.2 Results for GeLU and Quad activations We begin by reporting results for Huber loss with M = 1 and GeLU activation $$\beta = \frac { 1 } { 3 } , (74)$$ where Φ( · ) is the CDF of the standard normal distribution. For Quad activation we use α = 9, β = 1 / 3, We initialize θ (0) as a uniformly random unit vector, and run full-batch gradient descent with learning rate η = 1 . 5 for GeLU and η = 0 . 25 for Quad (these were selected by optimizing accuracy Figure 5: Phase transitions for Quad with single-neuron learner. Left: Success rate with sigmoid fits; each data point shows mean ± standard deviation across 60 independent runs. Right: Final correlation with median and 30th/70th percentile error bars (smooth splines fitted). Both metrics exhibit sharp phase transitions matching the predicted threshold δ ∗ ≈ 3 . 6 (dash-dot line). <details> <summary>Image 5 Details</summary>  ### Visual Description ## Line Graphs: Success Rate vs. δ = n/d and Final Correlation vs. δ = n/d ### Overview The image contains two side-by-side line graphs. The left graph plots "Success Rate (Threshold = 0.5)" against δ = n/d, while the right graph plots "Final Correlation" against the same δ = n/d. Both graphs include multiple data series (d = 1000, 2000, 3000, 4000) and a predicted threshold line at δ = 4. Error bars indicate variability in measurements. --- ### Components/Axes - **Left Graph**: - **Y-axis**: Success Rate (Threshold = 0.5), scaled from 0.0 to 1.0. - **X-axis**: δ = n/d, scaled from 2 to 7. - **Legend**: - Light blue: d = 1000 - Medium blue: d = 2000 - Dark blue: d = 3000 - Very dark blue: d = 4000 - **Predicted Threshold**: Dashed vertical line at δ = 4. - **Right Graph**: - **Y-axis**: Final Correlation, scaled from 0.0 to 1.0. - **X-axis**: δ = n/d, scaled from 2 to 7. - **Legend**: - Yellow: d = 1000 - Orange: d = 2000 - Red: d = 3000 - Pink: d = 4000 - **Predicted Threshold**: Dashed vertical line at δ = 4. --- ### Detailed Analysis #### Left Graph (Success Rate) - **Trends**: - All data series show a sigmoidal increase in success rate as δ increases. - Lower d values (e.g., d = 1000) exhibit gradual curves with slower convergence to 1.0. - Higher d values (e.g., d = 4000) show steeper rises, reaching near 1.0 by δ ≈ 5. - Error bars decrease in size as δ approaches 7, indicating reduced variability at higher δ. - **Key Data Points**: - At δ = 3: Success rates range from ~0.3 (d = 1000) to ~0.5 (d = 4000). - At δ = 5: Success rates exceed 0.9 for all d values. - Predicted threshold at δ = 4 aligns with the inflection point where success rates cross 0.5. #### Right Graph (Final Correlation) - **Trends**: - All data series peak sharply near δ = 4, then decline steeply. - Higher d values (e.g., d = 4000) achieve higher peak correlations (~0.95) but drop more abruptly post-δ = 4. - Lower d values (e.g., d = 1000) have lower peaks (~0.8) and gentler declines. - Error bars are largest near δ = 4, reflecting high variability at the threshold. - **Key Data Points**: - At δ = 3: Correlations range from ~0.2 (d = 1000) to ~0.5 (d = 4000). - At δ = 5: Correlations drop to ~0.1 (d = 1000) and ~0.3 (d = 4000). - Predicted threshold at δ = 4 coincides with the peak of all curves. --- ### Key Observations 1. **Threshold Behavior**: The predicted threshold at δ = 4 marks a critical transition point. Success rates and correlations peak here, after which performance declines. 2. **d-Dependent Performance**: - Larger d values (4000) achieve higher success rates and correlations but are more sensitive to δ > 4. - Smaller d values (1000) show slower convergence and lower peak performance. 3. **Error Patterns**: Variability increases near δ = 4, suggesting uncertainty in measurements at the threshold. --- ### Interpretation The data demonstrates that δ = n/d acts as a critical parameter influencing system performance. The predicted threshold at δ = 4 likely represents an optimal balance between resource allocation (n/d) and performance. Higher d values improve robustness and peak performance but amplify sensitivity to deviations beyond δ = 4. This suggests that system design should prioritize maintaining δ ≤ 4 while maximizing d to enhance reliability. The error bars highlight the need for larger sample sizes (d) to reduce uncertainty, particularly near the threshold. </details> over a grid of values of η ) for at most 1000 iterations. We terminate early once the correlation reaches ρ ( t ) ≥ 0 . 9, after which we run for 100 additional steps to ensure sufficient training. We also stop if the algorithm plateaus for too long, i.e., if the correlation changes by less than 0 . 01 over 500 consecutive steps. Simulation results. We consider d ∈ { 1000 , 2000 , 3000 , 4000 } with varying sample sizes n to explore different sampling ratios δ = n/d . For each configuration ( d, δ ), we perform 60 independent runs with different random seeds for θ ∗ , θ (0), and the data { x i , y i } n i =1 . We measure the correlation between the learned parameter and the ground truth via the normalized inner product $$\rho ( t ) := \left | \frac { \theta ( t ) } { 1 0 ( t ) } \right | ^ { 2 }$$ evaluated at the final time t = T (maximum 1000 steps). We report two key metrics: ( i ) the success rate , defined as the proportion of runs achieving ρ ( T ) ≥ 0 . 5, and ( ii ) the final correlation ρ ( T ) itself. For success rates, we fit sigmoid curves; for correlations, we report medians with 30th/70th percentile error bars. Figures 4 and 5 report simulation results for GeLU and Quad respectively. We plot empirical success rate and final correlation versus δ = n/d across four dimensions. Both metrics exhibit a sharp phase transition around δ = 6 (for GeLU ) and δ = 3 . 6 (for Quad ) : the success rate (left panel) jumps from near zero to near one, while the final correlation (right panel) transitions from ρ ( T ) ≈ 0 to ρ ( T ) ≈ 1. These behaviors match the predictions from high-dimensional asymptotics that are reported as vertical lines. We discuss the computation of these predictions next. Theoretical predictions. We compute the predicted thresholds δ ∗ ( t ) for t ∈ { 0 , 1 , . . . , 25 } (for GeLU ) or t ∈ { 0 , 1 , . . . , 25 } (for Quad ) as follows. For each t and each δ , we run the discrete-time DMFT recursion (Section 2.2) forward to time t , drawing fresh Monte Carlo samples (with sample size N = 9 × 10 5 ) at each step to estimate the expectations in (25) and update the covariance and response kernels ( C θ , C ℓ , R θ , R ℓ ). We then check whether an outlier eigenvalue emerges to the Figure 6: Predicted feature learning thresholds δ ∗ ( t ) with m = 1. Left: GeLU , m = 1. Right: Quad . Circles: mean over 3 or 20 seeds (respectively for GeLU and Quad ); error bars: standard deviation. Solid line and shaded band: spline interpolations of the mean and mean ± standard deviation. <details> <summary>Image 6 Details</summary>  ### Visual Description ## Line Graphs: δ*(t) vs Time Step t ### Overview Two line graphs are presented side-by-side, depicting the evolution of δ*(t) over time steps (t). Both graphs show a sigmoidal (S-shaped) curve with a plateau phase, accompanied by shaded regions representing variability (Mean ± Std). The left graph spans t=0 to t=25, while the right graph extends to t=30. --- ### Components/Axes - **X-axis**: "Time Step t" (integer values, 0–25 for left graph; 0–30 for right graph). - **Y-axis**: "δ*(t)" (continuous values). - Left graph: 4.5–6.0 - Right graph: 2.0–3.5 - **Legends**: - **Blue line**: "Mean" (solid line). - **Shaded blue region**: "Mean ± Std" (standard deviation around the mean). --- ### Detailed Analysis #### Left Graph (t=0–25) - **Trend**: δ*(t) increases rapidly from ~4.5 at t=0 to ~6.0 by t=10, then plateaus with minimal growth. - **Key Data Points**: - t=0: δ*(t) ≈ 4.5 - t=5: δ*(t) ≈ 5.3 - t=10: δ*(t) ≈ 5.8 - t=20: δ*(t) ≈ 5.95 - t=25: δ*(t) ≈ 6.0 - **Variability**: Shaded region (Mean ± Std) remains narrow (~±0.1) throughout, indicating low variance. #### Right Graph (t=0–30) - **Trend**: δ*(t) rises from ~2.0 at t=0 to ~3.5 by t=20, then plateaus. - **Key Data Points**: - t=0: δ*(t) ≈ 2.0 - t=5: δ*(t) ≈ 2.6 - t=15: δ*(t) ≈ 3.2 - t=25: δ*(t) ≈ 3.45 - t=30: δ*(t) ≈ 3.5 - **Variability**: Shaded region (~±0.1) is consistent, mirroring the left graph’s low variance. --- ### Key Observations 1. **Sigmoidal Behavior**: Both graphs exhibit a rapid initial increase followed by a plateau, suggesting a system approaching equilibrium. 2. **Plateau Values**: - Left graph plateaus at δ*(t) ≈ 6.0 (higher magnitude). - Right graph plateaus at δ*(t) ≈ 3.5 (lower magnitude). 3. **Consistency**: Shaded regions (Mean ± Std) are tightly clustered, indicating minimal variability in measurements. 4. **Time to Plateau**: - Left graph stabilizes by t=20. - Right graph stabilizes by t=25–30. --- ### Interpretation The graphs likely represent a dynamic system (e.g., biological, chemical, or engineering process) where δ*(t) evolves over time. The sigmoidal trend implies a transition from an initial state to a stable equilibrium. The higher plateau in the left graph suggests a system with greater capacity or different parameters compared to the right graph. The tight standard deviation bands indicate reliable, reproducible measurements. The extended time horizon in the right graph (t=30) may reflect a slower convergence rate or a different experimental setup. These results could inform optimization strategies, such as adjusting parameters to achieve desired δ*(t) values more efficiently. </details> left of the bulk spectrum by testing if the outlier equation (39) admits a solution z < c ( t ) -0 . 01, where c ( t ) is the left bulk edge (cf. (37)); the gap of 0 . 01 ensures numerical stability. We use binary search over δ (tolerance 0 . 01) to locate the critical threshold, and average results over 3 independent random seeds. The resulting thresholds are plotted in Figure 6. For GeLU , the threshold exhibits rapid initial growth from δ ∗ (0) ≈ 4 . 35 to δ ∗ (10) ≈ 5 . 85, and then appears to converge to a limiting value. To estimate this limit, we fit a degree-4 polynomial in 1 /t to the data, obtaining δ ∗ ( ∞ ) ≈ 6 . 05. We report the predictions for δ ∗ ( ∞ ) in Figures 4 and 5 and observe that they closely match the empirical phase transition. For Quad , the threshold exhibits rapid initial growth from δ ∗ (0) ≈ 2 . 0 to δ ∗ (10) ≈ 3 . 2, and then appears to converge to a limiting value. To estimate this limit, we fit a degree-4 polynomial in 1 /t to the data, obtaining δ ∗ ( ∞ ) ≈ 3 . 6. In summary, the empirical results confirm our expectation that, as n, d → ∞ , the threshold δ ∗ ( ∞ ) controls the learning of latent directions. However, the separation between the first phase of training ( t = O (1)), which in this case corresponds to pure overfitting, and the second one ( t = Θ(log d )) 7 , on which feature learning takes place, is relatively modest for moderate d . Since δ ∗ ( t ) < δ ∗ ( ∞ ) in this case, it is possible that feature learning occurs at smaller δ for moderate d . ## 4.3 Training dynamics and grokking We already illustrated the grokking phenomenon in Figure 2 of the introduction. We discussed its connection with the two stages learning dynamics in Section 3.2. Figure 7 illustrates how the learning dynamics change as a function of δ , for a GeLU neuron with d = 5000 and η = 0 . 5 on noiseless phase retrieval. Recall that the predicted feature learning threshold is δ NN ≈ 6 . 0. For δ = 5 (below threshold), no feature learning takes place: the correlation ρ ( t ) remains near zero throughout training. At δ = 10, slightly above threshold, grokking emerges but requires extended training time ( t ≈ 300, compared to t ≈ 80 in Fig. 2). As δ increases to 15, the transition accelerates (occurring around t ≈ 100), while grokking can still be observed. Finally, at δ = 50, learning becomes rapid with negligible generalization gap, and the grokking phenomenon becomes less pronounced. In total, these observations confirm our theoretical predictions: ( i ) no 7 The time Θ(log d ) is known to be the saddle point escape time; see, e.g., [AGJ21]. Figure 7: Learning dynamics of a GeLU neuron for varying sample ratios δ = n/d ∈ { 5 , 10 , 15 , 50 } , with d = 5000, η = 0 . 5, and noiseless phase retrieval target. Left: Correlation ρ ( t ) between learned and target parameters (log time scale is used for better visualization). Center: Training risk (solid) and test risk (dashed). Right: Generalization gap Risk test -Risk . Grokking is observed when δ is moderately above the threshold: training risk first drops while test risk remains high before eventually decreasing. The generalization gap peaks during the overfitting phase and vanishes once feature learning completes. <details> <summary>Image 7 Details</summary>  ### Visual Description ## Line Graphs: Model Performance Metrics vs Training Time ### Overview Three line graphs depict the evolution of three metrics (ρ(t), Train/Test Risk, Generalization Gap) across training time (t) for four hyperparameter values (δ = 5, 10, 15, 50). Each graph shows distinct trends with δ values represented by different colors. ### Components/Axes 1. **X-axis**: Training Time (t) ranging from 0 to 300 in all graphs 2. **Y-axes**: - Left graph: ρ(t) (0-1 scale) - Middle graph: Train/Test Risk (0.4-0.8 scale) - Right graph: Generalization Gap (0-0.4 scale) 3. **Legends**: - Bottom of each graph - Colors: δ=5 (yellow), δ=10 (green), δ=15 (cyan), δ=50 (blue) - Line styles: Solid for ρ(t), dashed for Risk, dotted for Gap ### Detailed Analysis **Left Graph (ρ(t))**: - All δ values start at 0 - Sharp rise begins ~t=40 for all δ - δ=50 reaches 1.0 by t=80 (highest) - δ=5 remains lowest (plateaus ~0.2 by t=160) - Approximate values: - δ=5: 0.2 (t=160) - δ=10: 0.3 (t=160) - δ=15: 0.4 (t=160) - δ=50: 1.0 (t=80) **Middle Graph (Train/Test Risk)**: - All δ start ~0.7 at t=0 - δ=50 drops fastest (0.5 by t=100) - δ=5 remains highest (0.6 by t=100) - Approximate values: - δ=5: 0.6 (t=100) - δ=10: 0.55 (t=100) - δ=15: 0.52 (t=100) - δ=50: 0.5 (t=100) **Right Graph (Generalization Gap)**: - All δ start at 0 - δ=5 peaks highest (~0.35 at t=100) - δ=50 peaks lowest (~0.15 at t=100) - All δ show decline after t=100 - Approximate values: - δ=5: 0.35 (t=100) - δ=10: 0.3 (t=100) - δ=15: 0.25 (t=100) - δ=50: 0.15 (t=100) ### Key Observations 1. **δ Inversely Correlated with Risk**: Higher δ values achieve lower risk faster (δ=50 reaches 0.5 risk at t=100 vs δ=5 at 0.6) 2. **Generalization Gap Peaks**: All δ show peak generalization gap ~t=100, with δ=5 having largest gap (0.35) 3. **ρ(t) Saturation**: δ=50 saturates ρ(t) earliest (t=80), while δ=5 remains unsaturated at t=160 4. **Consistent Trends**: All δ values show similar shape patterns across metrics, with δ=50 consistently outperforming others ### Interpretation The data suggests δ acts as a regularization parameter: - Higher δ values (50) reduce overfitting (lower risk, smaller generalization gap) - Lower δ values (5) show more pronounced overfitting (higher risk, larger gap) - The peak generalization gap at t=100 indicates a critical point where models begin overfitting before δ regularization takes full effect - ρ(t) saturation patterns correlate with risk reduction timing, suggesting ρ(t) may represent model confidence or density of learned features - The consistent δ ordering across all metrics implies a unified mechanism governing model behavior across training phases </details> learning for δ < δ NN ; ( ii ) grokking emerges for δ > δ NN ; ( iii ) grokking time increases as δ ↓ δ NN ; ( iv ) generalization gap diminishes for δ ≫ δ NN . The same scenario holds with other activation functions as shown in Appendix E. ## 5 Outline of the proofs This section outlines the proofs of Lemma 2.4 and Theorem 1. We organize the arguments in three steps: Step 1: Gaussian conditioning; Step 2: Detecting outlier eigenvalues via resolvent expansion and concentration; Step 3: Verifying that eigenvectors corresponding to specific outliers align with hard directions. Full proofs are given in Appendices B and C. ## 5.1 Step 1: Gaussian conditioning The starting point is the j -th Hessian diagonal block of Eq. (35), which we rewrite: $$H _ { j } ( t ) = - \sum _ { i = 1 } ^ { n } g ( \theta ( t ) X ,$$ where G j ( t ) = (1 /n ) diag ( g ( Θ ( t ) T x 1 , y 1 ; j ) , . . . , g ( Θ ( t ) T x n , y n ; j ) ) and we recall that $$g ( z , y ; j ) = f ( y , \frac { 1 } { m } z _ { i }$$ $$g ( z , y ) := e ^ { i ( y , o ( z ) ) }$$ Isolating the dependence on the GD trajectory. Although H j ( t ) depends on X in a complicated nonlinear way (through { Θ ( t ) T x i } i ≤ n and { y i } i ≤ n , where Θ ( t ) itself depends on X , y via the GD update rule), conditioning on a certain GD trajectory is equivalent to condition on a certain number of linear functions of X . Since X is Gaussian, its conditional distribution given a set of linear observations can be written explicitly. In the present case X decomposes into the sum of a low-rank component, which is measurable on the GD trajectory, and a high-rank component which is (aside from a projection step) independent from the GD trajectory. This decomposition is the foundation for both outlier detection (Step 2) and eigenvector analysis (Step 3). Gaussian conditioning decomposition. Define the parameter and gradient matrices $$\begin{aligned} \theta ( t ) := \left[ \theta _ { * } , \theta ( 0 ) , \ldots \right] \\ F ( 0 ; t ) := \frac { 1 } { \sqrt { n } } \left[ F ( X \theta ( t ) ) \right] \end{aligned}$$ By the gradient descent update rule, we have X T F (0; t ) ∈ span( Θ ( ∗ ; t )). Let F t be the σ -algebra generated Θ ( ∗ ; t ), F (0; t ). Conditional on F t , we have the distributional equivalence 8 $$X _ { | r _ { z } | } = P ^ { d } F ( 0 ; t ) X _ { new } P ^ { d } \theta ( s ; t ) + X _ { F }$$ where X new ∈ R n × d has i.i.d. N (0 , 1) entries independent of all other randomness, and P M denotes orthogonal projection onto span( M ). One of the crucial observations is that P F (0; t ) XP ⊥ Θ ( ∗ ; t ) = 0 due to the gradient descent structure. (Details in Appendix B, Lemma B.1.) Substituting (80) into (76) and expanding yields (using the shorthands Y ( t ) = P ⊥ F (0; t ) X new P ⊥ Θ ( ∗ ; t ) , and P t := P Θ ( ∗ ; t ) ) $$H _ { j } ( t ) = Y ( t ) ^ { d } G _ { j } ( t ) Y ( t ) = X ^ { T } _ { new } G _ { j } ( t ) X _ { new }$$ where the low-rank component has rank dependent uniquely on the number of GD steps t , the number of neurons m and the number of target directions k , but not on n, d . Bulk spectrum. By Cauchy interlacing, the low-rank perturbation in (81) do not affect the limiting spectral distribution. Since X new has i.i.d. Gaussian entries independent of G j ( t ), the asymptotic empirical spectral distribution of X T new G j ( t ) X new follows a general Marchenko-Pastur law. This asymptotic distribution depends on the aspect ratio δ = lim n,d →∞ n/d and the asymptotic distribution of G j ( t ) (given by discrete-time DMFT in Lemma 2.1). This establishes Lemma 2.4. We refer to Appendix B.3 for details omitted in this outline. ## 5.2 Step 2: Outlier matrix and concentration to deterministic limit We now locate outlier eigenvalues, namely eigenvalues that do not converge to the support of the asymptotic empirical spectral distribution µ j ∞ ( t ) (which, as we saw, follows a general MarchenkoPastur law) that from the bulk. Applying the Woodbury matrix identity to decomposition (81), we express the resolvent as $$\frac { ( H _ { j } ( t ) - z I ) ^ { - 1 } } { M _ { j } ( z ; t ) }$$ 8 For two random variables X , X ′ we write X ∣ ∣ F t d = X ′ if E [ φ ( X ) Z ] = E [ φ ( X ′ ) Z ] for any bounded measurable function φ and any bounded random variable Z measurable on F t . where R 0 ( z ) := ( Y ( t ) T G j ( t ) Y ( t ) -z I ) -1 is the bulk resolvent, Z ( t ) := X Θ ( ∗ ; t ) ( Θ ( ∗ ; t ) T Θ ( ∗ ; t ) ) -1 , and the matrix M j ( z ; t ) ∈ R 2( k + m ( t +1)) × 2( k + m ( t +1)) has the explicit form $$( t ) ^ { I } \left[ Z ( t ) ^ { T } G _ { j } ( t ) Z _ { I } \right]$$ We refer to M j ( z ; t ) (and its n, d →∞ limit defined below) as to the 'outlier matrix.' Since all eigenvalues of Y ( t ) T G j ( t ) Y ( t ) converge to the support of the general MarchenkoPastur law, it follows from (82) that z is an outlier eigenvalue of H j ( t ) only if det( M j ( z ; t )) = 0. Crucially, M j ( z ; t ) is low-dimensional. Namely, it has dimension 2( k + m + mt ) × 2( k + m + mt ), independent of the ambient dimensions n, d . As a consequence, it concentrates (for z in a suitable domain) around a deterministic matrix, and it is therefore sufficient to study the latter. Deterministic limit of the outlier matrix. Wenext describe the deterministic limit of M j ( z ; t ). Without loss of generality, we consider the case where the hard directions are the first r target directions. Lemma 5.1 (Restatement of Lemma C.1) . For z < min(0 , c j ( t )) , under the assumptions of Theorem 1, the outlier matrix (83) converges in probability to a deterministic limit $$\lim _ { n \rightarrow \infty } M _ { j } ( x ; t ) = M ^ { p - 1 } F ( z , t )$$ where O T t := [Θ T ∗ , Θ(0) T , . . . , Θ( t ) T ] collects the state-evolution variables. Furthermore, denoting by V ∗ ⩽ the first r coordinates of V ∗ , the matrix S j ( z ; t ) ∈ R ( k + m + mt ) × ( k + m + mt ) is block-diagonal : $$\frac { S _ { j } ( z ; t ) } { 0 } = \left[ - E \left[ \delta + c \right] \right] .$$ Lemma 5.1 reduces the problem of determining the outlier eigenvalues of the high-dimensional random matrix H j ( t ) ∈ R d × d to analyzing the deterministic matrix M ∞ j ( z ; t ), whose dimension is independent of n, d . The proof proceeds in two stages. First, we condition on the GD trajectory, i.e. on the σ -algebra F t . Using the decomposition (80), we can think as ( X , ε ) fixed and treat the independent Gaussian matrix X new from (80) as the only source of randomness. We apply leave-one-out expansions and Hanson-Wright inequalities to show that each entry concentrates to a conditional limit involving the Stieltjes transform of the bulk distribution. Second, we integrate over ( X , ε ). Namely, we use discrete-time DMFT (Lemmas 2.1 and 2.2), replacing to show that quantities measurable over F t arising from the first step concentrate around expectations over the state-evolution variables ( V ∗ , V (0) , . . . , V ( t ) , ε ) and (Θ ∗ , Θ(0) , . . . , Θ( t )). (Full proof in Appendix C, Sections C.2, C.3.) Block-diagonal structure and the hard block equation. The block-diagonal form (85) turns out to be a consequence of the fact that U H is the (maximal) hard subspace property and hence, by Lemma 2.1, hard directions V ∗ ⩽ are uncorrelated with all quantities observable from training. This decouples the hard block (upper-left r × r ) from the rest. Lemma 2.2 establishes that E [ O t O T t ] as a similar block-diagonal structure. Using this block-diagonal structure, and applying the Schur complement formula to (84), the equation det( M ∞ j ( z ; t )) = 0 factors as $$0 . \det ( M _ { j , i } H ( z ; t ) ) \cdot \det ( M _ { j , i } R ( z ; t ) ) =$$ for suitably defined matrices M ∞ j, R ( z ; t ), M ∞ j, H ( z ; t ). In particular, the equation det ( M ∞ j, H ( z ; t ) ) reads $$\det ( - z I _ { R } + E \frac { \delta G ^ { i } _ { T } } { \delta t } V _ { x } V ^ { i } _ { T } ) = 0 .$$ We will refer to this as to the 'hard block equation'. Note that this is precisely the outlier equation (39) in Theorem 1. Finally notice that outliers can converge not only to solutions of the hard block equation but also to other solutions of Eq. (86), namely to solutions of det ( M ∞ j, R ( z ; t ) ) = 0. However, in Step 3 we show that only solutions to (87) yield eigenvectors aligned with hard directions. Finally, we use Rouch´ e's theorem from complex analysis (Appendix C.4) to show that the neighborhood of solutions of det( M ∞ j ( z ; t )) must contain actual eigenvalues of the Hessian H j ( t ) (with matching multiplicity) converging to those solutions. ## 5.3 Step 3: Eigenvector alignment via residue calculation Finally, we determine which outlier eigenvectors align with the hard subspace Θ ∗ H = [ θ ∗ 1 , . . . , θ ∗ r ]. For simplicity, in this outline we assume that the zeros of det( M ∞ j ( z ; t )) = 0 are simple (multiplicity one). The general case follows analogously. Let z ∗ be a simple zero of det( M ∞ j ( z ; t )). By Rouch´ e's theorem (Appendix C.4), there exists a unique eigenvalue z n of H j ( t ) satisfying | z n -z ∗ | p - - → 0 with corresponding unit eigenvector ξ n . To access the eigenvector, we use the spectral projector $$= \sum _ { n = 1 } ^ { p _ { n , d } } - \frac { 1 } { 2 \pi i } \int _ { 0 } ^ { H } ( t ) - z I - 1 d z$$ where γ is a contour around z ∗ enclosing no other zero of det( M ∞ j ( z ; t )) and no part of the support of the general Marchenko-Pastur law µ ∞ j ( t ). The squared correlation with hard directions is $$( 8 9 )$$ The contour integral formulation converts the eigenvector alignment question into a residue calculation that can be analyzed using the limiting outlier matrix. Residue formula and block-diagonal structure. Substituting the resolvent expansion (82) into (88), we observe that the bulk resolvent R 0 ( z ) is analytic inside γ and contributes no residue. The only pole arises from M j ( z ; t ) -1 , reducing the contour integral to $$\begin{aligned} T_{R_0}(z) \Theta(x; t) Y^* GZ M_j(z; t) - 1 & = \frac{\partial T_R(e)}{\partial x} ( e ) \\ & = \frac{\partial T_R(e)}{\partial x} ( e ) \end{aligned}$$ By Lemma 5.1, we replace M j ( z ; t ) -1 by its limit M ∞ j ( z ; t ) -1 . The matrix Θ T ∗ H R 0 ( z )[ Θ ( ∗ ; t ) , Y T GZ ] has constant dimension and concentrates around a deterministic limit as n, d →∞ . First, by using Y Θ ( ∗ ; t ) = 0 and Lemma 2.2, we get $$\theta ^ { T } _ { H R _ { 0 } ( z ) \Theta ( \cdot ; t ) } = -$$ Second, using again Y Θ ( ∗ ; t ) = 0 , we get Θ T ∗ H R 0 ( z ) Y T GZ = 0. Summarizing, we proved $$\theta ^ { T } _ { H R _ { 0 } ( z ) | \theta ( x ; t ) , Y ^ { T } ( }$$ Consequently, using Lemma 5.1, the integrand in Eq. (90) converges to $$F _ { \alpha } ( z ) := A _ { \infty } ( z ) T ^ { M j } ( z ; t ) ^ { - 1 } A _ { \infty }$$ The only pole of F ∞ ( x ) inside γ arises because of the factor M ∞ j ( z ; t ) -1 at z = z ∗ . Let $$v = \left[ v _ { v } \right] \in k e r ( M ^ { 3 } g ( z )$$ Since z ∗ is a simple zero of det( M ∞ j ( z ; t )), the inverse M ∞ j ( z ; t ) -1 has a simple pole at z ∗ with residue $$\frac { v _ { 0 } T } { \sqrt { 2 } M g ( z , t ) w ^ { 3 / 2 } }$$ By the residue theorem, the contour integral equals $$\frac { 1 } { v ^ { T } \alpha ( z * ) } \int _ { v T \Delta } ^ { v T } \frac { F_{\infty}(z) d z } { A_{\infty}(z) } = - \frac { 1 } { ( z , ) ^ { 2 } }$$ where v H is the vector comprising the first r coordinates of v ≤ . The derivative ∂ z M ∞ j ( z ∗ ; t ) inherits a block structure from (84), and direct calculation yields $$\frac { v ^ { 1 } \partial _ { z } M _ { f } ( z ; t ) } { v } = - \frac { T } { v ^ { 2 } } ( E [ O _ { t } O _ { f } ] + a _ { 2 } S ( z ; t ) ) v < .$$ Substituting into (96), the ( z ∗ ) 2 factors cancel, yielding the following compact formula for eigenvector correlation: $$\rho \leq ( E [ O _ { t } O _ { t } ] + \alpha _ { 2 } S ( z _ { 2 } , t ) ) v \leq$$ By the factorization formula (86), z ∗ can be either a zero of det( M ∞ j, R ( z ; t )) or a zero of det( M ∞ j, H ( z ; t )), In the first case, we will show that v H = 0, whence the expression in Eq. (98) vanishes: such eigenvectors are asymptotically orthogonal to Θ ∗ H . On the other hand, if z ∗ is a zero of det( M ∞ j, H ( z ; t )), then v T ⩽ = [ v T H 0 T ] , and the block-diagonal structures of E [ O t O T t ] and S ′ ( z ∗ ; t ) simplify (98) to an expression involving only the hard block, yielding strictly positive correlation. This completes the eigenvector alignment analysis (full details in Appendix C.5). ## 6 Discussions and conclusions We studied feature learning in k -index models. For each t ≥ 0, we determine a threshold δ ∗ NN ( t ) in the number of samples per dimension, at which the spectrum of the Hessian of the empirical risk exhibits a qualitative change. For n/d above δ ∗ NN ( t ) the Hessian after t GD steps develops negative eigendirections which are correlated with 'hard' direction in the latent space. We expect the limit δ ∗ NN ( ∞ ) = lim t →∞ δ ∗ NN ( t ) to be the threshold for feature learning under GD. Within the conjectured scenario, feature learning under GD is essentially equivalent to feature learning in a spectral algorithm that computes the low-lying eigenvectors of a matrix of the form X T DX . The diagonal matrix D depends on the response variables y i through a preprocessing function implicitly determined by GD. Since this preprocessing is constrained by the network architecture rather than optimized for detection, the threshold δ NN is generally suboptimal compared to the best possible spectral algorithms. This scenario also provides a quantitative explanation of grokking: for δ > δ NN , the network first overfits to the training data, and subsequently learn the 'hard' features from the negative curvature of the Hessian, leading to a delayed drop in generalization error. Several problems remain open: - ( i ) Connect the spectral phase transition at δ NN ( ∞ ) with the GD dynamics. We expect that for 1 ≪ t ≲ log d , GD dynamics will be closely related to the dynamics in a saddle with the spectral properties that we characterize. Making this rigorous requires studying the dynamics beyond t = O (1). - ( ii ) Analysis for number of neurons m > 1, m = O (1). In this regime the reduction to the block diagonal form provided by Lemma 2.3 no longer holds and hence the analysis is more complicated. - ( iii ) Non-Gaussian and non-isotropic covariates. In particular, it would be important to understand the interplay between the low-dimensional structure in the function h and lowdimensional structures in the covariates distribution. ## Acknowledgments The authors thank Chen Cheng, Jason D. Lee and Basil Saeed for insightful discussions. The authors also thank Fan Nie and Jiatong Yu for the help with numerical simulations. This work was supported by the NSF through Award DMS-2031883, the Simons Foundation through Award 814639 for the Collaboration on the Theoretical Foundations of Deep Learning, and the NSF Award MFAI-2501597. ## References - [ABAM22] Emmanuel Abbe, Enric Boix-Adsera, and Theodor Misiakiewicz, The merged-staircase property: a necessary and nearly sufficient condition for sgd learning of sparse functions on two-layer neural networks , Advances in Neural Information Processing Systems, vol. 35, 2022, pp. 27814-27827. - [ABAM23] Emmanuel Abbe, Enric Boix-Adsera, and Theodor Misiakiewicz, Sgd learning on neural networks: leap complexity and saddle-to-saddle dynamics , Conference on Learning Theory, PMLR, 2023, pp. 2533-2635. - [ADK + 24] Luca Arnaboldi, Yatin Dandi, Florent Krzakala, Luca Pesce, and Ludovic Stephan, Repetita iuvant: Data repetition allows sgd to learn high-dimensional multi-index functions , arXiv:2405.15459 (2024). - [AGHJ25] Gerard Ben Arous, Reza Gheissari, Jiaoyang Huang, and Aukosh Jagannath, Local geometry of high-dimensional mixture models: Effective spectral theory and dynamical transitions , arXiv:2502.15655 (2025). - [AGJ21] Gerard Ben Arous, Reza Gheissari, and Aukosh Jagannath, Online stochastic gradient descent on non-convex losses from high-dimensional inference , Journal of Machine Learning Research 22 (2021), no. 106, 1-106. - [AMB + 19] Benjamin Aubin, Antoine Maillard, Jean Barbier, Florent Krzakala, Nicolas Macris, and Lenka Zdeborov´ a, The committee machine: computational to statistical gaps in learning a two-layers neural network , Journal of Statistical Mechanics: Theory and Experiment (2019), no. 12, 124023. - [AMS25] Kiana Asgari, Andrea Montanari, and Basil Saeed, Local minima of the empirical risk in high dimension: General theorems and convex examples , arXiv:2502.01953 (2025). - [Bac17] Francis Bach, Breaking the curse of dimensionality with convex neural networks , Journal of Machine Learning Research 18 (2017), no. 19, 1-53. - [BAP05] Jinho Baik, Gerard Ben Arous, and Sandrine Peche, Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices , The Annals of Probability 33 (2005), no. 5, 1643-1697. - [Bar96] Peter Bartlett, For valid generalization the size of the weights is more important than the size of the network , Advances in neural information processing systems 9 (1996). - [BBPV23] Alberto Bietti, Joan Bruna, and Loucas Pillaud-Vivien, On learning gaussian multiindex models with gradient flow , arXiv:2310.19793 (2023). - [BES + 22] Jimmy Ba, Murat Erdogdu, Taiji Suzuki, Zhichao Wang, Denny Wu, and Greg Yang, High-dimensional asymptotics of feature learning: How one gradient step improves the representation , Advances in Neural Information Processing Systems, vol. 35, 2022, pp. 37932-37946. - [BGN11] Florent Benaych-Georges and Raj Rao Nadakuditi, The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices , Advances in Mathematics 227 (2011), no. 1, 494-521. - [BH25] Joan Bruna and Daniel Hsu, Survey on algorithms for multi-index models , arXiv:2504.05426 (2025). - [BKM + 19] Jean Barbier, Florent Krzakala, Nicolas Macris, L´ eo Miolane, and Lenka Zdeborov´ a, Optimal errors and phase transitions in high-dimensional generalized linear models , Proceedings of the National Academy of Sciences 116 (2019), no. 12, 5451-5460. - [BMZ24] Rapha¨ el Berthier, Andrea Montanari, and Kangjie Zhou, Learning time-scales in twolayers neural networks , Foundations of Computational Mathematics (2024), 1-84. - [BS10] Zhidong Bai and Jack Silverstein, Spectral Analysis of Large Dimensional Random Matrices ( 2 nd edition) , Springer, 2010. - [BY93] Z. D. Bai and Y. Q. Yin, Limit of the Smallest Eigenvalue of a Large Dimensional Sample Covariance Matrix , The Annals of Probability 21 (1993), 1275-1294. - [CC17] Yuxin Chen and Emmanuel J Cand` es, Solving random quadratic systems of equations is nearly as easy as solving linear systems , Communications on pure and applied mathematics 70 (2017), no. 5, 822-883. - [CCM21] Michael Celentano, Chen Cheng, and Andrea Montanari, The high-dimensional asymptotics of first order methods with random data , arXiv:2112.07572 (2021). - [CM14] Benoit Collins and Camille Male, The strong asymptotic freeness of haar and deterministic matrices , The Annals of Probability 42 (2014), no. 3, 869-920. - [CM20] Sitan Chen and Raghu Meka, Learning polynomials of few relevant dimensions , Proceedings of Thirty Third Conference on Learning Theory, vol. 125, PMLR, 2020, pp. 1161-1227. - [CMW20] Michael Celentano, Andrea Montanari, and Yuchen Wu, The estimation error of general first order methods , Conference on Learning Theory, PMLR, 2020, pp. 1078-1141. - [DDM + 25] Leonardo Defilippis, Yatin Dandi, Pierre Mergny, Florent Krzakala, and Bruno Loureiro, Optimal spectral transitions in high-dimensional multi-index models , arXiv:2502.02545 (2025). - [DLB25] Alex Damian, Jason D. Lee, and Joan Bruna, The generative leap: Sharp sample complexity for efficiently learning gaussian multi-index models , arXiv:2506.05500 (2025). - [DLS22] Alex Damian, Jason D. Lee, and Mahdi Soltanolkotabi, Neural networks can learn representations with gradient descent , Conference on Learning Theory, Proceedings of Machine Learning Research, vol. 178, PMLR, 2022, pp. 5413-5452. - [DPVLB24] Alex Damian, Loucas Pillaud-Vivien, Jason D. Lee, and Joan Bruna, Computationalstatistical gaps in gaussian single-index models , arXiv:2403.05529 (2024). - [DTA + 24] Yatin Dandi, Emanuele Troiani, Luca Arnaboldi, Luca Pesce, Lenka Zdeborov´ a, and Florent Krzakala, The benefits of reusing batches for gradient descent in two-layer networks: Breaking the curse of information and leap exponents , arXiv:2402.03220 (2024). - [FW26] Zhou Fan and Leda Wang, High-dimensional learning dynamics of multi-pass stochastic gradient descent in multi-index models , 2026. - [FWNL25] Hengyu Fu, Zihao Wang, Eshaan Nichani, and Jason D. Lee, Learning hierarchical polynomials of multiple nonlinear features with three-layer networks , International Conference on Learning Representations, 2025. - [KY13] Antti Knowles and Jun Yin, The isotropic semicircle law and deformation of wigner matrices , Duke Mathematical Journal 162 (2013), no. 1, 25-78. - [KY17] Antti Knowles and Jun Yin, Anisotropic local laws for random matrices , Probability Theory and Related Fields 169 (2017), 257-352. - [KZM25] Filip Kovaˇ cevi´ c, Yihan Zhang, and Marco Mondelli, Spectral estimators for multi-index models: Precise asymptotics and optimal weak recovery , arXiv:2502.01583 (2025). - [LL20] Yue M. Lu and Gen Li, Phase transitions of spectral initialization for high-dimensional nonconvex estimation , IEEE Transactions on Information Theory 66 (2020), no. 2, 889-913. - [LM21] Zhenyu Liao and Michael W. Mahoney, Hessian eigenspectra of more realistic nonlinear models , Journal of Machine Learning Research 22 (2021), no. 163, 1-56. - [LOSW24] Jason D. Lee, Kazusato Oko, Taiji Suzuki, and Denny Wu, Neural network learns low-dimensional polynomials with sgd near the information-theoretic limit , arXiv:2406.01581 (2024). - [MM19] Marco Mondelli and Andrea Montanari, Fundamental limits of weak recovery with applications to phase retrieval , Foundations of Computational Mathematics 19 (2019), no. 3, 703-773. - [MP67] V. Marchenko and L. Pastur, Distribution of eigenvalues for some sets of random matrices , Mathematics of the USSR-Sbornik 1 (1967), no. 4, 457-483. - [MU25] Andrea Montanari and Pierfrancesco Urbani, Dynamical decoupling of generalization and overfitting in large two-layer networks , Advances in Neural Information Processing Systems, 2025. - [PBE + 22] Alethea Power, Yuri Burda, Harri Edwards, Igor Babuschkin, and Vedant Misra, Grokking: Generalization beyond overfitting on small algorithmic datasets , arXiv:2201.02177 (2022). - [RV18] Cynthia Rush and Ramji Venkataramanan, Finite sample analysis of approximate message passing algorithms , IEEE Transactions on Information Theory 64 (2018), no. 11, 7264-7286. - [SC95] Jack W Silverstein and Sang-Il Choi, Analysis of the limiting spectral distribution of large dimensional random matrices , Journal of Multivariate Analysis 54 (1995), no. 2, 295-309. - [SZ81] Haim Sompolinsky and Annette Zippelius, Dynamic theory of the spin-glass phase , Physical Review Letters 47 (1981), no. 5, 359. - [SZ82] , Relaxational dynamics of the edwards-anderson model and the mean-field theory of spin-glasses , Physical Review B 25 (1982), no. 11, 6860. - [Tao23] Terence Tao, Topics in random matrix theory , vol. 132, American Mathematical Society, 2023. - [TDD + 25] Emanuele Troiani, Yatin Dandi, Leonardo Defilippis, Lenka Zdeborov´ a, Bruno Loureiro, and Florent Krzakala, Fundamental computational limits of weak learnability in high-dimensional multi-index models , Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, vol. 258, PMLR, 2025, pp. 2467-2475. - [VDVW96] Aad W Van Der Vaart and Jon A Wellner, Weak convergence and empirical processes: with applications to statistics , Springer, 1996. - [Ver18] Roman Vershynin, High-dimensional probability: An introduction with applications in data science , vol. 47, Cambridge University Press, 2018. - [WNL24] Zihao Wang, Eshaan Nichani, and Jason D. Lee, Learning hierarchical polynomials with three-layer neural networks , International Conference on Learning Representations, 2024. - [ZWFL25] Bohan Zhang, Zihao Wang, Hengyu Fu, and Jason D. Lee, Neural networks learn generic multi-index models near information-theoretic limit , arXiv:2511.15120 (2025). ## A Technical tools We begin with some technical tools that will be used in the proof of the main results. Most of these are standard or variations of standard results. ## A.1 Rouch´ e's theorem for random functions Lemma A.1. Let D ⊂ C be a bounded connected open domain. For each n ≥ 1 , let f n : D → C be a random analytic function, and let f : D → C be a nonzero deterministic analytic function. Assume for every compact K ⊂ D , $$\sup _ { z \in K } | f ( z ) - f ( z ) | \rightarrow 0 .$$ Fix z ∗ ∈ D and suppose f has a zero of order k ≥ 1 at z ∗ ; that is, f ( z ) = ( z -z ∗ ) k g ( z ) near z ∗ with g analytic and g ( z ∗ ) = 0 . Then there exists r 0 > 0 such that for every 0 < r ≤ r 0 , the number of zeros N n ( r ) of f n in D ( z ∗ , r ) := { z ∈ C : | z -z ∗ | < r } (counted with multiplicity) satisfies $$P ( N _ { n } ( r ) = k ) - > 1 .$$ Furthermore, there exist random zeros z n, 1 , . . . , z n,k ∈ D ( z ∗ , r ) of f n such that $$\max _ { 1 \leq k \leq n } | z _ { n , k } - z ^ { * } _ { 1 } | \rightarrow 0 .$$ Proof. Since f is analytic and nonvanishing in D , its zeros are isolated. Choose r 0 > 0 such that D ( z ∗ , r 0 ) ⊂ D , z ∗ is the only zero of f in D ( z ∗ , r 0 ), and f has no zeros on the circle Γ r 0 := { z ∈ C : | z -z ∗ | = r 0 } . For 0 < r ≤ r 0 , define $$where F _ { r } := \{ z \in C : | z - z * | = 1 \} .$$ $$a _ { n , r } = \{ s \in F _ { p } ^ { | f _ { n } ( z ) - f ( z ) | < m _ { p } \}$$ Define the event By the uniform convergence (99), P ( A n,r ) → 1. On A n,r , we have | f n ( z ) -f ( z ) | < m r / 2 ≤ | f ( z ) | for all z ∈ Γ r . By Rouch´ e's theorem, f n and f have the same number of zeros (counted with multiplicity) in D ( z ∗ , r ). Since f has exactly k zeros in D ( z ∗ , r ) (all at z ∗ ), we conclude N n ( r ) = k on A n,r , proving (100). To prove (101), fix ε > 0 and choose 0 < ρ < r ≤ r 0 with ρ < ε . Define $$z ) < m _ { r } / 2 \} , ( 104 )$$ Since P ( A n,ρ ) → 1 and P ( A n,r ) → 1, we have P ( A n,ρ ∩ A n,r ) → 1. On A n,ρ ∩ A n,r , Rouch´ e's theorem applied to Γ ρ and Γ r shows that f n has exactly k zeros in D ( z ∗ , ρ ) and exactly k zeros in D ( z ∗ , r ), with no zeros on Γ ρ or Γ r . Therefore, all k zeros in D ( z ∗ , r ) lie in D ( z ∗ , ρ ). Denoting these zeros as z n, 1 , . . . , z n,k , we have $$\max _ { 1 \leq k } | z _ { n , k } - z ^ { * } _ { n , k } | \leq p < e _ { n , p }$$ Since P ( A n,ρ ∩ A n,r ) → 1 and ε is arbitrary, (101) follows. ̸ We also need the following straightforward consequence of uniform convergence. ̸ Lemma A.2. Let D ⊂ C be a bounded connected open set, and let f n : D → C be random analytic functions and f : D → C a deterministic analytic function satisfying (99) . Suppose f ( z ) = 0 for all z ∈ D . $$P ( inf | f _ { n } ( z ) \geq c x ) - 1 .$$ Then for every compact K ⊂ D , there exists c K > 0 such that ## A.2 From pointwise convergence to uniform convergence Lemma A.3. Let Γ ⊂ C be a compact set. Let f n : Γ → C be a sequence of random functions and f : Γ → C a deterministic function, both continuous on Γ . Suppose f n ( z ) p - - → f ( z ) for each z ∈ Γ . Assume there exist random variables L n satisfying L n = O P (1) uniformly in n such that $$\left | f _ { n } ( z ) - f _ { n } ( z ^ { \prime } ) \right | \leq L _ { n }$$ Then sup z ∈ Γ | f n ( z ) -f ( z ) | p - - → 0 and consequently ∫ Γ f n ( z ) d z p - - → ∫ Γ f ( z ) d z . Proof. This is a special case of general uniform convergence theory [VDVW96, Chapter 1.5]. The family of functions { f n : n ≥ 1 } is uniformly equicontinuous by (107), and for each z ∈ Γ, the family { f n ( z ) : n ≥ 1 } is tight because it converges in probability. Hence { f n : n ≥ 1 } is tight in sup norm by [VDVW96, Theorem 1.5.7] and converges weakly to f in sup norm by [VDVW96, Theorem 1.5.4]. Since the limit is deterministic, this is equivalent to sup z ∈ Γ | f n ( z ) -f ( z ) | p - - → 0. ## A.3 Asymptotics of minimum eigenvalues The following gives asymptotics of the minimum eigenvalue for a certain class of random matrices. Lemma A.4. Let X ∈ R n × d have i.i.d. N (0 , 1) entries and G = diag( g 1 , . . . , g n ) ∈ R n × n with uniformly bounded entries | g i | ≤ C , independent of X . Let δ n := n/d → δ ∈ (0 , ∞ ) . Assume the empirical distribution ̂ v g := 1 n ∑ n i =1 δ ·-g i converges weakly to a probability measure v g in probability. Let µ ⋆ denote the limiting spectral distribution of X T GX /n . Then, the minimal eigenvalue converges in probability to the left edge of the limiting spectrum: $$\lambda _ { \min } ( X ^ { T } G X / n ) - p _ { inf supp ( A ) }$$ We remark that this result is standard and expected; for instance, when G is positive definite, it follows from [KY17]. We provide a brief proof here just for completeness, extending to the general case where G may have negative entries. Proof of Lemma A.4. Denote Let a -:= inf supp( µ ⋆ ) denote the left edge of the limiting spectral measure. The proof proceeds by establishing matching lower and upper bounds for λ min ( S n ). $$S _ { n } = X ^ { T } G X / n \in R ^ { d \times d }$$ By standard random matrix theory, the empirical spectral distribution of S n converges in probability to a deterministic probability measure µ ⋆ on R , whose Stieltjes transform α : H → H is characterized by $$z + \frac { 1 } { a } = \delta \cdot E [ \frac { G } { \delta + a G } ] ^ { , }$$ where G ∼ v g . Sign-decomposition. To analyze the extreme eigenvalues, we decompose S n into positive and negative parts. Define $$S _ { n } = A _ { + } - A _ { - } , \quad S _ { n } = 1 _ { n } \sum _ { i : g _ { i } > 0 } g _ { i } x _ { i } ^ { T } , \quad A _ { + } := - 1 _ { n } \sum _ { i : g _ { i } < 0 } ( - g _ { i } ) x _ { i } ^ { T } ,$$ and x i ∈ R d denotes the i -th row of X . By Gaussian rotational invariance, there exist orthogonal matrices O , Q ∈ O( d ) such that O T A + O = B + and Q T A -Q = B -, with B ± diagonal. Therefore, $$S _ { n } = B _ { 4 } - U B - U ^ { \prime } T ,$$ where U is Haar-distributed on O( d ) and independent of the pair ( B + , B -). Strong convergence. Conditioning on G , we apply the anisotropic local law from [KY17] to the matrices A ± . Therefore, for any polynomial p , $$\int \limits _ { 0 } ^ { 1 } \int \limits _ { 0 } ^ { 1 } T _ { p } ( p | A + ) d t d p = \int \limits _ { 0 } ^ { 1 } p _ { 0 }$$ where µ ± are the limiting spectral distributions of A ± . Since B ± are diagonal matrices with the same spectra as A ± , the same strong convergence holds for B ± . By [CM14], since U is Haar-distributed, the pair ( B + , UB -U T ) converges strongly to a freely independent pair with marginal distributions µ + and µ -. Consequently, the empirical spectral distribution of S n converges in probability to the free additive convolution µ + ⊞ ( -µ -) = µ ⋆ . Lower bound. Suppose that along some subsequence we have λ min ( S n ) → a --ε for some fixed ε > 0 for contradiction. Set t := a --ε and consider the resolvent function f t ( x ) := 1 / ( x -t ). Since t < a -= inf supp( µ ⋆ ), the function f t is uniformly bounded on supp( µ ⋆ ). By polynomial approximation (via the Weierstrass theorem) on a compact set K ⊃ supp( µ ⋆ ) such that t / ∈ K , there exist polynomials p m ( x ) satisfying sup x ∈ K | p m ( x ) -f t ( x ) |→ 0. (113) then implies However, if λ min ( S n ) → t = a --ε , then ∥ ( S n -t I ) -1 ∥ op → ∞ , contradicting (114). Hence lim inf n →∞ λ min ( S n ) ≥ a -in probability. Upper bound. Fix ε > 0 and set t ′ := a -+ ε . Choose a smooth cutoff function ϕ ∈ C ∞ ( R ; [0 , 1]) satisfying ϕ ≡ 1 on ( -∞ , t ′ ] and supp( ϕ ) ⊂ ( -∞ , t ′ + ε/ 2]. By polynomial approximation and strong convergence, $$\int \limits _ { 0 } ^ { + \infty } \phi d u \geq \mu ( - \infty , a - t ) > 0 .$$ Since ϕ ( S n ) has nonnegative eigenvalues and Tr( ϕ ( S n )) > 0 with high probability for all sufficiently large n , at least one eigenvalue of S n must lie in ( -∞ , t ′ ]. Hence lim sup n →∞ λ min ( S n ) ⩽ a -in probability. Combining the lower and upper bounds yields (108). ## A.4 Concentration inequalities First we recall the following concentration bound for quadratic forms of Gaussian vectors. Lemma A.5 (Hanson-Wright) . Let x ∼ N ( 0 , I d ) and M ∈ R d × d be deterministic. Then for any t > 0 we have $$P ( x ^ { 1 } M _ { x } - T _ { r } ( M ) | p )$$ for a universal constant c > 0 . Second, we restate the following result from Lemma 21 of [AMS25] for reference convenience. Consider X ∈ R n × d with i.i.d. N (0 , 1) entries and G = diag( g 1 , . . . , g n ) ∈ R n × n with uniformly bounded entries | g i | ≤ C , sampled independently of X . Let δ n := n/d → δ ∈ (0 , ∞ ). Assume g i is i.i.d. with distribution v g . Let µ ⋆ denote the limiting spectral distribution of X T GX /n and α ⋆ denote the Stieltjes transform of µ ⋆ . Define the error quantity $$E _ { n } = c ^ { 1 } + ( l z ^ { 2 } + \frac { n } { n } ) ,$$ where C ′ depends only on C and L is a tunable parameter. Assume δ satisfies $$\frac { 1 } { 2 ^ { n } }$$ Then, with probability at least 1 -d -cL , we have $$| q | \leq C ^ { n } _ { n } 1 + ( | z | ^ { 2 } + n ^ { 2 } n ) ^ { 2 } E _ { n , 3 }$$ where C ′′ > 0 depends only on C . ## B Proofs for Sections 2.3 and 2.4 Throughout this section, we assume without loss of generality U ∗ H = [ e 1 . . . e r ] since the input data distribution and learning algorithm are rotationally invariant. ## B.1 Proofs of Lemma 2.1 and Lemma 2.2: Discrete time DMFT The DMFT equations (26) and (28) follow immediately from [CCM21] (see, e.g. Lemma 6.1). In order to prove Eqs. (27) and (29), we rewrite the DMFT equations (22) as $$V ( t ) = - \frac { 1 } { 6 } \sum _ { s = 0 } ^ { t - 1 } R _ { 0 } ( t , s ) E$$ $$\sum _ { s = 0 } ^ { t } \theta ( t + 1 ) = \theta ( t ) - n \sum _ { s = 0 } ^ { t } R _ { t } ( t * \theta _ { s } )$$ where we leverage the fact that the dependence of F ( Θ T x i , Θ T ∗ x i , ε i ) on Θ T ∗ x i , ε i is through y i = h ( Θ T ∗ x i , ε i ), namely F ( Θ T x i , Θ T ∗ x i , ε i ) = ̂ F ( Θ T x i , h ( Θ T ∗ x i , ε i )) for some ̂ F . We also write W ∗ H = W ∗ U H , Θ ∗ H = Θ ∗ U H and similarly for W ∗ E = W ∗ U E , Θ ∗ E = Θ ∗ U E . We are going to prove the following claims via induction over t : C1 ( t ): E [ ψ ( W (0) , . . . , W ( t -1) , ε, h ( W ∗ , ε ) , W ∗ E ) W ∗ H ] = 0 for every bounded measurable ψ . $$C 2 ( 1 ) : R _ { 1 } ( s , y ) ^ { \prime } = 0 , y _ { 1 } < 1 .$$ $$C 3 ( t : ( - 1 ) , \ldots , ( - 1 ) ^ { 1 0 } ) _ { n + 1 }$$ $$C H O : C _ { 6 } H _ { 5 } \times 7 H = 0 . 8 s < t$$ Note that C1 ( t +1) implies that $$e , h ( W _ { * } , e ) , W _ { * } E ) W _ { * } H ] = 0$$ for every bounded measurable ψ which is equivalent to Eq. (27) (general ψ can be obtained by taking monotone limits), and C3 ( t ) is equivalent to Eq. (29) (by Gaussianity of (Θ ∗ , (Θ( t ) : t ≥ 0))). The claims are obviously true for t = 0 (from the definition of U ∗ H ). We assume that it is true up to time t and prove it up to time t +1. C1 ( t +1). Note that from the induction hypothesis C4 ( t ), ( W (0) , . . . , W ( t )) is independent of W ∗ H , since W ( · ) and W ∗ H are jointly Gaussian with E [ W ( t ) W T ∗ ] = C θ ( t, ∗ ). Hence, with the shorthand W ≤ t := ( W (0) , . . . , W ( t )), we can write W ( t ) = g t W ( W ≤ t -1 , W ∗ E , ξ ) for some standard Gaussian vector ξ independent of everything else. Therefore $$= E [ \phi ( W ^{s}_{t-1} , e , h ( W *, e ) , W * E ) W * H ] = E [ \phi ( W ^{s}_{t-1} , W _ { s } E , ε , h ( W *, e ) , W * E , g t ( W ^{s}_{t-1} , W _ { s } E , ε , h ( W *, e ) , W * E , g t ] = 0$$ thus proving the induction claim. C2 ( t +1). By Stein's lemma where ( a ) follows from Eq. (23) and ( b ) from C4 ( t ). We thus get $$( t ) , h ( W _ { s } , e ) \cap W ^ { T } H = 0 ,$$ where the last step follows from C1 ( t +1). C3 ( t +1). Since (Θ( t ) : t ≥ 0) , Θ ∗ are jointly Gaussian, it is sufficient to prove E [ Θ( t +1)Θ T ∗ H ] = 0. Multiplying Eq. (121) with Θ T ∗ H and taking expectation yields $$\sum _ { s = 0 } ^ { t } R _ { e } ( t , s ) E [ \Theta ( s ) \Theta _ { H } ] - n \sum _ { ( t , * ) U _ { * } H U _ { * } T } ( b ) 0 ,$$ where ( a ) uses induction hypothesis C3 ( t ) and Θ ∗ ∼ N (0 , I k ) and ( b ) uses C2 ( t +1). C4 ( t +1). This follows immediately from C3 ( t +1). ## B.2 Proof of Lemma 2.3 Define the cross term $$H _ { cross } ( \Theta ) := - \sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { m } H _ { ij } ( x _ { i } ) \nabla f _ { e } ( x _ { i } ) \nabla f _ { e } ( x _ { j } ) .$$ For the lower bound, we recall that ℓ ′′ ≥ 0, so H cross ( Θ ) is positive semi-definite. Since m ∇ 2 Risk ( Θ ) = H diag ( Θ ) + H cross ( Θ ), we obtain For the upper bound, let H ∗ ( j, j ) ∈ R d × d denote the ( j, j ) block of H ∗ for ∗ ∈ { cross , diag } . We have $$\lambda _ { \min } ( n \downarrow ^ { 2 } Risk ( \Theta ) ) \geq \lambda _ { \min } ( H _ { diag } )$$ $$\begin{align} \lambda_{min}(m \nabla ^2 Risk(\Theta)) &= m \left[ u_i \right]_{j=1}^n \leq m_j \leq m_i = \lambda_m \\ By Assumptions 1, the coefficients l^*, a_i, σ^* are bounded. so \\ \end{align}$$ By Assumptions 1, the coefficients ℓ ′′ , a j , σ ′ are bounded, so $$\sum _ { i = 1 } ^ { n } \left[ \frac { 1 } { m } \sum _ { j = 1 } ^ { n } x_{i,j} f e ( x_{i}) a _{j}^{\sigma}( \theta ) x _ { i + b _ { j } } ^{2} x _ { i , j } m \right]$$ Combining (129) and (130) yields $$\lambda _ { \min } ( m \nabla ^ { 2 } Risk ( \Theta ) ) =$$ $$\begin{array}{ll} \xi ^ { 2 } Risk(\Theta) & = C \left[ 1 - \sum _ { n = 1 } ^ { m } x_{i,n}^T \right] \cdot (using (131)) \\ & = C \left[ 1 - \sum _ { n = 1 } ^ { m } x_{i,n}^T \right] \cdot (using (132)) \\ \end{array} Inequalities (128), (131), and (132) establish the lemma since \left| \frac{1}{n} \sum _ { i = 1 } ^ { n } x_{i,n}^T \right| \leq \lambda_{min}(H_d).$$ Inequalities (128), (131), and (132) establish the lemma since ∥ ∥ 1 n ∑ n i =1 x i x T i ∥ ∥ op = C (1+ o P (1)). ## B.3 Proof of Lemma 2.4: Bulk of the Hessian Recall that the gradient descent dynamics is given by Eq. (20). We define $$( 1 3 3 )$$ and $$F ( 0 ; t ) = - \frac { 1 } { \sqrt { n } } [ F ( X _ { 0 } ( t - 1 ) , X _ { 0 } ( t ) ]$$ where F = [ F 1 . . . F m ] . Note that we add the 1 √ n scaling in (134) for scaling consistency. Step I. Gaussian conditioning. We decompose X into orthogonal components via Gaussian conditioning. Define the σ -algebra $$F _ { t } = o ( e , \theta ( 0 ) , \ldots , F ( X \theta ( 0 ) , X$$ For two random variables X , X ′ we write X ∣ ∣ F t d = X ′ if E [ φ ( X ) Z ] = E [ φ ( X ′ ) Z ] for any bounded measurable function φ and any bounded random variable Z measurable on F t . Lemma B.1 (Gaussian conditioning) . We have $$X _ { | r _ { t } | } = P ^ { d } F ( 0 ; t ) X _ { n e u } ^ { d } P ^ { - 1 } \theta ( s ; t ) + X _ { P }$$ where P M := M ( M T M ) -1 M T denotes the orthogonal projection onto span( M ) , and X new ∈ R n × d is independent of X and ε , with entries drawn i.i.d. from N (0 , 1) . Proof. We first establish an algebraic decomposition of X : $$X = P ^ { 1 } F ( 0 ; t ) X P ^ { 1 } E ( t ; t ) + P ^ { 1 } F ( 0 ; t ) = P ^ { 1 } F ( 0 ; t ) X I P ^ { 1 } E ( t ; t ) + X I$$ The first equality is purely algebraic. The second equality follows from the fact that X T F (0; t ) has columns in span( Θ ( ∗ ; t )), This is in turn a consequence of the fact that by the gradient descent update (20), for each i ≤ t -1 we have It remains to show that, conditioned on F t , we may replace P ⊥ F (0; t ) XP ⊥ Θ ( ∗ ; t ) by P ⊥ F (0; t ) X new P ⊥ Θ ( ∗ ; t ) . This follows immediately from the fact that F t is generates by X Θ ( ∗ ; t ), X T F (0; t ), and P ⊥ F (0; t ) XP ⊥ Θ ( ∗ ; t ) is uncorrelated and hence independent of these quantities. $$\therefore x ^ { \prime } F ( X _ { 0 } ( i ) , X _ { 0 }$$ ## Step II. Low-rank decomposition. Define $$\begin{aligned} Y ( t ) := P _ { F } ^ { 1 } ( t ) \cdot X _ { new } ^ { P } _ { O } ( t ) \cdot e ^ { - i \theta t } \\ &= R ^ { n \times d } _ { t } \cdot Z ( t ) := X _ { E } ^ { 0 } ( t ) ( \theta * t ) \cdot T _ { O } ( s ) \cdot e ^ { - i \theta t } \end{aligned}$$ Using Lemma B.1 in the expression for the j -th block of the Hessian, cf. Eq. (76) we get $$The second term has rank 2 ( k + m )$$ The second term has rank at most 2( k + mt ). Step III. Removing low-rank perturbations via Cauchy interlacing. By Cauchy interlacing theorem, if A , B are symmetric matrices with ordered eigenvalues λ i ( A ), λ i ( B ), 1 ≤ i ≤ s and rank( B -A ) ≤ s , then λ i ( A ) ≤ λ i ( B ) ≤ λ i + s ( A ) for all 1 ≤ i ≤ d -s . As a consequence the limiting spectral distributions of the two matrices coincide. We state this as the following lemma (see e.g. [Tao23, Section 2.4.1] or [BS10, Theorem A.44]). Lemma B.2 (Interlacing for finite-rank perturbations) . Let H d be a sequence of d × d real symmetric random matrices, and denote their empirical spectral distributions by $$\mu H _ { d } = - \sum _ { i = 1 } ^ { d } \theta _ { i } - \lambda ( H _ { d } )$$ Suppose µ H d ⇒ µ in probability, where µ is a deterministic probability measure. Let W d be a sequence of real symmetric matrices with rank( W d ) ≤ s = O (1) . Define the perturbed matrices $$H _ { d } = H _ { d } + W _ { d }$$ and their empirical distributions $$\mu _ { H } ^ { d } = - \sum _ { i = 1 } ^ { d } \theta _ { i } ( H _ { i } )$$ Then we have $$\mu _ { H _ { a } } \Rightarrow \mu in probability as d \rightarrow$$ Step IV. Applying Marchenko-Pastur theory. By Lemma B.2 applied to the decomposition (140), the limiting spectral distribution LSD of H j ( t ) equals the LSD of $$H _ { j } ( t ) := Y ( t ) G _ { j } ( t ) Y ( t ) ,$$ since the rankO (1) perturbation does not affect the limit. We now further simplify ˜ H j ( t ) by removing the projection operators. Recall from (139) that Y ( t ) = P ⊥ F (0; t ) X new P ⊥ Θ ( ∗ ; t ) . Since F (0; t ) ∈ R n × mt and Θ ( ∗ ; t ) ∈ R d × ( k + mt ) have O (1) columns (by (133) and (134)), the projections P ⊥ F (0; t ) and P ⊥ Θ ( ∗ ; t ) differ from the identity by rankO (1) matrices. Therefore, Y ( t ) T G j ( t ) Y ( t ) differs from X T new G j ( t ) X new by a rankO (1) perturbation. Applying Lemma B.2 again, the LSD of ˜ H j ( t ) equals the LSD of $$H _ { j } ( t ) := X ^ { T } _ { new } G _ { j } ( t ) X _ { new } ,$$ where X new ∈ R n × d has i.i.d. N (0 , 1) entries. We now characterize the LSD of ̂ H j ( t ). Note that G j ( t ) is diagonal and independent of X new . By Lemma 2.1, the empirical spectral distribution of G j ( t ) converges in probability to the law of $$( 147 )$$ in distribution and in Wasserstein-2 distance. This falls into the framework of generalized Marchenko-Pastur laws for random matrices of the form X T GX , where X has i.i.d. Gaussian entries and G is independent of X and diagonal with a limiting spectral measure. The classical Marchenko-Pastur result states that the LSD µ j ∞ ( t ) of ̂ H j ( t ) is uniquely characterized by its Stieltjes transform α j t ( z ), which satisfies the fixed-point equation $$z + \frac { 1 } { a _ { f } ( z ) } = \delta ^ { i } \cdot E [ \frac { G _ { i } ^ { j } } { \sigma + C _ { i } a _ { f } ( z ) } ]$$ This completes the proof of Lemma 2.4. ## C Proofs of Theorems 1 and 2 This section proves Theorems 1 and 2, establishing the existence of a sharp phase transition for detecting hard directions via Hessian outliers. We refer to Section 5 for a roadmap. Throughout this appendix, as in the previous one, we assume without loss of generality U ∗ H = [ e 1 . . . e r ] . ## C.1 Outliers via resolvent analysis We fix j ∈ [ m ] and time t ≥ 0 throughout. Recall that $$H _ { j } ( t ) = X ^ { T } G _ { j } ( t ) X , G$$ where g is defined by Eqs. (77), (78). We also recall the augmented parameter and gradient matrices defined in Eq. (79) $$\theta ( t ) = [ \theta _ { 1 } , \theta ( 0 ) ]$$ $$F ( 0 ; t ) = \frac { 1 } { \sqrt { n } } [ F ( X _ { 0 } ^ { * } e ) . . . F ( X _ { 0 } ^ { * } ( t - 1 ) , X _ { 0 } ^ { * } ( 0 ) , X _ { 0 } ^ { * } e ) \cdots$$ In what follows, j ∈ [ m ], t ≥ 0 will be fixed and hence we drop these in Θ ( ∗ ; t ) , F (0; t ) , G j ( t ) , H j ( t ) whenever clear, writing Θ ∗ ; t , F 0; t , G , H . By Lemma B.1 in Appendix B, conditioned on the σ -algebra F t , the data matrix admits the representation X ∣ ∣ F t d = P ⊥ F 0; t X new P ⊥ Θ ∗ ; t + XP Θ ∗ ; t . We also recall the notations $$\begin{aligned} Y := P _ { 0 } ^ { l } X _ { new } P _ { 0 } ^ { l } e ^ { - R \times d } , Z := X _ { 0 } ^ { l } ( \sigma ^ { T } _ { t } \sigma _ { t } ) ^ { - 1 } e ^ { R } \end{aligned}$$ Let z ∈ C . We seek a formula for the resolvent ( H -z I ) -1 that isolates outlier contributions. Applying the Woodbury matrix identity to the resolvent, we obtain the representation (82), which we copy here $$( H - z I ) ^ { - 1 } = R _ { 0 } ( z ) - \left[ R _ { 0 } ( z ) , ( O _ { e + t } , Y ^ { T } G Z ) M ( z ) - 1 \right] \theta ^ { T } _ { * } R _ { 0 } ( z )$$ where R 0 ( z ) := ( Y T GY -z I ) -1 and the outlier matrix is defined by By Lemma 2.4 and Lemma A.4, for z < 0 in the left of the bulk of µ ∞ ( t ), the matrix Y T GY -z I is invertible with probability tending to one as n, d →∞ . From (153), z is an outlier eigenvalue of H if and only if det( M ( z )) = 0. Thus outlier eigenvalues satisfy $$det ( M ( z ) ) = 0 , for z < c ( t ).$$ where we recall that c ( t ) denotes the left edge of µ ∞ ( t ). Let D := { z ∈ C : Re( z ) < min( c ( t ) , 0) , | Im( z ) | ≤ 1 , | z | ≤ C } for C a sufficiently large constant. Further, let V ∗ ⩽ ∈ R r denotes the first r coordinates of V ∗ . We now analyze the asymptotic behavior of M ( z ) as n, d →∞ , which is the main result of this section. Lemma C.1. Under the assumptions of Theorem 1, for z ∈ D , we have $$M _ { j } ( z ; t ) \rightarrow M ^ { j } _ { f } ( z ; t ) ,$$ where the limiting outlier matrix is given by $$\begin{array}{ll} M ^ { j } ( z ; t ) = \left\{ \begin{array}{l} - 1_{z} E [ O , O ^ { T }] \\ I & \text{if } z \leq V * T \end{array} \right. \end{array}$$ Here S j R ( z ; t ) ∈ R ( k -r + m + mt ) × ( k -r + m + mt ) is a deterministic matrix obtained by taking expectations of certain matrix-valued functions of the discrete-time DMFT random variables (up to time t ). The proof of Lemma C.1 occupies the remainder of this section. We proceed by analyzing each block of M ( z ) separately. Step I. Simple blocks of outlier matrix. We now analyze the block structure of M ( z ) and establish concentration for its components. The first term in (154) has explicit inverse: $$\left | z ^ { i } G Z _ { 1 } \right | = \left | - I _ { 1 } - z ^ { i } G Z _ { 2 } \right |$$ For later use, we introduce concatenated vectors $$v _ { t } = \begin{bmatrix} v ( 0 ) \\ \vdots \\ v ( t ) \end{bmatrix},$$ Define the scalar random variable and consequently $$\left | \begin {matrix} z ^ { T } G Z _ { I } ^ { - 1 } \\ I ^ { - 1 } _ { 0 } \end {matrix} \right ]$$ $$(160)$$ where ( V ( t ) , V ∗ , ε ) are the state evolution random variables from Lemma 2.1 and h ( V ∗ , ε ) is the output function. In what follows we drop the subscript t and write G when the context is clear. Lemma C.2. Under the assumptions of Theorem 1, we have $$z ^ { T } G Z \rightarrow ( E [ O _ { t } O _ { t } ] , 1 ) ^ { - 1 } E [ G V _ { t + 1 } V _ { t + 1 } ] ( E [ O _ { t } O _ { t } ] ) ^ { - 1 } , 1 )$$ Proof. Expand the quadratic form: $$z ^ { T } G Z = ( \epsilon ^ { T } _ { i, t } o _ { i, t } ) ^ { - 1 } = ( \epsilon ^ { T } _ { i, t } o _ { i, t } ) ^ { - 1 } ($$ By Assumption 1, the function g ( · ) is uniformly bounded (by ∥ ℓ ′ ∥ ∞ ∥ σ ′′ ∥ ∞ < ∞ ) and Lipschitz continuous. Therefore, the test function $$\phi ( \theta ( 0 ) ^ { T } x _ { i } , \theta ( 1 ) ^ { T } x _ { j } , \ldots )$$ has at most quadratic growth and is locally Lipschitz. Applying Lemma 2.1 yields $$\sum _ { i = 1 } ^ { n } g ( \theta ( t ) , x _ { i } , y _ { i } )$$ where G is defined in (160). Similarly, Lemma 2.2 gives Θ T ∗ ; t Θ ∗ ; t p - - → E [ O t O T t ], which is positive definite under Assumption 2. By continuous mapping, $$\frac { ( \theta T _ { + } , \theta _ { + } ) ^ { - 1 } p _ { + } } { ( E [ O _ { t } , O _ { t } ] ) ^ { - 1 } } \rightarrow ( E [ O _ { t } , O _ { t } ] ) ^ { - 1 } .$$ Combining (163), (165), and (166) yields (161). Equation (162) follows by substituting into (158). Step II. Resolvent-dependent terms. We now analyze the second term in (154), which involves the resolvent R 0 ( z ) = ( Y T GY -z I ) -1 . This term has a 2 × 2 block structure: $$= \left[ z ^ { T } G Y R _ { 0 } ( z ) \right]$$ We analyze each block separately. Recall from (152) that Y = P ⊥ F 0; t X new P ⊥ Θ ∗ ; t . Hence $$( 1 6 8 )$$ Perform a spectral decomposition on the central sandwich: since P ⊥ F 0; t GP ⊥ F 0; t is symmetric, write $$P _ { F 0 , t } G P _ { F 0 , t } = O G O ^ { T } ,$$ where O ∈ R n × n is orthogonal and ̂ G = diag( ̂ g 1 , . . . , ̂ g n ) is diagonal. Similarly, write P ⊥ Θ ∗ ; t = UU T for some U ∈ R d × ( d -k -m -mt ) with U T U = I . Here U can be chosen as any orthonormal basis for the orthogonal complement of span( Θ ∗ ; t ). Define the rotated data matrix $$\overline { x } = 0 ^ { \prime } x _ { new } \in R ^ { n } \times d .$$ Since O is orthogonal and X new has i.i.d. N (0 , 1) entries, ˜ X also has i.i.d. N (0 , 1) entries (rotationinvariance of Gaussian). Substituting into (168), $$\begin{aligned} X ^ { \prime } G X U T - z I & = : U S U \\ - z I , & (by (168), (169), T _ { - z I } ) \\ ( 170 ) & (171) \end{aligned}$$ where Therefore $$R _ { 0 } ( z ) = ( Y ^ { T } G Y - z I )$$ This decomposition is valid for z = 0 and z not in the spectrum of S . Using (173), we compute the (1 , 1)-block of (167). Since Θ T ∗ ; t U = 0 by construction (the columns of U span the orthogonal complement of Θ ∗ ; t ), we have $$\theta _ { t } R _ { 0 } ( z ) \theta _ { x , t } = - \frac { 1 } { z } \theta _ { t } R _ { 0 } ( * t .$$ By Lemma 2.2, Θ T ∗ ; t Θ ∗ ; t p - - → E [ O t O T t ], so $$\theta _ { 1 } , R _ { 0 } ( z ) \theta _ { x _ { 1 } } - p _ { z } ^ { - 1 } [ E [ O _ { c } , T ] ]$$ For the cross-blocks (off-diagonal), note that $$\begin{aligned} \theta ^ { T } _ { R ( z ) } Y ^ { T } G Z &= \theta ^ { T } _ { t } R ( z ) Y ^ { T } G Z = - 1 _ { z } = 0 . \\ ( 1 ) P e _ { z } + U ( S - z I ) ^ { - 1 } U ^ { T } _ { 1 } y ^ { T } G Z + \theta ^ { T } _ { t } U ( S - z I ) ^ { - 1 } U ^ { T } _ { 1 } P e _ { z } Y ^ { T } G Z + \theta ^ { T } _ { t } U ( S - z I ) ^ { - 1 } U ^ { T } _ { 1 } y ^ { T } G Z = 0 . \end{aligned}$$ In the last step we used the fact that Y = P ⊥ F 0; t X new P ⊥ Θ ∗ ; t implies Y T P Θ ∗ ; t = 0 , and that Θ T ∗ ; t U = 0 . The (2 , 2)-block B lower := Z T GYR 0 ( z ) Y T GZ requires more care. We characterize the asymptotics of B lower as follows. Lemma C.3. Under the assumptions of Theorem 1, for z ∈ D , we have $$\begin{aligned} B _ { \infty } ^ { - p } B _ { 0 } ^ { p } ( z ; t ) , & B _ { 0 } ^ { \infty } ( z ; t ) = ( E [ O _ { t } O _ { t } ] ) ^ { - 1 } \left[ \frac { a ^ { i } ( z ) G ^ { j } _ { 2 } } { s } + G ^ { j } _ { 4 } d ^ { i } ( z ) V _ { * } < V _ { T } \right] \cdot 0 ^ { - 1 } \left[ \sum _ { R ( z ; t ) } S _ { R } ( z ; t ) \right] ^ { - 1 } (177)$$ Here ˜ S j R ( z ; t ) ∈ R ( k -r + m + mt ) × ( k -r + m + mt ) is some deterministic matrix depending on expectations over the discrete-time DMFT random variables. The proof of Lemma C.3 is carried out in Sections C.2 and C.3. Given Lemma C.3, we can now prove Lemma C.1. ̸ $$S = U ^ { \prime } X ^ { \prime } G X U$$ Proof of Lemma C.1. Recall from (154) that the outlier matrix M ( z ) is the sum of two terms. We analyze each term separately. First term. By Lemma C.2, the first term converges to $$\vert \overrightarrow { r } \vert ^ { 2 } = \vert - [ E ( O _ { t } , O T ) ] - E [ G ^ { i } _ { t } V _ { t } \cdot V _ { t + 1 } ] ( E ( O _ { t } ) - ( E [ O _ { t } ] ) ^ { 2 } ]$$ Second term. The second term in (154) has a 2 × 2 block structure given by (167). By (175), the (1 , 1)-block converges to -1 z E [ O t O T t ]. By (176), the cross-blocks equal to zero. By Lemma C.3, the (2 , 2)-block B lower converges to B ∞ lower ( z ; t ). Adding the two terms, we obtain $$- \left[ - 1 ^ { 2 } E [ O , O _ { t } ] + I ^ { 2 } B _ { oo } ( z ; t ) \right] \cdot ( 1 7 9 )$$ To simplify the (2 , 2)-block, we use the block decomposition V ∗ : t = [ V T ∗ ⩽ V T ∗ > : t ] T . By the second part of Lemma 2.1 (cf. (27)), the hard directions V ∗ ⩽ = P U ∗ H V ∗ has no correlation with others: $$E \left[ G ^ { i } V _ { * } \leq v T _ { * } \right ] = 0 .$$ Therefore, E [ G j t V ∗ : t V T ∗ : t ] is block diagonal: $$\sum _ { i = 1 } ^ { n } [ G ^ { i } V _ { \ast } \leq V _ { \ast } ] = \sum _ { i = 1 } ^ { n } E [ G ^ { i } V _ { \ast } > V _ { \ast } , t ] = \left[ E [ G ^ { i } V _ { \ast } < V _ { \ast } ] \right] .$$ Substituting the explicit form of B ∞ lower ( z ; t ) from Eq. (177) also using the block diagonal structure above, the (2 , 2)-block simplifies to where S j R ( z ; t ) := -E [ G j t V ∗ > : t V T ∗ > : t ] + ˜ S j R ( z ; t ). This gives the limiting outlier matrix (157). ## C.2 Resolvent concentration via leave-one-out B lower = Z T GYR 0 ( z ) Y T GZ can be written as $$\begin{aligned} B_{lower} &= z^\dagger GY ( - 1 _ { z } ) = z^\dagger GYU ( S - x _ { new } ) = z^\dagger GP _ { 0 , r } X _ { new } \\ &= \left( Y ^{z _ { 1 } } + U ( S - z I ) - 1 ^{u T } Y ^{T} GZ \right) (by eq. + U ( S - z I ) - 1 ^{u T } X ^{I}_{new} P _ { 0 , r } GZ . \end{aligned}$$ Define the Gaussian matrix $$R _ { 1 } = \overline { X U } = O ^ { T } X _ { new } U \in R ^ { n \times ( d - k ) }$$ Since U , O are measurable on F t , the entries of R are i.i.d. N (0 , 1) and independent of ( X , ε ). The Gaussian matrix R should not be confused with the resolvent R 0 ( z ) = ( Y T GY -z I ) -1 defined in Eq. (173). By Eq. (172), $$S = R ^ { \circ } G R _ { 1 }$$ $$S = U ^ { \prime } X _ { new } P _ { 0 , t } G P _ { 0 , t } X _ { new } U$$ Substituting this identity into Eq. (183) and using the fact P ⊥ F 0; t = P ⊥ F 0; t OO T , $$\rho _ { \text { lower } } = z ^ { 1 } G P _ { F _ { 0 , t } }$$ We regard this as a function of the random matrix R while conditioning on ( X , ε ). Denote the effective dimension after removing the conditioning subspace. For the Stieltjes transform α j t ( · ) of distribution µ j ∞ ( t ), when the context is clear, we drop the time subscript t and neuron index j , just denoting it by α ( · ). $$d _ { i } = d - k - m - mt .$$ The next lemma is the main result of this subsection. We establish concentration of the resolvent bilinear form a T R T ( R ̂ GR T -z I ) -1 Rb around its deterministic limit for random vectors a , b that are measureable in F t (thus independent of R ). Throughout, all statements and expectations are taken conditioning on { a , b , ̂ G } . We note that the following result is not surprising. For instance, if ̂ G is positive definite, the lemma can be derived via using its companion resolvent with the results from [KY17]. Lemma C.4. Consider a , b ∈ R n such that ∥ a ∥ , ∥ b ∥ = O P (1 / √ n ) and are measurable in F t . Then for z ∈ D , we have $$| \int _ { i = 1 } ^ { n } a _ { i b _ { i } } d _ { x } q ( z ) | - 1$$ where α ( z ) is the Stieltjes transform of distribution µ ∞ (cf. Lemma 2.4). Proof of Lemma C.4. We first state some basic facts about the spectrum of R T ̂ GR . Lemma C.5. Denote by ̂ µ n,d ( t ) the empirical spectral distribution of R T ̂ GR . Then $$\hat { \mu _ { n , t } } \rightarrow \hat { \mu _ { \infty } ( t ) } in probability.$$ where µ ∞ ( t ) is the limiting measure from Lemma 2.4. Proof. Recall from Eq. (186) that $$\begin{array}{ll} R ^ { \prime } G R = U ^ { \prime } X _ { new } ^ { T } P _ { 0 , t } G P _ { 0 , t } X _ { new } ^ { T } & R ^ { \prime } G R = U ^ { \prime } X _ { new } ^ { T } P _ { 0 , t } G P _ { 0 , t } X _ { new } ^ { T } \\ \end{array}$$ We compare this with X T new GX new , whose limiting spectral distribution is µ ∞ ( t ) by Lemma 2.4. Note that P ⊥ F 0; t GP ⊥ F 0; t -G has rank at most O (1) since F 0; t ∈ R n × mt . By Cauchy interlacing theorem (Lemma B.2), this rankO (1) perturbation does not affect the limiting spectral distribution. Similarly, the rotation by U preserves eigenvalues. Hence the ESD of R T ̂ GR converges to µ ∞ ( t ). In addition, applying Lemma A.4 to R T ̂ GR , we obtain λ min ( R T ̂ GR ) p - - → c ( t ). Step I. Regularization. For convenience, we consider an arbitrary compact subset K ⊂ D in this subsection and always assume that z ∈ K . For any such K , there exists a positive distance $$x ^ { - 1 } = \frac { 1 } { 2 z \epsilon K } ( \min ( c ( t ), 0 ) - R e ( z )$$ to the left edge. Without loss of generality, we assume that Im( z ) ≥ 0; the case Im( z ) < 0 follows by a similar reasoning. We further introduce an imaginary pertubation for technical convenience. Define the regularization parameter $$\frac { 1 } { n } = \frac { 1 } { \log n }$$ We emphasize that the concrete choice of η n does not matter as long as it converges to zero slowly enough. The following lemma shows that the regularization introduces negligible error. Lemma C.6. For any z ∈ K we have $$| \overrightarrow { R } ( R ^ { T } G R - z I ) | ^ { 2 } = 0 , (194)$$ for any sequence η n → 0 + . This lemma is proven in Section D.1. Denote the regularized resolvent $$( 1 9 5 )$$ Consequently, we work with the regularized resolvent Q in this subsection, and remove the regularization at the end. Step II. Diagonal-cross decomposition and leave-one-out. We now decompose the bilinear form into diagonal and off-diagonal parts. Let r i ∈ R d ⋆ denote the i -th row of R . Then we have $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } a _ { i , j } R _ { i } q _ { j } b = T _ { i } R _ { i } q _ { j }$$ ̸ The diagonal term S diag will be shown to concentrate around its deterministic limit, whereas the cross term S cross will be proved to vanish. To handle both terms, we employ the leave-one-out decomposition. For each i ∈ [ n ], define the leave-one-out resolvent $$q _ { i } = ( R ^ { T } G _ { i } - R - ( \frac { 1 } { q _ { e , i } } + i n ) I ) ^ { - 1 } , where G _ { i } := G - z + i n ) I ) ^ { - 1 } , where$$ and ̂ g i is the i -th diagonal entry of ̂ G . By the Sherman-Morrison formula applied to the rank-one perturbation R T ̂ GR = R T ̂ G -i R + ̂ g i r i r T i , we have $$Q = Q _ { i } - \frac { q _ { i } Q _ { i } r _ { i } T Q _ { i } } { 1 + q _ { i } r _ { i } Q _ { i } r _ { i } } .$$ Applying (198) to the bilinear form r T i Qr j , we obtain $$r _ { i } Q _ { r _ { j } } = \frac { r _ { i } Q - r _ { j } } { 1 + g _ { r _ { i } } ^ { 2 } r _ { j } ^ { 2 } } , for all i$$ Step III. Cross-term vanishing. Define for each i ∈ [ n ] ̸ $$\sum _ { j = 1 } ^ { n } h _ { i j } T Q _ { r _ { j } } , where h _ { i j } := n ^ { 2 } a _ { i j }$$ Then S cross = ∑ n i =1 T i . To bound this, we estimate the conditional second moment E [ | S cross | 2 ∣ ∣ ∣ F t ] . Lemma C.7. We have $$\int _ { \alpha } ^ { \beta } \{ f \} ^ { 2 } < \frac { 1 } { n } \sum _ { i = j } ^ { n } M _ { n } ^ { 1 } = \frac { 1 } { n } \sum _ { i = j } ^ { n } A _ { n } ^ { 2 }$$ The proof of this statement is given Section D.2. Since ∥ a ∥ , ∥ b ∥ = O P (1 / √ n ) by our assumption, we have M n = O P (1). Combining with η n = 1 / log n and d ⋆ = O ( n ), the bound (201) yields $$\sum _ { n = 1 } ^ { \infty } ( S _ { cross } | ^ { 2 } F _ { i } ) \leq \frac { 2 d _ { N } } { n ^ { 4 } } E$$ and hence by Markov's inequality S cross p - - → 0. Step IV. Diagonal-term concentration. To this end, we first prove the concentration of each diagonal term around its deterministic limit. Lemma C.8. For any z ∈ K , the following holds $$\sum _ { i \in [ n ] } | ( 1 / d _ { * } ) r _ { i } Q - i r _ { i } |$$ with probability at least 1 -d -C ⋆ . Proof. We establish uniform concentration of (1 /d ⋆ ) r T i Q -i r i around the Stieltjes transform α ( z + iη n ) via a two-step decomposition $$\begin{aligned} \{ ( 1 / d _ { x } ) r _ { i } T Q - r _ { i } - a ( z - \frac { i n _ { 0 } } { 2 } ) \} & = \{ ( 1 / d _ { x } ) r _ { i } T Q - i \} \\ \end{aligned}$$ The first term is controlled via the Hanson-Wright concentration bound (Lemma A.5). Since r i is independent of Q -i , using (116) with M = Q -i , and using the resolvent bounds ∥ Q -i ∥ op ≤ (1 /η n ) and ∥ Q -i ∥ F ≤ √ d ⋆ /η n , we set t = α n √ d ⋆ /η n with α n = C log n to obtain $$\rho ( r _ { i } T Q - r _ { i } - T r ( Q - n )$$ for some large enough universal constant C > 0. Via union bounds, we have for all i ∈ [ n ] $$\vert ( 1 / d _ { \ast } ) r ^ { T } Q - r _ { i } - ($$ with probability at least 1 -O ( n -C ). For the second term, we establish concentration of the normalized trace (1 /d ⋆ ) Tr( Q -i ) around the Stieltjes transform uniformly over all i ∈ [ n ] via invoking Lemma A.6. We specialize to our setting with η n = 1 / log n and choose L a sufficiently large constant. Since | z | ≤ C for z ∈ K , the error quantity (117) becomes $$\sum _ { n = 1 } ^ { \infty } ( \log n ) ^ { n } + \frac { 1 } { n } n n _ { m }$$ For n sufficiently large, the condition (118) is automatically satisfied due to E n = o (1). Using Lemma A.6 with probability 1 -O ( d -C ⋆ ), we obtain for all i ∈ [ n ] $$\frac { \vert ( 1 / d ) Tr ( Q _ { z } - i ) - a ( z + i ) \vert } { \sqrt { n } ^ { 2 } } \leq D _ { n } = \frac { ( \log n ) ^ { 4 } } { \sqrt { n } } = \frac { ( \log n ) ^ { 3 } } { ( \log n ) ^ { 5 } } \geq \frac { 1 + ( | z | ^ { 2 } + i n ) ^ { 2 } } { n ^ { 2 } }$$ By union bounds, this holds uniformly over i ∈ [ n ] with probability at least 1 -O ( d -C ⋆ ). Combining (206) and (208) via the triangle inequality decomposition (204), we obtain $$\begin{aligned} \left( 1 / ( d _ { n } ) \right)^{\frac { T } { q - r _ { i } } - a ( z - \frac { ( \log n ) ^ { 2 } } { \sqrt { n } } + \frac { ( \log n ) ^ { 1 9 / 2 } } { \sqrt { n } } < \left( \log n \right)^{\frac { T } { q - r _ { i } } } \end{aligned}$$ uniformly over all i ∈ [ n ] with probability at least 1 -O ( d -C ⋆ ). This establishes the stated result. We now combine the concentration bounds. Using the leave-one-out identity (199), we write $$r _ { i } Q _ { r _ { i } } = \frac { r _ { i } T Q _ { - i r _ { i } } } { 1 + g _ { r _ { i } } ^ { 2 } r _ { i } ^ { 2 } Q _ { - i r _ { i } } } .$$ Subtracting the target quantity d ⋆ α 1+ d ⋆ ̂ g i α , we obtain To control denominators, we invoke the following uniform bounds. Lemma C.9 (Denominator Bounds) . For all z ∈ K , η ∈ (0 , 1) , there exists an absolute constant C ε K > 0 such that $$\frac { 1 } { c _ { k } ^ { x } } \left | 1 + g _ { i } q ( z + i n ) r _ { i } \right | \geq \frac { 1 } { c }$$ uniformly and deterministically over all i ∈ [ n ] . The first bound follows from the gap ε K between z and the edge c ( t ). The second bound is proven in Section D.3. Substituting the bounds from Lemma C.9 and Lemma C.8 into (211), and using | α | ≤ C (since we are working in the compact subset K ), | ̂ g i | ≤ C/n , d ⋆ = O ( n ), we obtain $$\frac { \vert 1 + q _ { i } r ^ { T } Q - i r _ { i } \vert } { \vert 1 + q _ { i } r ^ { T } Q - i r _ { i } \vert } = \frac { C ^ { 4 } _ { n } C _ { k } ^ { x } D _ { n } = O P ( \left | \log _ { n } \frac { d _ { i } q _ { i } a } { \sqrt { n } } \right | )$$ Multiplying by d ⋆ and summing over i ∈ [ n ] with weights a i b i , and applying the triangle inequality, we obtain the following bound $$\begin{array}{ll} \sum _ { i = 1 } ^ { n } a _ { i } b _ { i } r _ { i } q _ { r _ { i } } - \sum _ { i = 1 } ^ { n } a _ { i } b _ { i } & \leq \sum _ { i = 1 } ^ { n } | a _ { i } b _ { i } | \cdot C ^ { 4 } _ { 4 } c _ { 8 } d _ { 8 } + D \\ \leq \left| a _ { 1 } || b _ { 1 } \right| \cdot C ^ { 4 } _ { 4 } c _ { 8 } d _ { 8 } + D \\ = O_P \left( \frac { 1 } { n } \cdot (log n)^{19/2} \sqrt{n} \right) \end{array}$$ Last, we remove the distortion iη n in the resolvent. Using Lipschitz continuity of α ( · ) in K (standard property for Stieltjes transforms), we have $$\left | a ( z + i n ) - a ( z ) \right | \leq C _ { e } x ^ { n }$$ for all z ∈ K . Repeating the denominator bounds and algebraic expansions as in (211), we obtain $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } a _ { ij } b _ { i j } d _ { i j } ( z + i n )$$ Combining (214) and (216) yields the desired diagonal-term concentration $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } a _ { i b _ { i } r _ { j } } q _ { r _ { i } } - \sum _ { i = 1 } ^ { n }$$ Combining the cross-term vanishing (202) and the diagonal-term concentration (217), and also Lemma C.6 to remove the regularization, we obtain the main result (Lemma C.4) of this subsection. ## C.3 Convergence of the outlier matrix This subsection completes the proof of Lemma C.3 by establishing the convergence of the intermediate form obtained from Lemma C.4 (Section C.2). We continue with the notation from Sections C.1 and C.2. Recall from Eq. (187) that $$B _ { \omega } = z ^ { T } G P _ { 0 , t } O B ( R ^ { T } G R - z I ) ^ { - 1 } R ^ { T } O ^ { T } P _ { 0 , t } G Z R ( R ^ { T } G R - z I ) ^ { - 1 } R ^ { T } O ^ { T } P _ { 0 , t } G Z R ( 2 1 8 )$$ and define ˜ f ( s ; z ) := d ⋆ α ( z ) 1+( d ⋆ /n ) sα ( z ) . Lemma C.4 implies that conditioning on ( X , ε ), we have $$\left\{ \begin{array}{ll} z ^ { - 1 } R ^ { T } O ^ { T } P _ { F o t } G Z - z ^ { T } G P _ { F o t } O f ( r c z I ) ^ { - 1 } & = 0. \\ g _ { z } \left( O ^ { T } P _ { F o t } G Z \right) | \frac{p}{2} & = 0. \\ \end{array} \right.$$ The main result of this subsection is the following convergence of the second term inside the absolute value in Eq. (219) to the deterministic limit B ∞ lower ( z ; t ) defined in Lemma C.3. Lemma C.10. For z ∈ D , we have $$z ^ { T } G P _ { 0 ; t } O f ( n G ; z )$$ where B ∞ lower ( z ; t ) is defined in (177) . The proof of Lemma C.10 occupies the remainder of this subsection. Combining (219) and Lemma C.10 immediately yields Lemma C.3. Proof of Lemma C.10. Throughout this subsection, for convenience, we work with z in an arbitrary fixed compact subset K ⊂ D as in Section C.2. Constants and bounds may depend on K , but are uniform over z ∈ K and n, d sufficiently large. We proceed in two steps. In Step I we reduce the intermediate form to a simpler one leveraging the specific form of ˜ f ( s ; z ). In Step II we establish convergence of the resulting expression via the discrete-time DMFT lemmas (Lemmas 2.1, 2.2). Step I. Replacement of ˜ f ( s ; z ) with f ( s ; z ) . Define $$f ( s ; z ) = \frac { 8 a ( z ) } { \delta + s a ( z ) } .$$ Recall that for a , b , we have the following representation $$\sum _ { i = 1 } ^ { n } a _ { i } f ( n g _ { i } ; z ) = a ^ { T } f ( \widehat { G } ; z )$$ To relate the left-hand side term in (220) to (222), we approximate d ⋆ by n/δ as follows. Lemma C.11. Let a , b ∈ R n with ∥ a ∥ , ∥ b ∥ = O P (1 / √ n ) . Then $$\sum _ { i = 1 } ^ { n } \frac { a _ { b _ { r } } } { 1 + d _ { q _ { j } } ( z ) }$$ We prove this technical lemma in Section D.4. With this lemma, we can work on the form d · Z T GP ⊥ F 0; t O f ( n ̂ G ; z ) O T P ⊥ F 0; t GZ instead of the more complicated one Z T GP ⊥ F 0; t O ˜ f ( n ̂ G ; z ) O T P ⊥ F 0; t GZ from now on, since the error is o P (1). Define Since f ( s ; z ) = δα ( z ) · ( δ + sα ( z )) -1 = -δ ( s z -s ) -1 , we obtain $$\frac { 1 } { n } - \frac { 1 } { n + 1 } = f ( z ) .$$ $$f ( \widehat { n } G ; z ) = - s ( \widehat { s } _ { 1 } I - \widehat { n } G ) ^ { - 1 }$$ Define the rescaled matrix ˜ G := n G = diag( g 1 , . . . , g n ), so that ˜ G has O (1) diagonal entries. Recall Z = X Θ ∗ ; t ( Θ T ∗ ; t Θ ∗ ; t ) -1 . Expanding Z and using G = ˜ G /n , we have $$= ( \theta T _ { * } \theta _ { * } t ) ^ { - 1 } d e$$ Since p-lim n,d →∞ ( Θ T ∗ ; t Θ ∗ ; t ) -1 = ( E [ O t O T t ] ) -1 , we can focus on the core term $$B _ { core } := \frac { d } { n ^ { 2 } } \Theta T _ { i , j } X ^ { i } T ^ { j } C$$ and add the ( E [ O t O T t ] ) -1 back in the end. By (169), O ̂ GO T = P ⊥ F 0; t GP ⊥ F 0; t , so $$( 228 ) .$$ Combining (228) with the identity (225), using d/n → 1 /δ , we obtain $$B _ { core } = - 1 _ { n } \phi _ { t } x ^ { T } G P _ { t }$$ Therefore it suffices to analyze the asymptotics of the right-hand side of (229). Step II. Convergence at s = s z . We analyze the expression $$\sum _ { n = 1 } ^ { \infty } ( s _ { I - P _ { F _ { 0 ; t } } } G P _ { F _ { 0 ; t } } ) ^ { - 1 } P _ { F _ { 0 ; t } } G P _ { F _ { 0 ; t } }$$ Lemma C.12. For s = s z , we have $$\left[ \begin{array}{c} G _ { 2 } ^ { s _ { z } } - V _ { T } ^ { s _ { z } } \\ 0 \\ C _ { \infty } ( s _ { z } ; z ) \end{array} \right]$$ where C ∞ ( s z ; z ) ∈ R ( k -r + m + mt ) × ( k -r + m + mt ) is a deterministic matrix depending on expectations over certain functions of the DMFT random variables. The remainder of this step is devoted to the proof of Lemma C.12. Proof of Lemma C.12. Throughout this proof, we write s = s z and define $$B = s _ { 2 } I - G .$$ By Lemma C.9, which provides a uniform lower bound on | 1+ d ⋆ ̂ g i α ( z ) | , we have ∥ B -1 ∥ op = O (1). By the Woodbury matrix identity, we have $$( s I - P F _ { 0 ; t } G P _ { 0 ; t } ) ^ { - 1 } = B ^ { - 1 } _ { 1 } - B ^ { - 1 } U ( I + K _ { p } U ^ { T } B ^ { - 1 } U ) ^ { - 1 } = B ^ { - 1 } _ { 1 } - B ^ { - 1 } U ( I + K _ { p } U ^ { T } B ^ { - 1 } U ) ^ { - 1 }$$ where $$where$$ Substituting (233) into (230), we obtain $$\begin{array}{ll} K _ { n } ( s ; z ) = K ^ { ( 1 ) } _ { n } ( s ; z ) - K ^ { ( 2 ) } _ { n } ( s ; z ) \\ &= \frac{1}{n!} \sum_{k=0}^{n-1} (-1)^k \binom{n}{k} s^k z^{n-k} \end{array}$$ where $$\begin{aligned} K^{(1)}( s ; z ) := 1 \overline{\sigma T}^t X^T G_n & K^{(2)}( s ; z ) := 1 \overline{\sigma T}^t X^T G_n \\ P_{F_0,t}^{-1} B^{-1} P_{F_0,t} GX\theta _{s,t} & P_{F_0,t}^{-1} B^{-1} U(I + K F_0 U^T B^{-1} U)^{-1} K F_0 U^T B^{-1} P_{F_0,t} \end{aligned}$$ Expand K (1) n using P ⊥ F 0; t = I -P F 0; t : $$\begin{aligned} K^{(1)}(s;z) &= - \epsilon T_{t+1} X^{T} G \\ &= - \epsilon T_{t} X^{T} G \\ &= - \epsilon T_{t-1} X^{T} G \\ &= T_{1} - T_{2} - T_{3} \end{aligned} Similarly, expand $K^{(2)}$ using $P_{F_0;t}$ on both sides:$$ Similarly, expand K (2) n using P ⊥ F 0; t = I -P F 0; t on both sides: $$\begin{equation} K^{(2)}_n(s;z)=1-\epsilon^{T}\chi^TG(1-t)=-\epsilon^{T}\chi^TG(1-t)-\epsilon^{T}\chi^TG(1-t)+\epsilon^{T}\chi^TG=:\tau_{5}-\tau_{6}-\tau_{7}. \end{equation}$$ Thus, K n ( s ; z ) = ∑ 8 ◦ =1 c ◦ T ◦ where c ◦ ∈ {-1 , +1 } are appropriate signs. The matrix K n ( s ; z ) inherits a block structure. We partition Θ ∗ ; t = [ Θ ∗ ; t, H Θ ∗ ; t, R ] where Θ ∗ ; t, H = Θ ∗ H ∈ R d × r contains the first r columns (hard directions θ ∗ 1 , . . . , θ ∗ r ) and Θ ∗ ; t, R ∈ R d × ( k -r + m + mt ) contains the remaining. Following the structure of Θ ∗ ; t , we write K n ( s ; z ) as $$K _ { n } ( s ; z ) = \left[ K _ { n } \vert H H \right] \left[ K _ { n } \vert R R \right] ^ { \cdot }$$ where we denote the hard-hard block as [ K n ] HH ∈ R r × r , the cross blocks as [ K n ] HR , [ K n ] RH , and the rest-rest block as [ K n ] RR ∈ R ( k -r + m + mt ) × ( k -r + m + mt ) . We first focus on the hard-hard block and the cross blocks. Lemma C.13. For any ◦ ∈ { 2 , 3 , . . . , 8 } , we have $$\begin{aligned} [ T _ { 0 } ] H H \rightarrow 0 _ { R } x r , & [ T _ { 0 } ] R H \rightarrow 0 _ { R } ( k - r + m + m t ) , \\ \end{aligned}$$ The proof is deferred to the end of this subsection. For the term T 1 , recall that B = s I -˜ G is diagonal with entries s -g i . Therefore, ˜ GB -1 ˜ G is also diagonal with entries g 2 i / ( s -g i ) and we have $$x _ { i } = - \theta ^ { 1 } _ { n } x ^ { T } G B ^ { 1 }$$ By Lemma 2.1, we have $$\begin{array}{c} x _ { 1 } - p E \left[ \frac { G ^ { 2 } t } { s - G _ { t } } O _ { t } O _ { t } \right] \\ = \left[ \frac { G ^ { 2 } t } { s - G _ { t } } V _ { * } V _ { T } \right] = 0 \end{array}$$ where the cross terms vanish because hard directions V ∗ ⩽ are uncorrelated with any function over the remaining coordinates V ∗ > : t := [ V T ∗ > V (0) T . . . V ( t ) T ] T (cf. Lemma 2.1). For the rest-rest block, each T ◦ do not vanish in general and converges to its own deterministic limit. By Lemmas 2.1 and 2.2, the full rest-rest block converges to a deterministic matrix: $$( s ; z ) \in R ^ { n } ( k - r + m + mt ) \times ( k - r + m + mt ) .$$ Combining the limit of T 1 in (242) with Lemma C.13, and the rest-rest block (243), we obtain the convergence stated in Lemma C.12. Concluding. By (229) and the definition (230), B core = -K n ( s z ; z )+ o P (1). Lemma C.12 implies $$- K _ { n } ( s ; z ) \rightarrow - K _ { \infty } ( s ; z ).$$ $$\frac { G _ { i } ^ { 2 } } { s _ { z } - G _ { t } } = \frac { a ( z ) G _ { i } ^ { 2 } } { \delta + G _ { t } a ( z ) } .$$ Combining Lemma C.11, (226), (244), (245), and recalling the normalization factor ( E [ O t O T t ] ) -1 on both sides, we obtain $$\sum _ { \sigma } ^ { T } G P _ { 0 , t } O f ( n G ; z ) = \sum _ { \sigma } ^ { T } ( E [ O ; σ ] ) ^ { - 1 } \left[ \frac { a ( z ) G _ { t } } { b + G _ { t } a ( z ) } V _ { * } < V _ { * } ≤ | S _ { R } ( z ; t ) | \right]$$ which is exactly B ∞ lower ( z ; t ) as defined in (177). This completes the proof of Lemma C.10. Define ˜ S R ( z ; t ) := -C ∞ ( s z ; z ). Since s z = -δ/α ( z ), we have Proof of Lemma C.13. Fix ◦ ∈ { 2 , 3 , . . . , 8 } and consider [ T ◦ ] pq where p ∈ [ r ]. Here [ T ◦ ] pq denotes the ( p, q )-th entry of T ◦ . The case q ∈ [ r ] follows analogously. By (237) and (238), the term T ◦ is a product of some matrices and contains at least one P F 0; t . Since P F 0; t projects onto span( F 0; t ), writing F 0; t = [ f 1 , . . . , f mt ] ∈ R n × mt , we have $$P _ { F 0 , t } = F _ { 0 ; t } ( F _ { 0 ; t } F _ { 0 ; t } )$$ where the sum contains mt = O (1) terms. Substituting (247) into T ◦ , the entry [ T ◦ ] pq decomposes into a sum of O (1) terms of the form $$( - \sum _ { i = 1 } ^ { n } \theta _ { p , x _ { i } } \cdot p _ { i } ) \times ( other terms )$$ where ρ i = ρ ( Θ (0) T x i , . . . , Θ ( t ) T x i , y i ) for some bounded Lipschitz function ρ ( · ) (using the explicit form of each f j , namely the formula of the gradient at each iteration). Since θ ∗ p is a hard direction, Lemma 2.1 gives $$\sum _ { i = 1 } ^ { n } \theta _ { i } p _ { i } - p _ { i } - E$$ because V ∗ p is uncorrelated with the other random variable in the product. The 'other terms' in (248) are O P (1) by Lemmas 2.1 and 2.2, hence the entire product vanishes in probability. Summing over the O (1) terms yields [ T ◦ ] pq p - - → 0. ## C.4 Eigenvalue localization via determinant zeros In Section C.1, we established that M j ( z ; t ) p - - → M ∞ j ( z ; t ) for z ∈ D . This subsection translates this convergence into a statement about convergence of Hessian eigenvalues. The main result is the following. Lemma C.14. We have the following two statements: - ( i ) Suppose det( M ∞ j ( z ; t )) = 0 has a solution z ∗ ∈ D of multiplicity k ≥ 1 . Then there exist random eigenvalues z n, 1 , . . . , z n,k of H j ( t ) such that $$\frac { \max _ { 1 < k \leq n } | z _ { n , e } - z _ { i } | - \rightarrow 0 . } { 1 < k \leq n }$$ ̸ - ( ii ) For any compact K ⊂ D , if det( M ∞ j ( z ; t )) = 0 for all z ∈ K , then with probability 1 -o (1) , H j ( t ) has no eigenvalues in K . Therefore, Hessian eigenvalues outside the support of the LSD concentrate near solutions of the determinant equation det( M ∞ j ( z ; t )) = 0. The proof of Lemma C.14 applies Rouch´ e's theorem to the determinant equation; we present the details below. Throughout this subsection, we work on compact subsets K ⊂ D . Lemma C.15. For any compact K ⊂ D and z, z ′ ∈ K , $$t ( M _ { j } ( z ; t ) ) \leq L _ { n } | z - z ' |$$ where L n = O P (1) uniformly over n . We prove this lemma in Section D.5. Next, we establish uniform convergence of the determinant. ## Lemma C.16. For any compact , $$\sup \{ det ( M _ { j } ^ { i } ( z ; t ) \} - p > 0 .$$ $$\sup \{ det ( M _ { j } ^ { z : t } ) - det ( M ^ { z } _ { j } ^ { x : t } ) \} - 1$$ Proof. By the Lipschitz bound (251) and the pointwise convergence from Lemma C.1, we apply Lemma A.3: equicontinuity combined with pointwise convergence implies uniform convergence. We now apply Rouch´ e's theorem to relate zeros of the finite-sample determinant to those of the limiting determinant. Proof of Lemma C.14. We apply Lemma A.1 and Lemma A.2 with f n ( z ) = det( M j ( z ; t )), f ( z ) = det( M ∞ j ( z ; t )), and D = D . Part ( i ) follows from Lemma A.1, using the fact that isolated eigenvalues of H j ( t ) outside the support of the limiting spectral distribution correspond to zeros of det( M j ( z ; t )) in D . Part ( ii ) follows similarly from Lemma A.2. We now analyze the limiting determinant equation using its block structure. Lemma C.17. det( M ∞ j ( z ; t )) = 0 holds if and only if one of the following two equations holds: $$\frac { \left | - z I _ { r } + E \right | } { \left | \sigma + C _ { r } ^ { 2 } \right | } = 0 ,$$ $$\det ( z \left[ E _ { R ; t } \theta _ { R ; t } \right] ) ^ { - 1 } + ( E [ \theta _ { R ; t } \theta _ { R ; t } ] ) ^ { - 1 } S d ( z ; t ) ( E )$$ In the second equation, $$\theta _ { R ; t } := [ \theta ^ { T } , \theta ( 0 ) ]$$ where Θ ∗ ,> ∈ R k -r denotes the last k -r coordinates of Θ ∗ . We note that the first equation is our master equation (39) in the main theorem. Proof of Lemma C.17. From Lemma C.1, the limiting outlier matrix has the form $$\begin{aligned} M ^ { 3 } ( z ; t ) = \left[ - I E [ O _ { t } ] \right] _ { I }$$ $$S ( z ; t ) = \left[ E ^ { - \frac { 1 } { 2 } v _ { * } ^ { 2 } + G _ { i } ^ { 2 } a _ { q } ^ { 2 } \right] .$$ We apply the Schur complement formula to Eq. (256). Since E [ O t O T t ] is positive definite, we have $$det ( M ^ { 2 } f ( z ; t ) ) = det ( - \frac { 1 } { z } [ O _ { t } o T ] ^ { - 1 } ) \cdot d$$ The first factor is nonzero since E [ O t O T t ] is invertible. $$( 2 5 4 )$$ For the second, by Lemma 2.2, E [ O t O T t ] is block-diagonal: $$\begin{aligned} E \left[ O _ { 4 } ^ { r } t \right ] = \left [ r _ { 0 } E \left | o _ { R _ { t } } \right | \right ] \\ & = ( 2 5 9 )$$ where the r × r block corresponds to hard directions and the ( k -r + m + mt ) × ( k -r + m + mt ) block corresponds to rest directions. Using (257) and (259), we expand the second as This matrix is block-diagonal. Therefore, det( M ∞ j ( z ; t )) = 0 if and only if $$\det ( - E \left[ \frac { s + G t d _ { q } ^ { i } ( z ) } { \sigma _ { q } ^ { i } } V _ { q } ^ { T } \right] ) = 0 .$$ This product equals zero if and only if one of the two equals zero, and the conclusion follows. ## C.5 Weak recovery of hard directions from eigenvectors This subsection analyzes the correlation between Hessian eigenvectors corresponding to outliers and the hard directions Θ ∗ H = [ θ ∗ 1 , . . . , θ ∗ r ] ∈ R d × r . The main result is the following. Lemma C.18. Let z ∗ ∈ D with det( M ∞ j ( z ∗ ; t )) = 0 and multiplicity β ≥ 1 . Let U = [ U T ⩽ U T > ] T ∈ R 2( k + m + mt ) × β be a basis for ker( M ∞ j ( z ∗ ; t )) , where U ⩽ ∈ R ( k + m + mt ) × β corresponds to the upper block and U > ∈ R ( k + m + mt ) × β to the lower block. Let U ⩽ , H ∈ R r × β denote the first r rows of U ⩽ . Consider the eigenvalues { z n,ℓ } β ℓ =1 of H j ( t ) converging to z ∗ (cf. Lemma C.14), with orthonormal eigenvectors { ξ n,ℓ } β ℓ =1 . Then $$\sum _ { ( t = 1 ) } ^ { \beta } \{ O _ { T } ^ { H } S _ { n , e } l \} ^ { 2 } - P _ { T r }$$ where S ′ ( z ∗ ; t ) := ∂ ∂z S ( z ; t ) ∣ ∣ z = z ∗ . By Lemma C.17, the determinant equation det( M ∞ j ( z ; t )) = 0 factors into hard and rest parts. If z ∗ arises from the rest block (i.e., from (254)), then U ⩽ , H = 0 and the correlation (262) equals zero. Instead, if z ∗ arises from the hard block (i.e., from (253)), then by the block-diagonal structure of E [ O t O T t ] and S ′ ( z ∗ ; t ), the formula (262) simplifies to $$\sum _ { e = 1 } ^ { \beta } \{ O _ { H } t _ { e , e } \} ^ { 2 } - T _ { r }$$ which is strictly positive and depends only on the hard block. Therefore, in order to have positive correlation, the corresponding eigenvalue must arise from the equation of the hard block. The proof of Lemma C.18 occupies the remainder of this subsection. Proof of Lemma C.18. Let z ∗ ∈ D be a solution of det( M ∞ j ( z ; t )) = 0 with multiplicity β ≥ 1. By Lemma C.14, there exist random eigenvalues { z n,ℓ } β ℓ =1 of H j ( t ) such that max ℓ ∈ [ β ] | z n,ℓ -z ∗ | p - - → 0. Let { ξ n,ℓ } β ℓ =1 denote the corresponding orthonormal eigenvectors. We measure the correlation via the following. Choose a contour γ centered at z ∗ with radius small enough that γ ⊂ D and γ encloses z ∗ but no other solutions of det( M ∞ j ( z ; t )) = 0 and no points in the support of the LSD. Define $$\rho ^ { n _ { d } } = - \frac { 1 } { 2 \pi i } \int _ { j } ^ { t } \phi ( H _ { j } ( t ) - z I ) ^ { n } d t$$ With probability 1 -o (1), the eigenvalues { z n,ℓ } β ℓ =1 are the only eigenvalues of H j ( t ) inside γ , hence $$p _ { n , d } = \sum _ { t = 1 } ^ { b } e _ { n , t } e _ { T , t }$$ The total squared correlation is $$\sum _ { ( t = 1 ) } ^ { \infty } | e ^ { i H \xi _ { n } a _ { n } } | ^ { 2 } = - T _ { r } ( \sigma ^ { * } H ( H _ { j } ( t ) - H ^ { p _ { n } d } O _ { s H } ) = - 1 _ { 2 m i } f _ { r } ( \sigma T _ { r } ( \sigma ^ { * } H ( H _ { j } ( t ) - H ^ { p _ { n } d } O _ { s H } ) )$$ From Section C.1, we have the following resolvent decomposition: $$( H _ { j } ( t ) - z I ) ^ { - 1 } = R _ { 0 } ( z )$$ where R 0 ( z ) := ( Y T GY -z I ) -1 is the resolvent and M j ( z ; t ) is the outlier matrix. Substituting (267) into (266), we obtain $$\sum _ { ( z = 1 ) } ^ { \beta } | e ^ { i h s _ { n , t } } | ^ { 2 } = - \frac { 1 } { 2 \pi i } + \frac { 1 } { 2 \pi i } f ( e ^ { i h R _ { 0 } ( z ) } )$$ The contour γ avoids the spectrum of Y T GY with probability 1 -o (1) due to Lemma A.4, hence R 0 ( z ) = ( Y T GY -z I ) -1 is analytic on and inside γ conditioned on the event. Therefore, $$f _ { r } ( \theta _ { h } R _ { 0 } ( z ) \theta _ { H } ) d z = 0 .$$ Thus, $$\sum _ { ( t = 1 ) } ^ { \beta } || e ^ { i T _ { H } R _ { 0 } ( z ) } | ^ { 2 } - 1 _ { 2 \pi i } f _ { r }$$ Next, we establish uniform convergence of the entire integrand in (270). ## Lemma C.19. Define $$\rho ( z ) \bigcirc \rho _ { H } ^ { T } R _ { 0 } ( z ) [$$ and the limiting matrix $$\sum _ { k = 1 } ^ { n } A _ { k } ( z ) M ^ { x } _ { j } ( z ; t ) - 1 ^ { | A _ { k } ( z ) | } e ^ { i t } A _ { k } ( z ) , \quad A _ { k } ( z ) = - \frac { 1 } { z } [ I _ { 0 } ] ^ { t }$$ Then $$\sup _ { z \in Y } | F _ { n , d } ( z ) - F _ { x 0 } ( z ) | o p \rightarrow 0 .$$ We prove this lemma in Section D.6. Next, as a consequence of Lemma C.19, $$f _ { r } F _ { n , d ( z ) } d z \rightarrow f _ { r } F _ { \infty ( z ) } d z .$$ It therefore suffices to compute the deterministic integral ∮ γ F ∞ ( z )d z by residue calculus. The only pole inside γ arises from M ∞ j ( z ; t ) -1 at z = z ∗ . By standard analytic perturbation theory, since det( M ∞ j ( z ; t )) has a zero of order β at z ∗ and U T ∂ ∂z M ∞ j ( z ∗ ; t ) U is positive definite on ker( M ∞ j ( z ; t )), the inverse M ∞ j ( z ; t ) -1 has a Laurent expansion near z ∗ with principal part $$z _ { z } ^ { M \dot { j } ( z ; t ) U } - 1 = \frac { U ( U _ { T } + R ( z ) , z - z * } ) } { z - z * }$$ where R ( z ) is analytic around z ∗ . Therefore, by the residue theorem, $$\frac { 1 } { 2 \pi i } \int _ { - \infty } ^ { + \infty } F _ { \alpha } ( z ) d z = \frac { 1 } { 2 \pi i } \int _ { - \infty } ^ { + \infty } F _ { \alpha } ( z )$$ where we abbreviate Λ := U T ∂ ∂z M ∞ j ( z ∗ ; t ) U . Note that A ∞ ( z ∗ ) extracts only the first r rows. Therefore, $$\left[ U _ { z } \right] = - \frac { 1 } { 2 } U _ { z } \leq H e R r \times B .$$ Therefore, $$\begin{aligned} T _ { r } ( A _ { \infty } ( z ^ { * } ) T _ { h } A - 1 U < h ) & = \frac { 1 } { ( z ^ { * } ) ^ { 2 } } Tr ( U ^ { T } _ { h } A - 1 U < h )$$ From the block structure of M ∞ j ( z ; t ) in (157), the kernel satisfies $$u _ { 2 } = - \frac { 1 } { 2 } [ o , o ^ { T } ] u _ { < } .$$ Differentiating M ∞ j ( z ; t ) with respect to z gives $$\frac { \partial } { \partial z } M _ { F } ^ { j } ( x , t ) = [ \frac { 1 } { 2 } E [ O _ { j } ^ { T } ]$$ where $$\frac { S ( x ; t ) } { \alpha ^ { 2 } _ { t } ( z ) } = \frac { S ( z ; t ) } { \alpha ^ { 0 } _ { t } ( z ) } = [ s _ { H } ( z ; t ) ] _ { 0 }$$ and α j ′ t ( z ) := ∂ ∂z α j t ( z ) is the derivative of the Stieltjes transform. For u ∈ ker( M ∞ j ( z ∗ ; t )) with u = [ u T ⩽ u T > ] T , we compute $$u ^ { \prime } - M x ( z ; t ) u = \frac { 1 } { ( z * ) ^ { 2 } }$$ Using (279), $$u ^ { \prime } = - M _ { z } ( z ^ { * } ; t ) u = - 1 _ { z } ( z ^ { * } )$$ Same calculation works for U . Substituting Λ = 1 ( z ∗ ) 2 U T ⩽ ( E [ O t O T t ] + S ′ ( z ∗ ; t )) U ⩽ from (283), the ( z ∗ ) 2 factors cancel: $$\frac { 1 } { 2 \pi i } \int _ { - \infty } ^ { + \infty } F _ { \alpha } ( z ) d z = T _ { r } ($$ Combining (274) and (284) yields the formula (262), completing the proof of Lemma C.18. ## C.6 Proof of Theorem 1: Finite time spectral phase transition Theorem 1 is a direct consequence of results in Sections C.4 and C.5. The statements of outlier eigenvalues follow from Lemma C.14. The weak recovery formulas in both regimes are established in Lemma C.18 and its subsequent discussion. ## C.7 Proof of Theorem 2: Large time limit This subsection establishes the results for t →∞ via leveraging the finite-time results of Theorem 1. For technical convenience, we work on arbitrary compact subsets K ⊂ D ∞ := { z ∈ C : Re( z ) < min(0 , c j ( ∞ )) } , where c j ( ∞ ) denotes the left edge of the limiting spectrum as t →∞ . Step I. Uniform convergence of determinants. For notational convenience, we define the determinant functions $$f ( z ) := \det ( - z I + E \frac { 8 G i ^ { 2 } } { δ + C _ { i } d _ { i } } )$$ $$\sum _ { i = 1 } ^ { n } \{ \frac { s G ^ { i } _ { o } } { s + G ^ { i } _ { o } d _ { o } ( z ) } V _ { * } ^ { i } T _ { * } \} - z I _ { r } + E$$ where G j t := g ( V ( t ) , h ( V ∗ , ε ); j ), G j ∞ := g ( V ( ∞ ) , h ( V ∗ , ε ); j ), and α j t ( z ), α j ∞ ( z ) are the Stieltjes transforms of the finite-time and infinite-time measures µ j ∞ ( t ), µ j ∞ ( ∞ ) with left edges c j ( t ), c j ( ∞ ). By construction, the outlier equations (39) and (47) from Theorems 1 and 2 correspond to f t ( z ) = 0 and f ∞ ( z ) = 0. Lemma C.20. For any compact K ⊂ D ∞ , we have sup z ∈K | f t ( z ) -f ∞ ( z ) |→ 0 as t →∞ . The proof is deferred to the end of this subsection. Step II. Regime δ > δ ∗ j ( ∞ ) . By definition (48), the equation f ∞ ( z ) = 0 has at least one negative solution z < c j ( ∞ ). Let z ∗ 1 , . . . , z ∗ p < c j ( ∞ ) denote the distinct negative solutions with multiplicities k 1 , . . . , k p ≥ 1. Since f ∞ is analytic (as the determinant of a matrix with analytic entries), its zeros are isolated. For each i ∈ { 1 , . . . , p } , choose a radius r i > 0 small enough that the open disk $$D _ { i } = \{ z \in C : | z - z _ { i } | < r _ { i } \}$$ satisfies D i ⊂ D ∞ , the disks { D i } p i =1 are pairwise disjoint, and f ∞ has exactly k i zeros (counted with multiplicity) in D i and no zeros on the boundary γ i := ∂ D i . By Lemma C.20, we have $$\sup \{ f ( z ) - f _ { \alpha } ( z ) \} \rightarrow 0 \text{ as } t \rightarrow \infty .$$ Define m i := min z ∈ γ i | f ∞ ( z ) | > 0. For all sufficiently large t , we have sup z ∈ γ i | f t ( z ) -f ∞ ( z ) | < m i / 2 by (288). Therefore, for all z ∈ γ i , $$\int f ( z ) - f _ { \infty } ( z ) | < m _ { i } / 2 ≤ | f _ { \infty } ( z )$$ By Rouch´ e's theorem, f t and f ∞ have the same number of zeros (counted with multiplicity) inside γ i . Since f ∞ has exactly k i zeros in D i , so does f t for all large t . Denote the zeros of f t in D i as { z ∗ ( t,q ) i } k i q =1 (counted with multiplicity). To show convergence z ∗ ( t,q ) i → z ∗ i , fix ε > 0 and choose 0 < ρ < r i with ρ < ε . By the same Rouch´ e argument applied to the circle { z : | z -z ∗ i | = ρ } , for all large t , f t has exactly k i zeros in { z : | z -z ∗ i | < ρ } and exactly k i zeros in D i . Therefore, all zeros of f t in D i lie in { z : | z -z ∗ i | < ρ } . Since ε is arbitrary, we obtain the convergence for t →∞ . Next, by Theorem 1, for each q ∈ [ k i ], there exist random eigenvalues z ( t ) n,q of H j ( t ) such that $$\lim _ { n \rightarrow \infty } q \in [ k ] ^ { 2 } - z _ { i } ^ { * } ( t ) = 0 .$$ Combining them via the triangle inequality, we obtain $$\lim _ { t \rightarrow \infty } n _ { d } \rightarrow \infty q ^ { [ k ] } \{ k \} = 0 ,$$ which is exactly the eigenvalue convergence statement in Theorem 2. Step III. Eigenspace convergence. Fix i ∈ { 1 , . . . , p } and consider a large enough time t . We use the same notation as in Step II. Let { ξ ( t ) n,ℓ } k i ℓ =1 denote orthonormal eigenvectors of H j ( t ) corresponding to eigenvalues inside γ i = ∂ D i . Denote the distinct zeros of f t ( z ) in D i as { z ( t,r ) i } q i ( t ) r =1 with multiplicities m i,r ( t ) satisfying ∑ r m i,r ( t ) = k i . For each r , let U ( t,i,r ) ⩽ , H ∈ R r × m i,r ( t ) be an orthonormal basis for the kernel $$\delta G ^ { i 2 } _ { t } - \frac { \vert z _ { t } \vert ^ { 2 } I _ { r } + E } { \left | V _ { * } \right | ^ { 2 } } .$$ By Theorem 1, for fixed t , $$\sum _ { ( t , i , r ) } \sum _ { ( t , i , r ) } | e ^ { i H S _ { H } ( t , i ) } | ^ { 2 } - \sum _ { ( t , i , r ) } T r = 1$$ as n, d →∞ , where and α j ′ t ( z ) := ∂ ∂z α j t ( z ). Since z ( t,r ) i → z ∗ i by Step II, dominated convergence yields $$S _ { H } ( t ; z ) = E \left[ \frac { \delta G ^ { 2 } _ { t } c ^ { i } _ { z } } { ( s + G _ { d } c _ { t } ) ^ { 2 } } V _ { < } ^ { T } \right]$$ We now take t →∞ . By Assumption 3, we have G j t → G j ∞ in L 2 (also in Wasserstein-2 sense), α j t → α j ∞ , and α j ′ t → α j ′ ∞ uniformly on K . $$\frac { S ^ { 1 } H ( t ; z _ { i } ) } { \left | + S ^ { 1 } H ( c ; z _ { i } ) \right | } \rightarrow S ^ { 1 } H ( t ; z _ { i } ) as t \rightarrow \infty .$$ ′ ∞ ′ H ∞ E δG ∞ α ( z ) ( δ + j j G ∞ j α j ∞ ( z )) ∗ Concatenating the bases, we form $$\begin{aligned} u ( t , i , q ( t ) ) & = [ u _ { \leq i } H ^ { t } ] \cdot \ldots \cdot u \\ & = R ^ { t } x _ { k _ { i } } , S _ { i } := \text{span} ( U ( t , i , H ^ { t } ) ) \end{aligned}$$ Since A t ( · ) → A ∞ ( · ) uniformly on K (Step I) and z ( t,r ) i → z ∗ i (Step II), we have A t ( z ( t,r ) i ) → A ∞ ( z ∗ i ) for each r . By the Davis-Kahan theorem (applied to the eigenspace of A ∞ ( z ∗ i ), which has dimension k i ), the corresponding k i -dimensional subspaces S i t converge as t →∞ with the limit $$E \left[ \frac { s G ^{2} _{\infty } } { s + G ^{2} _{\infty } d \alpha ( z _ { i } ) } \right] .$$ Let U i ∞ ∈ R r × k i be an orthonormal basis of S i ∞ . Since the trace formula U ↦→ Tr( U ( U T B i ∞ U ) -1 U T ) is continuous on the Stiefel manifold and basis-invariant, we obtain which matches the eigenvector correlation formula in Theorem 2. Step IV. Regime δ < δ ∗ j ( ∞ ) . By definition (48), f ∞ ( z ) = 0 has no negative solutions in the region z < c j ( ∞ ). Since f ∞ is continuous and nonvanishing on the compact set K , $$f : = inf \{ f _ { \alpha } ( z ) \} > 0 .$$ ] where S ( ; z ) := [ V ⩽ V T ⩽ ∗ . Define B i ∞ := I r + S ′ H ( ∞ ; z ∗ i ). By Lemma C.20, for all sufficiently large t , we have sup f ( z ) f ( z ) < m/ Therefore, for all z ∈ K z ∈K | t -∞ | 2. , $$( 3 0 0 ) - f _ { \infty } ( z ) \geq \beta - \beta / 2 = \beta / 2 > 0$$ which means f t ( z ) = 0 for all z ∈ K . By Theorem 1 (the case δ < δ ∗ j ( t )), for each large t , with probability 1 -o (1), the Hessian H j ( t ) has no eigenvalues in K arising from the hard block. The conclusion follows since K is arbitrary. ̸ The conclusion for eigenvectors also follows from the finite-time result. Taking t large enough and noting that the choice of K is arbitrary (as long as K ⊂ D ∞ ), we obtain $$\frac { p - \lim _ { n \rightarrow + \infty } | e ^ { i T } \sin \frac { 1 } { n } p || ^ { 2 } } { n d + \infty } = 0$$ for each fixed p ≥ 1 and all large enough t , where ξ p is the eigenvector associated with the p -th smallest eigenvalue of H j ( t ). This completes the proof of Theorem 2. Proof of Lemma C.20. Note that V ( t ) → V ( ∞ ) in L 2 by Assumption 3. Under Assumption 1, the map v ↦→ g ( v, h ( V ∗ , ε ); j ) is bounded and Lipschitz, therefore also G j t → G j ∞ in L 2 (in particular, G j t ⇒ G j ∞ ). This yields µ j ∞ ( t ) ⇒ µ j ∞ ( ∞ ) and c j ( t ) → c j ( ∞ ) (which is the statement of Lemma 2.5). Since K ⊂ D ∞ is at strictly positive distance from supp( µ j ∞ ( ∞ )), for all sufficiently large t we have K ⊂ D t := { z : Re( z ) < min(0 , c j ( t )) } . Therefore, $$\sup \{ \alpha ^ { i } ( z ) - \alpha ^ { j } _ { \infty } ( z ) \} \rightarrow 0 as t \rightarrow \infty$$ $$\frac { \delta G _ { i } ^ { 2 } } { \delta t } = - z I _ { r } + E$$ Define and similarly A ∞ ( z ) with t replaced by ∞ . By the determinant Lipschitz bound, $$\int _ { \alpha } ^ { \infty } f ( z ) - f _ { \infty } ( z ) | z - A | d p$$ From the definition of K , we have the uniform boundedness of G j t , G j ∞ , α j t ( z ), α j ∞ ( z ) on K for large enough t , and also the uniform boundedness of the operator norms ∥ A t ( z ) ∥ op and ∥ A ∞ ( z ) ∥ op over z ∈ K and large enough t . Therefore, it suffices to show sup z ∈K ∥ A t ( z ) -A ∞ ( z ) ∥ op → 0. Define ϕ ( g, α ) := δg 2 δ + gα , such that $$\begin{array}{ll} A _ { t } ( z ) - A _ { \infty } ( z ) = E [ \phi ( G _ { i } ^ { i } , a _ { i } ( z ) ) - \phi ( G _ { \infty } ^ { \infty } , a _ { \infty } ( z ) ) ] V _ { * } < V _ { 0 }$$ For fixed z , we have G j t → G j ∞ in L 2 and α j t ( z ) → α j ∞ ( z ). The function ϕ is Lipschitz continuous and absolutely bounded in K , and ∥ V ∗ ⩽ ∥ 2 is integrable. By dominated convergence, $$\begin{array}{ll} \lvert A _ { t } ( z ) - A _ { \infty } ( z ) \rvert_{op} & = 0 & \text{for each fixed } z \in K. \\ \end{array}$$ To obtain uniform convergence over z ∈ K , note that by (302), the family { α j t ( · ) } t ≥ 1 is equicontinuous on K for large enough t . Combined with the Lipschitz continuity of ϕ ( · ) in α uniformly over the possible values of g , the functions z ↦→ A t ( z ) are also equicontinuous for large enough t . Therefore, pointwise convergence (306) implies $$\sup \limits _ { z \in K } | A _ { t } ( z ) - A _ { \infty } ( z ) | \oplus 0 a s t$$ which yields the result by (304). ## D Proofs of supporting lemmas This appendix provides proofs of the supporting technical lemmas used in Appendix C. ## D.1 Proof of Lemma C.6 Proof. First, define the high-probability event $$\sum _ { i = 1 } ^ { n } ( R e ( z ) \leq \lambda _ { m i n } ( R G R ^ { T } ) - c )$$ where ε K > 0 is the gap defined in the beginning of Section C.2. By Lemma A.4, we have P ( E ) → 1 as n, d →∞ . On the event E , both resolvents are well-defined and satisfy the operator-norm bounds $$\begin{array}{ll} \vert ( R \overline { G } R ^ { T } - z I ) ^ { 1 } \vert _ { \oplus } & = 1 \\ \vert ( R \overline { G } R ^ { T } - ( z + i m ) I ) ^ { 1 } \vert _ { \oplus } & = - e x \end{array}$$ Second, define the bilinear forms $$S ( z ) := \sum _ { n = 0 } ^ { \infty } S ( z + i n ) := c$$ Using the fact that A -1 -B -1 = B -1 ( B -A ) A -1 with A = R ̂ GR T -z I and B = R ̂ GR T -( z + iη n ) I , we obtain $$\begin{array}{ll} \{ S ( z ) - S ( z + i n ) \} = n _ { m } \\ \leq n _ { e } \| a \| || b || R ^ { 2 } _ { lop } \{ ( } \begin{array}{l} a ^ { 2 } || b || R ^ { 2 } _ { lop } . \end{array} \right. \\ \leq \frac { n _ { e } } { 2 } \| a \| || b || R ^ { 2 } _ { lop } . \end{array} Since R \in R ^ { n \times d } has i.i.d. (bounded (309)) (311)$$ Since R ∈ R n × d ⋆ has i.i.d. N (0 , 1) entries (where we recall that d ⋆ = d -k -m -mt ), there exists a constant C such that, with probability 1 -o (1) [Ver18] $$\vert R \vert _ { \vert o p } \leq C ( \sqrt { n } + \sqrt { d } )$$ for some absolute constant C > 0. Conditioning on the events E and (312) (which together hold with probability 1 -o (1)), and substituting into (311), we obtain $$\left | S ( z ) - S ( z + i n ) \right |$$ for some universal constant C ′ > 0 independent of n, d . To conclude, we analyze the right-hand side of (313). By the assumption of Lemma C.6, we have ∥ a ∥ , ∥ b ∥ = O P (1 / √ n ), which implies ∥ a ∥∥ b ∥ = O P (1 /n ). In addition, we have ( √ n + √ d ) 2 = O ( n ). Combining these observations, the bound (313) becomes $$\frac { \vert S ( z ) - S ( z + i m ) \vert } { | e ^ { i k } | }$$ Since η n = 1 / log n → 0 as n →∞ , we conclude S ( z + iη n ) -S ( z ) p - - → 0. The convergence holds unconditionally since both events E and (312) occur with probability 1 -o (1). ## D.2 Proof of Lemma C.7 Proof. We establish the bound via decomposing the conditional second moment into diagonal and off-diagonal terms, then invoking Gaussian moment calculations and the leave-one-out technique. Throughout, we use the resolvent bound ∥ Q -i ∥ op ≤ (1 /η n ) (since Im( z + iη n ) ≥ η n ). ̸ First, recall from Section C.2 that $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { m } h _ { i j } r _ { i } q _ { j } , h _ { i j } := n$$ The second moment can be expanded as $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { i - 1 } [ T_{i , j } F_i ] + \sum _ { i = 1 } ^ { n } [ S_{cross} S_{cross} F_i ]$$ We bound the two terms respectively. For the off-diagonal terms, for i = j , we expand the product ̸ ̸ $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } E [ T_i \overline{T_j} | F_i ] = \frac { 1 } { n ^ { 4 } } \sum _ { k = 1 } ^ { n } E [ T_i \overline{T_j} | F_i ]$$ ̸ where all other cross-indexed terms have zero expectations. ̸ Conditioning on { r k } k = i (which determines both Q -i and r j since j = i ), and applying the leave-one-out formula r T i Qr j = r T i Q -i r j 1+ ̂ g i r T i Q -i r i , we obtain ̸ ̸ $$\begin{array}{ll} E \left[ h_{ij}^2 E \left[ r_i^T Q - r_j^T Q \right] \right]^2 & = E \\ & = E \\ \end{array}$$ ̸ ̸ where d ⋆ = d -k -m -mt . Summing over all i = j and defining M n := (1 /n 2 ) ∑ i,j h 2 ij = n 2 ∥ a ∥ 2 ∥ b ∥ 2 , we obtain the off-diagonal bound ̸ ̸ $$\sum _ { i = j } ^ { n } [ r _ { i } \overline { T _ { i } } q _ { r _ { i } } ] ^ { 2 } | F _ { i } | \leq d _ { i } M _ { n } C _ { i + j }$$ ̸ ̸ ̸ For the diagonal terms, we write T i using the leave-one-out formula as ̸ $$\frac { T _ { i } } { r _ { i } ^ { 1 } } = \frac { 1 } { n ^ { 2 } 1 + q _ { i r _ { i } } T _ { i } Q } .$$ ̸ Taking the expectation for | T i | 2 and conditioning on { r j } j = i , we obtain ̸ $$\sum _ { j = 1 } ^ { n } \sum _ { i = 1 } ^ { n } h_{ij} r_i$$ ̸ Summing over all i ∈ [ n ], we obtain $$\sum _ { i = 1 } ^ { n } \sum _ { j = 1 } ^ { n } E [ | T_i | ^ { 2 } | F_t | < \sum _ { i = 1 } ^ { n } d + C ^ { 2 } _ { n } n ^ { 4 } n A \sum _ { j = 1 } ^ { n } j ^ { 2 } i < \sum _ { i = 1 } ^ { n } d + M _ { n } C ^ { 2 } _ { n } n ^ { 2 } n A$$ ̸ Finally, combining the diagonal bound (322) and the off-diagonal bound (319) into the equation (316), we conclude $$\frac { \int _ { n } ^ { 2 } S_{n} d t } { \sum _ { n = 1 } ^ { 2 } F _ { i } } < - \frac { E [ | S_{cross } | ^ { 2 } ] } { M _ { n } C ^ { 2 } _ { n } } + \frac { d _ { 4 } M _ { n } C ^ { 2 } _ { n } } { n ^ { 2 } n ^ { 4 } } = \frac { M _ { n } C ^ { 2 } _ { n } } { n ^ { 2 } n ^ { 4 } } .$$ This completes the proof. ## D.3 Proof of Lemma C.9 Proof. Define the companion resolvent $$Q ( z ) = ( G R ^ { I } R - z I _ { n } ) ^ { - 1 } e ^ { i R ^ { n } }$$ By the Woodbury matrix identity (or direct verification), we have the following resolvent relation $$Q ( z ) = - 1 _ { n } + - G R _ { z }$$ This formula couples the two resolvents acting on spaces of different dimensions. Taking the ( i, i ) diagonal entry of both sides of (325), we obtain $$\frac { ( Q ( z ) ) _ { i } } { 1 - \frac { 1 } { z } + \frac { q _ { i } } { z ^ { r } } Q ( z ) _ { r } } ,$$ ̸ ̸ ̸ ̸ where r i ∈ R d ⋆ is the i -th row of R and ̂ g i is the i -th diagonal entry of ̂ G . We then invoke the leave-one-out formula, $$r _ { 2 } ^ { T Q ( z ) r _ { 1 } } = \frac { r _ { 1 } Q - ( z ) r _ { 1 } } { 1 + g _ { i } r _ { 1 } Q - ( z ) r _ { 1 } } ,$$ and substituting into (326), we obtain $$\frac { ( Q ( z ) ) _ { i } } { 1 - \frac { 1 } { z ^ { 2 } + g _ { r } ^ { 2 } r } Q _ { i } ( z ) _ { r } }$$ $$1 + q _ { i r } ^ { T } Q - i ( z ) r = - \frac { 1 } { z ( Q ( z ) ) i i }$$ For z + iη with η > 0, the matrix ̂ GR T R -( z + iη ) I n has eigenvalues bounded away from zero by η . Therefore, the following resolvent bound holds $$\begin{aligned} \lvert Q ( z + i n ) \rvert _ { \oplus } = \left\lVert \frac{1}{z + i n} \right\rVert = 1 . \end{aligned}$$ Rearranging yields From (329), we have for some constant C ′ > 0 depending only on the magnitude of z . Since the considered domain of z is bounded, C ′ is an absolute constant. ## D.4 Proof of Lemma C.11 Proof. Write ∆ := n/δ -d ⋆ . Since n/δ = δ n d/δ = d + o ( d ) and d ⋆ = d -C 0 , we have | ∆ | = o ( d ). In addition, we have $$\frac { 1 + d _ { g } i \alpha ( z ) } { 1 + d _ { g } i \alpha ( z ) - \frac { a ( z ) } { 1 + n g } } = \frac { a ( z ) } { 1 + d _ { g } i \alpha ( z ) }$$ By Lemma C.9, there exists c K > 0 such that | 1+ d ⋆ ̂ g i α ( z ) | ≥ c K and | δ + n ̂ g i α ( z ) | ≥ δc K uniformly for z ∈ K and i ∈ [ n ]. Since | ̂ g i | ≤ C/n , | α ( z ) | ≤ C K , and | ∆ | = o ( d ), we obtain from (332) $$\frac { \vert a ( z ) \vert } { 1 + d _ { 4 g i } ( a ( z ) ) } - \frac { \delta ( z ) } { c _ { k } } = \frac { a ( z ) } { 1 + d _ { 4 g i } ( a ( z ) ) } - \frac { \delta ( z ) } { c _ { k } }$$ Summing over i ∈ [ n ] and using ∥ a ∥∥ b ∥ = O P (1 /n ) we get $$\sum _ { i = 1 } ^ { n } \left( \frac { a ( z ) } { 1 + d g _ { i } a ( z ) } - \frac { b ( z ) } { 1 + d g _ { i } a ( z ) } \right)$$ which establishes (223). ## D.5 Proof of Lemma C.15 Proof. Recall the outlier matrix from Section C.1: $$M _ { 3 } ( z ; t ) = \left[ Z ^ { T } G Z I \right]$$ where Z = X Θ ∗ ; t ( Θ T ∗ ; t Θ ∗ ; t ) -1 , Y = P ⊥ F 0; t X new P ⊥ Θ ∗ ; t , G = G j ( t ), and Θ ∗ ; t , F 0; t are as defined in Section C.1. The resolvent R 0 ( z ) := ( Y T GY -z I ) -1 was introduced in Section C.1. By the resolvent identity, for z, z ′ ∈ K , $$R _ { 0 } ( z ) - R _ { 0 } ( z ^ { \prime } ) = ( z ^ { \prime } - z ) R _ { 0 } ( z ) R$$ From the Lemma A.4, for z ∈ K we have min i ∈ [ n ] | λ i ( Y T GY ) -z | ≥ ε K with probability 1 -o (1), where ε K > 0 depends on K but is uniform over n, d . Therefore, $$\sup \limits _ { x \in K } | R _ { 0 } ( z ) | \oplus \vert o p = - 1$$ Combining (336) and (337), we obtain $$\frac { \vert z - z ' \vert } { e ^ { 2 x } } with probability 1 - o ( 1 = \frac { | z - z ' | } { e ^ { 2 x } }$$ Define the block matrices $$a = \left[ \theta _ { t } , Y G Z \right] , C = \left[ z ^ { T } G Z , I \right]$$ From (335), the z -dependence enters only through R 0 ( z ): $$= A ( R _ { 0 } ( z ) - R _ { 0 } ( z ) ) B .$$ By submultiplicativity and (338), $$\vert M _ { j } ( z ; t ) - M _ { j } ( z ^ { \prime } ; t ) \vert o p$$ By Lemma 2.1 and the results established in Section C.1, $$\begin{array}{ll} || A || _ { \mathrm {op } } = O_P ( 1 ) , & || B || _ { \mathrm {op } } = O_P ( 1 ) , \\ z \in K. \end{array}$$ Therefore, from (341), $$\| M _ { j } ( z ; t ) - M _ { j } ( z ^ { \prime } ; t ) \|$$ Finally, invoke the standard determinant Lipschitz bound: for matrices X , Y ∈ R p × p , $$( | X || p ^ { - 1 } + || Y || p ^ { - 1 } ) | X - Y | _ { \sigma }$$ Here for M j ( z ; t ) the dimension p = O (1). Applying (344) with X = M j ( z ; t ) and Y = M j ( z ′ ; t ), and using (342) and (343), we obtain (251). ## D.6 Proof of Lemma C.19 We first establish that M j ( z ; t ) -1 converges to M ∞ j ( z ; t ) -1 uniformly over z ∈ γ . Lemma D.1. Choose the contour γ such that inf z ∈ γ | det( M ∞ j ( z ; t )) | > 0 . Then $$\sup _ { z \in T } M _ { j } ( z ; t ) - 1 - M ^ { p } _ { j } ( z ; t ) - 1 \uparrow$$ Proof. By Lemma C.1, we have sup z ∈ γ ∥ M j ( z ; t ) -M ∞ j ( z ; t ) ∥ op p - - → 0. Since inf z ∈ γ | det( M ∞ j ( z ; t )) | > 0 and det( M ∞ j ( z ; t )) is continuous on the compact set γ , there exists c γ > 0 such that sup z ∈ γ ∥ M ∞ j ( z ; t ) -1 ∥ op ≤ c -1 γ . By uniform convergence, with probability 1 -o (1), we have inf z ∈ γ | det( M j ( z ; t )) | ≥ c γ / 2 > 0, hence sup z ∈ γ ∥ M j ( z ; t ) -1 ∥ op ≤ C γ for some C γ < ∞ , and large enough n . Using the identity M j ( z ; t ) -1 -M ∞ j ( z ; t ) -1 = M j ( z ; t ) -1 ( M ∞ j ( z ; t ) -M j ( z ; t )) M ∞ j ( z ; t ) -1 and submultiplicativity, $$\| M _ { j } ( z ; t ) - 1 \|_{op} \| M ^ { p } _ { j } ( z ; t ) - M _ { j } ( z ; t ) \|_{op} \leq C _ { r } \cdot \| M _ { j } ( z ; t ) - M ^ { p } _ { j } ( z ; t ) \|_{op} \cdot c _ { r } ^ { 1 } ,$$ which converges to zero uniformly in probability over z γ by Lemma C.1. ∈ Proof of Lemma C.19. We decompose F n,d ( z ) as a product of three parts: $$F _ { n , d } ( z ) = A ( z ) ^ { T } \cdot M _ { j } ( z ; t ) ^ { - 1 } \cdot A$$ $$\begin{aligned} A ( z ) := \left[ O _ { t } ^ { t } R _ { 0 } ( z ) \Theta _ { t } H \right] = T ^ { - 1 } R _ { 0 } ( z ) G Y R _ { 0 } ( z ) \Theta _ { t } H \in I ^ { 2 } ( k + m + mt ) \times r . (348)$$ where We establish uniform convergence with respect to z ∈ γ for A ( z ). First, by Section C.1, specifically (174), we have $$( Y ^ { T } G Y - z I ) ^ { 1 } \theta _ { * , t } = - \theta ^ { T } _ { * , t } ( \theta _ { * } )$$ By Lemma 2.2, Θ T ∗ H Θ ∗ ; t p - - → [ I r 0]. Therefore, $$\rho _ { z } ^ { 1 } \Omega _ { * , t } ( z ) \Theta _ { * , H } = - \rho _ { z } ^ { 1 } [ r _ { 0 } ] ^ { T }$$ uniformly over z ∈ γ (since γ is compact and bounded away from z = 0). Second, recall that Z T GYR 0 ( z ) Θ ∗ ; t = 0. Since Θ ∗ H consists of the first r columns of Θ ∗ ; t , we have immediately $$Z ^ { \prime } G Y R _ { 0 } ( z ) \theta _ { 0 } \mu = 0$$ Combining (350) and (351), we obtain $$\sup _ { z \in Y } | A ( z ) - A _ { \infty } ( z ) | \to 0 .$$ $$A _ { \infty } ( z ) := - 1 [ r ^ { | z | } ( e ^ { i k + m + n t } ) ,$$ for all z ∈ γ . where By Lemma D.1, we have $$\sup _ { z \in Y } | M _ { j } ( z ; t ) - 1 - M ^ { p } _ { j } ( z ; t ) - 1 ||$$ Moreover, A ( z ), M j ( z ; t ) -1 , and their limits are uniformly O P (1) on γ . From the decomposition $$F _ { n , d } ( z ) - F _ { \infty } ( z ) = A ( z ) = ( A _ { 1 } + A _ { \infty } ) + A _ { \infty }$$ we obtain by submultiplicativity and (352), (354) that $$\sup _ { z \in Y } \left | F _ { n , d } ( z ) - F _ { \infty } ( z ) \right | o p \rightarrow 0 .$$ which completes the proof of (273). ## E Additional experiments ## E.1 Results for networks with moderately large width Numerically verifying our results and thresholds for large width networks is computationally prohibitive. Nevertheless, we provide numerical simulations for moderately large width networks to give insights into the phase transition phenomena. We extend our experiments to multi-neuron networks. For these experiments, we still use GeLU activation, Huber loss with parameter M = 1, and fixed learning rate η = 1 . 5. We consider widths m ∈ { 1 , 2 , 3 , 5 , 10 } with a j = 1 and b j = 0 for all neurons j ∈ [ m ]. The correlation is defined as the maximum over all m neurons: $$\rho ( t ) = \max _ { j \in [ m ] } | \theta _ { j } ( t ) | ^ { 2 }$$ Figure 8 presents the success rate and the final correlation for m = 5 across varying dimensions, exhibiting sharp phase transitions around δ = 3 . 5. We further investigate the proportional scaling regime where the network width grows proportionally with dimension, maintaining a fixed ratio m/d = 0 . 01. In these experiments, we consider pairs ( m,d ) ∈ { (4 , 400) , (6 , 600) , (8 , 800) , (10 , 1000) , (12 , 1200) } , with the same experimental setup as before: GeLU activation, Huber loss with parameter M = 1, learning rate η = 1 . 5, and all biases set to zero. Figure 9 shows the phase transition behavior under proportional width scaling. Interestingly, we observe that as both m and d increase proportionally, the phase transition threshold shifts a bit to the right. Figure 10 (left panel) shows success rate versus δ for different widths at fixed dimension d = 1000. As width increases, the phase transition threshold changes. Whether the threshold continues to decrease as width increases remains an open question due to computational constraints that limit us to moderately large widths and dimensions. Figure 10 (right panel) shows the final correlation for the same widths at fixed dimension d = 1000. <details> <summary>Image 8 Details</summary>  ### Visual Description ## Line Graphs: Success Rate vs. Final Correlation by Dataset Size (d) ### Overview The image contains two side-by-side line graphs comparing performance metrics (success rate and final correlation) across different dataset sizes (d = 1000, 2000, 3000, 4000). Both graphs plot metrics against δ = n/d, where n/d represents a normalized parameter. The left graph focuses on success rate (threshold = 0.5), while the right graph measures final correlation. Data points are plotted with markers and connected by lines, with error bars indicating variability. --- ### Components/Axes - **Left Graph (Success Rate)**: - **X-axis**: δ = n/d (ranges from 2 to 6) - **Y-axis**: Success Rate (Threshold = 0.5) (ranges from 0.0 to 1.0) - **Legend**: Located in the top-right corner, with colors corresponding to d values: - Light blue: d = 1000 - Medium blue: d = 2000 - Dark blue: d = 3000 - Very dark blue: d = 4000 - **Markers**: Circular data points with error bars. - **Right Graph (Final Correlation)**: - **X-axis**: δ = n/d (ranges from 2 to 6) - **Y-axis**: Final Correlation (ranges from 0.0 to 1.0) - **Legend**: Located in the top-right corner, with colors corresponding to d values: - Light green: d = 1000 - Medium green: d = 2000 - Dark green: d = 3000 - Very dark green: d = 4000 - **Markers**: Square data points with error bars. --- ### Detailed Analysis #### Left Graph (Success Rate) - **Trends**: - All d values show a sigmoidal increase in success rate as δ increases. - Success rate plateaus near 1.0 for δ ≥ 4.5 across all d values. - Lower d values (e.g., d = 1000) achieve plateau later than higher d values (e.g., d = 4000). - Error bars are smallest for d = 4000, indicating higher precision in larger datasets. - **Key Data Points**: - At δ = 3.0: - d = 1000: ~0.65 - d = 2000: ~0.70 - d = 3000: ~0.75 - d = 4000: ~0.80 - At δ = 4.0: - All d values exceed 0.90, with d = 4000 reaching ~0.95. #### Right Graph (Final Correlation) - **Trends**: - All d values show a gradual increase in final correlation as δ increases. - Correlation plateaus near 0.95 for δ ≥ 4.5 across all d values. - Lower d values (e.g., d = 1000) achieve plateau later than higher d values (e.g., d = 4000). - Error bars are smallest for d = 4000, indicating higher precision in larger datasets. - **Key Data Points**: - At δ = 3.0: - d = 1000: ~0.30 - d = 2000: ~0.35 - d = 3000: ~0.40 - d = 4000: ~0.45 - At δ = 4.0: - All d values exceed 0.70, with d = 4000 reaching ~0.85. --- ### Key Observations 1. **Dataset Size Impact**: Larger datasets (higher d) consistently outperform smaller datasets in both success rate and final correlation across all δ values. 2. **Threshold Behavior**: The success rate threshold of 0.5 is surpassed rapidly for δ ≥ 3.0, with larger datasets achieving it earlier. 3. **Diminishing Returns**: Both metrics plateau at high δ values, suggesting limited gains beyond δ = 4.5. 4. **Error Bars**: Smaller error bars for larger d values indicate more reliable measurements in bigger datasets. --- ### Interpretation The data demonstrates that increasing dataset size (d) improves both success rate and final correlation, particularly at lower δ values. This suggests that larger datasets enhance model robustness or generalization. The threshold of 0.5 in the success rate graph may represent a critical operational benchmark, where larger datasets achieve this threshold more efficiently. The plateauing trends imply that beyond a certain δ (n/d), additional data provides diminishing returns. The consistent performance across d values highlights the importance of dataset size in optimizing these metrics, with practical implications for resource allocation in data-driven systems. </details> Figure 8: Phase transitions for moderately large width learners with GeLU activation. Left: Success rate for m = 5, showing sharp transitions around δ ≈ 3 . 5. Right: Final correlation. Figure 9: Phase transitions for proportional width scaling with m/d = 0 . 01. Left: Success rate (threshold 0 . 5). Right: Final correlation (median with 30th/70th percentiles). <details> <summary>Image 9 Details</summary>  ### Visual Description ## Line Graphs: Success Rate vs. Final Correlation by δ = n/d ### Overview The image contains two side-by-side line graphs comparing performance metrics (Success Rate and Final Correlation) across varying values of δ = n/d (ranging from 2 to 6). Both graphs show multiple data series differentiated by system parameters (m and d), with error bars indicating measurement uncertainty. ### Components/Axes - **Left Graph (Success Rate)**: - **X-axis**: δ = n/d (2 to 6, linear scale) - **Y-axis**: Success Rate (Threshold = 0.5) (0.0 to 1.0, linear scale) - **Legend**: - m=4, d=400 (light blue circles) - m=6, d=600 (medium blue squares) - m=8, d=800 (dark blue diamonds) - m=10, d=1000 (cyan triangles) - m=12, d=1200 (navy pentagons) - **Legend Position**: Bottom-right corner - **Right Graph (Final Correlation)**: - **X-axis**: δ = n/d (2 to 6, linear scale) - **Y-axis**: Final Correlation (0.0 to 1.0, linear scale) - **Legend**: - m=4, d=400 (light green circles) - m=6, d=600 (medium green squares) - m=8, d=800 (dark green diamonds) - m=10, d=1000 (cyan triangles) - m=12, d=1200 (navy pentagons) - **Legend Position**: Bottom-right corner ### Detailed Analysis #### Left Graph (Success Rate) - **Trend**: All lines exhibit a sigmoidal curve, starting near 0.0 at δ=2 and plateauing near 1.0 by δ=5–6. - **Key Data Points**: - m=4, d=400: Reaches ~0.95 at δ=4.5, with error bars ±0.03. - m=12, d=1200: Plateaus at 1.0 by δ=4, error bars ±0.01. - **Notable**: Higher m/d values achieve higher success rates earlier (e.g., m=12,d=1200 plateaus at δ=4 vs. m=4,d=400 at δ=4.5). #### Right Graph (Final Correlation) - **Trend**: Similar sigmoidal shape but with shallower slopes. All lines plateau near 0.9–1.0. - **Key Data Points**: - m=4, d=400: Reaches ~0.85 at δ=5, error bars ±0.04. - m=12, d=1200: Plateaus at ~0.98 by δ=5.5, error bars ±0.02. - **Notable**: Correlation improves with higher m/d but plateaus below 1.0, suggesting diminishing returns. ### Key Observations 1. **Performance Scaling**: Both metrics improve with increasing m and d, but higher m/d values achieve plateaus faster. 2. **Threshold Sensitivity**: The left graph’s threshold (0.5) is exceeded earlier for larger systems (e.g., m=12,d=1200 reaches 1.0 by δ=4). 3. **Error Margins**: Larger systems (m=12,d=1200) exhibit tighter error bars, indicating more consistent measurements. ### Interpretation The data suggests that system performance (success rate and correlation) improves with increased m and d parameters, likely due to enhanced computational capacity or data granularity. The left graph’s threshold (0.5) is surpassed more efficiently by larger systems, while the right graph’s final correlation plateaus below 1.0, implying inherent limitations in the measured metric. The error bars highlight that larger systems yield more reliable results, reducing variability in outcomes. This could inform resource allocation decisions, where investing in higher m/d systems may optimize performance predictability. </details> Figure 10: Phase transitions for GeLU activation with varying width at dimension d = 1000. Left: Success rate with sigmoid fits showing how phase transitions change with increasing width. Right: Final correlation (median with 30th/70th percentiles) for varying widths. <details> <summary>Image 10 Details</summary>  ### Visual Description ## Line Graphs: Success Rate vs. Final Correlation by δ = n/d and m Values ### Overview The image contains two side-by-side line graphs comparing performance metrics (Success Rate and Final Correlation) across varying values of δ = n/d (x-axis) and m (legend parameter). Both graphs show sigmoidal trends with error bars, indicating statistical variability. The left graph focuses on success rates at a fixed threshold (0.5), while the right graph measures final correlation values. --- ### Components/Axes **Left Graph (Success Rate):** - **X-axis**: δ = n/d (ranges from 2 to 6, labeled in increments of 1) - **Y-axis**: Success Rate (Threshold = 0.5) (ranges from 0.0 to 1.0, labeled in increments of 0.2) - **Legend**: Positioned at bottom-right, with m values (1, 2, 3, 5, 10) mapped to light-to-dark blue gradients. - **Lines**: 5 curves (one per m value), all starting near 0 and plateauing near 1.0. **Right Graph (Final Correlation):** - **X-axis**: δ = n/d (same scale as left graph) - **Y-axis**: Final Correlation (ranges from 0.0 to 1.0, labeled in increments of 0.2) - **Legend**: Positioned at bottom-right, with m values (1, 2, 3, 5, 10) mapped to light-to-dark green gradients. - **Lines**: 5 curves (one per m value), all starting near 0 and plateauing near 1.0. --- ### Detailed Analysis **Left Graph (Success Rate):** - **m = 1**: Light blue line. Success rate rises sharply after δ ≈ 3.5, reaching ~0.95 by δ = 5.5. Error bars (±0.05–0.1) are largest at δ = 3–4. - **m = 2**: Medium blue line. Reaches ~0.9 by δ = 4, plateauing at ~0.95 by δ = 5. Error bars smaller than m = 1. - **m = 3**: Dark blue line. Achieves ~0.85 by δ = 3.5, plateauing at ~0.95 by δ = 4.5. Error bars minimal. - **m = 5**: Very dark blue line. Reaches ~0.9 by δ = 3, plateauing at ~0.98 by δ = 4. Error bars negligible. - **m = 10**: Darkest blue line. Peaks at ~0.95 by δ = 3, plateauing at ~0.99 by δ = 4. Error bars smallest. **Right Graph (Final Correlation):** - **m = 1**: Light green line. Correlation rises slowly, reaching ~0.8 by δ = 5.5. Error bars largest at δ = 3–4. - **m = 2**: Medium green line. Achieves ~0.75 by δ = 4.5, plateauing at ~0.9 by δ = 5. Error bars moderate. - **m = 3**: Dark green line. Reaches ~0.7 by δ = 4, plateauing at ~0.9 by δ = 5. Error bars smaller than m = 2. - **m = 5**: Very dark green line. Peaks at ~0.8 by δ = 4, plateauing at ~0.95 by δ = 5. Error bars minimal. - **m = 10**: Darkest green line. Achieves ~0.85 by δ = 4, plateauing at ~0.98 by δ = 5. Error bars smallest. --- ### Key Observations 1. **Inverse Relationship**: Higher m values achieve higher performance metrics (success rate/correlation) at lower δ values. For example: - m = 10 reaches ~0.95 success rate at δ = 3.5 (left graph), while m = 1 requires δ = 5.5. - m = 10 achieves ~0.98 final correlation at δ = 5 (right graph), while m = 1 only reaches ~0.8. 2. **Plateau Behavior**: All lines asymptotically approach 1.0, but higher m values plateau earlier. 3. **Error Trends**: Larger error bars for lower m values (e.g., m = 1–2) at δ = 3–4, suggesting greater variability in smaller systems. --- ### Interpretation The data demonstrates that increasing the parameter **m** improves system performance (success rate and final correlation) at lower δ values. This suggests: - **Efficiency Trade-off**: Higher m reduces the required δ to achieve high performance, but may involve computational or resource costs. - **Threshold Sensitivity**: The fixed threshold of 0.5 in the left graph highlights m's role in surpassing critical performance benchmarks. - **Scalability**: Systems with higher m values are more robust to variations in δ, as evidenced by smaller error bars and earlier plateaus. These trends imply that optimizing **m** is critical for balancing performance and resource constraints in systems governed by δ = n/d. </details> Figure 11: Grokking dynamics for different activations. Left: Quad with d = 5000, δ = 10, η = 0 . 15. Right: ReLU with d = 5000, δ = 17 . 5, η = 0 . 5. <details> <summary>Image 11 Details</summary>  ### Visual Description ## Line Graphs: Risk Reduction Over Training Time ### Overview Two side-by-side line graphs compare the evolution of risk metrics during training under two activation functions: "Quad" (left) and "ReLU" (right). Each graph tracks three metrics: 1. **1 - ρ(t)** (blue line) 2. **Risk(θ(t))/Risk(θ(0))** (orange line) 3. **Risk_test(θ(t))/Risk_test(θ(0))** (red line) ### Components/Axes - **X-axis (Left Graph)**: "Training Time t (Quad)" with values from 0 to 100. - **X-axis (Right Graph)**: "Training Time t (ReLU)" with values from 0 to 100. - **Y-axis (Both Graphs)**: "1 - ρ(t)" (blue line) and ratios of risk at time *t* to initial risk (orange/red lines). - **Legends**: Positioned on the left of each graph, with color-coded labels matching the lines. ### Detailed Analysis #### Left Graph (Quad Activation) - **Blue Line (1 - ρ(t))**: Starts at 1.0, decreases sharply to ~0.2 by t=60, then plateaus. - **Orange Line (Risk(θ(t))/Risk(θ(0)))**: Begins at ~0.8, drops to ~0.4 by t=40, then stabilizes. - **Red Line (Risk_test(θ(t))/Risk_test(θ(0)))**: Starts at ~0.6, declines to ~0.3 by t=60, then plateaus. #### Right Graph (ReLU Activation) - **Blue Line (1 - ρ(t))**: Similar to the left graph but with a slower initial decline, reaching ~0.2 by t=80. - **Orange Line (Risk(θ(t))/Risk(θ(0)))**: Starts at ~0.9, decreases to ~0.5 by t=60, then plateaus. - **Red Line (Risk_test(θ(t))/Risk_test(θ(0)))**: Begins at ~0.7, drops to ~0.4 by t=60, then stabilizes. ### Key Observations 1. **Blue Line (1 - ρ(t))**: Represents the inverse of the risk function. Both graphs show a monotonic decline, but the ReLU graph has a slower initial drop. 2. **Orange Line (Training Risk Ratio)**: Both graphs show a sharp decline in training risk, but the ReLU graph starts higher and declines more gradually. 3. **Red Line (Test Risk Ratio)**: Test risk decreases similarly in both graphs, but the ReLU graph exhibits a slightly slower decline. 4. **Divergence at Early Training**: The ReLU graph shows higher initial values for all metrics compared to the Quad graph, suggesting slower early progress. ### Interpretation - **Risk Reduction Dynamics**: The blue line (1 - ρ(t)) indicates that the model’s risk reduction accelerates early in training (Quad) but slows under ReLU. This suggests ReLU may require longer training to achieve similar risk reduction. - **Training vs. Test Risk**: The orange and red lines (training and test risk ratios) both decline, but the test risk (red line) lags slightly behind training risk (orange line) in both cases. This implies potential overfitting, though the gap narrows over time. - **Activation Function Impact**: The Quad activation achieves faster risk reduction initially, while ReLU shows more gradual progress. This could reflect differences in gradient dynamics or optimization landscapes between the two activation functions. ### Notable Trends - **Plateauing Behavior**: All metrics eventually plateau, indicating diminishing returns in risk reduction after a certain training duration. - **Consistency Across Metrics**: The blue line (1 - ρ(t)) consistently aligns with the orange and red lines, suggesting a direct relationship between the inverse risk function and the risk ratios. ### Conclusion The graphs demonstrate that training time significantly impacts risk reduction, with activation function choice influencing the rate of convergence. Quad activation achieves faster initial progress, while ReLU may require extended training for comparable results. The test risk (red line) lags behind training risk, highlighting the need for careful regularization to avoid overfitting. </details> ## E.2 Grokking across activation functions Figure 11 demonstrates that the same grokking phenomenon discussed in the main text extends to other activation functions. We show learning dynamics for Quad (left, δ = 10, η = 0 . 15) and ReLU (right, δ = 17 . 5, η = 0 . 5), both with d = 5000. Both experiments exhibit the two-stage pattern: an initial phase where the misalignment 1 -ρ ( t ) remains close to 1 with a constant gap between training and test losses, followed by a rapid transition where the network learns the signal direction and achieves near-perfect generalization.