Shallow Univariate ReLU Networks as Splines: Initialization, Loss Surface, Hessian, and Gradient Flow Dynamics

1. Introduction

Deep learning has revolutionized the field of machine learning (ML), leading to state of the art performance in image segmentation, medical imaging, machine translation, chess playing, reinforcement learning, and more. Despite being intensely studied and widely used, theoretical understanding of some of its fundamental properties remains poor. One critical, yet mysterious property of deep learning is the root cause of the excellent generalization achieved by overparameterized networks (Zhang et al., 2016).

In contrast to other machine learning techniques, continuing to add parameters to a deep network (beyond zero training loss) tends to improve generalization performance. This has even been observed for networks that are massively overparameterized, wherein, according to traditional ML theory, they should (over)fit the training data (Neyshabur et al., 2015). This overparameterization leads to a highly complicated loss surface, with increasingly many local minima. How does training networks with excess capacity lead to generalization? And how can generalization error decrease with overparameterization? How can gradient descent in a loss landscape riddled with local minima so frequently converge to near-global minima? One possible solution is a phenomenon known as implicit regularization (or implicit bias). Empirical studies (Zhang et al., 2016; Advani et al., 2020; Geiger et al., 2020) have shown that overparametrized NNs behave as if they possess strong regularization, even when trained from scratch with no explicit regularization—instead possessing implicit regularization or bias (IR or IB), since such regularization is not explicitly imposed by the designer.

We approach these issues by considering symmetries of the NN parameterization. Several recent works have brought attention to the importance of symmetries to deep learning (Badrinarayanan et al., 2015; Kunin et al., 2020; Tayal et al., 2020). One important advance is building symmetry aware architectures that automatically generalize, such as convolutional layers granting translation-invariant representations, or many other task-specific symmetries (Liu et al., 2016; Barbosa et al., 2021; Bertoni et al., 2021; Liu and Okatani, 2021). Understanding and exploiting symmetries may lead to similar improvements, resulting in better performance or models that require less training data due to implementation of expert-knowledge based symmetry groups.

We consider a new spline-based reparamaterization of a shallow ReLU NN, explicitly based on taking a quotient with respect to an existing weight symmetry. In particular, we focus on shallow fully connected univariate ReLU networks, whose parameters will always result in a Continuous Piecewise Linear (CPWL) output. We reparameterize the CPWL function implemented by the network in terms of the locations of the nonlinearities (breakpoints), the associated change in slope (delta-slopes), and a parameter that identifies which side of the breakpoint the associated neuron is active on (orientation). This parameterization is defined formally in section 2.1, and illustrated in Figure 2. We provide theoretical results for shallow networks, with experiments confirming these results.

1.1. Main Contributions

The main contribution of this work are as follows:

- Initialization: Increasingly Flat with Width. In the spline perspective, once the weight symmetry has been accounted for, neural network parameters determine the locations of breakpoints and their delta-slopes in the CPWL reparametrization. We prove that, for common initializations, these distributions are mean 0 with low standard deviation. Notably, the delta-slope distribution becomes increasingly concentrated near 0 as the width of the network increases, leading to flatter initial approximations.

- A characterization of the Loss Surface, Critical Points and Hessian, revealing that Flat Minima are due to degeneracy of breakpoints. We fully characterize the Hessian of the loss surface at critical points, revealing that the Hessian is positive semi-definite, with 0 eigenvalues occurring when elements of its Gram matrix set are linearly dependent. We show that many of these 0 eigenvalues occur due to symmetry artifacts of the standard NN parametrization. We characterize the ways this can occur, including whenever multiple breakpoints share the same active data—thus we expect that at any given critical point in an overparametrized network, the loss surface will be flat in many directions.

- Implicit Regularization is due to Flat Initialization in the Overparametrized Regime. We find that implicit regularization in overparametrized Fully Connected (FC) ReLU nets is due to three factors: (i) the very flat initialization, (ii) the curvature-based parametrization of the approximating function (breakpoints and delta-slopes, made explicitly clear by our symmetry-quotiented reparametrization) and (iii) the role of gradient descent (GD) in preserving initialization and providing regularization of curvature. In particular, each neuron has an global, rather than local, impact as each contributes an affine ReLU. Thus, backpropogation distributes the “work” of fitting the training data over many units, leading to regularized delta-slope and breakpoints. All else equal, these non-local effects mean that more overparametrization leads to each unit mattering less, thus typically resulting in better generalization due to implicit regularization (Neyshabur et al., 2015, 2018).

2. Theoretical Results

2.1. Spline Parametrization and Notation

Consider a fully connected ReLU neural net $\widehat{f}(x;\theta )$ with a single hidden layer of width H, scalar input x ∈ ℝ and scalar output ŷ ∈ ℝ, which we are attempting to match to a target function f(x), from which we have N sample input/output data pairs (x_n, y_n); the full vectors of inputs and outputs are denoted x and y, respectively. $\widehat{f}(·;\theta )$ is a continuous piecewise linear (CPWL) function since the ReLU nonlinearity is CPWL.

Here the NN parameters ${\theta }_{\mathrm{\text{NN}}}\triangleq \left\{b_0\right\}\cup {\left\{(w_i,b_i,v_i)\right\}}_{i=1}^H$ denote the global bias plus the input weight, bias, and output weight, respectively of neuron i, and (·)₊ ≜ max {0, ·} denotes the ReLU function.

An important observation is that the ReLU NN parametrization is redundant: for every function $\widehat{f}$ represented by Equation (1) there exists infinitely many transformations of the parameters ${\theta }_{\mathrm{\text{NN}}}^{\prime }\equiv {R}({\theta }_{\mathrm{\text{NN}}})$ s.t. the transformed function $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}}^{\prime })=\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}})$. These invariant transformations ${R}$ consist of (i) permutations of the hidden units and (ii) scalings of the weights and biases of the form $w_i\mapsto {\alpha }_i{w}_i,b_i\mapsto {\alpha }_i{b}_i,v_i\mapsto {\alpha }_i^{-1}{v}_i$ for α_i ∈ ℝ_>0 (Rolnick and Kording, 2019). The manifold of constant loss generated by the α-scaling symmetry is shown in Figure 1. The set ${G}$ of such function-invariant transformations together with function composition ◦ forms a group.

We want to understand the function implemented by this neural net, and so we ask: How do the CPWL parameters relate to the NN parameters? We answer this by taking the quotient with respect to the above α-scaling symmetry, transforming from the NN parametrization of weights and biases to the following CPWL spline parametrization:

where

and the Iversen bracket ⟦b⟧ is 1 when the condition b is true, and 0 otherwise. The CPWL spline parametrization is the Breakpoint, Delta-Slope, Orientation (BDSO) parametrization, ${\theta }_{\mathrm{\text{BDSO}}}\triangleq \left\{b_0\right\}\cup {\left\{({\beta }_i,{\mu }_i,s_i)\right\}}_{i=1}^H$, where ${\beta }_i\triangleq -\frac{b_i}{w_i}$ is (the x-coordinate of) the breakpoint (or knot) induced by neuron i, μ_i ≜ w_iv_i is the delta-slope contribution of neuron i, and s_i ≜ sgn w_i ∈ {±1} is the orientation of β_i (left for s_i = −1, right for s_i = +1). This terminology is illustrated in Figure 2.

Considering the BDSO parametrization provides a new, useful lens with which to analyze neural nets, enabling us to reason more easily and transparently about the initialization, loss surface, and training dynamics. The benefits of this approach derive from taking the quotient with respect to the α-scaling symmetry, leading to two useful properties: (1) the loss depends on the NN parameters θ_NN only through the BDSO parameters (the approximating function) θ_BDSO i.e., ℓ(θ_NN) = ℓ(θ_BDSO(θ_NN)), analogous to the concept of a minimum sufficient statistic in exponential family models; and (2) the BDSO parameters are more intuitive, allowing us to reproduce known results (Williams et al., 2019) more succinctly, and expand out to new results more easily. Much recent related work has also moved in this direction, analyzing function space (Balestriero and Baraniu, 2018; Hanin and Rolnick, 2019).

The BDSO parametrization makes clear that, fundamentally, $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}})$ is represented in terms of its second derivative. Examining

further illustrates the meaning of the BDSO parameters, especially μ_i and β_i: the second derivative is zero everywhere except at x = β_i for some i, where there is a Dirac delta of magnitude μ_i.

We note that the BDSO parametrization of a ReLU NN is closely related to but different than a traditional mth order spline parametrization ${\widehat{f}}_{\mathrm{\text{spline}}}(x)\triangleq \sum _{i=1}^K{\mu }_i{(x-{\beta }_i)}_{+}^m+\sum _{j=0}^m{c}_j{x}^j$ (Reinsch, 1967). The BDSO parametrization (i) lacks the base polynomial, and (ii) has two possible breakpoint orientations s_i ∈ {±1} whereas the spline is canonically all right-facing. Additionally, adding in the base polynomial (for the linear case m = 1) into the BDSO parametrization yields a ReLU ResNet parametrization.

It is also useful to consider the set of sorted breakpoints, ${\left\{{\beta }_p\right\}}_{p=0}^{H+1}$, where −∞ ≜ β₀ < β₁ < ⋯ < β_p < ⋯ < β_H+1 ≜ ∞, which induces two partitions: (i) a partition of the domain ${X}$ into intervals ${{X}}_p\triangleq [{\beta }_p,{\beta }_{p+1})$; and (ii) a partition of the training data ${{D}}_x\triangleq {\left\{x_n\right\}}_{n=1}^N$ into pieces ${\pi }_p\triangleq {{X}}_p\cap {{D}}_x$, so that ${{D}}_x={\cup }_{p=0}^{H+1}{\pi }_p$. Note that sorting the breakpoints (and hence the neurons) removes the permutation symmetry present in the standard NN parametrization.

Finally, we introduce some notation relevant to the training of the network: let $\widehat{ϵ}\in {ℝ}^{N\times 1}$ denote the vector of residuals $\text{y}-\widehat{\text{y}}$, i.e., ${\widehat{ϵ}}_n=f(x_n)-\widehat{f}(x_n)$, let $1_i\in {ℝ}^{N\times 1}$ denote the activation pattern of neuron i over all N inputs, $1_{i,n}=⟦w_i{x}_n+b_i>0⟧$, and let ${\widehat{ϵ}}_i\triangleq \widehat{ϵ}\odot 1_i\in {ℝ}^{N\times 1}$ and ${\text{x}}_i\triangleq \text{x}\odot 1_i\in {ℝ}^{N\times 1}$ denote the relevant residuals and inputs for neuron i, respectively. In words, relevant residuals and inputs for neuron i are those which are not zeroed out by the ReLU activation.

2.1.1. Basis Expansion, Infinite Width Limit

The expansion in Equation 2 can be further understood by viewing it as a basis expansion. Specifically, we can view the BDSO framework as a basis expansion of ${\widehat{f}}^{″}$ in a basis of Dirac delta functions, or, from Equation 1 as a basis expansion of $\widehat{f}$ in a basis of functions (x − β_i)_si. The form of these basis functions is fixed as the ReLU activation—analogous to the mother basis function in wavelets (Rao, 2002)—but can be translated via changes in β_i, reflected by s_i, or re-weighted via changes in μ_i.

Also, note that if the breakpoints β_i are fixed, and only the delta-slope parameters μ_i are optimized, then we effectively have a (kernel) linear regression model. Thus, the power of neural network learning lies in the ability to translate/orient these basis functions in order to better model the data. Intuitively, in a good fit, the breakpoints β_i will tend to congregate near areas of high curvature in the ground truth function, i.e., |f″(x)| ≫ 0, so that the sum-of-Diracs ${\widehat{f}}^{″}(x;{\theta }_{\mathrm{\text{BDSO}}})$ better approximates f″(x).¹

Consider the infinite width limit H → ∞, rewriting Equation 2 as

Using the Law of Large Numbers, ${𝔼}_{\beta }[\frac{c(\beta )}{p_{\beta }(\beta )}\delta (x-\beta )]=\frac{1}{H}{\sum }_i\frac{c({\beta }_i)}{p_{\beta }({\beta }_i)}\delta (x-{\beta }_i)+O(1/H)$ and so we may write μ_is_i ≜ c(β_i)/Hp_β(β_i). Then, from Equation 3, we can expand this as ${\mu }_i{s}_i={\widehat{f}}^{″}({\beta }_i;{\theta }_{\mathrm{\text{BDSO}}})/H{p}_{\beta }({\beta }_i)$.

2.2. Random Initialization in Function Space

We now study the random initializations commonly used in deep learning in function space. These include the independent Gaussian initialization, with $b_i~{N}(0,{\sigma }_b)$, $w_i~{N}(0,{\sigma }_w)$, $v_i~{N}(0,{\sigma }_v)$, and independent uniform initialization, with $b_i~{U}\left[-a_b,a_b\right]$, $w_i~{U}\left[-a_w,a_w\right]$, $v_i~{U}\left[-a_v,a_v\right]$. We find that common initializations result in flat functions, becoming flatter with increasing width.

Theorem 1. Consider a fully connected ReLU neural net with scalar input and output, and a single hidden layer of width H. Let the weights and biases be initialized randomly according to a zero-mean Gaussian or Uniform distribution. Then the induced distributions of the function space parameters (breakpoints β, delta-slopes μ) are as follows:

(a) Under an independent Gaussian initialization,

(b) Under an independent Uniform initialization,

Using this result, we can immediately derive marginal and conditional distributions for the breakpoints and delta-slopes.

Corollary 1. Consider the same setting as Theorem 1.

(a) In the case of an independent Gaussian initialization,

where $G_{pq}^{nm}(·|·)$ is the Meijer G-function and K_ν(·) is the modified Bessel function of the second kind.

(b) In the case of an independent Uniform initialization,

where Tri(·;a) is the symmetric triangular distribution with base [−a, a] and mode 0.

2.2.1. Implications

Corollary 1 implies that the breakpoint density drops quickly away from the origin for common initializations. As breakpoints are necessary to fit curvature, if the ground truth function f(x) has significant curvature far from the origin, then it may be far more difficult to fit. We show that this is indeed the case by training a shallow ReLU NN with an initialization that does not match the underlying curvature, with training becoming easier if the initial breakpoint distribution better matches the function curvature. This suggests that a data-dependent initialization, with more breakpoints near areas of high curvature, could potentially be faster and easier to train. Note that this simple curvature-based interpretation is made possible by our symmetry-quotiented reparametrization.

Next, we consider the typical Gaussian He (He et al., 2015) or Glorot (Glorot and Bengio, 2010) initializations. In the He initialization, we have ${\sigma }_w=\sqrt{2}$, ${\sigma }_v=\sqrt{2/H}$. In the Glorot initialization, we have ${\sigma }_w={\sigma }_v=\sqrt{2/(H+1)}$. We wish to consider their effect on the smoothness of the initial function approximation. From here on, we measure the smoothness using the roughness of $\widehat{f}$, $\rho \triangleq \int |{\widehat{f}}^{″}(x;{\theta }_{\mathrm{\text{BDSO}}}){|}^2\mathrm{\text{d}}x={\sum }_i{\mu }_i^2={\Vert \mu \Vert }_2^2$, where lower roughness indicates a smoother approximation.

Theorem 2. Consider the initial roughness ρ₀ under a Gaussian initialization. In the He initialization, we have that the tail probability is given by

where 𝔼[ρ₀] = 4. In the Glorot initialization, we have that the tail probability is given by

where $𝔼\left[{\rho }_0\right]=\frac{4H}{{(H+1)}^2}=O\left(\frac{1}{H}\right)$.

Thus, as the width H increases, the distribution of the roughness of the initial function ${\widehat{f}}_0$ gets tighter around its mean. In the case of the He initialization, this mean is constant; in the Glorot initialization, it decreases with H. In either case, for reasonable widths, the initial roughness is small with high probability, corresponding to small initial delta-slopes. This smoothness has implications for the implicit regularization phenomenon observed in recent work (Neyshabur et al., 2018), and studied later in section 2.6, in particular, Theorem 7.

Related Work. Several recent works analyze the random initialization in deep networks. However, there are two main differences; first, they focus on the infinite width case (Neal, 1994; Lee et al., 2017; Jacot et al., 2018; Savarese et al., 2019) and can thus use the Central Limit Theorem (CLT), whereas we focus on finite width case and cannot use the CLT, thus requiring nontrivial mathematical machinery (see Supplement for detailed proofs). Second, they focus on the activations as a function of input whereas we also compute the joint densities of the BDSO parameters i.e., breakpoints and delta-slopes. The latter is particularly important for understanding the non-uniform density of breakpoints away from the origin as noted above. The most closely related work is Steinwart (2019), which considers only the breakpoints, and suggests a new initialization with uniform breakpoint distribution.

2.3. Loss Surface in the Spline Parametrization

We now consider the squared loss $\frac{1}{2}\sum _{n=1}^N{(\widehat{f}(x_n;\theta )-y_n)}^2$ as a function of either the NN parameters ℓ(θ_NN) or the BDSO parameters $\tilde{\ell }({\theta }_{\mathrm{\text{BDSO}}})$. We begin with the symmetry-quotiented case:

Theorem 3. The loss function $\tilde{\ell }({\theta }_{\mathrm{\text{BDSO}}}=(\beta ,\mu ,\text{s}))$ is a continuous piecewise quadratic (CPWQ) spline. Furthermore, consider the evolution of the loss as we vary β_i along the x axis; this 1-dimensional slice $\tilde{\ell }({\beta }_i;{\beta }_{-i},\mu ,\text{s})$ is also a CPWQ spline in β_i with knots at datapoints ${\left\{x_n\right\}}_{n=1}^N$. Let p₁(β_i) (resp. p₂(β_i)) be the quadratic function equal to $\tilde{\ell }({\beta }_i;{\beta }_{-i},\mu ,\text{s})$ for β_i ∈ [x_n−1, x_n] (resp. [x_n, x_n+1)], which both have positive curvature, and let m_j ≜ arg min p_j(β_i). Then, with measure 1, the knots x_n fall into one of three types as shown in Figure 3:

- (Type I, Passthrough) m₁, m₂ < x_n, or x_n < m₁, m₂,

- (Type II, Repeller) m₁ < x_n < m₂

- (Type III, Attractor) m₂ < x_n < m₁.

We call Type I knots Passthrough knots because gradient flow starting from the higher side will simply pass through the knot. Similarly, we call Type II knots Repellers and Type III knots Attractors because β_i will be repelled by or attracted to the knot during gradient flow.

We now return to the loss ℓ(θ_NN) under the NN parametrization, and consider varying (w_i, v_i, b_i) such that the corresponding β_i changes but μ_i and s_i stay the same. Due to the α-scaling symmetry, this can be implemented in two distinct ways: (i) (w_i, v_i, b_i) ↦ (w_i, v_i, b_i + δ) or (ii) (w_i, v_i, b_i) ↦ ((1+δ)w_i, v_i/(1+δ), b_i). Version (ii) can be implemented by applying version (i) followed by the transformation (w_i, v_i, b_i) ↦ (αw_i, v_i/α, αb_i) for the appropriate α, which leaves the loss invariant. In fact, there is a continuum of transformations which could be achieved by version (i) followed by an arbitrary α-transformation. Thus, consider the 1-dimensional slice $\ell (b_i;\text{w},\text{v},{\text{b}}_{-i})$, which is equal to $\tilde{\ell }(-b_i/w_i;{\beta }_{-i},\mu ,\text{s})$, i.e., a horizontally reflected and scaled version of $\tilde{\ell }({\beta }_i;{\beta }_{-i};\mu ,\text{s})$. The 1-dimensional “slice” corresponding to version (ii) will similarly be a stretched version of $\tilde{\ell }({\beta }_i;{\beta }_{-i};\mu ,\text{s})$, with the caveat that this “slice” is along a hyperbola in the (v_i, w_i)-plane. Noting that any changes in (w_i, v_i, b_i) that implement the same change in β_i change the loss in the same way, let ℓ(β_i; θ_NN\β_i) denote the equivalence class of 1-dimensional slices of ℓ(θ_NN). This gives the following corollary:

Corollary 2. The 1-dimensional slices ℓ(β_i; θ_NN\β_i) of ℓ(θ_NN) that implement a change only in β_i are CPWQ with knots as in Theorem 3.

2.3.1. Critical Points of the Loss Surface

In addition to the fixed points associated with attractor datapoints, there are also critical points that lie in the quadratic regions, i.e., local (and global) minima which correspond to the minima of some parabolic piece for each neuron. We characterize these critical points here.

Theorem 4. Consider some arbitrary ${\theta }_{\mathrm{\text{NN}}}^{*}$ such that there exists at least one neuron that is active on all data. Then, ${\theta }_{\mathrm{\text{NN}}}^{*}$ is a critical point of $\tilde{\ell }(·)$ if and only if for all domain partitions ${{X}}_p$, the restriction $\widehat{f}(·;{\theta }_{\mathrm{\text{NN}}}^{*}){|}_{{{X}}_p}$ is an Ordinary Least Squares (OLS) fit of the training data contained in π_p.

An open question is how many such critical points exist. A starting point is to consider that there are C(N+H, H)≜(N+H)!/N!H! possible partitions of the data. Not every such partition will admit a piecewise-OLS solution which is also continuous, and it is difficult to analytically characterize such solutions.

Using Theorem 4, we can characterize growth of global minima in the overparametrized case. Call a partition Π lonely if each piece π_p ∈ Π contains at most one datapoint (see Figure 4). Then, we can prove the following results:

Lemma 1. For any lonely partition Π, there are infinitely many parameter settings θ_BDSO that induce Π and are global minima with $\tilde{\ell }({\theta }_{\mathrm{\text{BDSO}}})=0$. Furthermore, in the overparametrized regime H ≥ cN for some constant c ≥ 1, the total number of lonely partitions, and thus a lower bound on the total number of global minima of $\tilde{\ell }$ is $\left(\genfrac{}{}{0.0pt}{}{H+1}{N}\right)$ = O(N^cN).

Remark 1. Suppose that the H breakpoints are uniformly spaced and that the N datapoints are uniformly distributed within the region of breakpoints. Then in the overparametrized regime H ≥ dN² for some constant d ≥ 1, the induced partition Π is lonely with high probability 1 − e^−N2/(H+1) = 1 − e^−1/d.

Thus, with only O(N²) hidden units, we can almost guarantee a lonely partition at initialization. This makes optimization easier but is not sufficient to guarantee that learning will converge to a global minimum, as breakpoint dynamics could change the partition to crowded (see section 2.5). Note how simple and transparent the spline-based explanation is for why overparametrization makes optimization easier. For brevity, we include experiments testing our theory of loneliness in Figure 14 in Appendix.

2.4. The Gradient and Hessian of the Loss ℓ(θ_NN)

Previous works has revealed many key insights into the Hessian of the loss surface of neural networks. It has long been empirically known that GD tends to work surprisingly well, and that most or all local minima are actually near-optimal (Nguyen and Hein, 2017), although this is known to depend on the overparametrization ratio and the data distribution (Pennington and Bahri, 2017; Granziol et al., 2019). Flatter local minima have been shown to be connected to better generalization (Sankar et al., 2020), but the areas around critical points may be unusually flat, with the bulk of Hessian eigenvalues near 0, and only a few larger outliers (Ghorbani et al., 2019), meaning that many local minima are actually connected within a single basin of attraction (Sagun et al., 2017). Our theory can shed new light on many of these phenomena:

The gradient of ℓ(θ_NN) can be expressed as

From this, we can derive expressions for the Hessian of the loss ℓ(θ_NN) at any parameter θ_NN.

Theorem 5. Let θ^* be a critical point of ℓ(·) such that ${\beta }_i^{*}({\theta }^{*})\ne x_n$ for all i ∈ [H] and for all n ∈ [N]. Then the Hessian ${\text{H}}_{\ell }({\theta }^{*})$ is the positive semi-definite Gram matrix of the set of 3H + 1 vectors

as shown in eq. (4). Thus, ${\text{H}}_{\ell }({\theta }^{*})$ is positive definite iff the vectors of this set are linearly independent.

Remark 2. The zero eigenvalues of ${\text{H}}_{\ell }({\theta }^{*})$ correspond to

1. α-scaling: For each neuron i, the scaling transformation $(w_i,b_i,v_i)\mapsto ({\alpha }_i{w}_i,{\alpha }_i{b}_i,{\alpha }_i^{-1}{v}_i)$ leaves the approximating function $\widehat{f}(x)$ (and thus the training loss) invariant, thus generating a 1-dimensional hyperbolic manifold of constant loss; thus, the tangent vector of this manifold will be in the null space of ${\text{H}}_{\ell }(\theta )$ for all θ. This manifold is depicted in Figure 1.

2. 0 singularities: For each neuron i, if either of w_i or v_i are 0:

- if v_i = 0, then μ_i = 0, so that changes to the value of β_i (via either b_i or w_i) do not change $\widehat{f}(x)$, thus contributing two 0 eigenvalues

- if w_i = 0, then μ_i = 0 and β_i = ±∞, so any changes to b_i or v_i do not change $\widehat{f}(x)$, thus contributing two 0 eigenvalues

3. Functional changes: Any functional change that does not change the value of $\widehat{f}(x)$ at the training data will leave the loss unchanged. In the sufficiently overparametrized regime, there are many ways to make such changes. Some examples are shown in Figure 5.

The first two types are artifacts of the NN parametrization, caused by invariance to the α-scaling symmetry.

Additionally, we note that these are all also zero eigenvalues for non-critical points, but that there are other zero eigenvalues for non-critical points. While ${\text{H}}_{\ell }({\theta }^{*})$ does not depend on ϵ at the above critical points, at other points, ϵ terms can “cancel out” the corresponding term in eq. (4), yielding additional zero eigenvalues (as well as negative eigenvalues).

From this, we note that for any critical point θ^*, there is a large connected sub-manifold of parameter space of constant loss. Moving along this manifold consists of continuously deforming $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}}^{*})$ such that the output values at datapoints are left invariant or moving along the α-scaling manifolds which leave $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}}^{*})$ invariant. The 0 singularities lie at the intersection of every α-scaling manifold for neuron i, corresponding to the α = 0 and α = ∞ limits.

Consider the μ-only Hessian under the BDSO parametrization, $\frac{d2\ell }{d\mu 2}$, and rewrite our network as $\widehat{\text{y}}=\Phi (\text{x};\beta )\mu$, where $\Phi$ is the N × H feature matrix of activations (or basis functions) (x_n − β_i)_si, which is constant with respect to $\mu$. Then, we have $\frac{d2\ell }{d\mu 2}={\Phi }^{\top }\Phi$. Applying an SVD decomposition $\Phi =US{V}^{\top }$, then ${\Phi }^{\top }\Phi =V{S}^2{V}^{\top }$, such that the eigenvalues of the Hessian will be the singular values of $\Phi$ squared. Thus, we can decompose the vector of delta-slopes $\mu ={\mu }_r+{\mu }_n$ where ${\mu }_r$ is in the range of $\Phi$ while ${\mu }_n$ is in the nullspace. Thus, $\widehat{\text{y}}=\Phi (\text{x};\beta ){\mu }_r+\Phi (\text{x};\beta ){\mu }_n=\Phi (\text{x};\beta ){\mu }_r$, and the gradient with respect to $\mu$ is always constrained to the linear subspace ${\mu }_r$.

Figure 6A shows a 2-dimensional representation of the loss surface, which is quadratic along ${\mu }_r$, and constant along ${\mu }_n$. Figure 6B illustrates the effect of varying $\beta$, which changes $\Phi (\text{x};\beta )$ and hence ${\mu }_r$ and ${\mu }_n$.

2.4.1. The Flatness of the Hessian

An important question is: How does overparametrization impact (the spectrum of) the Hessian? The symmetry-quotiented spline parametrization enables us to lower bound the number of zero eigenvalues of the Hessian [e.g., the “flatness” (Li et al., 2018)] as follows.

Corollary 3. Let θ^* be a critical point of ℓ(·) such that its data partition is lonely, and either at least one neuron is active on all data or there is at least one pair of oppositely-faced neurons in the same data gap, so that $\text{x}\in span({B})$. Then, ℓ(θ^*) = 0 and ${\text{H}}_{\ell }({\theta }^{*})$ has exactly N non-zero eigenvalues, and thus 3H + 1 − N zero eigenvalues.

Intuitively, as overparametrization H/N increases, the number of neurons with shared activation patterns increases, which in turn means many redundant breakpoints between each pair of datapoints, which increases the number of flat directions (zero eigenvalues) of the Hessian. Additionally, each neuron is guaranteed one degree of freedom in the form of the α-scaling symmetry, leading to H intrinsic zero eigenvalues of the Hessian.

2.5. Gradient Flow Dynamics for Spline Parameters

Theorem 6. For a one hidden layer univariate ReLU network trained with gradient descent with respect to the neural network parameters ${\theta }_{\mathrm{\text{NN}}}={\left\{(w_i,b_i,v_i)\right\}}_{i=1}^H$, the gradient flow dynamics of the function space parameters ${\theta }_{\mathrm{\text{BDSO}}}={\left\{({\beta }_i,{\mu }_i)\right\}}_{i=1}^H$ are governed by the following laws:

Here i ∈ [H] indexes all hidden neurons and the initial conditions β_i(0), μ_i(0)∀i ∈ [H] must be specified by the initialization used (see Appendix B.8 for derivation).

2.5.1. Impact of Init Scale α

As mentioned previously, the standard NN parametrization has symmetries such that the function $\widehat{f}$ and the loss is invariant to α-scaling transformations of the NN parameters (section 2.1). However, such scalings do have a large effect on the initial NN gradients; specifically, applying the α-scaling transformation in eqs. (6) and (7), the initial learning rate for β_i scales as $1/{\alpha }_i^2$, while that of μ_i scales approximately as ${\alpha }_i^2$. Changing α at initialization has a large impact on the final learned solution. In Section 2.6 we will show how α determines the kind of (implicit) regularization seen in NN training (Woodworth et al., 2020).

2.5.2. Breakpoint Dynamics

Next, we wish to build some intuition of the conditions that lead to the different knot types and the effect the types have on training dynamics and implicit regularization. First, we rewrite eq. (6) as

and note that the dot product is a cumulative sum of weighted residuals (summing from the furthest datapoint on the active side of neuron i and moving toward β_i). Next, let ${{A}}_i=\left\{n|1_{i,n}=1\right\}$ denote the set of data indices which are active for neuron i, and ${{X}}_i$ be the active half-space of neuron i [i.e., (−∞, β_i) or (β_i, ∞)]. If we view ${\widehat{ϵ}}_i(t)$ as $\widehat{ϵ}({\beta }_i,t)\triangleq {\sum }_{n\in {{A}}_i}{\widehat{ϵ}}_n\delta ({\beta }_i-x_n)$, we can rewrite the dot product as an integral and write $\dot{{\beta }_i}$ as a function ϱ_i(β, t):

Finally, note that ψ₊(β, t) = ψ₋(∞, t)−ψ₋(β, t), i.e., ψ₊(·, ·) is just ψ₋(·, ·) reflected across the (finite) total integral; accordingly, we drop the subscript and refer to ψ(·, ·) from now on.

Then, to identify Repeller and Attractor knots, we are interested in the datapoints β_i = x_n where ${\dot{\beta }}_i$ changes sign, i.e., we are interested in the zero-crossings of ψ(·, t). Very roughly, at earlier stages of learning, large regions of the input space have high residual, so that many neurons get large, similar updates, leading to broad, smooth changes to $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}})$ and hence broad, smooth changes to $\widehat{ϵ}(x,t)$ and ψ(β, t). Applying such a change to ψ(β, t) will have the effect of “centering” ψ(β, t) such that it has an average value of 0 over large regions. Because the changes are broad at first, the steeper regions of ψ(β, t) will retain most of their steepness. Together with the “centering,” this leads to the conclusion that the local extrema of $\frac{\partial \psi (\beta ,t)}{\partial \beta }$ will be near zero-crossings of ψ(β, t).

Due to the broad centering effect, the second term ≈ 0,

Therefore, the local extrema of $\widehat{ϵ}(\beta ,t)$ will be local extrema of $\frac{\partial \psi (\beta ,t)}{\partial \beta }$ and hence near zero-crossings of ψ(β, t).

Note that the above analysis shows that ${\dot{\beta }}_i$ and hence the classification of each datapoint is determined by the same ψ(β, t), plus a scaling factor v_i(t)/w_i(t) and possibly a reflection based on sgn(w_i), i.e., groups of neurons with the same signs of v_i and w_i will yield the same datapoint classification. If datapoint persists in being an Attractor knot for a group of nearby neurons for a long enough time, and is surrounded by Passthrough knots, then those nearby neurons will converge on the datapoint and form a cluster. This leaves the question: what conditions lead to local extrema of $\widehat{ϵ}(\beta ,t)$ that persist long enough for this behavior?

Empirically, we find concentrations of breakpoints forming near regions of high curvature or discontinuities in the training data, leading to a final fit that was close to a linear spline interpolation of the data. In particular, breakpoints migrate toward nearby curvature which is being underfit and whose inflection is the same as that breakpoint's delta-slope. These clusters can remain fixed for some time, before suddenly shifting (all breakpoints in a cluster move together away from the shared attractor x_n) or splitting (two groups of β form two new sub-clusters that move away from x_n in each direction). These new movements can cause “smearing,” when the cluster loses coherence. Sometimes, this “smearing” effect is caused by the need to fit smooth curvature near a cluster: once the residual at the cluster falls below the residual of nearby curvature being underfit by the linear fit provided between clusters, the tendency is for the cluster to spread out; see Figure 7E.

As an illustrative example to explore this question, consider fitting a data set such as that shown in Figure 7. At initialization (Figure 7A), $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}})$ is much flatter than f(x). After a short time of training (Figure 7B), $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}})$ has matched the average value of f(x), but is still much flatter, leading to relatively broad regions where $\widehat{ϵ}(\beta ,t)$ is large and of the same sign, with zero crossings of ψ(β, t) near the first two inflection points. After a little more training (Figure 7C), we see that all three inflection points correspond to zero crossings in ψ(β, t), and that these zero crossings have persisted long enough that clusters of breakpoints have started to form, although the leftmost inflection point is already fit well enough that ψ(β, t) (and hence ${\dot{\beta }}_i$ is quite low in the region, so that those breakpoints are unlikely to make it “all the way” to a cluster before the data is well fit. Later (Figure 7D), the right two inflection points have nearby tight clusters. However, once $\widehat{f}(x;{\theta }_{\mathrm{\text{NN}}})$ reaches f(x) at these clusters, the local maxima of $\widehat{ϵ}(\beta ,t)$ move away from the clusters, and the clusters “spread out,” as shown in the converged fit Figure 7E.

2.5.3. Videos of Gradient Descent/Flow Dynamics

We have developed a BDSO spline parametrization-based visualization. For many of the experiments in this paper, the corresponding videos showing the learning dynamics are available at http://shorturl.at/tFWZ2.

2.6. Implicit Regularization

One of the most useful and perplexing properties of deep neural networks has been that, in contrast to other high capacity function approximators, overparametrizing a neural network does not tend to lead to excessive overfitting (Savarese et al., 2019). Where does this generalization power come from? What is the mechanism? One promising idea is that implicit regularization (IR) plays a role. IR acts as if a extra regularization term had been applied to the optimizer, despite no such regularization term being explicitly added by a designer. Instead, IR is added implicitly, brought about by some combination of architecture, initialization scheme, or optimization. Much recent work (Neyshabur et al., 2015, 2018) has argued that this IR is due to the optimization algorithm itself (i.e., SGD). Our symmetry-quotiented spline perspective makes the role of α-scaling symmetry apparent: IR depends critically on breakpoint and delta-slope learning dynamics and initialization. These results follow from section 2.5, showing that changes to the initialization scale result in dramatic changes to the relative speeds of learning dynamics for different parameters. In particular, in the kernel regime α → ∞, breakpoints move very little whereas delta-slopes move a lot, resulting in diffuse populations of breakpoints that distribute curvature. In stark contrast, in the adaptive regime α → 0 breakpoints move a lot whereas delta-slopes move very little, resulting in tight clusters of breakpoints that concentrate curvature.

2.6.1. Kernel Regime

We first analyze the so-called kernel regime, inspired by Chizat et al. (2019) where it is referred to as “lazy training.” In this section, we omit the overall bias b₀.

Lemma 2. Consider the dynamics of gradient flow on ℓ(·) started from ${\theta }_{\mathrm{\text{NN}},\alpha }\triangleq (\alpha {\text{w}}_0,\alpha {\text{b}}_0,{\text{v}}_0=0)$, where w_i ≠ 0∀i ∈ [H]. In the limit α → ∞, $\beta (t)$ does not change, i.e., each breakpoint location and orientation is fixed. In this case, the θ_NN model reduces to a (kernel) linear regression:

where $\mu \in {ℝ}^H$ are the regression weights and Φ(x; β) ∈ ℝ^N×H are the nonlinear features ϕ_ni≜(_xn−βi)si.

Using this, we can then adapt a well known result about linear regression (see e.g., Zhang et al., 2016):

Theorem 7. Let ${\mu }^{*}$ be the converged $\mu$ parameter after gradient flow on the BDSO model eq. (8) starting from some ${\mu }_0$, with $\beta$ held constant. Furthermore, suppose that the model perfectly interpolates the training data $\tilde{\ell }({\theta }_{\mathrm{\text{BDSO}}})=0$. Then,

Thus, the case where breakpoint locations and orientations are fixed, and ${\mu }_0=0$ reduces to ℓ₂-regularized linear regression on the delta-slope parameters $\mu$. A geometric representation of this process is shown in Figure 8. Recalling that ${\Vert \mu \Vert }_2^2$ is the roughness of $\widehat{f}(x;{\theta }_{\mathrm{\text{BDSO}}})$, this result demonstrates the importance of the initialization, and in particular the initial roughness (Theorem 2): with a high-roughness initialization minimizing ${\Vert \mu -{\mu }_0\Vert }_2^2$ will yield a high-roughness final fit.

In the overparametrized regime with an initialization that yields a lonely partition, the model will reach $\tilde{\ell }({\theta }_{\mathrm{\text{BDSO}}})=0$. If we consider the infinite-width limit, we get the following result, which is a specialization of Theorem 5 of Williams et al. (2019):

Corollary 4. Consider the setting of Theorem 7, with the additional assumptions that ${\mu }_0=0$ and the breakpoints are uniformly spaced, and let H → ∞. Then the learned function ${\widehat{f}}_{\infty }(x;{\mu }^{*},{\beta }^{*})$ is the global minimizer of

As such, ${\widehat{f}}_{\infty }(x;{\mu }^{*},{\beta }^{*})$ is a natural cubic smoothing spline with N degrees of freedom (Ahlberg et al., 1967).

Theorem 5 of Williams et al. (2019), is a generalization of Corollary 4 to the non-uniform case which replaces Equation (9) with

where p_β(β) is the initial density of breakpoints induced by the specific initialization used.² Deriving this result using the BDSO framework allows us to see that the initialization has the impact of weighting the curvature of certain locations more than others.

2.6.2. Kernel Regime Dynamics as PDE

Assume a flat (μ_i(0) = 0 ∀i ∈ [H]) kernel regime initialization. Specializing our spline parameterization dynamics results (Theorem 6), we have

Note that the terms r_{1, i}(t) and r_{x, i}(t) only depend on the index i through the mask $1_i$; in other words, they are the same for all breakpoints with the same activation pattern. This pattern is entirely determined by the orientation s_i and the data interval (x_n, x_n+1) that β_i falls into.

Thus, the possible activation patterns can be indexed by data point index n and orientation s. Letting ς_{n, s} denote the set of breakpoints with the activation pattern corresponding to (n, s), we have

Let ${\text{r}}_{1,n,s}(t)$ and ${\text{r}}_{x,n,s}(t)$ be vectors containing |ς_{n, s}| copies of r_{1, n, s}(t) (resp. r_{x, n, s}(t)), and let ${\text{r}}_1(t)$ and ${\text{r}}_x(t)$ be the concatenation of these over all n and s. This allows us to write

where $\mu$ and $\beta$ are the vectors of μ_i (resp. β_i) values for all i. Then, $\mu (t)$ can be viewed as a function $\mu :\left[H\right]\times \left\{+,-\right\}\times ℝ\to ℝ$, which is isomorphic to the function μ:ℝ × {+, −} × ℝ → ℝ given by $\mu (x,s,t)\triangleq \sum _{i=1}^H{\mu }_i(t)\delta (x-{\beta }_i)$, yielding the PDE

where r₁(x, s, t) and r_x(x, s, t) are piecewise constant functions of x with discontinuities at datapoints. Because r₁(x, s, t) is piecewise constant in x for all t, ${\int }_0^t{r}_1(x,s,\tau )\mathrm{\text{d}}\tau$ will also be piecewise constant, and likewise for r_x(x, s, t). Thus, μ(x, s, t) will be (discontinuous) piecewise linear in x for all t. Stepping back, we see that a finite-width ReLU network trained by gradient descent is a discretization (in both space and time) of this underlying continuous PDE.

Thus, $\widehat{f}$”_∞(x; t) = μ(x, +, t)+μ(x, −, t), and as μ(x, s, t) is (discontinuous) piecewise linear, this implies that ${\widehat{f}}_{\infty }(x;t)$ will be continuous piecewise cubic in x for all t, consistent with Corollary 4.

2.6.3. Adaptive Regime

Previous work (Arora et al., 2019a) has shown that the kernel regression regime is insufficient to explain the generalization achieved by modern DL. Instead, the non-convexity of the optimization (i.e., the adaptive regime, α → 0, called the “rich regime” in Woodworth et al., 2020) must be responsible. The spline parametrization clearly separates the two, with all adaptive regime effects due to breakpoint dynamics β_i(t): as α → 0, breakpoints become more mobile, and the breakpoint dynamics described in section 2.5 becomes more and more dominant. As shown in section 2.3, high breakpoint mobility can lead to clusters of breakpoints forming at certain specific data points. For low enough α, this clustering of breakpoints yielding a fit that is more linear spline than cubic spline; intermediate values of α interpolates these two extremes, giving a fit that has higher curvature in some areas than others.

Intriguingly, this suggests an explanation for why the kernel regime is insufficient: the adaptive regime enables breakpoint mobility, allowing the NN to adjust the initial breakpoint distribution (and thus, basis of activation patterns). This suggests that data-dependent breakpoint initializations may be quite useful (see Table 1 for a simple experiment in this vein).

Spline theory traditionally places a knot at each datapoint for smoothing splines (James et al., 2013). However, in circumstances where this is not feasible, either computationally or due to using a different spline, a variety of methods exist (MLE, at set quantiles of the regression variable, greedy algorithms, and more) (Park, 2001; Ruppert, 2002; Walker et al., 2005). More recent work has started to develop adaptive splines, where knots are placed to minimize roughness, length, or a similar metric, and knot locations are found via Monte Carlo methods or evolutionary algorithms (Miyata and Shen, 2005; Goepp et al., 2018). This adaptivity allows the spline to better model the function by moving knots to areas of higher curvature – remarkably similar to the behavior of a low-α NN.

In general, an individual breakpoint does not have a large impact on the overall function (O(1/H) on average); it is only through the cooperation of many breakpoints that large changes to the function are made. Another way of formulating this distinction is that the function (and hence the loss) depend on the breakpoints only through the empirical joint density ${\widehat{p}}_H(\beta ,\mu ,s)$. Similarly, training dynamics depend on the θ_NN joint density ${\widehat{p}}_H(w,v,b)$. This formulation is explored in Mei et al. (2018), where they derive a dynamics equation for the joint density (in θ_NN). Adapting this work to explore the dynamics of the θ_BDSO joint density is ongoing work.

2.6.4. Comparison With Concurrent Work

Independent of and concurrent with previous versions (Sahs et al., 2020a,b) of this work, Williams et al. (2019) has implicit regularization results in the kernel and adaptive regimes which parallel our results in this section rather closely. Despite the similarities, we take a significantly different approach. Comparing our results to those in Williams et al. (2019), the key differences are: (1) our BDSO parametrization has a clear geometric/functional interpretation whereas Williams et al.'s canonical parameters are opaque. (2) BDSO generalizes straightforwardly to high dimensions in a conceptually clean manner: oriented breakpoints become oriented breakplanes (see section 2.7); it is not clear what the multivariate analog of the canonical parameters would be. (3) Lemma 2 of Williams et al. (2019) and the surrounding text discuss the existence of what we call Attractor/Repulsor knot datapoints, but only gives an algebraic expression for their formation in terms of residuals and are unable to distinguish between the types, whereas we relate the persistence of such Attractors to the curvature of the ground truth function, particularly when breakpoints cross datapoints. (4) Our proof techniques are quite different from theirs, with different intuitions and intermediate results, and are also conceptually simpler. (5) Theorem 5 of Williams et al. (2019) is slightly more general form of our Corollary 4, extending to the case of non-uniform breakpoints. (6) Finally, we have many novel results regarding initialization, loss surface properties, Hessian and dynamics (everything outside of section 2.6). All in all, we believe our results are quite complementary to those of Williams et al., particularly as we extend our results to novel areas using our Hessian and loss surface analysis.

2.7. Extending to Multivariate Inputs

Throughout this work we chose to focus on the univariate case in order to build intuition and enable simpler theoretical results. However, in practice most networks have multivariate inputs; fortunately, our BDSO parametrization can be easily extended to multivariate inputs. For D-dimensional inputs, write

where the input weights are represented as D-dimensional vectors ${\text{w}}_i$. Using the reparametrization ${\eta }_i\triangleq v_i{\Vert {\text{w}}_i\Vert }_2$, ${\xi }_i\triangleq \frac{{\text{w}}_i}{{\Vert {\text{w}}_i\Vert }_2}$, ${\gamma }_i\triangleq \frac{-b_i}{{\Vert {\text{w}}_i\Vert }_2}$, the representation becomes

where η_i is a “delta-slope” parameter,³ and $({\xi }_i,{\gamma }_i)$ parametrize the orientation and signed distance from the origin of a D − 1-dimensional “oriented breakplane” (generalizing the 0-dimensional left-or-right oriented breakpoint represented by (s_i, β_i) in the 1-dimensional case). Generalizing our results to this parametrization is ongoing work.

3. Experiments

3.1. Suboptimality of Gradient Descent

Focusing on univariate networks allows us to directly compare against existing (near-) optimal algorithms for fitting Piecewise Linear (PWL) functions, including the Dynamic Programming algorithm (DP, Bai and Perron, 1998), and a very fast greedy approximation known as Greedy Merge (GM, Acharya et al., 2016) (Figure 9). Given a fixed budget of pieces (∝ to number of parameters e.g., network width), how well do these algorithms compare to SGD variants in fitting a quadratic (high curvature) function? DP and GM both quickly reduce training error to near 0 with order 100 pieces, with GM requiring slightly more pieces for similar performance. All the GD variants require far more pieces to reduce error below any target threshold, and may not even monotonically decrease with number of pieces. These results show how inefficient GD is w.r.t parameters, requiring orders of magnitude more for similar performance compared with PWL fitting algorithms.

3.2. Effect of Initial Breakpoint Distribution

We first ask whether the standard initializations will experience difficulty fitting functions that have significant curvature away from the origin (e.g., learning the energy function of a protein molecule). We train ReLU networks to fit a periodic function (sin(x)), which has high curvature both at and far from the origin. We find that the standard initializations do quite poorly away from the origin (Table 1, first row), consistent with our theory that breakpoints are essential for modeling curvature. Probing further, we observe empirically that breakpoints cannot migrate very far from their initial location, even if there are plenty of breakpoints overall, leading to highly suboptimal fits. In order to prove that it is indeed the breakpoint density that is causally responsible, we attempt to rescue the poor fitting by using a simple data-dependent initialization that samples breakpoints uniformly over the training data range [x_min, x_max], achieved by exploiting eq. (1). We train shallow ReLU networks on training data sampled from a sine and a quadratic function, two extremes on the spectrum of curvature. The data shows that uniform breakpoint density (Table 1, second row) rescues bad fits in cases with significant curvature far from the origin, with less effect on other cases, confirming the theory. We note that this could be a potentially useful data-dependent initialization strategy, one that can scale to high dimensions, but we leave this for future work.

3.3. Generalization: Implicit Regularization Emerges From Flat Init and Curvature-Based Parametrization

Our theory predicts that IR depends critically upon flatness of the initialization (Theorem 7 and corollary 4). Here, we test this theory for the case of shallow and deep univariate ReLU nets. We compare training with the standard flat initialization to a “spiky” initialization, and find that both fit the training data near perfectly, but that the “spiky” initialization has much worse generalization error (Figure 10 Top and Table 3 in Appendix).

It appears that generalization performance is not automatically guaranteed by GD, but is instead due in part to the flat initializations which are then preserved by GD. Our theoretical explanation is simple: integrating the dynamics in Equation (7) yields μ_i(t) = μ_i(0)+⋯  and so the initialization's impact remains.

3.4. Impact of Width and Init Variance

Next, we examine how smoothness (roughness) depends on the width H, focusing on settings with large gaps in the training data. We use two discontinuous target functions (shown in Figure 11), leading to a gap in the data, and test how increasing H (with α unchanged) affects the smoothness of $\widehat{f}$. We test this on a “smooth” data gap that is easily fit, as well as a “sharp” gap, where the fit will require a sharper turn. We trained shallow ReLU networks with varying width H and initial weight variance σ_w until convergence, and measured the total roughness of resulting CPWL approximation in the data gaps.

Figure 10 Bottom shows that roughness in the data gaps decreases with width and increases with initial weight variance, confirming our theory. A higher weight variance, and thus rougher initialization, acts in a similar fashion to the “spiky” initialization, and leads to increased roughness at convergence. In contrast, higher width distributes the curvature “work” over more units, leading to lower overall roughness.

3.5. Impact of Init Scale α

Finally, we explore how changing α impacts IR. Empirically, as α increases from 0 to ∞ we see three qualitative regimes: (i) an underfitting linear spline, (ii) an interpolating linear spline, and (iii) a roughness-minimizing natural cubic spline. This is quantified in Table 2, where we compare the NN fit to a linear interpolation and a natural cubic spline fit, for varying α. We first test in the special case that the initial function approximation is perfectly flat; we find excellent agreement with the linear interpolation and cubic spline fits for α = 3, 100 (Uniform initialization) and α = 10, 100 (He initialization). The impact of α on IR can be more easily visualized in our Supplementary Videos. In order to gauge the impact of the initialization breakpoint distribution, we also test with a standard He initialization (Cauchy distributed breakpoints, Corollary 1). In this case, we find that generalization error is uniformly higher for all α. More strikingly, the α regime (ii) above disappears, as a result of breakpoints being more concentrated around the origin and the initialization roughness being significantly nonzero. This supports the idea that the initial parameter settings, in particular the breakpoint distribution, has a critical impact on the final fit and its IR.

Taken together, our experiments strongly support that a smooth, flat initialization and overparametrization are both responsible for the phenomenon and strength of IR, while the initialization weight scale α critically determines the type of IR.

4. Discussion

We show that removing the α-scaling symmetry and examining neural networks in spline space enabled us to glean new theoretical and practical insights. The spline view highlights the importance of initialization: a smooth initial approximation is required for a smooth final solution. Fortunately, existing initializations used in deep learning practice approximate this property. In spline space, we also achieve a surprisingly simple and transparent view of the loss surface, its critical points, its Hessian, and the phenomenon of overparametrization. It clarifies how increasing width relative to data size leads with high probability to lonely data partitions, which in turn are more likely to reach global minima. The spline view also allows us to explain the phenomenon of implicit regularization, and how it arises due to overparametrization and the initialization scale α.

4.1. Related Work

Our approach is related to previous work (Frankle and Carbin, 2018; Arora et al., 2019b; Savarese et al., 2019) in that we wish to characterize parameterization and generalization. We differ from these other works by focusing on small width networks, rather than massively overparametrized or infinite width networks, and by using a spline parameterization to study properties such as smoothness of the approximated function. Previous work (Advani and Saxe, 2017) has hinted at the importance of low norm initializations, but the spline perspective allows us to prove implicit regularization properties in shallow networks. Finally, Williams et al. (2019) is closely related and is discussed in detail at the end of section 2.6.

4.2. Explanation for Correlation Between Flatness of Minima and Generalization Error

A key unexplained empirical observation has been that flatter local minima tend to generalize better (Li et al., 2018; Wei and Schwab, 2019). Our results above provide an explanation: as overparametrization ${O}=H/N$ increases, the flatness of the local minima (as measured by the number of zero eigenvalues) increases (Corollary 3) and the smoothness of the implicitly regularized function (as measured by inverse roughness $\rho ({\widehat{f}}_H)\ge \rho ({\widehat{f}}_{\infty })$) also increases. As previously shown (Wu et al., 2017), flatter and simpler local minima imply better generalization. Our work provides a parsimonious explanation for this: as we increase overparametrization, partitions become increasingly lonely, allowing for an increased redundancy in number of ways to exactly fit the training data (thus increasing the number of zero eigenvalues), while the inductive bias of gradient descent spreads the “work” (e.g., curvature changes due to delta-slopes) among many units, ensuring that each unit has a lesser effect and making the loss surface increasingly flat.

4.3. Future Work

Looking forward, there are still many questions to answer from the spline perspective: How does depth affect the expressive power, learnability, and IR? What kinds of regularization are induced in the adaptive regime and how do modern deep nets take advantage of them? How can data-dependent initializations of the breakpoints help rescue/improve the performance of GD? Can we design new global learning algorithms inspired based on breakpoint (re)allocation? We believe the BDSO perspective can help answer these questions.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JS led effort on theory section, with help from RP and AP. RP led effort on experiments, with help from all other authors. JS and RP were responsible for manuscript. All authors contributed to the article and approved the submitted version.

Funding

RP and AP were supported by Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/Interior Business Center (DoI/IBC) contract number D16PC00003. JS and AP were supported by NIH grant no. 1UF1NS111692-01. RP and AP were supported by funding from NSF NeuroNex grant no. 1707400.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frai.2022.889981/full#supplementary-material

Footnotes

References