Title: TOWARDS UNDERSTANDING THE TRUE LOSS SURFACE OF DEEP NEURAL NETWORKS USING RANDOM MATRIX THEORY AND ITERATIVE SPECTRAL METHODS Anonymous

Abstract: The geometric properties of loss surfaces, such as the local flatness of a solution, are associated with generalization in deep learning. The Hessian is often used to understand these geometric properties. We investigate the differences between the eigenvalues of the neural network Hessian evaluated over the empirical dataset, the Empirical Hessian, and the eigenvalues of the Hessian under the data generating distribution, which we term the True Hessian. Under mild assumptions, we use random matrix theory to show that the True Hessian has eigenvalues of smaller absolute value than the Empirical Hessian. We support these results for different SGD schedules on both a 110-Layer ResNet and VGG-16. To perform these experiments we propose a framework for spectral visualization, based on GPU accelerated stochastic Lanczos quadrature. This approach is an order of magnitude faster than state-of-the-art methods for spectral visualization, and can be generically used to investigate the spectral properties of matrices in deep learning.

1. INTRODUCTION:

The extraordinary success of deep learning in computer vision and natural language processing has been accompanied by an explosion of theoretical (Choromanska et al., 2015a; b; Pennington & Bahri, 2017) and empirical interest in their loss surfaces, typically through the study of the Hessian and its eigenspectrum (Ghorbani et al., 2019; Li et al., 2017; Sagun et al., 2016; 2017; Wu et al., 2017) .

Exploratory work on the Hessian, and its evolution during training (e.g., Jastrzębski et al., 2018) , attempts to understand why optimization procedures such as SGD can discover good solutions for training neural networks, given complex non-convex loss surfaces. For example, the ratio of the largest to smallest eigenvalues, known as the condition number, determines the convergence rate for first-order optimization methods on convex objectives (Nesterov, 2013) . The presence of negative eigenvalues indicates non-convexity even at a local scale. Hessian analysis has also been a primary tool in further explaining the difference in generalization of solutions obtained, where under Bayesian complexity frameworks, flatter minima, which require less information to store, generalize better than sharp minima (Hochreiter & Schmidhuber, 1997) . Further work has considered how large batch vs small batch stochastic gradient descent (SGD) alters the sharpness of solutions (Keskar et al., 2016) , with smaller batches leading to convergence to flatter solutions, leading to better generalization. These geometrical insights have led to generalization procedures, such as taking the Cesàro mean of the weights along the SGD trajectory (Garipov et al., 2018; Izmailov et al., 2018) , and algorithms that optimize the model to select for local flatness (Chaudhari et al., 2016) . Flat regions of weight space are more robust under adversarial attack (Yao et al., 2018) . Moreover, the Hessian defines the curvature of the posterior over weights in the Laplace approximation for Bayesian neural networks (MacKay, 1992; 2003) , and thus crucially determines its performance.

In this paper we use random matrix theory to analyze the spectral differences between the Empirical Hessian, evaluated via a finite data sample (hence related to the empirical risk) and what we term the True Hessian, given under the expectation of the true data generating distribution. 1In particular, we show that the differences in extremal eigenvalues between the True Hessian and the Empirical Hessian depend on the ratio of model parameters to dataset size and the variance per element of the Hessian. Moreover, we show that that the Empirical Hessian spectrum, relative to that of the True Hessian, is broadened; i.e. the largest eigenvalues are larger and the smallest smaller. We support this theory with experiments on the CIFAR-10 and CIFAR-100 datasets for different learning rate schedules using a large modern neural network, the 110 Layer PreResNet.

It is not currently known if key results, such as (1) the flatness or sharpness of good and bad optima, (2) local non-convexity at the end of training, or (3) rank degeneracy hold for the True Hessian in the same way as for the Empirical Hessian. We hence provide an investigation of these foundational questions.

2. RELATED WORK, CONTRIBUTIONS AND PAPER STRUCTURE:

Previous work has used random matrix theory to study the spectra of neural network Hessians under assumptions such as normality of inputs and weights, to show results such as the decreasing difference in loss value between the local and global minima (Choromanska et al., 2015a) and the fraction of negative eigenvalues at final points in training (Pennington & Bahri, 2017) . To the authors knowledge no theoretical work has been done on the nature of the difference in spectra between the True Hessian and Empirical Hessian in machine learning. Previous empirical work on neural network loss surfaces (Ghorbani et al., 2019; Papyan, 2018; Sagun et al., 2017; 2016; Jastrzębski et al., 2018) has also exclusively focused on the Empirical Hessian or a sub-sample thereof.

The work closest in spirit to our work is the spiked covariance literature, which studies the problem of learning the true covariance matrix given by the generating distribution, from the noisy sample covariance matrix. This problem is studied in mathematics and physics (Baik & Silverstein, 2004; Bloemendal et al., 2016b; a) with applications leading to extensively improved results in both sparse Principal Component Analysis (Johnstone & Lu, 2004) , portfolio theory (Laloux et al., 1999) and Bayesian covariance matrix spectrum estimation (Everson & Roberts, 2000) . A key concept from this literature, is the ratio of parameters to data samples and its effect on the observed spectrum; this ratio is mentioned in (Pennington & Bahri, 2017 ), but it is not used to determine the perturbation between the Empirical Hessian and True Hessian. (Kunstner et al., 2019) consider circumstances under which the Empirical Fisher information matrix is a bad proxy to the True Fisher matrix, but do not use a random matrix theory to characterise the spectral perturbations between the two. We list our main contributions below.

• We introduce the concept of the True Hessian, discuss its importance and investigate it both theoretically and empirically.

• We use random matrix theory to analytically derive the eigenvalue perturbations between the True Hessian and the Empirical Hessian, showing the spectrum of the Empirical Hessian to be broadened.

• We visualize the True Hessian spectrum by combining a GPU accelerated stochastic Lanczos quadrature (Gardner et al., 2018) with data-augmentation. Our spectral visualization technique is an order of magnitude faster than recent iterative methods (Ghorbani et al., 2019; Papyan, 2018) , requires one less hand-tuned hyper-parameter, and is consistent with the observed moment information. This procedure can be generically used to compute the spectra of large matrices (e.g., Hessians, Fisher information matrices) in deep learning.

3. PROBLEM FORMULATION:

For an input x ∈ R dx and output y ∈ R dy we have a given prediction function h(•; •) : R dx × R P → R dy , we consider the family of prediction functions parameterised by a weight vector w, i.e., H := {h(•; w) : w ∈ R P } with a given loss (h(x; w), y) : R dx × R dy → R. Ideally we would vary w such that we minimize the loss over our data generating distribution ψ(x, y), known as the true risk.

EQUATION

with corresponding gradient g true (w) = ∇R true (w) and Hessian H true (w) = ∇∇R true (w). If the loss is a negative log likelihood, then the true risk is the expected negative log likelihood per data point under the data generating distribution. However, given a dataset of size N , we only have access to the empirical or full risk

R emp (w) = 1 N N i=1 (h(x i ; w), y i ) (2)

and the gradients g emp (w) and Hessians H emp (w) thereof. The Hessian describes the curvature at that point in weight space w and hence the risk surface can be studied through the Hessian. Weight vectors which achieve low values of equation 2 do not necessarily achieve low values of equation 1.

Their difference is known as the generalization gap. We rewrite our Empirical Hessian as

EQUATION

where2 ε(w) ≡ H emp (w) -H true (w). A symmetric matrix with independent normally distributed elements of equal variance is known as the Gaussian Orthogonal Ensemble (GOE) and a matrix with independent non-Gaussian elements of equal variance is known as the Wigner ensemble.

3.1. FLATNESS AND GENERALIZATION WITH THE EMPIRICAL HESSIAN:

By the spectral theorem, we can rewrite H emp in terms of its eigenvalue, eigenvector pairs, [λ i , φ i ]:

EQUATION

The magnitude of the eigenvalues represent the magnitude of the curvature in the direction of the corresponding eigenvectors. Often the magnitude of the largest eigenvalue λ 1 , or the Frobenius norm (given by the square root of the sum of all the eigenvalues squared (

P i=1 λ 2 i ) 1/2

), is used to define the sharpness of an optimum (e.g., Jastrzębski et al., 2018; 2017) . The normalized mean value of the spectrum, also known as the trace P -1 P i λ 2 i , has also been used (Dinh et al., 2017) . In all of these cases a larger value indicates greater sharpness. Other definitions of flatness have looked at counting the number of 0 (Sagun et al., 2016) , or close to 0 (Chaudhari et al., 2016) , eigenvalues.

It is often argued that flat solutions provide better generalization because they are more robust to shifts between the train and test loss surfaces (e.g., Keskar et al., 2016) that exist because of statistical differences between the training and test samples. It is then compelling to understand whether flat solutions associated with the True Hessian would also provide better generalization, since the True Hessian is formed from the full data generating distribution. If so, there may be other important reasons why flat solutions provide better generalization.

Reparametrization can give the appearance of flatness or sharpness from the perspective of the Hessian. However, it is common practice to hold the parametrization of the model fixed when comparing the Hessian at different parameter settings, or evaluated over different datasets.

4. ANALYZING THE TRUE RISK SURFACE USING RANDOM MATRIX THEORY:

We show that the elementwise difference between the True and Empirical Hessian converges to a zero mean normal random variable, assuming a Lipshitz bounded gradient. Moreover, we derive a precise analytic relationship between the values of the extremal values for both the Empirical and True Hessians, showing the eigenvalues for the Empirical Hessian to be larger in magnitude. This result indicates that the true risk surface is flatter than its empirical counterpart.

4.1. REGULARITY CONDITIONS:

We establish some weak regularity conditions on the loss function and data, under which the elements of ε(w) converge to zero mean normal random variables.

Lemma 1. For an L-Lipshitz-continuous-gradient and almost everywhere twice differentiable loss function (h(x; w), y), the True Hessian elements are strictly bounded in the range -√ P L ≤ H j,k ≤ √ P L.

Proof. By the fundamental theorem of calculus and the definition of Lipshitz continuity λ max ≤ L

EQUATION

Lemma 2. For unbiased independent samples drawn from the data generating distribution and an L-Lipshitz loss the difference between the True Hessian and Empirical Hessian converges element-wise to a zero mean, normal random variable with variance ∝ 1/N .

Proof. The difference between the Empirical and True Hessian ε(w) is given as

EQUATION

By Lemma 1, the Hessian elements are bounded, hence the moments are bounded and using independence and the central limit theorem, equation 6 converges almost surely to a normal random variable P(µ jk , σ 2 jk /N ).

Remark. For finite P and N → ∞, i.e. q = P/N → 0, |ε(w)| → 0 we recover the True Hessian.

Similarly in this limit our empirical risk converges almost surely to our true risk, i.e. we eliminate the generalization gap. However, in deep learning typically the network size eclipses the dataset size by orders of magnitude.3 

4.2. A RANDOM MATRIX THEORY APPROACH:

In order to derive analytic results, we move to the large Dimension limit, where P, N → ∞ but P/N = q > 0 and employ the machinery of random matrix theory to derive results for the perturbations on the eigenspectrum between the True Hessian and Empirical Hessian. This differs from the classical statistical regime where q → 0. We primarily focus on the regime when q 1. In order to make the analysis tractable, we introduce two further assumptions on the nature of the elements ε(w). Assumption 1. The elements of ε(w) are identically and independently Gaussian distributed N (0, σ 2 ) up to the Hermitian condition.

Assumption 2. H true is of low rank r P .

Under assumption 1, ε(w) becomes a Gaussian Orthogonal Ensemble (GOE) and for the GOE we can prove the following Lemma 3. We discuss the necessity of the assumptions in Section 4.3.

Lemma 3. The extremal eigenvalues [λ 1 , λ P ] of the matrix sum A + B/ √ P , where A ∈ R P ×P is a matrix of finite rank r with extremal eigenvalues [λ 1 , λ P ] and B ∈ R P ×P is a GOE matrix with element variance σ 2 are given by

EQUATION

Proof. See Appendix.

Theorem 1. The extremal eigenvalues [λ 1 , λ P ] of the matrix sum H emp , where λ 1 ≥ λ 2 ... ≥ λ P (such that assumptions 1 and 2 are satisfied), where H true has extremal eigenvalues [λ 1 , λ P ], are given by

EQUATION

Proof. The result follows directly from Lemmas 2 and 3.

Remark. This result shows that the extremal eigenvalues of the Empirical Hessian are larger in magnitude than those of the True Hessian: Although we only state this result for the extremal eigenvalues, it holds for any number of well separated outlier eigenvalues. This spectral broadening effect has already been observed empirically when moving from (larger) training to the (smaller) test set (Papyan, 2018) . Intuitively variance of the Empirical Hessian eigenvalues can be seen as the variance of the True Hessian eigenvalues plus the variance of the Hessian, as we show in Appendix F. We also note that if the condition |λ i | > P/N σ is not met, the value of the positive and negative extremal eigenvalues are completely determined by the noise matrix.

4.3. NOTE ON GENERALIZING ASSUMPTIONS:

For Assumption 1, we note that the results for the Wigner ensemble (of which the GOE is a special case) can be extended to non-identical element variances (Tao, 2012) and element dependence (Götze et al., 2012; Schenker & Schulz-Baldes, 2005) . Hence, similarly to Pennington & Bahri (2017) , under extended technical conditions we expect our results to hold more generally. Assumption 2 can be relaxed if the number of extremal eigenvalues is much smaller than P , and the bulk of H true (w)'s eigenspectra also follows a Wigner ensemble which is mutually free with that of ε(w). Furthermore, corrections to our results for finite P scale as P -1/4 for matrices with finite 4'th moments and P -2/5 for all finite moments (Bai, 2008).

5. EFFICIENT COMPUTATION OF HESSIAN EIGENVALUES:

In order to perform spectral analysis on the Hessian of typical neural networks, with tens of millions of parameters, we avoid the infeasible O(P 3 ) eigen-decomposition and use the stochastic Lanczos quadrature (SLQ) algorithm, in conjunction with GPU acceleration. Our procedure in this section is a general-purpose approach for efficiently computing the spectra of large matrices in deep learning.

5.1. STOCHASTIC LANCZOS QUADRATURE:

The Lanczos algorithm (Meurant & Strakoš, 2006 ) is a power iteration algorithm variant which by enforcing orthogonality and storing the Krylov subspace, K m+1 (H, v) = {v, Hv, .., H m v}, optimally approximates the extremal and interior eigenvalues (known as Ritz values). It requires Hessian vector products, for which we use the Pearlmutter trick (Pearlmutter, 1994) with computational cost O(N P ), where N is the dataset size and P is the number of parameters. Hence for m steps the total computational complexity including re-orthogonalisation is O(N P m) and memory cost of O(P m). In order to obtain accurate spectral density estimates we re-orthogonalise at every step (Meurant & Strakoš, 2006) . We exploit the relationship between the Lanczos method and Gaussian quadrature, using random vectors to allow us to learn a discrete approximation of the spectral density.

A quadrature rule is a relation of the form,

b a f (λ)dµ(λ) = M j=1 ρ j f (t j ) + R[f ] (9)

for a function f , such that its Riemann-Stieltjes integral and all the moments exist on the measure dµ(λ), on the interval [a, b] and where R[f ] denotes the unknown remainder. The nodes t j of the Gauss quadrature rule are given by the Ritz values and the weights (or mass) ρ j by the squares of the first elements of the normalized eigenvectors of the Lanczos tri-diagonal matrix (Golub & Meurant, 1994) . For zero mean, unit variance random vectors, using the linearity of trace and expectation

E v Tr(v T H m v) = TrE v (vv N H m ) = Tr(H m ) = P i=1 λ m i = P λ∈D λ m dµ(λ) (10)

hence the measure on the LHS of equation 9 corresponds to that of the underlying spectral density.

The error between the expectation over the set of all zero mean, unit variance vectors v and the Monte Carlo sum used in practice can be bounded (Hutchinson, 1990; Roosta-Khorasani & Ascher, 2015) , but the bounds are too loose to be of any practical value (Granziol & Roberts, 2017) . However, in the high dimensional regime N → ∞, we expect the squared overlap of each random vector with an eigenvector of H, |v T φ i | 2 ≈ 1 P ∀i, with high probability (Cai et al., 2013) . This result can be seen intuitively by looking at the normalized overlap squared between two random Rademacher vectors, where each element is drawn from the distribution ±1 with equal probability, which gives 1/P . The Lanczos algorithm can also be easily parallelized, since it largely involves matrix multiplicationsa major advantage for computing eigenvalues of large matrices in deep learning. We exploit GPU parallelization in performing Lanczos (Gardner et al., 2018) , for significant acceleration.

5.2. WHY TO AVOID MULTIPLE RANDOM VECTORS AND KERNEL SMOOTHING:

In recent work using SLQ for evaluating neural network spectra (Papyan, 2018; Ghorbani et al., 2019) , the Hessian vector product is averaged over a set of random vectors, typically ≈ 10. The resulting discrete moment matched approximation to the spectra in equation 9 is smoothed using a Gaussian kernel N (λ i , σ 2 ). The use of multiple random vectors can be seen as reducing the variance of the vector overlap in a geometric equivalent of the central limit theorem. Whilst this may be desirable for general spectra, if we believe the spectra to be composed of outliers and a bulk as is empirically observed (Sagun et al., 2016; 2017) , then the same self-averaging procedure happens with the bulk. We converge to the spectral outliers quickly using the Lanczos algorithm, as shown in the Appendix by the known convergence theorem 3. Hence all the spectral information we need can be gleaned from a single random vector. Furthermore, it has been proven that kernel smoothing biases the moment information obtained by the Lanczos method (Granziol et al., 2018) . We hence avoid smoothing (and the problem of choosing σ), along with the ≈ 10× increased computational cost, by plotting the discrete spectral density implied by equation 9 for a single random vector, in the form of a stem plot.

5.3. WEAKNESSES OF STOCHASTIC LANCZOS QUADRATURE:

It is possible to determine what fraction of eigenvalues is near the origin by evaluating the weight of the Ritz value(s) closest to the origin. But due to the discrete nature of the Ritz values, which give the nodes of the associated moment-matched spectral density in equation 9, it is not possible to see exactly whether the Hessian of neural networks is rank degenerate. It is also not possible to exactly determine what fraction of eigenvalues are negative. Furthermore, smoothing the density and using quadrature has the problems discussed in Section 5.2. Instead, we determine the fraction of negative eigenvalues as the weight of negative Ritz values which are well separated from the smallest magnitude Ritz value(s), to avoid counting a split of the degenerate mass.

6. EXPERIMENTS:

In order to test the validity of our theoretical results in Section 4.2, we control the parameter q = P/N by applying data-augmentation and then running GPU powered Lanczos method (Gardner et al., 2018) on the augmented dataset at the final set of weights w f inal at the end of the training procedure. We use random horizontal flips, 4 × 4 padding with zeros and random 32 × 32 crops. We neglect the dependence of the augmented dataset on the input dataset in our analysis, but we would expect the effective number of data points N ef f < N due to degeneracy in the samples. For the Emprical Hessian, we use the full training and test dataset N = 50, 000 with no augmentation. We use PyTorch (Paszke et al., 2017) . We train a 110 layer pre-activated ResNet (He et al., 2016) neural network architecture on CIFAR-10 and CIFAR-100 using SGD with momentum set at ρ = 0.9 and data-augmentation. We include further plots and detail the training procedure for the VGG-16 (Simonyan & Zisserman, 2014) In order to investigate whether empirically observed phenomena for the Empirical Hessian such as sharp minima generalizing more poorly than their flatter counterparts also hold for the True Hessian, we run two slightly different SGD decayed learning rate schedules (equation 11) on a 110-Layer ResNet (Izmailov et al., 2018) .

EQUATION

For both schedules we use r = 0.01. For the Normal schedule we use α 0 = 0.1 and T = 225 whereas for the Overfit schedule we use α 0 = 0.001 and T = 1000.

We plot the training curves in the Appendices B, C and include the best training and test accuracies and losses, along with the extremal eigenvalues of the Empirical Hessian and Augmented Hessian in Table 1 for CIFAR-100 and We then run SLQ with m = 100 on the final set of training weights to compute an approximation to the spectrum for both schedules on the full 50, 000 training set (i.e the Empirical Hessian) and with m = 80 on an augmented 1, 500, 000 data-set as a proxy for the True Hessian, which we denote the Augmented Hessian. The 110-Layer ResNet has 1, 169, 972/1, 146, 842 parameters for CIFAR-100/CIFAR-10 respectively, hence for both we have q < 1. We primarily comment on the extremal eigenvalues of the Augmented Hessian and their deviation from those of the Empirical Hessian. We also investigate the spectral density at the origin and the change in negative spectral mass. A sharp drop in negative spectral mass, indicates due to Theorem 1 that most of the negative spectral mass is due to the perturbing GOE. We plot the Empirical Hessian spectrum for the end of training for CIFAR-100 in Figure 1 and the Augmented Hessian in Figure 4 . We plot the same discrete spectral density plots for the Empirical Hessian for CIFAR-10 in Figure 2 and the Augmented Hessian in Figure 3 . We note that the effects predicted by Theorem 1 are empirically supported by the differences between the Empirical and Augmented Hessian. For CIFAR-100 the extremal eigenvalues shrink from [-6.54, 20 .56] (Figure 1 ) for the full dataset to [-0.014, 3.026] (Figure 4 ) for the Augmented dataset. The Empirical Hessian has 43/100 negative Ritz values with a negative spectral mass total of 0.05. The Augmented Hessian has 5/80 negative Ritz values with total negative spectral mass of 0.00016. For CIFAR-10 the extremal eigenvalues shrink from [-1.374, 6 .821] (Figure 2 ) to [-0.0061, 0.5825] (Figure 3 ). The fraction of negative Ritz values shrinks from 35/100 with total weight 0.09 to 3/80 with total mass 0.03. As can be seen from the spectral plots of the full Hessians for both CIFAR-100 and CIFAR-10 in Figures 6 & 5 , and their respective Augmented Hessian 7 & 8, we also observe extreme spectral shrinkage, with the largest eigenvalue decreasing from [281.32, 184.4] to [45.68, 15.6 ] and the smallest increasing from [-7.33, -2.9 ] to [2×10 -5 , -9×10 -4 ] respectively for CIFAR-100/CIFAR-10. We note that despite reasonably extreme spectral shrinkage, the sharpest values of the Overfit schedule augmented Hessians are still significantly larger than those of the Normal schedule, hence the intuition of sharp minima generalizing poorly still holds for the True Hessian. The fraction of negative Ritz values for the Empirical Hessian is 43/100, with relative negative spectral mass of 0.032. For the Augmented Hessian, the number of negative Ritz values is 0. For the CIFAR-10 Overfit Augmented Hessian, there is one negative Ritz value very close to the origin with a large relative spectral mass 0.086, hence it could also be part of the spectral peak at the origin. 

6.2. OVERFITTING SGD:



6.3. RANK DEGENERACY AND RIGHT SKEW:

The vast majority of eigen-directions of the Augmented Hessian are, like those of the Empirical Hessian, locally extremely close to flat. We show this by looking at the 3 Ritz values of largest weight for all learning rate schedules and datasets, shown in Table 3 . The majority of the spectral bulk is carried by the 3 closest weights (≈ 0.99) to the origin for all learning rate schedules and datasets. In instances (CIFAR-100 SGD OverFit) where it looks like the largest weight has been reduced, the second largest weight (which is close in spectral distance) is amplified massively.

We further note that, as we move from the q 1 regime to q < 1, we move from having a symmetric to right skew bulk. Both have a significant spike near the origin. This experimentally corroborates the result derived in (Pennington & Bahri, 2017) , where under normality of weights, inputs and further assumptions the symmetric bulk is given by the Wigner semi-circle law and the right skew bulk is given by the Marcenko-Pastur. The huge spectral peak near the origin, indicates that the result lies in a low effective dimension.

For all Augmented spectra, the negative spectral mass shrinks drastically and we see for the Overfit Schedule for CIFAR-100 that there is no spectral mass below 0 and for CIFAR-10 there is one negative Ritz value (and 79 positive). This Ritz value λ n = -0.00093 has a weight of ρ = 0.085 compared to the largest spike at λ = 0.0001 with weight 0.889, hence we cannot rule out that both values belong to a true spectral peak at the origin.

7. DISCUSSION:

The geometric properties of loss landscapes in deep learning have a profound effect on generalization performance. We introduced the True Hessian to investigate the difference between the landscapes for the true and empirical loss surfaces. We derived analytic forms for the perturbation between the extremal eigenvalues of the True and Empirical Hessians, modelling the difference between the two as a Gaussian Orthogonal Ensemble. Moreover, we developed a method for fast eigenvalue computation and visualization, which we used in conjunction with data augmentation to approximate the True Hessian spectrum.

We show both theoretically and empirically that the True Hessian has smaller variation in eigenvalues and that its extremal eigenvalues are smaller in magnitude than the Empirical Hessian. We also show under our framework that we expect the Empirical Hessian to have a greater negative spectral density than the True Hessian and our experiments support this conclusion. This result may provide some insight as to why first order (curvature blind) methods perform so well on neural networks. Reported non-convexity and pathological curvature is far worse for the empirical risk than the true risk, which is what we wish to descend.

The shape of the true risk is particularly crucial for understanding how to develop effective procedures for Bayesian deep learning. With a Bayesian approach, we not only want to find a single point that optimizes a risk, but rather to integrate over a loss surface to form a Bayesian model average. The geometric properties of the loss surface, rather than the specific location of optima, therefore greatly influences the predictive distribution in a Bayesian procedure. Furthermore, the posterior representation for neural network weights with popular approaches such as the Laplace approximation has curvature directly defined by the Hessian.

In future work, one could also replace the GOE noise matrix ε(w) with a positive semi-definite white Wishart kernel in order to derive results for the empirical Gauss-Newton and Fisher information matrices, which are by definition positive semi-definite and are commonly employed in second order deep learning (Martens & Grosse, 2015) . Our approach to efficient eigenvalue computation and visualization can be used as a general-purpose tool to empirically investigate spectral properties of large matrices in deep learning, such as the Fisher information matrix. 

Jeffrey H Schenker and Hermann:



G LANCZOS ALGORITHM:

In order to empirically analyse properties of modern neural network spectra with tens of millions of parameters N = O(10 7 ), we use the Lanczos algorithm (Meurant & Strakoš, 2006) with Hessian vector products using the Pearlmutter trick (Pearlmutter, 1994) with computational cost O(N T m), where N is the dataset size and m is the number of Lanczos steps. The main properties of the Lanczos algorithm are summarized in the theorems 3,4 Theorem 3. Let H N ×N be a symmetric matrix with eigenvalues λ 1 ≥ .. ≥ λ n and corresponding orthonormal eigenvectors z 1 , ..z n . If θ 1 ≥ .. ≥ θ m are the eigenvalues of the matrix T m obtained after m Lanczos steps and q 1 , ...q k the corresponding Ritz eigenvectors then

λ 1 ≥ θ 1 ≥ λ 1 - (λ 1 -λ n ) tan 2 (θ 1 ) (c k-1 (1 + 2ρ 1 )) 2 λ n ≤ θ k ≤ λ m + (λ 1 -λ n ) tan 2 (θ 1 ) (c k-1 (1 + 2ρ 1 )) 2 (30)

where c k is the chebyshev polyomial of order k Proof: see (Golub & Van Loan, 2012) . Theorem 4. The eigenvalues of T k are the nodes t j of the Gauss quadrature rule, the weights w j are the squares of the first elements of the normalized eigenvectors of T k Proof: See (Golub & Meurant, 1994) . The first term on the RHS of equation 9 using Theorem 4 can be seen as a discrete approximation to the spectral density matching the first m moments v T H m v (Golub & Meurant, 1994; Golub & Van Loan, 2012) , where v is the initial seed vector. Using the expectation of quadratic forms, for zero mean, unit variance random vectors, using the linearity of trace and expectation

EQUATION

The error between the expectation over the set of all zero mean, unit variance vectors v and the monte carlo sum used in practice can be bounded (Hutchinson, 1990; Roosta-Khorasani & Ascher, 2015) . However in the high dimensional regime N → ∞, we expect the squared overlap of each random vector with an eigenvector of H, |v T φ i | 2 ≈ 1 N ∀i, with high probability. This result can be seen by computing the moments of the overlap between Rademacher vectors, containing elements P (v j = ±1) = 0.5. Further analytical results for Gaussian vectors have been obtained (Cai et al., 2013) . The number of Parameters in the V GG -16 network is P = 15291300 so our augmentation procedure is unable to probe the limit q < 1. The training procedure is identical to the PreResNet except for an initial learning rate of 0.05 and T = 300 epochs. Here we also see a reduction in both extremal eigenvalues. 

:

We consider loss surfaces that correspond to risk surfaces in statistical learning theory terminology.

For a fixed dataset, the perturbing matrix ε(w) can be seen as a fixed instance of a random variable.

CIFAR datasets, which have 50, 000 examples, are routinely used to train networks with about 50 million parameters.

We can always consider the transform A -ψ(A)I

the proof of Lemma1(Ledoit & Wolf, 2004) does not depend on the matrix being PSD

Figures and Tables:

Figure 2: Empirical Hessian spectrum for CIFAR-10 SGD.

Figure 3: Augmented 1, 500, 000 data Hessian spectrum for CIFAR-10 SGD.

Figure 6: Empirical Hessian Spectrum C10 OverFit SGD.

Figure 7: Augmented Data 1, 500, 000 Hessian Spectrum C100 OverFit SGD.

Figure 8: Augmented Data 1, 500, 000 Hessian Spectrum C10 OverFit SGD.

Figure 9: SGD CIFAR-100 training and test curves

Figure 20: CIFAR-100 PreResNet110 150'th epoch full dataspectrum

Figure 33: CIFAR-10 VGG16BN 300'th epoch full dataspectrum

in Appendix I.

Table content: <table><tr><td/><td/><td/><td colspan="2">225'th epoch -Full Data Spectrum</td></tr><tr><td>Spectral Density</td><td>10 8 10 5 10 2</td><td/><td/><td/></tr><tr><td/><td>6.54</td><td>2.49</td><td>Eigenvalue Size</td><td>11.53</td><td>20.56</td></tr><tr><td/><td/><td colspan="3">Figure 1: Empirical Hessian spectrum for CIFAR-100 SGD.</td></tr><tr><td/><td/><td/><td colspan="2">225'th epoch -Full Data Spectrum</td></tr><tr><td>Spectral Density</td><td>10 8 10 5 10 2</td><td/><td/><td/></tr><tr><td/><td>1.374</td><td>1.357</td><td>Eigenvalue Size</td><td>4.089</td><td>6.821</td></tr></table>

SCHEDULE -EMPIRICAL AND AUGMENTED HESSIAN ANALYSIS

Table content: <table><tr><td/><td>10 0</td><td/><td colspan="2">1000'th epoch -Full Data Spectrum</td></tr><tr><td>Spectral Density</td><td>10 8 10 6 10 4 10 2</td><td/><td/><td/></tr><tr><td/><td>7.3</td><td>88.9</td><td>Eigenvalue Size</td><td>185.1</td><td>281.3</td></tr><tr><td/><td/><td colspan="3">Figure 5: Empirical Hessian Spectrum C100 OverFit SGD.</td></tr><tr><td/><td/><td/><td colspan="2">1000'th epoch -Full Data Spectrum</td></tr><tr><td>Spectral Density</td><td>10 8 10 5 10 2</td><td/><td/><td/></tr><tr><td/><td>2.9</td><td>59.7</td><td>Eigenvalue Size</td><td>122.2</td><td>184.8</td></tr></table>

Schulz-Baldes. Semicircle law and freeness for random matrices with symmetries or correlations. arXiv preprint math-ph/0505003, 2005. Karen Simonyan and Andrew Zisserman. Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556, 2014. Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc., 2012. CIFAR-100 statistics table including extremal eigenvalues of the Empirical Hessian [λ 1 , λ n ], the True Hessian [λ 1 , λ n ] and test/train accuracies and losses.

Table content: <table><tr><td colspan="8">Lei Wu, Zhanxing Zhu, et al. Towards understanding generalization of deep learning: Perspective of</td></tr><tr><td colspan="6">loss landscapes. arXiv preprint arXiv:1706.10239, 2017.</td><td/></tr><tr><td colspan="8">Zhewei Yao, Amir Gholami, Qi Lei, Kurt Keutzer, and Michael W Mahoney. Hessian-based analysis</td></tr><tr><td colspan="8">of large batch training and robustness to adversaries. In Advances in Neural Information Processing</td></tr><tr><td colspan="4">Systems, pp. 4949-4959, 2018.</td><td/><td/><td/></tr><tr><td colspan="4">A TABLE SUMMARIES</td><td/><td/><td/></tr><tr><td colspan="2">Schedule λ 1</td><td>λ n</td><td>λ 1</td><td>λ n</td><td colspan="3">Test Acc Train Acc Test Loss Train Loss</td></tr><tr><td>Normal</td><td>20.56</td><td colspan="3">-6.54 3.025 -0.014</td><td>77.24</td><td>99.85</td><td>0.99</td><td>0.021</td></tr><tr><td>Overfit</td><td colspan="5">281.32 -7.33 45.68 2 × 10 -5 56.85</td><td>99.62</td><td>1.85</td><td>0.026</td></tr><tr><td colspan="2">Schedule λ 1</td><td>λ n</td><td>λ 1</td><td>λ n</td><td colspan="3">Test Acc Train Acc Test Loss Train Loss</td></tr><tr><td>Normal</td><td>6.82</td><td colspan="3">-1.37 3.025 -0.014</td><td>95.32</td><td>99.98</td><td>0.18</td><td>0.0034</td></tr><tr><td>Overfit</td><td colspan="2">184.4 -2.9</td><td>15.6</td><td colspan="2">-9 × 10 -4 86.46</td><td>99.91</td><td>1.85</td><td>0.026</td></tr></table>

CIFAR-10 statistics table including extremal eigenvalues of the Empirical Hessian [λ 1 , λ n ], the True Hessian [λ 1 , λ n ] and test/train accuracies and losses

Table content: <table><tr><td>Learning Schedule &amp; DataSet</td><td>[ρ 0 , λ 0 ]</td><td>[ρ 1 , λ 1 ]</td><td>[ρ 2 , λ 2 ]</td></tr><tr><td>CIFAR-100 SGD Normal</td><td>[0.92,0.0014]</td><td>[0.05,-0.03]</td><td>[0.018,0.06]</td></tr><tr><td colspan="2">CIFAR-100 SGD Normal (Augmented) [0.84,0.0017]</td><td>[0.1,-0.0007]</td><td>[0.045,0.003]</td></tr><tr><td>CIFAR-10 SGD Normal</td><td>[0.88,0.007]</td><td>[0.085,-0.01]</td><td>[0.023,0.02]</td></tr><tr><td>CIFAR-10 SGD Normal (Augmented)</td><td colspan="2">[0.88,-3.87e-6] [0.073,0.0006]</td><td>[0.03,-0.0007]</td></tr><tr><td>CIFAR-100 SGD Overfit</td><td>[0.65,-0.026]</td><td>[0.34,0.05]</td><td>[0.006,0.58]</td></tr><tr><td colspan="2">CIFAR-100 SGD Overfit (Augmented) [0.86,1.67e-5]</td><td>[0.1,0.0039]</td><td>[0.019,0.018]</td></tr><tr><td>CIFAR-10 SGD Overfit</td><td>[0.77,0.01]</td><td>[0.22,-0.04]</td><td>[0.005,0.32]</td></tr><tr><td>CIFAR-10 SGD Overfit (Augmented)</td><td colspan="3">[0.89,0.00013] [0.086,-0.00094] [0.016,0.0058]</td></tr></table>

Largest 3 weight Ritz values for different learning schedules and datasets. [ρ i , λ i ] denotes the weight and corresponding eigenvalue ρ 0 ≥ ρ 1 .... ≥ ρ P .

Table content: <table/>

Figure 22: Augmented CIFAR-100 VGG16BN Eigenspectrum Epoch 300, 5 million samples and m = 10 Lanczos steps epoch -Full Data Spectrum.pdf

Table content: <table><tr><td/><td>10 0 10 0</td><td/><td/><td colspan="3">50'th epoch -Full Data Spectrum 300'th epoch -Full Data Spectrum</td></tr><tr><td>Spectral Density Spectral Density</td><td>10 8 10 6 10 4 10 2 10 8 10 6 10 4 10 2</td><td/><td/><td/><td/><td/></tr><tr><td/><td/><td>0.330 3.03</td><td>0.489 0.54</td><td>1.308 1.95</td><td>Eigenvalue Size Eigenvalue Size</td><td>2.126 4.44</td><td>2.945 6.94</td><td>3.763 9.43</td></tr><tr><td colspan="8">6.59 0.344 0.635 0.312 K SPECTRAL PLOTS -CIFAR-10 VGG16BN 0.52 5.54 11.61 Eigenvalue Size 10 8 10 6 10 4 10 2 10 0 300'th epoch -Full Data Spectrum Figure 25: CIFAR-100 VGG16BN 50'th epoch full dataspectrum 17.67 10 0 75'th epoch -Full Data Spectrum Figure 31: CIFAR-10 PreResNet110 300'th epoch full dataspectrum Spectral Density Figure 21: CIFAR-100 VGG16BN 300'th epoch full dataspectrum 1.072 2.488 3.905 5.321 Eigenvalue Size 0.367 1.368 2.369 3.370 Eigenvalue Size Figure 27: CIFAR-100 VGG16BN 125'th epoch full dataspectrum 10 8 10 5 10 2 Spectral Density 10 6 10 4 10 2 10 0 0.254 0.820 1.386 1.952 Eigenvalue Size 10 9 10 6 10 3 000'th epoch -Full Data Spectrum Spectral Density Figure 26: CIFAR-100 VGG16BN 75'th epoch full dataspectrum 10 8 10 5 10 2 Spectral Density 5 million batch eigenvalue density 10 0 125'th epoch -Full Data Spectrum 1.647 0.295 1.057 2.410 3.762 Eigenvalue Size Spectral Density 10 9 10 6 10 3 Figure 32: CIFAR-10 VGG16BN 0'th epoch full dataspectrum Spectral Density 300'th epoch -Full Data Spectrum</td><td>23.73 4.372 6.737 2.518 5.114</td></tr><tr><td>Spectral Density Spectral Density</td><td>10 9 10 6 10 3 10 0 10 8 10 6 10 4 10 2 10 0</td><td>1.522 1.110</td><td>0.430 0.454</td><td colspan="3">2.382 000'th epoch -Full Data Spectrum 4.335 Eigenvalue Size 150'th epoch -Full Data Spectrum 2.017 3.581 Eigenvalue Size</td><td>6.287 5.145</td><td>8.239 6.709</td></tr><tr><td/><td/><td>2.799</td><td colspan="5">1.590 Figure 28: CIFAR-100 VGG16BN 150'th epoch full dataspectrum 0.382 0.827 2.036 Eigenvalue Size</td><td>3.245</td></tr><tr><td>Spectral Density Spectral Density</td><td>10 0 10 9 10 6 10 3 10 0 10 8 10 6 10 4 10 2</td><td>0.595 6.37</td><td colspan="5">150'th epoch -Full Data Spectrum Figure 23: CIFAR-100 VGG16BN 0'th epoch full dataspectrum 0.081 0.757 1.434 2.110 Eigenvalue Size 25'th epoch -Full Data Spectrum 0.60 5.17 10.93 16.70 Eigenvalue Size Figure 29: CIFAR-100 VGG16BN 150'th epoch full dataspectrum</td><td>2.786 22.47</td></tr><tr><td colspan="8">Figure 24: CIFAR-100 VGG16BN 25'th epoch full dataspectrum J SPECTRAL PLOTS -CIFAR-10 PRERESNET</td></tr><tr><td/><td>10 0</td><td/><td/><td colspan="3">50'th epoch -Full Data Spectrum</td></tr><tr><td>Spectral Density</td><td>10 8 10 6 10 4 10 2</td><td/><td/><td/><td/><td/></tr><tr><td/><td/><td>2.273</td><td>0.736</td><td>0.802</td><td>Eigenvalue Size</td><td>2.340</td><td>3.878</td><td>5.416</td></tr><tr><td/><td/><td/><td colspan="5">Figure 30: CIFAR-10 PreResNet110 50'th epoch full dataspectrum</td></tr></table>

Peer review (iclr2020):
The paper compares the true hessian and the empirical hessian of the loss function, showing the spectrum of the Empirical Hessian is generally broadened, and proposes a way to visualize the spectrum. The paper is well written with sound theoretical reasoning and empirical validation. However, the reviewer has the following concerns: As to the assumption in this paper, it is too strong to assume independence between different elements of Hessian, Wigner ensemble seems to be a better model, since the only independence only comes from samples instead of the elements of the Hessian. This would require a completely different analysis to deal with. Furthermore, it seems that Lemma 1 and Theorem 1 are asymptotic results. I don’t see how these are related to the empirical validation, since in experiments the compared dataset are CIFAR-10/CIFAR-100 and augmented CIFAR-10/CIFAR-100, which only differs in the number of samples. It makes more sense to have a non-asymptotic analysis here.