Deep Curvature Suite
Authors: Diego Granziol
Xingchen Wan
Timur Garipov
Venue:  (2021)
Keywords: Hessian computation
Deep Learning
Loss Curvature
Lanczos

Abstract: The curvature of the loss, provides rich information on the geometry underlying neural networks, with applications in second order optimisation and Bayesian deep learning. However, accessing curvature information is still a daunting engineering challenge, inaccessible to most practitioners. We hence provide a software package the \textbf{Deep Curvature Suite}, which allows easy curvature evaluation for large modern neural networks. Beyond the calculation of a highly accurate moment matched approximation of the Hessian spectrum using Lanczos, our package provides: extensive \emph{loss surface visualisation}, the calculation of the \emph{Hessian variance} and \emph{stochastic second order optimisers}. We further address and disprove many common misconceptions in the literature about the Lanczos algorithm, namely that it learns eigenvalues from the top down. We prove using high dimensional concentration inequalities that for specific matrices a single random vector is sufficient for accurate spectral estimation, informing our spectral visualisation method. We showcase our package practical utility on a series of examples based on realistic modern neural networks such as the VGG-$16$ and Preactivated ResNets on the CIFAR-$10$/$100$ datasets. We further detail $3$ specific potential use cases enabled by our software: research in stochastic second order optimisation for deep learning, learning rate scheduling using known optimality formulae for convex surfaces and empirical verification of deep learning theory based on comparing empirical and theoretically implied spectra.

Under review as a conference paper at ICLR 2021
DEEP CURVATURE SUITE
Anonymous authors
Paper under double-blind review
ABSTRACT
The curvature of the loss, provides rich information on the geometry underly- ing neural networks, with applications in second order optimisation and Bayesian deep learning. However, accessing curvature information is still daunting en- gineering challenge, inaccessible to most practitioners. We hence provide soft- ware package the Deep Curvature Suite, which allows easy curvature evalua- tion for large modern neural networks. Beyond the calculation of highly ac- curate moment matched approximation of the Hessian spectrum using Lanczos, our package provides: extensive loss surface visualisation, the calculation of the Hessian variance and stochastic second order optimisers. We further address and disprove many common misconceptions in the literature about the Lanczos algo- rithm, namely that it learns eigenvalues from the top down. We prove using high dimensional concentration inequalities that for specific matrices a single random vector is sufficient for accurate spectral estimation, informing our spectral visuali- sation method. We showcase our package practical utility on series of examples based on realistic modern neural networks such as the VGG-16 and Preactivated ResNets on the CIFAR-10/100 datasets. We further detail 3 specific potential use cases enabled by our software: research in stochastic second order optimisation for deep learning, learning rate scheduling using known optimality formulae for convex surfaces and empirical verification deep learning theory based on com- paring empirical and theoretically implied spectra.
INTRODUCTION
The success of deep neural networks trained with gradient based optimisers in speech and object recognition (LeCun et al., 2015), has led to an explosion in easy to high performance software implementations. Automatic differentiation packages such as TensorFlow (Abadi et al.] 2016) and PyTorch (Paszke al., 2017) have become widely adopted. Higher level packages, such Keras (Chollet (2155) allow practitioners users to state their model, dataset and optimiser in few lines of code, effortlessly achieving state of the art performance.
However, software for extracting second order information, representing the curvature of the loss point in weight space, has not kept abreast. Researchers aspiring to evaluate curvature infor- mation need to implement their own libraries, which are rarely shared or kept up to date. Naive implementations, which rely on full eigendecomposition (cubic cost the parameter count) are computationally intractable for all but the smallest of models. Hence, researchers typically ignore curvature information or use highly optimistic approximations. Examples in the literature include the diagonal elements of the matrix or of surrogate matrix Chaudhari et al. (2016); Dangel et al. (2019), which we show in Appendix can be very misleading.
MOTIVATION
The eurvature of the loss informs about the local conditioning of the problem (i.e the ratio of the largest to smallest eigenvalues Az). This determines the rate of convergence for first order methods and informs us about the optimal learning and momentum rates (Nesterov 2013). Hence easily accessible curvature information could allow practitioners to scale their learning rates in an optimal way throughout training, instead of relying expert scheduling, investigate this using our software in Section 5.2 Research areas where curvature information features most prominently analyses of the Loss Surface and Newton type optimization methods.

---

Under review as a conference paper at ICLR 2021
2.1 Loss SURFACES
Recent neural network loss surfaces using full eigendecomposition (Sagun et al. 2016 2017) have been limited to toy examples with less than five thousand parameters. Hence, loss surface visual- isation of deep neural networks have often focused on two dimensional slices of random vectors (Li et al., 2017) or the changes in the loss traversing set of random vectors drawn from the d- dimensional Gaussian distribution (Izmailov et al. 2018). It is not clear that the loss surfaces of modern expressive neural networks, containing millions or billions of dimensions, can be well cap- tured in this manner. Small experiments have shown neural networks have large rank degeneracy Sagun et al. (2016) with a small number large outliers. However high dimensional concentration theorems Vershynin (2088) guarantee that even large number of randomly sampled vectors are unlikely to encounter such outliers and hence have limited ability to discern the geometry between various solutions. Other works, that try to distinguish between flat and sharp minima have used the diagonal of the Fisher information matrix (Chaudhari et al. (0066), an assumption we will challenge in this paper. Specifically we show in Appendix E] that diagonal approximations do not capture key properties in synthetic and real neural network examples.
From a practical perspective, specific properties of the loss surface are not captured by the aforemen- tioned approaches. Examples include the flatness as specified by the trace, Frobenius and spectral norm. These measures have been extensively used to characterise the generalisation of solution found by SGD [Wu et al. 2018 Izmailov et al. 2018 He et al. 2019 Jastrzkebski | et al. 2017 2018 Keskar al. 2016). Under Bayesian and minimum description length argument (Hochreiter and Schmidhuber 1997) flatter minima should generalise hetter than sharper minima. The magnitude of these outliers have been linked to poor generalisation performance Keskar et al. (2016) and as a con- KerooraaneeKKkkarr
sequences the generalisation benefits of large learning rate SGD Wu et al. (2018 Jastrzkebski et al. (2017). These properties are extremely easy in principle to estimate, at a computational cost of a small multiple of gradient evaluations. However the caleulation of these properties not typically included in standard deep learning frameworks, which limits the ability of researchers to undertake such analysis.
Other important areas of loss surface investigation include understanding the effectiveness of hatch normalisation(Ioffe and Szegedy 2015). Recent convergence proofs (Santurkar et al., 2018) bound the maximal eigenvalue of the Hessian with respect to the activations and bounds with respect to the weights on per layer basis. Bounds on per layer basis do imply anything about the bounds of the entire Hessian and furthermore it has been argued that the full spectrum must be calculated to give insights on the alteration of the landscape (Kohler et al. 2018
2.2 SECOND ORDER OPTIMISATION METHODS
Second order optimisation methods solve the minimisation problem for the loss, L(w) E R asso- ciated with parameters w E Px11 and perturbation ow E PR1x to the second order in Taylor expansion,
8w argmin L(w 5w) ar 5w) - VL(w)
Sometimes, such as in deep neural networks when the Hessian H V2L(w) E PxP is not posi- tive definite, a positive definite surrogate is used. Note that Equation
[ is not lower bounded unless H is positive semi-defiiite. Often either a multiple of the identity added (known as damping) to the Hessian (H И + Dauphin et al. (2014) or surrogate positive definite approxima- tion to the Hessian H, such as the Generalised Gauss-Newton (GGN) (Martens 2010 Martens and Sutskever 2012) is employed. To derive the GGN, express the loss in terms of the activation L(w) = (f(w) of the output of the final layer f(w). Hence the elements of the Hessian can be written
_
dy dy H(w);
Pa(s(w) (fi(w)) afa(w) Of((w)ff (w) Owy
ao(s(w)) 8 (u((w) Ow,
Ow,
Sw,Ow;
k=0
The first term the LHS of Equation 2 is known the Generalised Gauss-Newtar matrix.
Despite the success of second order optimisation methods using the GGN for difficult problems on which SGD is known stall, such recurrent neural networks (Martens and Sutskever 2012). or
2

---

Under review as a conference paper at ICLR 2021
auto-encoders (Martens 2016). Researchers wanting to implement second order methods such as Vinyals and Povey] 2012 Martens and Sutskever, 2012; Dauphin et al. 2014) face the aforemen- tioned problems of difficult implementation. As small side note, Bayesian neural networks using the Laplace approximation feature the Hessian inverse multiplied by a vector (Bishop 2006).
2.3 CONTRIBUTIONS
We introduce our package, the Deep Curvature suite, a software package that allows analysis and visualisation of deep neural network curvature. The main features and functionalities of our package
are:
Eigenspectrum analysis of the curvature matrices Powered by the Lanczos Meurant and Strakos (2006) techniques implemented in GPyTorch (Gardner al. 2018) and outlined in Section 3, with a single random vector we use the Pearlmutter matrix-vector product trick Pearlmutter (1994) for fast inference of the eigenvalues and eigenvectors of the common curvature matrices of the deep neural networks. In addition to the standard Hessian matrix, we also include the feature for inference of the ingen-information of the Generalised Gauss- Newton matrix, a commonly used positive-definite surrogate to Hessiar
Advanced Statistics of Networks In addition to the commonly used statistics to evaluate network training and performance such as the training and testing losses and accuracy, we support computations of more advanced statistics: For example, we support squared mean and variance of gradients and Hessians (and GGN), squared norms of Hessian and GGN, L2 and L-inf norms of the network weights and etc. These statistics are useful and relevant for wide range of purposes such as the designs of second-order optimisers and network architecture.
Visualisations For all main features above, we include accompanying visualisation tools. In addition, with the sigen-informatte obtained, we also feature visualisations of the loss landscape by studying the sensitivity of the neural network to perturbations of weights.
Second Order Optimisation We include out of the box stochastic lanczos second order optimisers which use the absolute Hessian Dauphin al. (2014) or the Generalised Gauss
Newton Martens (2010) as curvature matrices
Batch Normalisation Support Our code allows for the use of batch normalisation Toffe and Szegedy 2015) in neural networks, which is integrated by default most modern convolutional neural networks. Since the way in watch batch normalisation is calculated has a big effect on the resulting curvature calculations and can be dependent on the order in which the samples are processed, we offer two models of evaluation. Train mode which depends on the order in which the samples are processed and similar to what the optimiser sees during training and evaluation mode, in which the statistics are pre-computed and hence invariant to sample shuffling. We expand on the differences between these methods in Appendix A
Our software makes calculating and visualising curvature information as simple as calculating the gradient at a saved checkpoint in weight-space. The computational complexity of our approach is O(mP), where is the number of Lanczos steps and P is the number of model parameters. This in stark contrast full exact eigendecomposition which has a numerical cost of O(P3), which is infeasible for large neural networks.
2.4 RELATED WORK
Recent work making curvature information more available using diagonal approximations, disallows the use of batch normalisation(Dange et al., 2019). Our software extends seamlessly to support
'The computation of the GGN-vector product is similar with the computational cost of two backward passes in the network. Also, GGN uses forward-mode automatic differentiation (FMAD) in addition to the commonly employed backward-mode automatic differentiation (RMAD). In the current PyTorch framework, the FMAD operation can be achieved using two equivalent RMAD operations. 2 whhic corresponds to the number of moments of the underlying Hessian/GGN density which we wish to match
3

---

Under review as a conference paper at ICLR 2021
Functionality
Deep Curvature Suite
PyHessian
Brih Manod
Tel Stooal
Butuh Manval
Full Naxval
230
Pert Macral
Part Vucval
Hession Density Estimation
GGN Deusity Estimation
202
200
Loss Surface Visualisation
150
152
Gradiemt Variance
120
Hessinn/GGN Variance
Second Order Optimisers
Table I: Table of Contribations: comparison between our released soft- ware and that of another Py Torch based Hessian tool PyHessiaa
Figure 1:
Devistion of the mini-batch Hessian deftyGGN (right)
speciral norm from the full dataset, Figure from Granziol 2020
batch normalization In Appendix E.5 we explicitly compare against their diagonal MC approxima- tions on synthetic and real neural network examples and find these approximations inadequate for spectral analysis. Independent work has also uses the Lanczos algorithm for Hessian computation Yao et al. (2019;; Ghorbani et al. 2019); Papyan (2018); Izmailoy et al. (2019). We compare the functionality tyt of our package to that of PyHsssia Yao al. (2018) in Table 1.
One major extra functionality we offer is the calculation of the GGN. For common activatios and loss functions, such as the cross-entropy loss and sigmoid activation, the generalised Gauss-Newton is equivalent to the Fisher information matrix (Pascanu and Bengio (013). Hence this approximation to the Hessian is interesting in its own right and many theoretical analyses, consider the generalised Gauss-Newton or Fisher information matrices opposed to the Hessian Pennington and Worah (2188; Karakida et al. (2019); Papyan (201)). Hence accessing the GGN spectrum could be major asset to researchers. Another extra functionality we offer is the calculation of the Hessian/GGN variance. This quantity has been used to derive optimal learning rates in Wu et al. (2018) and to predict the effect of minibatching curvature using a random matrix theory model [Granziol, 2020). We include an example figure where the theoretical prediction validated using our software in Figure 1.
2018)|Karakida
3 LANCZOS, MISCONCEPTIONS AND SPECTRAL DENSITY APPROXIMATIONS
The Lanczos Algorithm| Meurant and Strakoš (206)), on which our package is based, (Algorithm in the Appendix) is an iterative algorithm for learning a subset of the digenvalues/eigenveetoss of any Hermitian matrix, requiring only matrix vector products. It is an adaptation of the power iteration method, where the Krylov subspace X (H,v) v) spanfu, He, H is orthogonalised using Schmid, resulting in improved convergence. The Lanczos algorithm, is subject to many misconceptions in the literature, which limits its penetration in the research community. Two of which we state here and explicitly debunk in Appendix B.2
We can learn the negative and interior eigenvalues by shifting and inverting the matrix sign H + - H + pl
Lanezos learns the largest m [As u.) pairs of H E PxP with high probability (Dauphin ct al. 2014)
The key properties of the Lanczos algorithm, which we summarise here (and expand upon in Ap- pendix B are:
When seeded with a vector e the Lanezos algorithm gives discrete spectral density which matches the moments vHu, v, v u.., " v.
If the vector v is zero mean unit variance (such as a Gaussian or Rademacher Hutchinson (1990), this in expectation gives an m moment matched approximation to the spectral density ((A) of the underlying matrix H.
For large dimensional matrix P ->, we expect a single random vector to give an m moment matched approximation of the underlying spectral density. We give the proof for this in Appendix This informs our visualisation technique of using a single random vector.
When seeded with a random vector v, Lanczos converges quickly to the end points of the spectral support of p(A), the spectral density of H.
4

---

Under review as a conference paper at ICLR 2021
4 AN ILLUSTRATED EXAMPLE
example of using the Deep Curvature package. We train a VGG16 net- work on CIFAR-100 for 100 epochs. With the checkpoints generated, we may now compute analyse the eigenspectrum of the curvature matrix evaluated at the desired point in training.TD To evaluate the Hessian/GGN with 100 Lanczos iterations, we run the code given in Table
We give an illustration
Network Training
Eigenspectrum Computation
train_networf
compuee_cigenspecrrmn
dir='VGG16/
dataset="CIFARR000
datasct= latasct='CIFAR100"
data pattppa'a''daaa.
data "atth='datt''.
model="VGG16".
epochs=100. model= VGG166',
"ickppoitt.path==ccccpppppp pt'.
optimizer='SGD',
save.spectrum.pathhe "peccttsssserummm save_cigvec=/Truc
spectrum -ggn'
optimizer.kwargss
I'Ir'; 0.03, 'mom': 0.9, weight.deeay ; 5c-4))
lanczos_iters=
curvature_matrixs ggn lanczoss
Table 2: Training Noural Network and Calculating the Eigenspectrum using the Deep Curvature Suite package
This function call saves the spectrum results (including cigenvalues, cigenvectors and other related statistics) in the save spect rum path path string defined. To visualise the spectrum as stem plot simply the code given in Table with corresponding example Figure shown.
10
Eigenspectrum Visualization
Speciaal Dansity
10
plot.spectrum ('lanczos') path=' sgectrum-ggn path='spectrum-ggn npz') plt.
Eigenvalue Size
Table 3: Eigensspectrum plotting code and corresponding stem plot
Finally, with the eigenvalues and eigenvectors computed, might be interested in knowing how sensitive network is to perturbation along these directions. To achieve this, we first construct loss landscape by setting the number of query points and maximum perturbation to apply. To achieve that, we call the code given Table In this example, we set the maximum perturbation to be
Loss Surface
Train Loss
Test Loss
-000470
muild.loss__andscapp
0.5
40 04.0452
dataset=
"IFARet="CIFARRO0 CIFAR100
-0.0442
sata_path='datall VVGe1=
0.4
-0.0348
2.4
--00136
dist=1., points=21, spectrum path=hessian
0.3
-0 0.0000
2.2
0.0002
checkpoint.path== ct. sove.path='scape..
0.2
1.8656
plot_loss_landsee (' ('landscape -100. mpz plt show
0.1
2.0
2.3449
2.6695
1.8
6.88119
0.0
-0.2-0.1 0.0 0.1
0.2
-0.2-0.1
0.0 0.1 0.2
Table 4: Loss surface visualization along the sharpest Hessian eigenvectors with the Deep Curvature Suite
(dist argument) and number of query points along each direction to be 21 (n.poinnt argument). The corresponding plots in Training and Testing loss are also shown. Although we choose to show loss, the package returns the results for accuracy also, which show in Appendix D
5

---

Under review as a conference paper at ICLR 2021
5 EXAMPLE RESEARCH APPLICATIONS
In this section we list 3 example scientific applications of our software package. We investigate second order optimisation for deep neural networks using our inbuild stochastic Lanczos optimiser, online learning rate scheduling and verifying predictions made using random matrix theory in theo- retical deep learning.
5.1 SIMPLICITY OF SECOND ORDER OPTIMISATION WITH THE DEEP CURVATURE SUITE
To enable researchers experiment with stochastic second order optimisation algorithms for deep neural networks, we implement a Lanczos based optimiser (which takes the absolute Hessian Dauphin al. (2014), Generalised Gauss Newton Martens] (2010) input). The code for ning this optimiser is summarised in Table We plot the training error of the VGG-16 network on
Lanczos Training
SGD Training
train
(rain.networr(
dir= VVG16-CIFAR1000
di=="VGGI-CUFAARIOO/.
datasct= "CIasct="CIFARRO0... data data path='data/ "pppaaateepcccs=1100 epochs==100
dataset="CIFARI000
data_path='dan data/ sppcchs=epppppoo 110 model="VGG16", optimizer kwargs={
model="VGG16', optimizer= LLancGG LancGN optimizer.kwargs-2 'Ir': 1. damping 10,
'SGD',
'Ir': 0.01
"momentum", 0.9,
weight_decay batch 128
weeggh-deeay :
batch-size's 128))
curvature batch size': 128))
Table Comparison of SGD traniing Sccond Order Optimisation using the Deep Curvature Suite
CIFAR- 100 dataset, which has over 16 million parameters and hence is out reach of full inver- sion methods, in Figure 2b] We keep the ratio of damping constant to learning rate constant, where 5 = 10a, for a variety of learning rates in {1,0.1 0.1, 0.01, 0.001, 0.0001} with a batch size of 128 for both the gradient and the curvature, all of which post almost identical performance. We also com- pare against different learning rates of Adam and the best grid searched learning rate of SGD, both of which converge significantly slower per iteration compared to our stochastic Newton methods. We note that in terms of practical optimisation that the cost of running a stochastic second order optimiser is far greater than that SGD. For this example, where we use 20 Lanezos iterations, the per iteration cost scales as s 40 times that of SGD. Hence we expect these algorithms simply to serve as useful baselines against more computationally efficient approximations. One example future direction could be to invert the matrices less frequently as done in Martens and Grosse (2015), reducing computational cost.
5.2 SPECTRAL STOCHASTIC GRADIENT DESCENT WITH THE DEEP CURVATURE SUITE
SGD is incredibly sensitive to the choice of initial learning rate and its scheduling. Whilst adap- tive methods are considered more robust to the initial learning rate choice and its schedule, they have been shown to perform poorly on validation or held out test data Wilson al. (2017). Hence practitioners typically experiment with variety of schedules, such as the step, linear, exponential, polynomial and cosine annealing. One potential application of our codebase to use a stochastic estimate of the curvature and to use established optimal learning and momentum rates for the local quadratic approximation at the given point in weightspace Nesterov (2013), given Equation
2
V V\1 VAp
VAp
A A
V
VXp
Poolyak
VA1 + P
..Nestror
(3)
VA1
+ VAP
Where \1, Ap refer to the largest and smallest eigenvalues of the Hessian respectively. The learning rates and momenta a, P are given for both Polyak and Nesterov momentum respectively. There are several complications the above method. Firstly the Hessian of neural networks is in general not

---

Under review as a conference paper at ICLR 2021
Validation Error
1.0
SSGDMN
0.250
Training Error
1.0
Adam 1r-0.0003
0.8
SWA
Adam 1r000000
0.225
176
178
CN Ir-0.0001
0.5
1r=0.001
GN Lr=0.01
0.4
CN 1r-0.1
0.2
0.0
SGD 0.03
20
60
80
100
120
140
160
180
0
25
50
75
Epoch
Epoch
Emur CHIFAR-100 SSGDMIN) SWA PreResSet-
(b) Train Emor VGG-16 CIFAR-100 Second Order Optimisers
positive definite and hence Ap < (these optimality formulae are strictly for convex methods).
Secondly the loss surface at each point in weight space is likely to be different and it is impractical to recalculate the optimal step size using Equation 3 for every iteration, even if we sub-sample the data. We note from the work of Granziol (2020). that the loss surface when using mini-batch varies to that of the full surface (it is sharper) and hence the learning rates when training using a mini-batch must be optimal for the minibatch surface and not that of the full dataset. As an example use case for our software we consider whether the method of et al. (2018), which combines SGD with Polyak averaging Polyak and Juditsky (1992) at the end of training (which they denote SWA), to improve generalisation performance, could be improved. Instead of fine tuning the linear decay learning rate schedule by hand, we instead learn the learning rate using Equation at regular intervals by estimating the curvature using a sub-sampled Lanczos spectral estimate.
We the preactivated ResNet- 10 the CIFAR-100 dataset for 180 epochs and learn the Polyak and Nesterov learning rates and momenta every 20 epochs using the same mini-batch size as the gradient batch size 128 for 20 Lanezos iterations. Since the smallest Ritz values are very very close to O, which would result in momentum p 1, we use a heuristic to remove the smallest Ritz values, whereby the Ritz value of largest spectral mass has more than 50% of the spectral mass, it is removed and the resulting density renormalised, forming the new spectral density of interest. All methods employ Polyak Avearging at epoch 161 in Izmailov et al. (2018). We plot the results in Figure Pa and show that the Nesterov variant (SSGDMN) of online learning rate and momentum learning performs well against the baseline (SWA). Interestingly the Polyak variant (SSGDM) learns to drop the learning rate at regular intervals, similar to classical training methods, which does not perform as well as keeping the learning rate high and using Polyak averaging at the end of training.
5.3 VERIFYING RANDOM MATRIX THEORY IN DEEP LEARNING USING OUR PACKAGE
The full spectrum of deep neural networks can be used to validate/invaliddate novel theoretical con- tributions. As an example. analysis relating the loss surface to spin-glass models from condensed matter physics and random matrix theory (Choromanska al. 2015b a) rely on a number of unreal- istic assumptions, such as input independence (i.e. for an input x that the feature z; is independent of zy, for images which have spatial correlations, such assumptions are clearly violated). One po- tential verification of the practical applicability of these results, that does not involve checking that assumptions are satisfied, would be to visualise the spectra of large real networks and commonly used datasets. This follows because spin glass models have Hessians given by random matrices (Arous and Guionnet, 1997) with known analytical forms (Tao 2012 Akemann et al. [01I]. Hence statistical tests comparing the real world observed spectra and that underlying the specific spin glass model could be undertaken, showing that the results hold despite the assumptions not being met.
Two primary random matrix eigenvalue distributions which the spectra of many classes of ran- dom matrices converge (and which feature in theoretical deep learning papers) are the semi-circle and the Parchenko-Pasttre densities defined below.
elements X.,, i.d normal H X
elements Hi, j are normal
{+)(1 A o2 v
p(A) 2T
V 2 12 P(A) = 2 2g
(4)
Semi-Circke Law
Marchenko-Passtrr
there are N sammles and parameters and N < P the Marchenko- Baattenko-PPatuurr density has degeneracy fraction 1
x
7

---

Under review as a conference paper at ICLR 2021
10
10
-3
10
10
10
10
10
10
-7
10
10
15
-3.005
0.000
3.005
10
5.2
0.0
16.2
(a) (w) VGG1S8N spech aub
(b) Wigner semi-circle
(4) M (10) PResNes- 10 opoch 300
Figure 3: Neural network Hessian spectra end training for VGG-16BN and PreResNei-110 CIFAR-100 and example of the semi-circle law
In the example of [hhromanska et al. (2015a), the corresponding Hessian spectral denisty is the semi-circle law, shown in Figure 3b This differs greatly from neural network spectra found us- ing our package, such as that of the VGG-16 with BN shown in Figure 3a and the Preactived ResNet-110, for the CIFAR- 100 datasets, shown in Figure 3d Major differences include the rank- degeneracies (we find that over 90% of the spectral mass resides the peak nearest the origin) and outliers. The proof techniques used in Arous and Guionnet 1997), which underlie the claims in Choromanska et al. (2)1 specifically use the large deviation prineiple of the Gaussian Orthogonal Ensemble (which converges the semi-circle law) and hence forbid outliers.
The GGN has been modelled as the Marchenko- Pastur density in the literature (Pennington and Bahri 2017 Sagun et al. 2017 2016). The Residual matrix R(w), which is the difference of the Hessian and the GGN, has been modelled by the Wigner ensemble which gives the semi-circle law in Pennington and Bahri (2017). In Figure 4 we plot the GGN and Residual of VGG-16 network on the CIFAR-100 dataset. The GGN unlike the Marchenko-Pasiie has no spectral gap (the density tails off exponentially) and many outliers. Similarly the shape and rank-degeneracy of the Residual matrix varies significantly from the semi-cirele and contains outliers. Whilst single set of negative results does not invalidate an entire framework, implies that further research required to extend the claims of Choromanska et al. (2015a); Pennington and Bahri (2017) to practical deep learning. These results using our package demonstrate that the underlying spectral densities of the matrices in question vary considerably from those given by the theoretical predictions. This implies that the conclusions made may not necassarily hold in real neural networks and requires further investigation.
10
10"
10
10
10
10
10
10
10 -10
-5
40,1306
0.0000
0.1274
0.0
10.4
10
0.000
6,867
(6) (()))) VGG-16 Epoch 300
(b) Murchenko-
(c) ((w)) VGG-16 Fpoch 300
Figure 4: The Generalised Gauss. Newlon and Residual Matrix various Epochs of Training for the VGG- CIFAR- 100 and sample Marcenko-Past Pastur density
CONCLUSION
We introduce the Deep Curvature suite in PyTorch framework, based the Lanczos algorithm implemented in GPyTorrch (Gardner et al., 2018), that allows deep learning practitioners to learn spectral density information as well as cigenvalue/eigenveeorr pairs of the curvature matrices at specific points in weight space. Our package also allows the user to evaluate the Generalised Gauss Newton spectral density, the variance of both the Hessian and the Generalised Gauss Newton and includes two stochastic Lanczos based optimisers. Together with the software, we also include succinct summary of the linear algebra, iterative method theory including proofs of convergence and misconceptions and stochastic trace estimation that form the theoretical underpinnings of our work. We also give 3 example scientific applications of our work, namely second order optimisation on large scale neural networks. learning rate scheduling using spectral information and evaluation of thcoretical results which make predictions on spectral denisities of the Hessian of the neural networks which they are modelling.
8

---

Under review as a conference paper at ICLR 2021
REFERENCES
Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, al. Tensorflow: system for large- scale machine learning. In /2th {USENIX} Symposium on Operating Systems Design and Imple- mentation ((OSDI) 16), pages 265- 283, 2016.
Gernot Akemann, Jinho Baik, and Philippe Di Francesco. The Oxford handbook of random matrix theory. Oxford University Press, 2011.
Shun-ichi Amari and Hiroshi Nagaoka. Methods of information geometry, volume 191. American Mathematical Soc., 2007.
Anonymous. Towards understanding the true loss surface of deep neural networks using random matrix theory and iterative spectral methods. In Submitted to International Conference on Learn- ing Representations, 2020. /tR.https: pttctpa://ppeneviieee openreview net/ ft//foouuue Forummiddeilggaa22a2.. under review.
G Ben Arous and Alice Guionnet. Large deviations for wigner's law and voiculescu's non- commutative entropy. Probability theory and related fields, 108(4):517-542 1997.
Zhaojun Bai, Gark Fahey, and Gene Golub. Some large-scale matrix computation problems. Journal of Computational and Applied Mathematics, 74(1-2):71-89, 1996.
Christopher M Bishop. Pattern recognition and machine learning. springer, 2006.
Pratik Chaudhari, Anna Choromanska, Stefano Soatto, Yann LeCun, Carlo Baldassi, Christian Borgs, Jennifer Chayes, Levent Sagun, and Riccardo Zecehina. Entropy-SGD: Biasing gradi- ent descent into wide valleys, arXiv preprint arXiv:1611.01888 2016.
François Chollet. Keras. httpp::/ggttubb..cm///ooooleeeet http://hhholoeeeeeeee
2015.
Anna Choromanska, Mikael Henaff, Michael Mathieu, Gérard Ben Arous, and Yann LeCun. The loss surfaces of multilayer networks. In Artificial Intelligence and Statistics, pages 192-204, 2015a.
Anna Choromanska, Yann LeCun, and Gérard Ben Arous. Open problem: The landscape of the loss surfaces of multilayer networks. In Conference Learning Theory, pages 1756-1760, 2015b.
Felix Dangel, Frederik Kunstner, and Philipp Hennig. Backpack: Packing more into backprop.
arXiv preprint 19X1:::912.110888 2019.
Yann N Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio. Identifying and attacking the saddle point problem in nigh-dimensional non-convex op- timization. In Advances in neural information processing systems, pages 2933-2941 2014.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpy- torch: Blackbox matrix- matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, pages 7576-7586. 2018.
Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao. An investigation into neural net optimization via Hessian eigenvalue density. arXiv preprint 1rX1v:1901.10159, 2019.
Gene H Golub and Gérard Meurant. Matrices, moments and quadrature. Pitman Research Notes in Mathematics Series, pages 105-105, 1994.
Gene H Golub and Charles F Van Loan. Matrix computations, volume 3. JHU press, 2012.
Diego Granziol. Curvature is key: Sub-sampled loss surfaces and the implications for large batch training, arXiv preprint arXiv:2006.09099 2020.
Diego Granziol, Binxin Ru, Stefan Zohren, Xiaowen Dong, Michael Osbore,, and Stephen Roberts. Meme: An accurate maximum entropy method for efficient approximations in large-scale ma- chine learning. Entropy, (6)):55,,, 2019.

---

Under review as a conference paper at ICLR 2021
Haowei He, Gao Huang, and Yang Yuan. Asymmetric valleys: Beyond sharp and flat local minima.
arXiv preprint 9rXiv:192..00744 2019.
Sepp Hochreiter and Jtirenn Schmidhuber. Flat minima. Neural Computation, 9()::1-4,, 1997.
Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in SkatistcssSimulaan and Computation, 19(2):433 450, 1990.
Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. arXiv preprint
2015.
Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson. Averaging weights leads wider optima and better generalization. Uncertainty in Artificial In- telligence (UAI), 2018.
Pavel Izmailov, Wesley Maddox, Polina Kirichenko, Timur Garipov, Dmitry Vetrov, and An- drew Gordon Wilson. Subspace inference for bayesian deep learning. Uncertainty in Artificial Intelligence (UAI), 2019.
Stanislaw Jastrzkebski. Zachary Kenton, Devansh Arpit, Nicolas Ballas, Asja Fischer, Yoshua Bengio, and Amos Storkey. Three factors influencing minima in SGD. arXiv preprint arXiv:1711.04623, 2017.
Stanislaw Jastrzkebski, Zachary Kenton, Nicolas Ballas, Asja Fischer, Yoshua Bengio, and Amos Storkey. On the relation between the sharpest directions of DNN loss and the SGD step length. 2018.
Ryo Karakida, Shotaro Akaho, and Shun-ichi Amari. Universal statistics of fisher information in deep neural networks: Mean field approach. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1032-1041. PMLR, 2019.
Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Pe- ter Tang. On large-batch training for deep learning: Generalization gap and sharp minima, arXiv preprint arXiv:1609.048.0 2016.
Jonas Kohler, Hadi Daneshmand, Aurelien Lucchi, Ming Zhou, Klaus Neymeyr, and Thomas Hof- mann. Exponential convergence rates for batch normalization: The power of length-direction decoupling in non-convex optimization. arXiv preprint 1 10694, 2018.
Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. nature, 521(7553):36-444, 2015.
Hao Li, Zheng Xu, Gavin Taylor, and Tom Goldstein. Visualizing the loss landscape of neural nets. arXiv preprint arXiv: 1712.0991, 2017.
Lin Lin, Yousef Saad, and Chao Yang. Approximating spectral densities of large matrices. SIAM Review, 58(1):833)) ::44 65, 2016.
Vladimir Alexandrovich Marchenko and Leonid Andreevich Pastur: Distribution of eigenvalues for some sets of random matrices. Matematicheski Sbornik, 114(4)55775366, 1967.
James Martens. Deep learning via Hessian-free optimization. In ICML, volume 27, pages 735-742, 2010.
James Martens. Second-order optimization for neural networks. PhD thesis, University of Toronto. 2016. URL http:// www.ccc.trrooto..e... htt.://wtt..eeee/aaaaaate/////as
phd_ hnd_mmrreen.. .dd
James Martens and Roger Grosse. Optimizing neural networks with Kronecker-factored approxi-
mate curvature. In International conference on machine learning, pages 2408-2417, 2015.
James Martens and Ilya Sutskever. Training deep and recurrent networks with Hessian-free opti- mization. In Neural networks: Tricks of the trade, pages 479-535. Springer, Sprringr,,,, 2012.
10

---

Under review as a conference paper at ICLR 2021
Gérard Meurant and Zdenčk Strakoš. The Lanezos and conjugate gradient algorithms in finite pre- cision arithmetic. Acta Numerica, 15:471-542, 2006.
Yurii Nesterov. Introductory lectures convex optimization; A basic course, volume 87, Springer Science & Business Media, 2013.
Vardan Papyan. The full spectrum of deep net Hessians at scale: Dynamics with sample size. arXiv preprint ar1iv:1811 1.07062, 2018.
Vardan Papyan. Measurements of three-level hierarchical structure in the outliers in the spectrum of deepnet hessians. arXiv preprint arXiv: 002244, 2019.
Razvan Pascanu and Yoshua Bengio. Revisiting natural gradient for deep networks. arXiv preprint 1rX1v:130,.358, 2013.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary Detito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in Pytorch. 2017.
Barak A Pearlmutter. Fast exact multiplication by the Hessian. Neural computation, 6(1):147-160, 1994.
Jeffrey Pennington and Yasaman Bahri. Geometry of neural network loss surfaces via random matrix theory. In Proceedings of the 34th International Conference on Machine Learning-Volure 70, pages 2798-2806. JMLR. org, 2017.
Jeffrey Pennington and Pratik Worah. The spectrum of the fisher information matrix of a single- hidden-layer neural network. In Advances in Neural Information Processing Systems, pages 5410-5419, 2018
Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM journal on control and optimization, 30(4):838-855 855, 1992.
Levent Sagun, Léon Bottou, and Yann LeCun. Eigenvalues of the Hessian in deep learning: Singu- larity and beyond. arXiv preprint arXiv: /6: 1,00477, 2016.
Levent Sagun, Utku Evci, V Ugur Guney, Yann Dauphin, and Leon Bottou. Empirical analysis of the Hessian of over-parametrized neural networks. arXiv preprint arXiv:1706. (44.00040000444. 2017.
Shibani Santurkar, Dimitris Tsipras, Andrew Ilyas, and Aleksander Madry. How does batch nor- malization help optimization? 2483-2493, 2018.
In Advances Neural Information Processing Systems, pages
Terence Tao. Topics in random matrix theory, volume 132. American Mathematical Soc., 2012.
Shashanka Ubaru and Yousef Saad. Applications of trace estimation techniques.
Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge University Press, 2018.
Oriol Vinyals and Daniel Povey. Krylov subspace descent for deep learning, In Artificial Intelligence and Statistics, pages 1261-1268, 2012.
Eugene P Wigner. Characteristic vectors of bordered matrices with infinite dimensions i. In The Collected Works of Eugene Paul Wigner, pages 524-540. Springer, 1993.
Ashia C Wilson, Rebecca Roelofs, Mitchell Stern, Nati Srebro, and Benjamin Recht. The marginal value of adaptive gradient methods in machine learning. In Advances in Neural Information Processing Systems, pages 4148-4158, 2017.
Lei Wu, Chao Ma, and E Weinan. How sgd selects the global minima in over werr-parameterized learning: A dynamical stability perspective. In Advances in Neural Information Processing Systems, pages 8279-8228 2018.
11

---

Under review as a conference paper at ICLR 2021
Zhewei Yao, Amir Gholami, Qi Lei, Kurt Keutzer, and Michael W Mahoney. Hessian- anal- ysis of large batch training and robustness to adversaries. In Advances Neural Information Processing Systems, pages 4949 2018.
Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael Mahoney. Pyhessian: Neural networks through the lens of the hessian. arXiv preprint arXiv: 2019.
A SUBTLETIES OF BATCH NORMALISATION
Most modern convolutional neural networks use batch normalisation loffe and Szegedy 2015 Whilst refer the reader to loffe and Szegedy (2015) for full description of batch normalisation, a subtlety which is introduced that the output of a network for one sample depends on the samples in the mini-batch (as they are used to compute the statistics). Since these statistics are used to normalise and center the layers, they have a large effect on the curvature calculation. In our software package, we give the user two options to deal with this. We denote one Train Mode and the other Evaluation Mode. We detail them both below
A.1 TRAIN MODE
With batch norm in train mode the output of the network for one sample depends on other samples in the mini-batch. This effect influences the result of hessian-vector product computation on dataset. In particular, if the order of samples in the dataset is changed, then the set of mini-batches also changes and the hessian-vectue product will be different. We take this effect into account and fix the order of samples in the dataset, so that hessian-veede vector product for different networks is computed on the same set of mini-batches. Conceptually, one can use the following interpretation of neural network with batch norm layers in train mode: the input of such network is not a single sample but a mini-batch of samples. Instead of thinking about a dataset of samples we should think about a dataset of mini-batches. With the fixed order of samples, we guarantee that all models are evaluated on the same dataset of mini-batches.
A.2 EVALUATIONMODSS
One can also use BN layers in eval mode with BN statistics computed over the whole dataset. With BN layers in eval mode the output of a network on sample does not depend on other samples in the mini-batch. Our code the mode of BN layers can be specified as a parameter for hessian-vector product computation.
B LEARNING TO LOVE LANCZOS
The Lanezos Algorithm, which our package is based, (Algorithm() [ is iterative algorithm for learning a subset of the eigenvalues/ciieenvetors of any Hermitian matrix, requiring only matrix vector products. It can be regarded as far superior adaptation of the power iteration method, where the Krylov subspace X (H, v) = spanfu u, H H' is orthogonalised using Gram-Schmidt. Beyond having improved convergence the power iteration method (Bai et al.) 1996) by storing the intermediate orthogonal vectors in the corresponding Krylov subspace, Lanczos produces estimates of the eigenvectors and eigenvalues of smaller absolute magnitude, known Ritz vectors/values. Despite its known superiority to the power iteration method and relationship to orthogonal poly- nomials and hence when combined with random vectors the ability estimate the entire spectrum of a matrix, these properties are often ignored or forgotten by practitioners. We compare against diaognal approximations for random and synthetic matrices in Appendix
The Lanczos method be be explicitly derived by considering the optimization of the Rayleigh quo- tient (Golub and Van Loan 2012)
4 HU DTO
r(v)
(5)
over the entire Krylov subspace Kcm (H,,H v) as opposed to power iteration which is a particular H. Despite this, practitioners looking learn the leading
vector in the Krylov subspace
12

---

Under review as a conference paper at ICLR 2021
m 5
m
m 15
m
1.5
1x-10-1 10
3.9x
7.4x 10-2
3.9X 10
6.8x10-1
1,2x10
2.0x 10
1.1
2.7% 10
5.5x10
1.1x10
2.1x10
1.7x10-1
188x101
6.9x10-2
207x10-2
1.01
5.6x 10
10010--
15x10
2.0x
9.2x
8.4x10
7.6x10
6.
Table Lk-1/ Rk-1 For different values spectral gap A1/A2 and iteration number Table from Golub and Lean 2012
eigenvalue, very often resort to the power iteration, likely due to its implementational simplicity. We showcase the power iterations relative inferiority Lanezos with the following convergence theorems
Theorem 1. Let HPP be symmetric matrix with eigenvalues x 2. Xp and corresponding orthonormal eigenvectors 21p.pp. If 2 Om are the eigenvalues of the matrix Tmm obtained after Lanczos steps and 9m. the corresponding Ritz eigenvectors then
(A1 An) (01)
A1
81 2 A1
(1m-1(1+ 2p1))²
(6)
(λ1 - tan²
Ap S 0x I Am +
(cm-1 2p1))
where cn is the Chebyshev polyomial of order k. cos lhl ø1 pi (h1 A2)/(12 An)
Proof. see (Golub and Van Loan 2012).
Theorem 2. Assuming the same notation as in Theorem after m power iteration steps the corre- sponding extremal eigenvalue estimate is lower bounded by
X2 A1
2m-1
01 A (A1 -Ax) tan tan²(0:)
From the rapid growth of orthogonal polynomials such as Chebyshev, expect Lanezos superiority to significantly emerge for larger spectral gap and iteration number. To verify this experimentally. we collect the non identical terms in the equations 6 and of the lower bounds for Yι derived by Lanczos and Power iteration and denote them L and R respectively. For different values of \ and iteration number m we give the ratio of these two quatities in Table As can be clearly seen, the Lanczos lower bound is always closer to the true value, this improves with the iteration number m and relative edge is reduced if the spectral gap is decreased.
B. THE PROBLEM MOMENTS: SPECTRAL DENSITY ESTIMATION USING LANCZOS
In this section we show that that the Lanczos Tri- Tri-Diagonal matrix corresponds to an orthogonal polynomial basis which matches moments v" " and that when v is a zero mean random vector with unit variance, this corresponds to the moment of the underlying spectral density.
Stochastic trace estimation ance random vectors
Using the expectation of quadratic forms, for mean, unit vari-
E. Tr(u w) TrEE, (vo H ) = Tr(H
P
(8)
Σ Am = x"
AED
where have used the linearity of trace and expectation. Hence in expectation over the set of random vectors, the trace of the inner product of u and " equal to the m th moment of the spectral density of H.
13

---

Under review as a conference paper at ICLR 2021
The Lanczos tri tri-diagonal matrix T can be derived from the Moment matrix M, corresponding to the discrete measure (a( satisfying the moments μ4 = v² H'v = X'da(A) (Golub and Meurant 1994)
Lanczos-Stieltjes
1
0 Hu
H
v" HH
v" H
M
vdd Hm-lv
02 H 22m-2pp
and hence for a zero mean unit variance initial seed vector, the eigenvector/eigenvalue pairs of T contain information about the spectral density of H shown section B.1 This is given by the following Theorem
Theorem 3.
The eigenvalues T are the nodes 13 of the Gauss quadrature rule, the weights wy are the squares of the first elements of the normalized eigenvectors of TT
Proof. See Golub and Meurant 1994)
A quadrature rule is relation of the form,
M S F())d)(A) p(((tt)++ R[f]
j=1
for a function f. such that its Riemann-Stieltje integral and all the moments exist on the measure dp(A), on the interval [a, bj and where R[f] denotes the unknown remainder. The first term on the RHS of equation 9 using Theorem can be seen a discrete approximation to the spectral density
matching the first m moments v"" !T" (Golub and Meurant 1994 Golub and Van Loan 2012
For n, starting vectors, the corresponding discrete spectral density is given as
P(A)
=
1 Ne
Σ 1=1 k=1
(TH 1 ((2)
(10)
where T corresponds to the first entry of the eigenvector of the k--th eigenvalue, Ag. of the Lanczos tri-diagonal matrix, T, for the 1-th starting vector (Ubaru and Saad Lin et al. 2016))
B.2 COMPUTATIONAL COMPLEXITY
For large matrices, the computational complexity of the algorithm depends on the Hessian vector product, which for neural networks is O(mNP) where P denotes the number of parameters in the network, m is the number of Lanezos iterations and N is the number of data-points. The full re-orthogonalisation adds two matrix vector products, which is of cost O(m? where typically m² < N. Each random vector used can be seen another full run of the Lanczos algorithm, so for random vectors the total complexity is O(dm P(N + m))
Importance of keeping orthogonality The update equations of the Lanczos algorithm lead tri-diagonal matrix T = Rm whose eigenvalues represent the approximated eigenvalues of the matrix H and whose eigenvectors, when projected back into the the Krylov-subspace, X (H. e), give the approximated eigenvectors of H. In finite precision, it is known (Meurant and Strakoš) 2006) that the Lanezos algorithm fails to maintain orthogonality between Ritz vectors, with corre- sponding convergence failure. In order to remedy this, re-orthonormalise 1996) (as shown in line 9 of Algorithm| [] and observe a high degree of orthonormality between the Ritz eigenvectors. Orthonormality is also essential for achieving accurate spectral resolution as the Ritz value weights are given by the squares of the first elements of the normalised eigenvectors. For the practitioner wishing to reduce the computational cost of maintaining orthogonality, there exist more elaborate schemes (Meurant and Strakož 20066] Golub and Meurant 1994).
each step (Bai al.)
We also explicitly debunk some key persistent myths, given below.
14

---

Under review as a conference paper at ICLR 2021
Algorithm Lanezos Algorithm
1: Input: Hessian vector product ( (Hv), number of steps
2: Output: Ritz eigenvalueeeeigenvecoor pairs {Ae, u;} & quadrature weights 7
3: Set v := 4/ ((='v)
4: Set [3 := 0, Vold == e
5: Set (::
6: for j in I,.., m71 do
7:
Hv Woold T = a w 17
8:
9:
VVT w
10:
3 = v
11:
Vold
12:
U = w/B
V(:, + =
14: T(j +1) = Tÿ + =
B
end for
16: ( =Tig(T)
u, Ve;
18: T; {e 11,,0...0))
We can learn the negative and interior eigenvalues by shifting and inverting the matrix sign μI
H
Lanezos learns the largest m [}; u.] pairs of H PxP with high probability (Dauphin et al. 2014)
Since these two related beliefs are prevalent, we disprove them explicitly in this section, with Theo- and
Theorem 4. The shift and invert procedure H + - H + μI, changes the Eigenvalues of the Tri- diagonal matrix T (and hence the Ritz values) to >; ->; +
Proof. Following the equations from Algorithm 1
w (-H + μL)v & = = Hui + μI
U72 W) 1 (H +uI)e (w Hvi pl)
U72 [H o. Hvjjon 02 === -w2/llus
(1I)
02 o/ (-H H + 2T) = H 02 +
82 ||2]]
Assuming this for m - 1, and repeating the above steps for m we prove by induction and finally arrive at the modified tridiagonal Lanczos matrix
T= + III = ---i+unisis
(12)
Remark. No new Eigenvalues of the matrix H are learned. Although it is clear that the addition gf the identity does not change the Krylov subspace, such procedures are commonplace in code per- taining to papers attempting to find the smallest eigenvalue. This disproves the first misconception.
Theorem 5. For any matrix H E PxP such that X1 > A2 > >.... > Ap and m d < Eimmer As in expectation the of random vectors v the eigenvalues of the Lanezos Tridi- agonal matrix do not correspond to the top m eigenvalues of
Proof. Let us consider the matrix H = H
Am-1 1 Am I. 2
[1 > 0,
Vi <
(13)
A <, Vi m
15

---

Under review as a conference paper at ICLR 2021
Under the assumptions of the theorem, Tr(H) exist no W >> 0 such that
0 and hence by Theorem and Equation 8 there
m
P P
Σ
A
1 S m
(14)
P
satisfied for k = 1, as the LHS is manifestly positive and the RHS is negative. By Theorem this holds for the original matrix H.
Remark. Given that Theorem is satisfied over the expectation of the set of random vectors, which
by the CLT realised by Monte Carlo draws of random vectors as - oo the only way to really span the top eigenvectors is to have selected vector which lies in the dimensional subspace of the P dimensional problem corresponding to those vectors, which would correspond to knowing those vectors priori, defeating the point of using Lanczos all.
Remark. Given the relationship between the moments of the spectral density and the expectation over the set of random vectors, one may be curious to ask how we expect the deviation to be de- pending the number of random vectors actually used in practice. Whilst usually packages using Lanczos use a numher of random vectors (increasing the corresponding computational cost), in Appendix
C we demonstrate that we expect the difference to reduce we increase the problem dimension P. This informs our choice of only using single random vector in our package, we compare results for different random vectors in Appendix C Cnd find minimal differences.
C EFFECT OF VARYING RANDOM VECTORS
Given that the proofs for the moments of Lanczos matching those of the underlying spectral density, are true over the expectation over the set random vectors and in practice we only use a Monte Carlo average of random vectors, or in our experiments using stem plots, single random vector. We justify this with the following Lemma
Lemma 1. Let u E R Px1 random vector, where u, is zero mean and unit variance and finite 4 ''th moment E f = m4. Then for H E TR. PxP, then
i) Hu) Tr H
ii Varlu Hu} S (2 + m4) Tr(H" H)
Proof.
Elu" Hu)
Σ
H Bu, Σ H Tr
(15)
1,j=1
Ef Hul
Elu ugff
k
H.,
H (0 + 8,115, + K5, + madis.al
(16)
(J
k.l
(Tr H (2) Tr(+
Remark. Let us consider the signal to noise ratio for some positive definite > cI
Varlu" Hu)
1
Elu" Hu]
1 P- IAA)
A
(A2)
(17)
1 P-I
x2
where (..) denotes the arithmetic average. For the extreme case of all eigenvalues being identical, the condition number x = 1 and hence this reduces to 1/P - 0 in the P-> oo limit, whereas for a rank-1 1 matrix, this ratio remains 1. For the MP density, which well models neural network spectra, x is not a function of P as P - oo and hence we also expect this benign dimensional scaling to apply.
16

---

Under review as a conference paper at ICLR 2021
We verify this high dimensional result experimentally, by running spectral visualisation but using two different random vectors. We plot the results in Figure [ We find both figures 5a & 5b to be close to visually indistinguishable. There minimal differences in the extremal eigenvalues, with former giving {A1,An) An = {6.8885, 0.0455) } and the latter {6.8891, -0.0456), but the degeneracy O, bulk, triplet of outliers 2.27 and the large outlier at 6.89 is unchanged.
10
10
-
10
10
10
10
10
10
10
10
-10
10
-0.050
2.270
4.580
6.889
10
-0.05
2.27
4.58
6.89
(b) 3
Figure 5: VGG16 6 Epoch 300 end of training Lanezos for different random
C.I WHY WE DON'T KERNEL SMOOTH
Concurrent work, which has also used Lanezos with the Pearlmutter trick to learn the Hessian (Yao et al. 2018 Ghorbani et al. 2019). typically uses n. random vectors and then uses kernel smoothing, give final density. In this section argue that beyond costing factor of n, more computa- tionally, that the extra compute extended order to get more accurate moment estimates, which already argued in Section Carr asymptotically error free, is wasted due to the kernel smoothing (Granziol et al. 2019). The smoothed spectral density takes the form:
((A) = f 1o XX(A=X)(X))'' Σ wo
X₂)
(18)
We make some assumptions regarding the nature of the kernel function, Ais (A-Aa)):) order to prove our main theoretical result about the effect of kernel smoothing on the moments of the underlying spectral density. Both our assumptions are met by (the commonly employed) Gaussian kernel. Assumption 1. The kernel function k. (A >x) is supported the real line [-ox, oo).
Assumption 2. The kernel function k (A - -x)) is symmetric and permits all moments.
Theorem 6. The m-th moment of a Dirac mixture Ei= w.O(A - >d), which is smoothed by kernel kig satisfying assumptions 1 and & 2 is perturbed from its unsmoothed counterpart by an amouns Σi= w. E/ (2) ER, (2, (A25)) 23 where r = m if m is even and m 1 otherwise.
ER, (A)( (X21) denotes the 2j-th central moment of the kernel function ko (A).
Proof. The moments of the Dirac mixture are given as,
Σ w. f 661A-1)2 ((AA-A))A'a
n w.A!".
(X")
(19)
The moments of the modified smooth function (Equation equation 18 are
/
(Am) = Σ
(A)A
1=1
72
Σ
Wi (oo +
(20)
1=1
n
r/2
(1
Σ
wi Σ
(12))
2j
T=1
17

---

Under review as a conference paper at ICLR 2021
We have used the binomial expansion and the fact that the infinite domain is invariant under shift reparametarization and the odd moments of a symmetric distribution are
Remark. The above proves that kernel smoothing alters moment information, and that this process becomes more pronounced for higher moments. Furthermore, given that wi > O, (x)(2) > and (for the GGN lambbda; 0, the corrective term manifestly positive, the smoothed moment estimates are biased.
D LOCAL LOSS LANDSCAPE
The Lanczos algorithm with enforced orthogonality initialised with random vector gives moment matched discrete approximation to the Hessian spectrum. However this information is local to the point in weight space and the quadratic approximation may break down within the near vicinity. To investigate this, we use the loss landscape visualisation function of our package: We display this for the VGG-16 on CIFAR-100 in Figure 8 We for the training loss ball that the eigenvector corresponding to the largest eigenvalue increase in loss for the training, with other extremal eigenvectors, corresponding the eigenvalues = {2.67, }.35] overtaking it in loss change relatively rapidly. Interestingly for the testing loss, all the extremal eigenvectors change the loss much more rapidly, contradicting previous assertions that the test loss is a "shifted" version of the training loss (He et al. 2019) Izmailov et al. (018). We do however note some small asymettry between the changes in loss along the opposite ends of the eigenvectors. The flat directions remain flat locally and some of the eigen-vectors corresponding to negative values correspond to decreases test loss. We include the code in
6.88 only very locally corresponds to the sharpest
= -0.0470
0.5
-0.0452
2.6
-000442
0.4
-0.0348
2.4
-000136
0.3
= -0.0004
2.2
0.0002
0.0059
0.2
= 1.8656
2.0
2.3449
0.1
= 2.6695
1.8
6.8819
0.0
0.2 0.0 0.1 0.2
0.2 0.0 0.1 0.2
Trannig
Testing
Figure 6: VGG-16 1 CIFAR 100 Loss surface visunlised along 6 negative and position eigenvalues
10FAR-10 DATASET
To showcase the ability of our software to handle multiple datasets we display the Hessian of the VGG-16 trained in an identical fashion as its CIFAR-100 counterpart of CIFAR-10 Figure ?? along with the plot of selection of Ritz vectors traversing the training loss surface in Figure Sa and testing loss surface in Figure 8b along with also the training accuracy surface (Figure 9a and testing accuracy surface (Figure 9b).
18

---

Under review as a conference paper at ICLR 2021
100
A = 0.0470
72
À 0.0452
98
À = 0.0442
À = 0.0048
70
96
À 0.0136
À = 0.0004
94
68
À = 0.0002
À 0.0059
92
À 1.8656
66
A = 2.3449
90
A = 2.6695
64
À 6.8819
88
-01 0.1 0.0
0.1 0.2
-0. 0.0 0.1 0.2
(a) Tiinning Accuney
(b) Testing Accurac
Figure VGG-16 CIFAR- 100 Accuracy visualised along negative position eigenvalues
1.4
+ 000203
0.8
-000199
1.2
-0.0194
-0.0169
0.6
1.0
-000055
-0.0007
-0.0000
0.4
0.8
= 0.0016
0.8318
0.6
0.2
=1.0023
1.32265
0.4
1.9310
0.0
-0.2
0.0
0.2
0.2
0.0
0.2
(a) Training Loss
(b) Testing Loss
Figure 8: VGG-16 CIFAR-10 Loss surface visualised along negative and position eigenvalues
E EXAMPLES ON SMALL RANDOM MATRICES
In this section, we use some examples small random matrices to showcase the power of our package that uses the Lanczos algorithm with random vectors to learn the spectral density. Here, we look at known random matrices with elements drawn from specific distributions which converge to known spectral densities in the asymptotic limit. Here we consider Wigner Ensemble (Wigner 1993) and the Marcenko Pastur (Marchenko and Pastur 1967). both of which extensively used in simulations or theoretical analyses of deep neural network spectra (Pennington and Bahri 2017 Choromanska al. 2015a Anonymous 2020
19

---

Under review as a conference paper at ICLR 2021
100
A = 0.0203
À 0.0199
98
À = 0.0194
À = 0.0169
90
À 0.0055
96
À = 0.0007
88
À = 0.0000
94
À 0.0016
86
À = 0.8318
92
À = 1.1023
A 1.3265
84
À = 1.9310
90
82
-0.2
0.0
0.2
-0.2
0.0
0.2
(a) Tiiining Accuney
(b) Testing Accuracy
Figure 9: VGG-16 CIFAR-10 Accuracy visualised along negalive positive eigenvalues
E.I WIGNER MATRICES
Wigner matrices can be defined in Definition E. and their distributions of eigenvalues are governed by the semi-circle distribution law (Theorem
Definition E.1. Let {Y;) and [Zijasis Riisss be two real-valued families of zero mean, i.i. random variables, Furthermore suppose that E2222 1 and for each k E N
max (E) E|V1 <
(21)
Consider a P X P symmetric matrix Mp, whose entries are given by
Mp((i Y; (Mr(i,j) = Zij = Mp(j,),, if 2
(22)
The Matrix Mp is known real symmetric Wigner matrix.
Theorem 7. Let {Mp]P= be sequence of Wigner matrices, and for each P denote Xp Mp/VP. Then pxp, converges weakly, almost surely the semi circle distribution,
a(r)dr
V x= x|52
(23)
2T
For our experiments, we generate random matrices E PxP with elements drawn from the distribution (0,1 1) for P = (225, 10000) and plot histogram of the spectra found by eigende- composition, along with the predicted Wigner density (scaled by factor of /P) in Figures 10b & [0d]and compare them along with the discrete spectral density approximation learned by lanczos in m 30 steps using single random vector d in Figures T0a & 10cl It be seen that for a small number of steps m < P and a single random vector, Lanezos impressively captures not only the support of the eigenvalue spectral density but also its shape. We note as discussed in section
B.2 that the 30 Ritz values here do not span the top 30 eigenvalues even approximately.
E.2
MARCENKO-PASTUR
equally important limiting law for the limiting spectral density of many classes of matrices con- strained to be positive definite, such as covariance matrices, is the Marcenko- Pasturr law (Marchenko and Pastur 1967). Formally, given a matrix X E P with i.i.d zero mean entires with vari-
ance o²
0o. Let \ > X2, 2 Ap be eigenvalues of Yn +xx*. The random measure (P(A)) ( E A). AR
20

---

Under review as a conference paper at ICLR 2021
13
12
Stem
225
(b) Histogram. P 225
30,
10
Histogra,,
10
Figure 10: Lamezos stem plut for a single random vector with rys = 30 steps compared to actual eigenvalue histogram for matrices of the form
R PPx P.X.P. where drawn from distribution with unit variance, converging Wigner circle.
2
3
(a) Stem = 30, P = 10
(b) Histogram. 10
Figure 11: Lanezos stem plot for a single random vector with m. = 30 steps compared to actual eigenvalue histogram for matrices of the form
It Px where X /k, where each element PXA 0.51 is drawn from normal distribution with unil variance, converging Martenko-Pastur distribution q
Theorem 8. Assume that P,N - oc and the ratio P/N + q E (0, 00) (this is known as the Kolmogorov limit) then p. + y in distribution where
(1 ) (A), ifq>1
(24)
(A).),
dvg V(+ - =x)(x X
= o -
(25)
Ax2TO²2
Here, we construct random matrix X P with independently drawn elements from the dis- tribution (00, 1) and then form the matrix 7. X which is known to converge to the Marcenko- 2P and plot the associated histograms from full cigendecomposition in Figures 12b & 12d along with their
Pastur distribution. We use P (225. 10000 and
30, d 1 Lanczos stem coun- terparts in Figures| 12a 12c Similarly we see a faithful capturing not just of the support, but also of the general shape. We note that both for Figure 10 and Figure 12 the smoothness of the discrete spectral density for a single random vector increases significantly, even relative to the histogram.
We also the same experiment for P 10000 but this time with T 0.5P that exactly half of the eigenvalues will be O. We compare the Histogram of the eigenvalues in Figure Tib against its m 30, Lanezos stem plot in Figure Tla and find both the density the origin, along with the bulk and support to be faithfully captured.
E.3 COMPARISON TO DIAGONAL APPROXIMATIONS
As proxy for deep neural network spectra, often the diagonal of the matrix (Bishop) 2006) or the diagonal of surrogate matrix, such as the Fisher information, or that implied by the values of the Adam Optimizer (Chaudhari et al.) 2016) is used. We plot the true eigenvalue estimates for random matrices pertaining to both the Marcenko-Paste blue, along with the Lanczos estimate in red and the diagonal approximation in yellow. We see here that the diagonal approximation in both cases, fails to adequately the support or accurately model the spectral density. whereas the lanezos estimate is nearly indistinguishable from the true binned eigen- spectrum. This couree obvious from the mathematics of the un- normalised Wigner matrix. The diagonal elements simply draws from the normal distribution A00 1) and expect the
(Fig. 14a] and the Wigner density (Fig. [4b) in
21

---

Under review as a conference paper at ICLR 2021
13-
12
0.
2.
Stem
225
(b) Histogram. 225
30,
10
Histogram,
10
Figure 12: Lamezos stem plot For a single random vector with my = 30 steps compared to actual eigenvalue histogram for matrices of the form H E RPPP, where H Xx* /k, where cach cloment of X x = 2.P is drawn from normal distribation with unit variance, converging Marvenko- distribution q
115
27
30
13
80
30
33
20
10
Histogram
Diagonal
Figure 13: Generated matrices H
F 1000 with known eigenspecirum and Lanezos for different values of
(5, 15,30 30)
diagonal histogram plot to approximately follow this distribution (with variance 1). However the second moment of the Wigner Matrix can given by the Frobenius norm identity
-
1 +J==1
P A =
=
(26)
Similarly for the Marcenko-Pard Pastur distribution, We can easily that each element of H follows chi-square distribution of 1/Tx², with mean 1 and variance 2/T.
E.4 SYNTHETIC EXAMPLE
The curvature eigenspectrum of neural network often features a large spike at zero, a right-skewed bulk and some outliers (Sagun al. 2016 2017 In order simulate the spectrum of neural network, we generate a Matrix H E k 1000 with 470 eigenvalues drawn from the uniform distribution from [0, 15). 20 drawn from the uniform [0, 60) and 10 drawn from the uniform [-10, 0). The matrix is rotated through rotation matrix U, i. H = U UDU where D is the diagonal matrix consisting of the eigenvalues and the columns are gaussian random vectors which are orthogonalised using Gram-Schmin orthogonalisation. The resulting eigenspectrum is given in histogram in
Some examples of this can found later sections on real-life neural network experiments see Figures and
15
Manssako-Paatur
Wigner-Semi Clocle
Figure 14: Two randomly generated matrices H E & 500 with the histogram of the eigenvalues in blue, the Lanczos estimate
Tr 30, red and diagonal approximation in yellow
22

---

Under review as a conference paper at ICLR 2021
Figure 13a and then using the same random vector, successive Lanczos stem plots for different
number of iterations m [5, 30] are shown in Figure 13 Figure 13b] for a low number of steps, are the largest and smallest eigenvalues, some information
the degeneracy at À 0 is learned, is retained about the bulk density, but some of the outlier eigenvalues around À Rs 20 and ~ 30 completely missed out, along with all the negative outliers except the largest. For m 30 even the shape of the bulk is accurately represented, as shown in Figure [3d] Here, we would like to emphasise that learning the outliers is important in the neural network context, as they relate to important properties of the network and the optimisation process (Ghorbani et al. 2019).
On the other hand, we note that the diagonal estimate in Figure 13c gives absolutely no spectral information, with no outliers shown (maximal and minimal diagonal elements being 5.3 and 3.3 respectively and it also gets the spectral mass at O wrong. This builds on section F..3] as further- ing the case against making diagonal approximations in general. In neural networks, the diagonal approximation is similar to positing no correlations between the weights. This very harsh as- sumption and usually more reasonable assumption is to posit that the correlations between weights in the same layer are larger than between different layers, leading to block diagonal approximation Martens 2016), however often when the layers have millions of parameters, full diagonal approxi- mations are still used. (Bishop) 2006 Chaudhari et al. 2016).
E.5 COMPARISON TO BACKPACK
For 16-layer VGG network, on the CIFAR-100 dataset the GGN diagonal computations Daggue al. (2119) require over 125GB of GPU memory. Hence we use their Monte Carlo approximation to the spectral visualizations. We plot a histogram of the Monte Carlo approximation of the diagonal GGN (Diag-GGin) against both the Lanczos GGN (Lanc-GGN) and Lanezos Hessian (Lanc-Hess) Figure T3 Note that as the Lane-GGN and Lanc- Hess are displayed stem plots (with the discrete spectral density summing to 1 opposed to the histogram area summing to 1). We note that the Gauss-Newton approximation quite closely resembles its Hessian counterpart, capturing the majority of the bulk and the outlier eigenvectors at A R 6.88 and the triad near A Rs 2.29. The Hessian does still have significant spectral mass on the negative axis, around 37%. However most of this is captured by Ritz value -0.0003, with this removed, the negative spectral mass is only 0.05%. However the Diag-GGN gives a very poor spectral approximation. vastly overestimates the bulk region. which extends well beyond À implied by Lanczos and adds many spurious outliers between 3 and the misses the largest outlier of 6.88. Computational Cost Using a single NVIDIA GeForce GTX 1080 Ti GPU, the Gauss-Newton takes an average 26.5 seconds for each Lanczos iteration with the memory useage 2850Mb. Using the Hessian takes an average of 27.9 seconds for each Lanezos iteration with 2450Mb memory usage.
GGN diagonal against both our CGN-Hessian-LLneeze
10 2
10
10 0
10 -1
10 -2
10 -3
10 -
10
8
10 7
10
10 -8
10 9
-0.050
2.260
4.570
6.882
0.00
2.29
4.59
6.88
(a) Diag: GGN Lane GGN
(b) Diag GGN Hess
Figure 15: Diagonal Generalised Gauss-Newton carlo approximation (Ding-GGN) against 100 Lanezos using Gauss-Newton
vector products (Lanc. (Lamc-GGN) or Hessian vector products (Lance Hess)
23


Reviews:
[{"review_id": "i6G27-nNSCx", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "I think this paper needs to be desk-rejected since the supplementary material reveals the authors' identities in multiple places, e.g. the Readme file and the example notebook. Various .py files have an author attribute set.\n\nThis issue aside, both the package and the paper appear rushed to me and **not ready for publication**.\n\nPackage:\n* Most importantly, I could not find any tests whatsoever. Of course tests do not guarantee anything, especially if implemented poorly, but I hope that nobody would trust a completely untested third-party codebase for their own research. A thorough(!) test suite would be a strict requirement from my view to even consider accepting a code library paper at a top-tier conference.\n* There is no setup.py file to install the package and dependencies. I tried importing some of the modules from the root directory but got errors because of missing dependencies (e.g. backpack is not listed in the requirements.txt).\n* The package is split into multiples modules with generic names such as 'core'. I would strongly suggest to put them in a joint namespace.\n* What's the reasoning for including implementations of standard architectures that are available e.g. through torchvision and high-level training functions? Are those required? In general, I think that for a utility library it is best to assume that users have an established workflow of constructing and training their networks (unless the point of the library is to simplify that workflow, but this does not seem to be the case here).\n\nPaper:\n* I don't think that Section 4 fits into the flow of the paper. I would rather see that space used to go more into depth with the examples. At the moment you assume that the reader already knows why they would want to e.g. estimate the spectral density, but I think it would be worth both motivating the experiments and discussing the results.\n* I think basing the code examples on actual network and data loader objects would make them more accessible. Presumably passing string arguments for model and dataset into your compute_eigenspectrum and train_network functions is limited to the architectures and datasets that you implement? How do users provide their own datasets or networks? Or is this not supported? Section 2 seems to claim otherwise.\n* The font sizes and families around the code listings are extremely inconsistent, especially pages 5 and 6 are visually fairly unappealing.\n* There is some discussion around the relationship to Backpack, but what about PyHessian? That package seems more closely related since it also focuses on computing spectral densities.\n* Comparing first and second order optimisers on a per-epoch basis as in Figure 1 does not seem practically relevant to me, a wall clock time comparison would be of much higher interest to most potential users.\n* I would recommend carefully checking the paper for spelling and language. I'm not a native speaker, but to me it seems like there are quite a few inappropriate capitalisations (\"Training\", \"Testing\"), inconsistent capitalisations (\"eigenvalues\" vs \"Eigenvalues\") and hyphenations (\"eigen-decomposition\").\n* The reference pages are missing (after page 8).", "extra_scores": {}, "extra_text": {}}, {"review_id": "fteq0qfFfOg", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "###\nSummary:\nThis paper proposes a software package to ease and provide a standard way for Hessian-related computation, both for loss analysis and second order optimization. It also provides analysis on why Lanczos algorithm is a better choice to estimate Hessian eigenvalue compared to Power Iterations. Finally, it empirically shows the advantage of using Hessian approximation that goes beyond diagonal approximation for spectrum computation.\n\n###\nReasons for score: \nI lean toward acceptance. Implementing Hessian-related computation is complex and often a bottleneck for research. Easing the implementation of those computations lower the entry barrier for second order optimization/loss analysis research and has the potential to stimulate more works in those areas.\n \n###\n Pros:\n- Address an important implementation problem\n- Scale to medium size network used on CIFAR-10.\n \n###\nCons: \n- Comparisons with related work could be expanded. How does this library related to the other software that use Lanczos algorithm for Hessian computation?\n- How extensive is the library? How easy is to use it with an arbitrary Pytorch model? \n- It is not clear if the library would scale to neural network usually used for larger scale problem such as ImageNet.\n \n###\nQuestions:\n- It is not clear to me how to compute the Hessian of a Neural Network with Batch Norm using minibatch statistics due to the dependency on the other samples of the batch. Could you elaborate on this point?\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "Q4YyW0Lz33Z", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary:\n- This paper introduces a package for computing the second-order information of neural networks based on the Lanczos algorithm. The authors showcase the usages of the package with 1) visualizing the eigenspectrum of the curvature matrix; 2)visualizing the loss surface along with a specific direction, and 3) comparisons between different optimizers. Also, the authors claimed that they address several misconceptions about the Lanczos algorithm.\n\nOverall:\n\nThis paper is well-written and easy to follow. I believe that the developed package can be a good contribution to the machine learning community, especially for people working on second-order optimization and understanding the training dynamics and generalization performance of deep neural networks. However, I still have the following questions:\n\n-  Lanczos algorithm suffers from numerical instability. How do you deal with this issue? Can you provide more details about this?\n\n- I would suggest the authors move Section 4 to the appendix, considering the main scope of the paper is introducing a new machine learning package. In the meantime, the authors should highlight more on the difference between your package and the existing implementations. For example, you can add a table to summarize the features of all the related packages/implementations. Also, a table for comparing the running memory cost and running time cost will be very helpful.\n\nMinors:\n-  as well as the commonly used `Generalsed` Gauss Newton -> as well as the commonly used `Generalised` Gauss-Newton\n\n- The main interface functions are `organsed` as followed: -> The main interface functions are `organized` as following:\n\n- Krylov subspace K (H, v) = span{v,`H^v`,H2v...} is orthogonalised -> Krylov subspace K (H, v) = span{v,`Hv`,H2v...} is orthogonalised\n\n- the reference is missing in the main pdf.\n\n\nRating:\n- This paper did a good job of introducing the package with detailed examples. Also, this package will definitely ease the effort of researchers in related areas to compute the second-order information of the neural networks. In the meantime, I still have the concerns mentioned above. So, my rating is weak acceptance at the current stage. ", "extra_scores": {}, "extra_text": {}}, {"review_id": "65aJb_VQByr", "reviewer": "AnonReviewer5", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "### Summary\nThe paper presents a library for second-order analysis of the optimization of models implemented with PyTorch and potentially having millions of parameters. The library offers tools for easily computing and visualizing Hessian, curvature, loss landscape, and runnning simple statistics, for instance for studying the properties of local minima.\nCompared to existing tools, the proposed library is more complete and accurate, , as demonstrated by example analyses, and scales well with the number of model parameters and the dataset size.\nThe library itself seems very useful, however the paper needs major revision in order to be publishable.\n\n### Presentation\nThe paper is generally poorly written. Many concepts are not defined or not well enough, for instance:\n1. input independence, page 2, by which it is not clear if independence between features or between samples is meant;\n2. number of moments, page 3, not defined;\n3. curvature, by which it should be specified earlier it is meant the loss curvature.\n\nTo improve clarity, the paper should provide a comprehensive list of tools implemented in the library and give examples of why one would need to use each of them. The library's tools are described throughout the paper, but they are too scattered to have a complete view of the library.\nMoreover, the Lanczos algorithm, which is the central for the library, should be reported. \n\nFinally, the paper misses important parts, such as the reference page, a paragraph in page 5 and the cited and not reported Algorithm 1; it also has many typos; notations and domains are not provided for most equations and it is not clear why Equation (3) does not depend on the loss L.\n\nThe current presentation doesn't meet the standards for a conference paper.\n\n### UPDATE\nThe presentation of the paper is now convincing, with the background, contributions and concepts clearly stated. I hence increased my score. \n\nHowever, I agree with reviewer 2 that the library is not properly tested (I understand it can be hard to test the whole computation, but unit tests could be easily provided) and that it would have a higher impact if it were more modular, so that a user could easily add the loss analysis directly in her workflow.\nMoreover, I also think that the paper should be rejected given that the first submission wasn't anonymized and the paper wasn't in a state of being submitted.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a software package to ease and provide a standard way for Hessian-related computation, both for loss analysis and second order optimization. I think all reviewers agree with the usefulness of this work but differ in their assessment whether this work is ready for publication. Given the emphasize on providing a software package I share the view that careful testing and support of usability is important. While there is quite some spread in the scores I think in this case the average score is an appropriate way to compensate for subjective differences between reviewers. I think it is justified to encourage the authors to invest a bit more work to turn this into a fully convincing contribution.


Decision: Reject