Title: WIDE NEURAL NETWORKS ARE INTERPOLATING KERNEL METHODS: IMPACT OF INITIALIZATION ON GENERALIZATION

Abstract: The recently developed link between strongly overparametrized neural networks (NNs) and kernel methods has opened a new way to understand puzzling features of NNs, such as their convergence and generalization behaviors. In this paper, we make the bias of initialization on strongly overparametrized NNs under gradient descent explicit. We prove that fully-connected wide ReLU-NNs trained with squared loss are essentially a sum of two parts: The first is the minimum complexity solution of an interpolating kernel method, while the second contributes to the test error only and depends heavily on the initialization. This decomposition has two consequences: (a) the second part becomes negligible in the regime of small initialization variance, which allows us to transfer generalization bounds from minimum complexity interpolating kernel methods to NNs; (b) in the opposite regime, the test error of wide NNs increases significantly with the initialization variance, while still interpolating the training data perfectly. Our work shows that -contrary to common belief -the initialization scheme has a strong effect on generalization performance, providing a novel criterion to identify good initialization strategies.

1. INTRODUCTION:

Neural networks (NNs) have celebrated many successes over the past decade and achieved state of the art performance across various domains and tasks. From a theoretical standpoint, however, many aspects of neural networks are not well understood, as for example illustrated by Zhang et al. (2016) . Neural networks seem to contradict classical learning theory as, in many scenarios, they are able to fit random labels perfectly while still generalizing well when trained on the true labels. In addition, overparametrized neural networks frequently exhibit even improved test performance when the number of parameters is increased further (Belkin et al., 2019) . NN models thus often seem to avoid overfitting.

Very recently there have been advances in understanding the training and evaluation behavior of neural networks in the infinitely wide limit (Jacot et al., 2018; Hayou et al., 2019) and also in the strongly overparametrized regime (Du et al., 2018b; Li & Liang, 2018; Allen-Zhu et al., 2018a) , which is close to the infinite limit NN. Both lines of work express the behavior of the neural network in terms of the so-called neural tangent kernel (NTK). In particular, Du et al. (2018b) showed in this way that, under mild conditions, strongly over-parametrized neural networks converge to a global minimum of zero training error.

In another line of research, Belkin et al. (2018b) suggested a new picture of learning in the overparametrized regime, introducing an "interpolating kernel method" which successfully learns in this regime. This method selects the least complex function that interpolates all data points perfectly, as opposed to traditional methods which balance the function's complexity against its goodness-of-fit. It is suggested that this picture of overparametrized learning could help understand neural networks.

If overparametrized NNs are indeed linked to such kernel methods, then zero (or small) training error alone might not tell us much about the test performance, as the measure of function complexity could be determined by the NN architecture, initialization, and the optimization method. A bad NN design may lead to an unfavorable complexity measure that could result in a large generalization gap.

The question therefore remains in which way exactly are overparametrized neural networks linked to interpolating kernel methods? And how can this connection help us better understand neural networks and their puzzling features?

1.1 OUR CONTRIBUTIONS In this paper we answer the first question by making the link to interpolating kernel methods explicit for strongly overparametrized NNs. After that, we exhibit two implications for NNs which this connection allows us to draw.

First, for strongly overparametrized ReLU-NNs we make the bias explicit that is imposed by the NN initialization and by training with discrete gradient descent steps on the squared loss. We achieve this by decomposing the learned NN function as well as the test error into two terms: The first term corresponds to a minimum complexity interpolating kernel method, whereas the second term is proportional to the initialized weights and gives a non-zero contribution only if the test input is not part of the training set.

Second, we are able to bound the difference between a NN with standard initialization of small variance and the solution of a minimum complexity interpolating kernel method based on the NTK. This bound provides a new way to transfer existing test error bounds from interpolating kernel methods to NNs, and vice versa.

Our third contribution shows that without putting any constraints on the NN initialization, low training error does not imply low test error. This is because the second term in the above test error decomposition grows as a power of the initialization variance, as we show, whereas the first term is independent of this variance. Thus, the test error grows dramatically with increasing initialization variance, despite (close to) zero training error. Our findings provide an additional theoretical way to analyze initialization schemes, besides the existing heuristics and intention to prevent vanishing and exploding gradients (He et al., 2015; Glorot & Bengio, 2010) .

Our theoretical bounds and insights are nicely corroborated in several experimental settings.

1.2. NOTATION:

We consider a supervised learning setting with data

D = {(x i , y i )} N i=1 ⊂ B d 1 (0) × R consisting of N training points, i.e. the inputs come from the bounded d-dimensional ball B d 1 (0) = {x ∈ R d : ||x|| 2 ≤ 1} 1 .

We often view the training inputs X = (x 1 , . . . , x N ) T ∈ R N ×d and corresponding labels Y = (y 1 , . . . , y N ) T ∈ R N in matrix form. For any function g with domain R d , we define g(X) as the row-wise application of g. Similarly, for any set

D test ⊂ R d × R of N test test points we define X test ∈ R Ntest×d and Y test ∈ R Ntest . We call a dataset D consistent if x i = x j implies y i = y j .

We denote the output of a fully-connected NN with L layers and parameters θ by

EQUATION

where a(h) := max{h, 0} is the ReLU-activation function applied componentwise, m l denotes the number of neurons in layer l, and W l ∈ R m l ×m l-1 and b l ∈ R m l are the weights and biases.

For simplicity we choose the number of units in each hidden layer to be the same2 m l = m for l = 1, . . . , L -1 and m L = 1. θ l = [W l , b l ] denotes the vector of weights and biases in layer l, and θ = [θ 1 , . . . , θ L ] the vector of all parameters, whose total number is M = L l=1 (m l-1 + 1)m l . We keep the "initialization variance" σ 2 (equivalently, σ > 0) as an explicit scalar parameter, that can be varied.

Our parametrization (1) together with a standard normal initialization W l i,j , b l i ∼ N (0, 1) is the so-called NTK-parametrization (e.g. Du et al. (2018b) ). In a forward pass this is equivalent to the standard parametrization (setting σ/ √ m l → 1 in Eq. ( 1)) with the weights and biases initialized at variances σ 2 /m l and σ 2 , respectively (He et al., 2015) . During gradient training the two parametrizations are equivalent up to a parameter-dependent scaling factor of the learning rate (Lee et al., 2019) .

We consider the squared error loss (ŷ, y) = |ŷ -y| 2 /2 and train by minimizing the empirical risk

L D (θ) = (x,y)∈D (f NN θ (x), y) = f NN θ (X) -Y 2 2 /2 (2)

using gradient descent (GD). That is, starting from the initialization θ = θ 0 the weights are updated according to the discrete iteration θ t+1 = θ t -η∇ θ L(θ t ) for t = 0, 1, . . ., where η is a learning rate.

We connect NNs to interpolating kernel methods w.r.t. two specific kernel functions that are defined for any initialization θ = θ 0 of a NN (1). The first kernel Ψ : R d × R d → R is associated with training only the last layer θ L of a NN (Lee et al., 2019) ,

EQUATION

is the corresponding feature map. The second is the so-called neural tangent kernel (NTK) Φ : R d × R d → R associated with training all NN parameters (Jacot et al., 2018; Lee et al., 2019) ,

EQUATION

again represents the feature map. The "minimum complexity interpolating kernel method" fits a consistent training set D perfectly, using a function of minimum kernel norm w.r.t. a kernel K:

EQUATION

where H K is the reproducing kernel Hilbert space (RKHS) associated with K (Belkin et al. (2018b) ; see also App. B). Additionally, we denote the solution of a fully converged fully-connected ReLU-NN trained with gradient descent, squared loss and random initialization by f NN . We denote the test losses for these two predictors by L meth test = f meth (X test ) -Y test 2 /(N test ) for meth = NN, int. (2019b) developed data dependent generalization bounds based on the fact that the weights do not move far from their initialization, as shown in the above works. Furthermore, Arora et al. (2019a) experimentally evaluate the generalization behaviour of the infinite width solution of NNs in the NTK-limit.

2. RELATED WORK:

Most of these contributions focus on showing that the training loss converges to a global minimum or that the training dynamics is close to a trajectory that converges to zero training loss. We build upon these works, but in contrast our focus is the connection between the test error of a trained NN and the error of an interpolating kernel method. We achieve this by resolving the implicit bias that the initialization and gradient descent (GD) has on the test loss. The initialization bias on the training loss is minor compared to the test loss. The most related work to ours is by Zhang et al. (2019) . They investigate a similar issue as we do but use a different method to weaken the effects of non-zero initialization and also do not quantify the impact of initialization on generalization but give evidence for such an effect. The recent contribution Lee et al. (2019) , building upon Lee et al. (2017 ), Matthews et al. (2018 ) and Neal (1996) , explicitly solves the training dynamics and shows that the solution in the trained infinite model seems to be similar to a Gaussian Process, even though its covariance is not the Bayesian posterior of a simple prior. The solution of their continuous training time analysis is similar to the linearized model in the proof of our Thm. 2. We perform a discrete-time analysis, as in actual NN training, instead of the continuous-time (ODE-like) analysis of other works. Woodworth et al. (2019) uses a continuous time analysis as well to focus on the effect of noise in the gradients. Borovykh (2019) investigates the effect of the scaling between the "deep" and the "kernel" regime.

The second line of related work introduces interpolating kernel methods and the idea of an overparametrized learning regime which does not suffer from overfitting even though the number of parameters is be much larger than the number of training samples. This line of work was mainly introduced by Belkin et al. (2018b) and the follow-up works Ma et al. (2017) ; Belkin et al. (2018c; a) . The main idea is, instead of regularizing the function while trying to fit the data points well, to regularize (the complexity of) a function that perfectly interpolates all data points. This newly identified overparametrized learning regime challenges the conventional thinking about the biasvariance-tradeoff (Vapnik, 2013) from the underparametrized regime. Belkin et al. (2019) introduces a bouble-dip picture to describe the transition between the two regimes. Liang & Rakhlin (2018) have derived generalization bounds for these interpolating kernel methods using novel techniques.

Using our explicit link to NNs, this will provide an alternative to NN generalization bounds like margin or PAC-Bayes-based bounds (e.g., Bartlett et al. (2017) ; Neyshabur et al. ( 2017)). While a connection between overparametrized NNs and interpolating kernel methods has been suggested in the papers above, the link has not been made explicit in a rigorous way, to the best of our knowledge.

The implicit bias of (stochastic) gradient descent ((S)GD) and random initialization on the trained NN has been the topic of many previous works, e.g. 

3. LINKING WIDE NNS AND INTERPOLATING KERNEL METHODS:

We now develop the link between strongly overparametrized neural networks and minimum complexity interpolating kernel methods mathematically. We first illustrate the main idea with the simpler case where only the last NN layer is trained, before extending the result to general NN training.

3.1. PARADIGMATIC LINEAR CASE: TRAINING ONLY THE LAST NN LAYER:

As the NN function f NN θ (x) (1) is linear in the parameters θ L of the last layer, the associated features ψ(x) = ∇ θ L f NN θ (x) θ0 and corresponding kernel Ψ from (3) do not depend on the last-layer initialization θ L 0 (but depend on the other layers θ 1:L-1 0

), and f NN θ (x) = θ L , ψ(x) holds exactly for all θ and x. Starting NN gradient descent at θ L 0 = 0, we get the following simple link:

Theorem 1 Let f NN be the fully converged solution of a fully-connected strongly overparametrized ReLU-NN (1) with L layers, where only the last layer θ L = [W L , b L ] has been trained, using gradient descent under squared loss on a consistent dataset D, after it had been initialized at θ L 0 = 0, and all other layers have been kept at their random initialization ∼ N (0, 1). Then it holds that

EQUATION

) where f int Ψ is the minimum complexity interpolating kernel solution (5) of D w.r.t. kernel Ψ (3).

Proof. Abbreviating the last-layer parameters by ϑ = θ L and using the aforementioned linearity f θ (x) = ϑ, ψ(x) in ϑ, the gradient descent update rule reads:

EQUATION

This iteration can be solved explicitly using induction and the binomial theorem (Shah et al., 2018) :

EQUATION

where (ψ(X)ψ(X) T ) -1 denotes the pseudoinverse.

If the learning rate satisfies η < 2/λ max (ψ(X)ψ(X) T ) such that gradient descent converges, we get 1 -ηψ(X)ψ(X) T t → 0 in the limit t → ∞ as the spectral radius is smaller than one and the term in square brackets can be simplified to 1. To be precise this holds since D is consistent and we are in the strongly overparametrized regime, so that Y (and ψ(X)ϑ 0 ) lies in the range of ψ(X)ψ(X) T almost surely (over the random initialization θ 1:L-1 0 , on which ψ depends), see also Du et al. (2018b) .

Setting ϑ 0 = 0 as presupposed, we obtain the following predictor at full convergence (t = ∞):

f NN (x * ) = ϑ ∞ , ψ(x * ) = ψ(x * )ψ(X) T (ψ(X)ψ(X) T ) -1 Y.

(9) The same value ϑ ∞ can be obtained as the solution to the following least squares interpolation:

EQUATION

)

This can in turn be written as minimizing the RKHS norm of an interpolation function w.r.t. the kernel Ψ from Eq. ( 3) (for details on this see App. A and B):

EQUATION

The claim has thus been proven by deriving the explicit solution ( 9) and showing it to equal f int Ψ .

3.2. GENERAL CASE: TRAINING THE FULL NEURAL NETWORK:

In the case of training the full NN, we want to follow similar steps as in Sect. 3.1. For this, we first approximate f NN θ (x) affine-linearly by f lin θ around its initialization θ = θ 0 , using recent results on wide NNs (Du et al., 2018b) . We further drop the requirement that any components of θ 0 must vanish; note, e.g., that the suggestive requirement θ 0 = 0 (analogous to Sect. 3.1) would, for L ≥ 1, lead to vanishing features φ(x) = ∇ θ f NN θ (x) θ0 and a degenerate kernel Φ = 0. Nonzero θ 0 will however result in an initialization bias in f lin θ as we will see (Eq. ( 14)).

Theorem 2 Let f NN be the fully converged solution of a strongly overparametrized ReLU-NN (1) with L layers, each with m neurons, where all layers θ have been trained using gradient descent under squared loss on a consistent dataset D after random initialization θ 0 ∼ N (0, 1). Then there exists a N ∈ N such that for all NN with m ≥ N and ∀x * ∈ B d 1 (0) it holds with probability at least 1 -2δ over the initialization that

|f NN (x * ) -f int Φ (x * )| ≤ σ L δ 1/2 + O 1 m 1/2 (12)

where f int Φ is the minimum complexity interpolating kernel solution (5) of D w.r.t. kernel Φ (4).

For the proof of Thm. 2 we describe here especially those steps which are new compared to Thm. 1. The first property we need is the approximate linearity of the NN in its parameters θ during the whole gradient training. This is the reason why we restrict to strongly overparametrized fully-connected NNs. For such networks, a combination of the works (Jacot et al. (2018) 

(x) := f NN θ0 (x) + (∇ θ f NN θ (x) θ0 )(θ -θ 0 )

. Further, for t ≥ 0, let θt and θ t be the parameters obtained by gradient descent training starting from θ 0 with sufficiently small step size η of the full NN and its linearisation respectively. Then, there exists some N ∈ N such that for all m ≥ N it holds with probability at least 1-δ over the initialization that:

sup t |f NN θt (x)-f lin θt (x)| 2 ≤ O(1/m) for all x ∈ B d 1 (0).

Similar linearisation hold for strongly overparametrized RNNs and CNNs as well, not only for fully-connected NNs. The theorem proves that for every sufficiently wide NN the linearised solution is close to the trained full NN. For a more detailed analysis of the dependence of N on the width m, variance σ and δ we refer to Lee et al. (2019); Du et al. (2018b) . Using this result, we continue our analysis with f lin which we can rewrite as

f lin θ (x) = f NN θ0 (x) + θ -θ 0 , φ(x) = f NN θ0 (x) -θ 0 , φ(x) + θ, φ(x) were φ(x) = ∇ θ f NN θ (x)| θ0

is the feature vector defined in Eq. ( 4). Similar to the proof of Thm. 1, gradient descent on f lin θ with the loss function f lin θ (X) -Y 2 /2 and sufficiently small step size η < 2/λ max (ψ(X)ψ(X) T ) leads to the final (t → ∞) parameter

θ lin ∞ = θ 0 + φ(X) T φ(X)φ(X) T -1 (Y -f NN θ0 (X))

, resulting in the function

EQUATION

To connect Eq. ( 13), and thereby f NN , to the solution f int Φ of the interpolating model, two steps remain to be done. First we show that the two summands in Eq. ( 13) which are not proportional to the labels Y can be simplified into a single form for fully-connected ReLU-NNs; this can be done by showing f θ0 (x) = 1 L θ 0 , Φ(x) (see Lemma 2 for more details). Secondly, we can bound this single expression from above in terms of the initialization variance (Lemma 3; for proofs see Appendix. The term proportional to Y is then identified with f int .

Lemma 2 For a fully-connected ReLU-NN f NN θ (Eq. 1) with L layers, the trained linearized version f lin from Eq. ( 13) can be written as follows:

f lin (x) = φ(x)φ(X) T (φ(X)φ(X) T ) -1 Y + 1 L φ(x) 1 -φ(X) T (φ(X)φ(X) T ) -1 φ(X) θ 0 . (14)

Thus, the trained

f NN (x) = f lin (x) + O(1/ √ m)

decomposes essentially into the two terms of Eq. ( 14) (for x ∈ B d 1 (0)).

As |f NN (x) -f lin (x)| ≤ O(1/ √ m)

, the expression ( 14) holds for the trained f NN in the wide limit m → ∞ as well. This technical result enables all of our results in Thms. 2, 3, and 4.

Eq. ( 14) can be understood intuitively (in the case of a consistent D): The first term makes the model interpolate the training data D perfectly (f lin (X) = Y ), independently of the initialization θ 0 , whereas the second term vanishes on the training data x = X and thus only contributes to test inputs x = x * . The second term furthermore depends on the initialization θ 0 , more precisely on the component of θ 0 that is orthogonal to the feature manifold spanned by the data D, i.e. orthogonal to the range of φ(X) T . Note that 1 -φ(X) T (φ(X)φ(X) T ) -1 φ(X) is the projector onto the kernel of φ(X), about which the dataset D does not give any information.

Lemma 3 For a fully-connected ReLU-NN f NN θ , it holds for all x ∈ B d 1 (0) with probability at least 1 -δ over the initialization θ 0 that 1

L |φ(x) 1 -φ(X) T (φ(X)φ(X) T ) -1 φ(X) θ 0 | ≤ σ L √ δ .

As in the proof of Thm. 1 (see Eqs. ( 9)-( 11)) the first part of the decomposition in Eq. ( 14) is f int Φ . By Lemma 3, we see that for small σ 2 the distance between f lin and f int Φ becomes small with high probability, and f lin and f NN are close too because of m N . The triangle inequality

|f NN (x * ) -f int Φ (x * )| ≤ |f NN (x * ) -f lin /P hi (x * )| + |f lin (x * ) -f int (x * )|

gives the desired result. Thus, strongly overparametrized neural networks initialized with small variance and trained under gradient descent correspond to solutions of a minimum complexity interpolating kernel method.

4. IMPLICATIONS FOR NEURAL NETWORKS:

Now that we have established the link between NNs and interpolating kernel methods, we show two implications this has for the test loss of neural networks. First, we derive an upper bound on the NN test loss, which -in the regime of small NN initialization variance -allows us to transfer generalization bounds from interpolating kernel methods to neural networks. Second, a lower bound on the NN test loss gives -in the regime of large NN initialization variance -examples of NNs, trained under gradient descent after random initialization, that generalize badly despite (close to) zero training loss.

4.1. TRANSFERRING GENERALIZATION BOUNDS TO NEURAL NETWORKS:

From Thm. 2, it is now almost straightforward to connect the test error of the interpolating kernel method to the test error (w.r.t. squared loss) of the neural network.

Theorem 3 Let f NN be the fully-converged solution of a strongly overparametrized NN (1) with L layers and m hidden units per layer and initialization variance σ 2 , where all layers θ have been trained using gradient descent under squared loss on a consistent dataset D after random initialization θ 0 ∼ N (0, 1). Then it holds with probability at least 1 -2δ over the random initialization that

L NN test ≤ L int test + σ L / √ δ + O(1/ √ m).

To prove this theorem we use the fact that

|f lin (x) -f NN (x)| ≤ O(1/ √ m) holds uniformly for all x ∈ B d 1 (0) and that f lin (X test ) -f int (X test ) 2 /N test ≤ σ 2L

/δ holds with high probability; the last fact follows in the same way as in the proof of Lemma 3 by, instead of applying Markov's inequality to the squared deviation for a single test point x, applying it to averaged sum of squared deviations over all test points. Then, we make use of the triangle inequality and evaluate the square to prove the desired result. This inequality now enables us to transfer bounds for the test error from one formalism to the other and gives us a new grip on finding bounds on the test error for NNs. Good generalization bounds on minimum complexity interpolating models have been established in (Liang & Rakhlin (2018) ). On the other hand, Belkin et al. (2018b) showed that conventional bounds in the underparametrized regime might be of little use for minimum complexity kernel methods.

4.2. ZERO NN TRAINING ERROR ALONE DOES NOT IMPLY SMALL TEST ERROR:

In the past many people trained neural networks with the goal to achieve a low training error with the hope that this might also result in a low test error. The initialization was merely thought of as a necessity to converge faster to a low training error or to give a favourable starting point for the training dynamics. We now show with our analysis that random initialization alone is not enough, but instead we need certain constraints on the initialization to be able to give generalization bounds.

Before we prove this let us derive some simple properties. Let Φ(x, σ) be the feature vector of a NN initialized with NTK-parametrization with initialization variance σ. Then it is easy to see that Φ(x, σ) = σ L Φ(x, 1) because of the homogeneity of the derivative w.r.t constants. Further we define J(X test , σ)

:= 1 L Φ(X test ) 1 -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) θ 0 2 /

√ N test and J(X test ) := J(X test , 1). With this notation, we can show the following theorem:

Theorem 4 Let f NN be the fully-converged solution of a strongly overparametrized NN (1) with L layers and m hidden units per layer and initialization variance σ 2 , where all layers θ have been trained using gradient descent under squared loss on a consistent dataset D after random initialization θ 0 ∼ N (0, 1). Then it holds with probability at least 1 -δ over random initialization

L NN test ≥ |σ L J(X test ) -O(1/ √ m) -L int test |

For general datasets D, it is highly likely that J(X test , σ) = 0. To prove the statement we use the decomposition of Lemma 2. Now instead of making the variance small to bound the second term we can use the homogeneity of Φ(x, σ) to show that J(X test , σ) = σ L J(X test ). Now using twice the reverse triangle inequality we and the lemma 1 get the desired result. Now looking at the expression of Thm. 4 and using again the homogeneity of Φ we can see that L int test does not change when changing σ because the different factors cancel out and the term corresponding to the linear approximation should also not impact the result too much as long as we choose m large and do not increase σ significantly. Respecting this constraint we can now easily increase the variance to make the second term as large as we want. Thus, by increasing σ we are able to increase the test loss arbitrarily while not changing the training error at all because the second term in our decomposition does not contribute during training time. This shows that without considering an explicit initialization the training error of a NN might tell us very little about the actual test loss. Thus, our results underline the importance of finding good initialization strategies.

5. EXPERIMENTS:

In this section we give experimental evidence for our theoretical findings: the influence of NN initialization on the generalization performance of fully-trained wide NNs (Sect. 4) as well as the link to interpolating kernel methods in the first place (Sect. 3).

We perform experiments (Fig. 1 ) in three different settings3 : (1) A toy experiment, fitting a 2Ddataset with N = 10 datapoints sampled uniformly from f (x) = exp(-x 2 ), x ∈ [ -1, 1] 2 , with a NN of 1 hidden layer (L = 2) with m = 10000 neurons; here, the "wide NN" requirement m N is satisfied to the highest degree. (2) Fitting N = 100 MNIST digits 0 and 1 by a single-hidden-layer network with m = 4000 hidden units. In line with our framework and with other works on wide NNs (e.g., (Du et al., 2018b; Arora et al., 2019b )), we fit these classification labels using squared error loss. (3) Fitting MNIST digits 0 and 1 with 2 hidden layers (L = 3) of m = 1500 neurons each. The reason for using small training sets is to maintain the overparametrization limit m N .

In each setting we vary the variance σ 2 of the NN initialization, performing always 10 repetitions.

We stop training at the same "close to zero" training loss for all σ 2 , chosen in each setting so small that empirically the test loss remains almost constant (while the training loss is decreasing further).

The left panels (Fig. 1 ) show the test losses of the trained NN vs. the linearization f lin of the initialized NN (Lemma 1), at different NN initialization variances σ 2 . Both curves are very close, confirming the basic premise that the approximation of wide NNs by linear models is quite accurate for our NNs, even though we are not in the guaranteed overparametrized regime (e.g., m ≥ O(N 6 ) from Du et al. ( 2018b)). The approximation is better for the two MNIST experiments, where we are actually further from the this regime, whereas the larger repetition variances for the 2D bump are due to sampling the training inputs from a small volume of low dimension.

Most strikingly, the test loss of the trained NN depends heavily on its weight initialization θ 0 even though all NNs have been trained to the same low training error. This is shown by the increase at larger initialization variance σ 2 . Via the link to f lin , this can be understood from the test behavior of strongly overparametrized linear models f lin with big initialization θ 0 (Eq. ( 14) and Sect. 4.2).

The right panels underline this behavior predicted by our theory quantitatively: The test error of trained ReLU-NNs as well our our analytical lower bound on it (Thm. 4) both behave at big σ 2 like Θ(σ 2L ) (see also Eq. ( 12)), with L the number of layers. This in particular shows how one can find natural examples of (overparametrized) NNs, i.e. initialized randomly and trained by gradient descent, that interpolate the training data perfectly, yet with arbitrarily large test error.

In the regime of small σ 2 on the other hand, the right panels confirm Thm. 3, namely that the test error is close to the (low) test error of the interpolating kernel method f int Φ (Eq. ( 5)). Note that its test error in the figures is almost constant in σ 2 , as f int does not depend directly on the weight initialization θ 0 but on the (many and random) NN features φ(x), which via the kernel Φ only determine function smoothness. This is a more reasonable modeling behavior than that of wide NNs initialized at large σ 2 .

6. DISCUSSION:

In this paper we make the bias explicit that is imposed by gradient descent and random initialization on fully-connected strongly overparametrized NNs. We achieve this by decomposing a gradient descent-trained network with (close to) zero training error into an interpolating kernel term and a second term that vanishes on training inputs but not on test inputs.

Because this second term is independent of the training labels but depends strongly on the initialization, it is favorable for good generalization to reduce its impact. This can be achieved by choosing the variance of the NN initialization to be sufficiently small, as we have shown. On the other hand, too small an initialization variance limits the expressivity of the NN and its features.

Our work furthers the practical understanding of modern NNs by utilizing recent theoretical advances in strongly overparametrized NNs. In particular, we saw that the initialized weights are not merely a starting point for training, but significantly influence the performance of the NN trained to 4.2, as the linearized approximation of our NNs is quite accurate (see also Sect. 3). Right panels: Quantitatively, at big σ 2 , the NN test loss as well as its lower bound both grow as Θ(σ 2L ) (Thm. 4). At small σ 2 , on the other hand, the NN test error is close to that of the interpolating kernel method f int Φ (Eq. ( 5)), confirming Thm. 3. zero empirical error. Similarly, focussing on the first term of our decomposition, learning-theoretical results on interpolating kernel methods lead to a better understanding of puzzling features of NNs.

Our work provides intuition for the "essentially no barriers" phenomenon first observed in Draxler et al. (2018) ; Garipov et al. (2018) : When scaling all initialization weights θ 0 by a common continuous factor (like σ), all of the finally trained NNs f NN σ are connected continuously and their training error stays at (close to) zero. This is because Φ changes by only a scalar factor σ 2L , which continously changes the second (test) term in the linearization f NN ≈ f lin (Eq. 14) but cancels out of the first (training) term. It would be interesting for further work to investigate varying θ 0 more generally in the hypersurface of zero training loss, which induces non-scalar changes in Φ.

The connection of NNs with minimum complexity interpolating kernel methods opens up more directions for future research. An important next step would be to go from strongly overparametrized NNs to merely overparametrized NNs, where the (infinitesimal) kernel dynamics effected by NN training is to be understood in addition to fixed-kernel interpolation. Further, the connection remains to be made explicit for CNNs and RNNs, as well as for other (classification) loss functions, whose strict minimization often requires infinite weights. For SGD at least, instead of full gradient descent, the basic decomposition into a training part and its orthogonal component (Eq. 14) remains valid.

Another direction needed to fully leverage the link to NNs made explicit in Sect. 4.1, is to understand better the properties -e.g. the generalization behavior (Liang & Rakhlin, 2018 ) -of interpolating kernel methods, whose close connection to large modern-day NNs we have highlighted in this paper.

Taking the derivative with respect to θ and setting it to zero gives us θT = λ T Φ(X). Now plugging this into the Lagrangian we get

L(λ) = - 1 2 λ T Φ(X)Φ(X) T λ + λ T Y

Minimizing this term with respect to λ we get λ = Φ(X)Φ(X) T -1 Y . Inserting this relation we get

θ = Φ(X) T Φ(X)Φ(X) T -1 Y B CONNECTING THE INTERPOLATING MINIMUM 2-NORM SOLUTION TO THE MINIMUM RKHS-NORM INTERPOLATING KERNEL METHOD

Building upon the previous section we now want to show that a linear predictor where the weights are the solution of a minimum 2-norm interpolation with feature map Φ(x) is equivalent to a minimum complexity interpolating kernel method with

f = 1 2 arg min f ∈H K f H K subject to Y = f (X),

were the kernel of the RKHS is given by K(x, y) = Φ(x), Φ(y) and Φ is the same function used for the linear ansatz in the minimum 2-norm solution.

We can find the solution of the minimum complexity interpolating kernel method by minimizing the following Lagrangian:

L(f, λ) = 1 2 f H K + λ T (Y -f (X))

Because the RKHS-Norm is bounded by definition (for non degenerate data) we can use the representer theorem (f = N i=1 α i K(x i , .) with α i ∈ R) and the reproducibility of the kernel K to arrive at

L(α, λ) = 1 2 α T K(X, X)α + λ T (Y -K(X, X)α)

Now taking the derivative with respect to α we get α = λ. Inserting this into our Lagrangian we get

L(λ) = - 1 2 λ T K(X, X)λ + λ T Y

This Lagrangian is equivalent to the Lagrangian of the minimal two norm solution in the dual representation and due to the convexity of the problem the solutions must also be equivalent.

C LEMMA 2:

In this section we want to give the proof for lemma 2. Let f N N be a fully-connected NN as defined in the notation section. For simplicity we switch from the NTK initialization to the He et al. (2015) initialization. The proof is equivalent in both parametrizations. The only difference is that we can suppress some of the scaling factors in the initialization of the weights. We again assume a ReLUactivation function. Before showing the actual lemma we will first show the following expression which we will use to proof the lemma

EQUATION

were the scalar product is defined as the Hilbert-Schmid-product with A, B = tr A T B and W = diag(W L , ..., W 1 ) is the diagonal matrix with all the weigh matrices on the diagonal. We restrict ourselves to networks without biases but the proof with biases is almost identical.

To prove the lemma we can first proof that ∀l ∈ {1, ..., L} we have f

N N W l (x) = W l , ∇ W l f (x)

. If this statement is true it is easy to see that the formula above is also true.

A ReLU Network can be written in the following way

EQUATION

where l h l (x)≥0 is the diagonal matrix with the step function on the diagonal corresponding to the components of h l (x) ≥ 0. Now, we can define a(W L:l , x)

:= W L l h L-1 (x)≥0 W L-1 ...l h l (x)≥0 and b(W l:1 , x) := l h l-1 (x)≥0 W l-1 ... to get f N N W (x) = a(W L:l , x)W l b(W l:1 , x)

Using the fact that ∂aij ∂(W l )mn = 0 because we define

∂(I h l (x)=0 )ij ∂(W l )mn = 0.

Defining the derivative of the stepfunction at the step to be zero is also done in most of the standard frameworks like tensorflow. Additionally, the region where the stepfunction jumps is a d -1 dimensional subspace and thus the expression holds almost everywhere anyway. Now, since our weights do not move much and only a neglectable number change the step function we can just define this derivative to be zero. The derivative of b(W l:1 , x) with respect to W l is clearly zero because it does not depend on W l .

Using this we can see that almost everywhere

W l , ∇ W l f N N (x) = nm (W l ) nm ∂f N N W (x) ∂(W l ) nm = a(W L:l , x)W l b(W l:1 , x) = f N N (x),

Applying this result L-times we have get that almost everywhere f N N (x) = 1 L W, ∇ W f (x) . Thus, we get almost everywhere

f θ0 (x) = 1 L θ 0 , Φ(x)

and therefore also

f N N θ0 (x)-Φ(x)Φ(X) T (Φ(X)Φ(X) T ) -1 f N N θ0 (X) = 1 L Φ(x) l -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) θ 0 D LEMMA 3

In this section we want to show the statement of lemma 3. We want to show that

P 1 L |Φ(x) I -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) θ 0 | < σ 2L δ ≥ 1 -δ

First, let us focus on the operator I -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) . The second term has the properties of a projection operator (symmetric and equal to its square). Thus, all the eigenvalues are either zero or one. The operator projects M -dimensional vectors into a N -dimensional subspace and then back to the M -dimensional space. Thus, it has N eigenvalues which are one and all the rest are zero. Now, looking at the whole operator we see that it has M -N eigenvalues which are one and the rest are zero. Therefore, we can neglect the projection operator because this will only make our result larger due to the neglected N dimensions which would be projected out.

Using this and the formula we have derived in th proof for lemma 2 we get almost everywhere

1 L |Φ(x) I -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) θ 0 | ≤ | 1 L Φ(x)θ 0 | = |f NN θ0 (x)|

Now, we can use Markov's inequality to arrive at

P 1 L Φ(x) I -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) θ 0 2 < E |f NN θ0 (x)| 2 /δ ≥ 1 -δ

Thus, the only thing that is left to do is to show that E f NN θ0 (x) ≤ σ 2L . To start of we can again write the Neural Network as a product of matrices f NN θ0 (x) = W L l h L-1 (x)≥0 W L ...l h 0 (x)≥0 W 1 x. We use a similar approach as Du & Hu (2019) and generalize this approach to non-linear networks. We make use of the fact that for for a random matrix A ∈ R d1×d2 with i.i.d. N (0, 1) entries and a arbitrary non-zero vector v ∈ R d1 /0 the distribution of Av 2 v 2 is χ 2 d2 distributed. Additionally, it is easy to see that

l Av≥0 Av 2 Av 2 ≤ 1. Defining Z i =

Wil h i-1 (x)≥0 Wi-1...l h 1 (x)≥0 W1x 2 l h i-1 (x)≥0 Wi...l h 1 (x)≥0 W1x 2

we thus get

E |f NN θ0 (x)| 2 ≤ E[Z 1 . . . Z L ] = Π L l=1 E[Z l ] = σ 2L .

Inserting the result E [f 0 (x, θ 0 )] we get

P 1 L |Φ(x) I -Φ(X) T (Φ(X)Φ(X) T ) -1 Φ(X) θ 0 | < σ L / √ δ ≥ 1 -δ.

E DETAILS ON OUR EXPERIMENTS:

In this section we want to talk about how we conducted the experiments. Since the experimental settings are not very complicated this should be more than enough to easily reproduce our results.

As described in the experimental section we used three different setting to investigate the dependence of the Neural Network solution on the initialization variance. We used a standard He et al. (2015) initialization which is, as mention in the notation section, equivalent to the NTK-parametrization and standard normal initialization when we rescale the learning rate. Since we choose the theoretical value of η = 1.9 λmax(Ψ(X)Ψ(X) T ) for the learning rate the settings are equivalent. For all experiments we use Tensorflow with ReLU-activation functions, set the bias initially to zero and train with GD. For all networks we use a low number of training data points to make sure that we are in or close to the over-parametrized regime. This is necessary to assure the validity of the linearisation.

The first Network for the 2D exponential Bump is a single hidden layer NN with 10000 hidden units. We use 10 training data points and 300 testing points. We train 10 networks for each variance value and until the training loss drops first bellow a loss of 10 -5 . The final training loss of the curves plotted in the experimental section can be seen in figure 2 . We choose the stopping point for the training error in such a way that even though the training loss is still decreasing further the test loss is almost constant. 

:

We consider one-dimensional outputs in this paper, although this is easy to generalize.

Our arguments hold for general fully-connected NNs, too in the "wide limit" min l=1,...,L-1 {m l } → ∞.

More details on the settings needed to reproduce the experiments can be found in Appendix E.

Figures and Tables:

Soudry et al. (2018);Glorot & Bengio (2010);Daniely et al. (2016);Rahaman et al. (2018);Oymak & Soltanolkotabi (2018), which all elucidate certain aspects of the implicit bias. Our work is unique in that it investigates the effect (implicit bias of (S)GD and random initialization) on the test error of trained NNs, exhibiting NNs where -despite always vanishing training error -the test error can be either very good or arbitrarily large. Without gradient descent training and in a deterministic fashion, such examples have been constructed in(Mücke & Steinwart, 2019).

; Du et al. (2018b); Lee et al. (2019)) shows that the solution θ t of the training dynamics t = 0, 1, . . . stays close to the initialization θ 0 , and the deviation from a linear model vanishes with growing width m: Lemma 1 (Lee et al. (2019); Du et al. (2018b)) Denote the linearization of an NN (1) around its initialization θ 0 by f lin θ

Figure1: Experiments on the 2D bump f (x) = exp(-x 2 ) (top row) and on MNIST with L = 2 and L = 3 layers, respectively (last two rows). Left panels: The test error of trained wide NNs as well as their linearized versions increases significantly with the NN initialization variance σ 2 , despite training all NNs to the same very low training error. This is explained by the theory in Sect. 4.2, as the linearized approximation of our NNs is quite accurate (see also Sect. 3). Right panels: Quantitatively, at big σ 2 , the NN test loss as well as its lower bound both grow as Θ(σ 2L ) (Thm. 4). At small σ 2 , on the other hand, the NN test error is close to that of the interpolating kernel method f int Φ (Eq. (5)), confirming Thm. 3.

Figure 2: Benchmark experiment with a 2D-Gaussian Bump

Figure 3: Experiments with MNIST for a NN with 1 and 2 hidden layers

Peer review (iclr2020):
The paper considers the impact of initialization bias on test error in strongly overparameterized neural networks. The study uses tools from recent literature on the generalization of overparameterized neural networks, i.e. neural tangent kernels and interpolating kernel method, to provide useful insights on how the variance of weights initialization affects the test error. I have a few questions about theoretical results, but the paper has a convincing experiment that supports its theoretical claims. Addressing the following points will improve the exposition of the paper. 1. Please provide a little hint on how Lemma 2 rewrites the equation (13) for linearized function for easier readability without referring to the Appendix. 2. In the case of cross-entropy error, would the effect be similar? Could this be verified with a similar experiment as for MSE? 3. To what extent this result is observed in not as strongly overparameterized settings? In other words, it would be interesting to see what happens if you fix the architectural choice while increasing the number of training parameters, how long does the test error effect persist? Minor remark: - a few typos are present on pages 4, 5, 7, 8