Stable Minima Cannot Overfit in Univariate ReLU Networks: Generalization by Large Step Sizes

Authors: dan  qiao
kaiqi  zhang
esha  singh
daniel  soudry
yu-xiang  wang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 20:25:52

Abstract: We study the generalization of two-layer ReLU neural networks in a univariate nonparametric regression problem with noisy labels. This is a problem where kernels (\emph{e.g.} NTK) are provably sub-optimal and benign overfitting does not happen, thus disqualifying existing theory for interpolating (0-loss, global optimal) solutions. We present a new theory of generalization for local minima that gradient descent with a constant learning rate can \emph{stably} converge to.  We show that gradient descent with a fixed learning rate $\eta$ can only find local minima that represent smooth functions with a certain weighted \emph{first order total variation} bounded by $1/\eta - 1/2 + \widetilde{O}(\sigma + \sqrt{\mathrm{MSE}})$ where $\sigma$ is the label noise level, $\mathrm{MSE}$ is short for mean squared error against the ground truth, and $\widetilde{O}(\cdot)$ hides a logarithmic factor. Under mild assumptions, we also prove a nearly-optimal MSE bound of $\widetilde{O}(n^{-4/5})$  within the strict interior of the support of the $n$ data points. Our theoretical results are validated by extensive simulation that demonstrates large learning rate training induces sparse linear spline fits. To the best of our knowledge, we are the first to obtain generalization bound via minima stability in the non-interpolation case and the first to show ReLU NNs without regularization can achieve near-optimal rates in nonparametric regression.

Main paper: 
 summary of contributions: 
 In this paper, we present a new theory of generalization for solutions that gradient descent (GD) with a fixed learning rate can stably converge to. Specifically:1. We show that for 1D nonparametric regression (n data points with noisy labels), the solutions that GD can stably converge to must be regular functions with small (weighted) first-order total variation (Theorem 4.1 and Corollary 4.2), thus promoting sparsity in the number of linear pieces. This generalizes the result of Mulayoff et al. [2021] by removing the "interpolation" assumption.We also show that in the noisy case, there is no "flat" interpolating solution and gradient descent cannot converge to them unless the learning rate is O(1/n 2 ) no matter how overparameterized the two-layer ReLU network is (Theorem 3.1).2. We show that such solutions (stable local minima that GD converges to) cannot overfit, in the sense that the generalization gap vanishes as n → ∞ inside the strict interior of the data support (Theorem 4.3). Moreover, under a mild additional assumption on gradient descent finding solutions with training loss smaller than σ 2 , we prove that these solutions achieve near-optimal rate for estimating (the strict interior of) first-order bounded variation functions (Theorem 4.4) -provably faster than any kernel ridge regression estimators, including neural networks in the "kernel" regime.3. We conduct extensive numerical experiments to demonstrate our theoretical predictions, validate our technical assumptions, and illustrate the functional form of the ReLU NNs as well as the learned basis functions that gradient descent finds with different step sizes. These results reveal new insights into how gradient descent training aggressively learns representation and induces implicit sparsity.To the best of our knowledge, these results are new to this paper. We emphasize that (1) the training objective function is not explicitly regularized;(2) we do not early-stop training in favor of algorithmic-stability; and (3) the solutions that gradient descent converges to are not global optimal (interpolating) solutions unless the label noise is 0. Our approach is new in that we directly analyze the complexity of a superset of solutions that gradient descent can stably converge to, which enables us to prove end-to-end generalization bounds that are near-optimal in nonparametric regression tasks.Our analysis for gradient descent-training is categorically different from those in the "kernel" (a.k.a. "lazy") regime about interpolating solutions. Instead, we rigorously prove (and empirically demonstrate) that large step-size gradient descent do not behave this way and it does not converge to interpolating solutions. Our results fall into the non-kernel regime of neural network learning known as the "rich" (a.k.a. the "feature learning" or "sparse") regime [Chizat et al., 2019 , Woodworth et al., 2020] , in which the weights and biases can move arbitrarily far away from their initialization.Technical novelty. The main technical innovation in our analysis is in handling the 1 n i (y if θ (x i ))∇ 2 θ f θ (x i ) term that arises in the minima stability analysis when it was previously handled by Mulayoff et al. [2021] using interpolation, i.e., y i = f (x i ) ∀ i ∈ [n] . It turns out in the noisy-label case, we can decompose the term into a certain Gaussian complexity measure and a self-bounding style MSE of that f θ . A non-trivial step is to bound the largest eigenvalue of ∇ 2 f θ (•) by a constant which results in a uniform bound of both the Gaussian complexity and the MSE. These bounds themselves are vacuous in terms of the implied generalization error, but plugging them into the function-space constraint imposed by the noisy minima stability bound restricts the learned ReLU NN to be inside a weighted TV1 function class (Details in Theorem 4.1). This, in turn, allows us to amplify a vacuous MSE bound into a new MSE bound that is nearly optimal in the strict interior of the data support. The main technique in the last step involves bounding the metric entropy of the weighted TV1 class and then carefully working out a self-bounding (square loss) version of the Dudley's chaining argument [Wainwright, 2019] .It is important to note that while we analyze gradient descent training of overparameterized neural networks, the computational claim is very different from those in the kernel regime. The analysis in the kernel regime ensures that gradient descent finds interpolating solutions efficiently. We do not have a comparable efficiency claim. While computational guarantees on stationary point (and local minima) convergence in non-convex optimization problems are wellunderstood [Ghadimi and Lan, 2013, Jin et al., 2017] , we do not have guarantees on whether the solution that gradient descent finds satisfies our assumption on the training loss being "optimized" (smaller than label noise σ 2 ). Instead, our results provide a generalization gap bound for any stable solutions (Theorem 4.3) and a near-optimal excess risk (MSE vs ground truth) bound (Theorem 4.4) when the solution that GD finds happens to satisfy the assumption (empirically it does!). This is a meaningful middle ground between classical learning theory which does not concern optimization at all and modern theory that is fully optimization-dependent.While we focus on (full batch) gradient descent for a clean presentation, the minima stability for stochastic gradient descent is immediate under the stochastic definition of minima stability [Mulayoff et al., 2021] . The same reason applies to us focusing on univariate nonparametric regression. Our technique can be used to generalize the multivariate function space interpretation of minima stability from Nacson et al. [2022] to the noisy case, but it will take substantial effort to formalize the corresponding generalization bounds in the multivariate case, which we leave as a future work.Generalization in Overparameterized NNs, Interpolation, and Benign Overfitting. Most existing theoretical work on understanding the generalization of overparameterized neural networks focuses on the interpolation regime [Cao and Gu, 2019 , Frei et al., 2022 , Kou et al., 2023 , Buzaglo et al., 2024] . Handling label noise either requires explicit regularization [Hu et al., 2020, Zhang and Wang, 2022] , algorithmic stability through early stopping [Hardt et al., 2016, Richards and Kuzborskij, 2021] , or carefully crafted data distribution that leads to a phenomenon known as benign overfitting [Bartlett et al., 2020 , Frei et al., 2022 , Kou et al., 2023] . Benign overfitting could happen for nonparametric regression tasks [Belkin et al., 2019] , but there is well-documented empirical and theoretical evidence that benign overfitting does not occur for regression tasks with ReLU activation [Mallinar et al., 2022 , Haas et al., 2023 , Joshi et al., 2023] and that the excess risk is required to be proportional to the standard deviation of the label noise [Kornowski et al., 2024] . We are the first to go beyond the interpolation regime and show that gradient descent-trained neural networks generalize in noisy regression tasks without explicit regularization.Implicit bias of gradient descent. The implicit bias of gradient descent training of overparameterized NN is well-studied. It was shown that among the many globally optimal (interpolating) solutions, gradient descent finds the ones with the smallest norm in certain Hilbert spaces [Arora et al., 2019 , Mei et al., 2019] , classifiers with largest-margin [Chizat and Bach, 2020] , or the smoothest cubic spline interpolation [Jin and Montúfar, 2023] . None of these results, however, imply generalization bounds when the labels are noisy. Interestingly, our results show that gradient descent with a large step size induces an implicit bias that resembles sparse L1-regularization rather than the dense L2 regularization from gradient flow [Jin and Montúfar, 2023] .Implicit bias of minima stability. The closest to our work is the line of work on the implicit bias of minima stability [Ma and Ying, 2021 , Mulayoff et al., 2021 , Wu and Su, 2023] . We build directly on top of the function-space interpretation of minima stability established by Mulayoff et al. [2021] . However, these works critically rely on the minima interpolating the data, which makes their results inapplicable to settings with label noise. Mulayoff et al. [2021] also did not establish formal generalization bounds. Ma and Ying [2021] , Wu and Su [2023] do have generalization bounds, but (again) their results require interpolation and thus do not apply to our settings.Flat/Sharp Minima and generalization. Our work is also connected to the body of work on the hypothesis that "flat local minima generalize better". Despite compelling empirical evidence [Hochreiter and Schmidhuber, 1997, Keskar et al., 2017] , rigorous theoretical understanding of this hypothesis is still lacking [see, e.g., Wu and Su, 2023, and the references therein] . Our work contributes to this literature by formally proving that the hypothesis is real for two-layer ReLU NNs in a noisy regression task.Edge-of-Stability and Catapults. Empirical observations on how large learning rate training of NN finds solutions with Hessian's largest eigenvalue dancing around 2/η, i.e., "edge of stability" regime [Cohen et al., 2020] ; and that the loss may go up first before going down to a good solution ("catapult") [Lewkowycz et al., 2020] . Existing theoretical understanding of these curious behaviors of GD training is still limited to toy-scale settings (e.g., Arora et al. [2022] , Ahn et al. [2023] , Kreisler et al. [2023] ). Our work is complementary in that we provide generalization bound to the final solution GD stabilizes on no matter how GD gets there. Outside the context of GD and neural networks, "edge of stability" and the implicit bias of large step-size were observed for forward stagewise regression [see, e.g. Tibshirani, 2015, 2014, Section 4.4 and Page 42] albeit only empirically. Our results may provide a theoretical handle in formally analyzing these observations. Optimal rates of NNs in nonparametric regression. Finally, it was previously known that neural networks can achieve optimal rates for estimating TV1 functions [Suzuki, 2018 , Liu et al., 2021 , Parhi and Nowak, 2021 , Zhang and Wang, 2022] . Specifically, Savarese et al. [2019] , Ongie et al. [2020] , Parhi and Nowak [2021] , Zhang and Wang [2022] show that weight decay in ReLU networks is connected to total variation regularization. However, these works assume one can solve an appropriately constrained or regularized empirical risk minimization problem. Our work is the first to show that the optimal rate is achievable with gradient descent without weight decay. In fact, it was a pleasant surprise to us that both weight decay and large learning rate induce total variation-like implicit regularization in the function space. Two-layer neural network. We consider two-layer (i.e. one-hidden-layer) univariate ReLU networks, 
notations and problem setup: 
 EQUATIONwhere the network consists of k hidden neurons and ϕ(•) denotes the ReLU activation function.Training data and loss function. The training dataset is denoted byD = {(x i , y i ) ∈ R × R, i ∈ [n]}. {x i } n i=1is assumed to be supported by [-x max , x max ] for some constant x max > 0. We focus on regression problems with square loss ℓ(f, (x, y)) = 1 2 (f (x) -y) 2 . The training loss is defined asEQUATION1:k , w(2)1:k , b (2) ] ∈ R 3k+1 . As a short hand, we define ℓ i (θ) := ℓ(f θ , (x i , y i )) and L(θ) := 1 n i∈[n] ℓ i (θ).Gradient descent. We focus on the Gradient descent (GD) learner, which iteratively updates θ:EQUATIONwhere η > 0 is the step size (a.k.a. learning rate) and θ 0 is the initial parameter. Detailed calculation of gradient for two-layer ReLU networks is deferred to Appendix E. Below we define stability for local minima and discuss the conditions for a minimum to be stable.Twice differentiable stable local minima. Similar to Mulayoff et al. [2021] , we consider twice differentiable minima. According to Taylor's expansion around a twice differentiable minimum θ ⋆ ,EQUATIONwhere ∇ 2 L denotes the Hessian matrix and ∇L(θ ⋆ ) = 0. Therefore, as θ t gets close to θ ⋆ , the update rule for GD (2) can be approximated asθ t+1 ≈ θ t -η ∇L(θ ⋆ ) + ∇ 2 L(θ ⋆ )(θ t -θ ⋆ ). Such approximation motivates the definition of linear stability, which is first stated in Wu et al. [2018] . Definition 2.1 (Linear stability). With the update ruleθ t+1 = θ t -η ∇L(θ ⋆ ) + ∇ 2 L(θ ⋆ )(θ t -θ ⋆ ) , a twice differentiable local minimum θ ⋆ of L is said to be ϵ linearly stable if for any θ 0 in the ϵ-ball B ϵ (θ ⋆ ), it holds that lim sup t→∞ ∥θ t -θ ⋆ ∥ ≤ ϵ.Note that different from previous works [Wu et al., 2018 , Mulayoff et al., 2021] , we remove the expectation before ∥θ t -θ ⋆ ∥ since under GD everything is deterministic. Intuitively speaking, linear stability requires that once we have arrived at a distance of ϵ from θ ⋆ , we end up staying in the ϵ-ball B ϵ (θ ⋆ ). It is known that linear stability is connected to the flatness of the local minima. Lemma 2.2. Consider the update rule in Definition 2.1, for any ϵ > 0, a local minimum θ ⋆ is an ϵ linearly stable minimum of L if and only ifλ max (∇ 2 L(θ ⋆ )) ≤ 2 η .The implication is that the set of stable minima is equivalent to the set of flat local minima whose largest eigenvalue of Hessian is smaller than 2/η. The proof is adapted from Mulayoff et al. [2021] and we state the proof in Appendix C for completeness. When the result does not depend on ϵ (as above), we simply say "linearly stable". Throughout the paper, we overload the notation by calling a function f = f θ linearly stable function if θ is linearly stable."Edge of Stability" regime. Extensive empirical and theoretical evidence (Cohen et al. [2020] , Damian et al. [2024] , and see Section 1.2) have shown that the threshold of linear stability (from Lemma 2.2) is quite significant in GD dynamics: GD iterations initially tend to exhibit "progressive sharpening", where λ max (∇ 2 L(θ t )) is increasing, until finally GD reaches the "Edge of Stability", where λ max (∇ 2 L(θ t )) ≈ 2/η. We capture this phenomenon with the following definition. Definition 2.3 (Below Edge of Stability). We say that a sequence of parameters {θ t } t=1,2,••• generated by gradient descent with step-size η is ϵ-approximately Below-Edge-of-Stability (BEoS) for ϵ > 0 if there exists t * > 0 such that λ max (∇ 2 L(θ t )) ≤ 2e ϵ η for all t ≥ t * . Any θ t with t ≥ t * is referred to as an ϵ-BEoS solution.The BEoS regime provides a strong justification for the connection between the step size η and the largest eigenvalue of the Hessian matrix. It holds for all twice-differentiable solutions GD finds along the way -even if the GD does not converge to a (local or global) minimum. Empirically, BEoS is valid for both the "progressive sharpening" phase and the oscillating EoS phase for a small constant ϵ.Our goal in this paper is to understand generalization for both (twice differentiable) stable local minima (Definition 2.1) and any other solutions satisfying ϵ-BEoS (Definition 2.3), which are both subsets ofEQUATIONTo simplify the presentation, we focus on the case with ϵ = 0 w.l.o.g.1 and unless otherwise specified, a "stable solution" refers to an element of F(η, 0, D) in the remainder of the paper.For the data generation process, we will consider two settings of interest: (1) the fixed design nonparametric regression setting (with noisy labels) (2) the agnostic statistical learning setting. They have different data assumptions and performance metrics to capture "generalization".Nonparametric Regression with Noisy labels. In this setting, we assume fixed input x 1 , • • • , x n andy i = f 0 (x i ) + ϵ i for i ∈ [n], where f 0 : R → R is the ground-truth (target) function and {ϵ i } n i=1 are independent Gaussian noises N (0, σ 2 ). Our goal is find a ReLU NN f using the dataset to minimize the mean squared error (MSE):EQUATION)It is nonparametric because we do not require f 0 to be described by a smaller number of parameters, but rather satisfy certain regularity conditions. Specifically, we focus on estimating target functions inside the first order bounded variation classf 0 ∈ BV (1) (B, C n ) := f : [-x max , x max ] → R max x |f (x)| ≤ B, xmax -xmax |f ′′ (x)|dx ≤ C n ,where f ′′ denotes the second-order weak derivative of f and we define a short handTV (1) (f ) := xmax -xmax |f ′′ (x)|dx, which we refer to as the TV (1) (semi)norm of f throughout the paper. We refer readers to a recent paper [Hu et al., 2022, Section 1.2] for the historical importance and the challenges in estimating the BV functions. The complexity of such function class is discussed in Appendix D.Agnostic statistical learning and generalization gap. In this setting, we assume the n data points {(x i , y i )} n i=1 are drawn i.i.d. from an unknown distribution P defined on[-x max , x max ] × [-D, D]. The expected performance on new data points is called "Risk", R(f ) = E (x,y)∼P [ℓ (f, (x, y))].We define the absolute difference between training loss and the risk:Gen(f ) := GeneralizationGap(f ) = |R(f ) -L(f )|.We say that f generalizes if GeneralizationGap(f ) → 0 as n → ∞ with high probability. 
stable solutions cannot interpolate noisy labels: 
 A large portion of previous works studying minima stability assume the learned function interpolates the data. However, for various optimization problems, it is unclear whether there exists such an interpolating solution that is stable, especially when the number of samples n becomes large.For the nonparametric regression problem with noisy labels, we design an example where any interpolating function can not be stable. Before presenting the example, we first define the g function, which will be the weight function of the weighted TV (1) norm throughout the paper: forEQUATION)where X is drawn from the empirical distribution of the data (a sample chosen uniformly from {x j }).For various distributions of training data (e.g. Gaussian distribution, uniform distribution), most of g's mass is located at the center while g decays towards the extreme data points. The same g(x) is also applied as the weight function in Mulayoff et al. [2021] , where they derived an upper boundxmax -xmax |f ′′ (x)|g(x)dx ≤ 1 η -12 assuming f is linearly stable (Definition 2.1) and interpolating. We generalize the same upper bound to all stable solutions in F(η, 0, D) as in Theorem C.2. Below we construct a counter-example, where we can prove a contradicting lower bound of xmax -xmax |f ′′ (x)|g(x)dx for any interpolating f , thus disproving the assumption of interpolation.Counter-example. We fixx i = 2xmaxi n-1 -(n+1)xmax n-1for i ∈ [n] and f 0 (x) = 0 for any x, which implies that y i 's are independent random variables from N (0, σ 2 ).Theorem 3.1. For the counter-example, with probability 1 -δ, for any interpolating function f ,EQUATIONwhere the randomness comes from the noises {ϵ i }. Under this high-probability event, when n ≥Ω 1 ση log 1 δ, any stable solution f for GD with step size η will not interpolate the data, i.e.EQUATIONThe proof of Theorem 3.1 is deferred to Appendix F due to space limit. This result, together with Mulayoff et al. [2021, Theorem 1] , implies that gradient descent cannot converge to interpolating solutions unless η = O(1/n2 ). It also implies (when combined with Theorem 4.1) an intriguing geometric insight that all twice-differentiable interpolating solutions must be very sharp, i.e., its largest eigenvalue is larger than Ω(n 2 ) (see details in Appendix J ). Moreover, we highlight that the conclusion of Theorem 3.1 is consistent with our observation in Figure 2 (a), where the learned function tends to be smoother and would not interpolate the data as η becomes larger. Therefore, in the following discussion, we consider the case without assuming interpolation. 
main results: 
 In this section, we present the main results about stable solutions for GD (functions in F(η, 0, D)) from three aspects. Section 4.1 describes the implicit bias of stable solutions of gradient descent with large learning rate in the function space. Section 4.2 and 4.3 derive concrete generalization bounds that leverage the implicit biases in the distribution-free statistical learning setting and the non-parametric regression setting respectively. An outline of the proof of our main theorems is given in Appendix B. The full proof details are deferred to the appendix.We begin with characterizing the stable solutions for GD with step size η without the assumption of interpolation (there can be i ∈ [n] such that f (x i ) ̸ = y i ). Similar to the interpolating case, the learned stable function f enjoys a (weighted) TV (1) bound as below. Theorem 4.1. For a function f = f θ where the training (square) loss L is twice differentiable at θ,EQUATIONwhere g(x) is defined as (6). Moreover, if we assume y i = f 0 (x i ) + ϵ i for independent noise ϵ i ∼ N (0, σ 2 ), then with probability 1 -δ where the randomness is over the noises {ϵ i },xmax -xmax |f ′′ (x)|g(x)dx ≤ λmax(∇ 2 θ L(θ)) 2 - 1 2 + O σxmax • min 1, k n + xmax MSE(f ). (10)In addition, if f = f θ is a stable solution of GD with step size η on dataset D, i.e., f θ ∈ F(η, 0, D) as in (4), then we can replaceλmax(∇ 2 θ L(θ)) 2with 1 η in (9) and (10).Theorem 4.1, which we prove in Appendix G, associates the local curvature of the loss landscape at θ with the smoothness of the function f θ it represents as measured in a weighted TV (1) norm. In short, it says that flat solutions are simple. The result is a strict generalization of Theorem 1 in Mulayoff et al. [2021] which requires interpolation, i.e., L(θ) = 0. Observe that the number of neurons k does not appear in (9) and have no effect in (10) when k > n, thus the result applies to arbitrarily overparameterized two-layer NNs. Under the standard nonparametric regression assumption, ( 10) is a stronger bound that asymptotically matches the bound under interpolation [Mulayoff et al., 2021, Theorem 1] when k = o(n) and MSE(f θ ) = o(1)as the number of data points n → ∞.Another interesting observation when combining (10) with Theorem 3.1 is that for all interpolating solutions (observe that MSE(f θ ) = O(σ 2 ) w.h.p.)λ max (∇ 2 θ L(θ)) ≥ Ω(n 2 σ) -O(σ).To say it differently, all interpolating solutions are very sharp minima when the labels are noisy. This provides theoretical explanation of the empirical observation that noisy labels are harder to overfit using gradient training [Zhang et al., 2021] .Note that we leave the term L(θ) (or MSE(f )) in the TV (1) bound. Therefore, we can plug any upper bound for these terms into Theorem 4.1 for a concrete result, and below we instantiate the TV (1) bound with a crude MSE bound under the assumption that f is "optimized".Corollary 4.2. In the nonparametric regression problem with ground-truth function f 0 , for a stable solution f = f θ of GD with step size η where L is twice differentiable at θ, assume that f is optimized, i.e, the empirical loss of f is smaller than f 0 :EQUATIONwhere the randomness is over the noises {ϵ i } and O suppresses logarithmic terms of n, 1/δ.The assumption of optimized f is rather mild since, in practice, gradient-based optimizers are commonly quite effective in loss minimization (see also our experiments). With such assumption, we can derive an MSE upper bound (with high probability) of order O(σ 2 ) (details in Lemma G.5), and thus the TV (1) bound above.In some cases, we can decrease the MSE upper bound we assumed in Corollary 4.2, and use this to improve the resulting TV (1) bound ( 11). For instance, if the neural network is under-parameterized (i.e. k is smaller than n), in Appendix G.2 we derive a (high-probability) bound for MSE of order O(k/n), which implies that the last term in (11) becomes O(σx max k/n). The term vanishes if n/k is large enough, where the TV (1) bound reduces to the noiseless and interpolating case.Why would anyone care about the function space implications of stable solutions? The next theorem shows that these solutions cannot overfit (in the strict interior of the data support) without making strong assumptions on the shape of the data distribution. Theorem 4.3. Let P be a joint distribution of (x, y) supported on [-x max , x max ]×[-D, D]. Assume the dataset D ∼ P n i.i.d. For any fixed interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that with probability 1 -δ/2, g(x) ≥ c for all x ∈ I, if the function f = f θ is a stable solution of GD with step size η such that L is twice differentiable at θ and ∥f ∥ ∞ ≤ D, with probability 1 -δ (randomness over the dataset), the generalization gap restricted to I satisfiesEQUATIONwhere E I means that (x, y) is a new sample from the data distribution conditioned on x ∈ I and n I is the number of data points in D such that x i ∈ I.The proof of Theorem 4.3 is deferred to Appendix H. Briefly speaking, we show that in the strict interior of the data support, the generalization gap will vanish as the number of data points n increases. This vanishing generalization gap further implies that the expected performance on new data points is close to the (observable) training loss, i.e., the output stable solutions do not overfit.Regarding our assumptions, in addition to standard boundedness assumptions, we focus on the strict interior of the domain where g can be lower bounded (i.e. interval I). This is because for extreme data points, g(x) decays and thus imposes little constraint on the output function f . In Appendix H.1, we show that if the marginal distribution of x is the uniform distribution on [-x max , x max ], when n is sufficiently large, I can be chosen as [-2xmax 3 , 2xmax 3 ]. In this case, with high probability, I incorporates a large portion of the data points and n I = Ω(n). More illustrations about the choice of I under various data distributions are deferred to Appendix H.2.Meanwhile, the generalization gap bound has dependence 1 η on the learning rate η. Therefore, as we increase the learning rate (in a reasonable range), the learned stable solution tends to be smoother, which further implies better generalization performances.Finally, we zoom into the nonparametric regression task that we described in Section 2 where x 1 , • • • , x n are fixed and the noisy labels y i = f 0 (x i ) + N (0, σ 2 ) independently for i ∈ [n] for a ground truth function f 0 in the first order bounded variation class (see Section 2 for details). Similar to Theorem 4.3, we focus on the strict interior I of [-x max , x max ], but instead of a generalization gap bound, we prove an MSE bound against f 0 on I that nearly matches the theoretical limit. Theorem 4.4. Under the same conditions in Corollary 4.2, for any interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that g(x) ≥ c for all x ∈ I and f is optimized over I, i.e.EQUATION)where n I is the number of data points in D such that x i ∈ I.The proof of Theorem 4.4 is deferred to Appendix I. Below we discuss some interesting aspects of the result. First of all, Theorem 4.4 focuses on an interval I where g(x) can be lower bounded.In this way, we ignore the extreme data points and derive an MSE upper bound (restricted to I) of order O(n-4/5 I), which matches the minimax optimal rate for estimating BV (1) functions [see, e.g., Donoho and Johnstone, 1998, Theorem 1]. In contrast, it is well-known that neural networks in the kernel regime (and any other "linear smoothers") must incur a strictly suboptimal worst-case et al., 1990 et al., , Suzuki, 2018]] . According to the discussions in Appendix H.1, if the data follows a uniform distribution, the interval I can incorporate most of the data points where n I = Ω(n).MSE = Ω(n I -2/3 ) [DonohoMeanwhile, the MSE bound has dependence η -2/5 on the learning rate η. Such dependence is because with a larger learning rate, the learned function f will have smaller TV (1) bound, and therefore the set of possible output functions will contain fewer non-smooth functions, which implies a tighter MSE bound. However, this does not mean that a larger learning rate is always better. When η is too large, GD may diverge, and even if it does not, the set of stable solutions cannot approximate the ground truthf 0 well if xmax -xmax |f ′′ 0 (x)|dx ≫ 1 cη + σxmax c, thus failing to satisfy the "optimized" assumption. In our experiments, we verify the "optimized" assumption numerically for a wide range of η and demonstrate that by tuning learning rate η, we are adapting to the unknown TV (1) (f 0 ).Lastly, we remark that Theorem 4.4 holds for arbitrary k, even if the neural network is heavily over-parameterized (k ≫ n). The dependence 1 η + σx max on η results from the TV (1) bound in Corollary 4.2, where the term 1 η will dominate if η ≤ 1 σxmax , which is the case if the step size is not large. Furthermore, if n k is sufficiently large, we can improve the TV (1) bound in Corollary 4.2 as in Appendix G.2, and improve the mean squared errorMSE I (f ) to O σ 2 n I 4 5 xmax η 2 5accordingly.The assumptions and detailed statements for the improved MSE are deferred to Appendix I.1. 
experiments: 
 In this section, we empirically validate our claims by training a two-layer fully connected neural network with ReLU activation using gradient descent (GD) with varying step sizes. We focus on fitting a mildly overparameterized ReLU network to a simple nonparametric regression problem. The input dataset comprises of 30 equally spaced fixed design points {x i } n i=1 , where each x i ∈ [-0.5, 0.5] (n = 30). Label y i = f 0 (x i )+N (0, σ 2 ) with σ = 0.5 and f 0 (x) = (2x+1)1(x ≤ 0)+(-2x+1)1(x > 0). The two-layer ReLU network is parameterized by θ (see Section 2) with k = 100 neurons per layer. The network uses standard parameterization (scale factor of 1) and parameters are initialized randomly (see Figure 7 for the initial basis functions). Figure 3 provides further details on the learning curves and representation learning. We note that the learned representation is very different from the initialization, thus our experiments are clearly describing phenomena not covered by the "kernel" regime. In addition, it seems that all solutions that GD finds after a small number of iterations satisfy the "optimized" assumption as required in Theorem 4.4. In the appendix (Figure 8 ), we provide empirical justification for the other assumption we make about the twice-differentiability of the solutions. More experiments can be found in the appendix with more learning rate choices, as well as a discussion on the catastrophic and tempered overfitting of interpolating solutions when we adjust k. 
conclusion: 
 In this paper, we took a new look into how gradient descent-trained two-layer ReLU neural networks generalize from a lens of minima stability (and the closely related Edge-of-Stability phenomena). We focused on univariate inputs with noisy labels and showed GD with typical choice of learning rate cannot interpolate the data. We also established that local smoothness of the training loss functions implies a first order total variation constraint on the function the neural network represents, hence proving that all such solutions have a vanishing generalization gap inside the strict interior of the data support. In addition, under a mild assumption, we prove that these stable solutions achieve near-optimal rate for estimating first-order bounded variation functions. Future work includes generalization beyond 1D input, two-hidden layers, and understanding the choice of optimization algorithms.  
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [NA] Justification: The paper does not use existing assets.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
: 
 How do gradient descent-trained neural networks work? It is an intriguing question that depends on model architecture, data distribution, and optimization algorithms used for training [Zhang et al., 2021] . Specifically, in the overparameterized regime with specific random initialization of the weights, it was shown that gradient descent finds global optimal (0-loss or interpolating) solutions despite the non-convex objective function [Jacot et al., 2018 , Du et al., 2018 , Liu et al., 2022] . It was also shown that among the (many) global optimal solutions, the particular solutions that are selected by gradient descent often do not overfit despite having 0 training error [Chizat et al., 2019 , Arora et al., 2019 , Mei et al., 2019] , sometimes even if the data is noisy -a phenomenon known as "benign overfitting" [e.g., Belkin et al., 2019 , Bartlett et al., 2020 , Frei et al., 2022] . n i = 1 (f(xi) yi) 2 MSE(f) = 1 n n i = 1 (f(xi) f0(xi)) 2Figure 2: Empirical evidence of our claim. Constant step size gradient descent-trained two-layer ReLU neural networks generalize because of minima stability. The left panel shows that with increasing step size, gradient descent finds smoother solutions (linear splines) with a smaller number of knots. The middle panel illustrates our theoretical result with a numerically accurate upper bound using 1/η + O(1) of the curvature and TV1-complexity of the smooth solution. The right panel shows that tuning η gives the classical U-shape bias-variance tradeoff for overparameterized NN.What is less well-known is that interpolating solutions do overfit for ReLU neural networks (ranging from tempered to catastrophic) [Mallinar et al., 2022 , Joshi et al., 2023 , Haas et al., 2023] and the generalization bounds in the kernel regime are provably suboptimal for certain univariate nonparametric regression problems [Suzuki, 2018, Zhang and Wang, 2022] . These "exceptions" significantly limit the abilities of the kernel theory or "benign overfitting" theory in predicting the performance of an overparameterized neural network in practice.In many learning problems with noisy data, the best solutions are simply not among those that interpolate the data. For example, no interpolating solutions can be consistent in a fixed-design nonparametric regression problem. Even if interpolating solutions that satisfy benign overfitting can be found, they could be undesirable due to their "spikiness" [Haas et al., 2023] and lack of robustness [Hao and Zhang, 2024] . In addition, it was reported that when the label is noisy, it takes much longer for gradient descent to overfit [Zhang et al., 2021] . Most practical NN training would have entered the Edge-of-Stability regime [Cohen et al., 2020] or stopped before the interpolation regime kicks in.These observations motivate us to come up with an alternative theory for gradient-descent training of overparameterized neural networks that do not require interpolation. 
 
 Appendix: In this section, we illustrate a sequence of interpolating solutions that are the global optimal solutions, which is also the kernel limit in the "lazy" regime. The results are obtained by randomly initializing the weights, but solve the minimum norm solution by directly solving the least square problem (optimizing only the second layer weights.) Notice that the number of data points n = 30. When the model is barely able to interpolate (k = 30), the fitted function experiences the catastrophic overfitting. When the number of neurons k increases the interpolation solution becomes smoother and enters the tempered overfitting regime.A.3 Representation Learning of Large Learning Rate: Visualizing Learned Basis Functions.In this section, we visualize the basis functions at initialization and after training with different learning rate. We make several observations about Figure 7 . First, the learned basis functions are very different from the initialization, so a lot of representation learning is happening, in comparison to the "kernel" regime in which nearly no representation learning is happening. Second, as η gets smaller, the number of learned basis functions increases, hence increasing the number of knots in the fitted function. Third, the learned basis function displays a strong "clustering" effect in the sense that despite overparameterization, many learned basis functions end up being the same on the data support. Interestingly, they are not the same on R, we verified that they are still different outside the data support, e.g., one of the learned basis function has a knot at x < -800.A.4 Knots of the Learned ReLU NN (aka Linear Splines) and Their CoefficientsRecall that a linear spline is a continuous piecewise linear function and a two-layer ReLU NN with k neurons span the class of all linear splines with at most k knots. In Figure 8 , we visualize the locations of the knots of the linear spline that the learned ReLU neural networks represent. Recall that all ReLU NNs are linear splines, therefore the location of the knots (i.e., the change points of the linear pieces) describes the representation learning that happens. Each basis function is a ReLU function at the knot. The final ReLU NN is a linear combination of these learned basis functions. In the first panel, we plot the quantiles of the locations of the learned knots as the function of 1/η. In the second panel, we plot the sparsity of the learned coefficients in sparse L 1 and L p norm as the function of 1/η. The third panel plots the distance of the learned knots from the closest input data points. This empirically verifies that the solution that gradient descent finds at the end is a twice differentiable function w.r.t. the parameters in the sense that not a single learned knot is exactly at the input data point, thus ensuring the applicability of our theorems.In the first panel of Figure 8 , we find for large learning rate, most of the location of the learned knot is actually outside the data support on [-0.5, 0.5]. This is a somewhat surprising finding in that the mechanism of neural network learning may actually be "pushing the knot outside the data support" so they become inactive on the training data (and only the ReLU truncated 0s are active). This is a new (and very interesting) way to understand how sparsity arrives in gradient descent learning.The second panel describes the sparsification effect of the implicit biases from large learning rate, which again, indicates that the weighted TV1 constraint is indeed making the learned function sparse (in the coefficient vector). The third panel shows that despite that the learning rate gets as small as 1e -2, none of the learned basis function actually have knots coinciding with any of the input data, thus empirically justifying our assumptions on the twice-differentiability of the solutions GD finds.In this section, we outline proof ideas for the main theorems in Section 4.Proof overview of Theorem 4.1. According to direct calculation, the Hessian is given byEQUATION)Let v denote the unit eigenvector (∥v∥ 2 = 1) of1 n n i=1 (∇ θ f θ (x i ))(∇ θ f θ (x i ))T with respect to the largest eigenvalue, then it holds thatEQUATION)For the term (⋆), we connect the maximal eigenvalue at θ to the (weighted) TV (1) norm of the corresponding f = f θ . Let g(x) be defined as ( 6), Mulayoff et al. [2021, Lemma 4] shows that (details in Lemma G.1):EQUATION)For the term (#), we can bound it by the training loss L(θ) using Cauchy-Schwarz inequality and a somewhat surprising uniform upper bound of v T ∇ 2 θ f (x i )v in Lemma E.1:EQUATIONThe first conclusion (9) of Theorem 4.1 is derived by combining the inequalities above. For the second conclusion (10), note that the term (#) can be further decomposed as:EQUATION)Similar to (#), the term |(ii)| can be bounded by 2x max MSE(f θ ) using Cauchy-Schwarz inequality. Under the data-generating assumption y i -f 0 (x i ) = ϵ i ∼ N (0, σ 2 ) i.i.d., |(i)| can be bounded by a certain empirical Gaussian complexity term sup θ1 n i ϵ i h θ (x i ) with h θ (x i ) = v T ∇ 2 θ f (x i )v.The proof is complete by Lemma G.2 which bounds this Gaussian complexity term (w.h.p.) in two ways:(1) a dimension-free bound of O(σx max ) and ( 2 For a fixed ϵ-cover of the possible output function class, according to Hoeffding's inequality and a union bound, the uniform upper bound of Gen I (f ) can be bounded byO( log N ϵ /n I ) = O(ϵ -1 4 n -1 2 I) with high probability. Note that the approximation error is of order O(ϵ), therefore Gen I (f ) can be uniformly bounded over the possible output function class byEQUATIONwhere ϵ is chosen to minimize the bound. More details can be found in the proof of Theorem H.3.Proof overview of Theorem 4.4. Similar to Theorem 4.3, we can prove a constant upper bound for I |f ′′ (x)|dx, which implies that the metric entropy of the possible output function class (4) is of order ϵ -1/2 (Details in Lemma I.2). Therefore, the critical radius r is of order O(n -2 5 I ),(∇ 2 L(θ ⋆ )) ≤ 2 η .Proof of Lemma C.1. It holds thatEQUATIONwhere the first equation is from the update rule in Definition 2.1. The second equation holds because θ ⋆ is a local minimum and therefore ∇L(θ ⋆ ) = 0. As a result,EQUATIONOn one hand, ifEQUATIONwhere the second inequality is because all the eigenvalues of I -η∇ 2 L(θ ⋆ ) is bounded between [-1, 1]. Therefore, θ ⋆ is ϵ linearly stable for any ϵ.On the other hand, if θ ⋆ is ϵ linearly stable, we choose θ 0 such that θ0-θ ⋆ ∥θ0-θ ⋆ ∥ is the top eigenvector of ∇ 2 L(θ ⋆ ) and ∥θ 0 -θ ⋆ ∥ = ϵ. Then we haveEQUATIONwhich implies that lim sup t→∞ 1 -ηλ max ∇ 2 L(θ ⋆ ) t ≤ 1, and therefore λ max (∇ 2 L(θ ⋆ )) ≤ 2 η , which finishes the proof.The following Theorem C.2 is an extension of the main result in Mulayoff et al. [2021] . Recall that stable solutions refer to the functions in F(η, 0, D) as defined in Section 2. Theorem C.2 (Extension of Theorem 1 in Mulayoff et al. [2021] ). Let f = f θ be a stable solution for GD with step size η where the training loss L(f ) = 0 and L is twice differentiable at θ.EQUATION)where gEQUATION)Here X is drawn from the empirical distribution of the data (a sample chosen uniformly from {x j }).Proof of Theorem C.2. The proof of Theorem 1 in Mulayoff et al. [2021] first proves that λ max (∇ 2 L(θ)) ≤ 2 η according to the assumption that f is linearly stable, and then proves the conclusion above. Therefore, the same conclusion directly follows for the stable solutionsf = f θ in F(η, 0, D) satisfying λ max (∇ 2 L(θ)) ≤ 2 η .Some discussions about the result. Mulayoff et al. [2021] studied the problem assuming that the optimization converges to a global minimum and the learned function interpolates the training data, i.e. L(f ) = 0. In this way, they link the Hessian matrix to properties of the learned function f and show that stable solutions of GD correspond to functions whose second order derivative has a bounded weighted norm. Moreover, as the learning rate increases, the set of stable solutions contains less and less non-smooth functions. For a fixed learning rate, according to the curve of g(x), stable solutions tend to be smoother for instances near the center of the data distribution, and less smooth for instances near the edges. More discussions can be found in Mulayoff et al. [2021] .In this section, we first define the Besov function class. Then we recall the definition of bounded variation function class and discuss the connection between these two classes. Finally we bound the metric entropy of bounded variation function class using analysis about Besov class.Let Ω be the domain of the function class (which we omit in the definition) and ∥ • ∥ p denote the ℓ p norm. We first define the modulus of smoothness.Definition D.1. For a function f ∈ L p (Ω) where 1 ≤ p ≤ ∞, the r-th modulus of smoothness is defined by w r,p (f, t) = supEQUATION)where ∆ is defined asEQUATIONThen the norm of Besov space is defined as below.EQUATIONThen the norm of Besov space is defined as:EQUATIONFinally, a function f is in the Besov class B α p,q if ∥f ∥ B α p,q is finite. For more discussions and properties of Besov class, we refer interesting readers to Edmunds and Triebel [1996] .For the same domain Ω, recall that the bounded variation function class is defined asEQUATIONAccording to DeVore and Lorentz [1993] , bounded total variation class is closely connected to the Besov class. Specifically, for any constant B, C n , it holds thatEQUATIONNow we bound the complexity of BV (1) (1, 1), which will be helpful for bounding the complexity of BV (1) (B, C n ) in the future. We first define the metric entropy of a metric space (T, ρ).Definition D.3. For a set T with a corresponding metric ρ(•, •), let N (ϵ, T, ρ) denote the ϵ-covering number of T under metric ρ. Then the metric entropy of T with respect to ρ is log N (ϵ, T, ρ).More details and examples about covering and metric entropy can be found in Chapter 5 of Wainwright [2019] . Next we bound the metric entropy of a bounded subset of BV (1). Note that the ℓ ∞ metric over domainΩ is ρ ∞ (f, g) = sup x∈Ω |f (x) -g(x)|, which we denote by ∥ • ∥ ∞ for short. Lemma D.4. Assume the set T 1 = f : [-1, 1] → R 1 -1 |f ′′ (x)|dx ≤ 1, |f (x)| ≤ 1 and the metric is the ℓ ∞ distance ∥ • ∥ ∞ ,then there exists a universal constant C 1 > 0 such that for any ϵ > 0, the metric entropy of (TEQUATION)Proof of Lemma D.4. First of all, the domain Ω = [-1, 1] is a bounded set in R. Moreover, according to DeVore and Lorentz [1993] , we have BV (1) ⊂ B 2 1,∞ . Therefore, T 1 is a bounded subset of B 2 1,∞ (Ω) with both B 2 1,∞ and ℓ ∞ norm bounded by a universal constant. Therefore, the metric space (T 1 , ∥ • ∥ ∞ ) satisfies the assumptions in (the second point) of Corollary 2 in Nickl and Pötscher [2007] with r = ∞, d = 1, s = 2. Then combining the conclusion of Corollary 2 in Nickl and Pötscher [2007] and the fact that the metric entropy is upper bounded by bracketing metric entropy (Definition 4 in Nickl and Pötscher [2007] ), we finish the proof.Remark D.5. For our purpose of bounding the metric entropy of BV(1), we only consider the case where Ω is bounded. For more results regarding the metric entropy of (weighted) Besov space, please refer to Nickl and Pötscher [2007] .In this section, we calculate the gradient and Hessian matrix of f θ (x) with respect to θ. Re-call that f θ (x) = k i=1 w (2) i ϕ w (1) i x + b (1) i + b (2) where ϕ(x) = max{x, 0}. We denote θ = (w (1) 1 , • • • , w (1) k , b (1) 1 , • • • , b (1) k , w (2) 1 , • • • , w (2) k , b (2) ) T .According to direct calculation, for a given x ∈ [-x max , x max ] we have             ∇ w (1) i f θ (x) = xw (2) i 1 w (1) i x + b (1) i > 0 , ∀ i ∈ [k] ∇ b (1) i f θ (x) = w (2) i 1 w (1) i x + b (1) i > 0 , ∀ i ∈ [k] ∇ w (2) i f θ (x) = ϕ w (1) i x + b (1) i = w (1) i x + b (1) i 1 w (1) i x + b (1) i > 0 , ∀ i ∈ [k] ∇ b (2) f θ (x) = 1 (33)In this part, we calculate∇ 2 θ f θ (x) for a given x ∈ [-x max , x max ]. Below we calculate ∂ 2 f θ (x) ∂θi∂θj . First of all, if θ i = b (2) or θ j = b (2) , ∂ 2 f θ (x) ∂θi∂θj = 0. Then it remains to calculate ∂ 2 f θ (x) ∂θi∂θ ′ j where i, j ∈ [k] and θ, θ ′ ∈ {w (1) , b (1) , w (2) }. It is obvious that if i ̸ = j, ∂ 2 f θ (x)∂θi∂θ ′ j = 0. Therefore, we only calculate the case when j = i. Let δ denote the Dirac function, it holds that:                             ∂ 2 f θ (x) ∂w (1) i ∂w (1) i = w (2) i x 2 δ w (1) i x + b (1) i , ∀ i ∈ [k] ∂ 2 f θ (x) ∂w (1) i ∂b (1) i = ∂ 2 f θ (x) ∂b (1) i ∂w (1) i = xw (2) i δ w (1) i x + b (1) i , ∀ i ∈ [k] ∂ 2 f θ (x) ∂b (1) i ∂b (1) i = w (2) i δ w (1) i x + b (1) i , ∀ i ∈ [k] ∂ 2 f θ (x) ∂w (2) i ∂w (2) i = 0, ∀ i ∈ [k] ∂ 2 f θ (x) ∂w (1) i ∂w (2) i = ∂ 2 f θ (x) ∂w (2) i ∂w (1) i = x1 w (1) i x + b (1) i > 0 , ∀ i ∈ [k] ∂ 2 f θ (x) ∂b (1) i ∂w (2) i = ∂ 2 f θ (x) ∂w (2) i ∂b (1) i = 1 w (1) i x + b (1) i > 0 , ∀ i ∈ [k] (34)The gradient is generally not well-defined according to the existence of the Dirac function. However, under the assumption that f θ is twice differentiable with respect to θ (i.e. the knots of f do not coincide with x), all the Dirac functions take the value 0. In this case,         ∂ 2 f θ (x) ∂w (1) i ∂w (1) i = 0, ∀ i ∈ [k] ∂ 2 f θ (x) ∂w (1) i ∂b (1) i = ∂ 2 f θ (x) ∂b (1) i ∂w (1) i = 0, ∀ i ∈ [k] ∂ 2 f θ (x) ∂b (1) i ∂b (1) i = 0, ∀ i ∈ [k] (35)In this part, we upper bound the operator norm of the Hessian matrix. Equivalently, we upper bound v T ∇ 2 f θ (x)v under the constraint that ∥v∥ 2 = 1. We have the following lemma.Lemma E.1. Assume that f θ (x) is twice differentiable with respect to θ and x ∈ [-x max , x max ], for any v such that ∥v∥ 2 = 1, it holds thatEQUATION) Proof of Lemma E.1. Assume v = (α 1 , • • • , α k , β 1 , • • • , β k , γ 1 , • • • , γ k , ι) T ∈ R 3k+1 such that k i=1 (α 2 i + β 2 i + γ 2 i ) + ι 2 = 1.Note that the Hessian matrix ∇ 2 θ f θ (x) follows the structure:EQUATIONwhereA w (1) w (1) , A w (1) b (1) , A b (1) w (1) , A b (1) b (1) , A w (2) w (2) ∈ R k×k , A w (1) b (2) , A b (1) b (2) , A w (2) b (2) ∈ R k×1 , A b (2) w (1) , A b (2) b (1) , A b (2) w (2) ∈ R 1×k and A b (2) b (2) ∈ R are all zero matrices. Meanwhile, A w (1) w (2) , A b (1) w (2) , A w (2) w (1) , A w (2) b(1) ∈ R k×k are all diagonal matrices whose non-zero elements are between [-max{x max , 1}, max{x max , 1}]. Therefore, it holds that:EQUATIONwhere the second inequality holds because of Cauchy-Schwarz inequality. The last inequality results from x(1 -x) ≤ 1 4 .F Proof for the Counter-Example (Theorem 3.1)Theorem F.1 (Restate Theorem 3.1). For the counter-example, with probability 1 -δ, for any interpolating function f ,EQUATION)where the randomness comes from the noises {ϵ i }. Under this high-probability event, when n ≥Ω 1 ση log 1 δ, any stable solution f for GD with step size η will not interpolate the data, i.e.EQUATION)Proof of Theorem F.1. Consider any three consecutive points x j , x j+1 , x j+2 where j ∈ [n -2], note that their corresponding y's are y j , y j+1 , y j+2 which are i.i.d. Gaussian random variables N (0, σ 2 ).Then according to Mean Value Theorem, there exists a ∈ [x j , x j+1 ] and b ∈ [x j+1 , x j+2 ] such thatf ′ (a) = y j+1 -y j x j+1 -x j = n -1 2x max (y j+1 -y j ), f ′ (b) = y j+2 -y j+1 x j+2 -x j+1 = n -1 2x max (y j+2 -y j+1 ). (41)Therefore, it holds thatEQUATION)where the last equation means that y j+2 -2y j+1 + y j follows the distribution N (0, 6σ 2 ).We focus on the interval in the middle. For any x ∈ [x n/4 , x 3n/4 ], we haveEQUATION)Together with the definition of g(x) (6), we have for any x ∈ [x n/4 , x 3n/4 ], g(x) ≥ xmax 64 . Therefore, for any interpolating solutions f , it holds thatxmax -xmax |f ′′ (x)|g(x)dx ≥x max 64x 3n/4x n/4EQUATION)where G i 's are i.i.d samples from N (0, 6σ 2 ).Assume the median of |N (0, 1)| is c > 0, which is a universal constant. For any i ∈ [ n 6 ], defineEQUATIONThen we have H i = √ 6cσ with probability 1 2 . In addition, |G i | ≥ H i . According to Lemma K.1, with probability at least 1 -δ,EQUATION) for some universal constant c ′ .Together with the conclusion in Theorem C.2, for any interpolating and stable solution f , with probability 1 -δ,EQUATIONwhich does not hold when n ≥ Ω1 ση log 1 δ .G Proof for the TV (1) Bound (Theorem 4.1)We begin by calculating the gradient of empirical loss L at θ:EQUATIONwhere the detailed calculation of ∇ θ f θ (x) can be found in Appendix E. Then the Hessian is given byEQUATIONwhere∇ 2 θ f (x) is calculated in Appendix E. Let v denote the unit eigenvector (∥v∥ 2 = 1) of 1 n n i=1 (∇ θ f θ (x i ))(∇ θ f θ (x i ))T with respect to the largest eigenvalue, it holds thatλ max (∇ 2 θ L(θ)) ≥ v T ∇ 2 θ L(θ)v = λ max 1 n n i=1 (∇ θ f θ (x i ))(∇ θ f θ (x i )) T (⋆) + 1 n n i=1 (f θ (x i ) -y i )v T ∇ 2 θ f (x i )v (#) = λ max 1 n n i=1 (∇ θ f θ (x i ))(∇ θ f θ (x i )) T (⋆) + 1 n n i=1 (f 0 (x i ) -y i )v T ∇ 2 θ f (x i )v (i) + 1 n n i=1 (f θ (x i ) -f 0 (x i ))v T ∇ 2 θ f (x i )v (ii).(50)For term (⋆), we connect the maximal eigenvalue at θ to the (weighted) TV (1) norm of the corresponding f = f θ . Let the weight function g(x) be defined as ( 6), then the lemma below holds. Lemma G.1. Assume L is twice differentiable at θ and the corresponding function of θ is f , thenEQUATION)where g(x) is defined as (6).Proof of Lemma G.1. The proof of Lemma 4 in Mulayoff et al. [2021] directly proves the result for Mulayoff et al. [2021] when f is an interpolating solution.λ max 1 n n i=1 (∇ θ f θ (x i ))(∇ θ f θ (x i )) T , which is the λ max (∇ 2 θ L) in Lemma 4 ofFor the first inequality of Theorem 4.1, we directly handle the term (#) as below.EQUATIONwhere the first inequality results from Cauchy-Schwarz inequality. The second inequality is because of the uniform upper bound of v T ∇ 2 θ f (x i )v (Lemma E.1). For the second inequality of Theorem 4.1, we bound the two terms (i) and (ii). We begin with|(i)| = 1 n n i=1 v T ∇ 2 θ f (x i )v • ϵ i, which is a weighted sum of noises {ϵ i }. Lemma G.2. Assume ϵ i 's are independently sampled from N (0, σ 2 ) for some σ > 0, with probability at least 1 -δ, uniformly over all θ, v such that L is twice differentiable at θ and ∥v∥ 2 = 1,EQUATION)Proof of Lemma G.2. For the first part, according to Lemma E.1, we haveEQUATION)Since ϵ i 's are independently sampled from Gaussian distribution N (0, σ 2 ), according to concentration of Gaussian distribution and a union bound, with probability 1 -δ 2 , it holds that:EQUATION)Under this high-probability event, we haveEQUATION)For the second part, we bound the complexity ofv T ∇ 2 θ f (x i )v n i=1. Notice that θ is a function of the dataset, thus not independent to ϵ i . Also, v is a function of the dataset and θ. Our strategy is to apply an ϵ-net argument for both v and ∇ 2 θ f (x i ) for i = 1, • • • , n. We begin with considering the possibilities of {∇ 2 θ f (x i )} n i=1 . According to the detailed form of ∇ 2 θ f θ (x) in Appendix E, we have the set {∇ 2 θ f (x i )} n i=1 is fully determined by 1 w (1)j x i + b (1) j > 0 n,k i,j=1,1. Therefore it suffices to cover all the possibilities of 1 w (1)j x i + b (1) j > 0 n i=1 for all j ∈ [k].Without loss of generality, we can assume thatx 1 < x 2 < • • • < x n ,(1)j x i + b (1) j n i=1is also monotonic, which implies that there are + 1) possibilities of 1 w (1)j x i + b (1) j > 0 n i=1in total. As a result, the product space 1 w (1)EQUATIONwhere the third inequality results from the upper bound of operator norm of M (Lemma E.1). Therefore, the exact covering set of {∇ 2 θ f (x i )} n i=1 and an ϵ 6 max{xmax,1} -covering set of the unit Euclidean Ball with dimension 3k + 1 (which is exactly the domain of v) together provides an ϵ-cover of v T ∇ 2 θ f (x i )v n i=1 with respect to ∥ • ∥ ∞ . Meanwhile, an ϵ 6 max{xmax,1} -cover with respect to ∥ • ∥ 2 of the unit Euclidean Ball with dimension 3k + 1 has cardinality bounded byEQUATIONConsider a fixed stream (a i ) n i=1 in the covering set, we have|a i | ≤ 2 max{x max , 1} for all i ∈ [n].Then according to the concentration result of Gaussian distribution, with probability 1 -δ,EQUATIONWith a union bound over the ϵ-covering set of v T ∇ 2 θ f (x i )v n i=1 and conditioned on the high probability event of (55), with probability 1 -δ, uniformly over all possible ∥v∥ 2 = 1 and θ,EQUATIONFinally, by choosing ϵ = max{xmax,1} √ n , we haveEQUATIONRemark G.3. According to the Lemma G.2 above, there are two cases. When the neural network is over-parameterized, we can derive a constant upper bound for term (i). Meanwhile, if the number of neurons is smaller than the sample complexity, we can derive a tighter bound for (i) by coveringv T ∇ 2 θ f (x i )v n i=1 .For k, n such that k is polynomially smaller than n (e.g. k = n 1-α for some positive α), the term (i) will vanish if n converges to infinity.Meanwhile, the term (ii) can be bounded by the mean squared error MSE(f θ ) using Cauchy-Schwarz inequality and the uniform upper bound ofEQUATIONCombining the results above, we state the following weighted TV (1) bound of the learned function.We leave the training loss and mean squared error (MSE) in the bound, which will be handled later. Theorem G.4 (Restate Theorem 4.1). For a function f = f θ where the training (squareEQUATIONwhere g(x) is defined as (6). Moreover, if we assume y i = f 0 (x i ) + ϵ i for independent noise ϵ i ∼ N (0, σ 2 ), then with probability 1 -δ where the randomness is over the noises {ϵ i },EQUATION)In addition, if f = f θ is a stable solution of GD with step size η on dataset D, i.e., f θ ∈ F(η, 0, D) as in (4), then we can replaceλmax(∇ 2 θ L(θ)) 2with 1 η in (9) and (10).Proof of Theorem G.4. The first inequality results from plugging Lemma G.1 and ( 52) into (50). The second inequality holds by plugging Lemma G.1, Lemma G.2 and inequality (62) into (50). For the instantiation of λ max (∇ 2 θ L(θ)), the replacement holds due to the definition of F(η, 0, D) in (4).G.1 A Crude Bound for MSE and the Instantiated TV (1) Bound (Corollary 4.2)Now we prove a crude upper bound for MSE, which could instantiate a weighted TV (1) bound.Lemma G.5. Assume that the function f is optimized, then with probability 1 -δ, we haveEQUATIONProof of Lemma G.5. According to the assumption that f is optimized, it holds thatEQUATIONThen since (a + b) 2 ≤ 2 + 2b 2 always holds (AM-GM inequality), we haveEQUATIONRecall that ϵ i 's are independently sampled from N (0, σ 2 ), then according to the concentration of Gaussian distribution and a union bound, with probability 1 -δ,EQUATIONCombining the two results, with probability 1 -δ where the randomness is over the noises {ϵ i },EQUATIONFinally, Corollary 4.2 results from directly plugging Lemma G.5 into Theorem 4.1.G.2 An Improved MSE Bound and the Corresponding TV (1) BoundIn this part, we provide an improved upper bound of the mean squared error MSE(f θ ) and also the term (ii). We first make the assumption below that the parameters are from a bounded space. Assumption G.6. There exists some constant ρ > 0 such that gradient descent converges to some local minimum θ with ∥θ∥ ∞ ≤ ρ. In addition, we assume that ∥f 0 ∥ ∞ ≤ D where D > 0 is some universal constant satisfying that D ≤ kρ 2 (x max + 1).Assumption G.6 ensures that the parameter space is bounded while the ground truth function f 0 can be approximated well by some possible output function. The assumption will surely hold for some large enough constants ρ, D, which is without loss of generality. In the following analysis, we will replace D with its upper bound kρ 2 (x max + 1) to reduce the parameters in the logarithmic terms.To handle the mean squared error MSE(f θ ), we begin with an analysis on the complexity of the function class of two-layer ReLU networks with bounded parametersEQUATIONi ϕ w(1)i x + b (1) i + b (2) , ∥θ∥ ∞ ≤ ρ . Note thathere we assume that the input is from [-x max , x max ] and there exists an upper bound ρ > 0 on the parameter θ = (wEQUATIONk , b(1)EQUATIONk , wEQUATIONEQUATION)Proof of Lemma G.7. We consider the discrete function class belowEQUATIONwhere θ j is the j-th element of ( w(1)1 , • • • , w(1) k , b(1) 1 , • • • , b(1) k , w(2) 1 , • • • , w(2) k , b(2) ) T ∈ R 3k+1 and ε • Z is the set {ε • i, i ∈ Z}.Since for each element, the number of choices is bounded by 2ρ ε + 1, the total cardinality of Fε is bounded by 2ρ ε + 1 3k+1 .For each function f ∈ F ρ with corresponding parameter θ (∥θ∥ ∞ ≤ we choose function f from Fε with corresponding parameter θ such that | θj -θ j | is minimized for all j ∈ [3k + 1]. According to our definition of Fε , we have for all j ∈EQUATIONi ϕ w(1)EQUATIONi ϕ w(1)EQUATIONwhere the last inequality is by choosing ε = ϵ 5 max{xmax,1}kρ . Therefore, the ϵ-covering number N ϵ of function class F ρ with respect to ∥ • ∥ ∞ satisfiesEQUATIONwhich implies thatEQUATIONNow we provide an upper bound of MSE(f θ ) under the assumption that θ is optimized, which means that the empirical error of f θ is smaller than that of f 0 , i.e.EQUATION)Below we state the improved upper bound of mean squared error under such assumptions. Lemma G.8. Assume ϵ i 's are independently sampled from N (0, σ 2 ) for some σ > 0, if Assumption G.6 holds and f θ ∈ F ρ is optimized, then with probability at least 1 -δ, it holds thatEQUATIONEQUATION)For the optimized function f θ , we choose a function f from the ϵ-covering set in Lemma G.7 such that∥ f -f θ ∥ ∞ ≤ ϵ.Due to identical analysis as (55), with probability 1-δ 2 , max i {|ϵ i |} ≤ 2σ log 4n δ . Under this high-probability event, it holds thatEQUATION)According to Lemma G.7, the ϵ-covering number N ϵ of function class F ρ with respect to ∥ • ∥ ∞ satisfies log N ϵ ≤ 4k log 11 max{xmax,1}kρ 2 ϵ . Due to the concentration of Gaussian distribution and a union bound over the covering set, with 1 -δ 2 , for all f in the covering set,EQUATIONCombining the two high-probability events, with probability 1 -δ, we haveEQUATIONwhere the forth inequality results from Assumption G.6. Selecting ϵ = 1 k 2 ρ 2 n 2 max{xmax,1} and solving the second order inequality, it holds thatEQUATIONPlugging in the upper bound for MSE (Lemma G.8) to Theorem 4.1, we have the corollary below. Corollary G.9 (Improved version of Corollary 4.2). For a stable solution f = f θ of GD with step size η where L is twice differentiable at θ, assume Assumption G.6 holds and f ∈ F ρ is optimized, with probability 1 -δ, the function f satisfiesEQUATIONwhere the randomness is over the noises {ϵ i } and g(x) is defined as (6).Remark G.10. The additional term here O σx max k n improves over the constant additional term in Corollary 4.2 if k < n. In addition, if n k converges to infinity (e.g. k = n 1-α for some α > 0 and n is sufficiently large), the additional term could vanish.H Proof for the Generalization Gap Bound (Theorem 4.3) Recall that the generalization gap of function f is defined asEQUATIONwhere (x, y) is a new sample from the data distribution. The generalization gap measures the difference of the expected testing error and the training loss, and a small generalization gap implies that the model is not overfitting.For a stable solution f = f θ of GD with step size η where L is twice differentiable at θ, we first derive a corresponding analysis for the weighted TV (1) bound. Recall that the empirical loss is still defined asL(f ) = 1 2n n i=1 (f (x i ) -y i ) 2 .Then using the same calculation as (50), we haveEQUATIONwhere the last inequality results from Cauchy-Schwarz inequality.According to Lemma G.1, the term (⋆) satisfiesEQUATIONIn addition, the term (#) satisfies (w.l.o.g, we assume x max ≥ 1)EQUATION)Under the assumption that the learned functionEQUATIONCombining the results, the following TV (1) bound holds. Lemma H.1. Assume the data distribution satisfies that for all possible (x, y) from the distribution, |x| ≤ x max and |y| ≤ D, if the function f = f θ is a stable solution of GD with step size η such that L is twice differentiable at θ and ∥f ∥ ∞ ≤ D, then it holds thatEQUATIONwhere g(x) is defined as (6).Note that we assume there exists some interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that with probability 1 -δ/2 (randomness over the dataset), g(x) ≥ c for all x ∈ I (w.l.o.g. we assume c < 1), which further implies thatEQUATION)We base on the high-probability event above in the following discussions. Next we bound the metric entropy of the possible output function class.Lemma H.2. Define the setEQUATIONwhere O also absorbs the constant c.Proof of Lemma H.2. Define the set T 4 as:T 4 = f : [-x max , x max ] → R ∥f ∥ ∞ ≤ D, xmax -xmax |f ′′ (x)|dx ≤ 1 c 1 η - 1 2 + 2x max D .(91) Note that the metric entropy of T 3 is bounded by that of T 4 , therefore we directly prove the upper bound of log N (ϵ, T 4 , ∥ • ∥ ∞ ). Lemma D.4 ). For a fixed ϵ > 0, according to Lemma D.4, there exists a1 = f : [-1, 1] → R 1 -1 |f ′′ (x)|dx ≤ 1, |f (x)| ≤ 1 (as inEQUATIONWe defineg i (x) = 1 c 1 η -1 2 + 2x max D x max h i ( x xmax ) for all i ∈ [N ].Then g i 's are all defined on [-x max , x max ]. Obviously, we have {g i (x)} i∈[N ] also has cardinality N .For any f (x) ∈ T 4 , we define g(x) =1 1 c ( 1 η -1 2 +2xmaxD)xmax f (x • x max ) which is defined on [-1, 1].We now show that g(x) ∈ T 1 . First of all, for any x ∈ [-x max , x max ], we haveEQUATIONMeanwhile, it holds thatEQUATIONCombining the two results, we have g ∈ T 1 . Therefore, there exists some h i such that ∥g -h i ∥ ∞ ≤ ϵ 1 c ( 1 η -1 2 +2xmaxD)xmax . Since f (x) = 1 c 1 η -1 2 + 2x max D x max g( x xmax ), ∥g i -f ∥ ∞ = 1 c 1 η -1 2 + 2x max D x max ∥h i -g∥ ∞ ≤ ϵ.In conclusion, {g i } i∈[N ] is an ϵ-covering of T 4 with respect to ∥ • ∥ ∞ . Moreover, the cardinality of {g i } i∈[N ] is N , which finishes the proof. Now we provide our main result about the generalization gap (restricted to I). Below we define the generalization gap restricted to I:EQUATION)where E I means that (x, y) is a new sample from the data distribution conditioned on x ∈ I and n I is the number of data points in D such that x i ∈ I.Theorem H.3 (Restate Theorem 4.3). Let P be a joint distribution of (x, y) supported on [-x max , x max ] × [-D, D]. Assume the dataset D ∼ P n i.i.d. For any fixed interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that with probability 1 -δ/2, g(x) ≥ c for all x ∈ I, if the function f = f θ is a stable solution of GD with step size η such that L is twice differentiable at θ and ∥f ∥ ∞ ≤ D, with probability 1 -δ (randomness over the dataset), the generalization gap restricted to I satisfiesEQUATIONwhere E I means that (x, y) is a new sample from the data distribution conditioned on x ∈ I and n I is the number of data points in D such that x i ∈ I.Proof of Theorem H.3. We base on the following event that holds with probability 1 -δ/2:EQUATION)Then the output function f ∈ T 3 defined in Lemma H.2.For a fixed ϵ > 0, according to Lemma H.2, there exists an ϵ-covering set of T 3 (with respect to ∥ • ∥ ∞ ) whose cardinality N satisfies thatEQUATIONFor a fixed function f in the covering set, since the data set {(x i , y i )} xi∈I is still i.i.d. from the data distribution conditioned on x ∈ I, Hoeffding's inequality (Lemma K.3) implies that with probability 1 -δ, it holds thatEQUATIONTogether with a union bound over the covering set, we have with probability 1 -δ/2, for all f in the covering set,EQUATIONUnder such high probability event, for any f ∈ T 3 , let f be a function in the covering set such that ∥f -f ∥ ∞ ≤ ϵ. Then it holds that I Proof for the Refined MSE Bound (Theorem 4.4)EQUATIONIn this part, we base on the same conditions in Corollary 4.2, which is the weighted TV (1) upper bound. For an output stable solution f satisfying the conclusion of Corollary 4.2, we havexmax -xmax |f ′′ (x)|g(x)dx ≤ 1 η -1 2 + O(σx max )and we denote the right-hand side by S. In addition, according to the assumption that g(x) ≥ c for any x ∈ I, we further have I |f ′′ (x)|dx ≤ S c . Now we bound the complexity of the possible output function class.According to the definition of two-layer ReLU network and our assumption that ∥θ∥ ∞ ≤ ρ, it holds that:EQUATIONEQUATION) Define the set T = f : I → R |f (0)| ≤ kρ 2 + ρ, |f ′ (0)| ≤ kρ 2 , I |f ′′ (x)|dx ≤ S c .According to the inequalities above, the possible output function (if restricted to I) belongs to T. We begin with an analysis of the metric entropy of the intermediate function setT 2 . Lemma I.1. Assume the set T 2 = f : [-x max , x max ] → R f (0) = f ′ (0) = 0, xmax -xmax |f ′′ (x)|dx ≤ C 2for some constant C 2 > 0, and the metric is ℓ ∞ distance ∥ • ∥ ∞ , then there exists a universal constant C 1 > 0 such that for any ϵ > 0, the metric entropy of (TEQUATIONwhere C 1 can be chosen as the same C 1 in Lemma D.4.Proof of Lemma I.1. Let the setT 1 = f : [-1, 1] → R 1 -1 |f ′′ (x)|dx ≤ 1, |f (x)| ≤ 1 (as in Lemma D.4). For a fixed ϵ > 0, according to Lemma D.4, there exists aEQUATIONWe define gEQUATIONCombining the two results, we have g ∈ T 1 . Therefore, there exists some h i such that ∥g -h i ∥ ∞ ≤ ϵ C2xmax . Since f (x) = C 2 x max g( x xmax ), ∥g i -f ∥ ∞ = C 2 x max ∥h i -g∥ ∞ ≤ ϵ. In conclusion, {g i } i∈[N ] is an ϵ-covering of T 2 with respect to ∥ • ∥ ∞ . Moreover, the cardinality of {g i } i∈[N ] is N , which finishes the proof.With Lemma I.1, we are ready to bound the metric entropy of T.EQUATION)where S is the right-hand side of Corollary 4.2 and O also absorbs the constant c.Proof of Lemma I.2. For any function f ∈ T, it can be written as below:EQUATIONwhere g(x) = f (x) -f (0) -f ′ (0)x satisfies that g(0) = g ′ (0) = 0 and g ′′ (x) = f ′′ (x). Therefore, to cover T to ϵ accuracy, it suffices to cover the three parts to ϵ 3 accuracy with respect to ∥ • ∥ ∞ , respectively.For f (0), since |f (0)| ≤ kρ 2 + ρ, the covering number is bounded byEQUATIONFor f ′ (0)x, since |f ′ (0)| ≤ kρ 2 , the covering number is bounded byEQUATIONFinally, for g(x), since g(x) ∈ T 2 with C 2 = S c (g is extended linearly beyond the interval I), the covering number is bounded according to Lemma I.1 above.EQUATIONCombining the three parts, the metric entropy is bounded byEQUATIONwhere O also absorbs c, which is the constant lower bound of g(x).According to the metric entropy above, we are ready to provide a refined (high probability) bound for mean squared error (restricted to I). Note that we assume that the ground-truth function f 0 ∈ T, which is necessary for the mean squared error to vanish.Lemma I.3. Under the same conditions in Corollary 4.2, for any interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that g(x) ≥ c for all x ∈ I and f is optimized over I, i.e.xi∈I (f (x i ) -y i ) 2 ≤ xi∈I (f 0 (x i ) -y i ) 2, if the output stable solution θ satisfies ∥θ∥ ∞ ≤ ρ and the ground truth f 0 ∈ T, then with probability 1 -δ (over the random noises {ϵ i }), the functionEQUATIONwhere n I is the number of data points in D such that x i ∈ I.Proof of Lemma I.3. According to the assumption that f is optimized over I, we haveEQUATION)Similar to the calculation in Lemma G.8, it holds thatEQUATIONIt is obvious that the function class T is convex, together with the assumption that f 0 ∈ T, we have the function set T ⋆ := T -{f 0 } is star-shaped (details in Section 13 of Wainwright [2019] ).Note that the metric entropy of T satisfies that logN (ϵ, T, ∥ • ∥ ∞ ) ≤ O xmaxS ϵ. According to Corollary 13.7 of Wainwright [2019], the critical radius r satisfies thatEQUATION.Finally, according to Theorem 13.5 of Wainwright [2019], we have with probability 1 -δ,EQUATIONwhich finishes the proof.Finally, Theorem 4.4 is derived by plugging in the definition of S.Theorem I.4 (Restate Theorem 4.4). Under the same conditions in Corollary 4.2, for any interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that g(x) ≥ c for all x ∈ I and f is optimized over I, i.e. xi∈I (f (x i ) -y i ) 2 ≤ xi∈I (f 0 (x i ) -y i ) 2 , if the output stable solution θ satisfies ∥θ∥ ∞ ≤ ρ (for some constant ρ > 0) and the ground truth f 0 ∈ BV (1) (kρ 2 , 1 c O( 1 η + σx max )), then with probability 1 -δ (over the random noises {ϵ i }), the function f = f θ satisfiesEQUATIONwhere n I is the number of data points in D such that x i ∈ I.Proof of Theorem I.4. Note that BV(1) (kρ 2 , 1 c O( 1 η + σx max )) is a subset of T. Then the proof directly results from Lemma I.3 and S = O 1 η + σx max .I.1 The Improved Results for the Under-parameterized CaseWe assume that n/k is large enough such that the additional term in the TV (1) bound vanishes.Assumption I.5. We assume that n k is large enough such that the last term in ( 82) O(σx max k/n) ≤ 1 2 , which further implies that xmax -xmax |f ′′ (x)|g(x)dx ≤ 1 η , where g(x) is defined as (6).Assumption I.5 requires that n k is larger than some constant, which naturally holds if k = n 1-α for some α > 0 and n is sufficiently large. Under such assumption, we improve the MSE upper bound. Theorem I.6. Under the same conditions in Corollary G.9, assume that Assumption I.5 holds. For any interval I ⊂ [-x max , x max ] and a universal constant c > 0 such that g(x) ≥ c for all x ∈ I and f is optimized over I, i.e. xi∈I (f (x i ) -y i ) 2 ≤ xi∈I (f 0 (x i ) -y i ) 2 , if the output stable solution θ satisfies ∥θ∥ ∞ ≤ ρ (for some constant ρ > 0) and the ground truth f 0 ∈ BV (1) (kρ 2 , 1 cη ), then with probability 1 -δ (over the random noises {ϵ i }), the function f = f θ satisfiesMSE I (f ) = 1 n I xi∈I (f (x i ) -f 0 (x i )) 2 ≤ O σ 2 n I 4 5x max η 2 5EQUATION)where n I is the number of data points in D such that x i ∈ I.Proof of Theorem I.6. The proof is identical to Theorem I.4, with S replaced by 1 η . Remark I.7. Compared to Theorem 4.4, Theorem I.6 is better on the dependence of η by removing the additional term σx 2 max . Such improvement results from the improved TV (1) bound (Corollary G.9) and the fact that n/k is sufficiently large.J Twice-Differentiable Interpolating Solution with Noisy Labels must be "Sharp"Recall that in the counter-example, we fixx i = 2xmaxi n-1 -(n+1)xmax n-1for i ∈ [n] and f 0 (x) = 0 for any x, which implies that y i 's are independent random variables from N (0, σ 2 ). Proposition J.1. In the example above, assume f = f θ is an interpolating solution where L is twice differentiable at θ, then with probability 1 -δ, we have • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: This is a paper regarding theory of neural networks, where we believe there is no possible negative societal impact. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).To handle the case when ϵ > 0, just replace η with ηe -ϵ in all bounds in the remainder of the paper.W.l.o.g. we assume that xmax ≥ 1. If this does not hold, we can directly replace xmax in the bounds by 1.I ) using a self-bounding technique. More details about handling other parameters can be found in Appendix I. 


Decision: This paper provides new optimization and generalization results in an interesting univariate setting.  Reviewers are uniformly positive and this is an easy accept :)

Reviews: review: paper_summary: paper_summary: This paper studies the generalizability of stable local minima in univariate regression with shallow ReLU networks. Along the way, the authors provide a bound on a weighted total variation norm of networks corresponding to stable solutions which in turn provide a tighter generalization bound.
strength_weakness: strength_weakness: It is not clear to me why the experiments (e.g. figure 3) corroborate the theoretical results in the paper. 

To illustrate my point, suppose for a given $\eta > 0$, a solution $\theta^*$ is linearly stable. 
Now suppose one trains a ReLU NN with a different learning rate $\alpha > 0$ and it converges to $\theta^*$. When considering the generalization gap of the solution $\theta^*$ (e.g. theorem 4.4), only the constant $\eta$ is relevant. The point is that the actual learning rate of GD that was used to find stable solutions is irrelevant to the generalizability of the solution (at least in the context of theorem 4.4) . Hence, I feel that more compelling experiments would study generalization errors of stable local minima that are stable for varying $\eta$. 

In this vein, I feel that the statement, “Meanwhile, the dependence $\frac 1 {\eta}$ ...” (line 296 -298) to be inexact as it is only supported by the experiments.
questions: questions: In the context of theorem 4.3 and 4.4, should $n_I$ be the number of data points belonging to the interval $I \subset [-x_{max}, x_{max}|$ instead of the length of the interval $|I|$? Otherwise, the RHS of the bound in theorem 4.4 seems to be non-vanishing as $x_{max} $ is at least on the order of the length of the interval $|I|$.  

Can the authors also provide some intuition on how the weight function $g$ was chosen?
limitations: limitations: The authors have adequately addressed the limitations.
paper_summary: paper_summary: The paper studies uni-variate regression with two layer networks and builds on the following observation: If the basin (derived form a quadratic approximation) around a given minimum is too narrow, gradient descent with a fixed step size $\eta$ will escape it, only sufficiently wide basins can capture the iterates of GD with large step sizes. It is thus reasonable to restrict the function space attainable by GD to the set of functions that can be expressed by a minimum with a wide basin. By extending a previous result to the noisy-label setting, they show that the width of the basin (or the hessian of the loss) relates to a (weighted) measure of the TV norm of the function encoded by the minimum. The function space of interest thus becomes a space functions with low TV norm. By computing covering numbers of this space, uniform generalization bounds are derived for this restricted space of functions.
strength_weakness: strength_weakness: - *The "optimized" assumption*: This assumption which effectively assumes that there exists wide and low minima is somewhat justified in lines 267-270, but I believe it is reasonable to disagree with the statement that it is mild. The experiments provided as evidence for the assumption are conducted with really large noise levels sigmas that are of the same order as $f$. The authors empirical argument could be strengthened by considering noise levels that do not dwarf the signal. It seems difficult to believe that a smooth $f$ can attain a training loss that is lower than $f_0$ when no assumption on the smoothness of $f_0$ is made. The authors should clarify why a condition on $f_0$ is not necessary for this assumption to be mild.

- The refined bounds with underparametrization: The bounds are refined by considering under parametrized settings with the width $k$ being smaller than the number of data points but the authors do not discuss its impact on the "optimized" assumption. Why would an underparametrized network be able to attain train losses smaller than a ground truth on which no assumptions are made ? I believe there are trade-offs between $\eta$, the level of underparametrization and the attained train loss that are not mentioned sufficiently by the authors.

- A final minor weakness: The work extends [Mulayoff et al 2021] in limited ways but does not exploit the added generalization to derive new conclusions. That stepsize affects the number of knots was explored before. There is now a sigma appearing in the smoothness bounds and an MSE term, yet there are no comments on what is added by their extension. The interplay between sigma and eta is not discussed, for instance, does having a large eta  but high noise means that stable functions need no have a low TV norm or is the upper bound loose?  I believe such discussions could strengthen their work as some results, namely Thm 4.3, are straightforward applications of generalization bounds for  spaces with known metric entropy.
questions: questions: The weight function $g$ is inherited from [Mulayoff et al 2021]  but would really benefit from some more clarifications, even for completeness. It is clear that the interval $\mathcal{I}$ must be introduced to remove the weighing and allow the authors to use covering numbers for bounded TV norm spaces. Could the authors explain why it arises in theorem 4.1 ?
limitations: limitations: The limitations are discussed.
paper_summary: paper_summary: This paper studies the generalization properties of two-layer ReLU neural networks in a univariate nonparametric regression problem with noisy labels. It proposes a new theory for local minima to which gradient descent (GD) with a fixed learning rate $\eta$ stably converges. The paper shows that GD with a constant learning rate can only find stable local minima whose weighted TV(1) is bounded by $1/\eta-1/2+\tilde{O}(\sigma+\sqrt{\mathrm{MSE}})$. With this property of local minimas, they prove the generalization bound for univariate nonparametric regression. The theoretical results are validated by extensive simulations demonstrating that large learning rate training induces sparse linear spline fits.
strength_weakness: strength_weakness: 1. The analysis in this paper has little focus on optimization. There is no rigorous theoretical evidence that GD will find the solutions that satisfies the assumptions (though I know this is an open problem in the literature).
2. This paper is limited to univariate functions, which is not important in practice. The authors claims that the technique can be generated to multivariate functions but they did not show that.
questions: questions: How large is the family of target functions satisfying the BV(1) condition?
limitations: limitations: The authors have addressed some limitations of their work, acknowledging that their analysis is only for full batch gradient descent and univariate nonparamatric regression.
paper_summary: paper_summary: This paper studies the implicit regularization of large learning rates in gradient descent. The setting is univariate linear regression with two-layer ReLU neural networks. The authors show that, if GD converges to a local minimum, then the function implemented by the neural network at this local minimum has a bounded first-order total variation. This regularity avoids overfitting, as quantified by a generalization bound. Numerical experiments validate the results.
strength_weakness: strength_weakness: The paper considers a univariate case, some comments on the extension to the multivariate case are given at the end of Section 1.1.
questions: questions: Do the authors have any insight on the role of overparameterization in the setting they consider? It does not seem to have a large influence on the generalization bounds, but perhaps it would help the optimization (e.g., via a PL inequality)?

Minor remarks:
- lines 132: the connection between stepsize and L1/L2 regularization has been thoroughly analysed for diagonal linear networks in [2].
- lines 147-153: a relevant theoretical paper on the Edge of Stability is [1].

[1] Damin, Nichani, Lee, Self-Stabilization: The Implicit Bias of Gradient Descent at the Edge of Stability, ICLR 2023.

[2] Even, Pesme, Gunasekar, Flammarion. (S)GD over Diagonal Linear Networks: Implicit Regularisation, Large Stepsizes and Edge of Stability. NeurIPS 2023.
limitations: limitations: Limitations are properly addressed.

Comments: title: Thank you
Thanks a lot! We understand what you mean now!
Thank you!
Any further questions?
Any follow-up questions before discussion period expires?
Any further questions / comments?
Thank you for the rebuttal.
comment: Thanks for acknowledging our response and for your continued support!
I thank the authors for their rebuttal and I will keep my score.
Thanks for the follow-up questions!
> In Corollary 4.2, the optimized assumption, as far as I understand, is just assumed to be attainable without $f_0$ belonging to a BV space. My issue is that this assumption is unrealistic without constraints on $f_0$ and I believe it reasonable to disagree with your paragraph 267.

I see!  You are right that Corollary 4.2 is a valid statement without assuming regularity on f_0.  In fact, for Corollary 4.2 to be valid, all we need is MSE = $\tilde{O}(\sigma^2)$. The "optimized" condition was only used for convenience because we can bound $$\sqrt{MSE} = \\|f - f_0\\| \\leq \\|f - y\\| + \\|y - f_0\\| \leq 2\sigma.$$

We can instead just use $$\sqrt{MSE} \leq \sigma + \sqrt{\text{TrainingLoss}}$$ so Corollary 4.2 does not need the "optimized" assumption.  TrainingLoss being a constant is relatively mild, because all 0 initialization already have a training loss of $\frac{1}{2n}\sum_i y_i^2 = O(1)$ if label $y_i$ are all bounded. If gradient descent does not diverge, then it should find solutions with lower loss than initialization.

It is indeed unnatural to assume "optimized" for Corollary 4.2 and the discussion about it there is out of place.  We propose to defer that to Section 4.3 where it is actually used --- for solving a nonparametric regression task under the assumptions of the regularity of f_0.

> The interplays between $\sigma$ and $\eta$ (and underparametrization discussions): I believe it is very essential to include these discussions and experiments. 

Good call. We can add these experiments and discussions.   We focused on discussing the regime when $\sigma$ is a constant because such a low signal-to-noise ratio setting is the conventional setting for nonparametric regression. 

Our theorems are not asymptotic and they work for all $\eta>0, \sigma>0$. One more technical result we can explicitly state about the interplay between $\sigma$ and $\eta$  is to expose $\sigma$ as a parameter in Theorem 3.2.  The updated statement will say that interpolation is not possible unless $\eta < 1/(\sigma n^2)$.
Thank you for your clarifications.

*The optimized assumption*: In Corollary 4.2, the optimized assumption, as far as I understand, is just assumed to be attainable without $f_0$ belonging to a BV space. My issue is that this assumption is unrealistic without constraints on $f_0$ and I believe it reasonable to disagree with your paragraph 267. 

I will restate my concern in other words: there is, in my understanding, a problem of having an unconstrained $\eta$: your result states that for a given $\eta$, \emph{if the ground truth is sufficiently$\frac{1}{\eta}$ smooth, hence making the optimized assumption achievable, then Corollary 4.2 holds. In your result $\eta$ needs to have an upperbound set by the smoothness of $f_0$ otherwise, the optimized assumption is an act of faith. This need for an upperbound contrasts with the spirit of the paper which aims to analyze "large step-sizes". I would appreciate if the authors could clarify why they believe the optimized assumption in Corollary 4.2 is mild. From my perspective, the optimized assumption is implicitly disallowing large stepsizes.

*The interplays between $\eta$ and $\sigma$* (and underparametrization discussions): 
I believe it is very essential to include these discussions and experiments. The extension from [Mulayoff et al 2021] would be too limited otherwise.
Thanks again for your high-quality review and your support.
Thank you for the detailed response, and I have raised the score accordingly.
Thanks for acknowledging our response and for your continued support!
Thanks again for your time in reviewing our paper! We have addressed your technical questions above and shared some perspectives.  Could you kindly let us know if our rebuttal satisfactorily resolved your concerns? 

We would love an opportunity to address any further questions and comments you may have before the author discussion period expires.
Thanks again for your time in reviewing our paper! We have addressed your technical questions above.  Could you kindly let us know if our rebuttal satisfactorily resolved your concerns? 

We would love an opportunity to address any further questions and comments you may have before the author discussion period expires.
Thanks again for your time in reviewing our paper! We have addressed your technical questions above.  Could you kindly let us know if our rebuttal satisfactorily resolved your concerns? 

We would love an opportunity to address any further questions and comments you may have before the author discussion period expires.
I thank the authors for their rebuttal and keep my score.