Hessian Eigenspectra of More Realistic Nonlinear Models
Year: 2021

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 singleor multi-bulk, and with isolated eigenvalues on the leftor 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.

heading: 1 Introduction
text: The Hessian is ubiquitous in applied mathematics, statistics, and machine learning (ML). Given a (loss) functionL(w) of some parameters w ∈ Rp, 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) = ∂L(w)/(∂w∂wT). When a 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 H(w), e.g., by backpropagating the error to improve w [26]. Alternatively, once a ML model is trained, the Hessian (and the related Fisher information matrix [66, 68]) can be examined to identify outliers, perform diagnostics, and/or engage in other sorts of model validation [31, 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 a given point w in the parameter space is of central ∗Work done at ICSI and Department of Statistics, University of California, Berkeley, USA. †We refer the readers to an extended version of this article [43] for detailed proofs and more discussions.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
significance [18, 36, 13, 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 a 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, a quantity that is directly connected to the rates of convergence of various optimization algorithms [33].
For theoretical analysis of neural network (NN) models, Hessian eigenspectra are often assumed to follow well-known random matrix distributions such the Marc̆enko–Pastur law [45] or the Wigner’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, 15]. A somewhat more realistic setup involves using a so-called spiked model (or a spiked covariance model) [3, 5, 41]. In this case, the matrix follows a signal-plus-noise model and consists of full rank random noise matrix and low rank statistical information structure.3 The “signal” eigenvalues are generally larger than the noisy “bulk” eigenvalues; and the maximum eigenvalues, when isolated from the bulk, are referred to as the “spikes.” A substantial theory-practice gap exists, however. In both toy examples [27] and state-of-the-art 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 eigenspectra (and its dependence on the input data structure, the underlying response model and model parameters, as well as the loss function) for more realistic models is needed.

heading: 1.1 Our approach
text: 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 a few special cases (including linear least squares and logistic regression with homogeneous features), most Hessians of G-GLMs are not close to the Marc̆enko–Pastur and/or the semicircle law. Instead, the corresponding Hessian depends on the input feature structure, the underlying response model, and the loss function, in a more involved fashion that can be precisely characterized by the proposed analysis.
The G-GLM describes a generalized linear relation between the input feature xi ∈ Rp and the corresponding response yi, in the sense that there exists some parameters w∗ ∈ Rp such that for given wT∗xi, the response yi is independently drawn from
yi ∼ f(y | wT∗xi) (1)
for some conditional density function f(· | ·). This extends the classical GLM such as
logistic model: P(y = 1 | wT∗x) = (1 + e−w T ∗x)−1, y ∈ {−1, 1}, (2)
and covers a large family of models in applications in statistics and ML. Other examples include: (i) the (noisy) nonlinear factor model where y ∼ N (g(wT∗x), σ2) for some nonlinear g : R→ R and σ > 0 [14]; (ii) the (noiseless) phase retrieval model with y = (wT∗x)
2, in which case one wishes to reconstruct w∗ from its (squared) magnitude measurements [21]; and (iii) the single-layer NN model y = σ(wT∗x) for some nonlinear activation function σ(t) such as the tanh-sigmoid σ(t) = tanh(t).
For a given training set {(xi, yi)}ni=1 of size n, the standard approach to obtain/recover the parameter w∗ ∈ Rp of a G-GLM is to solve the following optimization problem
min w L(w) = min w
1
n n∑ i=1 `(yi,w Txi), (3)
for some loss function `(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 `(y, h) =
3Since the same informative pattern is repeated in each row or column of the matrix.
ln(1 + e−yh) in the case of logistic model (2). In many applications, however, the optimization problem in (3) may not be convex, for example to achieve superior robustness and/or accuracy [49, 71, 10], and can be NP-hard in general (the noiseless phase retrieval model y = (wT∗x)
2 with the square loss `(y, h) = (y−h2)2 as an example [9]). As we shall see, in such non-convex G-GLMs, the dominant Hessian eigenvector can be shown, in some cases, to positively correlate with the sought-for 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-eigenvector pairs.

heading: 1.2 Our main contributions
text: The main contribution of this work is the exact 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 1); and 2. the behavior of (possible) isolated eigenvalue-eigenvector pairs (Theorem 2 and 3),
as a function of the dimension ratio c = lim p/n, feature statistics, the loss function ` in (3), and the underlying response model in (1). Our results are based on a technical result of independent interest:
3. a deterministic equivalent (Theorem 4) of the random resolvent Q(z) = (H− zIp)−1, for z ∈ C not an eigenvalue of H, of the generalized sample covariance:4
H ≡ H(w) = 1n ∑n i=1 ` ′′(yi,w Txi)xix T i ≡ 1nXDX T, (4)
for X = [x1, . . . ,xn] ∈ Rp×n, D ≡ diag{`′′(yi,wTxi)}ni=1 ∈ Rn×n, and `′′(y, h) ≡ ∂2`(y, h)/∂h2, as n, p → ∞ with p/n → c ∈ (0,∞), under the setting of generic Gaussian feature xi ∼ N (µ,C), for µ ∈ Rp and positive definite covariance C ∈ Rp×p. We also demonstrate our results empirically by showing that:
4. for a given response model (1), 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 2; and
5. there may exist two qualitatively different spikes—one due to data signal µ and the other due to w∗ or w and thus the underlying model—which may appear on different sides of the main bulk, and their associated phase transition behaviors are characterized (Fig 4 versus 5).
To have a more clear picture of our contribution, we compare, in Fig 1a and 1b, the Hessian eigenvalues for the logistic model (2) with the logistic loss `(y, h) = ln(1 + e−yh) (which is the loss of choice within the maximum likelihood framework), for different choices of w 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, as confirmed by our theory. For phase retrieval model y = (wT∗x)
2 with square loss `(y, h) = (y − h2)2/4, the non-convex nature of the problem is reflected by a (relatively large) fraction of negative Hessian eigenvalues in Fig 1c. 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 scheme widely used in non-convex problems [8, 37, 34, 40, 2, 51].
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., [7], 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 Sec 3.1, logistic model with different choices of loss function in Fig 2, phase retrieval model in Figure 1c, and more in [43, Sec 4].
(ii) We extend the results in [57, 53, 44, 50] to G-GLMs by considering generic data statistics and loss function, whereas in previous efforts only much more homogeneous models are
4Here we consider w independent of xi, yi. With some additional effort, our results extend to the case where w takes some simple explicit function form of X and y, e.g., w is the least square solution to (X,y).
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.i.d. 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, 42, 53, 52], as well as on the associated eigenvectors that are shown closely connected to NN training dynamics [28]. Also, relative to [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. To our knowledge, these are not covered in the existing Hessian literature.
0.1 0.2 0.3 0.4 0
2
4
6
E ig
en va
lu es
Empirical
Theory
(a) Logistic, w = w∗ = µ
0.2 0.4 0
2
4 0.2 0.4 0
2
4
(b) Logistic, indep. w
−10 0 0
0.1
0.2
−10 0 0
0.1
0.2
−10 0 0
0.1
0.2
(c) Phase retrieval model
Top eigenvector
Theory
(d) Eigenvector in Fig 1c
Figure 1: Illustration of our main results: eigenspectral properties of the Hessian of G-GLMs with p = 800, n = 6 000 and C = Ip. Fig 1a versus 1b: absence versus presence of a right-hand side spike for different choices of w, logistic model (2) with logistic loss, and w∗ = µ ∼ N (0, Ip/p). Fig 1c versus 1d: the Hessian eigenspectra have a rather different shape (as opposed to the Marc̆enkoPastur-like in Fig 1a and 1b) for the (non-convex) phase retrieval model (1c) and the top eigenvector is known in this case to be a (noisy) estimate of w∗ (1d), as confirmed by our theory. With square loss `(y, h) = (y − h2)2/4, w∗ = [−2 · 1p/2; 2 · 1p/2]/ √ p, w ∼ N (0, Ip/p) and µ = 0.

heading: 1.3 Related work
text: 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 [81, 17] in the RMT literature, which is of the form C 1 2 ZDZTC 1 2 , for Z random and C,D deterministic or random but independent of Z. Our results generalize this, by allowing D to depend on Z, in a possibly nonlinear fashion, per (4). This is of direct interest for the Hessian of G-GLMs.
Hessian eigenspectra. The eigenspectra of Hessian matrices arising in ML models (in particular, for NNs) have attracted considerable interest recently [60, 61, 11, 22, 72, 24, 32, 20, 64, 79, 80, 25]. However, these investigations are often built upon somewhat unrealistic simplifying assumptions and reduce to the “mixed” behavior of Marc̆enko–Pastur and semicircle law [57, 15]. 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 on the data or underlying (response) model [30, 1, 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 [8, 37, 34, 40, 2, 51]. In [44], which was generalized in [50], the authors evaluated the eigenspectrum asymptotics of 1n ∑n i=1 f(yi)xix T i , for xi ∼ N (0, Ip). Their technical approach is, however, limited to the case of very homogeneous features with xi ∼ N (0, Ip). Here we generalize the analysis in [44, 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 [75, 69, 59], particularly for training NNs [79, 80, 64, 20, 73]. We
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 applications exhibit heavy-tailed (i.e., not Marc̆enko–Pastur) behavior, and that this can be used to define metrics to predict trends in the quality of state-of-the-art NNs without access to training or testing data [47, 48, 46]. It would be of interest to extend our approach to weight matrix analysis.

heading: 2 Main results
text: In the section, we present our main results: on the limiting Hessian eigenspectrum (in Sec 2.1); and on the behavior of the (possible) isolated eigenvalue-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-dimensional regime. Assumption 1 (High-dimensional asymptotics). As n, p → ∞ with p/n → c ∈ (0,∞), we have max{‖w‖, ‖w∗‖} = O(1) and xi i.i.d.∼ N (µ,C) with max{‖µ‖, ‖C‖} = O(1).

heading: 2.1 Limiting spectral measure
text: Our first result is the limiting Hessian eigenvalue distribution. This is a direct consequence of our main technical Theorem 4 and is proven in [43, Sec A.2]. Theorem 1 (Limiting spectral measure). Let Assumption 1 hold, we have, as n, p → ∞ with p/n → c ∈ (0,∞), the empirical spectral measure5 µH of the Hessian matrix H defined in (4) converges weakly and almost surely to a probability measure µ, defined through its Stieltjes transform m(z) = ∫ (t− z)−1µ(dt) as the unique solution to6
m(z) = 1
p tr Q̄b(z), δ(z) =
1 n tr ( CQ̄b(z) ) , Q̄−1b (z) ≡ E
[ g ·C
1 + g · δ(z)
] − zIp, (5)
where g ≡ ∂2`(y, h)/∂h2, for h = wTx ∼ N (wTµ,wTCw), (6)
and y and ` defined respectively in (1) and (3). Moreover, if we denote ν the law/distribution of g and assume the empirical spectral measure of C converges to ν̃ as p→∞, then (5) writes
m(z) = ∫ ( −z+t̃ ∫ t
1 + tδ(z) ν(dt)
)−1 ν̃(dt̃), δ(z) = ∫ ct̃
−z + t̃ ∫
t 1+tδ(z) ν(dt)
ν̃(dt̃). (7)
In the form of (7), the (Stieltjes transform of the) limiting Hessian spectral measure µ is determined by the ratio c = lim p/n and the two measures ν and ν̃. This formulation is closely connected to the separable covariance model [38, 6, 55, 17, 76] in RMT. Moreover, if ν(dt) = δ1(t) is a Dirac mass at one, this reduces to the popular sample covariance model [65]; taking further ν̃(dt) = δ1(t) gives the Marc̆enko-Pastur law. See Sec 3.1 for numerical evaluations of these special cases. In particular, the support of the (limiting) Hessian spectrum µ is directly linked to that of ν and ν̃. Remark 1 (Hessian eigen-support). Under Assumption 1, 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 µ is of compact support if and only if both ν and ν̃ have compact support [17, Porposition 3.4], µ may be of unbounded support, depending on the model and the loss.
An example of unbounded µ is the phase retrieval model with y = (wT∗x) 2 and square loss `(y, h) = (y − h2)2/4, for which we have g = 3(wTx)2 − (wT∗x)2 for x ∼ N (µ,C). As a 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 1c, where the Hessian spectrum has a “heavier” tail compared to Fig 1a (logistic model), and the empirically observed “isolated” eigenvalue is due to
5That is, the normalized counting measure of the eigenvalues λi(H) of H, i.e., µH ≡ 1p ∑p i=1 δλi(H).
6Uniqueness is ensured in such a way that =[m(z)] · =[z] > 0 for =[z] 6= 0 and zm(z) < 0 for =[z] = 0, so that (z,m(z)) is a valid Stieltjes transform couple, see more details in [29].
a 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 eigenvalue-eigenvector (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 f(t) = δ|t|≤ for some > 0.
Another example of unbounded µ is when the exponential loss [23] is used. Precisely, consider the logistic model (2) with `(y, h) = exp(−yh), we have that g = exp(−yh) for h ∼ N (wTµ,wTCw) which follows a log-normal distribution and has unbounded support. As such, the (limiting) Hessian eigenvalue distribution µ has also unbounded support. On the other hand, with logistic loss `(y, h) = ln(1 + e−yh), one has g ≤ 1/4 and µ is guaranteed to have bounded support. In Fig 2a and 2b, the empirical Hessian eigenvalues and the limiting distributions are compared for logistic and exponential losses, with a 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 Theorem 1 for different feature covariance structures7. In particular, one may observe a single main bulk with more “compact” Hessian spectrum as in Fig 2c or multiple bulks (two in the case of Fig 2d) with Hessian eigenvalues more “spread-out”, depending on the feature covariance structure ν̃. In the form of (7), the condition for the existence of multi-bulk eigenspectrum has been thoroughly discussed in [17, Sec 3.2–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.

heading: 2.2 Isolated eigenvalues and eigenvectors
text: As discussed in Remark 1, under Assumption 1, the Hessian has bounded (limiting) eigen-support if and only if ν, the law of g, has bounded support. Under this condition (or, after the application of some function f 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 1, assume that the law ν of the random variable g defined in (6) is of bounded support, define
G(z) = I3 + Λ(z)V TQ̄b(z)V ∈ R3×3, (8)
with Q̄b(z), δ(z) defined in (5), V ≡ [µ, Cw∗, Cw] ∈ Rp×3, U ≡ C 1 2 [w∗, w] ∈ Rp×2 and
Λ(z) ≡ E g 1 + g · δ(z)
[ 1 (U+z)T U+z U+z(U+z)T − (UTU)+ ] , z = C− 1 2 (x−µ) ∼ N (0, Ip), (9)
7Covariance 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 λ such that G(λ) has a zero eigenvalue (of multiplicity one), there exists an eigenvalue λ̂ of H such that λ̂− λ a.s.−−→ 0.
Theorem 2 provides an 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 δ(z) defined in Theorem 1. 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 Sec 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 n→∞ with p fixed), we have, by the strong law of large numbers, that H a.s.−−→ E[H], with
E[H] = E[`′′(y,wTx)xxT] = E[g] ·C + V [
1 E[g ·U+z]T E[g ·U+z] U+E[g · (zzT − Ip)](U+)T
] VT.
As a consequence, it is expected that in the large n, p→∞ limit, the top eigenvectors of H could also be related to the columns of V. This is the case in Fig 1d, where the top eigenvector is observed to be a “noisy” version of the model parameter w∗. More precisely, for (λ̂, û) an isolated eigenpair of H, the projection VTûûTV ∈ R3×3 can be shown to be asymptotically close to a deterministic matrix. This measures the “cosine-similarly” between the Hessian isolated eigenvector û with any column of V and consequently the performance of using û as an estimate of, for instance the model parameter w∗ for C = Ip.
This result is given in the following theorem, which is proven in [43, Sec A.4]. Theorem 3 (Isolated eigenvectors). In the setting of Theorem 2, for an isolated eigenvalue-eigenvector pair (λ̂, û) of H and λ the asymptotic position (of λ̂) given in Theorem 2, then
VTûûTV = −VTQ̄b(λ)V · Ξ(λ) + o(1), Ξ(λ) = (vTl,GG′(λ)vr,G)−1 · vr,GvTl,G, for Q̄b(z) and G(z) defined in (5) and (8), respectively, vl,G,vr,G ∈ R3 the left and right eigenvectors of G(λ) associated with eigenvalue zero, and G′(λ) the derivative of G(z) with respect to z evaluated at z = λ.

heading: 2.3 Technical tool: deterministic equivalent
text: Our main technical tool to derive Theorem 1, 2 and 3 is a so-called deterministic equivalent [29, 16] result for the Hessian resolvent Q(z) = (H−zIp)−1, that provides simultaneous access to the Hessian limiting eigenvalue distribution and the behavior of the possible isolated eigenpairs. Precisely, the normalized trace tr Q(z)/p gives the Stieltjes transformmH(z) = ∫ (t−z)−1µH(dt) of the empirical
spectral measure µH of H, from which µH can be recovered. Also, for (λ̂, û) an eigenpair of interest, with Cauchy’s integral formula we have |wTû|2 = − 12πı ∮ Γλ wTQ(z)w dz, for a deterministic vector w ∈ Rp and Γλ a positively oriented contour surrounding only λ̂. As such, for Q̄(z) a deterministic equivalent of Q(z), that is, Q(z) ↔ Q̄(z) with tr A(Q(z) − Q̄(z))/p → 0 and aT(Q(z) − Q̄(z))b → 0 almost surely as n, p → ∞, for A ∈ Rp×p 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 Q̄(z). This result is given as follows, with the proof deferred to [43, Sec A.1]. Theorem 4 (Deterministic equivalent). Let Q(z) ≡ (H− zIp)−1 be the resolvent of H defined in (4). Then, under Assumption 1, as n, p→∞ with p/n→ c ∈ (0,∞),
Q(z)↔ Q̄(z), with Q̄−1(z) = E [
g 1+g·δ(z) (C 1 2 (Ip −PU)C 1 2 + ααT) ] − zIp,
for random vector α ≡ µ + C 12 PUz ∈ Rp and g = `′′(y,wTµ + wTC 1 2 z) for z ∼ N (0, Ip) as defined in (6), y and δ(z) defined in (1) and (5), respectively, and PU ∈ Rp×p the projection onto the subspace spanned by the columns of U ≡ C 12 [w∗, w].

heading: 3 Evaluations and Discussions
text: In this section, we provide further discussions on the consequences of Theorem 1, 2 and 3, together with numerical evaluations. Implications of Theorem 1 on the Hessian eigenvalue distribution is
discussed in Sec 3.1. In Sec 3.2, we discuss the consequences of Theorem 2 and 3 on the possible isolated eigenpairs, for which two fundamentally different phase transitions are characterized.

heading: 3.1 Hessian eigenvalues distribution
text: For a better interpretation of Theorem 1 on the Hessian eigenspectrum, we consider here the special case of C = Ip, and start with the simple setting where the random variable g in (6) is constant, say g = 1: this happens, e.g., when the square loss `(y, h) = (y − h)2/2 is employed. In this case, the Hessian does not dependent on w,w∗ and the Stieltjes transform m(z) is the solution to zcm2(z)− (1− c− z)m(z) + 1 = 0 and corresponds to the Marc̆enko-Pastur law. As long as g is not constant, the limiting Hessian spectrum is, a priori, different from the Marc̆enkoPastur law, even in the C = Ip setting, since the associated Stieltjes transform m(z) is different from the solution to the Marc̆enko-Pastur 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) Marc̆enko-Pastur law. This can be understood with Theorem 1 and is due to the fact that, the distribution of g is more “concentrated” (around some constant, see Fig 3b versus 3d for a comparison between different cases). This is in sharp contrast to Fig 3c where with the exponential loss, the law of g has a much larger spread and the Hessian is therefore away from a Marc̆enko-Pastur-shape.
This “empirical fit” has been observed in [56, Fig 5], where acceleration methods proposed for a Marc̆enko-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 ML models. Nonetheless, it must be pointed out that this “visual approximation” by Marc̆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 (e.g., the phase retrieval model in Fig 1c), in the choice of loss function (e.g., the exponential loss in Fig 2b), or beyond the identity covariance setting (as in Fig 2c and 2d) would induce a Hessian spectrum that is very different from the Marc̆enko-Pastur law. In this vein, our Theorem 1 goes beyond such “loose” Marc̆enko-Pastur approximation and acts as a more accurate first example in the understanding of Hessian in more involved ML models beyond linear least squares that accounts for nonlinear transformations (such as activation function in NNs) and feature statistics.
While Theorem 1 is proven here only for Gaussian features, we conjecture, as is the case for many random matrix asymptotics, that it holds more generally beyond Gaussian distribution, see [43, Fig 5] and more discussions therein.

heading: 3.2 Isolated eigenvalues and their phase transitions
text: In this section, we discuss the implications of Theorem 2 and 3 on the possible isolated eigenvalueeigenvector pairs. More precisely, we show that, different from the classical spiked models extensively studied in RMT literature [3, 5, 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 as a function of the signal strength; here another type of Hessian spike arises due to the underlying G-GLM model (i.e., w∗ and w) and exhibits a rather different behavior.
0 0.2 0.4 0.6 0
2
4
E ig
en va
lu es
Eigenvalues
Limiting law
Right edge
0 0.2 0.4 0.6 0
2
4
E ig
en va
lu es
Eigenvalues Limiting law
Right edge
√ c 1 1.5
0
0.05
0.1
0.15
0.2
‖µ‖2
E ig
en ga
p Empirical Theoretical
√ c 1 1.5
0
0.5
1
‖µ‖2
E ig
en -a
lig n
Figure 4: Spike due to data signal in Corollary 1: theory versus practice of (left) Hessian eigenspectrum with ‖µ‖2 = 0.8, (middle) eigengap dist(λµ, supp(µ)), and (right) top eigenvector alignment α in (10), as a function of the signal strength ‖µ‖2, on logistic model with logistic loss, for µ ∝ [−1p/2, 1p/2], w = w∗ = 0, C = Ip, p = 512 and n = 2 048. Results averaged over 50 runs.
0 0.1 0.2 0.3 0
2
4
6
8
10
E ig
en va
lu es
Eigenvalues
Limiting law
Left edge
0 0.1 0.2 0.3 0
2
4
6
8
10
E ig
en va
lu es
Eigenvalues
Limiting law
Left edge
2 4 6 8 0
0.01
0.02
0.03
‖w‖
E ig
en ga
p Empirical Theoretical
2 4 6 8 0
0.2 0.4 0.6 0.8 1
‖w‖
E ig
en -a
lig n
Figure 5: Left-hand side spike due to response model in Corollary 2 in the absence of data signal: (left) Hessian spectrum for ‖w‖ = 2, with a left isolated eigenvalue λ̂w, (middle) eigengap dist(λw, supp(µ)), and (right) dominant eigenvector alignment (with w), as a function ‖w‖ with w ∝ [−1p/2, 1p/2], w∗ = µ = 0, C = Ip, p = 800 and n = 8 000. Results averaged over 50 runs.

heading: 3.2.1 Spike due to data signal
text: To study the spike due to data “signal” µ and its phase transition behavior, we focus here on the case w∗ = w = 0. This, in the case of logistic model (2) for example, gives rise to a much simpler form of limiting spectrum (per Theorem 1) and possible isolated eigenpairs (per Theorem 2 and 3), as summarized in the following corollary, with detailed derivations given in [43, Sec B.3].
Corollary 1 (Spike due to data signal: logistic model). Consider the logistic model in (2) with logistic loss, for w = w∗ = 0 and C = Ip, the limiting Hessian eigenvalue distribution is the Marc̆enko-Pastur law, but rescaled by a factor of g = 1/4. Moreover, there is at most one isolated eigenpair (λ̂µ, ûµ) of H and it satisfies
λ̂µ a.s.−−→
{ λµ = 1 4 (1 + ρ+ c · ρ+1 ρ ) ρ > √ c,
1 4 (1 +
√ c)2 ρ ≤ √ c;
, |µTûµ|2
‖µ‖2 a.s.−−→
{ α = ρ
2−c ρ2+cρ ρ >
√ c,
0 ρ ≤ √ c;
(10) with the signal strength ρ = limp→∞ ‖µ‖2 and c = lim p/n.
The behavior of the isolated eigen-pairs described in Corollary 1 follows the “classical” phase transition [4, 3, 54]: (i) the isolated eigenvalue always appears on the right-hand side of the main (Marc̆enko-Pastur) bulk and (ii) the eigenvalue amplitude and eigenvector alignment is “monotonic” with respect to the signal strength ‖µ‖2 in the sense that, for a fixed dimension ratio c, the largest Hessian eigenvalue is bound to become asymptotically isolated once ‖µ‖2 exceeds √ c and its value, as well as the eigenvector alignment, increase monotonically as ‖µ‖2 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.

heading: 3.2.2 Spike due to model
text: To investigate the spike due to the underlying model (i.e., w∗ and w), we position ourselves in the situation where µ = 0, that is, in the absence of data “signal”. This leads to 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 (2) with logistic loss, µ = 0 and C = Ip, then the Stieltjes transform m(z) satisfies m(z) = 1/(E[f(r, z)] − z) for f(r, z) = 1/(cm(z) + 2 + e−r + er) and r ∼ N (0, ‖w‖2) that depends on w but not on w∗. Moreover, there is at most one isolated eigenvalue λ̂w of H that is due to w and satisfies λ̂w − λw a.s.−−→ 0 with λw solution to 0 = det G(λw) = 1 +m(λw)E[f(r,λw)(r 2−‖w‖2)]
‖w‖2 .
The situation here is more subtle (than the spike due to data signal discussed in Sec 3.2.1): 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 5 illustrates the behavior of the spike due to w in the setting of Corollary 2. Note first that, different from the case of spike due to data signal µ, the spike in Fig 5-(left) 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.8 Also, note from Fig 5-(middle) that, different from the spike due to µ, the spike due to w exhibits here a “non-monotonic” behavior in the sense that, it is absent for small values of ‖w‖ (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 5-(right) that, the alignment between the associated isolated eigenvector and the parameter w is, unlike the eigengap in Fig 5-(middle), monotonically increasing as ‖w‖ grows large, as in the case of Fig 4-(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 eigenvector, which somehow goes against the conventional eigengap heuristic [67, 35]. It is worthy mentioning that, while, technically speaking, the proposed analysis is not capable of characterizing the behavior as ‖w‖ → ∞ under Assumption 1, empirical results suggest that for extremely large ‖w‖, the eigengap tends to “vanish”, the associated dominant eigenvector can still be used to recover w almost perfectly, see [43, Fig 8] as an example.

heading: 4 Conclusion
text: In this article, we provided a precise asymptotic characterization of Hessian eigenspectra (including the limiting eigenvalue distribution and the behavior of possible isolated eigenvalue-eigenvector 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 is challenging, e.g., one will need to first study the layer-wise Hessian and evaluate how the objective of interest (e.g., the eigenvalue distribution) of the full Hessian depends on each of its (layer-wise) blocks. Simplification can be made by specifying the choice of loss function and/or using different approximations (such as Gauss-Newton-type approximation [63] or working in the infinitely-wide but less-realistic NTK regime [32]), so as 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 (2021XXJS110) and CCF-Hikvision Open Fund (20210008) for providing partial support. MWM would like to acknowledge DARPA, IARPA (contract W911NF20C0035), 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.
8Depending 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-to-right transition. The eigenvector alignment can behave differently from
Fig 4-5 and establish a non-monotonic behavior as a function of ‖w‖ or ‖w∗‖; see [43, Sec 4.3].

Review summaries:
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.   Some other contributions:  (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.   (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$. 

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)$ where $f(\cdot | \cdot)$ is some conditional density function. Given a training set $(x_i, y_i)$ for $i=1, \dots, n$ a popular approach to estimate $w_*$ is to minimize the cost function $ L(w) = \sum_{i=1}^n \ell(y_i, x_i^T w) $ for some loss $\ell$ that is carefully chosen.  This framework is quite general and encompass many classical models such as the standard linear model, the logistic model or the phase retrieval problem.   The authors studies this model in the setting where - $n,p \to \infty$, while $n/p$ converges to some positive number. - $x_i \sim_{iid} \mathcal{N}(\mu, C)$ - $\|w\|, \|w_*\|, \| \mu \|, \|C\| = O(1)$   Under these assumptions, they provide a precise assymptotic caracterization of the eigenvalues and eigenvectors of the Hessian $ H = \frac{\partial^2 L}{\partial w \partial w}(w). $ More precisely:  1. 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. 2. 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$. 3. 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.   The main results are obtained through an so-called 'deterministic equivalent', that allows to study precisely the resolvant of the Hessian $H$.

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.