Analytic Insights into Structure and Rank of Neural Network Hessian Maps
Year: 2022

Abstract:
The Hessian of a neural network captures parameter interactions through secondorder derivatives of the loss. It is a fundamental object of study, closely tied to various problems in deep learning, including model design, optimization, and generalization. Most prior work has been empirical, typically focusing on low-rank approximations and heuristics that are blind to the network structure. In contrast, we develop theoretical tools to analyze the range of the Hessian map, which provide us with a precise understanding of its rank deficiency and the structural reasons behind it. This yields exact formulas and tight upper bounds for the Hessian rank of deep linear networks — allowing for an elegant interpretation in terms of rank deficiency. Moreover, we demonstrate that our bounds remain faithful as an estimate of the numerical Hessian rank, for a larger class of models such as rectified and hyperbolic tangent networks. Further, we also investigate the implications of model architecture (e.g. width, depth, bias) on the rank deficiency. Overall, our work provides novel insights into the source and extent of redundancy in overparameterized neural networks. 1

heading: 1 Introduction
text: Since the very infancy of neural networks, the Hessian matrix has been a central object of study. This is because the Hessian captures pairwise interactions of parameters via second-order derivatives of the loss function. As a result, the Hessian was productively employed, for instance, in (quasiNewton) optimization methods [1, 2], model design and pruning [3–5], generalization [6], network calibration [7], automatically tuning hyper-parameters [8]. But, from the outset the main practical challenge has been its size, scaling quadratically with the model dimensionality. This makes the problem severe for today’s DNNs which have millions or even billions of parameters [9, 10].
Consequently, most prior work has focused on designing scalable Hessian approximations, which either take the route of Hessian-vector products (R-operator) [11–13] or employ positive definite approximations by appealing to the Fisher information matrix. Additional approximations — without exception — are needed on top, such as diagonal approximations [3, 14] in the former or K-FAC [15– 17], restricted to layerwise or arbitrary blocks on the diagonal [18–20], in the latter.
∗Detailed list of contributions are: Sidak first discovered that the Hessian rank formula, in an early form, holds experimentally to high fidelity, thus kick-starting the project. Sidak came up with the proof technique and proved Theorems 3, 5, 9, 12. Sidak wrote essentially the entire paper and noted the rank-deficiency interpretation. Gregor proved Lemma 8, assisted in a part of Theorem 3, empirically observed the eventual formula for the Hessian rank, and essentially ran all the experiments for the final submission and made the corresponding figures. Correspondence to sidak.singh@inf.ethz.ch.
1The most recent version of our paper can be found at https://arxiv.org/abs/2106.16225 and the code corresponding to the experiments is located at https://github.com/dalab/hessian-rank.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
The goal of this paper is to advance the analytical understanding of the Hessian map of a neural network. We pursue the fundamental question of how the model architecture induces structural properties of the Hessian. In particular, we analyze the dimension of its range (i.e., the rank) and identify the sources of rank deficiency. Understanding the range of the Hessian map, in turn, delivers insights into the important aspect of how gradients change between iterations.
A reason why such an important direction currently remains sidelined is that non-linearities in a neural network result in an increased dependence on the data distribution, making a suitable theoretical analysis seem intractable. Following Saxe et al. [21], Kawaguchi [22], who delivered useful general insights on neural networks by looking at the linear case, we take a step back and rigorously characterize the range of the Hessian map and determine the resultant rank deficiency for deep linear networks. The key result of our paper is an exact formula along with tight upper bound on the rank of the Hessian — which effectively depend on the sum of hidden-layer widths. This stands opposed to the total number of parameters which are proportional to the sum of squared layer widths, thus implying a significant redundancy in the parameterization of neural networks (see Fig. 1).
The exact quantification of the Hessian rank gives a precise yet interpretable ballpark on the inherent complexity of neural networks since rank naturally measures the effective number of parameters (as best illustrated by the case of a quadratic minimization problem). This relationship is further reinforced by connection to the classical complexity measure of Gull [23], MacKay [24], which is equivalent to rank for a sufficiently small constant controlling the prior. Therefore, this sheds a novel perspective on the nature and degree of overparameterization in neural networks, and opens up interesting avenues for future investigation.
Contributions. The main contributions of our paper can be summarized as follows: (i) Section 3: We characterize the structure of the Hessian range by exhibiting it via matrix derivatives. (ii) Section 4: We prove tight upper bounds and provide exact formulas on the Hessian rank that are neatly interpretable. To the best of our knowledge, this is the very first time that such formulas and bounds are made available for neural networks. (iii) Section 5.2: In the non-linear case, we show that our (linear) rank formulae faithfully capture the numerical rank. (iv) Section 6: We demonstrate via experiments and theory that such rank bounds also hold throughout the course of training. (v) Section 7.1: In the non-linear case, we also provide a pessimistic yet non-trivial bound, which provably establishes degeneracy of the Hessian at the minimum. (vi) Section 7.2, 5.1: We extend our rank results to the case of bias and investigate the effects of architectural components (such as width, depth, bias) on rank. (vii) Appendix S8: As a by-product, our analysis also reveals interesting properties of the Hessian spectrum, and we prove the presence of additional redundancies due to repeated eigenvalue plateaus.
Related work. The study of the Hessian, in recent times, has been re-invigorated by the empirical observations of Sagun et al. [25, 26] who noted a high degree of degeneracy experimentally, and characterized the spectrum as being composed of a bulk around zero and few outlier eigenvalues. Since then, works such as [27, 28] have scaled the empirical analysis to bigger networks via efficient spectral density calculations and better explained the observations in [26]. However, a drawback is that due to their empirical nature, it remains hard to uncover the exact level of degeneracy via these methods [29], and neither are the exact factors that affect the rank deficiency discerned. Another approach has been to fit models from random matrix theory to match the observed Hessian spectrum in neural networks [29, 30], however, such methods do not properly capture the sharp peak at zero observed in empirical spectra [31]. Yet others have investigated the properties of the spectrum in the asymptotic regime [32, 33], but the bounds, if any, are quite coarse (see Fig. 4). Moreover, a common issue underlying most of these approaches is that they are blind to the layerwise compositional structure of neural networks, and are often motivated by a black-box decomposition of the Hessian.
In regards to indicating the parameter redundancy in neural networks, we are far from being the first work. Prior empirical studies have long reported similar observations, like in the form
of, post-training pruning [3, 4, 19, 20, 34] or inherently contained sub-networks (Lottery Ticket Hypothesis [35]) or intrinsic dimension of the loss landscape [36]. Recently, [37] have experimentally argued for the effective dimensionality from [23, 24] to be a good predictor of double descent [38]. Nevertheless, determining the precise extent of redundancy — and showing it provably — as well as the structural reasons behind it have remained illusory in such approaches.

heading: 2 Setup and Formalism
text: General notation. We use the shorthand, Wk:l, to refer to the matrix product chain Wk · · ·Wl, when k > l. When k < l, Wk:l will stand for the transposed product chain Wk
⊤ · · ·Wl⊤ . Besides, ⊗ denotes the Kronecker product of two matrices, vecr indicates the row-wise vectorization of matrices, rk refers to the rank of a matrix, and Ik is the identity matrix of size k × k. Deep Neural Networks (DNNs). A feedforward DNN is a composition of maps, i.e., F = FL ◦ · · · ◦ F 1, where the l-th layer map F l : RMl−1 → RMl , with input dimension d := M0, output dimension K := ML, total number of hidden neurons M := ∑L−1 l=1 Ml. Each layer map is parameterized by a weight matrix Wl and applies an elementwise activation function σl. So we have,

heading: F l = σl ◦Wl with Wl ∈ RMl×Ml−1 and σl : R→ R .
text: For the sake of tractability, we will often investigate linear DNNs where σl = id, and so F (x) = WL:1x. For compactness, we also represent the entire set of parameters by θ := {W1, · · · ,WL}, and where emphasis requires, we will subscript the DNN map F with it and write Fθ .
Next, assume that we are given a dataset S = {(xi,yi)}Ni=1 of N input-output pairs, drawn i.i.d from an underlying distribution px,y. Our focus will be on the squared loss (MSE), ℓx,y(θ) = 12∥y − ŷ∥ 2 and its residual δx,y := ŷ − y = ∂ℓx, y∂ŷ , where ŷ = Fθ(x) is the DNN prediction. The population loss, L is: L(θ) = E px,y [ℓx,y(θ)] . Finally, we will analyze the Hessian matrix HL = ∂2L ∂θ ∂θ .
Backpropagation in matrix derivatives. As all parameters are collected into matrices, we often work with matrix-matrix derivatives by vectorizing row-wise in the numerator (Jacobian) layout, i.e., ∂Y ∂X := ∂ vecr(Y) ∂ vecr(X)⊤ , see [39] and Appendix S2.5. Alongside this, we use the following rule:
∂AWB ∂W = A⊗ B⊤, e.g. ∂F ∂Wk = WL:k+1 ⊗ x⊤W1:k−1 ∈ RK×MkMk−1 . (1)
By the usual chain rule (backpropagation) one has for a linear DNN, at a sample (x,y):
∂ℓx,y ∂Wk
= [ Wk+1:Lδx,y ] ︸ ︷︷ ︸
backward ∈RMk
· [ Wk−1:1x ]⊤ ︸ ︷︷ ︸ forward ∈RMk−1 = Wk+1:L [ WL:1 xx⊤ − y x⊤ ] W1:k−1 . (2)
The above gradient with respect to Wk is of first order in Wk itself and second order in the other weight matrices. Lastly, let us setup the following shorthand, Ω := E [δx,y x⊤] = E [WL:1xx⊤ − yx⊤] , Σxx := E [xx⊤] , and Σyx := E [yx⊤], which we will use throughout the paper.

heading: 3 Hessian Maps of Linear DNNs

heading: 3.1 Hessian structure
text: The Hessian map has a natural block structure defined by the layers and their dimensionality. In order to leverage this structure, we directly take the derivative of the loss gradient in a matrix-by-matrix fashion. First, consider the k-th diagonal block of the Hessian, which is independent of y and is given by,
HkkL := ∂2L
∂Wk∂Wk = Wk+1:LWL:k+1 ⊗ Wk−1:1 Σxx W1:k−1 .
This follows from the matrix-derivative rule in Eq. (1) along with the Eq. (2) and taking expectation. The calculation of the off-diagonal Hessian blocks (kl-th block of size Mk Mk−1 × Ml Ml−1 )
involves the product rule. Note that the two occurrences of a weight matrix Wl, in Eq. (2), are once non-transposed (in WL:1xx⊤) and once transposed (in Wk+1:L or W1:k−1 respectively for k < l or k > l). For simplicity, let us express these parts without adding them and directly write the other Hessian contribution with respect to the transposed matrix, giving:
∀ k ̸= l, H̃klL := ∂2L
∂Wl∂Wk = Wk+1:LWL:l+1 ⊗ Wk−1:1 Σxx W1:l−1 . (3)
∀ k < l ĤklL := ∂2L
∂Wl⊤∂Wk = Wk+1:l−1 ⊗ Wk−1:1 Ω⊤ WL:l+1 . (4)
∀ k > l ĤklL := ∂2L
∂Wl⊤∂Wk = Wk+1:L ΩW1:l−1 ⊗ Wk−1:l+1 . (5)
Equivalence to Gauss-Newton Decomposition. A common approach is to look at the Hessian map from the perspective of the Hessian chain rule, where we have that, HL = Ho +Hf , with
Ho = E px,y [ ∇θF (x)⊤ [∂2ℓx,y] ∇θF (x) ] , and Hf = E px,y [ K∑ c=1 [∂ℓx,y]c ∇2θ Fc(x) ] .
For the MSE loss, the Gauss-Newton decomposition is in fact equivalent to what we discussed before (see details in Appendix S1), where Ho contains the blocks HkkL and H̃ kl L , while Hf consists of ĤklL (although with the non-transposed matrix). Henceforth, we will refer to the first term Ho as the outer-product Hessian, while we coin the second-term Hf as the functional Hessian.

heading: 3.2 Range of the Hessian map
text: Assuming a local Taylor-series approximation of the loss’ gradient, we have that for ∥∆θ∥ < ϵ , ∇θ+∆θ L ≈ ∇θL+HL ∆θ. This indicates how the gradients will change over any local perturbation ∆θ. Also, this holds over successive iterations of an optimization algorithm such as gradient descent. As a result, this serves to show the significance of the Hessian range as since its dimension (i.e., rank) will dictate the dimension of the space where gradients evolve — something which has been discussed only empirically before [40].
Moving on, to best highlight the layer-wise structure present in the Hessian range, let us multiply the two parts of the Hessian with a vector ∆θ, which we decompose as ∆θ =[ vecr(∆W 1)⊤ . . . vecr(∆W L)⊤ ]⊤ .
Range of the outer-product Ho The product of the k-th row-block Hk•o , corresponding to the kth layer, with ∆θ can be written succinctly as, Hk•o · ∆θ = vecr ( Wk+1:L · ∆̄ ·ΣxxW1:k−1 ) , where ∆̄ := ∑L l=1 W
L:l+1 ∆Wl Wl−1:1 ∈ RK×d. This essentially follows from the identity, vecr AXB = (A⊗B⊤) vecr X (see proof in Appendix S2.5). Note that ∆̄ represents the net change on the prediction map induced by changes to the weight matrices in the forward pass. We see that: (1) ∆θ is linearly compressed into a highly interpretable matrix ∆̄ with Kd entries. (2) The same compressed ∆̄ is shared across all result blocks as it is independent of the row block index k. Overall, this already hints that there is a significant intrinsic structure in the Hessian that constrains its range.
Range of the functional Hessian Hf Here we multiply with ∆̂θ, which is similar to ∆θ except that we consider ∆Wl⊤ instead of ∆Wl . Then the product corresponding to kth layer is, Hk•f · ∆̂θ = vecr ( Wk+1:L Ω [∆<k]⊤ + [∆>k]⊤ ΩW1:k−1 ) , where ∆<k :=∑k−1
l=1 W k−1:l+1 ∆Wl Wl−1:1 and ∆>k := ∑L l=k+1 W
L:l+1 ∆Wl Wl−1:k+1. Notice that the range of Hf is also inherently constrained like in the case of Ho, but there are two important differences. First, the data dependent part is now the covariance Ω = E [δx,y x⊤]. Clearly, if this matrix is low rank (say s), this will directly impact the rank of Hf . The other significant difference is that the weight matrix product is split at layer k. This reflects the fact that upstream and downstream layers have a different effect. We will see ahead that these factors will neatly gives rise to a dependence on the sum of hidden-layer widths.

heading: 4 Main result: Analysis of the Hessian rank
text: Preliminaries. Let us denote the rank of the uncentered covariance Σxx = E [xx⊤] by r. If r < d, then without loss of generality, consider Σxx := (Σxx) r×r, which is always possible by pre-processing the input. Thus, Σxx ≻ 0, always. Also, in such a case, we take W1 ∈ RM1×r. Further, the only assumption we make in our analysis is A1, which is in fact guaranteed at typical initialization almost surely (c.f. Appendix S2.6). In Section 6, we see what happens while training.
Assumption A1. Maximal Rank: ∀ l ∈ [L], Wl ∈ RMl×Ml−1 has rank equal to min(Ml,Ml−1).
Lastly, the omitted proofs in the coming subsections are located in Appendices S2, S3, S4 respectively.

heading: 4.1 Analytical tool
text: The key idea of our analysis technique is to reduce the rank of involved matrices to the rank of a certain special kind of matrix (or its variant), Z, as shown below:
Z = ( Iq ⊗C D⊗ In ) , with C ∈ Rm×n , D ∈ Rp×q . (6)
This row-partitioned matrix has a characteristic structure, where an identity matrix alternates between the two sides of the Kronecker product. Such matrices are in fact omnipresent in the Hessian structure, and importantly for our purpose, they possess additional properties on their rank. Inherent to these properties and our analysis, is the use of generalized inverse [41] and oblique (non-orthogonal) projector matrices. The following Lemma 1 from [42] details such a result: Lemma 1. Let Z be a matrix as in Eq. (6). Then, rk(Z) = q rk(C) + n rk(D)− rk(C) rk(D) .

heading: 4.2 Rank of the outer-product Hessian
text: Consider the following decomposition of Ho, which reveals its ‘outer-product’ nature: Proposition 2. For a deep linear network, Ho = AoBoAo⊤ , where Bo = IK ⊗Σxx ∈ RKd×Kd, and A⊤o = ( WL:2 ⊗ Id · · · WL:l+1 ⊗W1:l−1 · · · IK ⊗W1:L−1 ) ∈ RKd×p ,
where p is the number of parameters. A straightforward consequence is that if there is no bottleneck in between, i.e., no hidden-layer with width Mi < min(K, d), then the matrix Bo will control the rank of Ho. Hence, as a first upper bound we get, rk(Ho) ≤ min (rk(Ao), rk(Bo)) = rk(Bo) = rk(IK) rk(Σxx) = Kr . Such a decomposition is however not new (see [43]), but this only forms an initial step of our analysis and the current bound can be loose in the bottleneck case (e.g., an auto-encoder). Let us define the minimum dimension to be q := min(r,M1, · · · ,ML−1,K). Our main theorem can then be stated as: Theorem 3. Consider the matrix Ao mentioned in Proposition 2. Under the assumption A1, rk(Ao) = r rk(W 2:L) +K rk(WL−1:1)− rk(W2:L) rk(WL−1:1) = q (r +K − q) .
Now, from Theorem 3, it is evident that we can simply upper bound the rank of Ho, by the rank of Ao. But actually, it is possible to show an equality (using Lemma 18), as described below: Corollary 4. Under the setup of Theorem 3, the rank of Ho is given by rk(Ho) = q (r +K − q) .
Remark. It should be emphasized that this analysis not only delivers the rank of the outer-product Hessian but also of Neural Tangent Kernel [44] (see Appendix S9.5), Fisher information matrix, network Jacobian, due to their underlying intimate relationship, unaltered through the lens of rank.

heading: 4.3 Rank of the functional Hessian
text: For analyzing the rank of the functional Hessian Hf , we will continue operating by having one of the derivatives with respect to a transposed weight matrix (here, Wl ⊤ ), since rank does not change with
column or row permutations. We denote this modification of the functional Hessian by Ĥf . Then from Eqs. (4, 5), we can observe that there is a common structure within the blocks contained in the column l. Namely, all the weight matrix-chains have either l − 1 or l + 1 as their right indices. So our approach will be to bound the rank of the individual column blocks, as formalized below:
Theorem 5. For a deep linear network, the rank of l-th column-block, Ĥ•lf , of the matrix Ĥf , under the assumption A1 is given as rk(Ĥ•lf ) = q̂ Ml−1+ q̂ Ml− q̂ 2 , for l ∈ [2, · · · , L−1]. When l = 1, we have rk(Ĥ•1f ) = q̂ M1 + q̂ s− q̂ 2 . And, when l = L, we have rk(Ĥ•Lf ) = q̂ ML−1 + q̂ s− q̂ 2 . Here, q̂ := min(r,M1, · · · ,ML−1,K, s) = min(q, s) and s := rk(Ω) = rk(E [δx,y x⊤]).
The upper bound on the rank of Hf follows by combining the above result over all the columns: Corollary 6. Under the setup of Theorem 5, the rank of Hf can be upper bounded as, rk(Hf ) ≤ 2 q̂ M + 2 q̂ s− L q̂ 2 , where M = ∑L−1 i=1 Mi .
Although the above result is an upper bound, empirically this is precisely the formula at initialization, and thus we have the tightest upper bound on the rank of the functional Hessian. Besides, in general, we find that at initialization s = min(r,K), hence q̂ = q and we will just use q hereafter for simplicity.
Block-column independence. A surprising element of the above analysis is that just adding the ranks of the block-columns, corresponding to the respective layers, gives the rank of the entire Hf which is tight. This phenomenon is quite straightforward to see in a 2-layer network and there Corollary 6 is an equality. However, the more interesting observation is that this holds even for arbitrary depth and also extends to networks with non-linearities such as ReLU and Tanh. This implies that the column spaces associated with the layerwise block-columns do not overlap. Thus, Hf is similar to a block diagonal matrix, and it should be possible to uncover the similarity transformation. But this is beyond the current scope, and we leave it as an open question.

heading: 4.4 Overall bound on the Hessian Rank
text: Finally, in order to get an upper bound on the rank of the entire Hessian, we simply use rk(A+B) ≤ rk(A) + rk(B), along with the Corollary 4 and Corollary 6 (with q̂ = q as noted above), to obtain:
rk(HL) ≤ rk(Ho) + rk(Hf ) ≤ 2 qM − Lq2 + 2 q s+ q (r +K − q) .
Fact 7. The following equality holds empirically: rk(HL) = 2qM − Lq2 + q(r +K) .
This implies that our upper bound is off by just a constant additive factor of 2 q s − q2, which is another startling finding. Thus, revealing that the intersection of the Ho, Hf column spaces has a rather small dimension. E.g., take the typical case of no bottlenecks, q = s, then our upper bound exceeds the true rank by a small constant q2, i.e., minimum dimension squared. This suggests that the direct sum of their column spaces is not too far-fetched as an approximation to the column space of HL. Previously, [30] empirically noted a similar observation for 1-hidden layer networks and [33] showed a high degree of mutual orthogonality in the asymptotic regime. Our result shows a similar consequence in the finite case.
Alternative interpretation. Besides the above result, that the rank of the Hessian is proportional to the sum of widths, there is an alternate way of viewing this. Let us calculate the rank deficiency in the network when (uncentered) input-covariance has rank, i.e., r = d. Using Fact. (7) this comes
out to rank-deficiency(HL) = L−1∑ i=0 (Mi − q)(Mi+1 − q) , whereas the number of parameters is
equal to p = ∑L−1
i=0 Mi Mi+1. Hence this lends an elegant interpretation to our formula, whereby the amount of rank deficiency is equal to the number of parameters of a hypothetical network whose all layer-widths have been subtracted by the minimum layer-width of the original network.

heading: 4.5 Hessian Rank as effective # of parameters
text: Consider a quadratic optimization problem, argminθ∈Rp θ ⊤ Mθ, where assume the matrix M has a reduced rank r < p and is in rank normal form, M = [Mr;0]. As a result, we get an equivalent but reduced problem argminθr∈Rr θ ⊤ r Mr θr of size r < p, where θr is the vector containing the first r parameters. More generally, we can reparameterize θ 7→ U⊤θ where U contains the (top-r) eigenvectors of M. Thus, here it is strikingly clear that rank provides the effective number of parameters. Further, we can go to the case of neural networks by simply utilizing the 2nd order Taylor series expansion of the loss L— and so the Hessian HL plays the role of M. Overall, this suggests that the above results on the Hessian rank can provide insights into the effective # of parameters in neural networks, or in other words, into the inherent complexity of neural networks.

heading: 5 Empirical Results

heading: 5.1 Verification of Rank formulas and their behaviour
text: Setup. We test our results on a variety of datasets: MNIST [45], FashionMNIST [46], CIFAR10 [47]; for various loss types: MSE, cross-entropy, cosh; across several initialization schemes: Glorot [48], uniform, orthogonal [21]. Our theory extends to all these settings. However, due to space constraints we only show a subset of experiments here, but the rest can be found in the Appendix S9.
Procedure. To verify the prediction of our theoretical results, we perform an exact calculation of the rank by computing the full Hessian and the corresponding singular value decomposition (SVD). The available iterative schemes for rank approximation (c.f. [49]), although more memory-efficient, are only effective for well-separated spectra and thus cannot provide useful approximations in the case of neural network Hessians. Besides the exact Hessian computation, we also utilize FLOAT-64 precision to ensure ‘true’ rank calculation, resulting in increased memory costs. Hence, we downscale the image resolution to d = 64 to test on more realistic networks.
Results. We study how the rank varies as a function of the sample size N and the network architecture (for varying widths). Fig. 2 shows this for a linear network on CIFAR10 with MSE loss. First, in Fig 2a, we observe that our predictions match the true rank exactly across all sample sizes as the dependence on N is only exhibited in the rank of the empirical covariance Σ̂xx, confirming that rank is largely a distribution-independent quantity. So, for the rest of our experiments, we sufficiently subsample to ensure that the empirical and true covariance have the same rank. Next, in Fig. 2b, we see that our rank formulas hold for arbitrary-sized network architectures. To contextualize the growth of rank with increasing width, in Fig. 2c, we normalize it by # params p. We notice that rank/# params, which intuitively captures the fraction of effective parameters, saturates down to a small level — signalling the extent of redundancy present in the network parameterization (similar plots across depth and simulation of rank/# params on larger settings can be found in S9.1, S9.2).

heading: 5.2 The case of non-linearities
text: Although the linear nature of the neural network was crucial to our analysis, in this section we show experimentally that our results also extend to the non-linear setting — as numerical rank [50].
Visual comparison of the Hessian spectra. Let us first understand how non-linearities affect the Hessian spectrum relative to the spectrum of linear networks. Fig. 3 compares the spectra of HL, Hf , Ho in these two scenarios (linear vs Tanh), with a zoomed-in inset near the cut-off obtained from rank formulas corresponding to the linear case. We can observe the presence of numerous tiny, but not exactly zero, eigenvalues in the non-linear case. So, if we were to measure the rank with a threshold up to machine precision — as we did in the linear case — this would result in an inflation of the rank measurement for the non-linear scenario. From a practical point of view, tiny but non-zero eigenvalues hold little significance, so a more relevant quantity is the numerical rank [50] that uses a reasonable threshold to weed out such extraneous eigenvalues. The numerical rank, or alternatively the size of bulk around zero, for the non-linear case, indeed seems to be captured by our (linear) rank formulas to high fidelity, as evident from Fig. 3. Similar results for other non-linearities, loss functions, datasets can be found in the Appendix S9.4.
Quantitative measure of the fidelity of Rank formulas. To complement the visual grounds presented above, we now quantitatively measure the (in)significance of such spuriously tiny eigenvalues, in terms of the reconstruction error incurred by excluding them. Hence, we perform a low-rank approximation via the SVD and measure the (relative) reconstruction error at the value of rank from our formulas of the linear case. Fig. 4 displays the reconstruction error as a function of the # of top eigenvectors employed, for a ReLU network. As a reference, we consider the empirical rank measurement obtained at machine precision in this case.
We find that using the linear rank value provides an excellent reconstruction, hence demonstrating the fidelity of our rank formulas to serve as a measure of numerical rank in the non-linear case. The same observation extends to other non-linearities and losses, which we highlight in the Appendix S9.3. Consequently, iterative Hessian estima-
tion procedures, e.g. in second-order optimization methods [16, 51], could benefit from the linear rank as a guiding criterion for their design of Hessian approximation. Besides, these experiments also indicate that previous bounds, such as those by Jacot et al. [33], on the rank of outer-product Hessian and functional Hessian are quite coarse to be of much use (i.e., even > # parameters). This is because these bounds have a linear dependence on the product of: # of samples N and # of classes K. For the same reason, other works [16, 52] that bound the rank of the outer-product Hessian Ho, trivially, based on the # of samples N are of little use.

heading: 6 Evolution of Rank during training
text: The upper bounds on the Hessian rank detailed before inherently depend on the rank of the weight matrices. While initialization guarantees them to be of maximal rank, the rank of weight matrices might possibly decrease during training, thus bringing about a decrease in the Hessian rank. Under some additional assumptions, Lemma 8 shows that this does not happen and the rank remains constant.
Lemma 8. For a deep linear network, consider the gradient flow dynamics Ẇlt = −η∇WlLS(θ) ∣∣ θ =θt . Assume: (a) Centered classes: 1N ∑N i:yic=1 xi = 0, ∀ c ∈ [1, . . . ,K]. (b) Balancedness at initialization: Wk+1 ⊤
0 W k+1 0 = W k 0W
k⊤
0 . (c) Square weight-matrices: Wl ∈ RM×M , ∀ l and K = d = M . Then for all layers l, rk(Wlt) = rk(Wl0), ∀ t <∞ .
Balancedness is a common assumption that has been used in many previous works, like Arora et al. [53]. Centered classes can easily be enforced via a pre-processing step, however empirically this is not required. While the proof (see Appendix S5) holds for square case, empirically we also find this to be true for non-square matrices (Fig. 5, left) and non-linearities as shown in Fig. S2.
Consequence. An implication of this result is that our upper bounds on the Hessian rank remain valid throughout the training. Even if the rank of the weight matrices were to decrease, say in the rectangular case, the Hessian rank would only decrease and our bounds would still be valid. The other avenue for a decrease in rank is learning-driven, and as Ω = E [δx,y x⊤]→ 0, the functional Hessian
Hf → 0 (as it has an explicit dependence on the residual δx,y ) . So, if the learning architecture is powerful enough, the rank of Hf will approach 0, as well as leading to a decrease in the rank of HL, and so rk(HL)→ rk(Ho). In other words, the loss Hessian is completely captured through the outer product Hessian at convergence. We demonstrate this empirically via Fig. 5 (right) where we train a linear network using a linear teacher and visualize the rank dynamics as well as our upper bounds. We provide more details as well as experiments for the non-linear case in the Appendix S5.

heading: 7 Further results

heading: 7.1 Provable Hessian degeneracy with non-linearities
text: From Section 5.2, it is clear that our upper bounds faithfully estimate the numerical rank in the non-linear case, despite the numerous tiny eigenvalues which cause inflated rank measurements. This might raise the question if it is possible to establish, in this non-linear case, that the Hessian has provable degeneracy and not just approximate? The challenge is that here the data distribution manifests additionally via the activations at each layer, thus a theoretical analysis seems intractable without imposing strong assumptions. For a 1-hidden layer network, F (x) = W2σ(W1 x), we show it is still possible to get a pessimistic yet non-trivial upper bound with the following (mild) assumption:
Assumption A2. For each active hidden neuron i, the weighted input covariance has the same rank as the overall input covariance, i.e., rk(E [αx xx⊤]) = rk(Σxx) = r, with αx = σ′(x⊤ W1i •) 2 .
Theorem 9. Consider a 1-hidden layer network with non-linearity σ such that σ(z) = σ′(z)z and let M̃ be the # of active hidden neurons (i.e., probability of activation > 0). Then, under assumption A1 and A2, rank of Ho is given as, rk(Ho) ≤ rM̃ + M̃K − M̃ .
Assumption A2 is rather mild, i.e., in the finite-sample case, it holds as soon as the # of samples N > 2d, for typical initialization of parameter weights. Besides, the class of non-linearities which satisfy the above-mentioned condition includes e.g., ReLU, Leaky-ReLU. Further, this result extends to HL:
Corollary 10. At convergence to the minimum, the rank of the loss Hessian HL, for the same setup as Theorem 9, is upper bounded by: rk(HL) ≤ rM̃ + M̃K − M̃ .
Contrast this with the view from [54], who claim the spectrum to be generically non-degenerate. Or, unlike [55], we establish this without any assumptions on a particular kind of overparameterization.
Fact 11. For multiple hidden-layers, the following generalization of Theorem 9 holds empirically, rk(Ho) ≤ p−M1(d−r)− ∑L−1 i=1 Mi , where p is the # of parameters and assuming no dead neurons.
While these bounds are likely to be quite loose as noticeable from Section 5.2, but more importantly they help establish provable degeneracy of the Hessian at the minimum, with the number of ‘absolutelyflat’ directions (i.e., those in the Hessian null space) in proportion to the sum of hidden-layer sizes.

heading: 7.2 Effect of bias on the rank of Hessian
text: Now, we see how the Hessian rank changes when bias is enabled throughout a deep linear network. We make the following simplifying assumption, which is actually a standard convention in practice.
Assumption A3. The input data has zero mean, i.e., x ∼ px is such that E [x] = 0.
Theorem 12. Under the assumption A1 and A3, for a deep linear network with bias, the rank of Ho is upper bounded as, rk(Ho) ≤ q(r +K − q) +K , where q := min(r,M1, · · · ,ML−1,K).
The proof can be found in the Appendix S7. Empirically, we do not require the input to be mean zero and our upper bound actually holds with equality. Also, we list rank formulas for functional Hessian and the overall loss Hessian in Appendix S7, in the non-bottleneck case. E.g.,
Fact 13. rk (HL) = 2q′M + q′(r+K)− Lq′2 + Lq′ , where q′ := min(r+1,M1, · · · ,ML−1,K).
Here too, rank deficiency has a cleaner interpretation of being equal to the # of parameters in a hypothetical network with bias enabled, but with the minimum dimension q′ (that reflects the homogeneous
coordinate at input) subtracted: rank-deficiency(HL) = L−1∑ i=0 (Mi + 1 − q′)(Mi+1 − q′) .

heading: 8 Conclusion
text: Summary. Our paper provides a precise understanding of how the neural network structure constrains the Hessian range and the resulting rank deficiency. In contrast to the number of parameters which are proportional to layer-widths squared, we obtain that rank is proportional to layer-width. The proof strategy relies on bounding the rank of the two parts of the Hessian separately, i.e., the outerproduct Hessian Ho and the functional Hessian Hf , both of which are replete with the special Z-like structure. The analysis also reveals several striking properties of the Hessian, such as surprisingly small overlap in the column spaces of Ho and Hf , and independence of the layer-wise column blocks in Hf . While our results were derived assuming linear activations, we demonstrate that, even with non-linearities, our formulas faithfully capture the numerical rank.
All in all, by providing fundamental insights into the rank of the Hessian map — which ultimately point to the sources and extent of redundancy in overparameterized neural networks, as a result of its compositional structure — our work paves the way to exciting avenues for future research.
Discussion. In particular, our results merit discussion on some important aspects of deep learning:
(i) Overparameterization: Modern DNNs, with billions of parameters, are in stark contrast to the traditional statistical viewpoint of having # of parameters approximately equal to the # of samples. While several works have argued for measuring model complexity instead through weight norms [56], margins [57], compressibility [58], yet it remains hard to get an interpretable ballpark on the model complexity of neural networks. Since rank intuitively captures the notion of effective parameters, it could be a possible alternative to benchmark overparameterization, e.g. for double descent [38].
(ii) Flatness: A growing number of works [59–61] correlate the choice of regularizers, optimizers, or hyperparameters, with the additional flatness brought about by them at the minimum. However, the significant rank degeneracy of the Hessian, which we have provably established, also points to another source of flatness — that exists as a virtue of the compositional model structure —from the initialization itself. Thus, a prospective avenue of future work would be to compare different architectures based on this inherent kind of flatness.
(iii) Generalization: An interesting observation available from our work is that factors such as width, depth, enabling bias — commonly observed to improve generalization — also result in decreasing the rank/# parameters ratio, see Fig. S8b, S8c, S3. In a similar vein, recent work of [62] has provided a lower bound to the generalization error of statistical estimators in terms of the rank of the Fisher (which is intimately related to the Hessian) divided by # of parameters. Practically, one could use a further relaxation of rank as nuclear norm normalized by the spectral norm, in scenarios with spurious rank inflation. Overall, this suggests the relevance of studying the link between rank and generalization.
Limitations. Although empirically our (linear) rank results faithfully capture the behaviour of Hessian rank with non-linearities, we still lack a rigorous theoretical characterization of the same (which will likely require incorporating additional distribution-specific assumptions). Besides, our analysis is limited to the case of fully-connected networks, and it would be interesting to extend this to other architectures such as convolutional networks, residual networks, etc. Lastly, it still needs to be shown theoretically that our upper-bound on the rank of the functional Hessian holds with equality, despite that this is what we observe in practice.
Acknowledgments and Disclosure of Funding
We would like to thank Nicolò Ruggeri for reviewing a first draft of the paper. Further, we would like to thank the members of DA lab and Bernhard Schölkopf for useful comments. Sidak Pal Singh would also like to acknowledge the financial support from Max Planck ETH Center for Learning Systems.

Review summaries:
This paper provides several results on the rank of the Hessian of neural networks. The main results are:  (A) upper bounds on the rank of the Hessian of a linear neural network (i.e., identity activations). These are proved by separating the Hessian into the sum of the “outer-product Hessian” and the “functional Hessian”, and bounding the rank of these separately. (B) a conjectured formula for the exact Hessian rank that matches experiments, and is also quite close to the upper bounds. (C) experimental evidence that these results essentially hold for neural networks with nonlinear activations as well, even after training. (D) for 1-hidden layer nonlinear networks, under a mild assumption a very weak (but still nontrivial) bound on the rank is obtained. 

This paper studies the rank of Hessian for deep linear networks. Theoretically, with squared loss, and assuming the weight matrices are full-rank, this paper provides rank upper bound for outer-product, functional and full Hessian, and empirically verifies that the bounds are almost tight. This paper further shows empirically that the numerical rank with nonlinear networks can still be captured by the linear rank bound, and that the Hessian rank decreases during training. Finally, Hessian degeneracy results for 1-hidden-layer networks and extension to the case with bias are also provided.

This paper analyzes the rank of the Hessian matrix of fully connected neural networks. Theoretically, for a linear fully-connected neural network with square loss, the authors provided an upper bound for the rank of Hessian with respect to the population loss. They also did many experiments for small neural networks on MNIST and CIFAR-10 to show that their upper bound is almost tight, and works for different types of non-linear activations, losses, and initializations.

This paper derives tight upper bounds on the Hessian rank for deep linear networks. They show that these rank formulae (derived for the linear case) faithfully capture the numerical rank in the non-linear regime. They also demonstrate that these rank bounds hold during the course of training.