Stable Minima of ReLU Neural Networks Suffer from the Curse of Dimensionality: The Neural Shattering Phenomenon

Thanks for your support and detailed feedback! We will check and revise the minors. Responses to your questions below:

For example, this means that the lower bound of theorem 3.5 could be achieved for exotic configurations inside $\mathcal{F}_{\text{stable}}(\eta)$ that are not the solutions found by a typical optimization algorithm (as for theorem 3.7 see "Clarity" above). This fact limits the insights provided by the paper.

Thank you for mentioning this critical point. We argue that the configurations used in our lower bound are not exotic because they are precisely the structures that emerge from a typical training process, which our experiments empirically confirm.

Our work uses the set of stable minima as a prism to understand the dynamics of GD's representation leanring. Within this paradigm, our lower bound construction (see Append B.4 Figure 9 for the visualization of our lower bound construction) provides a intuition: it proves the existence of a "hard-to-learn" solution built from pseudo-sparse neurons, hinting at a dynamic where neurons pathologically shatter the dataset and get 'stuck' fitting local points. Moreover, the higher the dimension is, the easier for these neuron to shatter the dataset.

Furthermore, our experiments then provide direct empirical evidence for this theoretical construction. As shown in (Figure 1 right panel, Appendix A.1, A.2), vanilla Gradient Descent converges to a solution containing many neurons that exhibit these exact characteristics: a very low activation rate but a high weight norm. This demonstrates that the fundamental components of our lower bound represents a tangible failure mode produced by standard optimization algorithms in high dimensions.

Intuition on the weight function I interpret the weight function as a multiplication between (a) the expected activations of a neuron and (b) the probability of a neuron being activated times a data norm in the case the neuron is activated. In particular, what is the role of the (b) term ?

This is a great question! As derived in [Nacson et al. 2023, Appendix F], the form of $\tilde{g}( {u},t)$ is derived by brute-force computation of the loss function's Hessian matrix for a two-layer ReLU network:$$ \nabla_{\theta}^2\mathcal{L}=\frac{1}{n}\sum_{i=1}^n\nabla_{ {\theta}}f_{\theta}( x_i)\nabla_{\theta}f_{\theta}(x_i)^{\top} + \sum_{i=1}^n(f_{\theta}(x_i)-y_i)\nabla^2_{\theta}f_{ \theta}(x_i)

Thank you for your support and insightful comments! Responses to your questions below:

Proposition 2.2: I don't see how this can be true for any $\epsilon$, if this where true then wouldn't this mean that you get convergence to the minimizer from any initial condition? Perhaps another way of putting it is that the linearized dynamics can't be a good approximation for any $\theta$ right unless unless the loss is a quadratic form in the parameters?

Our proposition is not a claim about global convergence from an arbitrary starting point. Instead, it is about a local-stability condition. The purpose is to use this local analysis to derive a static condition that characterizes the stable minima. The "initial condition" in our stability definition refers to the moment an iterate enters the local quadratic region of a local minimum, $\theta^\star$. The proposition then shows that the curvature condition, $\lambda_{max}(\nabla^2\mathcal{L}(\theta^\star))\leq 2/\eta$, is necessary and sufficient for the linearized dynamics to remain trapped near $\theta^\star$. A rigorous proof can be found in [Qiao et al 24, Appendix C].

It is crucial to note that while this mathematical proposition about the linearized system holds for any $\epsilon$, we agree with the reviewer that this linearization is only a faithful approximation of the true GD dynamics for small $\epsilon$. This local regime is precisely where the concept of stability is meaningful, and it allows us to define the entire set of stable minima, which is the true object of our paper's generalization analysis. The key idea is that we are not analyzing the convergence from an arbitrary initialization, but rather characterizing the intrinsic properties of the entire set of stable minima that GD could possibly find.

Ref: [Qiao et al 24] Qiao et al Stable Minima Cannot Overfit in Univariate ReLU Networks: Generalization by Large Step Sizes, Neurips 2024

A small comment, you also say on line 29 that ReLU networks are thought not to benignly overfit, but there are quite a number of papers now showing this exact phenomenon for $d$ growing like $n^2$ (or even $n$ for leaky ReLU networks) instead of being fixed. In addition, the specialization of neurons (i.e., sparse activations) is often how one shows benign overfitting.

Our work is motivated by the observation that standard ReLU networks typically do not benignly overfit in fixed-dimensional, noisy regression. The recent work by [Haas et al. 2023] supports this, showing that benign overfitting in this setting requires modifying the ReLU activation with special high-frequency functions. In their paper, Theorem 4 explicitly confirms that standard ReLU NTKs are inconsistent when overfitting and Figure 1c empirically shows that the standard ReLU network exhibits harmful overfitting. This suggests the phenomenon is not intrinsic to standard ReLU. We would be grateful if the reviewer could share specific citations that show otherwise.

Ref: [Haas et al. 2023]: Haas et al, Mind the spikes: benign overfitting of kernels and dneural networks in fixed dimension. Neurips 2023

Can you provide some intuition as to the origin of the weighting function given in (4)..." & "What is special about ReLU here? Would similar results also hold for Leaky-ReLU, what is stopping you extend your techniques to other activation functions?

These are great questions! We would answer these two questions by explaining their common mechanism organically.

The weight function $g$ arises from two complementary perspectives: a static analysis of the loss landscape's curvature and a dynamic view of representation learning.

From a static viewpoint, as derived in [Nacson et al. 2023, Appendix F], the form of $g( {u},t)$ emerges directly from a careful analysis of the loss function's Hessian matrix for a two-layer ReLU network. Intuitively, a neuron defined by the hyperplane ${w}^{\top} {x} - b = 0$ contributes to the Hessian based on how many data points it activates. Hyperplanes that are "on the boundary" of the data support (i.e., activate very few data points) contribute less to the Hessian. Consequently, the weight $g( {u},t)$ is small in these regions, meaning the stability condition imposes a weaker regularization on neurons that are rarely active. ReLU is analytically special because its hard-sparsity property leads to a sparse loss Hessian, which allows for a clean derivation of the weighting function $g( {u},t)$ that is central to our analysis.

From a dynamic learning perspective, one can form an intuitive picture of how these pseudo-sparse neurons might emerge. High-dimensional geometry may make it easier for neurons' activation boundary to drift to shatter the dataset and isolate a very small subset of data points for each of them. Such a neuron only receives gradients from a few nearby points and if those points are already well-fit, the local gradient on this neuron's parameters can vanish, causing it to get "stuck" nearby the boundary. The small value of $g( {u},t)$ near the data boundary can be seen as the mathematical manifestation of this hypothesized dynamic: it creates "space" within the class of stable functions for these trapped, high-magnitude neurons to exist, a possibility our lower bound construction then formalizes.

Informally speaking, our work uses the set of stable minima as a prism to understand the dynamics of GD's representation leanring. The insight gained from ReLU is a general one that should be reasonably extrapolated to other activations like Leaky-ReLU, where weak negative gradients also fail to solve the underlying problem. Our intuition is that a small leaky coefficient would yield similar results (some light experiments confirmed that), as the weak gradient in the negative domain is likely insufficient to overcome this scattered signal problem. A large coefficient (e.g., an absolute value activation) might behave differently, but a rigorous analysis is challenging. We will conduct more experiments to investigate this behavior.

Ref: [Nacson et al. 2023]: Nacson et al. The Implicit Bias of Minima Stability in Multivariate Shallow ReLU Networks, ICLR 2023.

Does your class of $\epsilon$-stable functions only include those minimizers for which the loss is twice differentiable...

We thank the reviewer for this very insightful question that touches upon a key technical aspect of our work. We'd like to highlight three points. First, even though we are motivated by stable minima, the set of interest (below Prop 2.2) covers all NNs with Hessian's largest-eigenvalue being small. This set includes not just twice-differentiable stable local / global minima, but also many other low-curvature regions that are not local / global minima. Second, non-differentiable points in two-layer ReLU networks correspond to points where the hyperplane of a ReLU unit coincides with a data point. This is a measure 0 set. If we initialize the network randomly from a density and optimize using gradient descent (with constant stepsize), with probability 1, the whole trajectory will not run into such points. (see Footnote 2 in the paper). Third, analyzing twice-differentiable solutions is, in some sense, without loss of generality. This is because we can perturb the a non-differentiable solution infinitesimally to get a differentiable solution. Thus, by the Lipschitzness of the loss, it suffices that there is a differentiable solution with low-curvature near the non-differentiable solution of interest.

Also what is the main technical barrier moving to 3 layers as opposed to 2?

Moving from a two-layer to a three-layer network introduces two primary technical barriers that would require highly non-trivial extensions of our current analytical framework.

First, our analysis of the generalization gap and estimation error relies on calculating the metric entropy of some function class. For a two-layer ReLU network, the function space can be related to dictionary learning and ridge spline approximations, allowing us to leverage existing tools to bound its complexity [Parhi & Nowak 2023]. For a three-layer network, is is not even known what is the ``correct'' function class to consider. Even if the class were known, estimating the complexity of the corresponding function class would be significantly more involved and would require substantial new theoretical work.

Second, the derivation of the data-dependent weighting function is tied to a precise analysis of the loss function's Hessian matrix. For a three-layer network, the Hessian becomes much more complex.

Ref: [Parhi & Nowak 2023]: Rahul Parhi & Robert Nowak Near-minimax optimal estimation with shallow ReLU neural networks. IEEE Transcations on Information Theory, 2023

To confirm / clarify the main takeaway...

Yes, the reviewer's summary is an accurate takeaway within the setting of our paper. The implicit bias of gradient descent towards stable minima does not, on its own, prevent the curse of dimensionality, leading to a sample complexity that grows exponentially with the input dimension $d$. This is in sharp contrast to explicit regularization like weight decay.

However, we wish to add a crucial clarification: the driving force behind this phenomenon is the data distribution's interaction with the high-dimensional geometry, not just the ambient dimension in isolation. Our work highlights this "data-dependent regularity" and particularly examines the case of the uniform distribution on a ball. As mentioned in the discussion part of the paper, the induced weighted function $g$ inherits the full geometry of the data and becomes harder to describe and interpret for arbitary distributions, so we leave such an investivgation as future work.

Thank you for your support and insightful comments! Responses to your questions below:

The theoretical insight is interesting when compared to the results from relative flatness [1].

Thanks for bringing up the Relative Flatness paper [1]. We will cite and appropriately discuss that work. However, while both [1] and our paper study how flatness is connected to generalization, the settings of the two papers seem very different, and consequently, the nature of the results is different.

First and foremost, the definition of flatness in [1] seems to condition on the feature map $\phi$, i.e., the loss function is not a function of the parameter of $\phi$ and the first layer weights and biases are not learned. On the otherhand, in our paper, the loss takes both the first and second layer weights and biases as input. Our motivation is to study feature learning. Thus, in our setting, we do not consider a fixed feature map $\phi$.

Second, we focus on the overparameterized ($m>n$) regime with the hope of understanding whether "flatness/EoS/MinimaStability" can be used to explain the generalization of neural networks in this regime. In the generalization bound (Theorem 6) in [1], it seems that $m$ (this is $K$ in our paper) cannot be too large before the $n^{-2/(m+4)}$ bound becomes vacuous. Our upper bound is independent of $m$. Thus, when $m$ is sufficiently large, our bound is tighter than that of Theorem 6 of [1]. Furthermore, there was not a lower bound established in [1] that requires an exponential dependence on $m$.

Now if $m$ was kept fixed when $d$ increasing, the results from relative flatness would imply that generalization does not deteriorate, which would contradict this paper.

The reviewer is absolutely right here! Our lower bound does have an implicit assumption that $m$ is sufficiently large. If $m$ is fixed and $n\rightarrow \infty$, then we are in the classical regime and one can obtain an $O(\sqrt{dm/n})$-type uniform generalization gap bound (for all solutions and thus applicable to flat solutions as well). This bound does not have $\exp(d)$ dependence. In this under-parameterized regime, our upper bound (that depends exponentially on $d$) is still valid, but vacuous. Our exponential lower bound requires, naturally, that $m$ be sufficiently large to avoid the contradiction. A sufficient condition for our bound to kick in is that $m = \Omega(n)$.

In particular, we highlight that our Theorem 3.7 does not have this condition on $m$. This is because it deals with the weighted variation class directly, which can be viewed as a space of neural networks with arbitrary (even infinite) width and so the neural networks parameterization is abstracted away.

We will add the above discussion to the paper to increase readability.

How do the theoretical results relate to deep learning practice? There, input dimensions have grown enormously (LLMS and LVLMs), but while datasets have grown enormously as well, it seems their growth was not exponential in the input dimension.

This is a good point. We thought about this too. If the input data are truly somewhat uniformly distributed on a ball in $\mathbb{R}^d$ then our lower bound would apply directly there. However, images and texts are very structured, and we believe they are embedded in a much lower-dimensional manifold. The key to bridge the theory-practice gap could be to study how neural networks adapt to such low-dimensional structures.

Discussion on low-norm vs flat solutions is interesting. Does it only apply to MSE?

The function-space characterization of low-norm (i.e., weight decay) solutions is independent of the loss function of interest, thus the resulting generalization bounds apply to other losses. On the other hand, the function-space characterization of "flatness" in terms of the largest-eigenvalue of Hessian is sensitive to the choice of loss functions.

More specifically, we believe that the key features of low-norm vs flat solutions that we discovered in this paper do apply beyond MSE although there might be other factors to consider for other loss functions.

In particular, we studied MSE as a first step. Some parts of the techniques work beyond MSE but this would require a thorough study, e.g., if the Hessian decomposition of other twice-differentiable loss functions has a similar structure to that of the square loss. The papers that you brought up are interesting (thanks!), and we will look into them more closely to understand the difficulty of going beyond MSE as well as the merits of other weight-norm sensitive flatness definitions for future work.

Thank you for your feedback and the following question! The short answer is yes! More specifically, one localized bump functions is implemented by one ReLU atom, so the required number of bumps to satisfy Fano's lemma for the lower-bound risk is the required width of the NN, which is $K\asymp n$ (see Appendix B.4 lines 483-492 for more detail, see also Figure 9 for the visualization of the localized bump functions).

We thank the reviewer for their time and for providing the opportunity to clarify our contributions.

Clarification of the set-up of the experiment

The setup of our experiments was designed to test the precise claims of our generalization bounds. Our theory investigates the inductive bias of vanilla GD in the overparameterized regime, where generalization occurs even as network width scales with the sample size (e.g the scaling law). This is the reason why we use wide neural networks.

Furthermore, the decision to limit our experiments to $d\leq5$ is a direct and practical consequence of our theory. Our bounds predict that the required sample complexity for effective generalization scales exponentially with dimension $d$. The primary goal of our log-log plots is therefore not to achieve a low error in high dimensions, but to empirically demonstrate the trend that the generalization slope flattens as d increases. Our results show this unequivocally: the slope for $d=5$ is already exceptionally flat for a simple linear function, indicating severe sample inefficiency. Extending the experiments to higher dimensions would require a computationally infeasible number of samples merely to confirm an already established trend, providing no additional scientific insight for the immense cost involved.

Also, only linear ground truth is considered (it would be interesting to see at least one simple non-linear function, e.g. absolute value).

We have conducted experiments on non-linear functions (e.g., quadratic) and confirmed the same generalization trends, and we will add these experiments to the appendix in the final version.

Does network width affect the theoretical findings?

We focus on the overparameterized regime because the underparameterized regime does not require implicit bias for generalization. Our upper bound does apply to small $K$, but the bound could be bigger than the standard parametric rate $\sqrt{Kd/n}$. Our lower bound requires $K$ to be large to be relevant for neural networks. We will make the assumption that $K>n$ explicit in the final version.

a single ReLU hidden layer ($K=1$) should have no problem learning a noisy d-dimensional linear function, by shifting the ReLU bias to be outside of the data domain and effectively becoming a two-layer linear network in the data domain.

We agree that a 1-unit ReLU network is capable of learning a linear function. We ran new experiments and found that GD with sufficiently small learning rate does indeed finds good solutions.

Our experiments focused on the overparameterized setting (when $K$ is large) for the reason we described earlier (generalization bound via standard uniform convergence is readily available even without flatness or EoS). Our aim is not to find the best model / algorithm for fitting a linear function, but to illustrate how the modern regime of overparameterize + GD can be fragile even in such a simple setting.

Interestingly, we find that when $K = 16, d = 64$ and $n=1024$, the "pseudo-sparsity" phenomenon is still very prominent. We suspect that even a mild overparameterization (over the number of parameters to describe the target function) is rather harmful, if we solely rely on "flatness / implicit bias of GD stability".

This finding allows us to reiterate our core research paradigm. Our work uses the set of stable minima as a prism to understand the dynamics of GD's representation leanring, since directly analyzing large-step-size gradient dynamics is notoriously difficult. Within this paradigm, our lower bound construction (see Append B.4 Figure 9 for the visualization of our lower bound construction) provides a intuition: it proves the existence of a "hard-to-learn" solution built from pseudo-sparse neurons, hinting at a dynamic where neurons pathologically shatter the dataset and get 'stuck' fitting local points. For example, neurons whose activation boundary is near the edge of the data support only receives gradients from a few nearby points and if those points are already well-fit, the local gradient on this neuron's parameters can vanish, causing it to get "stuck" nearby the boundary. The higher the dimension is, the easier for these neuron to shatter the dataset.

We agree that the width of an NN plays a role on the probability of the occurance of this dynamic. As the reviewer astutely points out, such a "stuck" state is unlikely with just a single neuron. However, once representational redundancy is introduced (e.g K=16, d=64), the learning task may become distributed on various neurons and thus create the possibility for a pseudo-sparse neuron to become stable, as other neurons can compensate for the global loss, effectively shielding the "stuck" neurons from further updates. Moreover, our empirical results support this perspective (even with a narrow NN). We thank the reviewer for their insight and will clarify this perspective in the final version.

in recent years, GELU is increasingly used instead of ReLU.

We acknowledge the reviewer's point on GELU. Following the gradient analysis perspective, we believe this problem is not unique to ReLU. While GELU provides a non-zero gradient for negative inputs, this signal is weak and decays quickly. It is therefore plausible that this weak gradient is insufficient to pull a "stuck" neuron back from the data boundary. To test this, we conducted new experiments, and our empirical results show that the pseudo-sparsity and poor generalization phenomena indeed persist for GELU networks. We will add these supporting experiments to the appendix in the final version.