Theoretical Characterisation of the Gauss Newton Conditioning in Neural Networks

Authors: jim  zhao
sidak pal singh
aurelien  lucchi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 23:01:54

Abstract: The Gauss-Newton (GN) matrix plays an important role in machine learning, most evident in its use as a preconditioning matrix for a wide family of popular adaptive methods to speed up optimization. Besides, it can  also provide key insights into the optimization landscape of neural networks. 
In the context of deep neural networks, understanding the GN matrix involves studying the interaction between different weight matrices as well as the dependencies introduced by the data, thus rendering its analysis challenging.
In this work, we take a first step towards theoretically characterizing the conditioning of the GN matrix in neural networks. We establish tight bounds on the condition number of the GN in deep linear networks of arbitrary depth and width, which we also extend to two-layer ReLU networks.
We expand the analysis to further architectural components, such as residual connections and convolutional layers. 
Finally, we empirically validate the bounds and uncover valuable insights into the influence of the analyzed architectural components.

Main paper: 
 introduction: 
 The curvature is a key geometric property of the loss landscape, which is characterized by the Hessian matrix or approximations such as the Gauss-Newton (GN) matrix, and strongly influences the convergence of gradient-based optimization methods. In the realm of deep learning, where models often have millions of parameters, understanding the geometry of the optimization landscape is essential to understanding the effectiveness of training algorithms. The Hessian matrix helps identify the directions in which the loss function changes most rapidly, aiding in the selection of appropriate learning rates and guiding optimization algorithms to navigate the complex, high-dimensional space of parameters. However, in practice, the Hessian is not easily accessible due to the high computational cost and memory requirements. Instead, the GN matrix (or its diagonal form) is commonly employed in adaptive optimization methods such as Adam [Kingma and Ba, 2014] with the goal to improve the conditioning of the landscape. Although the Gauss-Newton matrix G O is only an approximation to the full Hessian matrix, it does seem to capture the curvature of the loss very well given the success of many second-order optimization methods based on approximations of the Gauss-Newton matrix, such as K-FAC [Martens and Grosse, 2020] , Shampoo [Gupta et al., 2018] or Sophia [Liu et al., 2023] . Particularly interesting is the last method, in which the authors observe that their optimizer based on the Gauss-Newton matrix performs even better than their optimizer based on the full Hessian matrix, implying that the Gauss-Newton matrix is a good preconditioner and captures the curvature of the loss landscape well. The prevalence of adaptive optimizers and their success in training neural networks illustrates the ability of the GN matrix to capture curvature information similar to that of the Hessian matrix.Although the non-convexity of neural loss landscapes can be considered as a given, some landscapes can be tougher to traverse than others. This can, for instance, arise due to, or be amplified by, aspects such as: (i) disparate scales of the input data and intermediate features, (ii) initializing the optimization process at a degenerate location in the landscape (imagine an entire set of neurons being dead), (iii) architectural aspects such as width, depth, normalization layers, etc.Given the close relation of the GN matrix to the network Jacobian (with respect to the parameters) and hence signal propagation [Lou et al., 2022] , its significance also reaches more broadly to problems such as model pruning and compression. For instance, it has been extensively observed empirically that training extremely sparse neural networks from scratch poses a significantly greater challenge compared to pruning neural networks post-training [Evci et al., 2019 ]. In an illustrative example shown in Fig. 1 , we observe that training sparse networks from scratch using stochastic gradient descent results in a slowdown of training, which is also reflected in an increase in the condition number of the GN matrix (similar experiments on Vision Transformers and other architectures can be found in Appendix I.1). This underscores the relevance of the conditioning of the GN matrix for understanding the behaviour of gradient-based optimizers and emphasizes the significance of maintaining a well-behaved loss landscape of the neural network throughout training.In fact, many fundamental components of deep learning frameworks, such as skip-connections as well as various normalization techniques [Ioffe and Szegedy, 2015 , Salimans and Kingma, 2016 , Ba et al., 2016 , Wu and He, 2018 , Miyato et al., 2018] , have been to some degree designed to mitigate the challenges posed by ill-conditioning in neural network optimization. This line of work continues to expand to this date, seeking out better normalization schemes, optimization maneuvers, and regularization strategies [Kang et al., 2016 , Wan et al., 2013] that allow for easier and faster optimization while avoiding directions of pathological curvature. But despite the numerous approaches to redress the poor conditioning of the landscape, understanding the precise factors of ill-conditioning within the network structure, and their relative contributions, has remained largely underdeveloped.Hence, our aim in this work is to carry out a detailed theoretical analysis of how the conditioning of the Gauss-Newton matrix is shaped by constituting structural elements of the network -i.e., the hidden-layer width, the depth, and the presence of skip connections, to name a few. We will shortly formally introduce the Gauss-Newton matrix and how it is connected to the Hessian of the loss function. Concretely, we would like to provide tight bounds for the condition number of these two terms as a function of the spectra of the various weight matrices and reveal the interplay of the underlying architectural parameters on conditioning. Furthermore, we aim to investigate the impact of both the dataset's structure and the initialization procedure on the conditioning of the loss landscape.Contributions. Taking inspiration from prior theoretical analyses of deep linear networks, we make a first foray into this problem by rigorously investigating and characterizing the condition number of the GN matrix for linear neural networks (the extension to the second term of the Gauss-Newton decomposition is touched upon in the Appendix). Our analysis holds for arbitrary-sized networks and unveils the intriguing interaction of the GN matrix with the conditioning of various sets of individual layers and the data-covariance matrix. We also complement our analysis in the linear case by studying the effect of non-linearities, via the Leaky-ReLU activation, albeit for two-layer networks. Importantly, as a consequence of our analysis, we show the precise manner in which residual networks with their skip connections or batch normalization can help enable better conditioning of the loss landscape. While our work builds on Singh et al. [2021] , our main contribution is the introduction of tight upper bounds for the condition number of the Gauss-Newton (GN) matrix for linear and residual networks of arbitrary depth and width. To the best of our knowledge, this has not been addressed in the literature before. Lastly, given that our bounds are agnostic to the specific values of the parameters in the landscape, we show the phenomenology of conditioning during the training procedure and the corresponding validity of our bounds.2 Setup and backgroundSetting. Suppose we are given an i.i.d. dataset S = {(x 1 , y 1 ), . . . ,(x n , y n )}, of size |S| = n, drawn from an unknown distribution p X,Y , consisting of inputs x ∈ X ⊆ R d and targets y ∈ Y ⊆ R k .Based on this dataset S, consider we use a neural network to learn the mapping from the inputs to the targets, F θ θ θ : X → Y, parameterized by θ θ θ ∈ Θ Θ Θ ⊆ R p . To this end, we follow the framework of Empirical Risk Minimization [Vapnik, 1991] , and optimize a suitable loss function L : Θ Θ Θ → R. In other words, we solve the following optimization problem,θ θ θ ⋆ = argmin θ θ θ∈Θ Θ Θ L(θ θ θ) = 1 n n i=1 ℓ (θ θ θ; (x i , y i )) ,say with a first-order method such as (stochastic) gradient descent and the choices for ℓ could be the mean-squared error (MSE), cross-entropy (CE), etc. For simplicity, we will stick to the MSE loss.Gauss-Newton Matrix. We analyze the properties of the outer gradient product of the loss function, G O , which we call the Gauss-Newton matrix, defined asEQUATIONwhere, ∇ θ θ θ F θ θ θ (x i ) ∈ R p×k is the Jacobian of the function with respect to the parameters θ θ θ, p is the number of parameters, k is the number of outputs or targets. This outer product of the gradient is closely related to the Hessian of the loss function via the Gauss-Newton decomposition [Schraudolph, 2002 , Sagun et al., 2017 , Martens, 2020 , Botev, 2020] , hence the chosen name, which decomposes the Hessian via the chain rule as a sum of the following two matrices:H L = H O + H F = 1 n n i=1 ∇ θ θ θ F θ θ θ (x i ) ∇ 2 F θ θ θ ℓ i ∇ θ θ θ F θ θ θ (x i ) ⊤ + 1 n n i=1 K c=1 [∇ F θ θ θ ℓ i ] c ∇ 2 θ θ θ F c θ θ θ (x i ),where ∇ 2 F θ θ θ ℓ i ∈ R k×k is the Hessian of the loss with respect to the network function, at the i-th sample. Note that if ℓ i is the MSE loss, thenH O = G O .Remark R1 (Difference between H L and G O ). When considering MSE loss, the difference between the Gauss Newton matrix G O and the Hessian of the loss function H L depends on both the residual and the curvature of the network F θ θ θ (x). Thus, close to convergence when the residual becomes small, the contribution of H F will also be negligible and G O is essentially equal to H L . Furthermore, Lee et al. [2019] show that sufficiently wide neural networks of arbitrary depth behave like linear models during training with gradient descent. This implies that the Gauss-Newton matrix is a close approximation of the full Hessian in this regime throughout training.Condition number and its role in classical optimization. Consider we are given a quadratic problem, argmin w∈R p 1 2 w ⊤ Aw, where A ≻ 0 is a symmetric and positive definite matrix. The optimal solution occurs for w * = 0. When running gradient descent with constant step size η > 0, the obtained iterates would be w k = (I -ηA)w k-1 . This yields a convergence rateof ∥w k -w * ∥ ∥w k-1 -w * ∥ ≤ max |1 -ηλ max (A)|, |1 -ηλ min (A)| . The best convergence is obtained for η = 2 (λ max (A) + λ min (A)) -1 , resulting in ∥w k ∥ ≤ κ-1 κ+1 ∥w k-1 ∥, where κ(A) = λmax(A)λmin(A) is called the condition number of the matrix. This ratio which, intuitively, measures how disparate are the largest and smallest curvature directions, is an indicator of the speed with which gradient descent would converge to the solution. When κ → ∞, the progress can be painfully slow, and κ → 1 indicates all the curvature directions are balanced, and thus the error along some direction does not trail behind the others, hence ensuring fast progress.Effect of condition number at initialization on the convergence rate As the condition number is a very local property, it is in general hard to connect the conditioning at network initialization to a global convergence rate. However, we would like to argue below that an ill-conditioned network initialization will still affect the rate of convergence for gradient descent (GD) in the initial phase of training. Let us denote the Lipschitz constant by L and the smoothness constant by µ. Furthermore, let the step size be such that η k ≤ 1 L . We present a modified analysis of GD for strongly convex functions, where we use local constants µ(k) and L(k) instead of the global smoothness and Lipschitz constant, respectively. Then by the definition of a single step of gradient descent and using the strong convexity and smoothness assumption2 we have:||θ k+1 -θ * || 2 ≤ (1 -η k µ)||θ k -θ * || 2 (2)So by recursively applying (2) and replacing µ by the local smoothness constants µ(k):||θ k -θ * || 2 ≤ k-1 i=0 (1 -η i µ(i))||θ 0 -θ * || 2 (3)One can clearly see the effect of µ(0) in the bound, which is even more dominant when µ(k) changes slowly. Of course, the effect of µ(0) attenuates over time, and that's why we are talking about a local effect. However, one should keep in mind that overparametrization leads the parameter to stay closer to initialization (at least in the NTK regime [Lee et al., 2019] ). 
related work: 
 Since the Gauss-Newton matrix is intimately related to the Hessian matrix, and the fact that towards the end of training, the Hessian approaches the Gauss-Newton matrix [Singh et al., 2021] , we carry out a broader discussion of the related work, by including the significance of the Hessian at large.There is a rich and growing body of work that points towards the significance of various Hessianbased measures in governing different aspects of optimization and generalization. One popular hypothesis [Hochreiter and Schmidhuber, 1997 , Keskar et al., 2016 , Dziugaite and Roy, 2017 , Chaudhari et al., 2019] is that flatter minima generalize better, where the Hessian trace or the spectral norm is used to measure flatness. This hypothesis is not undisputed [Dinh et al., 2017] , and the extent to which this is explanatory of generalization has been put to question recently [Granziol, 2021 , Andriushchenko et al., 2023] . Nevertheless, yet another line of work has tried to develop regularization techniques that further encourage reaching a flatter minimum, as shown most prominently in proposing sharpness-aware minimization [Foret et al., 2021] .(ii) Understanding Architectural and Training aspects of Neural Networks: Some other work has studied the challenges in large-batch training via the Hessian spectrum in Yao et al. [2018] . Also, the use of large learning rates has been suggested to result in flatter minima via the initial catapult phase [Lewkowycz et al., 2020] . The effect of residual connections and Batch normalization on Hessian eigenspectrum were empirically studied in Ghorbani et al. [2019] , Yao et al. [2020] , which introduced PyHessian, a framework to compute Hessian information in a scalable manner. More recently, the so-called edge of stability [Cohen et al., 2021] phenomenon connects the optimization dynamics of gradient descent with the maximum eigenvalue of the loss Hessian. Very recently, the phenomenon of Deep neural collapse was studied in Beaglehole et al. [2024] via the average gradient outer-product.(iii) Applications: There has also been a dominant line of work utilizing the Hessian for second-order optimization, albeit via varying efficient approximations, most notably via K-FAC [Martens and Grosse, 2020] but also others such as Yao et al. [2021] , Liu et al. [2023 ], Lin et al. [2023] . Given its versatile nature, the Hessian has also been used for model compression through pruning [LeCun et al., 1989 , Hassibi et al., 1993 , Singh and Alistarh, 2020] as well as quantization [Dong et al., 2019 , Frantar et al., 2023] , but also in understanding the sensitivity of predictions and the function learned by neural networks via influence functions [Koh and Liang, 2020, Grosse et al., 2023] , and countless more.Theoretical and empirical studies of the Hessian. There have been prior theoretical studies that aim to deliver an understanding of the Hessian spectrum in the asymptotic setting [Pennington and Bahri, 2017, Jacot et al., 2019] , but it remains unclear how to extract results for finite-width networks as used in practice. Besides, past work has analyzed the rank, empirically [Sagun et al., 2016 [Sagun et al., , 2017] ] as well as theoretically [Singh et al., 2021 [Singh et al., , 2023]] . Further, the layer-wise Hessian of a network can be roughly approximated by the Kronecker product of two smaller matrices, whose top eigenspace has shown to contain certain similarities with that of the Hessian [Wu et al., 2020] . In a different line of work by Liao and Mahoney [2021] , the limiting eigenvalue distribution of the Hessian of generalized generalized linear models (G-GLMs) and the behaviour of potentially isolated eigenvalue-eigenvector pairs is analyzed.Hessian and landscape conditioning. The stock of empirical repertoire in deep learning has been enriched by successful adaptive optimization methods plus their variants [Kingma and Ba, 2014] as well as various tricks of the trade, such as Batch Normalization [Ioffe and Szegedy, 2015] , Layer Normalization [Ba et al., 2016] , orthogonal initialization [Saxe et al., 2013 , Hu et al., 2020] , and the kind, all of which can arguably be said to aid the otherwise ill-conditioning of the landscape. There have also been theoretical works establishing a link between the conditioning of the Hessian, at the optimum and the double-descent like generalization behavior of deep networks [Belkin et al., 2019 , Singh et al., 2022] .Gauss-Newton matrix and NTK. In the context of over-parametrized networks, G O is for instance connected to the (empirical) Neural Tangent Kernel, which has been the focus of a major line of research in the past few years [Jacot et al., 2018 , Wang et al., 2022 , Yang, 2020] as the NTK presents an interesting limit of infinite-width networks. As a result the asymptotic spectrum and the minimum eigenvalue of the NTK has been studied in [Nguyen et al., 2022, Liu and Hui, 2023] , but the implications for finite-width networks remain unclear.Despite this and the other extensive work discussed above, a detailed theoretical study on the Gauss Newton conditioning of neural networks has been absent. In particular, there has been little work trying to understand the precise sources of ill-conditioning present within deep networks. Therefore, here we try to dissect the nature of conditioning itself, via a first principle approach. We hope this will spark further work that aims to precisely get to the source of ill-conditioning in neural networks and, in a longer horizon, helps towards designing theoretically-guided initialization strategies or normalization techniques that seek to also ensure better conditioning of the GN matrix. 
theoretical characterisation: 
 The main part of this paper will focus on analyzing the conditioning of G O in Eq. (1) as prior work [Ren and Goldfarb, 2019, Schraudolph, 2002] has demonstrated its heightened significance in influencing the optimization process. We will further discuss an extension to H F in Appendix F. Tying both bounds together yields a bound on the condition number of the overall loss Hessian in some simple setups.Pseudo-condition number. Since the GN matrix of deep networks is not necessarily full rank [Sagun et al., 2017 , Singh et al., 2021] , we will analyze the pseudo-condition number defined as the ratio of the maximum eigenvalue over the minimum non-zero eigenvalue. This choice is rationalized by the fact that gradient-based methods will effectively not steer into the GN null space, and we are interested in the conditioning of the space in which optimization actually proceeds. For brevity, we will skip making this distinction between condition number and pseudo-condition number hereafter.We will start with the case of linear activations. In this case a network with L layers can be expressed byF θ θ θ (x) = W L W L-1 • • • W 1 x, with W ℓ ∈ R a ℓ ×a ℓ-1 for ℓ = 1, . . . , L and a L = k, a 0 = d.To facilitate the presentation of the empirical work, we will assume that the widths of all hidden layers are the same, i.e. α ℓ = m for all ℓ = 1, . . . , L -1. Also, let us denote by Σ = 1 n n i=1 x i x T i ∈ R d×d the empirical input covariance matrix. Furthermore, we introduce the shorthand-notation Then, Singh et al. [2021] show that the GN matrix can be decomposed as G O = U(I k ⊗ Σ)U ⊤ , where I k is the identity matrix of dimension k and U ∈ R p×kd is given by:W k:ℓ = W k • • • W ℓ for k > ℓ and k < ℓ, W k:ℓ = W k ⊤ • • • W ℓ ⊤ .U = W 2:L ⊗ I d . . . W ℓ+1:L ⊗ W ℓ-1:1 . . . I k ⊗ W L-1:1 ⊤ . By rewriting U(I K ⊗ Σ)U ⊤ = U(I K ⊗ Σ 1/2 )(I K ⊗ Σ 1/2 )U ⊤ ,where Σ 1/2 is the unique positive semi-definite square root of Σ and noting that AB and BA have the same non-zero eigenvalues, we have that the non-zero eigenvalues of G O are the same as those ofG O = (I K ⊗ Σ 1/2 )U ⊤ U(I K ⊗ Σ 1/2 ) , (4) where U ⊤ U ∈ R Kd×Kd is equal to U ⊤ U = L ℓ=1 W L:ℓ+1 W ℓ+1:L ⊗ W 1:ℓ-1 W ℓ-1:1 .(5)Warm-up: the one-hidden layer case. In the case of one-hidden layer network, F θ θ θ (x) = WVx, we have that G O = WW ⊤ ⊗ Σ + I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 . To(W) + σ 2 min (V)).Then the condition number of GN for the one-hidden layer network with linear activations with m > max{d, k} is upper bounded byEQUATION)Proof sketch. Choosing A = WW ⊤ ⊗ Σ and B = I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 and i = j = 1 for the Weyl's inequality, i = j = n for the dual Weyl's inequality and using the fact that WW ⊤ ⊗ Σ andI k ⊗ Σ 1/2 V ⊤ VΣ 1/2 are positive semidefinite yield the result.It is important to point out that the convex combination in Eq. ( 6) is crucial and a more naive bound where we take the maximum of both terms is instead too loose, see details in Appendix E. Later on, we will observe that the general case with L layers also exhibits a comparable structure, wherein the significance of the convex combination becomes even more pronounced in deriving practical bounds.Remark R2 (Role of the data covariance). Besides the above dependence in terms of the convex combination of the bounds, we also see how the conditioning of the input data affects the conditioning of the GN spectra. This is observed in Figure 19 , where the condition number of the GN matrix is calculated on whitened and not whitened data. This observation might also shed light on why data normalization remains a standard choice in deep learning, often complemented by the normalization of intermediate layer activations.Following our examination of the single-hidden layer case, we now broaden our analysis to include Llayer linear networks. As before, we will first derive an expression of the GN matrix and subsequently bound the largest and smallest eigenvalue separately to derive a bound on the condition number.Obtaining the GN matrix involves combining (4) and ( 5), which yieldsG O = (I K ⊗ Σ 1/2 )U ⊤ U(I K ⊗ Σ 1/2 ) = L l=1 (W L:ℓ+1 W ℓ+1:L ) ⊗ (Σ 1/2 W 1:ℓ-1 W ℓ-1:1 Σ 1/2 ).By repeatedly applying Weyl's inequalities, we obtain an upper bound on the condition number.Lemma 2. Assume that m > max{d, k} and thatα ℓ := σ 2 min (W L:ℓ+1 ) • σ 2 min (W 1:ℓ-1 ) > 0 ∀ℓ = 1, . . . , L. Let γ ℓ := α ℓ L i=1 αi .Then the condition number of the GN matrix of a L-layer linear network can be upper-bounded in the following way:κ( G O ) ≤ κ(Σ) L ℓ=1 γ ℓ κ(W L:ℓ+1 ) 2 κ(W 1:ℓ-1 ) 2 ≤ κ(Σ) max 1≤ℓ≤L κ(W L:ℓ+1 ) 2 κ(W 1:ℓ-1 ) 2 .As mentioned earlier, we observe the same convex combination structure which, as we will soon see experimentally, is crucial to obtain a bound that works in practice.Empirical validation. The empirical results in Figure 2a show that the derived bound seems to be tight and predictive of the trend of the condition number of GN at initialization. If the width of the hidden layer is held constant, the condition number grows with a quadratic trend. However, the condition number can be controlled if the width is scaled proportionally with the depth. This gives another explanation of why in practice the width of the network layers is scaled proportionally with the depth to enable faster network training. Furthermore, Figure 3 empirically demonstrates the importance of bounding the condition number as a convex combination of the condition number of the weight matrices, as simply taking the maximum makes the bound vacuous. This is due to a 'self-balancing' behavior of the individual terms in Lemma 2, which is further elaborated in Appendix E. We show that the difference in magnitude between each condition number is largely driven by their smallest eigenvalue. Thus at the same time, they also have a smaller weight in the convex combination, making the overall condition number grow slower than the maximum bound would predict.Following the examination of linear networks, our focus will now shift to the analysis of linear networks where we introduce a skip connection with weight β between each layer, as illustrated in Figure 4 . This results in a residual network of the formF θ (x) = (W (L) + βI) • • • (W (1) + βI)x.Given the architecture of the residual network, a substantial portion of the analysis conducted in Section 4.2 for the L-layer linear neural network can be repurposed.The key insight lies in the realization that when the skip connection precisely bypasses one layer, it implies a modification of Eq. ( 5) toU ⊤ U = L ℓ=1 W L:ℓ+1 β W ℓ+1:L β ⊗ W 1:ℓ-1 β W ℓ-1:1 β , where we define W k:l β := (W k + βI) • • • (W l + βI) if k > l and W k:l β := (W k + βI) ⊤ • • • (W l + βI) ⊤ if k < l.I denotes the rectangular identity matrix of appropriate dimensions. Note, that we can apply the upper bounds derived in Lemma 2 analogously to arrive atEQUATIONwhereγ β ℓ := α β ℓ L i=1 α β i and analogously to the L-layer case α β ℓ := σ 2 min (W L:ℓ+1 β ) • σ 2 min (W 1:ℓ-1 β ) for ℓ = 1, . . . , L.Let us now analyze how skip connections affect the condition number. Denoting the SVD ofW ℓ = U ℓ S ℓ V ℓ ⊤ note that W ℓ +βI = U ℓ (S ℓ +βI)V ℓ T, and therefore σ i (W ℓ +βI) = σ i (W ℓ )+β for all singular values of W ℓ . Incorporating skip connections results in a spectral shift of each weight matrix to the positive direction by a magnitude of β. Furthermore, we will assume that the left and right singular vectors of the layers in-between coincide, that is V ℓ = U ℓ-1 . This assumption is fulfilled if the initialization values are sufficiently close to zero [Saxe et al., 2013] . Then for each ℓ we haveEQUATION) Since σmax(W i )+β σmin(W i )+β < σmax(W i ) σmin(W i )for all β > 0, we can conclude that the second upper bound in Eq. ( 7) for residual networks is smaller than the second upper bound for linear networks in Lemma 2. We verify that this observation also holds for the actual condition number and the tight upper bound in Eq. ( 7) as can be seen in Figure 5 . The latter shows that adding skip connections improves the conditioning of the network. Furthermore, Figure 20 in the appendix shows that the condition number is improved more for larger β because this pushes the condition number towards one as can be seen from Eq. ( 8). To conclude, the analysis of the condition number provides a rigorous explanation for the popularity of residual networks compared to fully connected networks. This is in line with existing empirical results in the literature [He et al., 2016 , Li et al., 2018 , Liu et al., 2019] .Remark: extension to convolutional layers. Note that the analysis on fully connected layers can also be extended to convolutional layers in a straightforward way by making use of the fact that the convolution operation can be reformulated in the form of a matrix-vector product using Toeplitz matrices, as discussed for instance in Singh et al. [2023] . In this case, we can apply the same analysis as for fully connected networks. See Appendix B for more details. 
extension to non-linear activations: 
 After the analysis of networks with a linear activation in the previous section, we will extend our results to non-linear activation functions σ that satisfy σ(z) = σ ′ (z)z for one-hidden layer networks F θ (x) = Wσ(Vx). Specifically, this extension includes all piece-wise linear activation functions such as ReLU or Leaky ReLU. For this purpose, we first rewrite the network output as a super-position of unit-networks F θi (x) over the number of hidden neurons m asEQUATIONwhere W •,i denotes the i-th column of W and V i,• denotes the i-th row of V. Using the following lemma, we can then derive an expression for the Gauss-Newton matrix.Lemma 3. (Lemma 25 in Singh et al. [2021]) Let F θi (x) = W •,i σ(V i,•x) be a unit-network corresponding to i-th neuron, with the non-linearity σ such that σ(z) = σ ′ (z)z. Let X ∈ R d×n denote the data matrix. Further, let Λ i ∈ R n×n be defined as (Λ i ) jj = σ ′ (V i,• x) j , for j = 1, . . . , n and zero elsewhere. Then the Jacobian matrix ∇ θ F θi (X) is given (in transposed form) by:∇ θ F θi (X) = 0 XΛ i ⊗ W ⊤ • i V i • XΛ i ⊗ I k 0 ⊤ .The Jacobian matrix of the full network is simply the sum over all units∇ θ F θ (X) = m i=1 ∇ θ F θi (X). Now that we have an expression for the Jacobian matrix, we can compute the GN matrix as G O = ∇F θ (X)∇F θ (X) ⊤ . Note that the non-zero eigenvalue of G O and G O := ∇F θ (X) ⊤ ∇F θ (X) ∈ R kn×kn are the same, for which we haveG O = m i=1 ∇F θi (X) ⊤ ∇F θi (X) = m i=1 G i O with G i O := Λ i X ⊤ XΛ i ⊗ W •i W ⊤ •i + Λ i X ⊤ V ⊤ i• V i• XΛ i ⊗ I k ∈ R kn×kn, since the mixed terms ∇F θi (X) ⊤ ∇F θj (X) with i ̸ = j are zero due to the block structure of the Jacobian of the unit-networks. Using the mixed product property of the Kronecker product we obtain 7), respectively.EQUATIONTo improve readability, we will writeΓ := m i=1 Λ i X ⊤ V ⊤ i• V i• XΛ i hereafter.Leaky ReLU activation. Let's now consider the case where the non-linear activation σ is LeakyReLU(x) = max{αx, x} for some constant α ∈ [0, 1). Typically α is chosen to be close to zero, e.g. α = 0.01. Again, the condition number κ( G O ) can be upper bounded by bounding the extreme eigenvalues independently. By observing that all terms in Eq. ( 10) are symmetric, we can again apply the Weyl's and dual Weyl's inequality. However to get a lower bound larger than zero forλ min ( G O ) a different approach is needed as W •,i W ⊤•,i is rank-deficient and the same steps to bound the smallest eigenvalue would lead to a vacuous value of zero. Instead, we will use the following observation that we can bound the extreme eigenvalues of the sum of a Kronecker product of PSD matrices (11) . Note that the y-axis is log scaled.A i , B i for i = 1, . . . , m by λ min ( m i=1 (A i ⊗ B i )) ≥ min 1≤j≤m λ min (A j ) • λ min ( m i=1 B i ).By first using Weyl's inequalities and subsequently applying the above bound with A i = Λ i X ⊤ XΛ i andB i = W •,i W ⊤•,i we achieve the following bound for the condition number for κ( G O ). Lemma 4. Consider an one-hidden layer network with a Leaky-ReLU activation with negative slope α ∈ [0, 1]. Then the condition number of the GN matrix defined in Equation (10) can be bounded as:EQUATIONThe proof can be found in Appendix H. We further empirically validate the tightness of the upper bounds on a subset of the MNIST dataset (n = 500), where we chose α = 0.01, which is the default value in Pytorch [Paszke et al., 2019] . As can be seen in Figure 6 , contrary to the linear setting the condition number increases with width. The upper bound also becomes tighter and more predictive with increasing width.Comparison to linear network. By running the same setup for varying values of α ∈ [0, 1] we can interpolate between ReLU (α = 0) and linear activation (α = 1) to see how the introduction of the non-linearity affects the conditioning. In Figure 23 , which had to be deferred to the appendix due to space constraints, we can see how reducing the value of α seems to consistently improve the conditioning of the GN matrix for the one-hidden layer case. Based on the observation that conditioning of the input data affects the conditioning of the GN spectrum, we also tested whether Batch normalization (BN) [Ioffe and Szegedy, 2015] , which is a very commonly used normalization scheme, has a similar effect on the condition number. For this, the condition number of the GN matrix was calculated for a one-hidden layer linear network with and without a Batch normalization layer. The experiment was run on a downsampled and subsampled version of Cifar-10, which was converted to grayscale with d = 64and n = 1000. The data was not whitened to see the effect of BN more strongly. As can be seen in Figure 7 , we indeed observe a clear improvement of the condition number when Batch normalization layers are added. Also, the trend of improved conditioning with increasing width remains after adding BN.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: Code to run the experiments is provided in the supplementary material together with instructions how to reproduce the results.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: For each experimental setting, the architectural setup, the initialization scheme and the optimizer together with its hyperparameters were specified. Since our work only evaluated the training setting, a data split was not necessary.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. Guidelines: 
discussion and conclusion: 
 Summary. In this work, we derived new analytical bounds for the condition number of neural networks, showcasing the following key findings: a) The conditioning of the input data has a nearly proportional impact on the conditioning of the GN spectra, underscoring the significance of data normalization. Empirical evidence further demonstrates that Batch Normalization similarly enhances the condition number. b) For linear networks, we showed that the condition number grows quadratically with depth for fixed hidden width. Also, widening hidden layers improves conditioning, and scaling the hidden width proportionally with depth can compensate for the growth. c) We showed how adding residual connections improves the condition number, which also explains how they enable the training of very deep networks. d) Preliminary experiments suggest that the ReLU activation seems to improve the conditioning compared to linear networks in the one-hidden layer case.Interesting use cases of our results. Through our analysis, we highlighted that the condition number as a tool from classical optimization is also an attractive option to better understand challenges in neural network training with gradient-based methods. Especially, knowing how different architectural choices will affect the optimization landscape provides a more principled way to design the network architecture for practitioners, for instance how to scale the width in relation to the depth of the network. The paper also gives a justification of why pruned networks are more difficult to train as they have worse conditioning. Although this is not our focus, it is possible that our analysis could inspire better techniques for pruning neural networks.Limitations and future work. We made a first step toward understanding the impact of different architectural choices on the conditioning of the optimization landscape. However, there are still many design choices, that have not been covered yet, such as an extension to non-linear networks for arbitrary layers, other architectures, such as transformers, and analytic bounds for normalization schemes, such as batch or layer normalization, which we will leave to future work. Another limitation of our current work is that the derived upper bounds are agnostic to the training dynamics. Therefore, they cannot distinguish the difference between random initializations and solutions of deep learning models after convergence. Incorporating training dynamics into the upper bounds would allow us to characterize solutions found by different architectures, which is left for future work. Another future direction is to extend the analysis to the Generalized Gauss-Newton matrix, which is particularly relevant for training with cross-entropy loss.• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: This paper is a theoretic contribution and as such we do not see any direct societal impact of the work performed.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This paper poses no such risks for misuse or dual-use as it is primarily theoretical in nature.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: The creators and original owners of the datasets used (MNIST and Cifar-10) the creators and original owners of models which were modified and used (ViT, ResNet20, ResNet32) were properly credited. In the case of the ResNet20, and ResNet32, the original owner Jiang Wenjie provided explicit permission to use the code.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [NA] Justification: This paper does not release new assets. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: A.1 Weyl's inequality and dual Weyl's inequalityIn this section, we first briefly review the Weyl's inequality and its dual form. We then demonstrate how these inequalities can be used to bound the extreme eigenvalues and thus the condition number of the sum of two Hermitian matrices A, B by the extreme eigenvalues of A and B. Lemma 5. Weyl [1912] Let A, B ∈ R n×n be two Hermitian matrices, and let λ 1 ≥ . . . ≥ λ n denote the ordered eigenvalues, then the Weyl's inequality and dual Weyl's inequality state thatλ i+j-1 (A + B) ≤ λ i (A) + λ j (B), i, j ≥ 1, i + j -1 ≤ n (Weyl) λ i+j-n (A + B) ≥ λ i (A) + λ j (B), i, j ≥ 1, i + j -n ≤ n.(Dual Weyl)In the following, we derive an upper and lower bound of the condition number of a symmetric matrix A, for which we haveκ(A) = |λ max (A)| |λ min (A)| .Recall that given a Hermitian n × n matrix A, we can diagonalize it by the spectral theorem. We obtain the following sequence of real eigenvaluesλ 1 (A) ≥ λ 2 (A) ≥ . . . ≥ λ n (A).We will make use of the Weyl's inequality and the dual Weyl's inequality to bound the eigenvalues of the sum of two Hermitian matrices A + B.From Lemma 5, we obtain the following inequalities: From Weyl's inequality we haveλ 1 (A + B) ≤ λ 1 (A) + λ 1 (B) λ n (A + B) ≤ λ n-j (A) + λ 1+j (B), for 0 ≤ j ≤ n -1.From the dual Weyl's inequality we haveλ 1 (A + B) ≥ λ n-j (A) + λ 1+j (B), for 0 ≤ j ≤ n -1 λ n (A + B) ≥ λ n (A) + λ n (B).Hence, we can find a lower and upper bound for the condition number of A + BEQUATIONA.2 Further bounds for extreme eigenvalues for positive semidefinite matricesFor square matrices A ∈ R p×p , B ∈ R q×q , denote by λ 1 , . . . , λ p the eigenvalues of A and by µ 1 , . . . , µ q the eigenvalues of B listed by multiplicity. Then the eigenspectrum of A ⊗ B ∈ R pq×pq consists of the eigenvalues λ i µ j i = 1, . . . , p, j = 1, . . . , q.Note that if A, B are additionally positive semidefinite (PSD), we have the following equality for the largest and smallest eigenvalues ofEQUATIONFurthermore, it follows from the sub-multiplicativity of the matrix norm that for square, PSD matrices A, B of same dimensions it holds thatλ max (AB) ≤ λ max (A) • λ max (B) λ min (AB) ≥ λ min (A) • λ min (B).The analysis on fully connected layers can be extended to convolutional layers in a straight forward way, by using the fact that the convolution operation can be reformulated in the form of a matrixvector product using Toeplitz matrices, as discussed for instance in Singh et al. [2023] . This would mean that we can apply the same analysis as for fully connected networks. For completeness we will discuss how below how the convolution operation can be rewritten as a matrix-vector product. For further details, the readers is referred to Singh et al. [2023] .Consider a convolution operation W * x of an input x ∈ R d with m filters of size k ≤ d, which are organized in the matrix W ∈ R m×k . For simplicity a stride of 1 and zero padding is assumed. Let z j:j+k-1 ∈ R k denote the vector formed by considering the indices j to j + k -1 of the original vector z ∈ R d . In this case, the convolution operation can be written asW * x =    ⟨W 1• , x 1:k ⟩ . . . ⟨W 1• , x d-k+1:d ⟩ . . . . . . ⟨W m• , x 1:k ⟩ . . . ⟨W m• , x d-k+1:d ⟩    ∈ R m×(d-k+1) ,where W i• denotes the i-th row of W.Let us further introduce Toeplitz matrices, T Wi• m i=1 , for each filter withT Wi• := toep(W i• , d) ∈ R (d-k+1)×d such that, T Wi• =      Wi1 . . . W ik 0 . . . 0 0 W i1 . . . W ik 0 . . . . . . 0 . . . . . . . . . 0 0 . . . 0 W i1 . . . W ik      .Finally, by stacking the Toeplitz matrices in a row-wise fashion, that isT W :=   T W1• . . . T Wm•   ∈ R m(d-k+1)×d ,we see that the matrix multiplication of T W with an input x gives the same result as vectorizing the convolution operation row-wise, i.e. vec r (W * x) = T W x.For a linear network with L hidden layers, each of which is a convolutional kernel, the network function can be formally represented asEQUATIONwhere the parameters of each hidden layer ℓ are denoted by W (ℓ) ∈ R m ℓ ×m ℓ-1 ×k l . m ℓ is the number of output channels, m ℓ-1 the number of input channels, and k ℓ the kernel size. As we assumed x to be one-dimensional, we have m 0 = 1.Similar to a single convolutional operation, the output of a single convolutional layer can be expressed via Toeplitz matrices. The Toeplitz matrix associated with the ℓ-th convolutional layer W (ℓ) can be expressed byT (ℓ) :=     T W (ℓ) (1,1)• . . . T W (ℓ) (1,m ℓ-1 )• . . . . . . T W (ℓ) (m ℓ ,1)• . . . T W (ℓ) (m ℓ ,m ℓ-1 )•     ∈ R m ℓ d ℓ ×m ℓ-1 d ℓ-1whereEQUATION(i,j)• , d ℓ-1 ) ∈ R d ℓ ×d ℓ-1 to itsassociated Toeplitz matrix. Now, the network function can be equivalently written asEQUATION)Figure 8 : Experiments on the condition number of the Gauss Newton matrix at initialization for a linear two-layer CNN at initialization on a random subsample of n=1000 of MNIST, whitened, with varying kernel size and number of filters. We can see a trend where the number of filters increases the condition number (in analogy to depth in MLPs) and the kernel size improves conditioning (in analogy to width in MLPs).We further provide experiments on linear CNNs in Figure 8 at initialization showcasing that we can observe similar trends where the number of filters increases the condition number of the Gauss Newton matrix (in analogy to depth in MLPs) and the kernel size improves conditioning (in analogy to width in MLPs).This work considers the pseudo condition number of some matrix A ∈ R n×n , which is the ratio of the largest eigenvalue over the smallest non-zero eigenvalue λmax(A) λnz,min(A) . If A is of full rank, the pseudo condition number is equivalent to the condition number. However, if A is rank deficient, it has to be typically estimated numerically if the rank is not known analytically. Based on the numerical estimation of the matrix rank, the smallest eigenvalue is chosen accordingly and the pseudo condition number is calculated.In this section we evaluate the sensitivity of choice of the smallest eigenvalue on the resulting pseudo condition number. This is illustrated in the example of a 1-hidden layer network with ReLU activation function and a hidden width of m = 20 at initialization and after 100 epochs of training. The details of this experimental setup are further specified in Appendix I.2. In Figure 9 we can see that depending on the eigenvalue distribution, the condition number can be quite robust (Epoch=0), but also very sensitive (Epoch=100) to the numerical matrix rank estimation. We conducted an experiment to track the condition number and evaluate how tight the upper bound is throughout training. For this we trained a 3-layer linear network with a hidden width of m = 500 with three different seeds for 5000 epochs with SGD with a mini-batch size of 256 and a constant learning rate of 0.2 on a subset of Cifar-10 (n = 1000) [Krizhevsky et al., 2009] , which has been downsampled to 3 × 8 × 8 images and whitened. The network was trained on a single NVIDIA GeForce RTX 3090 GPU and took around 5 minutes per run. The condition number was computed via the explicit formula in Lemma 2 on CPU and took around 4 minutes. The condition number and the upper bound are shown in Figure 10 together with the training loss. We make the following two observations: 1. The condition number takes values in the range of 6 to 12 throughout training, which corresponds to a maximal change of around 50% of the condition number at initialization. This indicates that the condition number at initialization can be indicative of how the conditioning is along the optimization trajectory. 2. The upper bound remains tight throughout training, which highlights that it can provide reliable information about the condition number of the Gauss-Newton matrix throughout training. In Lemma 2 we have seen that the upper bound of the condition number can be expressed in terms of a convex combination of the condition number of the weight matrices. In this section, we will elucidate the 'self-balancing' behavior of the individual terms in Lemma 2. For completeness, let us restate Lemma 2 here again: Lemma 2. Assume that m > max{d, k} and that α ℓ := σ 2 min (W L:ℓ+1 ) • σ 2 min (W 1:ℓ-1 ) > 0 ∀ℓ = 1, . . . , L. Let γ ℓ := α ℓ L i=1 αi . Then the condition number of the GN matrix of a L-layer linear network can be upper-bounded in the following way:κ( G O ) ≤ κ(Σ) L ℓ=1 γ ℓ κ(W L:ℓ+1 ) 2 κ(W 1:ℓ-1 ) 2 ≤ κ(Σ) max 1≤ℓ≤L κ(W L:ℓ+1 ) 2 κ(W 1:ℓ-1 ) 2 .By looking at the individual values of κ 2 (W L:ℓ+1 ) and σ 2 min (W L:ℓ+1 ) and similarly κ 2 (W 1:ℓ-1 ) and σ 2 min (W 1:ℓ-1 ) over ℓ = 1, . . . , L for networks of increasing depth L, we can see the 'selfbalancing' behavior in Figure 11 for downsampled MNIST (d = 196) over 10 runs and m = 300 neurons per hidden layer. Note that this behavior is also consistent for different widths m, input dimension d, and output dimension k. By first looking at the first row, we observe that κ 2 (W 1:ℓ-1 ) (top left) dominates the product of each term in Lemma 2 and follows an exponential growth rule in ℓ (note the log scale of the y-axis). By looking at the second row, we see that at the same time, the smallest singular value of κ 2 (W 1:ℓ-1 ) (bottom left) is decreasing exponentially, leading to the 'self-balancing' behavior of the condition number. The same observation also holds for W L:ℓ+1 , but since the number of weight matrices decreases with ℓ this leads to a mirrored plot. The 'self-balancing' behavior can be seen below in Figure 12 for ease of reading, where the individual terms of Lemma 2 are plotted with and without the weighting factor γ ℓ . Consider the case of a one-hidden layer linear network F θ θ θ (x) = WVx, with x ∈ R d , W ∈ R k×m , and V ∈ R m×d and some corresponding label y ∈ R k . Under the assumption of the MSE loss, the functional Hessian is then given as:EQUATIONwhere Ω = E[δ x,y x ⊤ ] ∈ R k×d is the residual-input covariance matrix, with δ x,y := F θ θ θ (x)y being the residual. Since we would like to estimate the spectrum of H F , let us suppose that it has aneigenvector v = a ⊗ b c ⊗ d. Assuming that λ is the corresponding eigenvalue, we have the following eigenproblem:H F • v = Ω c ⊗ d Ω ⊤ a ⊗ b = λ a ⊗ b c ⊗ dLet's start with the guess that b = d, then we obtain the following equations:λ a = Ω c and λ c = Ω ⊤ a.Solving this gives:Ω ⊤ Ω c = λ 2 c and ΩΩ ⊤ a = λ 2 a.Hence, a and c are simply the eigenvectors of ΩΩ ⊤ and Ω ⊤ Ω respectively. Both these have the same non-zero eigenvalues, and the resulting eigenvalues for the functional Hessian are nothing but the square root of corresponding eigenvalues of the matrix Ω ⊤ Ω -with both positive and negative signs. About b = d: We still have a degree of freedom left in choosing them. Any m-orthogonal vectors with unit norm would satisfy the above equation and, for simplicity, we can just pick the canonical basis vectors in R m .Teacher-Student setting Let us now further assume a Teacher-student setting, i.e. y = Zx for some Z ∈ R k×d . In this case note thatΩ = E[(WVx -y)x ⊤ ] = (WV -Z)E[xx ⊤ ] = (WV -Z)Σ.Now recall the definition of the condition number κ(H F ) = ∥H F ∥•∥H -1 F ∥ and using the observation from above we have that∥H F ∥ = σ max (H F ) = σ max ((WV -Z)Σ) ≤ σ max (WV -Z) • λ max (Σ) ∥H -1 F ∥ = σ min (H F ) = σ min ((WV -Z)Σ) ≥ σ min (WV -Z) • λ min (Σ).Therefore we haveEQUATION)Similar to bounding the condition number of the GN matrix, we can also clearly see a dependency on the conditioning of the input covariance.G Concrete bounds at initialization for one-hidden layer linear networkBased on the upper bound in Eq. ( 6) derived for the one-hidden layer linear network, we can look at how this bound behaves concretely at initialization. If we assume the entries of W and V to be i.i.d. Gaussian entries, then WW ⊤ and V ⊤ V will be Wishart matrices, for which asymptotically the distribution of the smallest and largest eigenvalues are known. Also there are known non-asymptotic bounds to bound the maximum and minimum singular value of a Gaussian matrix: Lemma 6 (Corollary 5.35 in Vershynin [2010] ). Let A be an N × n matrix whose entries are independent standard normal random variables N (0, 1). Then for every t ≥ 0, with probability at least 1 -2 exp(-EQUATIONUsing the Lemma above and a union bound, we can also derive a non-asymptotic bound for the condition number of the GN matrix under the assumption that we have a wide layer m > max(d, k):Lemma 7. Assume that m ≥ max(d, k) and that the entries of V ∈ R m×d , W ∈ R k×m are initialized Gaussian i.i.d. with V ij ∼ N (0, σ 2 v ) and W ij ∼ N (0, σ 2 w ), then with probability at least 1 -8 exp(-t 22 ) the condition number of the GN matrix κ( G O ) is upper bounded byκ( G O ) = λ max ( G O ) λ min ( G O ) ≤ κ(Σ) (σ 2 w ( √ m + √ k + t) 2 + σ 2 v ( √ m + √ d + t) 2 ) (σ 2 w ( √ m - √ k -t) 2 + σ 2 v ( m - √ d -t) 2 ) . Proof. Assume that V ∈ R m×d , W ∈ R k×m are initialized Gaussian i.i.d. with V ij ∼ N (0, σ 2 v ) and W ij ∼ N (0, σ 2w ), then the bounds in Lemma 6 are scaled with σ v and σ w respectively:σ v ( √ m - √ d -t) ≤ σ min (V) ≤ σ max (V) ≤ σ v ( √ m + √ d + t) σ w ( √ m - √ k -t) ≤ σ min (W) ≤ σ max (W) ≤ σ w ( √ m + √ k + t).Further noting thatλ min (V ⊤ V) = σ 2 min (V), λ max (V ⊤ V) = σ 2 max (V) λ min (WW ⊤ ) = σ 2 min (W) λ max (WW ⊤ ) = σ 2 max (W),we get the following non-asymptotic upper respectively lower bounds for the extreme eigenvalues with probability at least 1 -2 exp(-t 2 2 ) each:λ min (V ⊤ V) ≥ σ 2 v ( √ m - √ d -t) 2 , λ max (V ⊤ V) ≤ σ 2 v ( √ m + √ d + t) 2 λ min (WW ⊤ ) ≥ σ 2 w ( √ m - √ k -t) 2 , λ max (WW ⊤ ) ≤ σ 2 w ( √ m + √ k + t) 2Using a union bound, we get the upper bound for the condition number of the GN matrix.Based on the bound derived above, we can understand how the condition number behaves under initialization schemes used in practice. Consider the setting of a wide hidden layer m ≫ max{d, k}, in which the input dimension is much larger than the output dimension, i.e. d ≫ k. This is for instance the case for image classification, where the number of image pixels is typically much larger than the number of image classes. A popular initialization scheme is the one proposed in He et al. [2015] , commonly known as the Kaiming/He initialization, which corresponds to σ 2 w = 1 m , σ 2 v = 1 d in our case. Hence, under the assumption of m ≫ d, we will have σ 2 v ≫ σ 2 w and thereforeEQUATIONTo complement the analysis, we can also derive an asymptotic bound for the condition number when d, k, m go to infinity together.For this, we will use a result from Silverstein [1985] , in which it is shown that the smallest eigenvalue of the matrix M s = 1 s U s U ⊤ s , where U s ∈ R n×s composed of i.i.d. N (0, 1) random variables, converges to (1 -√ y) 2 as s → ∞, while n s → y ∈ (0, 1) as s → ∞. This corresponds to the case in which the hidden layer is larger than the input and output dimensions, that is m ≥ max(d, k). Similarly it can be shown that λ max → (1 + √ y) 2 for all y > 0 for s → ∞.Using this result we can bound the extreme eigenvalues of WW ⊤ and V ⊤ V at initialization. Assuming that the entries of W and V are entry-wise independent normally distributedW ij ∼ N (0, σ 2 w ), V ij ∼ N (0, σ 2 v ) (for instance σ 2 w = 1 m , σ 2 v = 1 d for Kaiming-initialization) we get for d, k, m → ∞: H Proofs Lemma 1. Let β w = σ 2 min (W)/(σ 2 min (W) + σ 2 min (V)).Then the condition number of GN for the one-hidden layer network with linear activations with m > max{d, k} is upper bounded byκ( G O ) ≤ κ(Σ) • σ 2 max (W) + σ 2 max (V) σ 2 min (W) + σ 2 min (V) = κ(Σ) • (β w κ(W) 2 + (1 -β w ) κ(V) 2 ). (6)Proof. Note that we have by the Weyl's in dual Weyl's inequality thatλ max ( G O ) = λ max (WW ⊤ ⊗ Σ + I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 ) ≤ λ max (WW ⊤ ⊗ Σ) + λ max (I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 ) λ min ( G O ) = λ min (WW ⊤ ⊗ Σ + I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 ) ≥ λ min (WW ⊤ ⊗ Σ) + λ min (I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 ).Then note that WW ⊤ ⊗ Σ and I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 are positive semidefinite and using the equality for the extreme eigenvalues of Kronecker products in Eq. ( 13), we have the boundλ max (WW ⊤ ⊗ Σ) + λ max (I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 ) = λ max (WW ⊤ )λ max (Σ) + λ max (Σ 1/2 V ⊤ VΣ 1/2 ) λ min (WW ⊤ ⊗ Σ) + λ min (I k ⊗ Σ 1/2 V ⊤ VΣ 1/2 ) = λ min (WW ⊤ )λ min (Σ) + λ min (Σ 1/2 V ⊤ VΣ 1/2 ).Using the submultiplicativity of the matrix norm we get the upper boundκ( G O ) ≤ κ(Σ) • λ max (WW ⊤ ) + λ max (V ⊤ V) λ min (WW ⊤ ) + λ min (V ⊤ V) .Finally, by noting that we haveλ i (WW ⊤ ) = σ 2 i (W) and λ i (V ⊤ V) = σ 2 i (V) for i ∈ {min, max} because we have m > max{d, k} yields the result.Lemma 4. Consider an one-hidden layer network with a Leaky-ReLU activation with negative slope α ∈ [0, 1]. Then the condition number of the GN matrix defined in Equation (10) can be bounded as:EQUATIONProof. As before we will bound the condition number by separately bounding the extreme eigenvalues first.First, recall the expression of the GN matrix asG O = m i=1 Λ i X ⊤ XΛ i ⊗ W •i W ⊤ •i + m i=1 Λ i X ⊤ V ⊤ i• V i• XΛ i ⊗ I kand note that we can apply the Weyl's and dual Weyl's inequality because(Λ i X ⊤ XΛ i )⊗(W •,i W ⊤ •,i ) and (Λ i X ⊤ V i,• V ⊤ i,• XΛ i ) ⊗ I kare symmetric, yielding the following bounds:λ min ( G O ) ≥ λ min m i=1 (Λ i X ⊤ XΛ i ) ⊗ (W •,i W ⊤ •,i ) + λ min m i=1 (Λ i X ⊤ V ⊤ i,• V i,• XΛ i ) ⊗ I k λ max ( G O ) ≤ λ max m i=1 (Λ i X ⊤ XΛ i ) ⊗ (W •,i W ⊤ •,i ) + λ max m i=1 (Λ i X ⊤ V ⊤ i,• V ⊤ i,• XΛ i ) ⊗ I k .Then we can bound the first term by using the fact that for the sum of a Kronecker product of PSD matrices A i , B i for i = 1, . . . , m we have λ min (EQUATIONby noting that Λ i only contains 1 or -α on its diagonals in case of Leaky-ReLU activation andm i=1 W •,i W ⊤ •,i = WW ⊤ . Analogously, we have λ max m i=1 (Λ i X ⊤ XΛ i ) ⊗ (W •,i W ⊤ •,i ) ≤ σ 2 max (X)σ 2 max (W).This yields the first upper bound.To get the second upper bound, we need to further bound the second term. For this, note thatEQUATIONAnalogously we haveEQUATIONwhich gives the second upper bound.We present here further details of the modified analysis of GD for strongly convex functions to study the effect of the condition number at initialization on the convergence rate, where we use local constants µ(k) and L(k) instead of the global smoothness and Lipschitz constant, respectively. Let L denote the Lipschitz constant and let the smoothness constant be denoted by µ. Furthermore, let the step size be such that η k ≤ 1 L . Then by the definition of gradient descent we get||θ k+1 -θ * || 2 = ||θ k -θ * -η k ∇f (θ k )|| 2 = ||θ k -θ * || 2 -2η k ∇f (θ k ) ⊤ (θ k -θ * ) + η 2 k ||∇f (θ k ) || 2 Strong convexity ≤ (1 -η k µ)||θ k -θ * || 2 -2η k (f (θ k ) -f (θ * )) + η 2 k ||∇f (θ k ) || 2 Smoothness ≤ (1 -η k µ)||θ k -θ * || 2 -2η k (f (θ k ) -f (θ * )) + 2η 2 k L(f (θ k ) -f (θ * )) = (1 -η k µ)||θ k -θ * || 2 -2η k (1 -η k L)(f (θ k ) -f (θ * ))Since we assumed that η k ≤ 1 L , the last term is negative. Therefore:EQUATIONSo by recursively applying ( 27) and replacing µ by the local smoothness constants µ(k):EQUATION)One can see the effect of µ(0) in the bound, which is even more dominant when µ(k) changes slowly.Of course, the effect of µ(0) attenuates over time, and that's why we are talking about a local effect.We repeated the experiments of pruning the weights at initialization for different architectures, including a small Vision Transformer (ViT)3 (Figure 13 ), a ResNet20 (Figure 14 ), a ResNet32 (Figure 15 ), a VGG5Net4 (Figure 16 ) and a Feed-forward network (Figure 17 ). The default initialization from PyTorch [Paszke et al., 2019] was used. In all setups, the weights were pruned layer-wise by magnitude and trained on a subset of Cifar-10 of n = 1000 images. The ViT was trained with AdamW with a learning rate of 1e -2 and weight decay 1e -2 . The ResNet20, ResNet32 and VGGnet were trained with SGD with momentum = 0.9 and weight decay of 10 -4 and a learning rate of 0.1 with a step decay to 0.01 after 91 epochs for ResNet20 and ResNet32 and after 100 epochs for the VGGnet. The Feed-forward network was trained with SGD with a constant learning rate of 0.01. All networks were trained with a mini-batch size of 64, the ResNet20, ResNet32 and the Feed-forward network were trained on a single NVIDIA GeForce RTX 3090 GPU and the ViT on a single GPU NVIDIA GeForce RTX 4090. The training time was around 1 hour for each setup. The computation of the condition number was run in parallel on GPU and took around 8 hours for each pruning rate of each network.For a one-hidden layer feed-forward network with ReLU activation, we empirically examined the effect of the hidden width on both the conditioning and the convergence speed when trained on a subset of MNIST (n = 1000). The networks were trained with SGD with a fixed learning rate, which was chosen via grid-search, on a single NVIDIA GeForce RTX 3090 GPU. The learning rate was 0.3 for the networks with width 15 and 20 and 0.5 for all the remaining networks. The networks were initialized via the Xavier normal initialization in PyTorch [Paszke et al., 2019] . The training time took less than one hour for all networks. The computation of the condition number was performed in parallel for each width and took around 4 hours in the largest setting. Figure 18 shows the trend that a larger width improves both the conditioning and the convergence speed. Additionally, the condition number is also more stable throughout training for larger widths.As we have seen previously in Figure 3 , the upper bound as a linear combination in Lemma 2 is crucial to get a tight bound on the condition number. Additionally, as illustrated in Figure 19 , whitening the data is essential to improve the conditioning. Therefore, if not otherwise specified, MNIST LeCun et al. [1998] and Cifar-10 Krizhevsky et al. [2009] will refer to whitened data and we refer to Lemma 2 or its equivalent in Eq. ( 7) for residual networks as the "upper bound". The shaded area in the figures corresponds to one standard deviation from the mean. If not other specified via the Kaiming normal initialization in PyTorch was used.As one can see from (8), the condition number for each term in the sum of the first upper bound in (8) improves when β increases. This can be seen by the fact that the ratio will be dominated by β and will go to 1 for β → ∞. This is also what we observe empirically in Figure 20 and Figure 22 , where the condition number is smaller for β = 1 compared to the other two settings, where β = 1 L < 1 √ L < 1 for deeper networks with L > 1.In Figure 23 the condition number of the GN matrix for a one-hidden layer network with Leaky-ReLU activation is plotted for different values of α, where α = 0 corresponds to ReLU and α = 1 corresponds to a linear activation. We can see how ReLU activation improves the conditioning in the one-hidden layer setting.The condition number was computed on CPU via the explicit formulas given in Lemma 2 and Eq. ( 7) and took approximately a few hours per experiment.The total compute for all experiments conducted amounts to approximately less than 100 hours. The main bottleneck is the computation of the condition number, which can be potentially further sped up by making full use of parallel computations on all available GPUs. The full research project required more compute than the experiments reported in the paper due to preliminary experiments. Nevertheless, the total compute of the full research project amounts to a run time of similar order of magnitude.In order to compute the condition number we need to compute the eigenspectrum of the GN matrix G O , which has dimension p × p, where p is the number of parameters or of the matrix ĜO , which has dimensions kd × kd, where d and k are the input, respectively output dimension of the network.The time complexity of calculating the eigenvalue decomposition has a cubic complexity. That is, in order to compute the condition number, we have a computational complexity of O(min(p, kd) 3 ). some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.For details of the derivation, we refer the reader to Appendix H.1The implementation of the ViT was based on Code taken from https://github.com/tintn/ vision-transformer-from-scratch/tree/main, which was published under the MIT license.The implementation of the ResNet20, ResNet32 and VGG5 was based on Code taken from https:// github.com/jerett/PyTorch-CIFAR10. The author granted explicit permission to use the code. 


Decision: This paper provides a theoretical characterization of the condition number of the Gauss-Newton (GN) matrix in neural networks, establishing tight bounds for linear networks and extending the analysis to some non-linear, residual, and convolutional architectures. This paper provides interesting analysis and will benefit the community, and therefore I recommend acceptance.

Reviews: review: paper_summary: paper_summary: This paper is dedicated to the theoretical characterization of the condition number of the Gauss-Newton (GN) matrix in neural networks. By studying deep linear networks and two-layer nonlinear networks, the authors establish tight bounds on the GN matrix's condition number and extend this analysis to architectures incorporating residual connections and convolutional layers. The methodology is rigorous, and the experimental validation is thorough, making significant contributions to understanding optimization processes in deep learning.
strength_weakness: strength_weakness: 1.The discussion in this paper on the relationship between the Gauss-Newton matrix condition number and the convergence rate of network optimization could be richer.

2.Convergence Rate Analysis in Figure 17: The network's optimal solution and the corresponding minimum loss differ under various settings, making it difficult to analyze the convergence rate from the loss changes. For instance, at epoch 300, the network with a width of 15 may have nearly converged, while the network with a width of 200 still shows a significant downward trend. Thus, it is hard to draw conclusions about the convergence speed. 

3.The caption of Figure 9 does not match the subfigures.

4.Figures 12-16 are not referenced or analyzed in the paper.

5.This paper lacks exploration of more general networks and recent network structures.
questions: questions: 1. Can this paper be extended to study the effect of network initialization on convergence rate?

2. Is it more appropriate to adjust the experiment on the effect of condition number and convergence rate in this paper to study different initializations or study the depth of the current network? They are more likely to keep the minimum loss close to each other, and make more use of the analysis of convergence rate.
limitations: limitations: As stated in the weaknesses.
paper_summary: paper_summary: This paper characterizes the conditioning of Gauss-Newton (GN) matrix. The contribution of this paper is clear and straightforward: for deep linear networks, it establishes a bound on the condition number of GN matrix, which is further extended to 2-layer ReLU networks. These bounds could be useful in certain scenarios. Numerical experiments are conducted to support the theoretical claims.
strength_weakness: strength_weakness: 1. My major concern is about the implications of the derived bounds. Specifically, I note that training/learning such as gradient descent learning dynamics, which is crucial in practice, is not involved in deriving the bound. Thus it is hard to see the implication of these derived bounds since they are the same for both before and after training of the model parameters, i.e., it provides no information about the effects of training and fails to characterize the properties of solutions of deep learning models in a specific task. This might even lead the bounds to be meaningless in certain scenarios. For example, in Lemma 1, let $k = 1$, then according to [1],  after enough iterations of gradient descent, the largest singular value $\sigma_{\max} \to \infty$ (as there is no formal definition of $\sigma_{\max}$ in Lemma 1 I assume that my understanding of its definition is correct) while $\sigma_{\min} \to 0$ for both $W$ and $V$, thus the bound blows up and becomes meaningless.

2. I find the introduction part of this paper a bit verbose, e.g., it spends about 4 pages (almost half of the main body) to discuss existing results and related works before presenting details of main results (starting from Lemma 1). I think it would be better to emphasize more about the technical contributions of the current work, which is not clear to me since it seems that many important steps have already been solved by previous works, e.g., Eq. (2), (3), and Lemma 3.

**Reference**

[1]  Ji and Telgarsky. Gradient descent aligns the layers of deep linear networks.
questions: questions: Could the authors give some implications of the proposed bounds and how we can better use it in practice?
limitations: limitations: The authors addressed their limitations in Section 7. In addition, in my view, the bounds fail to capture the distinctions between the solutions of deep learning models and random parameters such as the initialization, which limits its significance.
paper_summary: paper_summary: This paper examines the condition number of the Gauss-Newton matrix [1] in neural networks. It shows that normalization techniques, such as Batch Normalization [2], initial normalization, skip connections [3], and appropriate layer dimensions, reduce the condition number and therefore enhance the training stability. The objective is to provide insights into the training of deep neural networks by characterizing a new, tight upper bound for the condition number of the Gauss-Newton matrix. The paper primarily focuses on linear neural networks and shallow non-linear neural networks.

[1] SCHRAUDOLPH, Nicol N. Fast curvature matrix-vector products for second-order gradient descent. Neural computation, 2002, vol. 14, no 7, p. 1723-1738.
[2] IOFFE, Sergey et SZEGEDY, Christian. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In : International conference on machine learning. pmlr, 2015. p. 448-456.
[3] HE, Kaiming, ZHANG, Xiangyu, REN, Shaoqing, et al. Deep residual learning for image recognition. In : Proceedings of the IEEE conference on computer vision and pattern recognition. 2016. p. 770-778.
strength_weakness: strength_weakness: The paper is difficult to read and follow due to several typos and unclear ideas. The main contribution is not significant and heavily relies on Singh's work [4]. The authors discuss the K-FAC method [5], which uses the EFIM and not the Hessian, even though it can be an approximation in certain cases where the model's likelihood is in the exponential family [6]. The paper only addresses the Gauss-Newton matrix and not the Generalized Gauss-Newton matrix [6], which could be more relevant, especially when using a cross-entropy loss function instead of MSE. Additionally, the references are outdated (e.g., line 49 and most of the related work), and there are typographical errors, such as "generalized" being repeated on line 140 and "collapse" being misspelled on line 124. There is also a lack of references related to the Gauss-Newton matrix, Kronecker properties PyTorch library, etc...

[4] SINGH, Sidak Pal, BACHMANN, Gregor, et HOFMANN, Thomas. Analytic insights into structure and rank of neural network hessian maps. Advances in Neural Information Processing Systems, 2021, vol. 34, p. 23914-23927.
[5] MARTENS, James et GROSSE, Roger. Optimizing neural networks with kronecker-factored approximate curvature. In : International conference on machine learning. PMLR, 2015. p. 2408-2417.
[6] MARTENS, James. New insights and perspectives on the natural gradient method. Journal of Machine Learning Research, 2020, vol. 21, no 146, p. 1-76.
questions: questions: 1. Could you clarify the statement on line 31: “imagine an entire set of neurons being dead”? What implications does this scenario have for your study?
2. How good is the approximation of the Gauss Newton matrix $\mathbf{G}_O$ without considering the term $\mathbf{H}_F$?
3. Can you also conduct experiments on more challenging datasets like ImageNet [7] or TinyImageNet [8]?
4. Can you develop more about the time and space complexity of the computation of the condition number?
5. Can you explain and rephrase the lines 719 and 720 with the associated Figures 19 and 21?
6. Why did you compute the condition number on CPU rather than on GPUs (line 725)?
7. Can you also explore the relation between te condition number and the batch size?

[7] DENG, Jia, DONG, Wei, SOCHER, Richard, et al. Imagenet: A large-scale hierarchical image database. In : 2009 IEEE conference on computer vision and pattern recognition. Ieee, 2009. p. 248-255.
[8] Le, Ya and Xuan S. Yang. “Tiny ImageNet Visual Recognition Challenge.” (2015).
limitations: limitations: The limitations of the work are well addressed by the authors. I do not believe this work has any particular negative social impact.
paper_summary: paper_summary: The authors derived upper bounds for the condition number of the outer product of the jacobian of the neural network output (the Gauss Newton matrix) in the case of deep linear networks, and non-linear networks with a single hidden layer with piecewise-linear activation functions. They empirically evaluate these bounds using MNIST (with the exception of one plot in the main text that used CIFAR-10) and find their bounds reflect the general trend of the empirical condition number as depth (or in some cases width) are varied. The work is motivated by highlighting the importance of Hessian information in the optimisation of neural networks.
strength_weakness: strength_weakness: - ~~The paper generally uses the same dataset (MNIST) for almost all the plots, so it's hard to tell if these good behaviour of their bounds hold on other, potentially harder to optimise, problems such as CIFAR-10.~~ **It was pointed out in the rebuttal that CIFAR-10 was included in the appendix.**
- The precise motivation of studying the Gauss Newton is slightly weak. It's not clear to me that these results actually tell us anything substantial about optimising the neural network.
- The authors were not able to extend their bound to the deep non-linear network case.
questions: questions: - To be clear, on line 82 where we refer to "outer gradient product of the loss function", this is actually the outer gradient product of the network function and has no dependence on the loss function, right?
- Suggestion: In the equation below line 182 (which should probably be numbered) some extra brackets and the use of \cdots rather than \dots would make this expression clearer to read.
limitations: limitations: The authors briefly discuss some future directions their work could be taken in, but do not seriously critique the work they present here.

Comments: title: Official Comment by Reviewer 6bdM
Reply to author rebuttal
Answers to Questions 3-7
Answers to Question 1 and 2
comment: Thank you for your explanation, which has been generally helpful. Since my concerns have been resolved, I now hold a favorable view of this work and would like to raise my rating to 6.
Thank you for the authors' comprehensive rebuttal and the additional analyses provided in response to the feedback from myself and the other reviewers. 

The authors have well addressed my concerns, which has significantly improved my understanding of the paper. Based on this, I will raise my score to 6 (Weak Accept).
Thank you for your comprehensive response.

Yes, I see now that you have also run the main experiments on CIFAR-10, so this does address my concern. I wonder whether you could fit both MNIST and CIFAR-10 in figure 5, perhaps side by side? I believe the camera-ready version allows an extra page (please double check) so do consider this in case of acceptance.

I am convinced by your explanation of the motivation. On second inspection, it does appear that this is alluded to in the introduction, but perhaps some more explicit references to work like KFAC and Sophia might make the motivation clearer, as you said in your rebuttal.

I also agree that it is fair to leave the deep non-linear case to future work.

Thank you for clarifying my points of confusion. Please make modifications to wording if you believe it will help readability.

Overall, I think my initial score of 6 (weak accept) was slightly harsh. After realising that the experiments were also run on CIFAR-10, I am happy to raise my score to an 8 (strong accept) as I believe this is good work that has been carried out to an excellent standard.
**Response to rebuttal of weakness 1**

My concern still remains. Fig. 9 is not a practical setting: the model for Fig. 9 is a 3-layer linear model and cannot directly perform classification for data that is not linearly separable. The first plot of Fig. 9 only reveals that the loss converges, then 
1. if the model fits the data perfectly, then the data is linearly separable, which contradicts the second point of the author response.
2. if the model does not fit the data perfectly, then the first point of the author response is far from sufficient as the training dynamics is actually not a successful one. 

---
Overall, my point lies in that the proposed methods overlook many aspects of practical settings, therefore the range that the proposed methods can be applied to is rather limited.
**Question 3**:
> Can you develop more about the time and space complexity of the computation of the condition number?

**Answer**:
Thank you for this question. In order to compute the condition number we need to compute the eigenspectrum of the GN matrix $\mathbf{G}_O$, which has dimension $p \times p$, where $p$ is the number of parameters or of the matrix $\hat{\mathbf{G}}_O$, which has dimensions $kd \times kd$, where $d$ and $k$ are the input, respectively output dimension of the network. The time complexity of calculating the eigenvalue decomposition has a cubic complexity. That is, in order to compute the condition number, we have a computational complexity of $\mathcal{O}(\min(p, kd)^3)$.

**Question 4**:
> Can you also conduct experiments on more challenging datasets like ImageNet [7] or TinyImageNet [8]?

**Answer**:
Thank you for this suggestion. 
- **Only marginal gain in insight is expected from experiments on ImageNet.** We would like to argue that we expect only a marginal gain in insights from additional experiments on ImageNet compared to experiments on Cifar-10, which we have already conducted. 
In particular, note that the complexity of a given dataset appears only through the input covariance matrix $\boldsymbol{\Sigma}$, which is a separate factor in the upper bound. Furthermore, as we have mentioned in Remark R1, the effect of the conditioning of the input data on the conditioning of the GN spectra can be largely reduced by whitening or normalizing the input data. 
Thus, the effect that different datasets have on the condition number of the GN matrix can be largely alleviated through preprocessing, which is common practice.
- **Scaling up experiments is expensive and challenging.** An additional aspect that we want to emphasize is the additional effort to scale up the computation of the condition number, which is challenging and not straightforward, given the time and space complexity elaborated in the previous question.

Given the above explanations, we believe that experiments on ImageNet are less relevant to our work and currently also out of scope.

**Question 5**:
> Can you explain and rephrase the lines 719 and 720 with the associated Figures 19 and 21?

**Answer**:
Thank you for this question. As one can see from Eq. (6), the condition number for each term in the sum of the first upper bound in Eq.(6) improves when $\beta$ increases. This can be seen by the fact that the ratio will be dominated by $\beta$ and will go to 1 for $\beta \to \infty$. This is also what we observe empirically in Figure 19 an 21, where the condition number is smaller for $\beta = 1$ compared to the other two settings, where $\beta = 1/L < 1/\sqrt{L} < 1$ for deeper networks with $L > 1$.

**Question 6**:
> Why did you compute the condition number on CPU rather than on GPUs (line 725)?

**Answer**:
Thank you for this question. 
- In the case of linear and residual networks (with no activation function), where the GN matrix can be expressed analytically through Eq. (2), we did not find a significant time advantage of running the code on GPU compared to CPU.
- In the other case, where the GN matrix had to be computed numerically through automatic differentiation (for instance the experiments on pruning weights at initialization), we ran into memory problems with the GPUs. The way we resolved this was to build up the GN matrix through backpropagation on the GPU and move it to CPU to finally compute the condition number. Despite the slow down by moving the GN matrix between GPU and CPU, this still led to a speed up in time.

**Question 7**:
> Can you also explore the relation between the condition number and the batch size?

**Answer**:
Thank you for this interesting question! 

- Our work does not explicitly consider training dynamics, so generally speaking, our bounds consider $n$ to be the total number of sample points. However, one could also interpret $n$ as the batch size during training, in which case the condition number could be interpreted as the condition number of the "local loss landscape" of a single mini-batch. Note that $n$ only appears implicitly in our theoretical upper bounds through the condition number of the empirical input covariance matrix $\boldsymbol{\Sigma}$.

We kindly request the reviewer to let us know if there are any remaining
questions they may have. If they find that their queries have been sufficiently addressed, we would greatly appreciate it if they could reconsider their evaluation of our paper.
**Question 1**:
> Can this paper be extended to study the effect of network initialization on convergence rate?

**Answer**:
This is a very interesting question. As the condition number is a very local property, it is in general hard to connect the conditioning at network initialization to a global convergence rate. However, we would like to argue below that an ill-conditioned network initialization will still affect the rate of convergence for gradient descent (GD) in the initial phase of training. For this we will present a modified analysis of GD for strongly convex functions, where we use local constants $\mu(k)$ and $L(k)$ instead of the global smoothness and Lipschitz constant, respectively.
Let us denote the Lipschitz constant by $L$ and the smoothness constant by $\mu$. Furthermore, let the step size be such that $\eta_k \leq \frac{1}{L}$. Then by the definition of gradient descent, we have:

\begin{align*}
 ||\boldsymbol{\theta}\_{k+1} - \boldsymbol{\theta}^*  ||^2 &=  ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^* - \eta_k \nabla f(\boldsymbol{\theta}_k)  ||^2 \\\\
&=  ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^*  ||^2 -2 \eta_k \nabla f \left( \boldsymbol{\theta}_k \right)^\top \left( \boldsymbol{\theta}_k - \boldsymbol{\theta}^* \right) + \eta_k^2  ||\nabla f \left( \boldsymbol{\theta}_k \right)  ||^2  \\\\
& \stackrel{\text{Strong convexity}}{\leq} (1-\eta_k \mu)  ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^*  ||^2 - 2\eta_k (f(\boldsymbol{\theta}_k) - f(\boldsymbol{\theta}^*)) + \eta_k^2  ||\nabla f \left( \boldsymbol{\theta}_k \right)  ||^2  \\\\
& \stackrel{\text{Smoothness}}{\leq} (1-\eta_k \mu)  ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^*  ||^2 - 2\eta_k (f(\boldsymbol{\theta}_k) - f(\boldsymbol{\theta}^*)) + 2 \eta_k^2 L (f(\boldsymbol{\theta}_k) - f(\boldsymbol{\theta}^*)) \\\\
&= (1-\eta_k \mu)  ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^*  ||^2 - 2 \eta_k (1 - \eta_k L) (f(\boldsymbol{\theta}_k) - f(\boldsymbol{\theta}^*))
\end{align*}

Since we assumed that $\eta_k \leq \frac{1}{L}$, the last term is negative. Therefore:

\begin{equation}
 ||\boldsymbol{\theta}_{k+1} - \boldsymbol{\theta}^*  ||^2 \leq (1-\eta_k \mu)  ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^*  ||^2
\end{equation}


So by recursively applying Eq. (1) and replacing $\mu$ by the local smoothness constants $\mu(k)$:
\begin{equation}
 ||\boldsymbol{\theta}_k - \boldsymbol{\theta}^*  ||^2 \leq  \prod\_{i=0}^{k-1} (1-\eta_i \mu(i))  ||\boldsymbol{\theta}_0 - \boldsymbol{\theta}^*  ||^2
\end{equation}

One can clearly see the effect of $\mu(0)$ in the bound, which is even more dominant when $\mu(k)$ changes slowly. Of course, the effect of $\mu(0)$ attenuates over time, and that's why we are talking about a local effect. However, one should keep in mind that overparametrization leads the parameter to stay closer to initialization (at least in the NTK regime). We are happy to add a detailed discussion in the final version of the paper.

**Question 2**:
> Is it more appropriate to adjust the experiment on the effect of condition number and convergence rate in this paper to study different initializations or study the depth of the current network? They are more likely to keep the minimum loss close to each other, and make more use of the analysis of convergence rate.

**Answer**:
We are not sure what the Reviewer is asking precisely and kindly ask the Reviewer to rephrase the question, so we can dedicate an answer to it during the Author-Reviewer discussion phase.

We kindly request the reviewer to let us know if there are any remaining questions they may have. If they find that their queries have been sufficiently addressed, we would greatly appreciate it if they could reconsider their evaluationof our paper.