Hessian Eigenspectra of More Realistic Nonlinear Models
Authors: Zhenyu Liao
Michael W. Mahoney
Venue: NeurIPS 2021 Oral (2021)
Keywords: Generalized Linear Model
Hessian
Spectral Analysis
Random Matrix Theory.

Abstract: Given an optimization problem, the Hessian matrix and its eigenspectrum can be used in many ways, ranging from designing more efficient second-order algorithms to performing model analysis and regression diagnostics. 
When nonlinear models and non-convex problems are considered, strong simplifying assumptions are often made to make Hessian spectral analysis more tractable.
This leads to the question of how relevant the conclusions of such analyses are for realistic nonlinear models. 
In this paper, we exploit tools from random matrix theory to make a *precise* characterization of the Hessian eigenspectra for a broad family of nonlinear models that extends the classical generalized linear models, without relying on strong simplifying assumptions used previously. 
We show that, depending on the data properties, the nonlinear response model, and the loss function, the Hessian can have *qualitatively* different spectral behaviors: of bounded or unbounded support, with single- or multi-bulk, and with isolated eigenvalues on the left- or right-hand side of the main eigenvalue bulk. 
By focusing on such a simple but nontrivial model, our analysis takes a step forward to unveil the theoretical origin of many visually striking features observed in more realistic machine learning models.

Hessian Eigenspectra of More Realistic Nonlinear Models
Zhenyu Liao
School of Electronic Information and Communications
Huazhong University of Science and Technology, China zhenyu liao@hust. www.huuuutuuuu........omm cn
Michael W. Mahoney ICSI and Department Statistics University of California, Berkeley, USA mmahoney@stat. berkeley. edu
Abstract
Given an optimization problem, the Hessian matrix and its eigenspectrum can be used in many ways, ranging from designing more efficient eccond-order algorithms to performing model analysis and regression diagnostics. When nonlinear models and non-convex problems considered, strong simplifying assumptions made to make Hessian spectral analysis more tractable. This leads to the question of how relevant the conclusions of such analyses for realistic nonlinear models. In this paper, we exploit tools from random matrix theory to make precise characterization of the Hessian eigenspectra for broad family of nonlinear models that extends the classical generalized linear models, without relying on strong simplifying assumptions used previously. We show that, depending on the data properties, the nonlinear response model, and the loss function, the Hessian can have qualitatively different spectral behaviors: of bounded unbounded support, with single- or multi-bulk, and with isolated eigenvalues on the left- or right-hand side of the main eigenvalue bulk. By focusing on such simple but nontrivial model, our analysis takes step forward to unveil the theoretical origin of many visually striking features observed in more realistic machine learning models.
often
Introduction
The Hessian is ubiquitous in applied mathematics, statistics, and machine learning (ML). Given (loss) function L(w) of some parameters w E R", the Hessian H(w) RP*p is defined as the second derivative of the loss function with respect to the model parameter, i.e., H(w) aL (w) (0wōwT) When ML model is being trained, it is common to parameterize that model by w. and then train that model by minimizing some (smooth) loss function L(w), with the associated Hessian ([w). e.g., by backpropagating the error to improve w [26] A Alternativel, once MI. model is trained, the Hessian (and the related Fisher information matrix [66] (6)) can be examined to identify outliers, perform diagnostics, and/or engage in other sorts of model validation [3] 79 62).
For convex problems, the Hessian H(w) provides detailed information on how to adjust the gradient to achieve improved convergence, e.g., in Newton-like methods. For non-convex problems, the properties of the local loss "landscape" around given point w in the parameter space is of central
Work done at ICSI and Department Statistics, University California, Berkeley, USA.
We refer the readers to an extended version of this article 4 for detailed proofs and more discussions.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

39 77, 78 79]. In this case, most obviously, the signs of the smallest and largest Hessian eigenvalue can be used to test whether a given w is a local maximum, local minimum, or saddle point. More subtly, the Hessian eigenvalue distribution characterizes the local curvature of the loss function and provides direct access to, for instance, the fraction of negative Hessian eigenvalues that determines the number of (local) descent directions, quantity that is directly connected to the rates of convergence of various optimization algorithms [3
significance 18 36 13
For theoretical analysis of neural network (NN) models, Hessian eigenspectra are often assumed to follow well-known random matrix distributions such the Marčenko Pastur law [45] or the Wigner's s semicircle law [70]. This enables one to use Random Matrix Theory (RMT), but it involves (for NNs, at least) making relatively strong simplifying assumptions (e.g., the Hessian can be decomposed as the sum of the two freely independent matrices, the residual error, data feature, and weights are all
composed of i.i.d. zero mean normal random variables) [57 58 L5 A somewhat more realistic
K 5 41] In this case, the matrix follows a signal-plus cittnallpus-noise model and consists of full rank random noise matrix and low rank statistical information structure F The "signal" eigenvalues are generally larger than the noisy "bulk" " cigenvalues; and the maximum eigenvalues, when isolated from the bulk, are referred to as the "spikes." A substantial ineory-practice gap exists, however. In both toy examples [27] and ofate--f-tth-arr NN models [77] 78, (79, 80 64 20]. the strong simplifying assumptions are far from satisfactory. (A similar theory-practice gap has been observed for other NN matrices to which RMT has been applied, perhaps most notably weight matrices [47 48 46).) A more precise understanding of the Hessian cigenspectra (and its dependence the input data structure, the underlying response model and model parameters, as well as the loss function) for more realistic models is needed.
setup involves using so-called spiked model (or a spiked covariance model)
1.1 Our approach
In this article, we address these issues, in a setting that is simple enough to be analytically tractable but complex enough to shed light on realistic large-scale models. We consider a family of generalized generalized linear models (G-GLMs) that extends the popular generalized linear model (GLM) 19 31: and we show that, even for such simple models, the key simplifying assumptions used in previous theoretical analyses of Hessian can be very inexact. In particular, apart from few special cases (including linear least squares and logistic regression with homogeneous features), most Hessians of G-GLMs not close to the Marčenko Panttur and/or the semicircle law. Instead, the corresponding Hessian depends on the input feature structure, the underlying response model, and the loss function, more involved fashion that can be precisely characterized by the proposed analysis.
The G-GLM describes generalized linear relation between the input feature x; E Rp and the corresponding response y/+ in the sense that there exists some parameters w. E IP" such that for given w [ x,, the response ys is independently drawn from
y1 f(y w.xx)
(1)
for some conditional density function f({ 1). This extends the classical GLM such as
logistic model: PPyy = w.x) = (1 te w y (-1,
and covers large family of models in applications in statistics and ML. Other examples include: (i) the (noisy) nonlinear factor model where y/ ~ [[g(w].x x?) for some nonlinear IR - R and
a > [4]: (ii) the (noiseless) phase retrieval model with y = (w] [x)x)) in which case one wishes to reconstruct w. from its (squared) magnitude measurements (21:: and (iii) the single-layer NN model
y = af (wIx) x) for some nonlinear activation function (t) such as the tanh-sigmoid (t) = tanh(t).
For given training set ((x;, Jr) of size the standard approach to obtain/recover the parameter w, E R2 of a G-GLM is to solve the following optimization problem
1 Σ ely,, T xx),
min L(w)= min
for some loss function k(y.h) R R R. e.g., the negative log-likelihood of the observation model within the maximum likelihood estimation framework (31 such as the logistic loss ( h)
Since the same informative pattern is repeated in each row column of the matrix.

---

In(1 + e- ) in the case of logistic model 2 In many applications, however, the optimization
may not be convex, for example to achieve superior robustness and/or accuracy and can be NP-hard in general (the noiseless phase retrieval model y = (w, x)² with the square loss Eyy, h) = 2 as an example (9). As we shall see, in such non-convex G-GL.Ms, the dominant Hessian eigenvector can be shown, in some cases, to positively correlate with the sought-for
problem in 3
49
71 I0
parameter w, and therefore be used as the initialization of gradient descent methods 9 37 34]. This
particularly motivates our study of the possible isolated Hessian eigenvalue-eigenvecoor pairs.
1.2 Our main contributions
The main contribution of this work is the exacr characterization of Hessian eigenspectra for the family of G-GLMs, in the high-dimensional regime where the feature dimension p and the sample size n are both large and comparable. Precisely, we establish:
1. the limiting eigenvalue distribution of the Hessian matrix (Theorem I]: and 2. the behavior of (possible) isolated sigenvalue-eigenveeee pairs (Theorem and 3).
as function of the dimension ratio limp/n. n, feature statisties, the loss function E in and the underlying response model in (T. Our results are based on a technical result of independent interest:
3. deterministic equivalent (Theorem 4 of the random resolvent Q(z) (H x1p) for
2 E C not an eigenvalue of H, of the generalized sample covariance
H = (w) = 2/4 Ei= E" (y): w xxxxxxx LXDXT
(4)
for X = -x......, E Ixxn D diag (" w w }}} E RTX= and (" = ELy h)/Oh, as n. + with p/n - E (0, o0), under the setting of generic Gaussian feature x; ~ W(H, C), for μ E IR" and positive definite covariance C E RPXP, We also demonstrate results empirically by showing that:
4. for a given response model (T), the Hessian eigenvalue distribution depends on the choice of loss function and the data/feature statistics in an intrinsic manner, e.g., bounded versus unbounded support and single- versus multi-bulk in Fig and
5. there may exist two qualitatively different spikes one due to data signal μ and the other due w. or and thus the underlying model-which may appear on different sides of the main bulk, and their associated phase transition behaviors are characterized (Fig f versus 5
To have more clear picture of our contribution, we compare, in Fig Taland Tb the Hessian eigenvalues for the logistic model (2) with the logistic loss Eyy h) = In(1 + e Wh) (which is the loss of choice within the maximum likelihood framework), for different choices of in the parameter space. A nontrivial interplay between the response model, feature statistics and the parameter w is reflected by the range of the Hessian eigenvalue support and an additional right-hand spike in Fig 1b confirmed by our theory. For phase retrieval model y = w_'x²² with square loss C(y, h) = non-convex nature of the problem is reflected by (relatively large) fraction of negative Hessian eigenvalues in Fig Ic We also note that the top eigenvector (that corresponds to the largest eigenvalue) contains structural information of the underlying model, in the sense that it is positively correlated with w., as predicted by our theory. This is indeed connected to the Hessian- based initialization 5T).
h222// /4, the
scheme widely used in non-convex problems S 37
We conclude by emphasizing that, by focusing on the simple yet fundamental G-GLM, we obtain results that improve upon and are different than previous efforts in the following aspects:
(i) We provide precise asymptotic characterizations of the Hessian eigenspectra that go beyond, e.g., 77 where only Hessian lower bounds are given in the case of logistic model with logistic loss: our methodology and theoretical results hold much more generally for the family of G-GLM with arbitrary loss. As illustrating examples, we discuss linear least squares in See 3.1 3.1
logistic model with different choices of loss function in Fig phase retrieval model in Figure Tc and more in 43 Sec 4].
(ii) We extend the results in 57 53
50 to G-GLMs by considering generic data statistics
and loss function, whereas in previous efforts only much more homogeneous models are
'Here we consider w independent of x., Vs. With some additional effort, our results extend to the case where takes some simple explicit function form of X and y. e.g., w is the least square solution (X, y)
3

---

discussed, and sometimes under unrealistic assumption, e.g., the Hessian can be decomposed as the sum of two freely independent matrices, and the residual error, data feature, and weights are all composed of i.idd. zero mean normal random variables 57 58
(iii) Instead of focusing solely on the main eigenvalue bulk as in (57 58 our results also shed
novel light on the isolated eigenvalues (above and/or below the bulk) that are empirically
observed in the Hessian of modern NNs 60
22
2
53
52
well
the associated
cigenvectors that are shown closely connected to NN training dynamics (28). Also, relative
[57 [58) we show qualitatively different behaviors for the Hessian eigenspectra, e.g..
bounded versus unbounded support, single- versus multi-bulk as in Figure 2 knowledge, these are not covered in the existing Hessian literature.
To our
0.2
Eigenvalues
6
Fingiricd
Tup sigenvastior
4
Theory.
Theary
0.1
x
2
2
0.1 0.2 0.3 0.4
0.2
0.4
(a) Logistic, w w = (b) Logistic, indep. w
(c) Phase retrieval model (d) Eigenvector in Fig
Figure 1: Illustration of our main results: eigenspectral properties of the Hessian of G-GLMs with y 800, n 6 600 and C Ip. Fig la versus lb; absence versus presence of a right-hand side spike for different choices of w. logistic model (2) with logistic loss, and w. = fe ~ (O,.IIpp)) Ip/p). Fig versus versus 1d: the Hessian eigenspectra have rather different shape (as opposed to the Marčenko- Pastur-like in Fig la and [bb for the (non-convex) phase retrieval model (1c) and the top eigenvector is known in this case to be (noisy) estimate of w. (1d), as confirmed by theory. With square loss {v. h) = (y- h h2)²/4. W. = (-2 1p/2 2- 1 1/22/ VP. w ~ N(O, Ip/p) and μ = 0.
1.3 Related work
Here, we provide a brief review of related previous efforts, see more discussions in 43
Random matrix theory. Random matrices of the type 4 are related to the separable covariance model [7 in the RMT literature, which is of the form ZDZ C for Z random and C. D deterministic random but independent of Z. Our results generalize this, by allowing D to depend on Z, in possibly nonlinear fashion, per (4). This is of direct interest for the Hessian of GsGLMs.
Hessian eigenspectra. The cigenspectra of Hessian matrices arising in ML models (in particular, for NNs) have attracted considerable interest recently [60 [1] 22 72 24 32 20 64
79 80 25).
However, these investigations
often built upon somewhat unrealistic simplifying assumptions and reduce to the "mixed" behavior of Marčenko- -Paturr and semicirele law (57 [5] In contrast, here we focus on the more tractable example of G-GLM and provide precise results on the Hessian eigenspectra for structural feature on arbitrary loss. Efforts have also been devoted to efficiently computing the eigenspectra of Hessian matrices, which, despite not enabling analytic results, allows for the computation or approximation of Hessian spectra for large-scale ML models such as DNNs,
without any assumption the data or underlying (response) model 30
79
74 25
Spectral initialization in non-convex problems. A popular initialization scheme (of gradient- based methods) for non-convex problems is the spectral initialization, where the top eigenvectors of some
Hessian-type matrices are used as gradient descent initialization [S] 37 34 2 STI In 44, which was generalized in [50), the authors evaluated the eigenspectrum asymptotics of 1 Σ f(y, )x for x; ~ (00, (0, Ip). Their technical approach is, however, limited to the case of very homogeneous features with x; ~ (V(O,,)) I,). Here we generalize the analysis in [4 50 to the Hessian of G-GLM, by developing a systematic approach to account for both feature structures and loss functions.
Scalable second-order methods. Second-order methods are among the most powerful optimization methods that have been designed, and there have been several attempts to use their many advantages
for machine learning applications 775 69 59], particularly for training NNs [79 80 64 20 73
We
4

---

expect that our precise characterization of the Hessian sheds new light on the understanding and improved design of (e.g., computationally) more efficient second-order methods.
Spectra of realistic NN weight matrices. Recent work demonstrated that the eigenvalue distribution of weight matrices of state-of-the-art DNN models in computer vision and natural language applica- tions exhibit heavy-tailed (i.e., not Marčenko-Pastur behavior, and that this can be used to define metrics to predict trends in the quality of state-of-the-ar NNs without access to training testing It would be of interest to extend our approach to weight matrix analysis.
data 47 48 46..
2
Main results
In the section, we present our main results: on the limiting Hessian eigenspectrum (in 2ec and on the behavior of the (possible) isolated eigenvalue-eigen eigenvector(s) (in Sec 2.2). These two main results depend on a technical deterministic equivalent result for the Hessian resolvent (in Sec 2.3). which is of independent interest. We position ourselves in the following high high-dimensional regime.
Assumption 1 (High-dimensional asymptotics), As n, + co with p/n + c E (0, 00), we have max [|w]]] |, 1I3 = 0(1) and x,
i.i.d. Ni,, C) with max [| [[[|]] = O(1)
2.1 Limiting spectral measure
Our first result is the limiting Hessian eigenvalue distribution. This is a direct consequence of our main technical Theorem 4 and proven in [43 Sec A.2].
Theorem 1 (Limiting spectral measure). Let Assumption 7 hold, we have, as nt, - Dc with p/n - e. E (0, oo), the empirical spectral measure hh of the Hessian matrix H defined in (4) converges weakly and almost surely a probability measure p. defined through its Stieltjes transform m(z)= = [tt - 2) =u(dt) as the unique solution to
1 m(2)= = tr Qs(2), 6(2) tr (CQ+(=)), Q (x) p
1
9.C
=E
=Ipp
(5)
1 1 5(2)
where
forh= N(w H. wo Cw)), and y and defined respectively in ① and B}. Moreover, we denote v the law/distribution of and assume the empirical spectral measure of C converges as p - oc, then writes
(6)
m(z) / (-2++t (-2tt/ 1+ v(dt)) B(d)), 8(2)= / -z+tf [+10[=] v(d)) (df).
ct
(7)
In the form of (7), the (Stieltjes transform of the) limiting Hessian spectral measure is determined by the ratio c = Tim p/n and the two measures v and D. This formulation is closely connected the separable covariance model [38 6
76 in RMT. Moreover, if (ddt) = 81 (t) is Dirac mass at one, this reduces to the popular sample covariance model 655;; taking further (dt) 8, (t) gives the ManrennooPasttrr law. See See 3.1 for numerical evaluations of these special cases. In particular, the support of the (limiting) Hessian spectrum j is directly linked to that of v and i.
55
17
Remark (Hessian eigen-support). Under Assumption ! the (limiting) spectral measure of C has bounded support. However, this may not be the case for the law of g defined in (6). Since the Hessian eigenvalue distribution μ of compact support if and only if both v and have compact support NI7 Porposition 3.4], may be of unbounded support, depending on the model and the loss.
An example of unbounded μ the phase retrieval model with y = xw] and square loss f(y. h) (y- - 122²/4 /4, for which we have g 3(w Tx)2 (w]x)² for x .Ni,, C). As consequence, with say w. = w, g follows a chi-square distribution with one degree of freedom and has thus unbounded support. This corresponds to Fig [c. where the Hessian spectrum has a "heavier" tail compared to Fig La (logistic model), and the empirically observed "isolated" eigenvalue is due to
'That is, the normalized counting measure of the eigenvalues ).(H) of H. i.e., HN = = E ŌX((H)) Unnqqueeess is ensured in such way that 3[m())]- 3[z] > 0 for 3[2] o 0 and 2m(2) < 0 for 3[z] = O. so that (z, m())) is valid Stieltjes transform couple, more details in 29
5

---

Eigenvalues
0.6
4
0.4
2
0.2
0.1 0.2 0.3 0.4
2
4
0.2 0.4 0.6
0.5
1
(a) Logistic loss
(b) Exponential loss
D (6, + 62)
(d) 2 (8 84)
Figure 2: Impact of loss function: bounded 2a with logistic loss) versus unbounded 2b with exponential loss) Hessian eigenvalues, with 800, 6 000, logistic model in (2) with μ 0.
I= p. w. 0 and w = 1-1 p/2. 1//2// vp. Impact of feature covariance: Hessian spectrum of single- 2c with C diag p/23 2- 1//2 versus multi-bulk (2d with C = [liag[1p/2 4 1/22),
with 500, n = 6 000, logistic model with w
0p.. w = (O,,Ip/p))
finite-dimensional effect and will be "buried" in the noisy main bulk for larger values of n. p. Therefore, aiming for an (almost surely) isolated cigenvalue-cigenvecoo (e.g., to recover the model parameter w. using the top Hessian eigenvector), some pre-processing function f must be applied. This has been discussed in previous work 44 50 and corresponds to the so-called trimming strategy in phase retrieval 12, with for instance the truncation function (() d for some >>
Another example of unbounded μ is when the exponential loss 23 is used. Precisely, consider the logistic model (2) with Ely. h) exp(-yf) -yh), we have that y exp( yh) for h ~ N(w' fr, w Cw) which follows a log-normal distribution and has unbounded support. As such, the (limiting) Hessian eigenvalue distribution y has also unbounded support. On the other hand, with logistic loss ((y,)) one has g < 1/4 and μ is guaranteed to have bounded support. In Fig 2a and 2b the
In(1 + v wh), empirical Hessian eigenvalues and the limiting distributions are compared for logistic and exponential losses, with more "heavy-tailed" behavior observed for the exponential loss.
Clearly, depending on the measures (of g. which depends on feature statistics, loss and underlying model) and ū (of spectrum of feature covariance C), the Hessian spectrum can have very different forms. Here we compare the empirical Hessian eigenvalues with their limiting behaviors per Theo- rem | for different feature covariance structures In particular, with more "compact" Hessian spectrum as in Fig 2e multiple bulks (two the case of Fig 2d) with Hessian eigenvalues more "spread-out" depending
may observe single main bulk
the feature covariance structure i. In the form of 77). the condition for the existence of multi-bulk eigenspectrum has been thoroughly discussed in 17 See 3.2-3.4) 3.4] and can be numerically evaluated with ease.
As a side remark, the "multi-bulk behavior similar to Fig 2d has been empirically observed in Hessians of NNs in [42] 53 and is believed to be due to the classification structure within data (i.e., the data vectors are drawn from a mixture of distributions). Here, we provide an alternative explanation via feature covariance structure that holds beyond the classification setting.
2.2 Isolated eigenvalues and eigenvectors
As discussed in Remark under Assumption T. the Hessian has bounded (limiting) eigen-support if and only if v, the law of g. has bounded support. Under this condition (or, after the application of some function ff so that f(g) is bounded), we can then talk about the (possible) isolated Hessian eigenvalues, as in the following result, the proof of which is given in [43] Sec A.3].
Theorem 2 (Isolated eigenvalues). In the setting of Theorem assume that the law v of the random variable g defined in (6) of bounded support, define
G(2) = A(2)V Qo(:VV R3X3
with (=))))(()) defined in (5) V= ju, Cw., Cw] E RPX3
U C² (w.: w] E Rex2 and
A = F
g
U Utz(Utz)--U'T))))
2) -UTU)
z z= c c- (x- ~ 99)
1
1+
8(2)
"Covariance describes the joint variability or the "correlation" between entries of the feature vector and it is of particular significance in the analysis of image (with local structure) and time series data.

---

where we denote U the Moore-Penrose pseudoinverse of U. Then, for X such that G()) has zero eigenvalue (of multiplicity one), there exists an eigenvalue À of H such that X - X = 0.
Theorem 2 provides asymptotic characterization of the possible isolated Hessian eigenvalues by computing the determinant of the much smaller (three-by-three) deterministic matrix G closely related to the key quantity ((2) defined in Theorem Note that, Theorem 2 does not provide, at least explicitly, the phase transition condition under which these spikes become "isolated" from the main bulk. As we shall see in more details in See 3.2 two types of quantitatively different phase transitions can be characterized, due to the data "signal" and the underlying model, respectively.
We can also analyze the associated isolated eigenvectors. First note that, in the infinite data regime (i.e., for - Do with p fixed), we have, by the strong law of large numbers, that H 0.8, TE
[[H, with
Elg U+
E H E[[" (y, w x)xx B[g] C+V V
U+Elg (zz" I(U VT.
Eig Utz]
As consequence, is expected that in the large n. p oo limit, the top eigenvectors of H could also be related to the columns of V. This is the case in Fig Id where the top cigenvector is observed to be "noisy" version of the model parameter w.. More precisely, for (Å, ú) an isolated eigenpair of H. the projection V'uuV E R333 can be shown to be asympiotically close to a deterministic matrix. This measures the crosie-ssmillatee between the Hessian isolated eigenvector with any column of V and consequently the performance of using ú as an estimate of, for instance the model parameter w. for C C=Ip.
This result is given in the following theorem, which is proven in 43 Sec A.4].
Theorem 3 (Isolated eigenvectors). In the setling of (TThorrmm 2 for isolated eigenvalue-eigeenvcioo
pair (A, ú) of H and À the asymptotic position (of X) given in Theorem 2 then
Vūū V = (AVV E(A) + o(1), E(A) -v
for (2) and G(z) defined in (5) and (8) respectively, VI,G. V.G RR' the left and right eigenvec- tors of (A) associated with eigenvalue zero, and (A)) the derivative of G(z) with respect to evaluated at z X
2.3 Technical tool: deterministic equivalent
Our main technical tool to derive Theorem and is so-called deterministic equivalent 29 16 result for the Hessian resolvent Q(z) = (H
that provides simultaneous access to the Hessian limiting eigenvalue distribution and the behavior of the possible isolated eigenpairs. Precisely, the normalized trace tr Q (<2)/p gives the Stieltjes transform m=(2)) = f(t- " "A of the empirical
spectral measure μии of H, from which ичи can be recovered. Also, for (", ú) eigenpair of interest, 1 with Cauchy's 's integral formula we have |w = 2m fr, w (=())w dz, for deterministic
vector w E IR and I'A a positively oriented contour surrounding only Å. As such, for Q(2) deterministic equivalent of Q(z). that is, Q(z) Q(z) with tr AQ(=) - Q(=))/p ~ and
T (Q(=) Q(=))b 0 almost surely as n. + oc, for A RPxP and a, b Rp of bounded (Euclidean and spectral) norms, the limiting spectral measure (via the associated Stieltjes transform) and the isolated eigenpairs of H are directly accessible via the study of the deterministic equivalent (2).. This result is given as follows, with the proof deferred to [43 Sec A.1].
Theorem 4 (Deterministic equivalent). Let Q(=) = (H =Ip) -1 be the resolvent of H defined in Then, under Assumption I as n. + oa with p/n > c E (0, oo),
Q(z) + Q(z), with (=) E |
C (I, Pu) + +xx =Ipp
for random vector a = f C E RR and = C" w 44
TH + w (²z) for z ~ IO,II,) as defined in O y and 6(2) defined in ① and (5). respectively. and Pu RRXP the projection onto the subspace spanned by the columns of U- C [ w).
3 Evaluations and Discussions
In this section, we provide further discussions on the consequences of Theorem T2 and together with numerical evaluations. Implications of Theorem | on the Hessian eigenvalue distribution is
7

---

Eigenvalues
Rescaled MP
Distribution of
Empinical
Eigenvalues
Distribution of
Seale facter
0.2
0.4
0.1
6
10
(a) Logistic: Hessian
(b) Logistic: g:s
(c) Exp: Hessian
(d) Exp: ges
Figure 3: Empirical Hessian eigenspectra (blue, with limiting behavior in red per Theorem versus
rescaled and shifted Marčenko-Pastiee laws (green) in the setting of Fig with 800 and n 6 000. Fig[3a] versus [bb Marčenko-PPatur--Hn empirical distribution of the 9:8 versus the scaling factor (green in 3bb obtained by empirically
Hessian with logistic loss, Hessian eigenvalues (3a and
matching the minimal and maximal empirical Hessian eigenvalues
the Marčenko-Passe law).
Fig versus versus 3d] an example of Hon-Marrennoo-Paae Pasioo---a-urllliieeeeeeee Pake Hessian with exponential loss and the associated g:s (green). Note that the scales of the axes are different in different subfigures.
discussed in Sec 3.1 In Sec 3.2 discuss the consequences of Theorem 2 and on the possible isolated eigenpairs, for which two fundamentally different phase transitions are characterized.
3.1 Hessian eigenvalues distribution
For better interpretation of Theorem I on the Hessian eigenspectrum, consider here the special case of C = Ip. and start with the simple setting where the random variable in (6) is constant, say = 1: this happens, e.g., when the square loss (y, h) = (y - h)²/2 is employed. In this case, the Hessian does not dependent w, w. and the Stieltjes transform m(2) is the solution to =cm?(2) - (1 (1-C c 2))(2) ) +
and corresponds to the Marëenko-Pass law.
As long is not constant, the limiting Hessian spectrum is. priori, different from the Marčenko- Pastur law, even in the C I, setting, since the associated Stieltjes transform m(z) is different from the solution the Marčenko-Pastie
equation. However, we see in Fig 3a that, for the logistic model (2) with logistic loss, the Hessian spectrum is close, at least visually, to a (rescaled) Marčenko-Pastue law. This be understood with Theorem D and is due to the fact that, the distribution of g is more "concentrated" (around some constant, see Fig 3b versus 3d versus33dfo This is in sharp contrast to Fig 3c where with the exponential loss, the law of has much larger spread and the Hessian is therefore away from Marčenko- Mardenko-Pastur--haaee Pastur-shape.
a comparison between different cases).
This "empirical fit" has been observed in [56] Fig 5], where acceleration methods proposed for Marčenko-Pasiee Pastur distributed Hessian (in linear least squares) work reasonably well on logistic regression models. Our theory proposes a convincing theoretical explanation of this empirical observation on logistic regression, and possibly for others more involved MI. models. Nonetheless, must be pointed out that this "visual approximation" by Marčenko-Pastur law is not robust, in the sense that it "visually" holds only for, yet formally different from, the case of (i) logistic model with (ii) logistic loss and (iii) identity covariance C Ip: any change in the response model (c.g.. the phase retrieval model in Fig Ic in the choice of loss function (c.g.. the exponential loss in Fig 2b)., or beyond the identity covariance setting (as in Fig 2 and 2d) would induce Hessian spectrum that is very different from the Marčenko- Pastur law. this vein, our Theorem I goes beyond such "loose" Marčenko Parttur approximation and acts of Hessian in more involved ML models beyond linear least squares that accounts for nonlinear transformations (such activation function in NNs) and feature statistics.
more accurate first example in the understanding
While Theorem Tis is proven here only for Gaussian features, we conjecture, as the case for many random matrix asymptotics, that it holds more generally beyond Gaussian distribution, see 43 Fig 5] and more discussions therein.
3.2 Isolated eigenvalues and their phase transitions
In this section, we discuss the implications of Theorem and 3 on the possible isolated eigenvalue- eigenvector pairs. More precisely, we show that, different from the classical spiked models extensively studied in RMT literature 3 41,, for which (i) the isolated spike appears due to the presence

---

of some statistical "signal" in the data and (ii) a "monotonic" phase transition behavior can be characterized function the signal strength; here another type Hessian spike arises due the underlying G-GLM model (i.e., w. and w) and exhibits a rather different behavior.
0.2
Eigenvatues
Empirical
Eigen align
Eigenvalues
4
Limitng
Eigengap
0.15
0.1
Thooretical
0.5
Right odge
0.05
2
0 044aA44
04444
ve
1.5
Vc
1.5
0.2 0.4
0.6
||
||
Figure 4: Spike due to data signal in Corollary 1 theory versus practice of (left) Hessian eigen- spectrum with |||| 0.8, (middle) eigengap (... and (right) top eigenvector alignment a in (TO). as a function of the signal strength ||]]| on logistic model with logistic loss, for x 1p/p p/2.. W. 0, C Ip. 512 and 048. Results averaged over 50 runs.
10
Eigenvalues
8
Eigensalues
Eigengap
0.03
Empirical
Eigen-aligg align
0.8
Limijing
0.02
Theomtical
0.6
6
edge
0.4
0.01
0.2
2
2
4
6
8
0
||[]
||w]
0.1
0.2
0.3
Figure 5: Left-hand side spike due to response model in Corollary the absence of data sig- nal: (left) Hessian spectrum for ||w|| || = 2, with a left isolated eigenvalue dw: (middle) eigengap dist supp(u)), and (right) dominant eigenvector alignment (with w), function ||w|| with w x 1/2: 1p/2]. w. = μ 0, C Ip. = 800 and n 000. Results averaged over 50 runs.
3.2.1 Spike due to data signal
To study the spike due to data "signal" μ and its phase transition behavior, we focus here on the case w, w = O. This, in the case of logistic model (2) for example, gives rise to much simpler form of limiting spectrum (per Theorem I and possible isolated eigenpairs (per Theorem and 3). as
summarized in the following corollary, with detailed derivations given in 43 Sec B.3].
Corollary 1 (Spike due data signal: logistic model). Consider the logistic model in 2 with logistic loss, for w w. = 0 and C = 1p the limiting Hessian eigenvalue distribution is the Marčenko-Pastae law, but rescaled by factor of g = 1/4. Moreover, there is one isolated
eigenpair (Nys ún)) of H and it satisfies
{
A +11+p+c p+l) 111 2o)2
p
p> Vo,
JA F ||
3 --
p> Vc, p P Vc:
tep
PS Vc:
(10)
with the signal strength so |||l|2 and = lim p/n.
The behavior of the isolated
described in Corollary follows the "classical" phase transition 3 [4]: (i) the isolated eigenvalue always appears on the right-hand side of the main (Marčenko-Paseee bulk and (ii) the eigenvalue amplitude and eigenvector alignment is "monotonic" with respect to the signal strength ||a|2 the sense that, for fixed dimension ratio c, the largest Hessian eigenvalue is bound to become asympiotically isolated once ||da||² exceeds ve and its value, as well as the eigenvector alignment, increase monotonically as || grows. This is confirmed in Fig 4
As we shall see below, this is not the case for, e.g., the spike due to model parameter w.

---

3.2.2 Spike due to model
To investigate the spike due to the underlying model (i.e., w, and w), we position ourselves in the situation where μ O, that is, in the absence of data "signal". This leads the following corollary, the proof of which is given in [43] Sec B.4].
Corollary 2 (Spike due to model: logistic model). Consider the logistic model in ② with logistic loss, P 0 and C Ig. then the Stieltjes transform m(2) satisfies m(z) (E[(((r,))
2) for [T; =) = (/((((2) + 2 t + e") and NN(,, [[w])²) that depends on w but not on W.- Moreover, there is at most one isolated eigenvalue Row of H that is due to and satisfies
Aw Aw 0.5 O with Aw solution to 0 = det. G()w) =1 (tw)) B[f(r Yr² [w 12))
[[
The situation here is more subtle (than the spike due to data signal discussed in Sec 3.2. : as the model parameter w changes (e.g., as the "energy" ||w/| grows). both the Hessian (limiting) eigenvalue distribution and the possible spike location are impacted. Fig iillsstratss the behavior of the spike due to w in the setting of Corollary 2 Note first that, different from the case of spike due data signal μ. the spike in Fig (leff)t appears on the left-hand side of the main bulk: this particularly means that the Hessian may admit an eigenvalue that is significantly smaller than all the other eigenvalues Also, note from Fig 5ig (middle) that, different from the spike due to fs, the spike due to w exhibits here "non-monotonic" behavior the sense that, is absent for small values of |/| (as for small || in Fig 4 middle) and becomes "isolated" as ||w|| increases, but then again "merges into" the main bulk as w continues to increase, resulting an eigengap that falls back to zero.
It is perhaps even more surprising to observe in Fig (rightt that, the alignment between the associated isolated eigenvector and the parameter w is, unlike the eigengap in (middle), monotonically increasing as || grows large, as in the case of Fig (right). This suggests that, in the case of spike due to model, a smaller eigengap may not always imply less statistical "information" contained in the associated cigenvector, which somehow goes against the conventional eigengap heuristic [67] 35]. worthy mentioning that, while, technically speaking, the proposed analysis not capable of characterizing the behavior as ||w|| + оo under Assumption [ empirical results suggest that for extremely large ||w]], the eigengap tends to "vanish", the associated dominant eigenvector can still be
used to recover almost perfectly, see [43 Fig 8] as an example.
4 Conclusion
In this article, we provided precise asymptotic characterization of Hessian eigenspectra (including the limiting eigenvalue distribution and the behavior of possible isolated cigenvalue-eiienvvctor pairs) for the family of G-GLMs via a random matrix analysis. By assuming generic Gaussian input features, our results show that the Hessian can have a "heavy-tailed" or multi-bulk spectral behavior, and with isolated eigenvalues of the left- or right-hand side of the main bulk, depending on the underlying GLM and the choice of loss function.
Extending our approach to more involved ML models such as multi-layer neural networks chal- lenging, e.g., one will need to first study the layer-wise wise Hessian and evaluate how the objective of interest (e.g., the eigenvalue distribution) of the full Hessian depends on each of its (aser-wise) blocks. Simplification can be made by specifying the choice of loss function and/or using different approximations (such as NeuusssNettottyy approximation [63) or working in the infinitely but less-realistic NTK regime [32] )) so to facilitate analysis. We leave this for future work.
Acknowledgments and Disclosure of Funding
ZL would like to acknowledge the Fundamental Research Funds for the Central Universities of China (2021XXJSI 10) and CCF-Hikvision Open Fund (20210008) for providing partial support. MWM would like to acknowledge DARPA, IARPA (contract W91 INF20C0033) NSF, and ONR via its BRC on RandNLA for providing partial support. Our conclusions do not necessarily reflect the position or the policy of our sponsors, and no official endorsement should be inferred.
'Depending on the response model and loss, the spike due to model may also appear on the right-hand side of the bulk or even establish a left- teft-to-right transition. The eigenvector alignment can behave differently from Fig 4 S and establish a monooootoiic behavior as function of ||w]] or ||w=|]: see [43] Sec 4.3J.
10

---

References
[1] Ryan P Adams, Jeftrey Pennington, Matthew Johnson, Jamie Smith, Yaniv Ovadia, Brian Patton, and James Saunderson. Estimating the Spectral Density Large Implicit Matrices. arXiv, 2018.
Sanjeev Arora, Rong Ge, Tengyu Ma, and Ankur Moitra. Simple, efficient, and neural algorithms for sparse coding, Journal of Machine Learning Research, 40(2015), 2015.
[3] Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. The Annals of Probabiliiy, 3365):1643-1699 2005.
[4] Jinho Baik and Jack W Silverstein. Eigenvalues of large sample covariance matrices of spiked population models. Journal of multivariate analysis, 97(6): 1382-1408 2006.
[5] Florent Benaych-Georges and Raj Rao Nadakuditi. The eigenvalues and cigenvectors of finite, low rank perturbations of large random matrices. Advances in Mathematics, 227(1):494-4 2011.
[6] Zdzislaw Burda, Jerzy Jurkiewicz, and Bartlomiej Wactaw. Spectral moments of correlated wishart matrices, Physical Review E, 2,77112 71(2):026111 2005.
[7] Dankmar Böhning. Multinomial logistic regression algorithm. Annals of the Institute of Statistical Mathematics, 44(1):197-200, 1992.
[8] Emmanuel J Candes, Thomas Strohmer, and Vladislav Voroninski. Phaselift:
Exact and stable signal recovery from magnitude measurements via convex programming. Communications on Pure and Applied Mathematics, 66(8):1241-1274, 2013.
Emmanuel J. Candès, Xiaodong Li, and Mahdi Soltanolkotabi. Phase Retrieval via Wirtinger Flow: Theory and Algorithms. IEEE Transactions on Information Theary, 19855-77772 2015.
[10] Olivier Chapelle, Choon Teo, Quoe Le, and Alex Smola. Tighter bounds for structured estimation. Advances in neural information processing systems, 21:281-288, 2008
[11] Pratik Chaudhari, Anna Choromanska, Stefano Soatto, Yann LeCun, Carlo Baldassi, Christian Borgs, Jennifer Chayes, Levent Sagun, and Riccardo Zecchina. Entropy-SOO SGD: Biasing gradient descent into wide valleys. Journal of Statistical Mechanics: Theory and Experiment, (01(12))124018 2019.
[12] Yuxin Chen and Emmanuel Candes. Solving random quadratic systems of equations is nearly as casy as solving linear systems. In Advances in Neural Information Processing Systems, pages 739-747, 2015.
[13] Yuejie Chi, Yue M. Lu, and Yuxin Chen. Nonconvex Optimization Meets Low-Rank Matrix Factorization: An Overview. IEEE Transactions on Signal Processing, 67(20):5239 5269, 2018.
[14] Dennis Child. The essentials of factor analysis. Cassell Educational, 1990.
[15] 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 292-204, 2015.
[16] Romain Couillet and Merouane Debbah. Random matrix methods for wireless communications. Cambridge University Press, 1.
[17] Romain Couillet and Walid Hachem. Analysis of the limiting spectral measure of large random matrices of the separable covariance type. Random Matrices: Theory and Applications, 3(04): 1450016, 2014.
[18] Yann N Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho. Surya Ganguli, and Yoshua Bengio. Identifying and attacking the saddle point problem in high-dimensional non-convex optimization. In Advances in Neural Information Processing Systems, volume 27 of NIPS' 14, pages 2933-222222 29411 Curran Associates, Inc., 2014,
[19] Annette J Dobson and Adrian G Barnett. An introduction to generalized linear models. CRC press, 2018.
[20] Zhen Dong, Zhewei Yao, Amir Gholami, Michael W Mahoney, and Kurt Keutzer. HAWQ: Hessian aware quantization of neural networks with mixed -xexeddpreciioon.. In Proceedings of the IEEE International Conference on Computer Vision, pages 293 003-302,22 2019.
[21] James R. Fienup. Phase retrieval algorithms: comparison. Applied Optics, 21(15):2758, 1982.
[22] Stanislay Fort and Surya Ganguli. Emergent properties of the local geometry of neural loss landscapes. arXiv preprint 1rXiv.1910.05929. 20199

---

[23] Yoav Freund, Robert Schapire, and Naoki Abe. A short introduction to boosting. Journal- Suciety For Artificial Intelligence, 147111780) 1999.
[24] Mario Geiger, Stefano Spigler, Stéphane d' A'Asolli, Levent Sagun, Marco Baity-Jesi, Giulio Biroli, and Matthieu Wyart. Jamming transition paradigm understand the loss landscape of deep neural networks. Physical Review E, 100(11)012115, 2019.
[25]
Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao. An Investigation into Neural Net Optimization via Hessian Eigenvalue Density. 97:2232-2241, 2019.
[26] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016.
[27] Diego Granziol. Beyond random matrix theory for deep networks. arXiv preprint WrXiv:200000721,, 2020.
[28] Guy Gur-Ari, Daniel A. Roberts, and Ethan Dyer. Gradient descent happens in tiny subspace, 2018.
[29] Walid Hachem, Philippe Loubaton, and Jamal Najim. Deterministic equivalents for certain functionals of large random matrices. The Annals of Applied Probability, 17(3):75-9330, 2007.
[30] Insu Han, Dmitry Malioutov, Haim Avron, and Jinwoo Shin. Approximating Spectral Sums of Large- Scale Matrices using Stochastic Chebyshev Approximations. SIAM Journal on Scientific Computing. 1558))A 2017.
[31] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction. Springer Science & Business Media, 2009.
[32] Arthur Jacot, Franck Gabriel, and Clement Hongler. The asymptotic spectrum of the hessian of dnn throughout training. In International Conference on Learning Representations, 2019.
[33] Prateek Jain and Purushottam Kar. Non-convex optimization for machine learning. Foundations and Trends Machine Learning, 10(3-4):12-3366, 2017.
[34] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi. Low-rank matrix completion using alternating minimization. In Proceedings of the forty-fifth annual ACM symposium Theory of computing, pages 665-674, 2013.
[35] Antony Joseph and Bin Yu. Impact of regularization on spectral clustering. The Annals of Statistics. (765):776-1177 2016.
[36] Kenji Kawaguchi. Deep Learning without Poor Local Minima. In Advances in Neural Information Processing Systems, volume 29 of NIPS 16, pages 586-594. Curran Associates, Inc., 2016.
[37] Raghunandan II Keshavan, Andrea Montanari, and Sewoong Oh. Matrix completion from a few entries. IEEE transactions on information theory, 2666))2980-2998, 2010.
[38] A Khorunzhy, Eigenvalue distribution of large random matrices with correlated entries. Mat. Fiz, Anal. Geom, 3(1-2):80-101 101, 1996.
[39] Jason D. Lee, loannis Panageas, Georgios Piliouras, Max Simchowitz, Michael I. Jordan, and Benjamin Recht. First-order methods almost always avoid strict saddle points. Mathematical Programmingg 176(1- 2):311-3377 2019.
[40] Kiryung Lee, Yanjun Li, Marius Junge, and Yoram Bresler. Blind recovery of sparse signals from subsampled convolution. IEEE Transactions on Information Theory, 63(2):802-821, 2016.
[41] Thibault Lesieur, Florent Krzakala, and Lenka Zdeborová. Phase transitions in sparse pea. In 2015 IEEE International Symposium on Information Theory (ISIT). pages 1635-1639. IEEE, 2015.
[42] Xinyan Li, Qilong Gu, Yingxue Zhou, Tiancong Chen, and Arindam Banerjee. Hessian based analysis of sSGD for deep nets: Dynamics and generalization. In Proceedings of the 2020 SIAM International Conference on Data Mining, pages 190-198. STAM, 2020.
[43] Zhenyu Liao and Michael W. Mahoney. Hessian eigenspectra of more realistic nonlinear models. arXiv preprint arXiv:2 100.0:200000111 2021.
[44] Yue M Lu and Gen Li. Phase transitions of spectral initialization for high-dimensional non-convex estimation. Information and Inference: A Journal of the IMA, 2019.
[45] Vladimir Alexandrovich Marchenko and Leonid Andreevich Pastur. Distribution of eigenvalues for some sets of random matrices. Matematicheskii Sbornik, 114(4):007-33,, 1967.
12

---

[46] C. H. Martin and M. W. Mahoney. Post-mortem on deep learning contest: a Simpson's 's paradox and the complementary roles of scale metrics versus shape metrics. Technical Report Preprint: arXiv:2 100.04444. 2021.
[47] Charles Martin and Michael W. Mahoney, Implicit self regulaaiaaaiinnnnn deep neural networks: Evidence from random matrix theory and implications for learning. Journal of Machine Learning Research, 22(165):1-73, 2021.
[48] Charles H. Martin, Tongsu Serena Peng. and Michael W. Mahoney. Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data. Nature Communications,, 12(1):1-133 13, 2021.
[49] Llew Mason, Jonathan Baxter, Peter I. Bartlett, and Marcus R Frean. Boosting algorithms gradient descent. In Advances neural information processing systems, pages 512-518, 2000.
[50] Marco Mondelli and Andrea Montanari. Fundamental Limits of Weak Recovery with Applications Phase Retrieval. Foundations of Computational Mathematics, 19(3):003-733, 2019.
[51] Marco Mondelli, Christos Thrampoulidis, and Ramji Venkataramanan. Optimal Combination of Linear and Spectral Estimators for Generalized Linear Models. Foundations of Computational Mathematics, pages 1-54, 2021.
[52] Vardan Papyan. Measurements of three-level hierarchical structure in the outliers in the spectrum of deepnet hessians. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5012 10221, PMLR, 09-15 Jun 2019.
[53] Vardan Papyan. Traces of class/cross-class structure pervade deep learning spectra. Journal of Machine Learning Research, ( 2020.
[54] Debashis Paul. Asymptotics sample eigenstructure for a large dimensional spiked covariance model. Statistica Sinica, pages 1617-1642, 2007.
[55] Debashis Paul and Jack W Silverstein. No eigenvalues outside the support of the limiting empirical spectral distribution of a separable covariance matrix. Journal of Multivariate Analysis, (00(1):37-55) 57, 2009.
[56] Fabian Pedregosa and Damien Scieur. Average-case acceleration through spectral density estimation. arXiv preprint arXi::2002.047555 2020.
[57] 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- -oolumee 70, pages 2798 2806. JMLR, org, 2017.
[58] Jeffrey Pennington and Pratik Worah. The spectrum of the fisher information matrix of single-hidden-llyer neural network. Advances in Neural Information Processing Systems, 31:5410-5419, 2018.
[59] Farbod Roosta- Khorasani and Michael W. Mahoney. Sub-sampled Newton methods. Mathematical Programming, 1 (74(1-2):23,-326, 2019.
[60] Levent Sagun, Leon Bottou, and Yann LeCun. Eigenvalues of the hessian in deep learning: Singularity and beyond. arXiv preprint arXiv:1611 1.0776.. 2016.
[61] 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: (eee: 04454, 2017.
[62] Adepu Ravi Sankar, Yash Khasbage, Rahul Vigneswaran, and Vineeth N Balasubramanian. A deeper look at the hessian eigenspectrum of deep neural networks and its applications regularization. arXiv preprint arXiv:2012.0380) 2020.
[63] Nicol N. Schraudolph. Fast Curvature Matrix-Vector Products for Second-Order Gradient Descent, Neural Computation, 14(7):173317788, 2002.
[64] Sheng Shen, Zhen Dong, Jiayu Ye, Linjian Ma, Zhewei Yao, Amir Gholami, Michael Mahoney, and Kurt Keutzer. Q-BERT: Hessian Based Ultra Low Precision Quantization of BERT. In AAAI, pages 8815-8821, 2020.
[65] Jack W Silverstein and ZD Bai. On the empirical distribution of eigenvalues of a class of large dimensional random matrices. Journal of Mudtivariate analysis, 54(2):175 192, 1995.
13

---

[66] Aad W. Van der Vaart. Asymptotic Statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2000.
[67]
Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and computing, 17(4):395-416, 2007.
[68]
Martin J. Wainwright. High-Dimension Statisticx: A Non-Asympiotic Viewpoint. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press, 2019.
[69] Shusen Wang, Farhod Konsta-Khorassnn Peng Xu, and Michael W. Mahoney. GIANT: Globally Improved Approximate Newton Method for Distributed Optimization. In Advances in Neural Information Processing Systems, volume 31, pages 2332-2342. Curran Associates, Inc., 2018.
[70] Eugene P. Wigner. Characteristic vectors of bordered matrices with infinite dimensions. Annals of Mathematics, 62(3)5548-564, 1955.
[71] Yichao Wu and Yufeng Liu. Robust truncated hinge loss support vector machines. Journal of the American Statistical Association, 102(479):77--98, 2007.
[72] Yikai Wu, Xingyu Zhu, Chenwei Wu, Annie Wang, and Rong Ge. Dissecting bessian: Understanding common structure of hessian in neural networks. arXiv preprint arXiv:2010.04261, 2020.
[73] Peng Xu and Michael W. Mahoney. Second-order Optimization for Non-convex Machine Learning: an Empirical Study. In Proceedings of the 2020 SIAM International Conference on Data Mining, pages 199-207, 2020.
[74] Peng Xu, Farbod Roosta- Khora--aho and Michael W. Mahoney. Newton-type methods for non-convex optimization under inexact Hessian information. Mathematical Programming. 184(1-2):35-70, 2020.
[75] Peng Xu, Jiyan Yang, Farbod Roosta- Moista-Khorrsami, Christopher Ré, and Michael W. Mahoney. Sub-sampled Newton Methods with Non-uniform Sampling. Advances Neural Information Processing Systems, volume 29, pages 3000-3008. Curran Associates, Inc., 2016.
[76] Fan Yang. Edge universality of separable covariance matrices. Electronic Journal of Probability, 24, 2019.
[77] 7. Yan, A. Gholami, Q. Lei, K. Keutzer, and M. W. Mahoney. Hessian-based analysis of large hatch training and robustness to adversarics. Technical report, 2018. Preprint: arXiv: 180.::888000211
[78] Z. Yao, A. Gholami, P. Xu, Keutzer, and M. Mahoney, Trust region based adversarial attack on neural networks. Technical report, 2018. Preprint: arXiv:1812,06377..
[79] Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael W Mahoney. PyHessian: Neural networks through the lens of the Hessian. arXiv preprint arXiv: 19t.199922....... 2019.
880 Zhewei Yao, Amir Gholami, Sheng Shen, Mustafa Mustafa, Kurt Keutzer, and W. Michael Mahoney. Adahessian: An adaptive second order optimizer for machine learning. Proceedings of the AAAI Conference on Artificial Intelligence, (5(12):16655-000)) May 2021.
[81] Lixin Zhang. Spectral analysis of large dimensional random matrices. National University of Singapore PHD Thesis, 2006.
14


Reviews:
[{"review_id": "sx-8JBweQtG", "reviewer": "Reviewer_3BGa", "summary": "This paper studies generalized linear models in the large dimension, large sample size limit, obtains an analytic expression for the Hessian of the loss function, and discusses many qualitative features such as the relation to random matrix results and the appearance of isolated eigenvalues.", "questions": "", "limitations": "", "rating": 9, "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": "This paper studies a very broad class of generalized linear models including the standard ones, phase retrieval, and essentially any model with a response distribution depending on $w_* \\cdot x$ and a loss function of the form $ \\sum_{i\\in\\mbox{samples}} L(y_i,w\\cdot x_i) $.  An analytic expression for the Hessian is obtained for the large dimension p and large sample size n limit with fixed ratio p/n.  One assumes that the data distribution is approximately normal in this limit.  The derivation is an explicit computation, more or less decomposing the Hessian along the  $w_*,w$ and orthogonal directions and thus obtaining a spiked covariance model.  The result has a similar analytic form to those obtained using random matrix theory and reduces to them in special cases, but takes into account the full nonlinearity of the loss function and the resulting spectra can be quite different.\nComparisons are made to studies of the Hessian for particular models: logistic, limits of neural networks, many empirical studies.\nThe qualitative features seen there, in particular isolated eigenvalues and \"multi-bulk\" spectra, are reproduced.\nPrevious observations that empirical results are not fit by RMT results such as Marcenko-Pastur are explained. \nA new isolated eigenvalue (or \"spike\") whose location is a function of the model parameter w is sometimes present.\n\nThis looks like quite a general and powerful result and I am enthusiastic about the paper.  Although similar techniques are used in many recent papers, so far as I know the result in the generality stated here is new.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No apparent societal impact.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.5"}, "extra_text": {}}, {"review_id": "_7jtAhHU6m9", "reviewer": "Reviewer_rM2S", "summary": "This paper studies the spectrum of empirical cost functions used with data from generalized generalized linear models (G-GLM). Such models depend on some parameters $w^* \\in \\mathbb{R}^p$ that link the input features $x_i \\in \\mathbb{R}^p$ to the responses $y_i \\in \\mathbb{R}$ by the relation $y_i \\sim  f(y |  (w^*)^T x_i)$\nwhere $f(\\cdot | \\cdot)$ is some conditional density function.\nGiven a training set $(x_i, y_i)$ for $i=1, \\dots, n$ a popular approach to estimate $w_*$ is to minimize the cost function\n$\nL(w) = \\sum_{i=1}^n \\ell(y_i, x_i^T w)\n$\nfor some loss $\\ell$ that is carefully chosen.\n\nThis framework is quite general and encompass many classical models such as the standard linear model, the logistic model or the phase retrieval problem.\n\n\nThe authors studies this model in the setting where\n- $n,p \\to \\infty$, while $n/p$ converges to some positive number.\n- $x_i \\sim_{iid} \\mathcal{N}(\\mu, C)$\n- $\\|w\\|, \\|w_*\\|, \\| \\mu \\|, \\|C\\| = O(1)$\n\n\nUnder these assumptions, they provide a precise assymptotic caracterization of the eigenvalues and eigenvectors of the Hessian\n$\nH = \\frac{\\partial^2 L}{\\partial w \\partial w}(w).\n$\nMore precisely:\n\n1. They show that the empirical distribution of the eigenvalues of $H$ converges weakly almost surely to a probability measure $\\mu$ defined by an equation on its Stieltjes transform.\n2. When $\\mu$ has a bounded support, they caracterize precisely the assymptotic postitions of potential isolated eigenvalues. When such isolated eigenvalues exists, the authors also compute the limiting angle between the correspondeing eigenvectors and the vectors $w,w_*,\\mu$.\n3. The authors illustrate their results on popular models and losses, showing the existence of interesting phenomena. In particular, they show that isolated eigenvalues can not only appears due to 'signal' in the data (large $\\|\\mu\\|$) but also appear because of the model itself.\n\n\nThe main results are obtained through an so-called 'deterministic equivalent', that allows to study precisely the resolvant of the Hessian $H$.", "questions": "", "limitations": "", "rating": 9, "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": "I found this paper very interesting and very relevant to a machine learning audience.\n\nObtaining a theoretial understanding of (non-convex) loss functions is an outstanding challenge that has recieved a lot of attention over the recent years. This paper provides a complete and precise analysis of the spectrum of the Hessian of G-GLM allowing to understand better the optimization of such models.\nEven if this model remains simple it is very relevant because very popular and more realistic than the ones studied in the literature.\nIt would be interesting to see if one could extend this analysis to multilayer models, but one would probably have difficulties to obtain exploitable analytic formulas.\n\nHence I think that the approach taken by the authors is completely relevant. The results they obtain are a clear improvement of the ones present in the literature. It is interesting to see that some assumptions made in earlier papers are perhaps too simplistic (always predict a Marcenko Pastur eigenvalue density) for more realistic models such as G-GLM.\n\nThe results obtained in this paper are strong.\nThe paper is very well written, with many examples and comments that allow to understand the subtelties of the results.\n\nThe results of this paper can also be of practical interest, since one could use the top eigenvector of $H$ as initialization for iterative algorithm. On that point, I would be interested to know if one could achieve a correlation with $w_*$ by minimizing $H$ for some carefully chosen $\\ell$ that is better than the spectral estimators proposed in [46].\n\nI would be also curious to know if one could potentially obtain results about the Hessian at point $w$ that depend on the data $x_i, y_i$. One could for instance study the spectrum of $H$ for $w$ for which $L$ is small, but that may be very difficult...\n\n\nTo sum up, I think that this is a paper with strong results that answer and raises many interesting questions.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2h30"}, "extra_text": {}}, {"review_id": "RJ-hN6Cs48f", "reviewer": "Reviewer_MLx1", "summary": "This papers develops new results in random matrix theory applied to characterize the limiting spectral distribution of Hessian matrices of common loss functions.  Specifically,  the setting is what the authors call 'generalized GLMs',  i.e. a class of models that extends classical generalized linear models.   The paper addresses a shortcoming in the literature where much of the analysis of Hessian matrices uses fairly strong simplifying assumptions (e.g. such as Marcenko-Pastur).  Even though the paper doesn't consider deep NNs,  even within the context of G-GLMs it finds rich predictable spectral structure not yet accounted for in the literature, and should have important consequences.", "questions": "", "limitations": "", "rating": 8, "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": "The main contribution of the paper is a new mathematical result describing the empirical spectral measure of Hessian matrices for G-GLMs.  The result employs random matrix theory (RMT) tools: Stieltjes transforms, deterministic equivalents of random resolvents, free probability, e.t.c.), and extends existing results for so-called 'separable covariance models' to capture G-GLMs in full generality.  The RMT line of attack has been considered in earlier papers e.g. Benaych-Georges and Nadakuduti [4],  and Pennington and Worah [52]. However, the simplifying assumptions used in these papers to make the analysis more mathematically tractable lose important detail of the model.  This paper goes through the derivations not making these assumptions.  Similar to Marcenko-Pastur it provides remarkably accurate predictions of  empirical spectral distributions, and is able to describe important salient features of such spectra (existence of isolated eigen-values,  the shape of the bulk of the distribution,  whether the distribution is a single-bulk or multi-bulk (multi-modal),  and whether it has finite or infinite support).  \n\nI believe that the paper offers the ML community a new powerful tool which would be valuable in getting new more accurate insight into both G-GLMs,  and (eventually) deep NNs, and perhaps help pave the way for more efficient second-order SGD algorithms.  The paper is well written, and shows important implications of their results (from initialization of non-convex models, to inadequacy of currently popular assumptions on Hessian analysis).  I would have liked to see some discussion / hints on how to extend these results to NNs where they'd be of most interest,  and some comments on large-scale approximate spectral computations.  Furthermore, since RMT tools may be out of reach for many ML researchers, I would recommend to put together a toolbox that can predict/estimate Hessian spectra for moderate-size models. \n\nComments: \n1) The term 'generalized generalized linear model'  sounds a bit awkward (I thought it's a typo at first).  I don't have great suggestions here -- but since it's a natural generalization of GLMs -- it's worth thinking about a better name.  E.g. Extended GLMs,  or wide-sense GLMs...   \n\n2) Since you extend the RMT analysis from papers such as [52] which ultimately aim to analyze deep NNs,  it may be worth summarizing a pathway of how to use your results for multi-layer NNs.\n\n3) It's worth adding more context for computational tools from randomized-linear algebra that can be used to compute approximate eigen-spectra of Hessian matrices for very large models.  You cite [72] among many citations (there are also related jax tools), but provide little context -- e.g. while they do not enable analytic results, they are indeed able to compute accurate Hessian spectral densities for large deep networks without any assumptions and can be used for analysis of specific models. \n* Adams, Ryan P., et al. \"Estimating the spectral density of large implicit matrices.\" arXiv preprint arXiv:1802.03451 (2018).\n* Han, Insu, et al. \"Approximating spectral sums of large-scale matrices using stochastic chebyshev approximations.\" SIAM Journal on Scientific Computing 39.4 (2017): A1558-A1585.\n\n4) As you're working in the high-dimensional setting (n,p large), Hessian spectra can be thought of as being 'smoothed-out' due to finite-sample noise.  Could there be value in applying 'deconvolution' tools as proposed in\nEl Karoui, Noureddine. \"Spectrum estimation for large dimensional covariance matrices using random matrix theory.\" The Annals of Statistics 36.6 (2008): 2757-2790.\nto the G-GLM setting? \n\n5) Minor:  you use loss and objective in the first few sentences -- which is a bit confusing (e.g. perhaps objective could involve regularization terms in addition to loss).  Worth using consistent language. \n\n6) Are there any success stories in second-order SGD in practice?  It's worth mentioning if this is still an active research area, or if there are indeed practical existing algorithms.  E.g. LBFGS is not exactly a Hessian-based method. \n\n7) page 1: typo faction -> fraction. \n\n8) page 3:  \"(maximum likelihood) logistic loss\".  Unclear -- do you mean Hessian at the ML-estimate of the logistic loss?  Then it doesn't make sense to compute it for 'different choices of w in the parameter space'.  Rewrite to make it clearer. \n\n9) It's worth adding a reference / definition for \"freely independent matrices\".  Also not clear what you mean by 'very homogeneous features'.  Also the term 'the law of g' can be unclear.\n\n10) For random Z and C, D independent of Z.  It's unclear what's random -- if D is not random -- then what does it mean to be 'independent' of Z? \n\n11) Is the setting for figure 3 again p=800, n=6000?  The meaning of the green bar in (d) is confusing. \n\n12)  Page 9:  are isolated eigenvalues above the bulk and below the bulk equally interesting? ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A. ", "review_text": "", "extra_scores": {"ethical_concerns": "N/A. ", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "2P37NVFrb5", "reviewer": "Reviewer_d8k3", "summary": "The paper explores the spectrum of Hessians for a class of generalized generalized linear models which includes convex problems such as logistic loss and non-convex problems such as the phase retrieval. The authors make some simplifying assumptions such as data $x_i \\sim N(\\mu, C)$ so that random matrix theory tool can be applied. The result is a complete characterization of the spectrum of the Hessian via the Stieltjes transform which shows how the covariance in the data, the loss function, and the distribution of the feature vector each effect the density of the eigenvalues. In particular, they show that the Hessian $H(w)$ is deterministically equivalent to a generalized sample covariance matrix $H(w) = 1/n X D X^T$ where $D$ is a diagonal matrix with the second derivatives of the loss evaluated at the current $w$. The authors go on to analyze various properties of $H$ for logistic loss and phase retrieval in detail and show that for many cases the Hessian does not resemble a pure Marchenko-Pastur. \n\nSome other contributions:\n\n(1). They show that the spectrum of the Hessian $H$ is unbounded due to either (i) the ESD of the covariance in the data converging to an unbounded measure (this is known even for the classes spiked models) and (ii) if $g := \\frac{\\partial^2 \\ell}{\\partial h^2} |_{x^T w} \\sim \\nu$ is unbounded. In the case of phase retrieval, the $g$ follows a chi-distribution and therefore unbounded; this accounts for the heavy-tailed features in the spectrum. \n\n(2). The authors show that there are two different reasons for an outlier eigenvalue to occur. The first is that the data itself is not normalized to be 0 mean and so the $\\mu$ creates a spike (usually to the left of the bulk). The second outlier occurs because of the parameters $w$. ", "questions": "", "limitations": "", "rating": 8, "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": "The paper extends the models for Hessians in novel and practical way to generalized generalized linear models (e.g. nonconvex objective functions such as phase retrieval) and adds covariance structure. The paper is well-written with some simulated data to show how their theoretical results match the empirical spectral distributions. Although I did not check all the details for the mathematics, I am confident that results are correct.\n\nI have a few clarifying questions for the authors, but otherwise I am happy to accept this paper. It has some interesting random matrix theory results for Hessians. \n\n(1). If the largest eigenvalue converges to the outlier, does the second largest eigenvalue converge to the edge of the bulk\n(2). Could the authors clarify on the assumptions on $w$ versus $w_*$? For instance, is $w^*$ or $w$ chosen independent of $X$? More importantly how are the targets $y$ chosen and does this choice effect anything? Does adding noise to the targets effect anything?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The work is purely theoretical. The authors made clear their assumptions (e.g. the data matrices distributed as Gaussian and the class of objective functions they are considering). ", "review_text": "", "extra_scores": {"ethical_concerns": "There are no ethical concerns; the result is purely theoretical. ", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper develops new mathematical results in random matrix theory that characterize the limiting spectral distribution of Hessian matrices corresponding to so-called generalized generalized linear models (G-GLMs) in the high-dimensional regime. Using the technique of deterministic equivalents, the authors are able to relax a number of the simplifying distributional assumptions that prior work has used to pursue analyses of this type. The result is an asymptotically exact characterization of the limiting spectrum that exhibits the types of structures that often show up in practice, including isolated outlier eigenvalues, multimodal densities, etc. As such, this paper offers the NeurIPS community a new and powerful perspective for understanding the spectra of high-dimensional models, and potentially paving the way for new developments in second-order optimization and beyond.

The phrase in the title of the paper, "more realistic," actually provides a good characterization of the scope and generality of the results: the focus here does not actually touch on realistic models, NNs, etc., but the G-GLM setup and the general distributional assumptions move the needle significantly in the "realistic" direction. While this paper does observe some phenomena that do tend to appear in practical configurations (such as large outliers), it remains unclear whether the explanations offered here are in fact the same ones that underlie the phenomena more generally. An expanded version of this paper would surely benefit from some detailed (empirical) analysis of this question. Still, absent these additions, the techniques and insights derived in this paper will nevertheless be of interest to the community, and this paper will make a great addition to NeurIPS.

Decision: Accept (Oral)