S HARP DESCRIPTION OF LOCAL MINIMA IN THE LOSS LANDSCAPE OF HIGH - DIMENSIONAL TWO - LAYER R E LU NEURAL NETWORKS
arXiv:2604.09412v1 [stat.ML] 10 Apr 2026
A P REPRINT Jie Huang1 , Bruno Loureiro2 , and Stefano Sarao Mannelli* 3,4 1
Physics, Chalmers University of Technology and University of Gothenburg 2 Departement d’Informatique, École Normale Supérieure, PSL & CNRS 3 Data Science and AI, Computer Science and Engineering, Chalmers University of Technology and University of Gothenburg 4 School of Computer Science and Applied Mathematics, University of the Witwatersrand
A BSTRACT PK We study the population loss landscape of two-layer ReLU networks of the form k=1 ReLU(wk⊤ x) in a realisable teacher–student setting with Gaussian covariates. We show that local minima admit an exact low-dimensional representation in terms of summary statistics, yielding a sharp and interpretable characterisation of the landscape. We further establish a direct link with one-pass SGD: local minima correspond to attractive fixed points of the dynamics in summary statistics space. This perspective reveals a hierarchical structure of minima: they are typically isolated in the well-specified regime, but become connected by flat directions as network width increases. In this overparameterised regime, global minima become increasingly accessible, attracting the dynamics and reducing convergence to spurious solutions. Overall, our results reveal intrinsic limitations of common simplifying assumptions, which may miss essential features of the loss landscape even in minimal neural network models.
1
Introduction
Modern neural networks are trained by first-order algorithms on loss functions that are highly non-convex, structured, and shaped by both the architecture and the data distribution. Nevertheless, despite the inherent difficulty posed by non-convexity, day-a-day practice suggests that this is not a hindrance in the successful training of neural networks. Explaining this discrepancy requires a principled understanding of the geometry of these landscapes: the structure of their critical points, the nature of their basins of attraction, and the global organisation of minima and saddles. A principled understanding of these landscapes and how descent-based algorithms navigate them is therefore a central challenge for developing a mathematical theory of deep learning. Recent years have seen significant progress in this direction, most notably through results showing that the optimisation landscape of ReLU-based infinitely wide neural networks becomes benign in the mean-field limit, where global convergence guarantees can be established [Chizat and Bach, 2018, Mei et al., 2018, Rotskoff and Vanden-Eijnden, 2022, Sirignano and Spiliopoulos, 2020]. However, such asymptotic regimes remain largely non-quantitative and do not capture the finite-width mechanisms by which landscape trivialisation emerges as overparametrisation is increased [Bach and Chizat, 2021]. In particular, they offer limited guidance on how wide a network must be for these benign properties to hold, how landscape geometry changes as a function of model size, and which structural features of the data and target model control the onset of a benign optimisation regime. In this work, we address these questions by studying the population risk landscape of teacher–student two-layer neural networks with ReLU activation under Gaussian inputs. In this setting, Safran and Shamir [2018] established the existence of spurious local minima and provided compelling evidence that increasing the network width leads to a simplification of the landscape ruling out the hypothesis of a uniformly benign finite-width landscape with a computer-assisted construction. In a follow-up work, Safran et al. [2021] formally proved that mild overparameterisation fundamentally alters the local Hessian, converting non-global minima into saddle points. Yet, empirical evidence clearly ∗
Corresponding author: [email protected]
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
shows that the overparameterised landscape is not entirely benign, as new, higher-order local minima can form, leaving the geometry of these remaining traps largely unexplained by local analyses. By contrast, in this work we develop a refined description of the population landscape by leveraging a summarystatistics representation of the population loss landscape. This approach allows us to explicitly characterise the geometry of minima and saddle points, to analyse their stability, and to track how the landscape deforms as a function of overparameterisation. In doing so, we provide a detailed account of how landscape trivialisation emerges, complementing existing existence results with a geometric and quantitative understanding of the underlying optimisation problem. More precisely, our main contributions are: • We give an exact characterisation of the population loss landscape of two-layer ReLU networks via a reduced set of self-consistent equations for low-dimensional summary statistics, see Result 1. This reduced formulation enables an efficient analytical description of the minima and clarifies their relationships as fixed points of gradient-flow on the population loss. • Leveraging this description, we show that, in the well-specified setting where the teacher and student have equal width, the population loss of two-layer ReLU networks exhibits a hierarchy of isolated local minima separated by high-loss barriers. These minima are stable attractors for gradient-based algorithms and therefore constitute a significant obstacle to optimisation; see Fig.1A. (left) for an illustration. • We demonstrate that overparameterisation fundamentally alters this landscape topology, as illustrated in Fig.1A. (right). Using the string method [Weinan et al., 2002, Ren et al., 2007], we show that the previously isolated spurious minima expand into continuous, connected manifolds linked by flat paths, greatly facilitating convergence to global minima. • As shown in Fig.1B., these geometric properties directly govern the dynamics, which converge to quantised loss values characterised by the discrete hierarchy. Furthermore, the transition from the well-specified to the over-parameterised regime is reflected in the training statistics, where increasing the number of hidden units leads to a significantly larger fraction of trajectories reaching the global minimum. Taken together, our results deliver a geometric characterisation of the population loss landscape of two-layer ReLU networks and its consequences for optimisation. We show how overparametrisation controls the emergence, stability, and eventual disappearance of spurious minima, thereby providing a principled explanation for the optimisation benefits of overparameterisation.
Loss
Loss
A. Representation of the minima in the well-specified (left) and overparamterised (right) cases.
B. Minima density and associated value of the order parameters.
Figure 1: Geometry and statistics of local minima. (a) Schematic comparison of the loss landscape. In the wellspecified regime (left), minima are isolated points (marked ‘x’), whereas in the over-parameterised regime (right), they form continuous connected manifolds (red segments). (b) Validation of the theoretical predictions. The histogram shows the distribution of population risk reached by gradient flow (104 runs) in the high-dimensional standard normal teacher initialisations. Vertical dashed lines correspond to the analytical loss levels of the families derived in Result 2. Insets show heatmaps of the student-student (Q) and student-teacher (R) overlaps for a representative solution within the dominant family.
2
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Further related works High-dimensional Dynamics and Mean-Field Theory. The study of learning dynamics in two-layer neural networks using statistical physics techniques was pioneered by Biehl and Schwarze [1995], Saad and Solla [1995], who introduced the order-parameter description for soft committee machines. Recent works have rigorously established the validity of these mean-field descriptions in the high-dimensional limit [Goldt et al., 2019, 2020]. This framework has been extended to classify dynamical regimes based on the scaling of the learning rate and batch size [Ben Arous et al., 2022, Veiga et al., 2022, Arnaboldi et al., 2024], and to unify the description of single-pass, multi-pass SGD, and large-width [Arnaboldi et al., 2023a]. While these works focus on the evolution of the risk, our work leverages this formalism to characterise the fixed points of these dynamics, specifically analysing the hierarchical structure of minima that emerges from the non-convex landscape. Optimisation Landscape of Two-Layer Networks. A central question in deep learning theory is identifying conditions under which the loss landscape is benign. In the infinite-width limit, mean-field analyses [Mei et al., 2018, Chizat and Bach, 2018, Sirignano and Spiliopoulos, 2020, Rotskoff and Vanden-Eijnden, 2022] have shown that the loss landscape becomes convex-like, guaranteeing global convergence. However, this asymptotic result does not preclude the existence of spurious local minima in finite-width networks, as demonstrated empirically and theoretically by Safran and Shamir [2018], Safran et al. [2021] and Venturi et al. [2018]. Our results bridge this gap by providing a quantitative description of these finite-width spurious minima and explicating the geometric mechanism by which overparameterisation simplifies the landscape, a phenomenon qualitatively discussed in Draxler et al. [2018], Simsek et al. [2021] and rigorously analysed for phase retrieval in Sarao Mannelli et al. [2020a], Davis et al. [2020]. While Safran et al. [2021] relied on local Hessian evaluations to show that overparameterisation introduces negative eigenvalues at spurious minima, our summary-statistics approach provides a complementary, global view. We explicitly characterise the topology of these saddle-point transitions and track the full continuous manifolds connecting solutions. Symmetry, Connectivity, and Saddle Points Dynamics. The permutation symmetry of hidden units plays a crucial role in shaping the loss landscape. Fukumizu and Amari [2000] originally identified that hierarchical plateaus in learning curves arise from symmetry-breaking transitions. This hierarchical learning manifests as a “staircase” profile, where the dynamics are governed by a separation of timescales between the learning of distinct features [Abbe et al., 2022, Jain et al., 2024, Berthier et al., 2025, Montanari and Urbani, 2025]. In the specific context of ReLU networks, Arjevani and Field [2019] used group theory to show that spurious minima exhibit a “principle of least symmetry breaking,” resulting in weight matrices with highly symmetric, block-like transitivity partitions. However, their discrete algebraic approach does not easily translate to generalisation dynamics. Our framework naturally recovers these symmetric blocks via the macroscopic order parameters (Eq. 8), linking them directly to the test error and learning trajectories. More recent geometric analyses [Brea et al., 2019, Simsek et al., 2021] have argued that overparameterisation connects isolated minima through permutation saddles, effectively merging distinct basins of attraction. Similar connectivity properties have been explored in the context of the binary perceptron Baldassi et al. [2019, 2021], Annesi et al. [2023], Barbier et al. [2024], Barbier [2025]. Our use of the string method [Weinan et al., 2002, Ren et al., 2007, Samanta and Weinan, 2013] confirms this geometric picture in the context of ReLU networks, showing that the “flat directions” predicted by Baldassi et al. [2020] are realised as equipotential manifolds connecting symmetry-equivalent solutions. While these saddle points do not affect the long-time asymptotic performance, they play a pivotal role in the transient, often slowing down optimisation. This dynamical role of saddles has been investigated in polynomial networks [Zhang et al., 2025] and classical statistical-physics models [Cugliandolo and Kurchan, 1993, 1994]. Dynamics of Specialised and Simplified Architectures. Beyond the standard teacher-student setup, similar dynamical barriers and landscape complexities have been analysed in generalised linear models and multi-index models [Arous et al., 2021, Bietti et al., 2025, Şimşek et al., 2024], where the empirical risk can exhibit intricate spurious valleys [Maillard et al., 2020, Asgari et al., 2025]. Furthermore, the mechanism of saddle-to-saddle dynamics has been largely explored in deep linear networks [Saxe et al., 2014, 2019]. Recent results suggest that ReLU networks can exhibit similar behaviour: after an initial alignment phase [Boursier et al., 2022], the dynamics may follow a linear-network-like evolution, connecting non-convex training with the deep linear network formalism [Saxe et al., 2022, Jarvis et al., 2025]. Finally, in models like phase retrieval [Chen et al., 2019, Mannelli et al., 2019, Arnaboldi et al., 2023b] and tensor PCA [Sarao Mannelli et al., 2019, 2020b], the structure of the Hessian [Bonnaire et al., 2024, 2025], the initalisatisation scheme [Jarvis et al., 2025, Annesi et al., 2025], and the temperature of the dynamics [Baity-Jesi et al., 2018] are known to dictate the escape from spurious traps. These traps are often organised in complex topologies, such as the “canyon” structures observed in high-dimensional chaos [Fournier et al., 2025] or the triplets of minima found in spiked tensor models [Ros et al., 2019, Pacco et al., 2025]. Our work complements these findings by focusing on how these mechanisms manifest specifically in the ReLU nonlinearity via the compensation mechanism. 3
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
2
Problem Formulation and Summary Statistics Description
We consider the teacher-student framework, a convenient way of generating structured but tractable supervised learning tasks. In this framework, the target function y = f⋆ (x) — refereed as the teacher — is chosen from the same hypothesis class as the statistical model of interest — refereed as the student. In this work, we focus on the class two-layer neural networks with uniform readout layer and unbounded width. More precisely, each training example (xµ , y µ ) consists of an input vector xµ ∈ Rd and an output scalar y µ , where the components of xµ are drawn i.i.d. from the standard normal distribution N (0, 1). The target output is generated by a fixed teacher network with M hidden units: y µ = ϕ(xµ , W ∗ ), (1) ∗ where W ∗ ∈ RM ×d comprises the teacher weights wm (for m = 1, . . . , M ) and the network function is defined as ! M ∗ ⊤ X x w m √ ϕ(x, W ∗ ) = ReLU . (2) d m=1
Similarly, the student network approximates the target with K hidden units: ⊤ K X wk x √ ϕ(x, W ) = ReLU , d k=1
(3)
where the parameter matrix is W ∈ RK×d with rows corresponding to the weight vectors wk . We formally define the parameterisation regimes based on the relative width of the student network compared to the teacher. Definition (Parameterisation Regimes). Let M and K denote the number of hidden units in the teacher and student networks, respectively. 1. Well-specified: The regime where K = M . In this case, the function space representable by the student coincides exactly with that of the teacher. 2. Overparameterised: The regime where K > M . In this case, the function space representable by the student includes and exceeds that of the teacher. Note that this definition can be different from the common usage of overparametrisation by practitioners, which typically compares the number of parameters in the model class to the amount of available data. Indeed, since the complexity of the target function is well-defined in our teacher-student task, it is natural to define overparametrisation with respect to this complexity. The objective of the student network is to approximate the teacher function by minimising the discrepancy between their outputs. We define the population risk (or generalisation error) as the expected squared difference with respect to the input distribution: i 1 h 2 L(W ; W ∗ ) = Ex (ϕ(x, W ) − ϕ(x, W ∗ )) . (4) 2 The student learns by adjusting its parameters W to minimise L(W ; W ∗ ) via gradient flow: ẇk = −ηEx [Gk ] , ∗
Gk = (ϕ(x, W ) − ϕ(x, W )) H
wk ⊤ x √ d
x √ , d
(5)
where H(·) denotes the Heaviside step function. To analyse the system’s dynamics, we introduce the relevant order parameters (or sufficient statistics) in the form of weight overlaps: 1 1 1 ∗⊤ ∗ ∗ wn , Qij = wi⊤ wj , Rim = wi⊤ wm , Tmn = wm (6) d d d with Q ∈ RK×K , R ∈ RK×M , and T ∈ RM ×M . For analytical clarity, we assume orthonormal teacher weights, such that T = IM (the identity matrix). This assumption simplifies the discussion without loss of generality, with the same phenomenology as long as T is full-rank. Our analysis extends straightforwardly to arbitrary teacher configurations, with additional supporting results provided in Appendix C. 4
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Result 1 (Necessary Conditions for Minima). The stationary points of the population risk correspond to the zeros of the gradient flow. These satisfy a system of coupled non-linear equations governing the order parameters Q and R: FR (Q, R) = 0,
FQ (Q, R) = 0.
(7)
The explicit, analytical forms of the functionals FR and FQ are detailed in Appendix B. Notably, these equations are exact and independent of the input dimension d. We discuss the broader generality of these results in Sec.3.2. Result 2 (Classification of Local Minima). The local minima are not unique; they are organised into distinct families characterised by a structural parameter k1 ∈ [0, M ], representing the specific count of student units that become anti-aligned with the teacher vectors. By imposing a specific ansatz based on k1 , the general conditions in Result 1 reduce to a simplified set of equations (derived in Appendix B). By imposing a block-symmetric ansatz based on this symmetry-breaking arrangement, the general conditions in Result 1 reduce to a tractable set of equations. Where real solutions exist, this allows us to analytically determine the order parameters Q(k1 ) and R(k1 ) and directly compute macroscopic observables like the exact generalisation error. This structural parameter k1 directly maps to the transitivity partitions of the isotropy subgroups identified by Arjevani and Field [2019]. Our macroscopic reduction confirms that their discrete algebraic classes of spurious minima dictate the quantised levels of the population risk. Fig.1B. illustrates the geometric and dynamical consequences of these theoretical findings, serving as a visual guide to the results presented in Sec.3. Landscape Topology (Fig.1A.): The organisation of minima in the different families (Result 2) depends on the parameterisation regime and changes drastically as the student becomes over-parameterised. In the well-specified case, families consist of isolated, point-like minima separated by loss barriers. In the over-parameterised case, however, these barriers dissolve, and minima within a family merge to form continuous, connected manifolds. We explore this geometric transition in depth in Sec.3.1. Dynamical Selection (Fig.1B.): The existence of these families strongly affects the dynamics as discussed in Sec.3.2. The figure shows the density of solutions found by gradient flow across 104 initialisations. This highlights the sharp concentration of solutions around specific loss levels validates our theoretical predictions (vertical dashed lines). The insets display the structure of the order parameters Q and R for typical solutions, highlighting the symmetry-breaking patterns characteristic of each family.
3
Landscape geometry and consequences for optimisation
We leverage our theoretical framework to dissect the optimisation landscape. To do so, we rely on the numerical integration of the derived ODEs for the order parameters Q and R. This allows us to simulate the mean-field trajectories, which we use to complement, explain, and validate the behaviour of the simulations. Landscape Characterisation
Density
Density
200 100 00.00
0.01
0.02
0.03 Loss
0.04
A. K = 17, M = 17.
0.05
750 500 250 00.00
2000 Density
3.1
0.01
0.02
0.03 Loss
0.04
B. K = 18, M = 17.
0.05
1000 00.00
0.01
0.02
0.03 Loss
0.04
0.05
C. K = 19, M = 17.
Figure 2: Loss families and theoretical value. Histogram of final loss values obtained from ODE dynamics at M = 17 for K = M (left), K = M + 1 (center), and K = M + 2 (right) starting from 104 initialisations in orthonormal teacher configuration. Dashed vertical lines indicate the corresponding theoretical loss values obtained from the Result 2 For K = M and K = M + 1, the theoretical values capture the locations of the local minima families and the minor discrepancies may be attributed to numerical effects. In contrast, for K = M + 2, local minima are strongly suppressed in M = 17 and the loss distribution concentrates near zero. 5
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Spurious minima exist and are reachable. We start our investigation of the population loss landscape by the wellspecified case where both teacher and student have the same width (K = M ) — the minimal case for which the model has the capacity required to approximate the target. The evolution of the weights follows the gradient flow described in Eq. (5). By adopting the point of view of the order parameters, we can write a closed set of ordinary differential equations (ODEs) governing the time evolution of Q(t) and R(t), see Appendix A. In Fig.1b and Fig.2, we simulate these dynamics numerically starting from 104 independent random initialisations drawn from the standard normal distribution. Consistent with the empirical results by Safran and Shamir [2018], we observe that in the well-specified setting the gradient flow frequently becomes trapped in suboptimal stationary points with non-vanishing generalisation error. As shown in the Figures, the distribution of final loss values exhibits a quantised structure: steady states do not populate a continuum but rather concentrate tightly around discrete loss levels, indicating that the landscape is dominated by distinct families of spurious minima. To understand the geometric origin of these families, we examine the spectral properties of the order parameters. Due to the permutation symmetry of the hidden units, the network and loss are invariant under reordering of the student neurons. As exemplified in the insets of Fig.1b and discussed in detail in Appendix B.1, we observe that the discrete loss levels correspond to specific symmetry-breaking configurations where a subset of k1 student neurons becomes anti-aligned with the teacher vectors (i.e., Rim < 0), a structure made explicit in Fig. 8. This anti-alignment induces a large local error, which triggers a compensation mechanism in the remaining aligned neurons. Specifically, the aligned units adjust their orientation to offset the negative contribution of the anti-aligned units, effectively zeroing out the gradient. This cooperative effect stabilises the system in an isolated local minimum, preventing the anti-aligned units from flipping back to the correct orientation. As the number of anti-aligned units increases, the required compensation grows until the configuration becomes unstable; thus, these stable minima can be naturally ordered by the integer count of anti-aligned neurons, defining an hierarchy between them. This empirical regularity can be analytically formalised by imposing a block-symmetric ansatz on the order parameters (R, Q). We partition the student hidden units into two groups: a set I1 of size k1 (anti-aligned) and a set I2 of size k2 (aligned), such that k1 + k2 = K. We visualise this ansatz in Eq. (8), where we adopt the notation B(x, y) to represent a block matrix with diagonal elements x and off-diagonal elements y (i.e., xI + y(J − I)). cross cross B(r1diag , r1off ) r12 B(q1diag , q1off ) q12 R= , Q= . (8) cross cross r21 B(r2diag , r2off ) q12 B(q2diag , q2off ) Here, the blue blocks correspond to the k1 anti-aligned units, the green blocks to the k2 aligned units, and the grey blocks capture the interaction between the two populations. This specific block-symmetric structure of the order parameters serves as the macroscopic, statistical-mechanics equivalent of the discrete isotropy subgroups identified by Arjevani and Field [2019], confirming that the ”principle of least symmetry breaking” dictates the clustering of the network’s weights. In Appendix B, we report additional supporting validation of the ansatz. The Role of Over-parameterisation. Having characterised the spurious minima in the well-specified case, we turn our attention to how overparametrisation, i.e. increasing the student width to K > M , alters the landscape and facilitates convergence. Global Convergence and Connectivity. To quantify the benefit of overparameterisation, we track the optimisation trajectories of 103 random initialisations under gradient flow across three regimes: well-specified (K = M ) and mildly overparameterised (K = M + 1 and K = M + 2). As illustrated in Fig.2a, in the well-specified setting a significant fraction of trajectories remain trapped in the high-loss families identified in Sec.3.1. In contrast, the addition of even a single extra neuron (K = M + 1), Fig.2b-c, dramatically expands the basin of attraction of the global minimum, with nearly all trajectories converging to zero loss. These empirical results indicate that overparameterisation effectively destabilises the spurious minima that plague the K = M landscape. To understand the geometric origin of this improvement, we analyse the connectivity between solutions using the zero-temperature string method (see Appendix D for the precise definition). We first consider the K = M case. Due to permutation symmetry, there exist multiple discrete global minima. We construct a path between two such symmetric solutions: (R∗ , Q∗ ) and a permutation thereof. As shown by the orange curve in Fig.3A., the minimum loss path connecting them must traverse significant loss barriers, illustrating that in the well-specified regime, symmetric solutions are isolated in weight space. This behavior is consistent across all observed permutations. We then repeat the analysis for K = M + 1. We embed the K = M solutions into the larger space by appending a zero-row to the weight matrices (representing an inactive extra neuron). Connecting these embedded solutions via the string method reveals a qualitatively different geometry (purple curve, Fig.3A.). This embedding is the macroscopic analogue to the “neuron splitting” operation analysed by Fukumizu and Amari [2000], Safran et al. [2021]. While 6
0.05 0.04 0.03 0.02 0.01 0.00
distance from initial point
Loss
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
K=M K=M+1
0.0
0.2
0.4 0.6 0.8 String Coordinate
1.0
A. Loss along the string path.
100 10 1 10 2 10 3 10 4
K=17, M=17 K=18, M=17 K=19, M=17 10 6
10 4 10 2 perturbation
100
B. Perturbative analysis of fixed points.
Figure 3: Left: Loss along the string path. We show the loss along the string path connecting two two local minima with different weights in K = M case and K = M + 1 case obtained from the string method in GD. For K = M , the path crosses several different loss barriers, while for K = M + 1 the string remains at constant loss, indicating a flat direction connecting the minima. Right: Perturbative analysis of fixed points. Distance from the initial configuration obtained by adding a Gaussian perturbation of standard deviation σ in weight space. The plot consider the average distance after the network has been allowed to relax with gradient descent for 1, 000 steps and learning rate 0.01. The mean (solid line) and standard deviation (shaded area) are evaluated by collecting 20 independent perturbations. The different line represent different level of over-parameterisation. The exact-parametirisation case appears to be less than 0.001 from the original configuration, a distance that is compatible with the limited epochs allowed by gradient descent, and it is therefore considered stable from small enough σ. On contrary, the overparameteritdsed case show instability even at small σ.
Safran et al. [2021] showed this operation induces a negative Hessian eigenvalue, our string method approach reveals the global consequence: the loss barrier vanishes entirely, and the loss remains approximately constant along the path. This implies that the additional degree of freedom introduces a flat direction—a continuous manifold of zero-loss solutions—that connects previously isolated minima, thereby allowing the optimiser to traverse the landscape without overcoming potential barriers. While the string of fixed points originated by adding a zero-weight neuron connects the landscape, it is crucial to determine if it is a stable attractor or not. Standard Hessian analysis is ill-suited here due to the non-differentiability of the ReLU activation. Instead, we rely on a perturbative dynamical analysis. We initialise the system at a fixed point and apply a Gaussian perturbation to the weights: W → W + ξ, where ξ ∼ N (0, σ 2 I). Note that we perturb the weights rather than the order parameters directly, as arbitrary perturbations to Q and R may violate the geometric constraints (details in Appendix E). We then allow the system to relax via gradient flow. The results, reported in Fig.3B., demonstrate that, while in the well-specified case spurious fixed points are stable to perturbations, i.e. the dynamics consistently return to the attractor. Conversely, in the over-parameterised case (K ≥ M + 1), these same fixed points become unstable. The perturbation excites the extra neuron, allowing the system to escape the saddle point. Landscape Geometry in the Overparameterised Regime We now extend our topological analysis to the overparameterised setting (K = M + 1). We generalise the ansatz (Result 2) to accommodate the additional degrees of freedom, identifying stationary points where the extra neurons carry non-trivial weight (details in Appendix B). This revealed that in the overparameterised landscape, the lowest-order spurious minima are destabilised. Specifically, for the family characterised by a single anti-aligned neuron (k1 = 1), we are not able to find stable solutions to the fixed-point equations. This theoretical prediction is corroborated by numerical integration of the ODEs reported in Table 1: across all initialisations, no trajectory converged to a state exhibiting a single anti-aligned direction. However, higher-order spurious minima (with k1 ≥ 2) persist. These families, which correspond to more complex symmetrybreaking configurations, remain stable even in the presence of the additional neuron. This paints a more complex picture with respect to Safran et al. [2021]. Indeed, while their results established that overparameterisation can destabilise local minima via simple neuron splitting (provided the weight norms are bounded—a condition naturally satisfied by P our macroscopic order parameters, where i ∥wi ∥ ≈ M − k1 ), our global analysis reveals a much more intricate landscape. We find that while simple, uncoupled defects (k1 = 1) are indeed annihilated by the addition of a single neuron, higher-order spurious minima (k1 ≥ 2) actively persist. Crucially, these surviving local minima are not merely artifacts of adding a zero-weight neuron to an exact-parameterisation solution; they are complex, coupled structures that emerge uniquely in the overparameterised space. Because they cannot be derived through simple zero-padding or 7
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
k1 = 5
Loss
0.030 0.025
k1 = 4
0.020
k1 = 3
0.015 k1 = 2
0.0
0.2
0.4 0.6 0.8 String Coordinate
1.0
Figure 4: Connectivity of solution manifolds in the over-parameterised regime (K = M + 1). Loss profiles along the minimum loss paths (computed via the string method) connecting two symmetric realisations of the same local minimum. Different colorus correspond to different families indexed by k1 (the number of anti-aligned units). The paths are perfectly flat, indicating that the isolated fixed points of the well-specified regime have expanded into continuous, connected manifolds.
local Hessian perturbations, identifying these surviving traps fundamentally requires solving the exact macroscopic fixed-point equations. M inima Order k1 = 0 k1 = 1 k1 = 2 k1 = 3 k1 = 4
K = 17 M = 17 13.09% 27.52% 29.05% 18.94% 7.55%
K = 18 M = 17 59.29% 0.00% 2.10% 10.83% 8.99%
K = 19 M = 17 99.63% 0.00% 0.05% 0.31% 0%
Table 1: Percentage of 104 random initialisations that converge to minima with k1 anti-aligned units in different (K, M ); k1 = 0 corresponds to the global minimum. While high-loss minima still exist for k1 ≥ 2, their geometry is fundamentally altered. In the well-specified case, distinct permutations of the same solution were isolated by loss barriers. To test if this isolation persists in the overparameterised regime, we apply the string method to find the minimum loss path connecting two symmetric configurations within the same k1 family. The results, illustrated in Fig.4, reveal a topological phase transition. For all surviving families, the computed paths exhibit constant loss profiles (flat strings). This indicates that the loss barriers separating symmetric solutions have vanished. Consequently, the discrete, isolated minima observed at K = M have merged into continuous, flat manifolds. While additional phenomenons appear for k1 ≥ 2, they do not modify this conclusion; for completeness, all string method results for the case k1 = 2 are reported in Appendix F. 3.2
Dynamical Characterisation
As discussed in the problem formulation Sec.2, the evolution of the weights is governed by a closed set of ODEs for the order parameters, Q̇ = FQ (Q, R) and Ṙ = FR (Q, R). While the full explicit forms are deferred to Appendix A, the dynamics can be written in terms of Gaussian expectation terms whose form depends on the activation function. For ReLU, these admit closed-form expressions, given in Appendix A.4; corresponding derivations for Leaky ReLU and sigmoidal (erf) activations are provided in Appendix G. In Fig.5, we numerically integrate the gradient flow ODEs to validate our theoretical landscape analysis. The simulations confirm that the minima identified with our ansatz are important attractors of the problem. Indeed, the loss trajectories do not settle into a continuum, but rather concentrate into discrete, quantised bands. These plateaus align precisely with the loss levels of the spurious families predicted by Result 2. Furthermore, comparing the regimes (K = M vs. K > M ) illustrates the dynamical advantage of overparameterisation, which destabilises some of the intermediate bands and facilitates the cascade toward the global minimum. 8
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
100
M=20, K=20 M=20, K=21 M=20, K=22
Loss
10 1
10 2
10 3
0.0
0.5
1.0
1.5
Steps
2.0
2.5
3.0
1e6
Figure 5: Training dynamics and loss quantisation across parameterisation regimes. Evolution of population risk under gradient flow for M = 20 teacher units (1,000 random initialisations per condition, learning rate η = 0.1 , orthonormal teacher configuration). All regimes exhibit entrapment in discrete high-loss plateaux. However, while the well-specified case (K=20, orange) is dominated by these suboptimal attractors, mild over-parameterisation (K = 21, purple; K = 22, blue) progressively destabilises intermediate plateaus, allowing a larger fraction of trajectories to escape towards the global minimum. The observed loss bands correspond to the theoretical loss levels of the spurious families predicted by Result 2. Alternative Optimisation Schemes. While our primary analysis focuses on gradient flow, our framework is flexible and can be readily adapted to alternative descent-based algorithms often studied in theoretical literature or employed in practice. We derive the corresponding modified ODEs for the following settings in Appendix A and Appendix H: GD (nGD). Optimisation is performed on the hypersphere with a fixed norm constraint, ∥wk ∥2 = • Normalized √ d. The dynamics project the gradient onto the tangent space of Sd−1 : η ẇk = −ηE[Gk ] + ⟨E[Gk ], wk ⟩wk . (9) d • Orthonormalized GD (onGD). The weights are constrained to the Stiefel manifold, ensuring the student weight matrix remains orthonormal (W W ⊤ = IK ). The update rule projects the gradient via PW (Z) = Z − d1 ZW ⊤ W : Ẇ = −ηPW (E[G]). (10) • Two-Layer GD (2L-GD). Both student weights wk and readout coefficients vk are trained simultaneously. The continuous-time dynamics are coupled: ẇk = −ηE[Gkw ] and v̇k = −ηE[Gkv ]. (11) Consequently, the flow equations must be augmented to include the evolution of the second layer, i.e., Q̇ = FQ (Q, R, v). • Online Stochastic Gradient Descent (SGD). Optimisation is performed on the empirical loss using new samples at each step: ∆wk = −ηGk . While formally a discrete process, in the limit of high-dimensional inputs (d → ∞), the dynamics can be described by a deterministic set of ODEs (see Appendix A.6). A key question concerns the generality of the landscape properties derived under gradient flow. Result 3 (Landscape Equivalence for SGD). If the learning rate scales as η = od (1), the additional diffusive terms in the ODEs for SGD vanish (see Eqs. 30-32 Appendix A.6). In this regime, the trajectory of SGD converges to that of the gradient flow. Consequently, the topological analysis of the landscape presented in previous sections applies directly to SGD in the appropriate scaling limit. Empirical Comparison and Convergence. We compare the convergence properties of these optimisers in Table 2, aggregating statistics from 104 initialisations. The results highlight a sharp distinction between unconstrained (GD, 2L-GD) and constrained (nGD, onGD) dynamics. First, constrained dynamics exhibit significantly slower convergence timescales, often failing to reach steady states even after extended training (denoted as N/C in Table 2). To provide a complete picture of these failure modes, Appendix I.1 provides the full classification tables for nGD across all network widths, alongside the unquantised final loss distributions characteristic of onGD. Second, the nature of the landscape for onGD is fundamentally different. The compensation mechanism identified in Sec.3.1 relies on the ability of aligned neurons to adjust their magnitude to cancel the error from anti-aligned units. The orthonormality constraint forbids such magnitude adjustments. Therefore, 9
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Optimiser GD 2L-GD nGD onGD
K = 17 M = 17 13.25% 13.24% 14.12% N/C
K = 18 M = 17 64.18% 67.91% 58.35% N/C
K = 19 M = 17 77.50% 99.48% N/C N/C
Table 2: Convergence frequency across optimisation schemes. The table reports the percentage of runs reaching the global minimum out of 104 initialisations. Columns represent different parameterisation regimes (exact vs. overparameterised). Consistent with Safran and Shamir [2018], over-parameterisation aids convergence in unconstrained dynamics. A more thorough comparison is reported in Appendix I. However, constrained dynamics (nGD, onGD) show severe slowdowns, resulting in non-convergence (N/C) within the allocated computational budget (1.2 × 107 steps), a comprehensive analysis is detailed in Appendix I.1. the specific families of spurious minima discussed previously cannot exist in the onGD setting, rendering its landscape distinct from the standard ReLU network landscape. Conversely, nGD recovers a minima structure analogous to GD. However, in the overparameterised regime (K = M +2), we observe a high frequency of “mixed” states characterised by combinations of minima from different families. While these states are theoretically unstable, their gradients suggest proximity to saddle points Zhang et al. [2025], causing severe dynamical slowdowns. This accounts for the lower convergence rates compared to GD and 2L-GD.
4
Conclusions
In this work, we provided a sharp characterisation of the population loss landscape for high-dimensional two-layer ReLU networks, connecting perspectives from dynamical approaches and landscape analysis and shedding new light on open problems left by Safran and Shamir [2018]. Our analysis reveals a fundamental topological phase transition driven by over-parameterisation. In the well-specified regime, the landscape is dominated by a hierarchy of isolated spurious minima. Conversely, over-parameterisation destabilises low-order minima and merges families of high-order ones into connected manifolds. This geometric mechanism—the emergence of flat directions—allows the optimiser to circumvent loss barriers, providing a principled explanation for the trainability of wide networks. Leveraging mean-field theory, we confirmed that these landscape features act as strong attractors for the gradient flow, governing the “quantised” loss distributions observed empirically. We further established the generality of these findings, extending the validity of this landscape description to SGD under appropriate scaling. By doing so, our macroscopic framework bridges two distinct perspectives in the literature: the group-theoretic symmetries of spurious minima [Arjevani and Field, 2019] and the local Hessian instabilities triggered by neuron splitting [Safran et al., 2021]. We show that these static, finite-width geometric properties govern the global, high-dimensional learning trajectories. Key questions remain for future investigation. While we have established the connectivity of solutions, the precise measure of the basins of attraction for different initialisations remains to be quantified. Additionally, extending this analysis to a broader class of non-linearities and optimisation protocols offers a promising future direction. Overall, this work establishes a rigorous, interpretable framework for understanding how over-parameterisation mitigates the inherent non-convexity of neural network training.
Acknowledgments The authors thank Ashkan Panahi for precious feedback. BL would like to acknowledge discussion with Berfin Simsek, Francesca Mignacco, Ludovic Stephan, Sebastian Goldt and Courtney Paquette concerning this problem. S.S.M. was supported by the Wallenberg AI, AutonomousvSystems, and Software Program (WASP). B.L. was supported by the French government, managed by the National Research Agency (ANR), under the France 2030 program with the project references “ANR-23-IACL-0008” (PR[AI]RIE-PSAI) and “ANR-25-CE23-5660” (MAPLE), as well as the Choose France - CNRS AI Rising Talents program. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at Alvis (Chalmers Centre for Computational Science and Engineering, C3SE) partially funded by the Swedish Research Council through grant agreement no. 2022-06725, under project NAISS 2024/22-1082.
References Lenaic Chizat and Francis Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. Advances in neural information processing systems, 31, 2018. 10
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Song Mei, Andrea Montanari, and Phan-Minh Nguyen. A mean field view of the landscape of two-layer neural networks. Proceedings of the National Academy of Sciences, 115(33):E7665–E7671, 2018. Grant Rotskoff and Eric Vanden-Eijnden. Trainability and accuracy of artificial neural networks: An interacting particle system approach. Communications on Pure and Applied Mathematics, 75(9):1889–1935, 2022. Justin Sirignano and Konstantinos Spiliopoulos. Mean field analysis of neural networks: A law of large numbers. SIAM Journal on Applied Mathematics, 80(2):725–752, 2020. Francis Bach and Lénaı̈c Chizat. Gradient descent on infinitely wide neural networks: Global convergence and generalization. arXiv preprint arXiv:2110.08084, 2021. Itay Safran and Ohad Shamir. Spurious local minima are common in two-layer relu neural networks. In International conference on machine learning, pages 4433–4441. PMLR, 2018. Itay M Safran, Gilad Yehudai, and Ohad Shamir. The effects of mild over-parameterization on the optimization landscape of shallow relu neural networks. In Conference on Learning Theory, pages 3889–3934. PMLR, 2021. E Weinan, Weiqing Ren, and Eric Vanden-Eijnden. String method for the study of rare events. Physical Review B, 66 (5):052301, 2002. Weiqing Ren, Eric Vanden-Eijnden, et al. Simplified and improved string method for computing the minimum energy paths in barrier-crossing events. The Journal of chemical physics, 126(16), 2007. Michael Biehl and Holm Schwarze. Learning by on-line gradient descent. Journal of Physics A: Mathematical and general, 28(3):643, 1995. David Saad and Sara Solla. Dynamics of on-line gradient descent learning for multilayer neural networks. Advances in neural information processing systems, 8, 1995. Sebastian Goldt, Madhu Advani, Andrew M Saxe, Florent Krzakala, and Lenka Zdeborová. Dynamics of stochastic gradient descent for two-layer neural networks in the teacher-student setup. Advances in neural information processing systems, 32, 2019. Sebastian Goldt, Marc Mézard, Florent Krzakala, and Lenka Zdeborová. Modeling the influence of data structure on learning in neural networks: The hidden manifold model. Physical Review X, 10(4):041044, 2020. Gerard Ben Arous, Reza Gheissari, and Aukosh Jagannath. High-dimensional limit theorems for sgd: Effective dynamics and critical scaling. Advances in neural information processing systems, 35:25349–25362, 2022. Rodrigo Veiga, Ludovic Stephan, Bruno Loureiro, Florent Krzakala, and Lenka Zdeborová. Phase diagram of stochastic gradient descent in high-dimensional two-layer neural networks. Advances in Neural Information Processing Systems, 35:23244–23255, 2022. Luca Arnaboldi, Yatin Dandi, Florent Krzakala, Bruno Loureiro, Luca Pesce, and Ludovic Stephan. Online learning and information exponents: The importance of batch size Time/Complexity tradeoffs. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp, editors, Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 1730–1762. PMLR, 21–27 Jul 2024. Luca Arnaboldi, Ludovic Stephan, Florent Krzakala, and Bruno Loureiro. From high-dimensional & mean-field dynamics to dimensionless odes: A unifying approach to sgd in two-layers networks. In The Thirty Sixth Annual Conference on Learning Theory, pages 1199–1227. PMLR, 2023a. Luca Venturi, Afonso S Bandeira, and Joan Bruna. Spurious valleys in two-layer neural network optimization landscapes. arXiv preprint arXiv:1802.06384, 2018. Felix Draxler, Kambis Veschgini, Manfred Salmhofer, and Fred Hamprecht. Essentially no barriers in neural network energy landscape. In International conference on machine learning, pages 1309–1318. PMLR, 2018. Berfin Simsek, François Ged, Arthur Jacot, Francesco Spadaro, Clément Hongler, Wulfram Gerstner, and Johanni Brea. Geometry of the loss landscape in overparameterized neural networks: Symmetries and invariances. In International Conference on Machine Learning, pages 9722–9732. PMLR, 2021. Stefano Sarao Mannelli, Eric Vanden-Eijnden, and Lenka Zdeborová. Optimization and generalization of shallow neural networks with quadratic activation functions. Advances in Neural Information Processing Systems, 33:13445–13455, 2020a. Damek Davis, Dmitriy Drusvyatskiy, and Courtney Paquette. The nonsmooth landscape of phase retrieval. IMA Journal of Numerical Analysis, 40(4):2652–2695, 2020. Kenji Fukumizu and Shun-ichi Amari. Local minima and plateaus in hierarchical structures of multilayer perceptrons. Neural networks, 13(3):317–327, 2000. 11
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
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. In Conference on Learning Theory, pages 4782–4887. PMLR, 2022. Anchit Jain, Rozhin Nobahari, Aristide Baratin, and Stefano Sarao Mannelli. Bias in motion: Theoretical insights into the dynamics of bias in sgd training. Advances in Neural Information Processing Systems, 37:24435–24471, 2024. Raphaël Berthier, Andrea Montanari, and Kangjie Zhou. Learning time-scales in two-layers neural networks. Foundations of Computational Mathematics, 25(5):1627–1710, 2025. Andrea Montanari and Pierfrancesco Urbani. Dynamical decoupling of generalization and overfitting in large two-layer networks. arXiv preprint arXiv:2502.21269, 2025. Yossi Arjevani and Michael Field. On the principle of least symmetry breaking in shallow relu models. arXiv preprint arXiv:1912.11939, 2019. Johanni Brea, Berfin Simsek, Bernd Illing, and Wulfram Gerstner. Weight-space symmetry in deep networks gives rise to permutation saddles, connected by equal-loss valleys across the loss landscape. arXiv preprint arXiv:1907.02911, 2019. Carlo Baldassi, Enrico M Malatesta, and Riccardo Zecchina. Properties of the geometry of solutions and capacity of multilayer neural networks with rectified linear unit activations. Physical review letters, 123(17):170602, 2019. Carlo Baldassi, Clarissa Lauditi, Enrico M Malatesta, Gabriele Perugini, and Riccardo Zecchina. Unveiling the structure of wide flat minima in neural networks. Physical Review Letters, 127(27):278301, 2021. Brandon Livio Annesi, Clarissa Lauditi, Carlo Lucibello, Enrico M Malatesta, Gabriele Perugini, Fabrizio Pittorino, and Luca Saglietti. Star-shaped space of solutions of the spherical negative perceptron. Physical Review Letters, 131 (22):227301, 2023. Damien Barbier, Ahmed El Alaoui, Florent Krzakala, and Lenka Zdeborová. On the atypical solutions of the symmetric binary perceptron. Journal of Physics A: Mathematical and Theoretical, 57(19):195202, 2024. Damien Barbier. How to escape atypical regions in the symmetric binary perceptron: a journey through connectedsolutions states. SciPost Physics, 18(3):115, 2025. Amit Samanta and E Weinan. Optimization-based string method for finding minimum energy path. Communications in Computational Physics, 14(2):265–275, 2013. Carlo Baldassi, Fabrizio Pittorino, and Riccardo Zecchina. Shaping the learning landscape in neural networks around wide flat minima. Proceedings of the National Academy of Sciences, 117(1):161–170, 2020. Yedi Zhang, Andrew Saxe, and Peter E Latham. Saddle-to-saddle dynamics explains a simplicity bias across neural network architectures. arXiv preprint arXiv:2512.20607, 2025. Leticia F Cugliandolo and Jorge Kurchan. Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model. Physical Review Letters, 71(1):173, 1993. Leticia F Cugliandolo and Jorge Kurchan. On the out-of-equilibrium relaxation of the sherrington-kirkpatrick model. Journal of Physics A: Mathematical and General, 27(17):5749, 1994. 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(106):1–51, 2021. Alberto Bietti, Joan Bruna, and Loucas Pillaud-Vivien. On learning gaussian multi-index models with gradient flow part i: General properties and two-timescale learning. Communications on Pure and Applied Mathematics, 78(12): 2354–2435, 2025. Berfin Şimşek, Amire Bendjeddou, and Daniel Hsu. Learning gaussian multi-index models with gradient flow: Time complexity and directional convergence. arXiv preprint arXiv:2411.08798, 2024. Antoine Maillard, Gérard Ben Arous, and Giulio Biroli. Landscape complexity for the empirical risk of generalized linear models. In Mathematical and Scientific Machine Learning, pages 287–327. PMLR, 2020. Kiana Asgari, Andrea Montanari, and Basil Saeed. Local minima of the empirical risk in high dimension: General theorems and convex examples. arXiv preprint arXiv:2502.01953, 2025. Andrew M. Saxe, James L. McClelland, and Surya Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. In Yoshua Bengio and Yann LeCun, editors, 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014. URL http://arxiv.org/abs/1312.6120. Andrew M Saxe, James L McClelland, and Surya Ganguli. A mathematical theory of semantic development in deep neural networks. Proceedings of the National Academy of Sciences, 116(23):11537–11546, 2019. 12
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Etienne Boursier, Loucas Pillaud-Vivien, and Nicolas Flammarion. Gradient flow dynamics of shallow relu networks for square loss and orthogonal inputs. Advances in Neural Information Processing Systems, 35:20105–20118, 2022. Andrew Saxe, Shagun Sodhani, and Sam Jay Lewallen. The neural race reduction: Dynamics of abstraction in gated networks. In International Conference on Machine Learning, pages 19287–19309. PMLR, 2022. Devon Jarvis, Sebastian Lee, Clémentine Carla Juliette Dominé, Andrew M. Saxe, and Stefano Sarao Mannelli. A theory of initialisation’s impact on specialisation. In The Thirteenth International Conference on Learning Representations, ICLR 2025, Singapore, April 24-28, 2025. OpenReview.net, 2025. URL https://openreview.net/forum?id= RQz7szbVDs. Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma. Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval. Mathematical Programming, 176(1):5–37, 2019. Stefano Sarao Mannelli, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborova. Passed & spurious: Descent algorithms and local minima in spiked matrix-tensor models. In international conference on machine learning, pages 4333–4342. PMLR, 2019. Luca Arnaboldi, Florent Krzakala, Bruno Loureiro, and Ludovic Stephan. Escaping mediocrity: how two-layer networks learn hard generalized linear models with sgd. arXiv preprint arXiv:2305.18502, 2023b. Stefano Sarao Mannelli, Giulio Biroli, Chiara Cammarota, Florent Krzakala, and Lenka Zdeborová. Who is afraid of big bad minima? analysis of gradient-flow in spiked matrix-tensor models. Advances in neural information processing systems, 32, 2019. Stefano Sarao Mannelli, Giulio Biroli, Chiara Cammarota, Florent Krzakala, Pierfrancesco Urbani, and Lenka Zdeborová. Marvels and pitfalls of the langevin algorithm in noisy high-dimensional inference. Physical Review X, 10 (1):011057, 2020b. Tony Bonnaire, Giulio Biroli, and Chiara Cammarota. From zero to hero: How local curvature at artless initial conditions leads away from bad minima. arXiv preprint arXiv:2403.02418, 2024. Tony Bonnaire, Giulio Biroli, and Chiara Cammarota. The role of the time-dependent hessian in high-dimensional optimization. Journal of Statistical Mechanics: Theory and Experiment, 2025(8):083401, 2025. Brandon Livio Annesi, Dario Bocchi, and Chiara Cammarota. Overparametrization bends the landscape: Bbp transitions at initialization in simple neural networks. arXiv preprint arXiv:2510.18435, 2025. Marco Baity-Jesi, Levent Sagun, Mario Geiger, Stefano Spigler, Gérard Ben Arous, Chiara Cammarota, Yann LeCun, Matthieu Wyart, and Giulio Biroli. Comparing dynamics: Deep neural networks versus glassy systems. In International Conference on Machine Learning, pages 314–323. PMLR, 2018. Samantha J Fournier, Alessandro Pacco, Valentina Ros, and Pierfrancesco Urbani. Non-reciprocal interactions and high-dimensional chaos: comparing dynamics and statistics of equilibria in a solvable model. arXiv preprint arXiv:2503.20908, 2025. Valentina Ros, Gerard Ben Arous, Giulio Biroli, and Chiara Cammarota. Complex energy landscapes in spiked-tensor and simple glassy models: Ruggedness, arrangements of local minima, and phase transitions. Physical Review X, 9 (1):011003, 2019. Alessandro Pacco, Alberto Rosso, and Valentina Ros. Triplets of local minima in a high-dimensional random landscape: correlations, clustering, and memoryless activated jumps. Journal of Statistical Mechanics: Theory and Experiment, 2025(3):033302, 2025.
Appendix This appendix provides detailed derivations, algorithmic descriptions, and additional experimental results supporting the Main Text. The content is organised as follows: • Appendix A: Dynamics. Full derivation of the mean-field ODEs for Gradient Descent, the calculation of relevant Gaussian integrals (I2 , I3 ), and the extension to Stochastic Gradient Descent. • Appendix B: Equations of Fixed Points. Derivation of the implicit equations for stationary points, the block-symmetric ansatz formulation, and numerical verification of the solutions. 13
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
• Appendix C: Robustness to Finite Input Dimensionality. Analysis of the impact of finite input dimensions (d) on the loss distribution and the structural robustness of the macroscopic order parameters. • Appendix D: String Method. Details on the zero-temperature string method used to probe connectivity, including the definition of the induced metric in order-parameter space and the path-finding algorithm. • Appendix E: Stability of fixed points from higher order minima. Perturbative dynamical analysis examining the stability of fixed points and the escape mechanisms from saddle points. • Appendix F: Additional String Method Results in the Overparameterised Regime. Further empirical examples of connectivity and minimum energy paths connecting symmetric solutions when K > M. • Appendix G: Alternative Activation Functions. Analytical derivation of integrals and landscape analysis for Leaky ReLU and Sigmoidal (erf) activations. • Appendix H: Constraint Dynamics. Derivations of the mean-field equations for Normalized GD (nGD) and Orthonormalised GD (onGD), along with extended empirical results (Appendix I.1) analysing their convergence properties and unquantised loss distributions. • Appendix I: Spurious minima in two-layer ReLU networks. Tabulated statistics of global versus local minima across different parameterisation regimes for unconstrained gradient descent, validating agreement with prior empirical works.
A
Dynamics
As shown in Fig. 6, the mean-field ODEs provide an accurate description of the loss dynamics observed in simulations of finite-width neural networks trained by gradient descent. We now provide the derivation of these equations. 100
1.0
Simulated ODEs
0.8
10 1
1.0
R[0,0] R[0,1] R[1,0] R[1,1]
0.8 0.6 Value
Loss
Value
0.6
10 2
0.4
0.5
1.0
1.5 Steps
A. Training loss
2.0
2.5 1e5
0.4
0.2
0.2
0.0
0.0
10 3 0.0
Q[0,0] Q[0,1] Q[1,1]
0.0
0.5
1.0
1.5 Steps
2.0
2.5 1e5
B. Entries of the order parameter R
0.0
0.5
1.0
1.5 Steps
2.0
2.5 1e5
C. Entries of the order parameter Q
Figure 6: Comparison between mean-field ODE predictions and network simulations. Panel (a) shows the training loss, while panels (b) and (c) display representative entries of the order parameters R and Q as functions of training steps. Solid lines denote simulation results and markers indicate ODE predictions. A.1
Population Gradient
We provide additional details on the ODEs for our learning problems. The network is trained using Gradient Descent (GD), Normalized GD (nGD), or Orthonormalized GD (onGD) on the MSE loss function. In this section, we focus on 14
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
the derivation for standard GD with ReLU activation functions. The derivations for nGD and onGD are detailed in Appendix H, and integrals for other activations (e.g., Leaky ReLU, erf) are provided in Appendix G. Following the notation used by Goldt et al. [2019], we use indices i, j, k, l = 1, . . . , K for student units and n, m = 1, . . . , M for teacher units. For the sake of completeness, we discussed the general two-layer network case, were both first layer (W ) and second layer (V ) are trained. Eventually, we will ignore the update of the second layer for the case where this is relevant. Additionally, we will consider a generic activation function g and specialise on individual activation differences later on. We express the loss function between student and teacher networks as: *" K #2 + M X X 1 ∗ L= vk g(λk ) − vm g(ρm ) , 2 m=1
(12)
k=1
w⊤ x
w∗ ⊤ x
√ where λk = √kd and ρm = m are the pre-activations. The continuous-time dynamics of the student weights wi d under gradient descent are given by: x dwi ′ = −η vi ∇wi L = −η vi ∆ g (λi ) √ , (13) dt d dvi = −η∇vi L = −η ⟨∆ g(λi )⟩ , (14) dt P P ∗ where ∆ = ϕ(x) − ϕ∗ (x) = j vj g(λj ) − m vm g(ρm ) is the error signal, and η is the learning rate for weights.
A.2
Derivation of Dynamics
The pre-activations λ and ρ are joint multivariate normal distributions with zero mean. The covariance matrix Σ is determined by the order parameters already reported in Eqs.6 of the Main Text: Qik =
wi⊤ wk = ⟨λi λk ⟩, d
Rin =
wi⊤ wn∗ = ⟨λi ρn ⟩, d
Tnm =
∗ wn∗ ⊤ wm = ⟨ρn ρm ⟩. d
To derive the closed-form ODEs, we apply the chain rule to the definition of the order parameters. For the student-student overlap Qik = d1 wi⊤ wk : ! 1 dwi ⊤ dw dQik k = wk + wi⊤ . (15) dt d dt dt Substituting the weight dynamics from Eq. (13) into the first term, we obtain: 1 dwi ⊤ η x⊤ wk η wk = − vi ∆ · g ′ (λi ) √ = − vi ⟨∆ · g ′ (λi )λk ⟩ . d dt d d d
(16)
The remaining factor 1/d represents the intrinsic time scale of the mean-field dynamics. By absorbing this factor into the continuous time variable (i.e., defining the macroscopic time step dt ∼ d1 ), we obtain the O(1) dynamical equations: dQik = −ηvi ⟨∆g ′ (λi )λk ⟩ − ηvk ⟨∆g ′ (λk )λi ⟩, dt
(17)
where η denotes the effective learning rate 2 . The second term follows by symmetry (i ↔ k). Expanding the error signal ∆, we arrive at terms of the form ⟨g ′ (λi )λk g(λj )⟩, which motivates the definition of the integral I3 . This formulation applies generally to two-layer networks with activation function g(·) and hidden layer width not larger than O(d). See Arnaboldi et al. [2023a] for a detailed discussion. Similarly, for the student-teacher overlap Rin : √ dRin 1 dwi ⊤ ∗ η = wn = − vi ⟨∆ · g ′ (λi ) dρn ⟩ =⇒ −ηvi ⟨∆ · g ′ (λi )ρn ⟩. dt d dt d 2
(18)
Alternatively, one can absorb the entire η/d term into the differential time, showing that it is enough to have η or 1/d to obtain the ODEs. More discussion on this can be found in Arnaboldi et al. [2023a].
15
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Finally, for the second-layer weights vi , the gradient descent dynamics are given by: dvi = −η ⟨∆ · g(λi )⟩ . dt PK PM ∗ Substituting the explicit form of the error signal ∆ = j=1 vj g(λj ) − m=1 vm g(ρm ), we obtain: K M X X dvi ∗ = −η vj ⟨g(λj )g(λi )⟩ − vm ⟨g(ρm )g(λi )⟩ . dt m=1 j=1
(19)
(20)
Here, the expectations are taken over the joint Gaussian distribution of the fields. We observe terms involving the product of two activation functions, ⟨g(·)g(·)⟩, which motivates the definition of the two-variable integral I2 (x1 , x2 ) ≡ ⟨g(x1 )g(x2 )⟩. A.3
Full Equations
The closed-form ODEs training for two layers are: M K X X dRin ∗ = ηvi vm I3 (i, n, m) − vj I3 (i, n, j) , dt m=1 j=1 M K M K X X X X dQik ∗ ∗ = ηvi vm I3 (i, k, m) − vj I3 (i, k, j) + ηvk vm I3 (k, i, m) − vj I3 (k, i, j) , dt m=1 m=1 j=1 j=1 K M X X dvi ∗ vj I2 (i, j) . vm I2 (i, m) − = η dt m=1 j=1
(21)
(22)
(23)
where we slightly abuse notation in the arguments of I2 , I3 to explicitly show which variables are involved. For example, I3 (i, n, m) implies C00 = Qii , C11 = Tnn , C22 = Tmm , C01 = Rin , etc. A.4
ReLU Activation Function Case
For ReLU, the Gaussian integrals have closed-form expressions. We collect here the explicit formulas used in the dynamics; corresponding expressions for other activations are given in Appendix G. Let (x0 , x1 , x2 ) be jointly Gaussian variables with zero mean and covariance matrix Σ, where Cij = ⟨xi xj ⟩. Two-Variable Integral (I2 ).
Used for the dynamics of the second-layer weights v. I2 (x1 , x2 ) = ⟨g(x1 )g(x2 )⟩.
For ReLU, explicit forms are given by: p I2 (Σ) =
2 C11 C22 − C12 C12 C12 + arccos − √ . 2π 2π C11 C22
(24)
(25)
Three-Variable Integral (I3 ). Used for the dynamics of Q and R. We define I3 specifically as the expectation involving a derivative, a multiplier, and an activation input: I3 (x0 , x1 , x2 ) = ⟨g ′ (x0 ) x1 g(x2 )⟩. For ReLU, g ′ (x) = H(x) = 1x>0 (Heaviside step function). The analytic solution is: q 1 C02 C01 2 I3 (Σ) = C12 arccos − √ C22 C00 − C02 , + 2π C00 C22 C00 In the full equations below, the notation I3 (i, k, j) maps to x0 = λi , x1 = λk , and x2 = λj . 16
(26)
(27)
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
A.5
Loss Function in Terms of Order Parameters
Finally, utilising the definition of the two-variable integral I2 , we can express the population loss L in a compact closed form. Expanding the squared error term and exchanging the summation and expectation operations, we obtain: * K !2 + M X X 1 ∗ L= vk g(λk ) − vm g(ρm ) 2 m=1 k=1
=
K X
M K X M X 1 X ∗ ∗ 1 ∗ vi vj ⟨g(λi )g(λj )⟩ + vn vm ⟨g(ρn )g(ρm )⟩ − vk vm ⟨g(λk )g(ρm )⟩. 2 i,j=1 2 n,m=1 m=1
(28)
k=1
Substituting the integral definition I2 (x, y) = ⟨g(x)g(y)⟩, the loss function becomes: L=
K M K X M X 1 X ∗ ∗ 1 X ∗ vi vj I2 (i, j) + vn vm I2 (n, m) − vk v m I2 (k, m). 2 i,j=1 2 n,m=1 m=1
(29)
k=1
A.6
Stochastic Gradient Descent
In the case of Stochastic Gradient Descent, as shown in Saad and Solla [1995], Biehl and Schwarze [1995], the ODE for Q presents an additional term resulting from the stochasticity of the update. This results in the equations below: M K X X dRin ∗ = ηvi vm I3 (i, n, m) − vj I3 (i, n, j) , (30) dt m=1 j=1 M K M K X X X X dQik ∗ ∗ = ηvi vm I3 (i, k, m) − vj I3 (i, k, j) + ηvk vm I3 (k, i, m) − vj I3 (k, i, j) dt m=1 m=1 j=1 j=1 # " M K M K X X X X ∗ ∗ ∗ vj vl I4 (i, k, j, l) , (31) + η 2 vi vk vm vj vm I4 (i, k, m, j) + vn I4 (i, k, n, m) − 2 m,n=1
M X
j=1 m=1 K X
j,l=1
dvi = η vn∗ I2 (i, n) − vj I2 (i, j) . dt n=1 j=1
(32)
Where the new term I4 is given by I4 (x0 , x1 , x2 , x3 ) = ⟨g ′ (x0 ) g ′ (x1 ) g(x2 ) g(x3 )⟩.
B
Equations of Fixed Points
The implicit equations described in Result 1 describing the fixed points are the following: 0 = FR (Q, R) =
M K X X dRin = I3 (i, n, m) − I3 (i, n, j), dt m=1 j=1
0 = FQ (Q, R) =
M K M K X X X X dQik = I3 (i, k, m) − I3 (i, k, j) + I3 (k, i, m) − I3 (k, i, j), dt m=1 m=1 j=1 j=1
where I3 for the case of ReLU is provided in the Appendix Eq.27. These equations are derived from the gradient flow dynamics. It is easy to see, and discuss in detail in Appendix A, that the equations are described by the order paramters Eqs.6 introduced in the Main Text. Since the gradient depends on the weights only through their overlaps Q and R, any stationary point of the high-dimensional population risk Eq.28 is necessarily a stationary point of the summary statistics equations FR (Q, R) = 0, FQ (Q, R) = 0 for Gaussian inputs, therefore the flow must be zero giving Result 1. As discussed in Sec.3.2, our results are more general than just gradient flow on two-layer ReLU networks. In Appendix A, H, and G we discuss how the flow is derive more in general, including other activation function and the stochastic gradient dynamics case. 17
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Gaussian Mixture Model Fit for Q diag
Gaussian Mixture Model Fit for T diag
Empirical (Q diag) GMM total (n=2) Component 1 Component 2
5
Density
4
8
Empirical (T diag)
7
GMM total (n=1) Component 1
6 5
3
4 3
2
2 1 1 0
0.6
0.7
0.8
0.9
1.0
1.1
1.2
0
1.3
Element value
0.8
0.9
1.0
1.1
1.2
Element value
Figure 7: Gaussian mixture analysis of diagonal order parameters.The left panel shows the distribution of diagonal elements of the matrix Q, which is well described by a two-component Gaussian mixture model. The right panel shows the distribution of diagonal elements of the matrix T , which follows a single Gaussian distribution centred at 1. The data are obtained by grouping order parameters within the same minima family, re-arranging the corresponding (Q, R) matrices to maximise the maximise the absolute values of diagonal of Q and extracting the diagonal entries.
B.1
Ansatz Formulation
To characterise the structure of the order parameters, we employ a Gaussian Mixture Model (GMM) analysis, which represents a distribution as a weighted sum of Gaussian components and provides a simple tool for identifying clustered structure in empirical data. We focus on the family of minima empirically observed by the results from large-scale cluster runs (e.g., the same structure shown in Fig.1b). In order to show the order parameter matrices R and Q in a consistent way, we have to consider the permutational invariance of the units in the solution. Therefore, we re-arrange the student’s units to maximise the absolute values of the diagonal elements of the student-teacher overlap R. After alignment, we analyse the distribution of the diagonal elements of the matrices Q and T across the ensemble of minima by fitting these distributions using a GMM to identify distinct populations of hidden units. In Fig.7, we show an example of this procedure for elements in Q. In Fig. 7 (Left), the diagonal elements of Q exhibit a clear bimodal distribution. A two-component GMM reveals two well-separated Gaussian components, whose relative weights closely match the ratio k1 /k2 (with k2 denoting the size of the aligned group). This indicates that student units separate into two distinct populations during training, supporting the block-structured ansatz used in our theoretical analysis. As a control, in Fig.7 (Right) we apply the same GMM analysis to the diagonal elements of the teacher-teacher overlap matrix T . In contrast to Q, the distribution of T is well described by a single Gaussian component (n = 1), as shown in Fig. 7 (Right). This intuition is used to derive the ansatz reported in Sec.3.1. The block-symmetric ansatz proposed in Eq. 8 can be formally connected to the discrete group-theoretic analysis of Arjevani and Field [2019]. They demonstrate that critical points in two-layer ReLU networks under Gaussian inputs retain large isotropy subgroups (e.g., Sk1 × Sk2 ), leading to weight matrices with highly symmetric transitivity partitions. In our macroscopic framework, these discrete weight-space partitions manifest exactly as the uniform block-diagonal and off-diagonal constants in the Q and R order parameters. Therefore, the dynamical compensation mechanism we observe is the statistical-mechanics dual to their “principle of least symmetry breaking.” Examples of the Ansatz. We report here some examples for the case k1 = 2, k2 = 3 in the well-specified regime (K = M ): ς s φ φ φ q e u u u s ς φ φ φ e q u u u R = f f b τ τ , Q = u u p µ µ . (33) f f τ b τ u u µ p µ u u µ µ p f f τ τ b 18
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
and over-parametrised case (i.e. K = M + 1): ι ι κ κ κ ς s φ φ φ s ς φ φ φ R= , f f b τ τ f f τ b τ f f τ τ b
Ω v v q v e Q= γ u γ u γ u
v e q u u u
γ γ γ u u u u u u . p µ µ µ p µ µ µ p
(34)
To provide a direct visual check of the structural assumptions underlying the ansatz, we report a representative example of the order parameters (R∗ , Q∗ ) obtained from the simulations for K = 18 and M = 17, shown in Fig. 8. Loss=0.0001
Loss=0.0117
Loss=0.0195
Loss=0.0263
Loss=0.0319
Loss=0.0380
1.00 0.75
R
0.50
0.00
Value
0.25
0.25 0.50
Q
0.75 1.00
Figure 8: Different order parameters at K = 18 and M = 17. The plots show the values of R (first row) and Q (second row) for minima representative of the different families. From left to right, we see results for k1 = 0, 2, 3, 4, 5, 6.
B.2
Solution of the Fixed Point Equation Under the Ansatz
The fixed point equations of dQik /dt, dRin /dt (Eqs. 22 and Eqs. 21) simplify after substituting the ansatz defined in Eqs. 8. Different entries of Qik and Rin (e.g., e and ς) satisfy distinct equations, denoted by Ee , Eς , and so on. In the K = M setting, this procedure results in 11 coupled equations characterising fixed points reported below: ! p p p 1 e u [k1 ς q 2 − e2 − ς q 2 − e2 + (k1 − 1) q s cos−1 − p + f k2 q cos−1 − √ − k1 ς q − s2 2πq pq q2 ! i p p p p q ς + k2 ς pq − u2 − k2 ς q − ϕ2 + q ς cos−1 − p + ς q − s2 − ς q − ς 2 − q cos−1 − √ = 0, q q2
Eς = −
(35)
! p p p 1 e u [k1 s q 2 − e2 − s q 2 − e2 + q cos−1 − p (k1 − 2) s + ς + f k2 q cos−1 − √ − k1 s q − s2 2πq pq q2 ! p p p p q s i + k2 s pq − u2 − k2 s q − φ2 + q s cos−1 − p + s q − s2 − s q − ς 2 − q cos−1 − √ = 0, q q2 (36)
Es = −
!! !! p 1 u e q −1 −1 −1 2 2 Eφ = [−q b + (k2 − 1)τ cos −√ + φ −(k1 − 1) q − e + q cos −p − q cos −p 2πq pq q2 q2 i p p p p φ + (k1 − 1) φ q − s2 − k2 φ pq − u2 + k2 φ q − φ2 + φ q − ς 2 + q cos−1 − √ = 0, q (37)
19
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks ! p p p s 1 h e −1 −1 2 2 2 2 2 −√ Ee = q − s − q − e + qs cos −p k1 e + e q − e − q e(k1 − 2) + q cos πq q q2 ! p p p p u q (38) + k2 e q − φ2 − pq − u2 − qu cos−1 − √ − e q − s2 + e q − ς 2 − eq cos−1 − p pq q2 φ s ς s i + k2 φ q cos−1 − √ + q ς cos−1 − √ + qs cos−1 − √ − 2qs cos−1 − √ = 0, q q q q !! p e −1 2 2 ! p (k − 1) q − e + q cos − p 1 q2 b φ q 1 h u p − b2 −1 p + φ cos−1 − √ + b cos−1 − √ − u + cos − Eu = 2 2π p p q q q p ! p u p − f2 1 u f −1 −1 2 − p e(k1 − 1) + q cos −√ + s cos −√ + k1 u pq − u ) + (k1 − 1) p pq p p ! p p u p − f2 u q − s2 s f ς + + (k1 − 1) f cos−1 − √ + ς cos−1 − √ + + f cos−1 − √ p q q p q p p ! p2 − µ2 + p cos−1 − √µ 2 q (k2 − 1)µ + p cos−1 − √upq + k2 u pq − u2 (k2 − 1) p p − − u + cos−1 − p p q p2 ! p u p − τ2 τ + (k2 − 1) φ cos−1 − √ + + (k2 − 1) p p
p p ! q − φ2 u q − φ2 u q − ς2 i φ + τ cos−1 − √ + + = 0, q q q q (39) ! p p p p p 1 µ Ef = − − f k2 p − τ 2 − f p − b2 − f k1 p − f 2 + f k1 pq − u2 + f k2 p2 − µ2 + f (k2 − 1) p cos−1 − p 2πp p2 ! p p p f u u − f p2 − µ2 + f p cos−1 − p + f p − τ 2 − p cos−1 − √ + k1 ps cos−1 − √ − ps cos−1 − √ p pq pq p2 u = 0, + pς cos−1 − √ pq (40) ! p p p p p p 1 p Eb = − [−b p − b2 − bk1 p − f 2 + bk1 pq − u2 + bk2 p2 − µ2 − bk2 p − τ 2 − b p2 − µ2 + bp cos−1 − p 2πp p2 ! i p b u µ + b p − τ 2 − p cos−1 − √ + k1 pφ cos−1 − √ + (k2 − 1)pτ cos−1 − p = 0, 2 p pq p (41) ! p p p µ u 1 − k1 τ p − f 2 + k1 pφ cos−1 − √ Eτ = − [−τ p − b2 + p b + (k2 − 2)τ cos−1 − p + k1 τ pq − u2 2 2πp pq p ! p p p p p τ + k2 τ p2 − µ2 − k2 τ p − τ 2 − τ p2 − µ2 + pτ cos−1 − p + τ p − τ 2 − p cos−1 − √ ] = 0, p p2 (42) ! p p 1 e s Eq = [−(k1 − 1) q 2 − e2 − e(k1 − 1) cos−1 − p + (k1 − 1) q − s2 + s cos−1 − √ π q q2 ! p p p u φ q − k2 pq − u2 + u cos−1 − √ + k2 q − φ2 + φ cos−1 − √ − q cos−1 − p + q − ς 2 (43) pq q q2 i ς + ς cos−1 − √ = 0, q
20
u
p
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks p p 1 p b f u Ep = [ p − b2 + b cos−1 − √ p − f 2 + f cos−1 − √ pq − u2 + u cos−1 − √ + k1 − k1 π p p pq ! ! p p µ τ p −1 −1 −1 −p −√ − (k2 − 1) p2 − µ2 + (k2 − 1)(−µ) cos + k2 p − τ 2 + (k2 − 1)τ cos − p cos −p p p2 p2 p − p − τ 2 ] = 0, (44) p h p p b 1 τ f Eµ = µ p − b2 + pτ cos−1 − √ + bp cos−1 − √ + k1 µ p − f 2 − pq − u2 + f p cos−1 − √ πp p p p !! p p µ τ u p2 − µ2 + p cos−1 − p −pu cos−1 − √ − k2 µ + k2 µ p − τ 2 + k2 pτ cos−1 − √ pq p p2 (45) ! ! ! p µ µ p 2 −1 −1 −1 + µ p2 − µ2 − p cos −p + 2µp cos −p − µp cos −p p2 p2 p2 i p τ − µ p − τ 2 − 2pτ cos−1 − √ = 0. p
To validate the theoretical ansatz proposed in the Main Text, we solve the reduced fixed-point equations derived above. Specifically, we substitute the block-diagonal parameterization into the general mean-field ODEs (derived in Appendix A.3), transforming the high-dimensional optimization problem into a set of coupled non-linear algebraic equations. We used the solver ”FindRoot” implemented in Wolfram Mathematica to solve the fixed point equations. The function implements a Newton’s method implementation. The resulting solutions (R∗ , Q∗ ) constitute the theoretical fixed points of the mean-field dynamics. By evaluating the population loss function at these theoretical coordinates, we obtain the predicted energy barriers. These values correspond to the vertical dashed lines shown in Fig. 2 (Main Text), which demonstrate a close agreement with the peaks of the empirical loss histograms. Following the same procedure as in the case K = M , the above analysis can be straightforwardly generalized to the overparameterised setting K = M + 1. Due to their length, we do not report these equations explicitly. Instead, Fig. 9 displays the order parameters (R∗ , Q∗ ) obtained directly from the numerical solver for K = 18 and M = 17 in different k1 . The resulting structures and losses closely match those shown in Fig. 8. Loss=0.0006
Loss=0.0119
Loss=0.0194
Loss=0.0262
Loss=0.0326
Loss=0.0388
1.00 0.75
R
0.50
0.00
Value
0.25
0.25 0.50
Q
0.75 1.00
Figure 9: Order parameters obtained from solvers at K = 18 and M = 17. The plots show the values of R (first row) and Q (second row) for minima computed by the solver. Columns are ordered from left to right by the number of anti-aligned units, k1 = 0, 2, 3, 4, 5, 6.
C
Robustness to Finite Input Dimensionality
In the main text, our analysis of the fixed points we assume orthonormal teacher configuration (T = IM ) for simplifying the equations. While this is exact in the thermodynamic limit (d → ∞), in this section, we investigate the robustness of our predictions to finite input dimensions d. We simulate the dynamics starting from actual neural network weight initialisations at dimensions d ∈ {196, 392, 784}. 21
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Fig. 10 illustrates the impact of finite dimensionality on the final population loss distribution. We compare the empirical histograms obtained from finite d simulations against the idealised infinite-dimensional case. While the discrete, quantised hierarchy of the local minima is strictly preserved, finite dimensions introduce variance. This results in a progressive broadening of the loss distribution peaks around the theoretically predicted values. As expected, the empirical distributions tightly converge toward the analytical infinite-dimensional limit as d increases. To further validate that the structural symmetries discussed in our ansatz (Sec.3.1) persist beyond the idealised orthonormal teacher scenario, Fig.11 visualises the student-teacher overlap matrix R. We uniformly sample 10 configurations at the end of the dynamics from the first-order local minimum (i.e., the second quantised group corresponding to k1 = 1, as observed in Fig.10). Despite the noise introduced by the finite-dimensional teacher, the distinct macroscopic block structure—including the aligned and anti-aligned unit clusters—remains clearly identifiable. As d decreases, the off-diagonal fluctuations become more pronounced, yet the fundamental geometric organisation in blocks remains robust. This confirms that our ansatz accurately reflects the empirical behaviour of finite-dimensional networks.
300
orthonormal (infinity) d=784 d=392 d=196
250
Density
200 150 100 50 0 0.00
0.02
0.04 Population loss
0.06
0.08
orthonormal (infinity)
Figure 10: Impact of finite dimensionality on the loss distribution. Histograms of the final population loss obtained by simulating the ODEs starting from actual neural network weight initialisations at dimensions d ∈ {196, 392, 784}, compared against the infinite-dimensional limit (orthonormal teacher, T = I). Dashed vertical lines indicate the corresponding theoretically predicted loss values. While smaller dimensions introduce a broadening effect due to variance, the discrete families of minima remain robustly centred around the predicted macroscopic levels.
1.0
d=392
0.0
Value
d=784
0.5
d=196
0.5
1.0
Figure 11: Student-teacher overlap matrix (R) across different dimensions. Heatmaps of the R matrix for 10 randomly sampled configurations converging to the first-order local minimum (k1 = 1). The top row represents the idealised d → ∞ limit, while subsequent rows correspond to finite dimensions d = 784, 392, and 196. The blocksymmetric structure predicted by our theoretical ansatz clearly emerges across all regimes, with finite-size fluctuations increasing at lower dimensions.
22
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
D
String method
(a) current string discrete images on a path
(c) reparameterize reparameterize {xi } uniform in arclength
(b) one update step update interior images only xi ← xi + ∆xi (e.g., − η∇f )
B
B
A
A
endpoints fixed update i = 1, . . . , N − 2
B
A
Figure 12: Illustration of the string method. In this section, we briefly introduce the string method [Weinan et al., 2002, Ren et al., 2007, Samanta and Weinan, 2013], which is used to investigate whether distinct minima are separated by energy barriers or connected by flat valleys. This algorithm seeks the Minimum Energy Path (MEP) connecting two fixed configurations in the order parameter space. Conceptually similar to the Nudged Elastic Band method [Draxler et al., 2018], this technique enables the study of low-loss trajectories between distinct points. We illustrate the string method in Fig. 12. Given two fixed endpoints, the string method approximates a continuous path in parameter space by a discrete set of images. These images evolve under gradient descent, and are periodically reparameterised to ensure approximately uniform spacing with respect to the chosen metric. D.1
The Induced Metric and Interpolation
A crucial requirement for the string method is a proper definition of distance to discretise the path. Since the gradient descent dynamics are defined on the weights W ∈ RK×d , the geometry of the landscape is intrinsically Euclidean in the weight space. Therefore, the string evolution in the reduced order parameter space must rely on an induced metric that preserves this geometry. Given the usage of the mean-field description based on the order parameters (sufficient statistics) Θ = (R, Q), we lose the information of the individual weights. We define the distance between two macroscopic states Θ1 and Θ2 as the expected squared Euclidean distance between their typical microscopic realizations W1 and W2 : 1 2 2 D (Θ1 , Θ2 ) = EW1 ,W2 ∥W1 − W2 ∥F , (46) d where the expectation is taken over the independent ensembles of weights consistent with the respective order parameters. Expanding the Frobenius norm yields: 1 ∥W1 − W2 ∥2F = Tr d
W1 W1⊤ W2 W2⊤ W1 W2⊤ W2 W1⊤ + − − d d d d
By definition of the order parameters, the self-overlap terms converge to their expectations: 1 ⊤ E Wi Wi = Qi . d
.
(47)
(48)
To evaluate the cross-term W1 W2⊤ , we decompose the weights into a signal component (parallel to the teacher subspace) and a noise component (orthogonal to the teacher): W = W∥ + W⊥ = RW ∗ + Ξ,
(49)
where W ∗ ∈ RM ×d are the teacher weights and Ξ represents the orthogonal noise. Substituting this decomposition into the cross-term: ! ∗ ∗⊤ 1 1 W W W1 W2⊤ = (R1 W ∗ + Ξ1 )(R2 W ∗ + Ξ2 )⊤ = R1 R2⊤ + noise cross-terms. (50) d d d 23
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
We assume that W1 and W2 are drawn as independent samples from their respective macroscopic ensembles. Under this assumption, the noise components Ξ1 and Ξ2 are statistically uncorrelated. In the limit d → ∞, the inner product of independent random vectors vanishes: d1 Ξ1 Ξ⊤ 2 → 0. Furthermore, assuming orthonormal teacher weights 1 ∗ ∗⊤ (T = d W W = I), the expectation simplifies to: 1 (51) E W1 W2⊤ ≈ R1 R2⊤ . d Substituting these expectations back into the distance equation, we obtain the metric employed in our analysis: q dist(Θ1 , Θ2 ) = Tr Q1 + Q2 − R1 R2⊤ − R2 R1⊤ .
(52)
We emphasise that the independence assumption used to derive Eq. (52) implies that the microscopic noise components Ξ are uncorrelated. While valid for random initialisations in high dimensions, this assumption does not strictly hold for consecutive points along a continuous trajectory, where the noise components are necessarily correlated. However, the Minimum Energy Path is a geometric object that is invariant under reparameterisation of the curve. The specific choice of metric serves only to distribute the discretisation points (images) along the string. Therefore, Eq. (52) constitutes a valid parameterisation choice that respects the global Euclidean scaling of the problem, even if it overestimates local distances in the small-step limit. D.2
Algorithm
The string is discretized into P + 1 images {Θ0 , Θ1 , . . . , ΘP }, where Θi = (Qi , Ri ). The endpoints Θ0 and ΘP are fixed at the known minima in our case. The algorithm proceeds iteratively: 1. Evolution: Each interior image Θi (i = 1, . . . , P − 1) evolves independently for a fictitious time step ∆τ according to the mean-field dynamics derived in Appendix A.3: Θi ← Θi + ∆τ · F(Θi ),
(53)
where F(Θ) = (Q̇, Ṙ) denotes the vector field of the order parameters (defined by Eqs. 21 and Eqs. 22). This step drives the string towards the valley floor (local stationarity). 2. Reparameterization: To prevent the images from collapsing into the minima, we redistribute them to maintain uniform spacing. We first compute the total arc-length L of the current string using the metric in Eq. (52). Then, we construct a cubic spline interpolation of the path and resample P + 1 new images at equidistant arc-length coordinates sk = Pk L (for k = 0, . . . , P ). 3. Termination: The algorithm is iterated for a sufficient number of steps to ensure the string relaxes to a stationary state. Upon convergence, the transverse forces vanish, and the resulting curve represents the Minimum Energy Path.
E
Stability of fixed points from higher order minima
As noted in the Main Text, the non-differentiability of the ReLU activation function at the origin precludes a standard Hessian-based stability analysis. We employ a perturbative approach to probe the local landscape geometry. Given a fixed point configuration W̄ derived by solving the fixed point equations in the well-specified setting (Appendix B.2), we apply an isotropic Gaussian perturbation in the weight space to see if it acts as a stable attractor or a saddle point: Winit = W̄ + ξ, with ξ ∼ N (0, σ 2 I) (54) where σ represents the strength of the perturbation. In practice, the perturbation is done in the mean-field description using Q and R. Notice that directly perturbing Q and p √ R may result in non-realisable configurations, in particular Q ⪰ 0, Qij = Qji ≤ Qii Qjj , and Rim ≤ Qii Tmm . The system is then allowed to relax for 100, 000 steps following the gradient flow dynamics with η = 0.01. We monitor the discrepancy in loss value between the unperturbed configuration and the perturbed one. The perturbation is repeated 20 times with different seeds, and the average and the variance of the loss gap are recorded and shown in Fig.3B. of the Main Text. 24
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Our empirical observation that spurious fixed points become unstable for K ≥ M + 1 is supported by the exact local Hessian analysis of Safran et al. [2021]. They mathematically proved that taking a non-global minimum in the K = M landscape and “splitting” a neuron to embed it in the K = M + 1 landscape necessarily introduces a strictly negative eigenvalue in the Hessian block corresponding to the duplicated units. When we introduce the Gaussian perturbation ξ to the extended weight matrix, the system is pushed off the exact critical point. Because the Hessian now possesses a direction of negative curvature, the gradient flow does not return to the spurious state, but instead escapes the newly formed saddle point and descends toward the continuous global minimum manifold.
F
Additional String Method Results in the Overparameterised Regime
To complement the representative examples shown in the main text, we report here additional string-method results connecting equal-loss minima within the k1 = 2 family. These configurations are constructed by exchanging the roles of extra, aligned, and anti-aligned neurons. More precisely, Fig. 13 shows four different endpoint constructions, where in each case the left endpoint corresponds to the original configuration of (Q, R), while the right endpoint is obtained by permuting the student neurons, as indicated by the legend in Fig. 14. The corresponding loss values along the strings are shown in Fig. 14. These additional experiments reveal a richer phenomenology than the representative flat strings discussed in the main text. In particular, when k1 = 2, the transition between neurons of the same kind (red) remains essentially flat, while the transitions between an extra neuron and an anti-aligned neuron (blue), and between an aligned neuron and an anti-aligned neuron (green), collapse towards zero loss along the path, indicating directions that connect these configurations to the global minimum. By contrast, the path connecting an extra neuron and an aligned neuron (orange) retains a clearly visible barrier. This observation does not contradict the existence of flat manifolds in the overparameterised regime. Rather, it shows that connectivity among minima in the landscape is not fully connected, but—in the case of k1 = 2—it is still clustered: while some pairs of equal-loss minima are connected by flat directions, other interpolation paths may still involve excursions to higher loss. Interestingly, for k1 > 1, all the minima appear to be fully connected. As a final note, we stress that the existence of flat directions is already sufficient to show that the isolated minima of the well-specified regime (K = M ) merge into continuous manifolds when K = M + 1. extra neuron & anti-aligned neuron
0
extra neuron & aligned neuron
0
aligned neuron & anti-aligned neuron
0
same kind of neurons
0
Origin
0
1.00 0.75
2
2
2
2
2
4
4
4
4
4
0.50
6
6
6
6
6
0.25
8
8
8
8
8
0.00
10
10
10
10
10
12
12
12
12
12
14
14
14
14
14
16
16
16
16
16
0.0
2.5
5.0
7.5
10.0 12.5 15.0
0.0
2.5
5.0
7.5
10.0 12.5 15.0
0.0
2.5
5.0
7.5
10.0 12.5 15.0
0.0
2.5
5.0
7.5
10.0 12.5 15.0
0.25 0.50 0.75 0.0
2.5
5.0
7.5
10.0 12.5 15.0
1.00
Figure 13: R of different settings. Examples of endpoint configurations within the k1 = 2 family in the overparameterised regime, shown here through the corresponding R matrices. For the strings in Fig. 14, the left endpoint corresponds to the origin R, while the right endpoint is obtained by permuting the student neurons.
G
Alternative Activation Functions
In this section, we provide the analytical expressions for the Gaussian integrals I2 and I3 for Leaky ReLU and Error Function (erf) activations. Based on these derivations, we also present numerical observations of the learning landscape and spurious minima structures for these activations in this section. The notation follows the definitions in Appendix A: I2 (x1 , x2 ) = ⟨g(x1 )g(x2 )⟩, ′
I3 (x0 , x1 , x2 ) = ⟨g (x0 )x1 g(x2 )⟩, where the expectations are taken over a multivariate Gaussian distribution with covariance matrix Σ. 25
(55) (56)
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
0.025
extra neuron & anti-aligned neuron extra neuron & aligned neuron aligned neuron & anti-aligned neuron Same kind of neurons
0.020
Loss
0.015 0.010 0.005 0.000 0.0
0.2
0.4
Coordinate
0.6
0.8
1.0
Figure 14: Strings of different settings. Loss values along strings connecting the endpoint configurations shown in Fig. 13. The red curve remains essentially flat, whereas the blue and green curves collapse to zero loss, revealing directions that connect these configurations to the global minimum. By contrast, the orange curve displays a visible barrier. This illustrates that in the overparameterised regime the surviving k1 = 2 minima belong to a connected but direction-dependent landscape geometry.
G.1
Leaky ReLU
The Leaky ReLU activation function with a leakage parameter α ∈ [0, 1] is defined as: x if x ≥ 0, g(x) = LReLU(x; α) = αx if x < 0.
(57)
Its derivative is the generalized step function g ′ (x) = I(x ≥ 0) + αI(x < 0). Note that setting α = 0 recovers the standard ReLU results, while α = 1 corresponds to a linear network. Two-Variable Case
The explicit form for the correlation between two Leaky ReLU units is: q C12 (1 − α)2 LReLU 2 C11 C22 − C12 + C12 arccos − √ I2 (Σ) = αC12 + , 2π C11 C22
(58)
where C11 , C22 are the variances and C12 is the covariance of (x1 , x2 ). Three-Variable Case
The integral involving the derivative, a multiplier, and the activation is given by: q (1 − α)2 C02 C01 LReLU 2 I3 (Σ) = αC12 + C12 arccos − √ + C22 C00 − C02 . 2π C00 C22 C00
(59)
Here, the covariance matrix for the joint variables (x0 , x1 , x2 ) is: Σ=
G.2
C00 C01 C02
C01 C11 C12
C02 C12 C22
! .
Leaky ReLU Results
Based on the derived integrals, we numerically investigate the stationary points of networks with leaky ReLU activation. With the Leaky ReLU parameter set to α = 0.01, the Fig. 15 shows the results for K = M = 12 under nGD. It can be clearly observed that anti-aligned units appear, similar to the ReLU case. 26
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Loss=0.0014
Loss=0.0014
Loss=0.0179
Loss=0.0184
Loss=0.0334
Loss=0.0340
1.00 0.75
R
0.50
0.00
Value
0.25
0.25 0.50
Q
0.75 1.00
Figure 15: Order Parameters of Leaky ReLU in K = M = 12
G.3
Sigmoidal (erf)
Consider the sigmoidal activation function defined by the error function: x g(x) = erf √ . 2 q 2 The derivative is proportional to a Gaussian: g ′ (x) = π2 e−x /2 . Two-Variable Case (I2 )
(60)
Based on the arcsine law for Gaussian integrals, the result is: ! 2 C12 erf I2 (Σ) = arcsin p . π (1 + C11 )(1 + C22 )
(61)
Three-Variable Case (I3 ) The integral involving the derivative of the error function allows for an analytical solution: 2 C01 C02 I3erf (Σ) = p C − . (62) 12 2 1 + C00 π (1 + C00 )(1 + C22 ) − C02 G.4
Results for Sigmoidal (erf) Activation
Based on the derived integrals, we numerically investigate the stationary points of networks with erf activation by running the ODEs for a few seeds. Fig. 16 illustrates the structure of the order parameters for K = M = 12 across different local minima found by nGD. Notice the absence of anti-aligned neurons in sharp contrast with what is observed in the Main Text for ReLU and in Fig.15 for Leaky ReLU. Loss=0.0014
Loss=0.0090
Loss=0.0148
Loss=0.0161
Loss=0.0380
Loss=0.0309
1.00 0.75
R
0.50
0.00
Value
0.25
0.25 0.50
Q
0.75 1.00
Figure 16: Order Parameters of erf in K = M = 12
27
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
H
Constraint Dynamics
In this section, we derive the mean-field ODEs for Normalised Gradient Descent (nGD) and Orthonormalized Gradient Descent (onGD), which impose distinct geometric constraints on the student weights during training. We also show the resulting learning dynamics for both methods and a brief discussion of stability, convergence properties, and the learnability of nGD and onGD. Consistent with previous sections, we denote the input dimension by d. H.1
Normalised Student Setting (nGD)
In Normalised Gradient Descent, we enforce the constraint that the norm of each student weight vector remains fixed at initialisation, typically ∥wi ∥2 = d (implying Qii = 1). The continuous-time update rule subtracts the radial component of the gradient: dwi η = −η∇wi L + (∇wi L · wi⊤ )wi . (63) dt d This projection modifies the order parameter dynamics by introducing a decay term proportional to the self-overlap. The resulting ODEs are: K K M M X X X X dRin ∗ ∗ = ηvi vj I3 (i, n, j) − ηRin vi vj I3 (i, i, j) , (64) vm I3 (i, n, m) − vm I3 (i, i, m) − dt m=1 m=1 j=1 j=1 M K M K X X X X dQik ∗ ∗ = ηvi vm I3 (i, k, m) − vj I3 (i, k, j) + ηvk vm I3 (k, i, m) − vj I3 (k, i, j) dt m=1 m=1 j=1 j=1 M K M K X X X X ∗ ∗ − ηQik vi vm I3 (i, i, m) − vj I3 (i, i, j) − ηQik vk vm I3 (k, k, m) − vj I3 (k, k, j) . m=1
m=1
j=1
j=1
(65) ii Note that for the diagonal elements, substituting k = i and Qii = 1 confirms that dQ dt = 0, satisfying the constraint.
H.2
Empirical Observations for Normalized Gradient Descent
We investigate the loss distribution of Normalized Gradient Descent (nGD) in the overparameterised regime. Specifically, we consider a student network with K = 19 hidden units learning from a teacher with M = 17 units (K = M + 2). Fig. 17A. displays the histogram of the final population loss. Contrary to the expectation that overparameterisation allows the network to achieve zero loss, nGD exhibits a fundamental barrier. As shown in the histogram, the system fails to reach the rigorous global minima (Loss ≈ 0). Even in the best-found minima (Loss ≈ 0.355), a significant error persists. This is a direct consequence of the constraint Qii = 1. In unconstrained GD, the K − M excess students would decay to zero norm to eliminate interference, as shown in Fig. 8. In nGD, however, the order parameters shown in Fig. 18 (left columns) indicate that these excess units are constrained to maintain unit norm. They cannot be silenced and effectively act as intrinsic noise sources, creating an irreducible error even when the M target features are perfectly retrieved. The high density of solutions in the suboptimal interval Loss ∈ [0.365, 0.370] (approx. 48%) arises from the combinatorial diversity of ”blocking” defects. As visualised in the right four columns of Fig. 18, these states exhibit student groupings of varying sizes (e.g., k1 = 3, 4, 5). Unlike the unique diagonal alignment required for the best solution, there are combinatorially many ways to form these suboptimal clumps. The normalisation constraint stabilises these high-loss configurations, effectively trapping the dynamics in a high-energy part. H.3
Orthonormalized Student Setting (onGD)
In Orthonormalized Gradient Descent, we enforce the stronger constraint that the student weights remain orthonormal throughout training, i.e., W W ⊤ = dIK (implying Q = IK ). The update rule projects the gradient onto the tangent space of the Stiefel manifold: dW η = −η∇W L + (∇W L W ⊤ )W, (66) dt d 28
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
K=19, M=17, nGD
300
K=17, M=17, onGD 25 20
200
Density
Density
250 150 100
15 10 5
50 0 0.355
0.360
0.365
0.370
0.375
Population loss
0.380
0.385
0 0.18
0.390
A. nGD Loss Distribution
0.20
0.22
0.24
0.26
Population loss
0.28
0.30
B. onGD Loss Distribution
Figure 17: Loss distribution in different constraints These histograms show the empirical density of the final population loss after 20,000,000 running steps. (a) Normalized Gradient Descent (nGD) with K = 19, M = 17. Despite overparameterisation, the loss distribution is strictly multimodal, exhibiting discrete peaks at high loss values. (b) Orthonormalised Gradient Descent (onGD) with K = 17, M = 17. The distribution is unimodal and broad, concentrated entirely in a high-error regime (centered at 0.22). Loss=0.3557
Loss=0.3559
Loss=0.3676
Loss=0.3695
Loss=0.3715
Loss=0.3714
1.00 0.75
R
0.50
0.00
Value
0.25
0.25 0.50
Q
0.75 1.00
Figure 18: Order parameters of local minima in nGD (K = 19, M = 17). The columns correspond to different final losses. Left Two Columns (Best Minima): The R matrices exhibit a near-perfect diagonal structure, indicating successful retrieval of the 17 teacher features. However, the irreducible error persists because the 2 excess students (constrained to Qii = 1) cannot decay to zero. Right Four Columns (Suboptimal Minima): These states represent the dominant local minima cluster. They exhibit ”blocking” defects (e.g., 2-to-1 assignments seen as red blocks in R), where extra units are misallocated, further elevating the loss. where W ∈ RK×d . Under this setting, the student-student covariance matrix is fixed, so: dQik = 0. (67) dt The dynamics of the student-teacher overlap R are modified by a mixing term involving all other student units: M K K M K X X X X X dRin ∗ ∗ = ηvi vm I3 (i, n, m) − vj I3 (i, n, j) − ηvi Rkn vm I3 (i, k, m) − vj I3 (i, k, j) . dt m=1 m=1 j=1 j=1 k=1
(68) H.4
Empirical Observations for Orthonormalized Gradient Descent
We next examine the learning dynamics of Orthonormalized Gradient Descent (onGD). We focus on the student-teacher setting with K = M = 17. Fig. 17B. presents the distribution of the final population loss. As shown in Fig. 17B., onGD does not achieve meaningful learning. The final loss distribution is concentrated in the high-error regime of [0.18, 0.30], with a peak around 0.22. The structural reason for this failure is evident in Fig. 19: 29
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
Loss=0.2061
Loss=0.2143
Loss=0.2140
Loss=0.2167
Loss=0.2176
Loss=0.2053
1.00 0.75
R
0.50
0.00
Value
0.25
0.25 0.50
Q
0.75 1.00
Figure 19: Final order parameters for onGD (K = M = 17). The columns correspond to different final losses. Bottom Row: The student-student overlap matrices Q retain a strict identity structure (Q = I), satisfying the orthogonality constraint. Top Row: The student-teacher overlap matrices R exhibit a disordered, noise-like pattern. Unlike successful learning scenarios shown in Fig. 8, the entries remain diffuse and small. This lack of structural alignment confirms that the student units fail to retrieve the teacher vectors in this constrained setting. the entries of the student-teacher overlap matrix R remain disordered, confirming that the network is unable to retrieve the teacher configuration. The failure of onGD appears to be closely related to the excessive rigidity imposed by the orthogonality constraint. In a typical successful learning trajectory (as seen in unconstrained GD), student units often pass through intermediate phases where they become correlated (Qik ̸= 0) to ”sense” the teacher’s structure. By enforcing Qik = 0 at i ̸= k, onGD effectively forbids these cooperative intermediate states. Our results suggest that for K = 17, the gradient flow on this constrained manifold lacks accessible paths connecting the random initialization to the teacher configuration, effectively locking the system in a high-loss state. To comparison, we present representative learning trajectories in a setting with few hidden units (K = M = 2). As shown in Fig. 20, successful learning is possible in this regime, indicating that the obstruction induced by the orthogonality constraint is dimension-dependent: systems with only a few hidden units can still navigate the Stiefel manifold to find the solution. (b) Final R
(a) Training loss
Loss
Value
10 1
0.6 0.4 0.2 0.0 0.2 0.4 0.6
10 2 0.0
0.2
0.4
Training step
0.6
0.8
1.0 1e6
Figure 20: Dynamics of onGD in a setting with small hidden layers (K = M = 2). (a) The training loss (log scale) decreases monotonically and converges to a near-zero value (∼ 10−3 ), indicating successful optimisation. (b) The final student-teacher overlap matrix R exhibits a clear diagonal structure (red blocks indicate high positive correlation). This confirms that, unlike the case with larger hidden layers (K = M = 17), the orthogonality constraint perfectly retrieves the teacher weights in low dimensions.
I
Spurious minima in two-layer ReLU networks
We reproduce the gradient descent results for ReLU soft committee machines using mean-field dynamics over 10,000 random initialisations. Table 3a reports the proportions of global minima, spurious local minima, and non-minima in the balanced case K = M , while Table 3b presents the corresponding statistics in the overparameterised setting K = M + 1. The observed proportions are in close agreement with those reported in Safran and Shamir [2018].
30
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
a Results for case K = M .
b Results for the overparameterised case K = M + 1.
K
M
Local (%)
Global (%)
None (%)
K
M
Local (%)
Global (%)
None (%)
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.91 5.30 13.16 21.15 33.27 43.74 54.62 63.87 71.69 77.99 82.48 86.79 90.42 92.54 94.76
98.82 94.37 86.09 78.53 66.38 55.92 45.12 35.93 28.08 21.81 17.44 13.18 9.51 7.38 5.21
0.27 0.33 0.75 0.32 0.35 0.34 0.26 0.20 0.23 0.20 0.08 0.03 0.07 0.08 0.03
9 10 11 12 13 14 15 16 17 18 19 20
8 9 10 11 12 13 14 15 16 17 18 19
0.07 0.75 1.93 3.83 5.74 8.37 13.50 18.12 28.58 36.88 42.25 51.46
99.74 98.95 97.26 94.25 91.42 87.46 81.16 75.36 65.47 57.50 52.62 44.24
0.19 0.30 0.81 1.92 2.84 4.17 5.34 6.52 5.95 5.62 5.13 4.30
Table 3: Statistics of minima under mean-field dynamics. Proportions of local, global, and non-minima found by integrating the mean-field ODEs across 10,000 random seeds for each (K, M ) pair.
31
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
I.1
Extended Empirical Results for Constrained Dynamics
In this section, we provide a more granular look at the empirical convergence properties of the constrained dynamics, specifically Normalized Gradient Descent (nGD) and Orthonormalized Gradient Descent (onGD). Tables 4 and 5 report the statistical breakdown of the final states reached by nGD for the well-specified (K = M ) and mildly overparameterised (K = M + 1) regimes, respectively. Following the methodology of Safran and Shamir [2018], we define an empirical threshold ε to classify a run as having reached a global minimum, and a gradient norm threshold δ (evaluated on dQ/dt and dR/dt) to identify non-convergent trajectories. As reported in the rightmost columns of the tables, these thresholds were selected via visual inspection of the loss histograms and vary depending on the network width. Notably, we observe a general slowdown in the dynamics for the normalised case compared to unconstrained gradient descent. This slowdown is particularly pronounced in lower-dimensional configurations, where the system is more heavily restricted by the constraint. Because the spherical constraint shrinks the feasible space and alters the local geometry, navigating to the global minimum becomes progressively harder, requiring careful tuning of the threshold. Table 4: Statistics of minima in nGD under mean-field dynamics for case K=M. Proportions of local, global, and non-minima obtained by integrating the mean-field ODEs under Normalized GD over 10,000 random seeds for each (K, M ) with K = M . A result is classified as a global minimum when its error is smaller than ε. For the result whose errors are larger than ε, if is gradients (dQij /dt or dRin /dt) are larger than δ, it will be classified as None. Otherwise it will be classified as a local minimum. From the analysis in Section H, the global minima of nGD are very large. K
M
Local (%)
Global (%)
None (%)
ε
δ
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.00 0.15 0.02 0.01 0.00 0.62 4.87 12.93 24.15 36.45 49.95 58.57 67.92 75.37 81.49 85.88 89.77 92.01 94.38
99.77 99.85 99.98 99.99 100.00 99.38 95.13 87.07 75.85 63.55 50.05 41.43 32.08 24.63 18.51 14.12 10.23 7.99 5.62
0.23 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2
1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2 1e-2
Table 5: Statistics of minima in nGD under mean-field dynamics for case K=M+1. Proportions of local, global, and non-minima obtained by integrating the mean-field ODEs under Normalized GD over 10,000 random seeds for each (K, M ) with K = M + 1. A result is classified as a global minimum when its error is smaller than ε. For the result whose errors are larger than ε, if is gradients (dQij /dt or dRin /dt) are larger than δ, it will be classified as None. Otherwise it will be classified as a local minimum. From the analysis in Section H, the global minima of nGD are very large. K
M
Local (%)
Global (%)
None (%)
ε
δ
3 4 5 6
2 3 4 5
2.73 48.95 1.34 1.40
97.27 51.05 98.66 98.60
0.00 0.00 0.00 0.00
0.123 0.11 0.105 0.105
1e-1 1e-1 1e-1 1e-1
Continued on next page
32
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
K
M
Local (%)
Global (%)
None (%)
ε
δ
7 8 9 10 11 12 13 14 15 16 17 18 19 20
6 7 8 9 10 11 12 13 14 15 16 17 18 19
2.87 1.80 1.49 2.57 4.25 6.61 10.58 15.64 21.55 27.54 34.73 41.66 48.76 55.83
97.12 98.19 98.48 97.42 95.74 93.39 89.42 84.36 78.45 72.46 65.27 58.34 51.24 44.17
0.01 0.01 0.03 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1
1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1 1e-1
In contrast to nGD, we observe a complete failure of convergence when enforcing the stricter Stiefel manifold constraints of onGD. Because the network never successfully reaches the global minimum in our simulations, tabulating the convergence statistics gives no additional insights. Instead, we visualise the full empirical distribution of the final population loss in Figures 21 and 22.
Density
The top rows of both figures display the well-specified case (K = M ), while the bottom rows show the overparameterised case (K = M + 1). Notice that the final loss distributions for onGD do not exhibit the discrete, quantised band structure characteristic of both unconstrained GD and nGD. This lack of quantisation is a direct consequence of the orthogonality constraint. Indeed, the compensation mechanism—which relies on aligned neurons adjusting their partial overlaps and magnitudes to offset the error of anti-aligned units—is strictly forbidden on the Stiefel manifold. Without the ability to form these cooperative partial overlaps, the network cannot settle into the structured spurious families. K=2, M=2
800 700 600 500 400 300 200 100 0
K=3, M=3
50
50
50
40
40
40
30
30
30
20
20
20
10
10
10
0
0
Density
K=3, M=2
K=4, M=3
500
60
400
50 40 30
200
20
100 0 0.00
10 0.05
0.10 Loss
0.15
0.20
K=5, M=5 60
60
60
300
K=4, M=4
70
70
0 0.00
0.05
0.10 Loss
0.15
80 70 60 50 40 30 20 10 0.20 0 0.00
0
K=5, M=4
0.05
0.10 Loss
K=6, M=5
0.15
70 60 50 40 30 20 10 0.20 0 0.00
0.05
0.10 Loss
0.15
0.20
Figure 21: Distribution of loss for small networks in onGD under mean-field dynamic. Histograms of the final population loss obtained after long-time integration of the mean-field ODEs under Orthonormalised GD over 10,000 random initialisations. Each panel corresponds to a small student-teacher size configuration (K, M ). The top row shows the equal case K = M and the bottom row shows the overparameterised case K = M + 1. All panels share the same horizontal axis for direct comparison.
33
Sharp description of local minima in the loss landscape of high-dimensional two-layer ReLU neural networks
K=7, M=7
70 60 Density
K=15, M=15
K=19, M=19
30
20
40
25
40
30
20
15
30
20
15
10
50
20
10
10
10 0
0
0
K=12, M=11
60
40
50
30
K=20, M=19
25 20 15
20
10
20
10
10 0
0
K=16, M=15
30
40 30
5
5
K=8, M=7
70
Density
K=11, M=11
50
0.1
0.2
Loss
0.3
0.4
0.5
0
5 0.1
0.2
Loss
0.3
0.4
0.5
0
0.1
0.2
Loss
0.3
0.4
20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.5 0.0
0.1
0.2
Loss
0.3
0.4
0.5
Figure 22: Distribution of loss for wider networks in onGD under mean-field dynamic. Histograms of the final population loss obtained after long-time integration of the mean-field ODEs under Orthonormalised GD over 10,000 random initialisations. Each panel corresponds to a student-teacher size configuration (K, M ) with larger widths than those shown in Fig. 21. The top row shows the equal case K = M and the bottom row shows the overparameterised case K = M + 1. All panels share the same horizontal axis for direct comparison.
34