The Hessian of a neural network captures parameter interactions through second-order 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, providing us with a precise understanding of its rank deficiency as well as 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 networks.

§ INTRODUCTION

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 (quasi-Newton) optimization methods [Zhu et al., 1997, Setiono and Hui, 1995], model design and pruning [LeCun et al., 1990, Hassibi and Stork, 1992, Hochreiter and Schmidhuber, 1995], generalization [MacKay, 1992], network calibration [Denker and LeCun, 1990], automatically tuning hyper-parameters [LeCun et al., 1992].
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 [He et al., 2016, Devlin et al., 2019].

Consequently, most prior work has focused on designing scalable Hessian approximations, which either take the route of Hessian-vector products (R-operator) [Pearlmutter, 1994, Martens and Sutskever, 2011, Bastien et al., 2012] or employ positive definite approximations by appealing to the Fisher information matrix. Additional approximations — without exception — are needed on top, such as diagonal approximations [LeCun et al., 1990, Schaul et al., 2013] in the former or  K-FAC [Martens and Grosse, 2020, Botev et al., 2017, Goldfarb et al., 2020] restricted to layerwise or arbitrary blocks on the diagonal [Laurent et al., 2018, Dong et al., 2017, Singh and Alistarh, 2020], in the latter.

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.

Hessian spectrum of linear and ReLU networks at initialization. Dashed lines indicate our rank predictions. Results have been averaged over 3 runs.

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., 2014, Kawaguchi, 2016, 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 Figure <ref>).

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. This relationship is further reinforced by connection to the classical complexity measure of Gull, 1989, MacKay, 1992, 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:

* Section <ref>: We characterize the structure of the Hessian range by exhibiting it in the form of matrix derivatives.

* Section <ref>: 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.

* Section <ref>: In the non-linear case, we show that our (linear) rank formulae faithfully capture the numerical rank.

* Section <ref>: We demonstrate via experiments and theory that such rank bounds also hold throughout the course of training.

* Section <ref>: In the non-linear case, we also provide a pessimistic yet non-trivial bound, which provably establishes degeneracy of the Hessian at the minimum.

* Section <ref>, <ref>: We extend our rank results to the case of bias and investigate the effects of architectural components (such as width, depth, bias) on rank.

* Appendix <ref>: 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 (e.g., see Fig. <ref>). =-1

Other kinds of redundancies in the Hessian structure due to the occurrence of repeated eigenvalues. In Appendix <ref>, we discuss the conditions when they arise and the reasons behind it. Additional such results can also be found therein.

Related work.
The study of the Hessian, in recent times, has been re-invigorated by the empirical observations of Sagun et al., 2016, Sagun et al., 2017 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  [Ghorbani et al., 2019, Papyan, 2019] have scaled the empirical analysis to bigger networks via efficient spectral density calculations and better explained the observations in [Sagun et al., 2017]. However, a drawback is that due to their empirical nature, it remains hard to uncover the exact level of degeneracy via these methods [Granziol et al., 2020], 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 [Pennington and Bahri, 2017, Granziol et al., 2020], however, such methods do not properly capture the sharp peak at zero observed in empirical spectra [Granziol, 2020]. Yet others have investigated the properties of the spectrum in the asymptotic regime  [Karakida et al., 2019, Jacot et al., 2019], but the bounds, if any, are quite coarse (see Fig. <ref>). 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 reference to 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 [LeCun et al., 1990, Hassibi and Stork, 1992, Han et al., 2016, Dong et al., 2017, Singh and Alistarh, 2020] or inherently contained sub-networks (Lottery Ticket Hypothesis [Frankle and Carbin, 2019]). Recently, [Maddox et al., 2020] have experimentally argued for the effective dimensionality from [Gull, 1989, MacKay, 1992] to be a good predictor of double descent [Belkin et al., 2019]. Nevertheless, determining the precise extent of redundancy and the structural reasons behind it, have remained illusory in such approaches.

§ SETUP AND FORMALISM

General notation. We use the shorthand, $\wM{k:l}$, to refer to the matrix product chain $\wM k \cdots \wM l$, when $k>l$. When $k < l$, $\wM{k:l}$ will stand for the transposed product chain $\wMt{k} \cdots \wMt{l}$.
Besides, $\kro$ denotes the Kronecker product of two matrices, $\vect_r$ indicates the row-wise vectorization of matrices, $\rank$ refers to the rank of a matrix, and $\Im_{k}$ is the identity matrix of size $k\times k$.

Deep Neural Networks (DNNs). A feedforward DNN is a composition of maps, i.e., $F = \Fmap L \circ \dots \circ \Fmap 1$, where the $l$-th layer map  $\Fmap l: \Re^{M_{l-1}} \to \Re^{M_l}$, with  input dimension $d:=M_0$, output dimension $K:=M_L$, total number of hidden neurons $M:=\sum_{l=1}^{L-1} M_l$. Each layer map is parameterized
by a weight matrix $\wM{l}$ and applies an elementwise activation function $\sigma^{l}$. So we have,
\[\,\Fmap l = \sigma^{l} \circ \wM{l}\,\,\text{with}\,\, \wM l \in \Re^{M_{l} \times M_{l-1}}\,\,\text{and}\,\,\sigma^{l}: \Re \to \Re\,. \]

For the sake of tractability, we will often investigate linear DNNs where $\sigma^l = \text{id}$, and so $F(\x) = \wM{L:1}\x$.
For compactness, we also represent the entire set of parameters by $\btheta:=\{\wM{1}, \cdots, \wM{L}\}$, and where emphasis requires, we will subscript the DNN map $F$ with it and write $\Fn$.

Next, assume that we are given a dataset $S=\{(\x_i,\y_i)\}_{i=1}^{N}$ of $N$ input-output pairs, drawn i.i.d from an underlying distribution $p_{\x, \y}$.
Our focus will be on the squared loss (MSE),  $\ell_{\x,\,\y}(\btheta) = \frac 12 \|\y - \hat\y\|^2$ and its residual $ \pmb \delta_{\x,\,\y} :=\hat\y-\y = \jacobi{\ell_{\x,\,\y}}{\hat\y}$,
where $\hat\y=F_{\btheta}(\x)$ is the DNN prediction.
The population loss, $\mathcal{L}$ is:  ${\mathcal{L}}(\btheta) = \E_{p_{\x, \y}}\left[\ell_{\x,\,\y}(\btheta)\right]\,$.
Finally, we will analyze the Hessian matrix $\HL
= \dfrac{\partial^2 \mathcal{L}}{\partial  \btheta \, \partial  \btheta}$.

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., $\dfrac{\partial \mb Y}{\partial \mb X} := \dfrac{\partial \vect_r(\mb Y)}{\partial \vect_r(\mb X)^\top}\,$, see [Magnus and Neudecker, 2019] and  Appendix <ref>. Alongside this, we use the following  rule:

\begin{align}\label{eq:matrix-derivative}
\frac{\partial \mb A \mb W \mb B}{\partial \mb W} = \mb A \otimes \; \mb B^\top, \quad \text{e.g.} \quad
\jacobi{F}{\wM \idx} &=
\wM{L:\idx+1} \otimes \x^\top \wM{1:\idx-1} \,\, \in \Re^{K \times M_{\idx} M_{\idx-1}} \,.
\end{align}

By the usual chain rule (backpropagation) one has for a linear DNN, at a sample $(\x, \y)$:

\begin{align}
\jacobi{\ell_{\x, \y}}{\wM \idx}
& = \underbrace{\Big[ \wM{\idx+1:L}\pmb \delta_{\x, \y} \Big]}_{\text{backward } \in \Re^{M_{\idx}}} \,\cdot \, \underbrace{\Big[  \wM{{\idx}-1:1} \x\Big]^\top}_{\text{forward } \in \Re^{M_{{\idx}-1}}}
= \wM{\idx+1:L} \Big[  \wM{L:1}\,\x\,\x^\top - \y\,\x^\top \Big] \wM{1:\idx-1}\,.
\label{eq:dnn-grad-linear}
\end{align}

The above gradient with respect to $\wM \idx$ is of first order in $\wM \idx$ itself and second order in the other weight matrices. Lastly, let us setup the following shorthand, $\Om:=\E[\pmb \delta_{\x, \y} \, \x^\top]=\E[\wM{L:1} \x \x^\top - \y \x^\top]\,$, $\covx:=\E[\x\x^\top]\,$, and $\covyx:=\E[\y\x^\top]$, which we will use throughout the paper.

§ HESSIAN MAPS OF LINEAR DNNS

§.§ Hessian structure

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 $\idx$-th diagonal block of the Hessian, which is independent of $\y$ and is given by,

\[
\HLsup{\idx\idx} := \frac{\partial^2 \mathcal{L}}{\partial \wM \idx \partial \wM \idx} \, \,
= \,\,\, \wM{\idx+1:L} \wM{L:\idx+1} \; \otimes \;
\wM{\idx-1:1} \,\covx\, \wM{1:\idx-1} \,.
\]

This follows from the matrix-derivative rule in Eq. (<ref>) along with the Eq. (<ref>) and taking expectation.
The calculation of the off-diagonal Hessian blocks ($kl$-th block of size $M_k \,M_{k-1}\times M_l \,M_{l-1}\,$) involves the product rule. Note that the two occurrences of a weight matrix $\wM{l}$, in Eq. (<ref>), are once non-transposed (in $\wM{L:1}\x \x^\top$) and once transposed (in $\wM{k+1:L}$ or $\wM{1: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:

\begin{align}
\forall \,k \neq l, \quad \HLtildeSup{kl} := \frac{\partial^2 \Loss}{\partial \wM l \partial \wM k}  &= \wM{k+1:L}\wM{L:l+1} \; \otimes \; \wM{k-1:1} \,\covx\, \wM{1:l-1}\,.\label{eq:outer-kl}\\
\forall \,k<l \quad \widehat{\mb H}^{kl}_{\mathcal{L}} := \frac{\partial^2 \Loss}{\partial \wMt l {\partial \wM k}} & = \wM{k+1:l-1} \;\otimes\; \wM{k-1:1} \,{\Om}^\top \,\wM{L:l+1} \,.\label{eq:func-hess-eqn1}\\
\forall \,k>l \quad \widehat{\mb H}^{kl}_{\mathcal{L}} := \frac{\partial^2 \Loss}{\partial \wMt l {\partial \wM k}} & = \wM{k+1:L} \,\Om\, \wM{1:l-1} \;\otimes\;  \wM{k-1:l+1} \,.\label{eq:func-hess-eqn2}
\end{align}

Equivalence to Gauss-Newton Decomposition.$\quad$ 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_{p_{\x, \y}}\left[\nabla_{\btheta} {F(\x)}^\top \; [\partial^2 \ell_{\x,\y} ]\; \nabla_{\btheta} F(\x)\right]\,, \text{\,and} \quad \HF=\E_{p_{\x, \y}}\left[\sum_{c=1}^{K} \,[\partial \ell_{\x, \y}]_{c}\; \nabla^2_{\btheta}\, F_c(\x)\right]\,.\]

For the MSE loss, the Gauss-Newton decomposition is in fact equivalent to what we discussed before (see details in Appendix <ref>), where $\HO$ contains the blocks $\HLsup{kk}$ and $\HLtildeSup{kl}$, while $\HF$ consists of $\HLhatSup{kl}$ (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.

§.§ Range of the Hessian map

Assuming a local Taylor-series approximation of the loss' gradient, we have that for $\|\Delta \btheta\|<\epsilon\,$, $\nabla_{\btheta + \Delta \btheta} \,\mathcal{L} \approx \nabla_{\btheta} \mathcal{L} + \HL \, \Delta \btheta$. This indicates how the gradients will change over any local perturbation $\Delta \btheta$. As a result, this also holds over successive iterations of an optimization algorithm such as gradient descent, and serves to show the significance of the Hessian range. Let us multiply the two parts of the Hessian with a vector $\Delta \btheta$, which we decompose as $\Delta \btheta = \left[\vect_r(\Delta \wM 1)^\top \,\, \dots \,\, \vect_r(\Delta \wM L)^\top\right]^\top \,$ to best reflect the layerwise structure.

Range of the outer-product $\HO$ The product of the $k$-th row-block $\HOsup{k\bullet}$, corresponding to the $k^{\text{th}}$ layer, with $\Delta \btheta$ can be written succinctly as,

\[
\HOsup{k\bullet} \,\cdot\, \Delta \btheta
= \vect_r\left(\wM{k+1:L} \cdot \bar \Delta \cdot \covx \wM{1:k-1}\right)\,, \,\,\text{where}\,\,
\bar\Delta  := \sum_{l=1}^L \wM{L:l+1} \,\,\Delta \wM  l\,\,
\wM{l-1:1} \,\in \Re^{K \times d}.\]

This essentially follows from the identity, $\vect_r \Am \Xm \Bm = (\Am \kro \Bm^\top) \vect_r\Xm$ (see proof in Appendix <ref>).
Note that $\bar \Delta$ represents the net change on the prediction map induced by changes to the weight matrices in the forward pass. We see that: (1) $\Delta \btheta$ is linearly compressed into a highly interpretable matrix $\bar\Delta$ with $Kd$ entries. (2) The same compressed $\bar\Delta$ 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 $\widehat{\Delta \btheta}$, which is similar to $\Delta \btheta$ except that we consider ${\Delta \W^l}^\top$ instead of $\Delta \W^l\,$.
Then the product corresponding to $k^{\text{th}}$ layer is,

\[
\HFsup{k\bullet} \cdot \widehat{\Delta \btheta}
= \vect_r \left(\wM{k+1:L} \,\Om\, [\Delta^{<k}]^\top
\,+\, [\Delta^{>k}]^\top\,\Om\,\wM{1:k-1}\right)\,,\]

where $\Delta^{<k}:=\sum_{l=1}^{k-1}
\wM{k-1:l+1}\,\,\Delta \wM{l}\,\, \wM{l-1:1} $ and $\Delta^{>k}:=\sum_{l=k+1}^L \wM{L:l+1} \,\,\Delta \wM{l}\,\, \wM{l-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 $\Om=\E[\pmb \delta_{\x,\y}\, \x^\top]$. 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.

§ MAIN RESULT: ANALYSIS OF THE HESSIAN RANK

Let us denote the rank of the uncentered covariance $\covx = \E[\x \x^\top]$ by $r$. If $r<d$, then without loss of generality, consider $\covx := (\covx)_{\,r\times r}$,
which is always possible by pre-processing the input. Thus, $\covx \succ 0$, always. Also, in such a case, we
take $\wM{1}\in\mathbb{R}^{M_1 \times r}$. Further, the only assumption we make in our analysis is <ref>, which is in fact guaranteed at typical initialization with high probability (c.f. Appendix <ref>). In Section <ref>, we see what happens while training. =-1

Maximal Rank: $\,\forall \, l \in [L]$, $\wM{l} \in \mathbb{R}^{M_l \times M_{l-1}}$ has rank equal to $\min(M_l, M_{l-1})$.

Lastly, the omitted proofs in the coming subsections are located in Appendices <ref>, <ref>, <ref> respectively.

§.§ Analytical tool
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), $\Zm$, as shown below:

\begin{align}
\Zm = \begin{pmatrix}\Im_q \otimes \Cm \\ \Dm \otimes \Im_n\end{pmatrix}\,,\quad \text{with} \,\, \Cm\in\mathbb{R}^{m\times n}\,,\,\,\Dm\in\mathbb{R}^{p\times q}\,. \label{eq:z-structure}
\end{align}

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 [Rao, 1972] and oblique (non-orthogonal) projector matrices. The following  Lemma <ref>  from [Chuai and Tian, 2004] details such a result:

Let $\Zm$ be a matrix as in Eq. (<ref>). Then,
\rank(\Zm)=q\, \rank(\Cm)+n\, \rank(\Dm)-\rank(\Cm)\, \rank(\Dm)\,.$

§.§ Rank of the outer-product Hessian

Consider the following decomposition of $\HO$, which reveals its `outer-product' nature:

For a deep linear network, $\,\HO = \Am_o  \Bm_o {\Am_o}^\top\,$, where $\,\Bm_o = \Im_K \kro \covx\, \in \mathbb{R}^{Kd\times Kd}$,

\[\text{and}\,\,\,\Am_o^\top=\begin{pmatrix}
\wM{L:2} \kro \Im_d \quad \cdots \quad
\wM{L:l+1} \kro \wM{1:l-1}
\quad \cdots \quad
\Im_K \kro \wM{1:L-1}
\end{pmatrix}\,\, \in \mathbb{R}^{Kd\times 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 $M_i<\min(K, d)$, then the matrix $\Bm_o$ will control the rank of $\HO$.
Hence, as a first upper bound we get,
$\,\rank(\HO)\leq \min\left(\rank({\Am_o}), \rank(\Bm_o)\right) = \rank(\Bm_o) = \rank(\Im_K) \, \rank(\covx) = K r\,$. =-1

Such a decomposition is however not new (see [Bernacchia et al., 2018]), 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, M_1, \cdots, M_{L-1}, K)$. Our main theorem can then be stated as: =-1

Consider the matrix $\Am_o$ mentioned in Proposition <ref>.
Under the assumption <ref>,
\rank(\Am_o) =
r \rank(\wM{2:L}) + K \rank(\wM{L-1:1}) - \rank(\wM{2:L})\rank(\wM{L-1:1}) = q \,( r + K  - q)\,.

Now, from Theorem <ref>, it is evident that we can simply upper bound the rank of $\HO$, by the rank of $\Am_o$. But actually, it is possible to show an equality (using Lemma <ref>), as described below:

Under the setup of Theorem <ref>, the rank of $\HO$ is given by $\rank(\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 [Jacot et al., 2018] (see Appendix <ref>), Fisher information matrix, network Jacobian, due to their underlying intimate relationship, unaltered through the lens of rank.

§.§ Rank of the functional Hessian

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, $\wMt{l}$), since rank does not change with column or row permutations. We denote this modification of the functional Hessian  by $\HFhat$. Then from Eqs. (<ref>, <ref>), 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:

For a deep linear network, the rank of $l$-th column-block, $\HFhatsup{\bullet l}$, of the matrix $\HFhat$, under the assumption <ref> is given as
$\rank(\HFhatsup{\bullet l}) = q\, M_{l -1} + q\,M_{l}  - q^2\,,\,$ for $l \in [2, \cdots , L-1]$.
When $l=1$, we have $
\rank(\HFhatsup{\bullet 1}) =q\, M_{1} + q \, s - q^2\,.$ And, when $l=L$, we have $
\rank(\HFhatsup{\bullet L}) =q\, M_{L-1} + q \, s - q^2\,.$
Here, $q := \min(r, M_1, \cdots, M_{L-1}, K, s)$ and $s:=\rank(\Om)=\rank(\E[\pmb \delta_{\x,\y}\,\x^\top])$.

The upper bound on the rank of $\HF$ follows by combining the above result over all the columns:

Under the setup of Theorem <ref>, the rank of $\HF$ can be upper bounded as,
\qquad\qquad\qquad\qquad\rank(\HF) \leq 2\, q\, M  + 2 \,q\, s - L\, q^2\,, \quad \text{where} \quad M=\sum_{i=1}^{L-1} M_i\,.$

Although the above corollary states an upper bound, empirically we find that this is precisely the formula at initialization, and thus we have the tightest upper bound on the rank of the functional Hessian. Also, we find that usually at initialization $s=\min(r, K)$, although this is not needed for the proof. =-1

Block-column independence.$\quad$ 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 <ref> 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.

§.§ Overall bound on the Hessian Rank

Finally, in order to get an upper bound on the rank of the entire Hessian, we simply use $\rank(\Am + \Bm) \leq \rank(\Am) + \rank(\Bm)$, along with the Corollary <ref> and Corollary <ref>, to obtain:

\[
\rank(\HL) \leq \rank(\HO) + \rank(\HF) \,\leq\, 2\, q\, M  - L\, q^2  + 2 \,q\, s + q\,( r + K - q)\,.
\]

The following equality holds empirically:
$\, \colorboxed{gray!50}{\rank(\HL) = 2q \, M - L\, q^2  + q(r + K)}\,.$

This implies that our upper bound is off by just a constant additive factor of $2\,q\,s - q^2$, 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 $q^2$, 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, [Pennington and Bahri, 2017] empirically noted a similar observation for 1-hidden layer networks and  [Jacot et al., 2019] showed a high degree of mutual orthogonality in the asymptotic regime. Our result shows a similar consequence in the finite case.=-1

Alternative interpretation.$\quad$ 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. (<ref>) this comes out to

\[\rankdef(\HL) = \sum\limits_{i=0}^{L-1} (M_i \,- \,q)(M_{i+1}\,-\,q)\,,\]

whereas the number of parameters is equal to $p=\sum_{i=0}^{L-1} M_i \, M_{i+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.

§ EMPIRICAL RESULTS

§.§ Verification of Rank formulas and their behaviour

Rank vs # samples $N$

Rank vs minimal width $M_{*}$

Rank vs depth $L$

Behaviour of rank on CIFAR10 using MSE, with hidden layers: $50, 20, 20, 20$ (Fig. <ref>), $M_{*}, M_{*}$ (Fig. <ref>) and $L$ layers of width $M=25$ (Fig. <ref>). The lines indicate the true value and circles denote our formula predictions.

Setup.$\quad$ We test our results on a variety of datasets: MNIST [LeCun and Cortes, 2010], FashionMNIST [Xiao et al., 2017], CIFAR10 [Krizhevsky, 2009]; for various loss types: MSE, cross-entropy, cosh; across several initialization schemes: Glorot [Glorot and Bengio, 2010], uniform, orthogonal [Saxe et al., 2014]. 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 <ref>.

Procedure.$\quad$ 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. [Ubaru and Saad, 2016]), 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.$\quad$We study how the rank varies as a function of the sample size $N$ and the network architecture (for varying widths). Fig. <ref> shows this for a linear network on CIFAR10 with MSE loss. First, in Fig <ref>, 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 $\samplecovx$, 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. <ref> and <ref>, we see that our rank formulas hold for arbitrary-sized network architectures, across varying width and depth.

$\effparams$ vs minimal width $M_{*}$

$\effparams$ vs depth $L$

Behaviour of rank/#params on CIFAR10 using MSE, with hidden layers: $M_{*}, M_{*}$ (Fig. <ref>) and $L$ layers of width $M=25$ (Fig. <ref>). The lines indicate the true value and circles denote our formula predictions.

To contextualize the growth of rank with increasing architecture sizes, in Fig. <ref>, we normalize it by # parameters $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.

§.§ Simulation of rank formulas for large settings

In the previous subsection, we have established how our formulas hold exactly in practice. An added benefit of these formulas is that they allow us to visualize how the rank will growth in relation to the number of parameters for bigger architectures — without actually having to do the Hessian computations. In Fig. <ref> we show such a simulation for increasing width and depth.  The simulations make the limiting behaviour of the fraction $\frac{\text{rank}}{\# params}$ even more apparent, as the fraction decreases with more and more overparametrization (in terms of both depth and width), until it reaches a threshold.

The left subfigure, which is the width-simulation plot, also shows an interesting behaviour. In the early phase, the outer-product Hessian $\HO$ dominates the functional Hessian $\HF$ in terms of $\frac{\text{rank}}{\# params}$, but after a certain width $\HF$ starts to dominate $\HO$ and continues to do so throughout. It would be of relevance for future work to further investigate the interaction between $\HO$ and $\HF$, and provide an understanding of these two phases.

Simulating rank/#params over width (left) and depth (middle). For the width plot: $L=4, r=2352$, and for the depth plot: $r=2352$ and hidden layer size $5000$.

§.§ The case of non-linearities

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 [Golub et al., 1976].

Spectrum of the loss Hessian $\HL$ (left), functional Hessian $\HF$ (middle)  and outer product $\HO$ (right), for linear and non-linear networks. Black dashed lines are the predictions of the bulk size via our rank formulas. We use 2 hidden layers of size $30,20$ with Tanh activation on MNIST.

Visual comparison of the Hessian spectra.$\quad$
Let us first understand how non-linearities affect the Hessian spectrum relative to the spectrum of linear networks. Fig. <ref> 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 [Golub et al., 1976] 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. <ref>. Similar results for other non-linearities, loss functions,  datasets can be found in the Appendix <ref>.

Hessian reconstruction error for a ReLU network of hidden layer sizes $30, 20$ as the rank of the approximation is increased. As before, the three sub-figures reflect this for the loss Hessian $\HL$ (left), functional Hessian $\HF$ (middle),  and outer product $\HO$ (right) respectively. The x-axis represents the number of top eigenvectors that form the low-rank approximation.
The dashed vertical lines indicate the cut-off at various values of the rank: first line at the prediction based on the linear model, second line at the empirical measurement of rank, and third line based on upper bounds from [Jacot et al., 2019], which become too coarse to be of any use (actually even $>$ # parameters).

Quantitative measure of the fidelity of Rank formulas.$\quad$
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. <ref> 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 ($0.002\%$ error in case of $\HL$), 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 <ref>.
Consequently, iterative Hessian estimation procedures, e.g. in second-order optimization methods [Botev et al., 2017, Lacotte et al., 2021], 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., 2019, on the rank of outer-product Hessian and functional Hessian are quite coarse to be of much use. 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 [Shen, 2016, Botev et al., 2017] that bound the rank of the outer-product Hessian $\HO$, trivially, based on the # of samples $N$ are of little use.

§ EVOLUTION OF RANK DURING TRAINING

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 <ref> shows that this does not happen and the rank remains constant.=-1

For a deep linear network, consider the gradient flow dynamics $\dot{\Wm}^{l}_t = -\eta \nabla_{\Wm^{l}}\mathcal{L}_S(\bm{\btheta})\big{|}_{\bm{\btheta} \,=\, \bm{\btheta}_t}$. Assume: (a) Centered classes: $ \frac{1}{N} \sum_{i:y_{ic} = 1}^{N}\x_i= \bm{0}, \hspace{2mm} \forall\, c \in [1,\dots,K]$. (b) Balancedness at initialization: $\Wm^{{k+1}^\top}_0\Wm_0^{k+1} = \Wm^{k}_0\Wm^{k^\top}_0$. (c) Square weight-matrices: $\Wm^{l} \in \mathbb{R}^{M \times M}\,,$ $\,\forall \,l\,$ and $\,K=d=M$. Then for all layers $l$,
$\,\, \rank(\Wm^{l}_t) = \rank(\Wm^{l}_0), \hspace{3mm} \,\forall \, t<\infty\,.$

Balancedness is a common assumption that has been used in many previous works, like Arora et al., 2019. Centered classes can easily be enforced via a pre-processing step, however empirically this is not required. While the proof (see Appendix <ref>) holds for square case, empirically we also find this to be true for non-square matrices and non-linearities as shown in Fig. <ref> (and more in Appendix <ref>).

Rank dynamics for a ReLU network with hidden layer sizes $25, 20, 15$, trained on a subset of MNIST. We show the evolution of the rank of the weights (left) with usual Glorot initialization (guaranteeing maximal rank) and the rank of the Hessians as a function of training time (right).

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 $\Om=\E[\pmb \delta_{\x, \y} \, \x^\top] \to 0$, the functional Hessian $\HF\to \mb 0$ (as it has an explicit dependence on the residual $\pmb \delta_{\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 $\rank(\HL) \to \rank(\HO)$. In other words, the loss Hessian is completely captured through the outer product Hessian at convergence. These observations also extend to the non-linear case, as shown in Fig. <ref>.=-1

§ FURTHER RESULTS

§.§ Provable Hessian degeneracy with non-linearities

From Section <ref>, 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) = \wM{2}\sigma(\wM{1}\,\x)$, we show it is still possible to get a pessimistic yet non-trivial upper bound
with the following (mild) assumption:=-1

For each active hidden neuron $i$, the weighted input covariance has the same rank as the overall input covariance, i.e.,  $\rank(\E[\alpha_{\x}\,\x \x^\top]) = \rank(\covx) = r$, with $\,\alpha_{\x} = {\sigma^\prime(\x^\top \, \wM{1}_{i\,\bullet})}^2$. =-1

Consider a 1-hidden layer network with non-linearity $\sigma$ such that $\sigma(z) = \sigma^{\prime}(z) z$ and let $\widetilde{M}$ be the # of active hidden neurons (i.e., probability of activation $>0$). Then, under assumption <ref> and <ref>,  rank of $\HO$ is given as,
$\rank(\HO) \leq r\widetilde{M} \,+\, \widetilde{M} K \,-\, \widetilde{M}\,$.

Assumption <ref> is rather mild, i.e., in the finite-sample case, it holds as soon as the # of samples $N> 2 d$, 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$: =-1

At convergence to the minimum, the rank of the loss Hessian $\HL$, for the same setup as Theorem <ref>, is upper bounded by: $\rank(\HL) \leq r\widetilde{M} \,+\, \widetilde{M} K \,-\, \widetilde{M}\,$.

Contrast this with the view from [Pennington and Worah, 2018], who claim the spectrum to be generically non-degenerate. Or, unlike [Poggio et al., 2018], we establish this without any assumptions on a particular kind of overparameterization.

For multiple hidden-layers, the following generalization of Theorem <ref> holds empirically, $\rank(\HO) \leq p - M_1 (d-r) - \sum_{i=1}^{L-1} M_i\,$, where $p$ is the # of parameters and assuming no dead neurons. =-1

While these bounds are likely to be quite loose as noticeable from Section <ref>, but more importantly they help establish provable degeneracy of the Hessian at the minimum, with the number of `absolutely-flat' directions (i.e., those in the Hessian null space) in proportion to the sum of hidden-layer sizes.=-1

§.§ Effect of bias on the rank of Hessian

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. =-1

The input data has zero mean, i.e., $\x \sim p_{\x}$ is such that $\E[\x]=0$.

Under the assumption <ref> and <ref>,  for a deep linear network with bias, the rank of $\HO$ is upper bounded as,
$\rank(\HO) \leq q(r+K-q) + K\,$, where $q:= \min(r, M_1, \cdots, M_{L-1}, K)$.

The proof can be found in the Appendix <ref>. 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 <ref>, in the non-bottleneck case. E.g.,

\rank\left(\HL\right) =
2q^\prime M \,+\, {q^\prime}(r + K) \,-\, L {q^\prime}^2 \,+\, L{q^\prime}
\,$, where $q^\prime := \min(r+1, M_1, \cdots, M_{L-1}, K)$.

Here as well, rank deficiency has a cleaner interpretation of being equal to the # of parameters in a hypothetical network with bias enabled, albeit with the minimum dimension $q^\prime$ (that reflects the homogeneous coordinate at input) subtracted:

\[\quad\rankdef(\HL) = \sum\limits_{i=0}^{L-1} (M_i + 1 \,- \,q^\prime)(M_{i+1}\,-\,q^\prime)\,.\]

Empirical verification. We briefly discuss some empirical experiments that verify the accuracy of our Hessian rank formulas in case of bias. Fig. <ref> illustrates this for network architectures of arbitrary depth using the CIFAR10 dataset and MSE loss. Thus, it demonstrates that our rank formulas — also in the case of bias — exactly predict the empirically observed Hessian rank. Further results, e.g., across arbitrary width as well as number of samples are located in the Appendix <ref>.

$\effparams$ vs depth $L$

Effect of bias on $\effparams$ for varying depth

The fraction of effective parameters, i.e., rank/#params as a function of depth $L$ for MSE loss on down-scaled CIFAR10, using bias. For the left Fig. <ref>, we use network architectures with $L$ hidden layers each of width $M=25$. In the case of Fig. <ref> on the right, we simulating this ratio over depth for $\HL, \HLsup{\text{bias}}$ — in other words, with bias disabled and bias enabled. For this, we use $M=1000$, $d=r=784$ and $K=100$.

Next, in Fig. <ref> we showcase the resulting effect of enabling bias in the network on the Hessian rank, by simulating the ratio $\frac{\text{rank}}{\# params}$ across increasing depth, for the loss Hessian $\HL$ with and without bias. We find that in both cases the $\frac{\text{rank}}{\# params}$ curve saturates to a small threshold. But interestingly, we see that enabling bias further results in a decrease in this ratio.

§ CONCLUSION

Summary.$\quad$ 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 outer-product Hessian $\HO$ and the functional Hessian $\HF$, both of which are replete with the special $\Zm$-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, our work delivers important insights into the nature and degree of parameterization contained in a neural network as a result of its compositional structure.

Discussion.$\quad$ Our results merit discussion on some of the fundamental 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  [Neyshabur et al., 2015], margins [Bartlett et al., 2017], compressibility [Arora et al., 2018], 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 [Belkin et al., 2019].

(ii) Flatness:
A growing number of works [Keskar et al., 2016, Jastrzębski et al., 2018, Baldassi et al., 2020] 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. <ref>, <ref>, <ref>. In a similar vein, recent work of [Seroussi and Zeitouni, 2021] 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. =-1

Besides the above aspects, there are many other interesting questions that naturally arise. These include, to list a few: (a) finding the Hessian rank for convolutional networks, (b) using the rank formulae to build informed models of the Hessian spectrum [Pennington and Bahri, 2017, Granziol et al., 2020], (c) obtaining an equality on the bounds — which is what we observe empirically at initialization, (d) better understanding the interaction between $\HO$ and $\HF$. To conclude, by providing fundamental insights into the range of the Hessian map, our work paves the way to exciting avenues for future research.

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.

[Zhu et al., 1997]
Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal.
Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale
bound-constrained optimization.
ACM Transactions on Mathematical Software (TOMS), 230
(4):0 550–560, 1997.

[Setiono and Hui, 1995]
Rudy Setiono and Lucas Chi Kwong Hui.
Use of a quasi-newton method in a feedforward neural network
construction algorithm.
IEEE Transactions on Neural Networks, 60 (1):0
273–277, 1995.

[LeCun et al., 1990]
Yann LeCun, John S Denker, and Sara A Solla.
Optimal brain damage.
In Advances in neural information processing systems, pages
598–605, 1990.

[Hassibi and Stork, 1992]
Babak Hassibi and David G. Stork.
Second order derivatives for network pruning: Optimal brain surgeon.
In NIPS, pages 164–171, 1992.

[Hochreiter and Schmidhuber, 1995]
Sepp Hochreiter and Jürgen Schmidhuber.
Simplifying neural nets by discovering flat minima.
In Advances in neural information processing systems, pages
529–536, 1995.

[MacKay, 1992]
David J. C. MacKay.
A Practical Bayesian Framework for Backpropagation Networks.
Neural Computation, 40 (3):0 448–472, 05
ISSN 0899-7667.
URL <https://doi.org/10.1162/neco.1992.4.3.448>.

[Denker and LeCun, 1990]
John S. Denker and Yann LeCun.
Transforming neural-net output levels to probability distributions.
In Proceedings of the 3rd International Conference on Neural
Information Processing Systems, NIPS'90, page 853–859, San Francisco, CA,
USA, 1990. Morgan Kaufmann Publishers Inc.
ISBN 1558601848.

[LeCun et al., 1992]
Y. LeCun, P. Simard, and Barak A. Pearlmutter.
Automatic learning rate maximization by on-line estimation of the
hessian's eigenvectors.
In NIPS 1992, 1992.

[He et al., 2016]
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.
Deep residual learning for image recognition.
IEEE Conference on Computer Vision and Pattern Recognition
(CVPR), 2016.

[Devlin et al., 2019]
Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.
Bert: Pre-training of deep bidirectional transformers for language
Proceedings of NAACL-HLT, 2019.

[Pearlmutter, 1994]
Barak A Pearlmutter.
Fast exact multiplication by the hessian.
Neural computation, 60 (1):0 147–160, 1994.

[Martens and Sutskever, 2011]
James Martens and Ilya Sutskever.
Learning recurrent neural networks with hessian-free optimization.
In ICML, 2011.

[Bastien et al., 2012]
Frédéric Bastien, Pascal Lamblin, Razvan Pascanu, James Bergstra, Ian
Goodfellow, Arnaud Bergeron, Nicolas Bouchard, David Warde-Farley, and Yoshua
Theano: new features and speed improvements.
arXiv preprint arXiv:1211.5590, 2012.

[Schaul et al., 2013]
Tom Schaul, Sixin Zhang, and Yann LeCun.
No more pesky learning rates, 2013.

[Martens and Grosse, 2020]
James Martens and Roger Grosse.
Optimizing neural networks with kronecker-factored approximate
curvature, 2020.

[Botev et al., 2017]
Aleksandar Botev, Hippolyt Ritter, and David Barber.
Practical gauss-newton optimisation for deep learning.
In International Conference on Machine Learning, pages
557–565. PMLR, 2017.

[Goldfarb et al., 2020]
Donald Goldfarb, Yi Ren, and Achraf Bahamou.
Practical quasi-newton methods for training deep neural networks.
arXiv preprint arXiv:2006.08877, 2020.

[Laurent et al., 2018]
César Laurent, Thomas George, Xavier Bouthillier, Nicolas Ballas, and Pascal
An evaluation of fisher approximations beyond kronecker
factorization, 2018.
URL <https://openreview.net/forum?id=ryVC6tkwG>.

[Dong et al., 2017]
Xin Dong, Shangyu Chen, and Sinno Jialin Pan.
Learning to prune deep neural networks via layer-wise optimal brain
surgeon, 2017.

[Singh and Alistarh, 2020]
Sidak Pal Singh and Dan Alistarh.
Woodfisher: Efficient second-order approximation for neural network
compression, 2020.

[Saxe et al., 2014]
Andrew M. Saxe, James L. McClelland, and Surya Ganguli.
Exact solutions to the nonlinear dynamics of learning in deep linear
neural networks, 2014.

[Kawaguchi, 2016]
Kenji Kawaguchi.
Deep learning without poor local minima, 2016.

[Gull, 1989]
Stephen F. Gull.
Developments in Maximum Entropy Data Analysis, pages 53–71.
Springer Netherlands, Dordrecht, 1989.
ISBN 978-94-015-7860-8.
URL <https://doi.org/10.1007/978-94-015-7860-8_4>.

[MacKay, 1992]
David MacKay.
Bayesian model comparison and backprop nets.
In J. Moody, S. Hanson, and R. P. Lippmann, editors, Advances
in Neural Information Processing Systems, volume 4. Morgan-Kaufmann,

[Sagun et al., 2016]
Levent Sagun, Leon Bottou, and Yann LeCun.
Eigenvalues of the hessian in deep learning: Singularity and beyond.
arXiv preprint arXiv:1611.07476, 2016.

[Sagun et al., 2017]
Levent Sagun, Utku Evci, V Ugur Guney, Yann Dauphin, and Leon Bottou.
Empirical analysis of the hessian of over-parametrized neural
arXiv preprint arXiv:1706.04454, 2017.

[Ghorbani et al., 2019]
Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao.
An investigation into neural net optimization via hessian eigenvalue
In International Conference on Machine Learning, pages
2232–2241. PMLR, 2019.

[Papyan, 2019]
Vardan Papyan.
Measurements of three-level hierarchical structure in the outliers in
the spectrum of deepnet hessians.
arXiv preprint arXiv:1901.08244, 2019.

[Granziol et al., 2020]
Diego Granziol, Timur Garipov, Dmitry Vetrov, Stefan Zohren, Stephen Roberts,
and Andrew Gordon Wilson.
Towards understanding the true loss surface of deep neural networks
using random matrix theory and iterative spectral methods, 2020.
URL <https://openreview.net/forum?id=H1gza2NtwH>.

[Pennington and Bahri, 2017]
Jeffrey Pennington and Yasaman Bahri.
Geometry of neural network loss surfaces via random matrix theory.
In Doina Precup and Yee Whye Teh, editors, Proceedings of the
34th International Conference on Machine Learning, volume 70 of
Proceedings of Machine Learning Research, pages 2798–2806. PMLR,
06–11 Aug 2017.
URL <http://proceedings.mlr.press/v70/pennington17a.html>.

[Granziol, 2020]
Diego Granziol.
Beyond random matrix theory for deep networks.
ArXiv, abs/2006.07721, 2020.

[Karakida et al., 2019]
Ryo Karakida, Shotaro Akaho, and Shun ichi Amari.
Universal statistics of fisher information in deep neural networks:
Mean field approach, 2019.

[Jacot et al., 2019]
Arthur Jacot, Franck Gabriel, and Clément Hongler.
The asymptotic spectrum of the hessian of dnn throughout training.
arXiv preprint arXiv:1910.02875, 2019.

[Han et al., 2016]
Song Han, Huizi Mao, and William J. Dally.
Deep compression: Compressing deep neural networks with pruning,
trained quantization and huffman coding, 2016.

[Frankle and Carbin, 2019]
Jonathan Frankle and Michael Carbin.
The lottery ticket hypothesis: Finding sparse, trainable neural
arXiv: Learning, 2019.

[Maddox et al., 2020]
Wesley J. Maddox, Gregory Benton, and Andrew Gordon Wilson.
Rethinking parameter counting: Effective dimensionality revisted.
arXiv preprint arXiv:2003.02139, 2020.

[Belkin et al., 2019]
Mikhail Belkin, Daniel J. Hsu, Siyuan Ma, and Soumik Mandal.
Reconciling modern machine-learning practice and the classical
bias–variance trade-off.
Proceedings of the National Academy of Sciences, 116:0
15849 – 15854, 2019.

[Magnus and Neudecker, 2019]
Jan R Magnus and Heinz Neudecker.
Matrix differential calculus with applications in statistics
and econometrics.
John Wiley & Sons, 2019.

[Rao, 1972]
C Radhakrishna Rao.
Generalized inverse of a matrix and its applications.
In Vol. 1 Theory of Statistics, pages 601–620. University of
California Press, 1972.

[Chuai and Tian, 2004]
Jianjun Chuai and Yongge Tian.
Rank equalities and inequalities for kronecker products of matrices
with applications.
Applied Mathematics and Computation, 1500
(1):0 129–137, 2004.
ISSN 0096-3003.

[Bernacchia et al., 2018]
Alberto Bernacchia, Mate Lengyel, and Guillaume Hennequin.
Exact natural gradient in deep linear networks and its application to
the nonlinear case.
In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi,
and R. Garnett, editors, Advances in Neural Information Processing
Systems, volume 31. Curran Associates, Inc., 2018.

[Jacot et al., 2018]
Arthur Jacot, Franck Gabriel, and Clément Hongler.
Neural tangent kernel: Convergence and generalization in neural
arXiv preprint arXiv:1806.07572, 2018.

[LeCun and Cortes, 2010]
Yann LeCun and Corinna Cortes.
MNIST handwritten digit database.

[Xiao et al., 2017]
Han Xiao, Kashif Rasul, and Roland Vollgraf.
Fashion-mnist: a novel image dataset for benchmarking machine
learning algorithms, 2017.

[Krizhevsky, 2009]
Alex Krizhevsky.
Learning multiple layers of features from tiny images.

[Glorot and Bengio, 2010]
Xavier Glorot and Yoshua Bengio.
Understanding the difficulty of training deep feedforward neural
In Yee Whye Teh and Mike Titterington, editors, Proceedings of
the Thirteenth International Conference on Artificial Intelligence and
Statistics, volume 9 of Proceedings of Machine Learning Research,
pages 249–256, Chia Laguna Resort, Sardinia, Italy, 13–15 May 2010. PMLR.
URL <http://proceedings.mlr.press/v9/glorot10a.html>.

[Ubaru and Saad, 2016]
Shashanka Ubaru and Yousef Saad.
Fast methods for estimating the numerical rank of large matrices.
page 468–477, 2016.

[Golub et al., 1976]
Gene H. Golub, Virginia C. Klema, and Gilbert W. Stewart.
Rank degeneracy and least squares problems.
Technical report, Stanford, CA, USA, 1976.

[Lacotte et al., 2021]
Jonathan Lacotte, Yifei Wang, and Mert Pilanci.
Adaptive newton sketch: Linear-time optimization with quadratic
convergence and effective hessian dimensionality, 2021.

[Shen, 2016]
Hao Shen.
Towards a mathematical understanding of the difficulty in learning
with feedforward neural networks, 2016.

[Arora et al., 2019]
Sanjeev Arora, Nadav Cohen, Noah Golowich, and Wei Hu.
A convergence analysis of gradient descent for deep linear neural
networks, 2019.

[Pennington and Worah, 2018]
Jeffrey Pennington and Pratik Worah.
The spectrum of the fisher information matrix of a
single-hidden-layer neural network.
In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi,
and R. Garnett, editors, Advances in Neural Information Processing
Systems, volume 31. Curran Associates, Inc., 2018.

[Poggio et al., 2018]
Tomaso Poggio, Kenji Kawaguchi, Qianli Liao, Brando Miranda, Lorenzo Rosasco,
Xavier Boix, Jack Hidary, and Hrushikesh Mhaskar.
Theory of deep learning iii: explaining the non-overfitting puzzle,

[Neyshabur et al., 2015]
Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.
Norm-based capacity control in neural networks.
In Peter Grünwald, Elad Hazan, and Satyen Kale, editors,
Proceedings of The 28th Conference on Learning Theory, volume 40 of
Proceedings of Machine Learning Research, pages 1376–1401, Paris,
France, 03–06 Jul 2015. PMLR.
URL <http://proceedings.mlr.press/v40/Neyshabur15.html>.

[Bartlett et al., 2017]
Peter Bartlett, Dylan J. Foster, and Matus Telgarsky.
Spectrally-normalized margin bounds for neural networks, 2017.

[Arora et al., 2018]
Sanjeev Arora, Rong Ge, Behnam Neyshabur, and Yi Zhang.
Stronger generalization bounds for deep nets via a compression
approach, 2018a.

[Keskar et al., 2016]
Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy,
and Ping Tak Peter Tang.
On large-batch training for deep learning: Generalization gap and
sharp minima.
arXiv preprint arXiv:1609.04836, 2016.

[Jastrzębski et al., 2018]
Stanisław Jastrzębski, Zachary Kenton, Devansh Arpit, Nicolas Ballas, Asja
Fischer, Yoshua Bengio, and Amos Storkey.
Finding flatter minima with sgd, 2018.
URL <https://openreview.net/forum?id=r1VF9dCUG>.

[Baldassi et al., 2020]
Carlo Baldassi, Fabrizio Pittorino, and Riccardo Zecchina.
Shaping the learning landscape in neural networks around wide flat
Proceedings of the National Academy of Sciences, 1170
(1):0 161–170, 2020.

[Seroussi and Zeitouni, 2021]
Inbar Seroussi and Ofer Zeitouni.
Lower bounds on the generalization error of nonlinear learning
arXiv preprint arXiv:2103.14723, 2021.

[Matsaglia and Styan, 1974]
George Matsaglia and George P. H. Styan.
Equalities and inequalities for ranks of matrices.
Linear and Multilinear Algebra, 20 (3):0
269–292, 1974.
URL <https://doi.org/10.1080/03081087408817070>.

[Jain, 1989]
Anil K. Jain.
Fundamentals of Digital Image Processing.
Prentice-Hall, Inc., USA, 1989.
ISBN 0133361659.

[Arora et al., 2018]
Sanjeev Arora, Nadav Cohen, and Elad Hazan.
On the optimization of deep networks: Implicit acceleration by
International Conference on Machine Learning (ICML),

[Magnus and Neudecker, 1999]
Jan R. Magnus and Heinz Neudecker.
Matrix Differential Calculus with Applications in Statistics
and Econometrics.
John Wiley, second edition, 1999.
ISBN 0471986321 9780471986324 047198633X 9780471986331.

Supplementary Material

§ BACKPROPAGATION IN MATRIX-DERIVATIVES FOR THE GENERAL CASE

In the general case, we can represent the gradient in analogy as to Eq. (<ref>),

\begin{align}
\jacobi{F}{\wM{k}}= \Jmap{L:k+1} \Lmap{k}\kro \Big[ \Fmap{k-1:1}(\x) \Big]^\top,\,\, \text{where}\,\, \Fmap{k-1:1} = \Fmap{k-1} \circ \cdots \circ \Fmap{1}\,.
\end{align}

And, further $\mb J$ denotes the Jacobian map across the indexed layers, which is itself a composition of elementary Jacobians.

\begin{align}
\Jmap{k} & = \Lmap{k} \wM{k}, \quad \Lmap{k} = \diag(\dot \sigma^{(k)}), \quad  \Jmap{k+1:L} = {\Jmap{L:k+1}}^\top \,.
\end{align}

The Jacobian maps depend on functions of the input as each $\Lmap{k}$ depends on the pre-activation $\wM{k}\, \x_{k-1}$. By the usual chain rule (backpropagation) one has for a linear DNN, at a sample $(\x, \y)$:

\begin{align}
\jacobi{\ell}{\wM k} := \jacobi{\ell_{\x, \y}}{\wM k}
& = \underbrace{\Big[ \wM{k+1:L}\pmb \delta_{\x, \y} \Big]}_{\text{backward } \in \Re^{M_k}} \;\cdot \; \underbrace{\Big[  \wM{k-1:1} \x\Big]^\top}_{\text{forward } \in \Re^{M_{k-1}}} %\\[2mm]
= \wM{k+1:L} [\wM{L:1}\x \x^\top - \y \x^\top] \wM{1:k-1}\,.
\label{eq:supp-dnn-grad-linear}
\end{align}

The gradient with regard to $\wM k$ is first order in $\wM k$ and second order in the other matrices. In the general case, we get

\begin{align}
\jacobi{\ell}{\wM k}
& = \Big[\Jmap{k+1:L}\pmb \delta_{\x, \y}\Big] \cdot \Big[ \Fmap{k-1:1}(\x) \Big]^\top
\label{eq:dnn-grad}
\end{align}

Clearly, the partial forward maps are non-linear, whereas the backward maps are linearized at an argument determined by the current input.

§.§ Equivalence with Gauss-Newton decomposition

Remember, $\btheta \in \mathbb{R}^p$ denotes the (vectorized) parameters of the neural network map $F$, which then feeds into the loss $\ell$. Then, the Hessian of the composition of $\ell$ and $F$ with respect to $\btheta$ (computed over a sample $(\x, \y)$, but we omit specifying it for brevity) is,

\[
\small
\nabla^2_{\btheta} (\ell \circ F) = \nabla_{\btheta} F^\top \; [\partial^2 \ell ]\; \nabla_{\btheta} F \;+\; \sum_{c=1}^{K} \, [\partial \ell]_{c}\; \nabla^2_{\btheta}\, F_c  \,.
\]

where $\partial \ell$ and $\partial^2 \ell$ are respectively the gradient and Hessian  of the loss $\ell$ with respect to the network function, $F$. Also, $\nabla_{\btheta} F \in \mathbb{R}^{K\times p}$ is the Jacobian map of the network function $F(\x)$ with respect to the parameters $\btheta$. Let us employ the shorthand $\nabla_{k} F$ to denote the Jacobian of the network function with respect to the weight matrix $\wM{k}$ in the numerator-layout style as mentioned earlier. Similarly, let $\nabla^2_{kl} \,F_c$ be the Hessian of $c$-th component of the network function with respect to weight matrices $\wM{k}, \wM{l}$. Hence, we obtain,
\[\partial^2 \ell = \Im_K\,, \quad \nabla_{\idx} F =  \wM{L:\idx+1} \,\otimes \, \x^\top \wM{1:\idx-1}\,,\partial^2 \ell = \Im_K\,,\]

\begin{align*}
\nabla_{\idx} F &=  \wM{L:\idx+1} \,\otimes \, \x^\top \wM{1:\idx-1}\,,\\[1mm]
\nabla_{k} F^\top  \mb I_{K}
\nabla_{l} F & =\wM{k+1:L} \wM{L:l+1} \, \otimes \, \wM{k-1:1} \x\x^\top \wM{1:l-1}\,.
\end{align*}

Now, the other term has a reduction (or contraction) with the components of the residual, i.e., $\partial \ell = \pmb \delta_{\x, \y}$. Importantly, we notice that it has a block-hollow structure since,

$\nabla^2_{\idx \idx}\, F_c  = \mb 0, \,\, \forall c \in [1\cdots K]$.

Thus all the diagonal blocks come from the first term in the Gauss-Newton decomposition. Now, comparing the above expressions with the Eqns. (<ref>, <ref>, <ref>), it is evident that the two approaches yield the same structure of the Hessian.

As a side-remark, note that in Eqns. (<ref>, <ref>, <ref>) we express the $kl$-th block as $\dfrac{\partial^2 \Loss}{\partial \wM l \partial \wM k}$
instead of $\dfrac{\partial^2 \Loss}{\partial \wM k \partial \wM l}$,
only to ensure consistent shape as per matrix derivative convention (but corresponding entries are ofcourse equal).

§ TOOLS FOR THE ANALYSIS

§.§ General notation

We employ the shorthand notation, $\wM{k:l}$, to refer to the matrix product chain $\wM k \cdots \wM l$, when $k>l$. When $k < l$, $\wM{k:l}$ will stand for the transposed product chain $\wMt{k} \cdots \wMt{l}$.
In the edge case $k=l$, this will imply either $\wM{k}$ or $\wMt{k}$ depending on the context. Although, we will make the notation explicit on occasions where it might not be evident.
Besides, we use the $\kro$ to denote the Kronecker product of two matrices, $\vect_c$ and $\vect_r$ to denote the column-wise and row-wise vectorization of matrices respectively. $\Im_{k}$ denotes the identity matrix of size $k$, while $\mathds{1}_{k}$ denotes the all ones vector of length $k$, $\mb 0_{k}$ denotes an all-zeros matrix of size $k$. The generalized inverse [Rao, 1972] [This is like a general version of pseudoinverse which only satisfies the first Moore-Penrose condition.] of a matrix $\Am$ is given by a matrix $\Am^{-}$ which obeys $\Am \Am^{-} \Am = \Am$. The notation $\Am^{\bullet i}$ or $\Am_{\bullet i}$ denotes the $i$-th column of the matrix $\Am$, while $\Am^{i \bullet}$ or $\Am_{i \bullet}$ refers to its $i$-th row. We place the column and row indices into the subscript or superscript depending on the context they are used.

§.§ Helper Lemmas

Let $\Am\in\mathbb{R}^{m\times n} $ and $\Bm\in\mathbb{R}^{p\times q}$. Then the row-partitioned matrix $\left[\begin{array}{c}\Im_q \otimes \Am \\ \Bm \otimes \Im_n\end{array}\right]$ has the rank,

\[
\rank\left[\begin{array}{c}\Im_q \otimes \Am \\ \Bm \otimes \Im_n\end{array}\right]=q \rank(\Am)+n \rank(\Bm)-\rank(\Am) \rank(\Bm)
\]

The proof relies on the following rank formula due to [Matsaglia and Styan, 1974], and are based on using the generalized inverse of a matrix. Besides, the following proof closely follows Chuai and Tian, 2004.

$$\rank\left[\begin{array}{c}\Am \\ \Cm\end{array}\right]=\rank(\Am)+\rank\left(\Cm-\Cm \Am^{-} \Am\right)=\rank(\Cm)+\rank\left(\Am-\Am \Cm^{-} \Cm\right)$$

Here, $ \Cm^{-}$ denotes the weak (generalized) inverse of $\Cm$, i.e., any solution such that $\Cm  \Cm^{-} \Cm = \Cm$. Then, we have that,

\begin{align*}
\rank\left[\begin{array}{c}\Im_q \otimes \Am \\ \Bm \otimes \Im_n\end{array}\right]&= \rank(\Im_q \otimes \Am)+ \rank\big(\left(\Bm \otimes \Im_n\right) - \left(\Bm \otimes \Im_n\right){(\Im_q \otimes \Am)}^{-} \left(\Im_q \otimes \Am\right) \big) \\
&= \rank(\Im_q) \rank(\Am)+ \rank\big(\left(\Bm \otimes \Im_n\right) - \left(\Bm \otimes \Im_n\right){(\Im_q \otimes \Am^{-})} \left(\Im_q \otimes \Am\right) \big) \\
&= q \rank(\Am)+ \rank\big(\left(\Bm \otimes \Im_n\right) - \left(\Bm \otimes  A^{-} A\right) \big)\\
&= q \rank(\Am)+ \rank\big(\Bm \otimes \left(\Im_n -  A^{-} A\right) \big) \\
&= q \rank(\Am)+ \rank(\Bm)\rank\left(\Im_n -  A^{-} A\right) \\
&\overset{(a)}{=} q \rank(\Am)+ \rank(\Bm)(n-\rank(\Am)) \\
&= q \rank(\Am)+ n\rank(\Bm) - \rank(\Am)\rank(\Bm)
\end{align*}

where, in (a) we have used that $\rank\left(\Im_n -  A^{-} A\right)= n - \rank(\Am)$, which follows from the fact that column space of the matrix $\Im_n -  A^{-} A$, satisfies $\mathcal{C}\left(\Im_n -  A^{-} A\right)  \subset \mathcal{N}(A)$ (because the null space $\mathcal{N}(A)$ is the set of column vectors $\alpha$ for which $A\alpha=0$ and since any $x = \left(\Im_n -  A^{-} A\right) y \, \implies Ax=0 $) which means $\rank\left(\Im_n- A^{-} A\right)\leq \operatorname{dim} \mathcal{N}(A)$ and $\rank\left(\Im_n -  A^{-} A\right) \geq \rank\left(\Im_n\right)- \rank\left(A^{-} A\right) = n - \rank\left( A\right) = \operatorname{dim}\mathcal{N}(A)$.

Let $\Am_1\in\mathbb{R}^{m_1\times n_1} $, $\Am_2\in\mathbb{R}^{m_2\times n_2}$, $\Bm\in\mathbb{R}^{p\times n_1}$. Then the column block matrix $\left[\begin{array}{c} \Am_1\otimes \Am_2 \\ \Bm \otimes \Im_{n_2}\end{array}\right]$ has the rank,

\rank\left[\begin{array}{c}\Am_1\otimes \Am_2 \\ \Bm \otimes \Im_{n_2}\end{array}\right]= \rank(\Am_2)
\left(
\rank\left[\begin{array}{c}\Am_1 \\ \Bm\end{array}\right] - \rank(\Bm)
\right)
+ n_2 \, \rank(\Bm)

\begin{align*}
\rank\left[\begin{array}{c}\Am_1\otimes \Am_2 \\ \Bm \otimes I_{n_2}\end{array}\right] &\overset{(a)}{=} \rank(\Bm \kro \Im_{n_2}) + \rank\left((\Am_1\otimes \Am_2) - (\Am_1\otimes \Am_2) (\Bm \kro \Im_{n_2})^{-}(\Bm \kro \Im_{n_2})\right) \\
&= n_2 \, \rank(\Bm) + \rank\left((\Am_1\otimes \Am_2) - (\Am_1\otimes \Am_2) (\Bm^{-}\Bm \kro \Im_{n_2})\right)\\
&= n_2 \, \rank(\Bm) + \rank\left((\Am_1 - \Am_1 \Bm^{-}\Bm)\otimes \Am_2)\right)\\
&= n_2 \, \rank(\Bm) + \rank(\Am_1 - \Am_1 \Bm^{-}\Bm)\rank(\Am_2) \\
&\overset{(b)}{=} n_2 \, \rank(\Bm) +\left( \rank\left[\begin{array}{c}\Am_1 \\ \Bm \end{array}\right] - \rank(\Bm)\right)\rank(\Am_2) \\[1mm]
\left(
\rank\left[\begin{array}{c}\Am_1 \\ \Bm\end{array}\right] - \rank(\Bm)
\right)
+ n_2 \, \rank(\Bm)
\end{align*}

Step (a) and (b) are due to the rank formula for the column block matrix [Matsaglia and Styan, 1974].

$$\rank\left[\begin{array}{c}\Xm \\ \Ym\end{array}\right]=\rank(\Ym)+\rank\left(\Xm- \Xm\Ym^{-} \Ym\right)$$

And, the other steps follow from the basic properties of Kronecker product.

Let $\Am \in\mathbb{R}^{p\times q} $, $\Bm\in\mathbb{R}^{m\times p}$, then we have that:

\rank\left[\begin{array}{c} \Am\\ \Bm\Am\end{array}\right]=\rank(\Am)

This follows simply from the definition of generalized inverse ($\Am \Am^{-} \Am = \Am$) and the rank formula for column-block matrix.

\begin{align*}
\rank\left[\begin{array}{c}\Am \\ \Bm\Am\end{array}\right]&=\rank(\Am)+\rank\left((\Bm\Am)- (\Bm\Am) \Am^{-} \Am\right)\\
&=\rank(\Am)+\rank(\Bm\Am- \Bm\Am) = \rank(\Am)
\end{align*}

Let $\Bm\in\mathbb{R}^{m\times p}$, then we get the elementary identity:

\rank\left[\begin{array}{c} \Im_p\\ \Bm\end{array}\right]=\rank(\Im_p) = p\,.

Given a matrix $\mb M:= \mb A \mb B \mb A^T$, with $\Bm \succ 0$ symmetric, then $\rank(\Mm) = \rank(\Am)$.

Since $\Bm\succ 0$, we can write $\Mm = (\Am \Bm^{\frac{1}{2}}) (\Am \Bm^{\frac{1}{2}})^\top$, where $\Bm^{\frac{1}{2}}$ denotes the matrix square root of $\Bm$. This implies $\rank(\Mm) = \rank(\Am\Bm^{\frac{1}{2}})$, since the null space of $\Xm^\top \Xm$ is the same as the null space of any arbitrary matrix $\Xm$, i.e.,  $\mathcal{N}(\Xm^\top\Xm) =\mathcal{N}(\Xm)$, and additionally using $\rank(\Xm^\top)=\rank(\Xm)$. Next, as $\Bm \succ 0$, we have $\Bm^{\frac{1}{2}}\succ 0$, which further implies $\Bm^{\frac{1}{2}}$ is full rank. Hence $\rank(\Am\Bm^{\frac{1}{2}}) = \rank(\Am)$, which at last gives, $\rank(\Mm) = \rank(\Am)$.

§.§ {Left, Right, Pseudo}- Inverses

For a matrix $\Am$ with full column rank, the left inverse is defined to be a matrix $\Am^{-L}$ such that $\Am^{-L} \Am = \Im$. Likewise, when the matrix $\Am$ has full rank, we can define a right inverse which is a matrix $\Am^{-R}$ such that $\Am \Am^{-R} = \Im$. The left and right inverses need not be unique. But often a nice or convenient choice for the left inverse is $\Am^{-L} = {(\Am^\top \Am)}^{-1} \Am^\top$, while that for the right inverse is $\Am^{-R} = \Am^\top {(\Am \Am^\top)}^{-1}$.

The (Moore-Penrose) pseudoinverse $\Am^{\dagger}$ of a matrix $\Am \in \mathbb{R}^{m\times n}$ is a unique matrix that satisfies the following properties:

\begin{align*}
\Am \Am^\dagger \Am &= \Am \\
\Am^\dagger \Am \Am^\dagger &= \Am^\dagger \\
{(\Am \Am^\dagger)}^\top &= \Am \Am^\dagger \\
{(\Am^\dagger \Am)}^\top &= \Am^\dagger \Am  \\
\end{align*}

When the matrix has full column rank or full row rank, then the pseudoinverse agrees with the particular choice of left and right inverse we mentioned above. In such scenarios of full column or row rank, when we want to refer to this choice of left or right inverse, we will simply denote them by $\Am^{\dagger}$.

§.§ Block row and column operations

In our proofs, we make use of row and column operations, taken altogether on blocks of matrices rather than just individual rows or columns. So, here we clarify what we actually mean by such block row and column operations, and how their usage does not affect the rank of matrix to which these are applied. Essentially, we will look at the corresponding “elementary matrices” that get formed and argue that they multiplying with them does not change rank.

Let us consider that we have with us the following matrix, $\Pm$, with its row blocks labeled as $R_1, \cdots, R_L$.

\begin{equation}
\Pm =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
\\[2mm]
R_1& \Om^1\kro\Nm^1   \\[1mm]
\Vdots & \vdots \\[1mm]
R_{i} & \Om^i\kro\Nm^i  \\[2mm]
\Vdots & \vdots \\[1mm]
R_{L}  &   \Om^L\kro\Nm^L  \\
\end{pNiceArray}.
\end{equation}

§.§.§ Factoring-out block operations

Assume that the matrix $\Om^i$ has full column rank and then $\Om^i\kro \Im$ also has full column rank and is left-invertible. That means we can write the following decomposition $\Pm = \Pm_1 \Pm_2$.

\begin{equation}
\Pm = \underbrace{
\begin{pNiceArray}[first-row,nullify-dots]{ccccc}
& \\[2mm]
\Im & \cdots & \bm 0 & \cdots& \bm 0   \\[1mm]
\vdots &\ddots &\vdots & \ddots& \vdots\\[1mm]
\bm 0 & \cdots& \Om^i \kro \Im& \cdots & \bm 0 \\[2mm]
\vdots & \ddots & \vdots& \ddots & \vdots\\[1mm]
\bm 0 &\cdots & \bm 0 &\cdots & \Im \\
\end{pNiceArray}}_{\Pm_1}
\underbrace{
\begin{pNiceArray}[first-row,nullify-dots]{c}
\\[2mm]
\Om^1\kro\Nm^1  \\[1mm]
\vdots \\[1mm]
\Im \kro\Nm^i \\[2mm]
\vdots \\[1mm]
\Om^L\kro\Nm^L   \\
\end{pNiceArray}}_{\Pm_2}\,.
\label{eq:blockop-example}
\end{equation}

Since even $\Pm_1$ is left-invertible (identity matrix in all diagonals except for $\Om^i\kro\Im$ but which is left-invertible), we have that,

\[\rank(\Pm) = \rank(\Pm_1 \Pm_2) = \rank(\Pm_2)\,.\]

This is what we actually mean when applying the row operation,
\[R_i \leftarrow {(\Om^i\kro \Im)}^{\,\dagger} R_i\,.\]

Subsequently, we start working with the matrix $\Pm_2$, although we may not explicitly update the name of the matrix.

Remark. Likewise, we could have instead assumed $N^i$ to be left invertible, and factored out a corresponding $\Pm_1$ matrix with $\Im\kro N^i$ as one of its diagonal blocks. Going further, in a similar manner, we can define the column operation analogue of this by factoring out from the right a matrix which is right-invertible, and will thus preserve rank.

§.§.§ Deletion block operations

The block operations in this section are nothing but the analogue of the usual row operations, except they are carried out at the level of blocks.
Say that we have done the above factoring-out operation in Eq. (<ref>). Now, we are given that $N^1=N^i$, i.e, we have the matrix:

\begin{equation}
\Pm =
\begin{pNiceArray}[first-row,nullify-dots]{c}
\\[2mm]
\Om^1\kro\Nm^1  \\[1mm]
\vdots \\[1mm]
\Im \kro\Nm^1 \\[2mm]
\vdots \\[1mm]
\Om^L\kro\Nm^L   \\
\end{pNiceArray}\,.
\end{equation}

Now, consider we multiply from the left with the matrix $\Qm$ to yield $\Pm^\prime$:

\begin{equation}
\Pm^\prime = \underbrace{
\begin{pNiceArray}[first-row,nullify-dots]{ccccc}
& \\[2mm]
\Im & \cdots & - \Om^1 \kro \Im & \cdots& \bm 0   \\[1mm]
\vdots &\ddots &\vdots & \ddots& \vdots\\[1mm]
\bm 0 & \cdots&  \Im& \cdots & \bm 0 \\[2mm]
\vdots & \ddots & \vdots& \ddots & \vdots\\[1mm]
\bm 0 &\cdots & \bm 0 &\cdots & \Im \\
\end{pNiceArray}}_{\Qm}
\begin{pNiceArray}[first-row,nullify-dots]{c}
\\[2mm]
\Om^1\kro\Nm^1  \\[1mm]
\vdots \\[1mm]
\Im \kro\Nm^1 \\[2mm]
\vdots \\[1mm]
\Om^L\kro\Nm^L   \\
\end{pNiceArray} \quad = \quad
\begin{pNiceArray}[first-row,nullify-dots]{c}
\\[2mm]
\bm 0 \\[1mm]
\vdots \\[1mm]
\Im \kro\Nm^1 \\[2mm]
\vdots \\[1mm]
\Om^L\kro\Nm^L   \\
\end{pNiceArray}\,.
\label{eq:blockop-example}
\end{equation}

Since $\Qm$ is a upper-triangular matrix with ones on the diagonal, it is invertible. As a result, rank of $\Pm$ does not change when multiplied by $\Qm$ from the left.

This is what we actually mean when applying the row operation,
\[R_1 \leftarrow R_1 - {(\Om^1\kro \Im)} R_i\,.\]

Subsequently, we start working with the matrix $\Pm^\prime$, although we may not explicitly update the name of the matrix.

Remark. Notice, we could have done a similar thing on other side of Kronecker factors. Going further, in a similar manner, we can define the column operation analogue of this by multiplying on the right such a matrix which is invertible, and which will thus preserve rank.

Aliter. One can perhaps intuit these block operations from the point of view of inclusion of subspaces, but here we are being a bit pedantic.

§.§ Matrix derivatives

Let us start by discussing some simple facts on vectorization. Consider a matrix $\Am \in \mathbb{R}^{m\times n}$, and recall that $\vect_r$ and $\vect_c$ denote rowwise and columnwise vectorization respectively. Then firstly we have the following simple relation between them:

\begin{equation}\label{eq:vects}
\vect_c(\Am) = \vect_r(\Am^{\top})\,.
\end{equation}

Now, we give the proof of the commonly-used identity, $\vect_c(\Am\Xm\Bm) = \left(\Bm^{\top} \kro \Am\right) \vect_c(\Xm)$, where $\Am \in \mathbb{R}^{m\times n}\,,\, \Xm \in \mathbb{R}^{n\times p}\,,\, \Bm \in \mathbb{R}^{p\times q}$. For more details on this, refer to  [Magnus and Neudecker, 2019, Jain, 1989].

The main idea is to write $\Xm$ in the form of canonical basis vectors $\e_i$, i.e., $\Xm = \sum_{i=1}^{p} \Xm_{\bullet i} \,\e_{i}^{\top}$, where $\Xm_{\bullet i}$ denotes the $i$-th column of $\Xm$.

\begin{align*}
\vect_c\left(\Am  \sum_{i=1}^{p} \Xm_{\bullet i} \,\e_{i}^{\top} \Bm \right) &= \vect_c\left(\sum_{i=1}^{p} \left(\Am \Xm_{\bullet i}\right) \,\left(\Bm^{\top} \e_{i}\right) \right)\\
&\overset{(a)}{=} \sum_{i=1}^{p} \left(\Bm^{\top} \e_{i}\right) \kro \left(\Am \Xm_{\bullet i}\right) \\
&\overset{(b)}{=} \left(\Bm^{\top} \kro \Am\right) \, \sum_{i=1}^p \e_i \kro \Xm_{\bullet i} = \left(\Bm^{\top} \kro \Am\right) \vect_c(\Xm) \quad \QEDB
\end{align*}

In step (a), we have used that for two vectors $\ab, \bb \,$, the following basic fact $\vect_c(\ab \bb^{\top}) = \bb \kro \ab$ holds. And, in (b) we have employed the mixed-product property of Kronecker products.

Since, we utilize row-wise vectorization in our paper, let use find the equivalent relation in terms of that:

\begin{align}\label{eq:axb-kro}
\vect_r\left(\Am\Xm \Bm\right) \overset{\textrm{Eq.~\eqref{eq:vects}}}{=} \vect_c\left(\Bm^{\top} \Xm^\top \Am^\top\right) = \left(\Am\kro \Bm^\top\right) \vect_c\left(\Xm^\top\right) \overset{\textrm{Eq.~\eqref{eq:vects}}}{=} \left(\Am\kro \Bm^\top\right) \vect_r(\Xm)\,.
\end{align}

Recall, we use the numerator (Jacobian) layout to express matrix-by-matrix derivatives, i.e.,
\[\dfrac{\partial \mb Y}{\partial \mb X} := \dfrac{\partial \vect_r(\mb Y)}{\partial \vect_r(\mb X)^\top}\,.\]

Thus, when $\mb Y = \Am \Xm \Bm$, we use the above property along the first identification theorem of vector calculus [Magnus and Neudecker, 2019] as mentioned below:

(first identification theorem):

\mathrm{d} f=A(\x) \,\mathrm{d} \x \Longleftrightarrow \frac{\partial f(\x)}{\partial \x^{\top}}=A(\x)\,,

where, $\mathrm{d}$ denotes the differential.

Finally, this yields that:

\[\dfrac{\partial \Am\Xm\Bm}{\partial \mb X} = \Am \kro \Bm^\top \,.\]

Remark. If we were using the column vectorization, we would have instead obtained $\Bm^\top \kro \Am$.

§.§ Rank of weight matrices at initialization

Here we study the rank of random matrices to understand how the weight matrices of a neural network at initialization influence the rank.

Consider a random matrix $\Wm \in \mathbb{R}^{m \times n}$ for $m,n \in \mathbb{N}$ where each entry is sampled i.i.d. w.r.t. to some continuous (i.e. not discrete) probability distribution $p$, i.e. $W_{ij} \sim p $. Then it holds that

\begin{equation*}
\operatorname{rank}(\Wm) = \operatorname{min}(m, n) \hspace{2mm} \text{a.s.}
\end{equation*}

Assume w.l.o.g. that $n\leq m$ (otherwise consider the transposed matrix) and enumerate the columns of
$\Wm$ as $\Wm_1, \dots, \Wm_{n} \in \mathbb{R}^{m}$.
We need to show the linear independence of $\{\Wm_1, \dots,
\Wm_n\}$ over $\mathbb{R}^{m}$.
Let us show this inductively.
First $\Wm_1 \not= \bm{0}$ with probability $1$ since $p$ is continuous. Consider now $\Wm_2$.
Conditioned on the previously sampled column, $\Wm_1$, due to the continuous nature of the distribution,
the probability of $\Wm_2$ being contained in the span of $\Wm_1$ is zero:

\begin{equation*}
\mathbb{P}\left(\Wm_{2} \in \operatorname{span}(\Wm_1 )\right) = \int_{\mathbb{R}}\underbrace{\mathbb{P}\left(\Wm_{2} \in \operatorname{span}(\Wm_1 )|\Wm_1\right)}_{=0}p(\Wm_1)d\Wm_1  = 0
\end{equation*}

Finally, consider $\Wm_m$ and condition on the previously sampled vectors $\Wm_1, \dots, \Wm_{m-1}$. By the induction hypotheses, they span an $m-1$-dimensional space. Again, due to the independence of $\Wm_m$ from the previous vectors and the fact that an $m-1$-dimensional subspace has Lebesgue measure $0$ in $\mathbb{R}^{m}$, it holds in a similar fashion that

\begin{equation}
\mathbb{P}\left(\Wm_{m} \in \operatorname{span}(\Wm_1, \dots, \Wm_{m-1})\right) = 0
\end{equation}

and the matrix hence has full rank.

It turns out that we can apply a similar argument for the case of the product of two random matrices:

Consider random matrices $\Vm \in \mathbb{R}^{m \times n}$ and $\Wm \in \mathbb{R}^{n \times k}$, both drawn with i.i.d. entries according to some continuous probability distribution $p$. Define $\bm{Z} = V W \in \mathbb{R}^{m \times k}$. Then it holds that

\begin{equation*}
\operatorname{rank}(\bm{Z}) = \operatorname{min}(m, n, k)
\end{equation*}

First, notice that by Lemma <ref>, both matrices have full rank, i.e. $\operatorname{rank}(\Vm) = \operatorname{min}(m,n)$ and $\operatorname{rank}(\Wm) = \operatorname{min}(n, k)$.
By standard linear algebra results (not involving the fact that we have random matrices), we get that for $n \leq m$, $\operatorname{rank}(\bm{Z}) = \operatorname{rank}(\Wm) = \operatorname{min}(m, k)$ and for $n \leq k$, $\operatorname{rank}(\bm{Z}) = \operatorname{rank}(\Vm) = \operatorname{min}(m, n)$.

Thus it remains to show the case where $n \geq k, m$, i.e. the contracting dimension is the biggest. Assume w.l.o.g. that $k \leq m$ (otherwise study the transposed matrix). The columns of $\bm{Z}$ are given by $\bm{z}_i = \Vm\Wm_i$ for $i=1, \dots, k$. It thus suffices to show the linear independence of  $\{\bm{z}_1, \dots, \bm{z}_k\}$.

Notice that $\operatorname{rank}(\Vm) = m$ from the assumptions, thus $\{\Vm\x: \x \in \mathbb{R}^{n}\}$ is a $m$-dimensional subspace of $\mathbb{R}^{n}$. We will apply  a similar argument as in Lemma <ref>. Consider $\bm{z}_1 = \Vm\Wm_1$. Due to the independence and the fact that $\operatorname{rank}(\Vm)=m>0$, $\bm{z}_1 \not= 0$ a.s. Assume that $\bm{z}_1, \dots, \bm{z}_{k-1}$ are linearly independent, they thus span a $k-1$ dimensional space.

Conditioned on $\bm{z}_1, , \dots, \bm{z}_{k-1}$, $\bm{z}_k = \Vm\Wm_k$ is a random vector in $\operatorname{Im}(\Vm)$ (think of $\Vm$ as a fixed linear map since we are conditioning on it). Since $\operatorname{span}(\bm{z}_1, \dots, \bm{z}_{m-1})$ forms a $k-1$ dimensional subspace of $\operatorname{Im}(\Vm)$, which has dimension $m$, the subspace again has Lebesgue measure zero and we conclude that

\begin{equation*}
\mathbb{P}\left(\bm{z}_k \in \operatorname{span}(\bm{z}_1, \dots, \bm{z}_{k-1})\right) = \int_{\mathbb{R}}\underbrace{\mathbb{P}\left(\bm{z}_k \in \operatorname{span}(\bm{z}_1, \dots, \bm{z}_{k-1})|\bm{z}_1, \dots, \bm{z}_{k-1}\right)}_{=0} p(\bm{z}_k)d\bm{z}_k
\end{equation*}

We can now easily use this result for an arbitrary sized matrix product:

Consider random matrices $\Wm^{i} \in \mathbb{R}^{m_i \times m_{i+1}}$ for $i=1, \dots, n$ where each entry is initialized i.i.d. w.r.t. a continuous distribution $p$. Define the product matrix $\Wm = \Wm^1\dots \Wm^{n}$. Then it holds that

\begin{equation*}
\operatorname{rank}(\Wm) = \operatorname{min}(m_1, \dots, m_n)
\end{equation*}

Apply Lemma <ref> recursively, i.e. for $\Wm = \Wm^{1}\Wm^{2:n}$, then for $\Wm^{2:n} = \Wm^{2}\Wm^{3:n}$ up until $\Wm^{n}$.

§ RANK OF THE OUTER-PRODUCT TERM

We begin by discussing the proof of Proposition <ref>. Then, we briefly discuss the example of two-layer networks to motivate the proof, and after that we discuss the proof of the Theorem <ref>. Subsequently, we present the proof of Corollary <ref>.

§.§ Proof of Proposition <ref>

For a deep linear network, $\,\HO = \Am_o  \Bm_o {\Am_o}^\top\,$, where $\,\Bm_o = \Im_K \kro \covx\, \in \mathbb{R}^{Kd\times Kd}$,

\begin{equation*}
\text{and}\,\,\,\Am_o^\top=\begin{pmatrix}
\wM{L:2} \kro \Im_d \quad \cdots \quad
\wM{L:l+1} \kro \wM{1:l-1}
\quad \cdots \quad
\Im_K \kro \wM{1:L-1}
\end{pmatrix}\,\, \in \mathbb{R}^{Kd\times p}\,,
\end{equation*}

From Eq. (<ref>) we can notice that any block, say $kl$-th, of $\HO$ can be re-written as,

\begin{align*}
\HOsup{kl} =\left(\wM{k+1:L} \kro \wM{k-1:1} \right) \; \left(\Im_K \kro \covx\right) \; {\left(\wM{L:l+1} \otimes \wM{1:l-1}  \right)}\,,
\end{align*}

where, we have used the mixed-product property of Kronecker products i.e., $\Am\Bm\kro\Cm\Dm=(\Am\kro\Cm)(\Bm\kro\Dm)$. Now, it is clear from looking at the terms which are on the left and right of $\Im_K \kro \covx$, that we get the required decomposition.

As mentioned in the preliminaries, we consider, without loss of generality, that when the (uncentered) input covariance $\covx$ has rank $r<d$, then we take it to be

\[\covx = \begin{pmatrix}
(\covx)_{\,r\times r}& \mb{0}_{\,r\times (d-r)} \\[2mm]
\mb{0}_{\,(d-r)\times r}& \mb{0}_{\,(d-r)\times (d-r)}
\end{pmatrix}\,.\]

which is always possible by pre-processing the input (although this is not needed in practice). Thus, in such a scenario we can equivalently work with $\covx := (\covx)_{\,r\times r}$ and just consider the first $r$ columns of $\wM{1}$ (whose shape will then be $\mathbb{R}^{M_1\times r}$). Otherwise, if $r=d$, then we just continue with $\covx$ and $\wM{1}\in \mathbb{R}^{M_1\times d}$ as usual. To simplify our discussion ahead, we will always write $r$ in place of $d$, however the meaning of it should be clear from this remark.

Thus, in our presentation of the proof of Theorem <ref>, we will make a similar adaptation to Proposition <ref>, and so

\begin{equation*}
\Am_o^\top=\begin{pmatrix}
\wM{L:2} \kro \Im_r & \cdots \quad
\wM{L:l+1} \kro \wM{1:l-1}
\quad \cdots &
\Im_K \kro \wM{1:L-1}
\end{pmatrix}\,\, \in \mathbb{R}^{Kr\,\times \, p^\prime}\,,
\end{equation*}

where $p^{\prime} = p+ K(r-d)$.

§.§ Example for two-layer networks

Let us first illustrate Theorem <ref> via the example of $L=2$, i.e., we have a 2-layer network $\Fn(\x)= \wMb \wMa \x$, with weight matrices $\wMb \in \mathbb{R}^{K\times M_1}$ and $\wMa \in \mathbb{R}^{M_1\times r}$. Applying Proposition (<ref>), we obtain $\Am_o$ with the familiar structure:

\[\Am_o = \begin{pNiceArray}{c}
\mb W^{2^\top}  \kro \Im_r \\[1mm]
\Im_K \,\kro \wMa \\[1mm]
\end{pNiceArray}\]

Applying Lemma <ref> on $\Am_o\,$ thus yields:
$\rank(\Am_o) = r \, \rank(\mb W^{2^\top} ) \,+\, K \, \rank({\wMa})\, - \,\rank(\mb W^{2^\top} )  \rank({\wMa})$. If we assume that the hidden layer is the bottleneck, i.e., $q:=\min(r, M_1, K)=M_1$, then we get $\rank(\Am_o) = r \, M_1 + K \, M_1 - {M_1}^2$, keeping in mind the Assumption <ref>.

While here the special $\Zm$-like structure is apparent at the outset, the general case of $L$-layers is more involved and requires additional work to reduce to this structure, as illustrated in our proof ahead.

§.§ Proof of Theorem <ref>

Let us restate the Theorem <ref> from the main text,
Consider the matrix $\Am_o$ mentioned in Proposition <ref>.
Under the assumption <ref>,

\[
\rank(\Am_o) =
r \rank(\wM{2:L}) + K \rank(\wM{L-1:1}) - \rank(\wM{2:L})\rank(\wM{L-1:1}) = q \,( r + K  - q)\,.
\]

The proof is divided into two parts: (1) Bottleneck case and (2) Non-bottleneck case. For more details about the block-row operations that we employ here, please refer to the Section <ref>.

Part 1: Bottleneck case.
We assume without loss of generality that the layer $\ell-1$ has the minimum layer width out of all hidden layers, besides what is known that $M_{\ell-1} < \min(d, K)$. We will therefore have that $\Wm^{\ell}\in \mathbb{R}^{M_{\ell}\times M_{\ell-1}}$ will have full column rank and will be left-invertible. Due to random initialization of weight matrices (see Section <ref>) we will also have $\wM{k:\ell}, \, k\geq \ell$ to also have a left inverse.

Then let us write the block matrix $\Am_o$ in the column manner and label the row blocks corresponding to layer $\ell$ as $R_\ell$:
code-for-first-row = ,
code-for-first-col = ,

\begin{equation}\label{eq:startA}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\[2mm]
R_1& \wM{2:L}\kro  \Im_r   \\[1mm]
\Vdots & \vdots \\[1mm]
R_{\ell-1} & \wM{\ell:L} \kro  \wM{\ell-2:1} \\[2mm]
R_{\ell} & \wM{\ell+1:L} \kro  \wM{\ell-1:1} \\[2mm]
R_{\ell+1} & \wM{\ell+2:L} \kro  \wM{\ell:1} \\[2mm]
\Vdots & \vdots \\[1mm]
R_{L}  &   \mathbf{I}_K \kro \wM{L-1:1}  \\
\end{pNiceArray}.
\end{equation}

Consider the following (block) row operations:

R_{k} \leftarrow \left(\bm I_{M_{k}} \kro {\wM{k-1:\ell}}\right)^{\,\dagger} R_{k}, \quad \forall k \in {\ell +1, \cdots, L}

These row-operations are valid as the pre-factor has full-column rank, and are rank preserving, as discussed in Section <ref>. In this way, we have

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\[2mm]
R_1& \wM{2:L}\kro  \Im_r   \\[1mm]
\Vdots & \vdots \\[1mm]
R_{\ell-2} & \wM{\ell-1:L} \kro  \wM{\ell-3:1} \\[2mm]
R_{\ell-1} & \wM{\ell:L} \kro  \wM{\ell-2:1} \\[2mm]
R_{\ell} & \wM{\ell+1:L} \kro  \wM{\ell-1:1} \\[2mm]
R_{\ell+1} & \wM{\ell+2:L} \kro  \wM{\ell-1:1} \\[2mm]
\Vdots & \vdots \\[1mm]
R_{L}  &   \mathbf{I}_K \kro \wM{\ell-1:1}  \\
\end{pNiceArray}.
\end{equation}

Similarly, $\wM{{\ell-1}^\top} \in \mathbb{R}^{M_{\ell-2} \times M_{\ell-1}}$ as well as $\wM{k\,:\,\ell-1}$ for  $k\leq \ell-1$ is also full-column rank and thus left-invertible. Then apply the following row operations,

R_{k} \leftarrow \left(\wM{k +1\,:\, \ell-1} \kro I_{M_{k-1}}\right)^{\,\dagger} R_{k} , \quad \forall k \in {\ell -2, \cdots, 1}

And we get,

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}[code-before =\rowcolor{blue!15}{1,9}]
& \\[2mm]
R_1& \wM{\ell:L}\kro  \Im_r   \\[1mm]
\Vdots & \vdots \\[1mm]
R_{\ell-2} & \wM{\ell:L} \kro  \wM{\ell-3:1} \\[2mm]
R_{\ell-1} & \wM{\ell:L} \kro  \wM{\ell-2:1} \\[2mm]
\hline\\
R_{\ell} & \wM{\ell+1:L} \kro  \wM{\ell-1:1} \\[2mm]
R_{\ell+1} & \wM{\ell+2:L} \kro  \wM{\ell-1:1} \\[2mm]
\Vdots & \vdots \\[1mm]
R_{L}  &   \mathbf{I}_K \kro \wM{\ell-1:1}  \\
\end{pNiceArray}.
\end{equation}

Now using $R_1$, we can apply the deletion block operations to remove $\{R_2, \cdots, R_{\ell-1}\}$ as the left term in the Kronecker is identical. Next, via $R_{L}$ we can apply the deletion block operations to get rid of $\{R_\ell, \cdots, R_{L}\}$ as now the right term in the Kronecker is identical.  We are left with:

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\[2mm]
R_1& \wM{\ell:L}\kro  \Im_r   \\[2mm]
R_{L}  &   \mathbf{I}_K \kro \wM{\ell-1:1}  \\
\end{pNiceArray}.
\end{equation}

Finally, we can apply Lemma <ref> to obtain that, when $\Wm^{\ell}$ has full column rank or ${M_{\ell-1}}$ is the minimum layer-width (i.e., the bottleneck dimension):

\begin{align*}
\rank(\Am_o)
&= r \rank(\wM{\ell:L}) + K \rank(\wM{\ell-1:1}) - \rank(\wM{\ell:L})\rank(\wM{\ell-1:1}) \\
&= r {M_{\ell-1}} + K {M_{\ell-1}} - {M_{\ell-1}}^2 \\
&=  q \,( r + K  - q)\,,
\end{align*}

where, in the last step, we have used the definition of $q:=\min(r, M_1, \cdots, M_{L-1}, K)=M_{\ell-1}$ which gives rise to the equivalent expression.

Part 2: Non-bottleneck case.

This part is very similar and we will see that it uses just one set of row operations (like for the layers $> \ell$ and $< \ell$). In the non-bottleneck case, there are further two possibilities:

When $K$ is the minimum: This means that $\wMt{L}$ has full column-rank as $M_{L} = K$ is the minimum of all layer widths and input dimensionality. We start from the same $\Am_o$ matrix as in Eq. (<ref>). Now consider the following factoring-out operations:

R_{k} \leftarrow \left({\wM{k+1:L}} \kro \bm I_{M_{k-1}} \right)^{\,\dagger} R_{k}  \quad \forall k \in {1, \cdots, L-1}

These row-operations are valid (i.e., rank-preserving) as the pre-factor of (block-)row $R_k$ has full column rank, as discussed in Section <ref>. This results in,

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}[code-before =\rowcolor{blue!15}{1}]
& \\[2mm]
R_1&\mathbf{I}_K\kro  \Im_r   \\[1mm]
\Vdots & \vdots \\[1mm]
R_{\ell-1} & \mathbf{I}_K\kro  \wM{\ell-2:1} \\[2mm]
R_{\ell} & \mathbf{I}_K \kro  \wM{\ell-1:1} \\[2mm]
R_{\ell+1} & \mathbf{I}_K \kro  \wM{\ell:1} \\[2mm]
\Vdots & \vdots \\[1mm]
R_{L}  &   \mathbf{I}_K \kro \wM{L:1}  \\
\end{pNiceArray}.
\end{equation}

This is then followed by,

R_{k} \leftarrow R_{k} - \left( \bm I_{K}  \kro {\wM{k-1:1}} \right) R_{1}, \quad \forall \, k \in {2, \cdots, L}\,.

This results in,

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\
R_1& \Im_K \kro  \Im_r   \\[1mm]
\end{pNiceArray}.
\end{equation}

When $r$ is the minimum: This means that $\wM{1}$ has full column-rank as $M_{0} = r$ is the minimum of all layer widths, input and output dimensionality.  We start from the same $\Am_o$ matrix as in Eq. (<ref>). Now consider the following factoring-out operations:

R_{k} \leftarrow \left(\bm I_{M_{k}} \kro {\wM{k-1:1}}\right)^{\,\dagger} R_{k} \quad \forall \, k \in {2, \cdots, L}

As before, these row-operations are valid (i.e., rank-preserving) as the pre-factor of (block-)row $R_k$ has full column rank, as discussed in Section <ref>. This results in,

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}[code-before =\rowcolor{blue!15}{7}]
& \\[2mm]
R_1& \wM{2:L}\kro  \Im_r   \\[1mm]
\Vdots & \vdots \\[1mm]
R_{\ell-1} & \wM{\ell:L} \kro   \Im_r \\[2mm]
R_{\ell} & \wM{\ell+1:L} \kro   \Im_r \\[2mm]
R_{\ell+1} & \wM{\ell+2:L} \kro   \Im_r \\[2mm]
\Vdots & \vdots \\[1mm]
R_{L}  &   \mathbf{I}_K \kro  \Im_r  \\
\end{pNiceArray}.
\end{equation}

This is then followed by,

R_{k} \leftarrow  R_{k} - \left({\wM{k+1:L}}  \kro \bm I_{r} \right) R_{L}, \quad \forall k \in {1, \cdots, L-1}\,.

This results in,

\begin{equation}
\Am_o =
\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\
R_L& \Im_K \kro  \Im_r   \\[1mm]
\end{pNiceArray}.
\end{equation}

Resulting rank:  Thus, for either scenario of $K$ or $r$ being the minimum we get,
$$\rank(\Am_o) = \rank(\bm I_{K} \kro  \Im_r) = K r. $$

Final note. It is easy to check that for both the first and second case,  we can summarize the obtained rank in the form of the following equality:

$$\rank(\Am_o) = r \rank(\wM{L:2}) + K \rank(\wM{L-1:1}) - \rank(\wM{L:2})\rank(\wM{L-1:1}) =  q \,( r + K  - q).$$

§.§ Proof of Corollary <ref>

Let us first recall the Corollary,

Under the setup of Theorem <ref>, the rank of $\HO$ is given by
\[\rank(\HO) = q \,( r + K  - q)\,.\]

It is quite evident that we can use the rank of $\Am_o$ to bound the rank of $\HO$ due to the decomposition from Proposition <ref>. However, as mentioned in the main text, we can show an equality using the Lemma <ref> since $\HO$ is also of the form $\Am \Bm \Am^{\top}$ with $\Bm = (\Im_K \kro \covx) \succ \mb{0}$.

\begin{align*}
\rank(\HO) \overset{\textrm{Prop.~\ref{eq:outer-decomp}}}{=} \rank(\Am_o\Bm \Am_o^{\top}) = &\overset{\textrm{Lemma~\ref{lemma:outer-equality}}}{=} \rank(\Am_o)\\ &\overset{\textrm{Thm.~\ref{theorem:ub-outer}}}{=}
q \,( r + K  - q)\,.
\end{align*}

§ RANK OF THE FUNCTIONAL HESSIAN TERM

§.§ Proof of Theorem <ref>

Similar to the outer-product case, here as well the proof is divided into two cases. However, there is a subtle difference in that the input-residual covariance matrix ($\Om=\E[\pmb \delta_{\x, \y} \, \x^\top]$) can also dictate the rank, besides the weight matrix with minimum dimension. Especially since during training, as the residual approaches zero, the matrix $\Om\rightarrow0$. To abstract this, we update our definition of $q$ to include $s=\rank(\Om)$, and thus we get

\[q=\min(r, M_1, \cdots, M_{L-1}, K, s)\,.\]

Notice the additional $s$ as the last argument in the minimum. Further, since $\Om\in \mathbb{R}^{K\times r}$, we have that $s\leq\min(K, r)$.

Given these considerations, we split our analysis to the case where $q=s$ (referred to as the non-bottleneck case) and where $q=\min(M_1, \cdots, M_{L-1})$ (referred to as the bottleneck case). Note, if $q=r$ or $q=K$, then these cases are already subsumed by $q=s$ case, since $s\leq\min(K,r)$.

Further, in each of the two cases, we will analyse the rank of the block-columns of the functional Hessian formed with respect to a transposed weight matrix, i.e., $\HFhatsup{\bullet \ell}$. This makes it easier to deal with underlying structure of the actual block-columns $\HFcol{\ell}$, and does not affect the rank as it is invariant to row or column permutations. Finally, before we proceed into the details of the two parts, let us recollect the theorem statement,

For a deep linear network, the rank of $\ell$-th column-block, $\HFhatsup{\bullet \ell}$, of the matrix $\HFhat$, under the assumption <ref> is given as
$\rank(\HFhatsup{\bullet \ell}) = q\, M_{\ell -1} + q\,M_{\ell}  - q^2\,,\,$ for $\ell \in [2, \cdots , L-1]$.
When $\ell=1$, we have $
\rank(\HFhatsup{\bullet 1}) =q\, M_{1} + q \, s - q^2\,.$ And, when $\ell=L$, we have $
\rank(\HFhatsup{\bullet L}) =q\, M_{L-1} + q \, s - q^2\,.$
Here, $q := \min(r, M_1, \cdots, M_{L-1}, K, s)$ and $s:=\rank(\Om)=\rank(\E[\pmb \delta_{\x,\y}\,\x^\top])$.

§.§.§ Non-bottleneck case ($q = s$)

In this case, the Theorem boils down to showing the following:

* The rank of the $\ell^\text{th}$ column ($\ell \in [2, \cdots , L-1]$) of the functional Hessian, i.e. $\HFhatsup{\bullet \ell}$, is given by:

\[
\rank(\HFhatsup{\bullet \ell}) = s\, M_{\ell-1} + s\,M_{\ell}  - s^2\]

* The rank of the first column of the functional Hessian, i.e. $\HFhatsup{\bullet 1}$, is given by:

\[
\rank(\HFhatsup{\bullet 1}) = s\, M_{1}
\]

* The rank of the last column of the functional Hessian, i.e. $\HFhatsup{\bullet L}$, is given by:

\[
\rank(\HFhatsup{\bullet L}) = s\, M_{L-1}
\]

In this case, $\Om$, which shows up in every single block will dictate the rank. As before, for more details about the factoring-out and deletion block-row operations that we employ here, please refer to the Section <ref>. Let us now look at each of the parts:

Part 1. Let us first consider the case of the inner-columns, with $\ell>1$ and $\ell<L$. Now, the expression for $\HFhatsup{\bullet \ell}$ is given by,

code-for-first-row = ,
code-for-first-col = ,

\small
\HFhatsup{\bullet \ell}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wMt{\ell}}) \\[2mm]
\vect_{r}(\wM{1}) & \wM{2:\ell-1} \kro \Om^\top  \wM{L:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{j}) & \wM{j +1:\ell-1} \kro \wM{j-1:1} \Om^\top \wM{L:\ell+1}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{\ell-1}) & \Im_{M_{\ell-1}}\kro \wM{\ell-2:1} \Om^\top  \wM{L:\ell+1}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell}) & \bm{0}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell+1}) &  \wM{\ell+2:L} \Om \wM{1:\ell-1} \kro \Im_{M_{\ell}} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k}) & \wM{k+1:L} \Om \wM{1:\ell-1} \kro \wM{k-1:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L}) & \Om \wM{1:\ell-1} \kro \wM{L-1:\ell +1} \\[1mm]

\end{pNiceArray}

Now, since $\rank(\Om) = s$, we can express it as $\Om = \Cm \Dm$, where $\Cm \in \mathbb{R}^{K\times s}$ and $\Dm \in \mathbb{R}^{s \times d}$. Then consider the following factoring-out block operations:

R_{j} \leftarrow  \left(\Im_{m_j} \kro \wM{j-1:1}\Dm^\top \right)^\dagger R_{j}, \quad \forall j \in {1, \cdots, \ell-1}

where, $^\dagger$ denotes the left-inverse. When $j=1$, $\wM{j-1:1}:=\Im_r$. These row-operations are valid (i.e., rank-preserving) as $\wM{j-1}\Dm^\top$ has full column rank and is left-invertible. Similarly, consider the following factoring-out block operations:

R_{k} \leftarrow  \left(\wM{k+1:L} \Cm \kro \Im_{M_{k-1}} \right)^\dagger R_{k}, \quad \forall k \in {\ell+1, \cdots, L}

Note, when $k=L$, $\wM{k+1:L}:=\Im_K$. These row-operations are also valid (i.e., rank-preserving) as $\wM{k+1:L}\Cm$ has full column rank and is left-invertible. Then we can express the $\HFhatsup{\bullet \ell}$ as follows:

\small
\HFhatsup{\bullet \ell}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wMt{\ell}}) \\[2mm]
\vect_{r}(\wM{1}) & \wM{2:\ell-1} \kro \Cm^\top  \wM{L:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{j}) & \wM{j +1:\ell-1} \kro \Cm^\top \wM{L:\ell+1}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{\ell-1}) & \Im_{M_{\ell-1}}\kro \Cm^\top  \wM{L:\ell+1}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell}) & \bm{0}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell+1}) & \Dm \wM{1:\ell-1} \kro \Im_{M_{\ell}} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k}) & \Dm \wM{1:\ell-1} \kro \wM{k-1:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L}) & \Dm \wM{1:\ell-1} \kro \wM{L-1:\ell +1} \\[1mm]

\end{pNiceArray}

Now, it is clear that we can use row block $\ell-1$ to eliminate all row blocks prior to it and row block $\ell+1$ to eliminate all row blocks after it by the following deletion block operations:

R_{j} \leftarrow R_{j} - \left(\wM{j +1:\ell-1} \kro \Im_{q}\right) R_{\ell-1}, \quad \forall j \in {1, \cdots, \ell-1}

R_{k} \leftarrow R_{k} - \left(\Im_{q} \kro \wM{k-1:\ell+1}\right) R_{\ell+1}, \quad \forall k \in {\ell+1, \cdots, L}

Hence, we have that
\begin{align*}
\rank(\HFhatsup{\bullet \ell}) = \rank\left(\begin{array}{c}
\Im_{M_{\ell-1}}\kro \Cm^\top  \wM{L:\ell+1} \\[1mm]
\Dm \wM{1:\ell-1} \kro \Im_{M_{\ell}}
\end{array}\right)  =  s\, M_{\ell-1} + s\,M_{\ell}  - s^2 \end{align*}

where, we used the Lemma from [Chuai and Tian, 2004] in the last step.

Part 2. The procedure for this part will follow Part 1 procedure for blocks after the zero block.

\begin{align*}
\rank(\HFhatsup{\bullet \ell}) = \min(s\, M_1,  d\, M_1) = s\, M_1
\end{align*}

Part 3. The procedure for this part will follow Part 1 procedure for blocks before the zero block.

\begin{align*}
\rank(\HFhatsup{\bullet \ell})
=\min(s\, M_{L-1},  K\, M_{L-1}) = s\, M_{L-1}
\end{align*}

§.§.§ Bottleneck case ($q \neq s$)

Here, we need to prove the following:

* The rank of the $\ell^\text{th}$ column ($\ell \in [2, \cdots , L-1]$) of the functional Hessian, i.e. $\HFhatsup{\bullet \ell}$, is given by:

\rank(\HFhatsup{\bullet \ell}) =q\, M_{\ell-1} + q\,M_{\ell}  - q^2$$

* The rank of the first column of the functional Hessian, i.e. $\HFhatsup{\bullet 1}$, is given by:

\rank(\HFhatsup{\bullet 1}) = q\, M_{1} + q \, s - q^2

* The rank of the last column of the functional Hessian, i.e. $\HFhatsup{\bullet L}$, is given by:

\rank(\HFhatsup{\bullet L}) =q\, M_{L-1} + q \, s - q^2

Let us assume that $M_k$ is the bottleneck width, so it will “dictate" the rank now. Like in previous parts, for more details about the block-row operations that we employ here, please refer to the Section <ref>. Let us now look at each of the parts:

Part 1. Let us first consider the case of the inner-columns, with $\ell>1$ and $\ell<L$. Further, let us take $k>\ell$, and the procedure in the other scenario of $k<\ell\,$ is similar.

We have that $\HFhatsup{\bullet \ell}$ is as follows:

code-for-first-row = ,
code-for-first-col = ,

\small
\HFhatsup{\bullet \ell}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wMt{\ell}}) \\[2mm]
\vect_{r}(\wM{1}) & \wM{2:\ell-1} \kro \Om^\top  \wM{L:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{j}) & \wM{j +1:\ell-1} \kro \wM{j-1:1} \Om^\top \wM{L:\ell+1}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{\ell-1}) & \Im_{M_{\ell-1}}\kro \wM{\ell-2:1} \Om^\top  \wM{L:\ell+1}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell}) & \bm{0}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell+1}) &  \wM{\ell+2:L} \Om \wM{1:\ell-1} \kro \Im_{M_{\ell}} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k}) & \wM{k+1:L} \Om \wM{1:\ell-1} \kro \wM{k-1:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L}) & \Om \wM{1:\ell-1} \kro \wM{L-1:\ell +1} \\[1mm]

\end{pNiceArray}

Notice, we can write $\wM{i+1:L} = \wM{i+1:k} \wM{k+1:L}, \, \forall i \leq k -1 $. The matrix $\wM{i+1:k}$ is left-invertible as $M_k$ is the bottleneck width. Then, we consider the following factoring-out block operations (see Section <ref>) for the below zero part:

R_{i} \leftarrow  \left(\wM{i+1:k} \kro \Im_{M_{i-1}} \right)^\dagger R_{i}, \quad \forall i \in {\ell +1, \cdots, k-1}

where, $^\dagger$ denotes the left-inverse. Now, for a layer $i$ between $\{k+2, \cdots,L\}$, we notice that $\wM{i-1:\ell+1} = \wM{i-1:k+1} \wM{k:\ell+1}$ with $\wM{i-1:k+1}$ being left-invertible. Hence consider the factoring-out operations,

R_{i} \leftarrow  \left(\Im_{M_{i}} \kro \wM{i-1:k+1} \right)^\dagger R_{i}, \quad \forall i \in {k +2, \cdots, L}

Now, for the row-blocks in the part above zero, we have that $\wM{j-1:1}\Om^\top\wM{L:\ell+1} = \wM{j-1:1}\Om^\top\wM{L:k+1}\wM{k:\ell+1}, \, \forall j \leq \ell-1$. Notice, that $\wM{j-1:1}\Om^\top\wM{L:k+1}$ is left invertible due to the bottleneck $M_k$. Hence, consider the following block-row operations:

R_{j} \leftarrow  \left(\Im_{m_j}\kro \wM{j-1:1}\Om^\top\wM{L:k+1} \right)^\dagger R_{j}, \quad \forall j \in {1, \cdots, \ell-1}

Here, it does not matter to us what the actual value of $s$ is, since we know that $M_k$ is the minimum width which guarantees that the above expression containing $\Om^\top$ is left invertible.

Overall, we can thus express the $\HFhatsup{\bullet \ell}$ as follows:

\small
\HFhatsup{\bullet \ell}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wMt{\ell}}) \\[2mm]
\vect_{r}(\wM{1}) & \wM{2:\ell-1} \kro \wM{k:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{j}) & \wM{j +1:\ell-1} \kro \wM{k:\ell+1}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{\ell-1}) & \Im_{M_{\ell-1}}\kro \wM{k:\ell+1}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell}) & \bm{0}  \\[2mm]
\hline\\
\vect_{r}(\wM{\ell+1}) & \wM{k+1:L}\Om \wM{1:\ell-1} \kro \Im_{M_{\ell}} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k}) & \wM{k+1:L} \Om \wM{1:\ell-1} \kro \wM{k-1:\ell+1} \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L}) & \Om \wM{1:\ell-1} \kro \wM{k:\ell +1} \\[1mm]

\end{pNiceArray}

Now, it is clear that we can use row $\ell-1$ to eliminate all rows prior to it as well as the rows from $k+1$ to $L$. While the row $\ell+1$ can be used to eliminate all rows from $\ell+2$ until $k$. Hence, we have that,

\begin{align*}
\rank(\HFhatsup{\bullet \ell}) = \rank\left(\begin{array}{c}
\Im_{M_{\ell-1}}\kro \wM{k:\ell+1} \\[2mm]
\wM{k+1:L} \Om \wM{1:\ell-1} \kro \Im_{M_{\ell}}
\end{array}\right)  =  q\, M_{\ell-1} + q\,M_{\ell}  - q^2 \end{align*}

where, we used the Lemma from [Chuai and Tian, 2004] in the last step.

Part 2. Now we deal with the column block corresponding to first layer, which is:

\small
\HFhatsup{1}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wM{1}}^\top) \\[2mm]
\vect_{r}(\wM{1}) & \bm{0}  \\[2mm]
\vect_{r}(\wM{2}) &  \wM{3:L} \Om  \kro \Im_{M_1}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k-1}) &  \wM{k:L} \Om \kro \wM{k-2:2}  \\[2mm]
\vect_{r}(\wM{k}) & \wM{k+1:L} \Om \kro \wM{k-1:2} \\[2mm]
\vect_{r}(\wM{k+1}) &  \wM{k+2:L} \Om \kro \wM{k:2}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L}) & \Om \kro \wM{L-1:2} \\[1mm]

\end{pNiceArray}

Basically, we have to follow the same procedure for row blocks before $k$ and for the ones after $k$ as done in the part 1.
In other words, for the row blocks prior to $k$, we can write the $\wM{i+1:L} = \wM{i+1:k}\wM{k+1:L}$. Since $\wM{i+1:k}$ is left-invertible due to the bottleneck, we consider the following factoring-out operations:

R_{i} \leftarrow  \left(\wM{i+1:k} \kro \Im_{M_{i-1}} \right)^\dagger R_{i}, \quad \forall i \in {2, \cdots, k-1}

Now, for a layer $i$ between $\{k+2, \cdots,L\}$, we notice that $\wM{i-1:2} = \wM{i-1:k+1} \wM{k:2}$ with $\wM{i-1:k+1}$ being left-invertible. Hence the factoring-out operations will be,

R_{i} \leftarrow  \left(\Im_{M_{i}} \kro \wM{i-1:k+1} \right)^\dagger R_{i}, \quad \forall i \in {k +2, \cdots, L}

This resulting $\HFhatsup{1}$ is as follows:

\small
\HFhatsup{1}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wM{1}}^\top) \\[2mm]
\vect_{r}(\wM{1}) & \bm{0}  \\[2mm]
% \hline\\
\vect_{r}(\wM{2}) &  \wM{k+1:L} \Om  \kro \Im_{M_1}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k-1}) &  \wM{k+1:L} \Om \kro \wM{k-2:2}  \\[2mm]
\vect_{r}(\wM{k}) & \wM{k+1:L} \Om \kro \wM{k-1:2} \\[2mm]
\vect_{r}(\wM{k+1}) &  \wM{k+2:L} \Om \kro \wM{k:2}  \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L}) & \Om \kro \wM{k:2} \\[1mm]

\end{pNiceArray}

Now, we can use row $2$ to eliminate rows $3$ until $k$, and similarly we can use row $L$ to eliminate rows $k+1$ to $L-1$. Thus, we have that:

\begin{align*}
\rank(\HFhatsup{1}) =  \rank\left(\begin{array}{c}
\wM{k+1:L} \Om  \kro \Im_{M_1} \\ [1mm]
\Om \kro \wM{k:2}
\end{array}\right)  = q\, M_1 + q\, s - q^2\end{align*}

Different from the previous analysis, in the above matrix we did not have a “naked" Identity matrix on the left, so we could not directly use the Lemma <ref>. But, we used its generalized version contained in Lemma <ref> along with the Lemma <ref> to obtain the rank in the final step.

Part 3. Finally, we have the last column:

\small
\HFhatsup{\bullet L}= \begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \vect_{r}({\wM{L}}^\top) \\[2mm]
\vect_{r}(\wM{1}) & \wM{2:L-1} \kro \Om^\top \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{k-1}) & \wM{k:L-1} \kro \wM{k-2:1} \Om^\top  \\[1mm]
\vect_{r}(\wM{k}) & \wM{k +1:L-1}\kro \wM{k-1:1} \Om^\top   \\[2mm]
\vect_{r}(\wM{k+1}) & \wM{k +2:L-1}\kro \wM{k:1} \Om^\top   \\[1mm]
\Vdots & \vdots \\[1mm]
\vect_{r}(\wM{L-1}) & \Im_{M_{L-1}} \kro \wM{L-2:1} \Om^\top  \\[1mm]
\vect_{r}(\wM{L}) & \bm{0}  \\[2mm]
\end{pNiceArray}

We can follow a similar strategy as carried out in Part 2 to get,

\begin{align*}
\rank(\HFhatsup{\bullet L}) = \rank\left(\begin{array}{c}
\wM{k+1:L-1}\kro \Om^\top \\ [1mm]
\Im_{M_{L-1}} \kro \wM{k:1} \Om^\top
\end{array}\right)  = q\, M_{L-1} + q\, s - q^2\end{align*}

§.§ Proof of Corollary <ref>

Let us remember the Corollary from the main text,
Under the setup of Theorem <ref>, the rank of $\HF$ can be upper bounded as,

\[\rank(\HF) \leq 2\, q\, M  + 2 \,q\, s - L\, q^2\,, \quad \text{where} \quad M=\sum_{\ell=1}^{L-1} M_{\ell}\,.\]

By using the above Theorem <ref> and applying the fact that $\rank(\left[\Am \,\, \Bm\right])\leq \rank(\Am) + \rank(\Bm)$ on the column-blocks of $\HFhat$, we get the desired upper bound on the rank of entire functional Hessian.

\begin{align*}
\rank(\HF) &= \rank(\HFhat) \,\,\leq \,\, \sum\limits_{\ell=1}^L \rank(\HFhatsup{\bullet \ell}) \\
&= q\, M_1 + q\, s - q^2 + \sum\limits_{\ell=2}^{L-1} \left(q\, M_{\ell -1} + q\,M_{\ell}  - q^2\right) + q\, M_{L-1} + q\, s - q^2 \\
&= 2\, q\, \left(\sum\limits_{\ell=1}^{L-1}M_{\ell}\right)  + 2 \,q\, s - L\, q^2 \,\,=\,\, 2\, q\, M  + 2 \,q\, s - L\, q^2
\end{align*}

§ EVOLUTION OF RANK DURING TRAINING

Here we prove Lemma <ref>, stating that the rank of the individual weights remains invariant under gradient flow dynamics.
For sake of readability, let us restate Lemma <ref>:
For a deep linear network, consider the gradient flow dynamics $\dot{\Wm}^{l}_t = -\eta \nabla_{\Wm^{l}}\mathcal{L}_S(\bm{\btheta})\big{|}_{\bm{\btheta} \,=\, \bm{\btheta}_t}$. Assume: (a) Centered classes: $ \frac{1}{N} \sum_{i:y_{ic} = 1}^{N}\x_i= \bm{0}, \hspace{2mm} \forall\, c \in [1,\dots,K]$. (b) Balancedness at initialization: ${\Wm^{{k+1}^{\top} }_0}\Wm_0^{k+1} = \Wm^{k}_0\\Wm^{k^{\top}}_0$. (c) Square weight-matrices: $\Wm^{l} \in \mathbb{R}^{M \times M}\,,$ $\,\forall \,l\,$ and $\,K=d=M$. Then for all layers $l$,
$\,\, \rank(\Wm^{l}_t) = \rank(\Wm^{l}_0), \hspace{3mm} \,\forall \, t<\infty\,.$

To prove Lemma <ref>, we first need some helper lemmas. We start with a well-known result:

Under assumption b) in Lemma <ref> and gradient flow dynamics, it holds that

$$\left(\Wm^{k+1}_t\right)^{T}\Wm_t^{k+1} = \Wm^{k}_t\left(\Wm^{k}_t\right)^{T}$$

We refer to [Arora et al., 2018] (Theorem 1) for a proof.

This essentially guarantees that balancedness is preserved throughout training if we guarantee it at initialization. Let us assume that the weights $\Wm^{l}$ are full rank, which for the squared matrix case is equivalent to $\operatorname{det}\left(\Wm^{l}\right) \not = 0$. The Jacobi formula allows us to extend the dynamics of $\Wm^{l}$ to the determinant:

Given a dynamic matrix $t \mapsto \Am(t)$, the determinant follows:

$$\frac{d}{dt}\operatorname{det}\left(\Am(t)\right) = \operatorname{tr}\left(\operatorname{adj}\left(\Am(t)\right)\frac{d \Am(t)}{dt}\right)$$

where $\operatorname{adj}\left(\Am(t)\right)$ is the adjugate matrix satisfying
$$\Am(t) \operatorname{adj}\left(\Am(t)\right) =  \operatorname{adj}\left(\Am(t)\right)\Am(t)= \operatorname{det}(\Am(t)) \bm{I}$$

For a proof of this standard result we refer to Magnus and Neudecker, 1999 for instance.

We are now ready to proof the main claim:

Proof of Lemma <ref>:
Let us first simplify the right-hand side of the gradient flow equation using the assumptions, before applying Lemma <ref> :

\begin{equation*}
\begin{split}
\nabla_{\Wm^{l}}\mathcal{L}(\Wm)\big{|}_{\Wm = -\eta \Wm_t} &= -\eta\Wm_t^{l+1:L}\bm{\Omega}\Wm_t^{1:l-1} \\&= -\eta\Wm_t^{l+1:L}\Wm_t^{L:l+1}\Wm_t^{l}\Wm_t^{l-1 : 1}\bm{\Sigma}\Wm_t^{1:l-1} \\
&\stackrel{L1}{=} -\eta\left(\Wm_t^{l}\Wm_t^{l  T}\right)^{L-l}\Wm_t^{l}\Wm_t^{l-1 : 1}\bm{\Sigma}\Wm_t^{1:l-1}
\end{split}
\end{equation*}

Applying Lemma <ref>, gives

\begin{equation*}
\begin{split}
\frac{d}{dt}\operatorname{det}\left(\Wm^{l}_t\right) &= \operatorname{tr}\left(\operatorname{adj}\left(\Wm^{l}_t\right)\dot{\Wm}^{l}_t\right) \\
&= -\eta\operatorname{tr}\left(\operatorname{adj}\left(\Wm^{l}_t\right)\left(\Wm_t^{l}\Wm_t^{l  T}\right)^{L-l}\Wm_t^{l}\Wm_t^{l-1 : 1}\bm{\Sigma}\Wm_t^{1:l-1}\right) \\
&= -\eta\operatorname{tr}\left(\operatorname{det}\left(\Wm^{l}_t\right)\left(\Wm_t^{l T}\Wm_t^{l  }\right)^{L-l}\Wm_t^{l-1 : 1}\bm{\Sigma}\Wm_t^{1:l-1}\right) \\
&= -\eta\operatorname{det}\left(\Wm^{l }_t\right)\operatorname{tr}\left(\left(\Wm_t^{l T}\Wm_t^{l }\right)^{L-l}\Wm_t^{l-1 : 1}\bm{\Sigma}\Wm_t^{1:l-1}\right) \\
&= -\eta\operatorname{det}\left(\Wm^{l}_t\right) \operatorname{tr}(\Bm(t))
\end{split}
\end{equation*}

We can write the solution of this differential equation as

$$\operatorname{det}\left(\Wm^{l }_t\right) = e^{-\eta \int_{0}^{t} \operatorname{tr}(\Bm(s)) ds}\operatorname{det}\left(\Wm^{l}_0\right)$$

Of course, we have no idea how to solve the integral in the exponential, but the exact solution does not matter as it is always positive and thus not zero. Thus, as long as $\operatorname{det}\left(\Wm^{l}_0\right) \not= 0$, also $\operatorname{det}\left(\Wm^{l}_t\right) \not = 0$ holds, at least for a finite time horizon $t < \infty$.

We illustrated this theoretical finding through non-linear networks. To complete the picture we also show a linear network with its corresponding rank dynamics for weights and Hessians in Figure <ref>. We use a linear teacher setting, i.e. $\y = \Wm \x$ where $\x \sim \mathcal{N}(\bm{0}, \Im)$ and $\Wm \in \mathbb{R}^{d \times K}$, allowing the model thus to achieve a perfect training error. We train the model for $2000$ epochs with SGD. Again we observe how the weights remain full-rank throughout training and as a consequence, the rank of the Hessians remains constant initially. Due to the linear nature of the teacher, an exact training error of zero is achievable. As soon as the error reaches a certain threshold (around $10^{-30}$), the contribution of the functional Hessian starts to vanish as its eigenvalues, one-by-one become too small to count towards the rank. As a consequence, the rank of the loss Hessian also decreases as it is composed of both the functional and the outer Hessian. Also observe how the outer Hessian remains constant throughout the entire optimization, as it does not depend on error. As expected, at the end of training, the loss Hessian collapses onto the outer Hessian.

Evolution of Rank:  We display the rank of the weights $\wM{i}$ (left) and the rank of the Hessians (right) throughout training. We use a $3$ hidden-layer linear network trained on Gaussian data which is labeled by a linear teacher.

§ PESSIMISTIC BOUNDS IN THE NON-LINEAR CASE

§.§ Proof strategy for the $\HO$ upper bound

The key idea behind our upcoming proof strategy is to do the analysis à la super-position of “unit-networks”, i.e., networks with one-hidden neuron. This is because we can reformulate the network function $F_{\btheta}$ as sum of network functions of $M$ unit-networks, one per each hidden neuron. Mathematically,

\begin{equation*}
F_{\btheta}(\x) = \sum\limits_{i=1}^M F_{\theta_i}(\x)\,,\quad \text{where,}\quad F_{\theta_i}(\x) = \wMb_{\bullet \, i} \, \sigma(\mb W^{1^\top}_{i \,\bullet} \, \x).
\end{equation*}

where, $\wMa_{i \,\bullet} \in \mathbb{R}^{r}$ is the $i$-th row of $\wMa$ and $\wMb_{\bullet \, i}\in \mathbb{R}^{K}$ is the $i$-th column of $\wMb$.
To ease notation, we will henceforth use the shorthand: $\wMasupp:=\wMa$, and $\wMbsupp:=\wMb$. The column and row indexing follow the convention as before.

Our approach will be to analyze the rank of the outer-product Hessian $\HO$ in the case of empirical loss, or in other words with finitely many samples. Then, we will do a limiting argument to extend it to case of population loss. The benefit of such an approach is that we can bound the rank of $\HO$ via the rank of the Jacobian of the function, $\nabla_{\btheta} F_{\btheta}(\Xm)$ formed over the entire data matrix $\Xm\in \mathbb{R}^{d\times N}$. We will, in turn, bound the rank of such a Jacobian matrix $\nabla_{\btheta} F_{\btheta}(\Xm)$ via the ranks of the Jacobian matrices for the unit-network functions $\nabla_{\btheta} F_{\theta_i}(\Xm)$.

Similar to Remark <ref>, when the rank $r$ of the data matrix (or alternatively the empirical (uncentered) input covariance) is less than $d$, we will consider without loss of generality that $\Xm \in \mathbb{R}^{r\times N}$.

In other words, for the sake of analysis, we consider that the data matrix has been pre-processed to take into account rank lower than input dimension.

§.§ Lemmas on the structure and rank of Jacobian matrix of the unit-network

As a first step, the following Lemma <ref> shows the structure of the Jacobian matrix for the unit-network, formed over the entire data matrix.

Consider the unit-network $\Fni{i}(\x) = \wMbsupp_{\bullet \, i} \, \sigma(\wMasupp^\top_{i \,\bullet} \, \x)$ corresponding to i-th neuron, with the non-linearity $\sigma$ such that $\sigma(z) = \sigma^\prime(z) z$. Let the data-matrix be denoted by $\Xm \in \mathbb{R}^{r \times N}$. Further, let us denote the matrix which has on its diagonal — the activation derivatives $\sigma^\prime(\wMasupp_{i \,\bullet}^\top \, \x)$ over all the samples $\x$, and zero elsewhere, by $\bm \Lambda^i \in \mathbb{R}^{N\times N}$.  Then the Jacobian matrix $\nabla_{\btheta} \Fni{i}(\Xm)$ is given (in transposed form) by,

code-for-first-row = ,
code-for-first-col = ,

$$\nabla_{\btheta} F_{\theta_i}(\Xm)^\top = \hspace{-1cm}\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\
\Vdots &        \bm{0} \\[1mm]
\wMasupp_{i \,\bullet}     & \Xm \, \bm \Lambda^i \kro \wMbsupp^\top_{\bullet \,i}  \\[2mm]
\wMbsupp_{\bullet \,i}    &  \wMasupp_{i \,\bullet}^\top  \Xm \, \bm \Lambda^i \kro \Im_{K}  \\[1mm]
\Vdots &        \bm{0} \\
\end{pNiceArray}

Let us start simple by computing the gradient with respect to $k$-th component of the function, i.e. $F^k$, which comes out as follows:

\[
\frac{\partial F^k}{\partial w_{ij}} = \mathds{1}\lbrace k=i\rbrace \, {\sigma(\Vm\x)}_j
\]

\[
\frac{\partial F^k}{\partial v_{ij}} = w_{ki} \,\, \sigma^\prime(\wMasupp_{i \,\bullet}^\top\,\x)\,\, x_j
\]

where ${\sigma(\Vm\x)}_j = \sigma( \Vm_{j \,\bullet}^\top \,\x)$, with $\Vm_{j \,\bullet}$ being the $j$-th row of $\Vm$, in a column vector format as mentioned above. And, for example in case of ReLU, $\sigma^\prime(z) = \partial\, \sigma(z)/\partial z = \mathds{1}\lbrace z>0\rbrace$.

For a fixed sample $\x$, consider that the activation derivatives for all hidden-neurons are stored in a diagonal matrix $\bm\Lambda^{\x}\in\mathbb{R}^{M\times M}$, i.e., $\bm \Lambda^{\x}_{jj} =  \sigma^\prime( \Vm_{j \,\bullet}^\top \,  \x)$. Then we can rewrite $\sigma(\Vm \x) = \bm \Lambda^{\x} \, \Vm \x$ for all non-linearities that satisfy $\sigma(z) = \sigma^\prime(z) z$. So, we have that the Jacobian of the function with respect to the parameters comes out to be,

code-for-first-row = ,
code-for-first-col = ,

\[\nabla_{\bm\theta} F_{\btheta}(\x)^\top= \hspace{-1cm}\begin{pNiceArray}[first-row,first-col,nullify-dots]{ccccc}
& \bm{F}^{1} & \cdots &  \bm{F}^{k} &  \cdots &  \bm{F}^{K}\\[2mm]
\wMasupp_{1 \bullet}   & w_{11} \, \bm \Lambda^{\x}_{11} \,\, \x & \cdots &  w_{k1}  \, \bm \Lambda^{\x}_{11} \,\, \x & \cdots &  w_{K1}  \, \bm \Lambda^{\x}_{11} \,\, \x\\
\Vdots    & \vdots &  & \vdots &  & \vdots \\[1mm]
\wMasupp_{i \bullet}   & w_{1i} \,  \bm \Lambda^{\x}_{ii}  \,\, \x & \cdots &  w_{ki} \,   \bm \Lambda^{\x}_{ii}  \,\, \x & \cdots &  w_{Ki}  \,  \bm \Lambda^{\x}_{ii}  \,\, \x\\
\Vdots    & \vdots &  & \vdots &  & \vdots \\[1mm]
\wMasupp_{M \bullet}   & w_{1M} \,   \bm \Lambda^{\x}_{MM} \,\, \x & \cdots &  w_{kM} \,   \bm \Lambda^{\x}_{MM} \,\, \x& \cdots &  w_{KM}  \,  \bm \Lambda^{\x}_{MM} \,\, \x\\[2mm]
\hline
& & & &  \\[-1mm]
\wMbsupp_{1 \bullet}   &  \bm \Lambda^{\x} \, \Vm \x & \cdots & \bm{0} & \cdots &  \bm{0}\\
\Vdots    & \vdots & \ddots &  &  \ddots& \vdots \\[1mm]
\wMbsupp_{k \bullet}    &  \bm{0}&  &  \bm \Lambda^{\x} \, \Vm \x & &  \bm{0}\\
\Vdots    & \vdots &  \ddots &  &  \ddots & \vdots \\[1mm]
\wMbsupp_{K \bullet}   &  \bm{0}& \cdots & \bm{0} & \cdots &  \bm \Lambda^{\x} \, \Vm \x \\
\end{pNiceArray}\]

Hence, from above we can we write the Jacobian of the $i$-th unit-network function $\wMbsupp_{\bullet \, i} \, \sigma(\wMasupp^\top_{i \,\bullet} \, \x)$ with respect to the entire set of parameters and at a given input $\x$, as follows:
code-for-first-row = ,
code-for-first-col = ,

\[\nabla_{\btheta} \, F_{\theta_i}(\x)^\top= \hspace{-1cm}\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\
\Vdots &        \bm{0} \\[1mm]
\wMasupp_{i \,\bullet}   & \wMbsupp^\top_{\bullet \,i} \kro \bm \Lambda^{\x}_{ii}  \,\, \x\\[2mm]
\wMbsupp_{\bullet \,i}    &  \Im_{K}\kro  \wMasupp^\top_{i \,\bullet}  \big(\bm \Lambda^{\x}_{ii}  \,\x\big) \\[1mm]
\Vdots &        \bm{0} \\
\end{pNiceArray} \,\,\,\,\, \overset{(a)}{=} \hspace{-6mm}\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\
\Vdots &        \bm{0} \\[1mm]
\wMasupp_{i \,\bullet}   &  \bm \Lambda^{\x}_{ii}  \,\, \x \kro \wMbsupp^\top_{\bullet \,i}\\[2mm]
\wMbsupp_{\bullet \,i}    & \wMasupp^\top_{i \,\bullet}  \big(\bm \Lambda^{\x}_{ii}  \,\x\big)  \kro  \Im_{K}\\[1mm]
\Vdots &        \bm{0} \\
\end{pNiceArray}\,,\]

where, in (a) we have used the fact that for vectors $\ab,\,\bb$ we have that $\ab^\top \,\kro\, \bb = \bb \,\kro \,\ab^\top = \bb\, \ab^\top$, as well as the fact that $\wMasupp^\top_{i \,\bullet} \big(\bm \Lambda^{\x}_{ii}  \,\x\big)$ is a scalar which allows us to commute the factors in the corresponding Kronecker product.

Finally, we can express the above Jacobian across all the samples in the data matrix, as stated in the lemma:
code-for-first-row = ,
code-for-first-col = ,

\[\nabla_{\btheta} F_{\theta_i}(\Xm)^\top = \hspace{-1cm}\begin{pNiceArray}[first-row,first-col,nullify-dots]{c}
& \\
\Vdots &        \bm{0} \\[1mm]
\wMasupp_{i \,\bullet}     & \Xm \, \bm \Lambda^i \kro \wMbsupp^\top_{\bullet \,i}  \\[2mm]
\wMbsupp_{\bullet \,i}    &  \wMasupp_{i \,\bullet}^\top  \Xm \, \bm \Lambda^i \kro \Im_{K}  \\[1mm]
\Vdots &        \bm{0} \\
\end{pNiceArray}\,.
\]

Here, we utilized that $\Am \kro \Bm = \left[\Am_{\bullet\, 1} \kro \Bm,\, \cdots,\, \Am_{\bullet\, n} \kro \Bm\right]$ for some arbitrary matrix $\Am$ containing $n$ columns. Besides, we have collected the activation derivatives for the $i$-th neuron, i.e., $\bm \Lambda^\x_{ii}\,=\sigma^\prime( \Vm_{i \,\bullet}^\top \,  \x)\,$, across all samples $\x$, into the diagonal matrix $\bm \Lambda^i \in \mathbb{R}^{N\times N}$.

From the above Lemma, we can also see that the benefit of analyzing via the unit-networks is that we only have to deal with the activation derivatives of a single neuron at a time. Besides, now that we know the structure of the unit-network Jacobian, we will  analyze its rank. But before, let's recall the assumption <ref> from the main text, in our current notation:
For each active hidden neuron $i$, the weighted input covariance has the same rank as the overall input covariance, i.e.,  $\rank(\E[\alpha_{\x}\,\x \x^\top]) = \rank(\covx) = r$, with $\,\alpha_{\x} = {\sigma^\prime(\x^\top \, \wMasupp_{i\,\bullet})}^2$.

This assumption can be translated into finite-sample case as follows. First, note that the (uncentered) input covariance $\covx$ corresponds to $\frac{1}{N} \Xm \Xm^\top$, while the weighted covariance $\E[\alpha_{\x}\,\x \x^\top]$ corresponds to the matrix $\frac{1}{N}\Xm\bm\Lambda^i \bm\Lambda^i \Xm^\top$. This is straightforward to check, and notice $\alpha_\x = (\bm \Lambda^\x_{ii})^2$. Then, the equivalent assumption is to require $\rank\left(\Xm\bm\Lambda^i \bm\Lambda^i \Xm^\top\right) = \rank\left(\Xm\Xm^\top\right) = r$, ignoring the constant $\frac{1}{N}$ which does not affect rank. Further, since for any arbitrary matrix $\Am$, we have that $\rank(\Am\Am^\top) =\rank(\Am)$. Thus, our equivalent assumption can be simplified to as follows:
assump:data(finite-sample equivalent)
For each active hidden neuron $i$, assume that $\rank\left(\Xm\bm \Lambda^i\right)=\rank\left(\Xm\right)=r$ , where $\bm \Lambda^i$, as detailed before, contains the activation derivatives across all samples for this neuron $i$.

Under the same setup as Lemma <ref> and Assumptions <ref>, <ref> (or equivalently <ref>), the rank of the Jacobian matrix, $\nabla_{\btheta} \Fni{i}(\Xm)$, of the $i$-th unit-network is given by:

\[\rank({\nabla_{\btheta} \Fni{i}(\Xm)}) = r + K - 1\,.\]

From Lemma <ref>, the Jacobian matrix is given by (ignoring the zero blocks which do not matter for the analysis of rank),

$$\nabla_{\btheta} \Fni{i}(\Xm)^\top = \begin{pmatrix}
\Xm \, \bm \Lambda^i\kro \wMbsupp^\top_{\bullet \,i}  \\[2mm]
\wMasupp^\top_{i \,\bullet}  \Xm \, \bm \Lambda^i \kro \Im_{K} \\[1mm]
\end{pmatrix} = \underbrace{\begin{pmatrix}
\Im_{r} \kro \wMbsupp^\top_{\bullet \,i}\\[2mm]
\wMasupp^\top_{i \,\bullet} \kro \Im_{K}  \\[1mm]
\end{pmatrix}}_{\Am_i \,\in\, \mathbb{R}^{(r+K)\times Kr}} \underbrace{\left(\Xm \, \bm \Lambda^i\kro \Im_{K} \right)}_{\in \, \mathbb{R}^{Kr\times KN}}

Now this factorization reveals the familiar $\Zm$-like structure, and so the matrix labelled $\Am_i$ in the above factorization has rank equal to $r + K - 1$ by Lemma <ref>. And, $\rank(\Xm \, \bm \Lambda^i\kro \Im_{K}) = K \rank(\Xm \, \bm \Lambda^i)= Kr$, by employing assumption <ref>. Thus, this matrix $\Xm \, \bm \Lambda^i\kro \Im_{K}$ is right invertible. Hence, we have:

\[\rank({\nabla_{\btheta} \Fni{i}(\Xm)}) =\rank({\nabla_{\btheta} \Fni{i}(\Xm)}^\top) = \rank(\Am_i) = r + K-1\,.\]

§.§ Proof of Theorem <ref>

Now, that we are equipped to prove the Theorem, and let us recall its statement from the main text:
Consider a 1-hidden layer network with non-linearity $\sigma$ such that $\sigma(z) = \sigma^{\prime}(z) z$ and let $\widetilde{M}$ be the # of active hidden neurons (i.e., probability of activation $>0$). Then, under assumption <ref> and <ref>,  rank of $\HO$ is given as,