Deep linear networks for regression are implicitly regularized towards flat minima

Authors: pierre  marion
lénaïc  chizat
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 07:17:41

Abstract: The largest eigenvalue of the Hessian, or sharpness, of neural networks is a key quantity to understand their optimization dynamics. In this paper, we study the sharpness of deep linear networks for univariate regression. Minimizers can have arbitrarily large sharpness, but not an arbitrarily small one. Indeed, we show a lower bound on the sharpness of minimizers, which grows linearly with depth. We then study the properties of the minimizer found by gradient flow, which is the limit of gradient descent with vanishing learning rate. We show an implicit regularization towards flat minima: the sharpness of the minimizer is no more than a constant times the lower bound. The constant depends on the condition number of the data covariance matrix, but not on width or depth. This result is proven both for a small-scale initialization and a residual initialization. Results of independent interest are shown in both cases. For small-scale initialization, we show that the learned weight matrices are approximately rank-one and that their singular vectors align. For residual initialization, convergence of the gradient flow for a Gaussian initialization of the residual network is proven. Numerical experiments illustrate our results and connect them to gradient descent with non-vanishing learning rate.

Main paper: 
 introduction: 
 Neural networks have intricate optimization dynamics due to the non-convexity of their objective. A key quantity to understand these dynamics is the largest eigenvalue of the Hessian or sharpness S(W) (see Section 3 for a formal definition), in particular because of its connection with the choice of learning rate η. Classical theory from convex optimization indicates that the sharpness should remain lower than 2/η to avoid divergence (Nesterov, 2018) . The relevance of this point of view for deep learning has recently been questioned since neural networks have been shown to often operate at the edge of stability (Cohen et al., 2021) , where the sharpness oscillates around 2/η, while the loss still steadily decreases, albeit non-monotonically. Damian et al. (2023) explained the stability of gradient descent slightly above the 2/η threshold to be a general phenomenon for non-quadratic objectives, where the third-order derivatives of the loss induce a self-stabilization effect. They also show that, under appropriate assumptions, gradient descent on a risk R L implicitly solves the constrained minimization problemEQUATIONTrainability of neural networks initialized with a sharpness larger than 2/η is also studied by Lewkowycz et al. (2020) , which describes a transient catapult regime, which lasts until the sharpness goes below 2/η. These results beg the question of quantifying the largest learning rate that enables successful training of neural networks. For classification with linearly separable data and logistic loss, Wu et al. (2024) show that gradient descent converges for any learning rate. In this work, we address the case of deep linear networks for regression. As illustrated by Figure 1a , the picture then differs from the classification case: when the learning rate exceeds some critical value, the network fails to learn. We further remark that this critical value does not seem to be related to the initial scale: when the learning rate is under the critical value, learning is successful for a wide range of initial scales. In Figure 1b , we see that the large initialization scales correspond to initial sharpnesses well over the 2/η threshold, confirming that training is possible while initializing beyond the 2/η threshold. Figure 1 : Training a deep linear network on a univariate regression task with quadratic loss. The weight matrices are initialized as Gaussian random variables, whose standard deviation is the x-axis of plots 1a and 1b. Experimental details are given in Appendix C.To understand where the critical value for the learning rate comes from, we characterize the sharpness of minimizers of the empirical risk. We show that there exist minimizers with arbitrarily large sharpness, but not with arbitrarily small sharpness. Indeed, the sharpness of any minimizer grows linearly with depth, as made precise next. Theorem 1. Let X ∈ R n×d be a design matrix and y ∈ R n a target. Then the minimal sharpness S min of any linear network x → W L . . . W 1 x of depth L that implements the optimal linear regressor w⋆ ∈ R d satisfies 2∥w ⋆ ∥ 2-2 L 2 La ⩽ S min ⩽ 2∥w ⋆ ∥ 2-2 L 2 (2L -1)Λ 2 + (L -1) 2 a 2 ,where Λ is the largest eigenvalue of the empirical covariance matrix Σ := 1 n X ⊤ X, and a := (w ⋆ /∥w ⋆ ∥) ⊤ Σ(w ⋆ /∥w ⋆ ∥).We note that this bound is similar to a result of Mulayoff and Michaeli (2020) , although we alleviate their assumption of data whiteness ( Σ = I). In particular, we do not assume that data covariance matrix is full rank. Furthermore, our proof technique differs, since we do not require tensor algebra, and instead exhibit a direction in which the second derivative of the loss is large.This result shows that it is not possible to find a minimizer of the empirical risk in regions of low sharpness. Combined with (1), this suggests an interpretation of the critical value for the learning rate: gradient descent should not be able to converge to a global minimizer as soon asS min > 2 η ⇔ η > 2 S min ≃ (∥w ⋆ ∥ 2-2 L 2 La) -1 .This is confirmed experimentally by Figure 2 , which shows that the critical value of the learning rate matches our theoretical prediction. We note that this gives a quantitative answer to the observations of Lewkowycz et al. (2020) , which shows the existence of a maximal architecture-dependent learning rate beyond which training fails. The dependence of learning rate on depth (namely, constant over depth) also matches other papers that study scaling of neural networks (Chizat and Netrapalli, 2024; Yang et al., 2024) . La) -1 of Theorem 1 (dashed vertical line).To deepen our understanding of the training dynamics, we aim at quantifying the sharpness of the minimizer found by gradient descent (when it succeeds): we know that it has to be larger than S min , but is it close to S min or is it much larger? Inspecting Figure 1b , we see that the answer empirically depends on the interplay between initialization scale and learning rate. For small initialization scale, the sharpness after training is equal to a value relatively close to S min and independent of the learning rate. As the initialization scale increases, the sharpness of the trained network also increases, and plateaus at the value 2/η. The plateauing for large initialization scales can be explained by the edge of stability analysis (1), which upper bounds the sharpness of the minimizer by 2/η (see Figure 1d ). However, this gives no insight on the value of the sharpness when the learning rate is sufficiently small so that the network does not enter the edge of stability regime (see Figure 1c ).In this paper, we study the limiting case for vanishing η, i.e., training with gradient flow. As our main finding, we bound the sharpness of the minimizer found by gradient flow in the case of overdetermined regression, meaning that the sample size n is larger that the data dimension d and the data empirical covariance matrix is nonsingular. In this case, we prove that the ratio between the sharpness after training and S min is less than a constant depending mainly on the condition number of the empirical covariance matrix Σ. In particular, the ratio does not depend on the width or depth of the network. This shows an implicit regularization towards flat minima. Note that the phenomenon we exhibit is different from the well-studied implicit regularization towards flat minima caused by stochasticity in SGD (Keskar et al., 2017; Smith and Le, 2018; Blanc et al., 2020; Damian et al., 2021; Li et al., 2022; Liu et al., 2023) . In the present study, the dynamics are purely deterministic, and the low sharpness is due to the fact that the weight matrices found by gradient flow have (approximately) the same norm across layers, and that this norm is (approximately) the smallest possible one so that the network can minimize the risk.Link with generalization. Flatter minima have been found to generalize better (Hochreiter and Schmidhuber, 1997; Jastrzębski et al., 2017; Keskar et al., 2017; Jiang et al., 2020) , although the picture is subtle (Dinh et al., 2017; Neyshabur et al., 2017; Andriushchenko et al., 2023) . However, in this paper, we focus on the link of sharpness with (non-convex) optimization dynamics, rather than with generalization abilities. Indeed, our implicit regularization result holds in the overdetermined setting, where all minimizers of the empirical risk implement the same function thus have the same generalization error, although they differ in parameter space and in particular have different sharpnesses. We leave to future work extensions to more complex settings where our approach may link sharpness and generalization, beginning with deep linear networks for underdetermined regression. We refer to Appendix C for more comments on the link with generalization.We investigate two initialization schemes, quite different in nature: small-scale initialization and residual initialization. Let us explain both settings, by presenting our approach, contributions of independent interest, and related works.Small-scale initialization. In this setting, we consider an initialization of the weight matrices W k with i.i.d. Gaussian entries of small variance. Initialization scale is known to play a key role in training of neural networks: small-scale initialization corresponds to the "feature learning" regime where the weights change significantly during training, by opposition to the "lazy" regime (see, e.g., Chizat et al., 2019) . We show convergence of the empirical risk to zero, then characterize the structure of the minimizer found by gradient flow, a novel result of interest independently of its connection with sharpness. At convergence, the weight matrices are close to being rank-one, in the sense that all their singular values but the largest one are small. Furthermore, the first left singular vector of any weight matrix aligns with the first right singular vector of the next weight matrix. From this specific structure, we deduce our bound on the sharpness of the trained network. The bound is illustrated on Figure 1b , where our lower and upper theoretical bounds on the sharpness are plotted as dotted black lines. We observe that the sharpness after training, when starting from a small-scale initialization, is indeed situated between the black lines.The result and proof extend the study by Ji and Telgarsky (2020) for classification, although the parameters do not diverge to infinity contrarily to the classification case, thus requiring a finer control of the distance to the rank-one aligned solution. In regression, implicit regularization towards lowrank structure in parameter space was also studied by Saxe et al. (2014) ; Lampinen and Ganguli (2019) ; Gidel et al. (2019) ; Saxe et al. (2019) ; Varre et al. (2023) for two-layer neural networks and in Timor et al. (2023) for deep ReLU networks. This latter paper assumes convergence of the optimization algorithm and show that a solution with minimal ℓ 2 -norm has to be low-rank. In our linear setting, we instead show convergence. As detailed below, we impose mild requirements on the structure on the initialization beyond its scale; they are satisfied for instance by initializing one weight matrix to zero and the others with i.i.d. Gaussian entries. In particular, we do not require the so-called "zero-balanced initialization" as is common in the literature on deep linear networks (see, e.g., Arora et al., 2018; Advani et al., 2020; Li et al., 2021) or a deficient-margin initialization as in Arora et al. (2019a) . Finally, the limit when initialization scale tends to zero has been described for deep linear networks in Jacot et al. (2021) for multivariate regression. It consists in a saddle-to-saddle dynamics, where the rank of the weight matrices increases after each saddle. The present study considers instead a non-asymptotic setting where the initialization scale is small but nonzero, and shows convergence to a rank-one limit because univariate and not multivariate regression is considered.We note that sharpness at initialization can be made arbitrarily small, since it is controlled by the initialization scale, while sharpness after training scales as Θ(L). This therefore showcases an example of sharpening during training (although we make no statement on monotonicity).Residual initialization. Architectures of deep neural networks used in practice often present residual connections, which stabilize training (He et al., 2015) . A simple non-linear residual architecture writesh k+1 = h k + σ(N k h k ).Removing the non-linearity, we get h k+1 = (I + N k )h k , which prompts us to consider deep linear networks with square weight matrices W k ∈ R d×d that are initialized asEQUATIONwhere the N k ∈ R d×d are filled with i.i.d. standard Gaussian entries and s ⩾ 0 is a hyperparameter tuning the initialization scale. The scaling of the residual branch in 1/ √ d is common and corresponds for instance to the so-called Glorot and He initializations respectively from Glorot and Bengio (2010) and He et al. (2015) . It ensures that the variance of the residual branch is independent of the width d.Similarly, as studied in Arpit et al. (2019) ; Marion et al. (2022) , the scaling in 1/ √ L is the right one so that the initialization noise neither blows up nor decays when L → ∞. Note that, in practice, this scaling factor is often replaced by batch normalization (Ioffe and Szegedy, 2015) , which has been shown empirically to have a similar effect to 1/L scaling (De and Smith, 2020) .In this setting which we refer to as residual initialization, we show global convergence of the empirical risk. To our knowledge, it is the first time that convergence is proven for a standard Gaussian initialization of the residual network outside the large width d ⩾ n regime. Previous works considered either an identity initialization (Bartlett et al., 2018; Arora et al., 2019a; Zou et al., 2020) or a smooth initialization such that ∥W k+1 (0) -W k (0)∥ F = O(1/L) (Sander et al., 2022; Marion et al., 2024) . The extension to standard Gaussian initialization leverages sharp bounds for the singular values of the product of W k (0). Our main assumption, in alignment with the literature, is that the risk at initialization should not be larger than a constant (depending on Σ and s).We then show that the weights after training can be writtenW k (∞) = I + s √ Ld N k + 1 L θ k ,where the Frobenius norm of the θ k is bounded by a constant (depending only on s). This structure finally enables us to bound the sharpness of the trained network. Remark that, to connect this analysis with our discussion of sharpness in univariate regression, we add to the residual network a final fixed projection vector p ∈ R d , so that our neural network writes x → p ⊤ W L . . . W 1 x, but the proof of convergence also holds without this projection.Experimentally, we give in Appendix C plots in the residual case that are qualitatively similar to Figure 1 . The main difference is that the initial sharpness is less sensitive to the initialization scale s.Organization of the paper. Section 2 details our setting and notations. Section 3 studies the sharpness of minimizers of the empirical risk and proves Theorem 1. Dynamics of gradient flow starting from small-scale initialization and residual initialization are respectively presented in Sections 4 and 5. The Appendix contains proofs, additional plots, experimental details, and related works. 
setting: 
 Model. We consider linear networks of depth L from R d to R, which are linear mapsx → p ⊤ W L . . . W 1 x (3) parameterized by weight matrices W 1 , . . . , W L , where W k ∈ R d k ×d k-1 , d 0 = d and p ∈ R d Lis a fixed vector. This definition includes both fully-connected networks by setting d L = 1 and p = 1, and residual networks by settingd 1 = • • • = d L = d, the W k closeto the identity, and p to some fixed (potentially random) vector in R d . We let W = (W 1 , . . . , W L ) and w prod = W ⊤ 1 . . . W ⊤ L p. Given X ∈ R n×d a design matrix and y ∈ R n a target, we consider the empirical risk for regressionR L (W) := 1 n ∥y -XW ⊤ 1 . . . W ⊤ L p∥ 2 2 = 1 n ∥y -Xw prod ∥ 2 2 =: R 1 (w prod ) .The notations R L (W) and R 1 (w prod ) may seem redundant, but are actually practical to define gradients of the risk both with respect to a single matrix W k and to the product w prod . Let Σ := 1 n X ⊤ X the empirical covariance matrix, and λ and Λ respectively its smallest nonzero and largest eigenvalue. For now, we do not assume that Σ is full rank, so there is more than one solution to the regression problem min wprod∈R d R 1 (w prod ). We denote by w ⋆ ∈ R d the smallest norm solution, and we let R min = R 1 (w ⋆ ) be the minimum of R 1 (and R L ). In all the following, we assume that w ⋆ ̸ = 0. Note that, due to the overparameterization induced by the neural network, there exists an infinity of parameterizations of the mapping x → w ⋆⊤ x.Gradient flow. We consider that the neural network is trained by gradient flow on the empirical risk R L , that is, the parameters evolve according to the ordinary differential equationEQUATION)An application of the chain rule givesEQUATIONwithEQUATION)where the second equality is a consequence of ∇R 1 (w ⋆ ) = 0.Notations. For k ∈ {1, . . . L}, we denote respectively by σ k , u k , and v k the first singular value (which equals the ℓ 2 operator norm), the first left singular vector and the first right singular vector of W k . The ℓ 2 operator norm of a matrix M is denoted by ∥M ∥ 2 and its Frobenius norm by ∥M ∥ F . Its smallest singular value is denoted by σ min (M ). For a vector v, we let ∥v∥ 2 its Euclidean norm.Finally, for quantities that depend on the gradient flow time t, we omit for concision their explicit dependence on t when it is dispensable. 
estimates of the minimal sharpness of minimizers: 
 To define the sharpness of the model, we let D = L k=1 d k d k-1 and identify the space of parameters with R D , which amounts to stacking all the entries of the weight matrices in a large D-dimensional vector. Then the norm of the parameters seen as a D-dimensional vector can be related to the Frobenius norm of the matrices by∥W∥ 2 2 = L k=1 ∥W k ∥ 2 F .This allows us to define the Hessian of the risk H : R D → R D×D , and we denote by S(W) its largest eigenvalue for some parameters W, or sharpness. We note that there exists alternative definitions of the sharpness, but this one is the most relevant to study optimization dynamics. Our results are specific to this definition. The following result gives estimates on the minimal sharpness of minimizers of the empirical risk (and is a strictly stronger statement than Theorem 1).Theorem 2. Let S min = inf W∈arg min R L (W) S(W) and a := (w ⋆ /∥w ⋆ ∥) ⊤ Σ(w ⋆ /∥w ⋆ ∥). We haveS min ⩾ 2a∥w ⋆ ∥ 2-1 L 2 ∥p∥ 1 L L k=1 1 ∥W k ∥ F , and 2∥w ⋆ ∥ 2-2 L 2 ∥p∥ 2 L La ⩽ S min ⩽ 2∥w ⋆ ∥ 2-2 L 2 ∥p∥ 2 L (2L -1)Λ 2 + (L -1) 2 a 2 .The proof of the result relies on the following variational characterization of the sharpness as the direction of the highest change of the gradientEQUATION)Lower bounds are proven by considering well-chosen directions Wk , for instance Wk = (1+ξβ k )W k for the first lower bound. The upper bound is proven by constructing a specific minimizer and bounding its sharpness. The first lower bound shows that the sharpness of minimizers can be arbitrarily high if one of the matrices has a low-enough norm. More precisely, take any minimizer W = (W 1 , . . . W L ) and consider W C = (CW 1 , W 2 /C, W 3 , . . . , W L ), for some C > 0. Then W C is still a minimizer, andS(W C ) ⩾ 2a∥w ⋆ ∥ 2-1 L 2 ∥p∥ 1 L ∥W 2 /C∥ F = 2a∥w ⋆ ∥ 2-1 L 2 ∥p∥ 1 L C ∥W 2 ∥ F C→∞ ----→ ∞.The fact that a reparameterization of the network can lead to arbitrarily high sharpness is consistent with a similar result by Dinh et al. (2017) for two-layer ReLU networks.Note that the first lower bound is arbitrarily small for minimizers such that the norms ∥W k ∥ F are large. On the contrary, the second lower bound is uniform and asymptotically matches the upper bound when L → ∞: we have S min ∼ 2∥w ⋆ ∥ 2 2 La. As already noted by Mulayoff and Michaeli (2020) , the intuition behind the linear scaling of the bound with depth can be seen from a one-dimensional example: take f (x 1 , . . . , x L ) = L k=1 x i . Then an easy computation shows that the sharpness of f at (1, . . . , 1) is equal to L -1. This showcases a simple example where the output of f is constant with L while its sharpness grows linearly with L. 
analysis of gradient flow from small initialization: 
 In this section, we characterize the structure of the minimizer found by gradient flow starting from a small-scale initialization. The proof is inspired by the one of Ji and Telgarsky (2020) for linearlyseparable classification, with a finer analysis due to the finiteness of minimizers in our setting.We consider the model (3) with d L = 1 and p = 1. Denoting by R 0 the empirical risk when the weight matrices are equal to zero, we can state our assumption on the initialization.(A 1 ) The initialization satisfies that R L (W(0)) ⩽ R 0 and ∇R L (W(0)) ̸ = 0.It is satisfied for instance if one of the weight matrices W k is equal to zero at initialization while the others have i.i.d. Gaussian entries, so that R L (W(0)) = R 0 and ∇ k R L (W(0)) ̸ = 0 (almost surely).Linear networks trained by gradient flow possess the following remarkable property (Arora et al., 2018) that shall be useful in the remainder. Lemma 1. For any time t ⩾ 0 and any k ∈ {1, . . . , L -1},W ⊤ k+1 (t)W k+1 (t) -W ⊤ k+1 (0)W k+1 (0) = W k (t)W ⊤ k (t) -W k (0)W ⊤ k (0) .Define nowε := 3 max 1⩽k⩽L ∥W k (0)∥ 2 F + 2 L-1 k=1 ∥W k (0)W ⊤ k (0) -W ⊤ k+1 (0)W k+1 (0)∥ 2 .Note that ε only depends on the initialization, and can be made arbitrarily small by scaling down the initialization. The following key lemma connects throughout training three key quantities to ε. Lemma 2. The parameters following gradient flow satisfy for any t ⩾ 0 that• for k ∈ {1, . . . , L}, ∥W k (t)∥ 2 F -∥W k (t)∥ 2 2 ⩽ ε , • for j, k ∈ {1, . . . , L}, |σ 2 k (t) -σ 2 j (t)| ⩽ ε , • for k ∈ {1, . . . , L -1}, ⟨v k+1 (t), u k (t)⟩ 2 ⩾ 1 - ε σ 2 k+1 (t).The first identity of the Lemma bounds the sum of the squared singular values of W k (t), except the largest one. In other words, it quantifies how close W k (t) is to the rank-one approximation given by the first term in its singular value decomposition. The second statement bounds the distance between the spectral norms of any two weight matrices. The last bound quantifies the alignment between the first left singular vector of W k (t) and the first right singular vector of W k+1 (t). In particular, if ε is small and σ 2 k+1 (t) is of order 1, then v k+1 (t) and u k (t) are nearly aligned. We next use this Lemma to show that the neural network satisfies a Polyak-Łojasiewicz (PL) condition, which is one of the main tools to study non-convex optimization dynamics (Rebjock and Boumal, 2023) . A well-known result, recalled in Appendix A for completeness, shows that this implies exponential convergence of the gradient flow to a minimizer W SI of the empirical risk. Theorem 3. Under Assumption (A 1 ), the network satisfies the PL condition for t ⩾ 1, in the sense that there exists some µ > 0 such that, for t ⩾ 1,L k=1 ∇ k R L (W(t)) 2 F ⩾ µ(R L (W(t)) -R min ) .The proof leverages the structure of the gradient of the risk with respect to the first weight matrix, which relies on the linearity of the neural network and the fact that we consider a univariate output. More precisely, recall that, by (5),∇ 1 R L (W(t)) = (W L (t) . . . W 2 (t)) ⊤ d1×1 ∇R 1 (w prod (t)) ⊤ 1×d0 .Therefore the Frobenius norm of the gradient decomposes as the product of two vector norms∇ 1 R L (W(t)) 2 F = ∥W L (t) . . . W 2 (t)∥ 2 2 ∥∇R 1 (w prod (t))∥ 2 2 ⩾ 4λ∥W L (t) . . . W 2 (t)∥ 2 2 (R L (W(t)) -R min ), where the lower bound unfolds from a straightforward computation. The delicate step is to lower bound ∥W L (t) . . . W 2 (t)∥ 2 , which we approach by distinguishing depending on the magnitude of σ 1 (t) = ∥W 1 (t)∥ 2 . If σ 1 (t) is large, we use Lemma 2 to deduce both that all σ k (t) are large and then that the first singular vectors of successive weight matrices are aligned. This implies that the product of weight matrices has a large norm. To analyze the case where σ 1 (t) is small, we use Assumption (A 1 ) to bound away R L (W(t)) from R 0 for t ⩾ 1, and therefore w prod (t) from 0. The fact that w prod (t) cannot be too close to 0 while σ 1 (t) is small implies that ∥W L (t) . . . W 2 (t)∥ 2 is large. All in all, this allows us to lower bound ∥W L (t) . . . W 2 (t)∥ 2 , and the PL condition follows.To characterize the weights at the end of the training, we make the following assumption.(A 2 ) The data covariance matrix Σ is full rank, and we haveε ⩽ 1, ∥w ⋆ ∥ 2 ⩾ 1, and 32L √ ε ⩽ 1.The first statement ensures unicity of the minimizer w ⋆ of R 1 , and thus, given that the risk goes to 0, we have w prod → w ⋆ . The last condition means that the initialization has to be scaled down as the depth increases, so that ε = O(1/L 2 ). Intermediates conditions are technical. We can then show the following corollary.Corollary 1. Under Assumptions (A 1 )-(A 2 ), there exists T ⩾ 1, such that, for all t ⩾ T and k ∈ {1, . . . , L},∥w ⋆ ∥ 2 2 1/L ⩽ σ k (t) ⩽ 2∥w ⋆ ∥ 2 1/L .Together with Lemma 2, this result gives a precise description of the structure of the weights at the end of the gradient flow trajectory. Up to the small factor ε, the weights are rank-one matrices, with equal norms and aligned singular vectors. Since the product of weights aligns with w ⋆ , this means that the first right singular vector of W 1 has to align with w ⋆ , and then the weight matrices align with their neighbors in order to propagate the signal in the network.Combining this specific structure of the weights with the variational characterization (7) of the sharpness and the explicit formulas ( 5)-( 6) for the gradients, we derive the following upper bound on the sharpness of the found minimizer.Corollary 2. Under Assumptions (A 1 )-(A 2 ), the following bounds on the sharpness of the minimizer W SI hold:1 ⩽ S(W SI ) S min ⩽ 4 Λ λ .This result shows that the sharpness of the minimizer is close to S min in the sense that their ratio is bounded by a constant times the condition number of Σ. For example, in the case of white data ( Σ = I), we obtain that S(W SI ) ⩽ 4S min . 
analysis of gradient flow from residual initialization: 
 We now study the case of residual initialization. We consider a linear network of the form (3) withW k (t) = I + s √ Ld N k + 1 L θ k (t) ,where each matrix is a d×d matrix and the N k are filled with i.i.d. Gaussian entries N (0, 1). We refer to Section 1 for a discussion of the scaling factor in front of N k . Following the standard initialization of residual networks, we assume that the θ k are initialized to zero. Note that the scaling factor 1/L in front of the θ k has no impact on the dynamics; it is convenient for exposition and computations, since we show that, with this scaling factor, the θ k are of order O(1) after training.Before stating the main result of this section on the convergence of the gradient flow, recall that p ∈ R d is a fixed vector appended at the end of the residual network to project its output back in R. Theorem 4. There exist C 1 , . . . , C 5 > 0 depending only on s such that, if L ⩾ C 1 and d ⩾ C 2 , then, with probability at least1 -16 exp(-C 3 d) , if R L (W(0)) -R min ⩽ C 4 λ 2 ∥p∥ 2 2 Λ ,the gradient flow dynamics (4) converge to a global minimizer W RI of the risk. Furthermore, the minimizer W RI satisfiesEQUATIONTo our knowledge, Theorem 4 is the first result showing convergence of gradient flow for standard Gaussian initialization of residual networks without assuming overparameterization. The main requirement is that the loss at initialization be not too large, as is standard in the literature analyzing gradient flow for deep linear residual networks (Bartlett et al., 2018; Arora et al., 2019a; Zou et al., 2020; Sander et al., 2022) . Note that our bound on the loss at initialization does not depend on the width d, depth L, or sample size n. We emphasize that the same proof holds for multivariate regression, in the absence of the projection vector p. We focus here on univariate regression to connect the result with the analysis of sharpness for univariate regression in Section 3. Details on adaptation to multivariate regression are given in Appendix B.7. Finally, the precise dependence of C 1 to C 5 on s can be found in the proof.The proof is a refinement of the analysis for identity initialization of residual networks (Zou et al., 2020; Sander et al., 2022) . From the expression of the gradients (5), we can show that4Λ∥p∥ 2 2 ∥Π L:k+1 ∥ 2 2 ∥Π k-1:1 ∥ 2 2 (R(W) -R min ) ⩾ ∥∇ k R L (W)∥ 2 F ⩾ 4λ∥p∥ 2 2 σ 2 min (Π L:k+1 )σ 2 min (Π k-1:1 )(R(W) -R min ) , with Π L:k := W L . . . W k and Π k:1 := W k . . . W 1 . Letting t * = inf t ∈ R + , ∃k ∈ {1, . . . , L}, ∥θ k (t)∥ F > C 5 ,the crucial step is to lower bound σ 2 min (Π L:k+1 ) and σ 2 min (Π k-1:1 ) uniformly for t ∈ [0, t ⋆ ], in order to get a PL condition valid for t ∈ [0, t ⋆ ]. Then, the condition on the loss at initialization is used to prove that t ⋆ = ∞, thereby the PL condition holds for all t ⩾ 0. We deduce both convergence and the bound on the norm of θ RI k . The lower bound on σ 2 min (Π L:k+1 ) and σ 2 min (Π k-1:1 ) is straightforward in the case of an identity initialization. In the case of Gaussian initialization, the proof is more intricate, and leverages the following high probability bounds on the singular values of residual networks. Lemma 3. There exist C 1 , . . . , C 4 > 0 depending only on s such that, ifL ⩾ C 1 , d ⩾ C 2 , u ∈ [C 3 , C 4 L 1/4 ] ,then, with probability at least1 -8 exp - du 2 32s 2 , it holds for all θ such that max 1⩽k⩽L ∥θ k ∥ 2 ⩽ 1 64 exp(-2s 2 -4u) and all k ∈ {1, . . . , L} thatEQUATIONσ min I + s √ Ld N k + 1 L θ k . . . I + s √ Ld N 1 + 1 L θ 1 ⩾ 1 4 exp - 2s 2 d -u .The proof of this result goes in three steps. We first study the evolution of the norm of the activations across the layers of the residual network when θ = 0, and prove a high-probability bound by leveraging concentration inequalities for Gaussian and χ2 distributions. Then an ε-net argument allows to bound the singular values. Finally, the extension to θ in a ball around 0 is done via a perturbation analysis. The proof technique is related to previous works (Marion et al., 2022; Zhang et al., 2022) , but provides a crisper and sounder bound. More precisely, Marion et al. (2022) show a bound on the norm of the activations with a probability of failure that decays polynomially with the width d, which is not sufficient to apply the ε-net argument that requires an exponentially decreasing probability of failure. As for Zhang et al. (2022) , they provide a similar bound with the purpose of showing convergence of wide residual networks, however with a less sharp probability of failure that increases polynomially with depth.1 Finally, as previously, the dependence of C 1 to C 4 on s can be found in the proof.The characterization (8) of the minimizer in Theorem 4 allows to bound its sharpness, as made precise by the following corollary. It holds under the same assumptions and high-probability bound as the conclusion of Theorem 4. Corollary 3. Under the assumptions of Theorem 4, and if the data covariance matrix Σ is full rank, there exists C > 0 depending only on s such that the following bounds on the sharpness of the minimizer W RI hold:1 ⩽ S(W RI ) S min ⩽ C Λ λ .As for Corollary 2, the proof relies on the fact that the norms of the weight matrices are close to each other and to the smallest possible norm to minimize the risk. This result shows again an implicit regularization towards a low-sharpness minimizer. Experimental illustration connecting the result with gradient descent with non-vanishing learning rate is provided in Appendix C. 
conclusion: 
 This paper studies dynamics of gradient flow for deep linear networks on a regression task. For small-scale initialization, we prove that the learned weight matrices are approximately rank-one and that their singular vectors align. For residual initialization, convergence of the gradient flow for a Gaussian initialization is proven. In both cases, we obtain that the sharpness of the solution found by gradient flow is close to the smallest sharpness among all minimizers. Interesting next steps include studying the dynamics at any initialization scale, for non-vanishing learning rates, as well as extension to non-linear networks. We refer to Appendix C for additional comments and preliminary experimental results regarding possible extensions. 
 
 Appendix: Organization of the Appendix. Appendix A presents some useful preliminary lemmas. The proofs of the results of the main paper are presented in Appendix B, while additional plots, discussion, and experimental details are given in Appendix C. Finally, Appendix D discusses some additional related work.Lemma 4. For α > 0 and x ∈ [0, 1 /2],(1-x) α ⩾ 1 -2αx .For α > 0 and x > 0 such that αx ⩽ 1,(1 + x) α ⩽ 1 + 2αx .Proof. Regarding the first inequality of the Lemma, we have(1 -x) α = exp(α log(1 -x)) ⩾ exp(α(-2x)) ⩾ 1 -2αx ,where the first inequality holds for x ∈ [0, 1 /2]. The second inequality of the Lemma is proven by(1 + x) α = exp(α log(1 + x)) ⩽ exp(αx) ⩽ 1 + 2αx ,where the second inequality holds when αx ⩽ 1.Lemma 5. There exists an absolute constant C > 0 such that, for L ⩾ C and x ⩾ 1,L exp(- √ Lx) ⩽ 4 exp(-x) .Proof. For any x ∈ R,L exp(- √ Lx) ⩽ 4 exp(-x) ⇔ exp(( √ L -1)x) ⩾ L 4 .Then, for L ⩾ 1 and x ⩾ 1, we have exp((√ L -1)x) ⩾ 1 + ( √ L -1)x + 1 2 ( √ L -1) 2 x 2 = 1 + ( √ L -1)x + 1 2 (L + 1 -2 √ L)x 2 ⩾ ( √ L -1)x + L 4 + 1 2 ( L 2 + 1 -2 √ L)x 2 For L large enough, L 2 + 1 -2 √ L ⩾ 0. Thus, since x ⩾ 1, exp(( √ L -1)x) ⩾ ( √ L -1)x + L 4 + 1 2 ( L 2 + 1 -2 √ L)x = L 4 + ( L 4 - 1 2 )x ⩾ L 4 ,where the last inequality holds for L large enough. This concludes the proof.Lemma 6. Let h ∈ R d , N ∈ R d×d with i.i.d. standard Gaussian entries, andY 1 = ∥N h∥ 2 2 ∥h∥ 2 2 , Y 2 = h ⊤ N h ∥h∥ 2 2 , Then Y 1 ∼ χ 2 (d) and Y 2 ∼ N (0, 1) .Proof. We have, for i ∈ {1, . . . , d},(N h) i = d j=1 N ij h j .By independence of the (N ij ) 1⩽j⩽d , we deduce that (N h) i follows a N (0, ∥h∥ 2 ) distribution. Furthermore, by independence of the rows of N , the (N h) i are independent. ThusY 1 = 1 ∥h∥ 2 d i=1 (N h) 2 i follows a χ 2 (d) distribution. Moving on to Y 2 , we have Y 2 = 1 ∥h∥ 2 d i,j=1 N ij h i h j .Thus Y 2 follows a centered Gaussian distribution, and by independence of the (N ij ) 1⩽i,j⩽d , its variance is equal to1 ∥h∥ 4 d i,j=1 h 2 i h 2 j = 1 ,which concludes the proof.The next lemma shows that the PL condition implies exponential convergence of the gradient flow. It is a well-known fact (see, e.g., Rebjock and Boumal, 2023 , for an overview of similar conditions), proved here for completeness. Lemma 7. Let f : R D → R be a differentiable function lower bounded by f min ∈ R, and consider the gradient flow dynamics dx dt = -∇f (x(t)) .If f satisfies the Polyak-Łojasiewicz inequality for t ⩾ 0 ∥∇f (x(t))∥ 2 2 ⩾ µ(f (x(t)) -f min ) , then x(t) converges to a global minimizer x ∞ , and, for t ⩾ 0,f (x(t)) -f min ⩽ (f (x(0)) -f min )e -µt .Proof. By the chain rule,d dt f (x(t)) = ∇f (x(t)), dx dt = -∥∇f (x(t))∥ 2 2 .Plugging in the Polyak-Łojasiewicz inequality,d dt f (x(t)) ⩽ -µ(f (x(t)) -f min ) . Thus d dt (f (x(t)) -f min ) ⩽ -µ(f (x(t)) -f min ) .To solve this differential inequality, one can for instance use the comparison theorem (Michel Petrovitch, 1901) , which states that f (x(t)) -f min is smaller that the solution g of the initial value problemg(0) = f (x(0)) -f min , d dt g(t) = -µg(t) . This shows that f (x(t)) -f min ⩽ (f (x(0)) -f min )e -µt .To show convergence of the iterates, letF (t) = √ µ t 0 ∥∇f (x(s))∥ds + 2 f (x(t)) -f min .We have by the Polyak-Łojasiewicz inequalityF ′ (t) = ∥∇f (x(t))∥ √ µ - ∥∇f (x(t))∥ f (x(t)) -f min ⩽ 0 . Then 0 ⩽ t 0 dx ds ds = t 0 ∥∇f (x(s))∥ds ⩽ F (t) √ µ ⩽ F (0) √ µ ,showing that x(t) converges.Lemma 8. With the notation introduced in Section 2, the following identities holdsR L (W) = R min + 1 n ∥X(w ⋆ -w prod )∥ 2 2 , and 4λ(R L (W) -R min ) ⩽ ∥∇R 1 (w prod )∥ 2 2 ⩽ 4Λ(R L (W) -R min ) .Proof. We haveR L (W) = 1 n ∥y -Xw prod ∥ 2 2 = 1 n ∥y -Xw ⋆ + X(w ⋆ -w prod )∥ 2 2 = 1 n ∥y -Xw ⋆ ∥ 2 + 1 n ∥X(w ⋆ -w prod )∥ 2 2 + 2 n y -Xw ⋆ , X(w ⋆ -w prod ) = R min + 1 n ∥X(w ⋆ -w prod )∥ 2 2 + 2 n X ⊤ (y -Xw ⋆ ), w ⋆ -w prod ,where the scalar product is equal to zero because ∇R 1 (w⋆ ) = -2 n X ⊤ (y -Xw ⋆ ) = 0.This gives the first identity of the Lemma. Next, denoting π(w prod ) the projection of w prod on the orthogonal subspace to the kernel of X, we haveR L (W) -R min = 1 n ∥X(w ⋆ -π(w prod ))∥ 2 2 = 1 n (w ⋆ -π(w prod )) ⊤ X ⊤ X(w ⋆ -π(w prod )) ⩽ 1 n ∥w ⋆ -π(w prod )∥ 2 ∥X ⊤ X(w ⋆ -π(w prod ))∥ 2 .By the formula (6) for the gradient of R 1 ,∥∇R 1 (w prod )∥ 2 = 2 n ∥X ⊤ X(w ⋆ -w prod )∥ 2 = 2 n ∥X ⊤ X(w ⋆ -π(w prod ))∥ 2 ⩾ 2λ∥w ⋆ -π(w prod )∥ 2 ,where the last lower bound holds because both w ⋆ and π(w prod ) are in the orthogonal to the kernel of X. Plugging in the formula above, we obtain thatR L (W) -R min ⩽ 1 2λ ∥∇R 1 (w prod )∥ 2 • 1 2 ∥∇R 1 (w prod )∥ 2 = 1 4λ ∥∇R 1 (w prod )∥ 2 2 .Finally, to obtain the upper bound on the gradient, note that∥∇R 1 (w prod )∥ 2 2 = 4 n 2 ∥X ⊤ X(w ⋆ -w prod )∥ 2 2 = 4 n 2 (w ⋆ -w prod ) ⊤ X ⊤ XX ⊤ X(w ⋆ -w prod ) ⩽ 4Λ n (w ⋆ -w prod ) ⊤ X ⊤ X(w ⋆ -w prod ) = 4Λ(R L (W) -R min ) . Lemma 9. Take W 1 , . . . W L ∈ R d×d such that, for all k ∈ {1, . . . , L}, ∥W k . . . W 1 ∥ 2 ⩽ M , and σ min (W k . . . W 1 ) ⩾ m ,where M ⩾ 1 and m ∈ (0, 1). Then, for all θ such that max1⩽k⩽L ∥θ k ∥ 2 ⩽ m 2 4M 2 , letting Wk = W k + θ k L , we have ∥ Wk . . . W1 ∥ 2 ⩽ 2M , and σ min ( Wk . . . W1 ) ⩾ m 2 .Proof. First note that the assumptions imply that, for anyEQUATIONIt shall come in handy to extend this formula to the case where j = k, where we define the empty matrix product to be equal to the identity matrix, which has an operator norm of 1 ⩽ M m . Now, take any θ as in the Lemma and any h 0 ∈ R d . Let, for k ⩾ 0,h k = W k . . . W 1 h 0 , and hk = Wk . . . W1 h 0 . Then h k = W k h k-1 and hk = W k + θ k L hk-1 . Thus hk -h k = θ k L hk-1 + W k ( hk-1 -h k-1 ) .Since h0 -h 0 = 0, we get by recurrence that, for k ⩾ 1,EQUATIONFrom there, let us prove by recurrence thatEQUATION)This equation holds for k = 0 since M ⩾ 1. Next, assume that it holds up to a certain rank k -1, and let us prove it at rank k. From (10), we get thatEQUATION)Since M ⩾ 1 and m < 1, the bound on ∥θ j ∥ 2 from the assumptions of the Lemma implies in particular that ∥θ j ∥ 2 ⩽ m 2M . Utilizing this, as well as (9) and the recurrence hypothesis (11) up to rank k -1, we get∥ hk -h k ∥ 2 ⩽ 1 L k j=1 M m • m 2M • 2M ∥ h0 ∥ 2 ⩽ M ∥ h0 ∥ 2 . Then ∥ hk ∥ 2 ⩽ ∥h k ∥ 2 + ∥ hk -h k ∥ 2 ⩽ ∥W k . . . W 1 ∥ 2 ∥h 0 ∥ 2 + M ∥ h0 ∥ 2 ⩽ 2M ∥ h0 ∥ 2 .This concludes the proof of the recurrence hypothesis at rank k. We therefore get that∥ Wk . . . W1 ∥ 2 ⩽ 2M .To prove the lower bound on the smallest singular value of Wk . . . W1 , observe that∥ hk ∥ 2 ⩾ ∥h k ∥ 2 -∥ hk -h k ∥ 2 ⩾ σ min (W k . . . W 1 )∥h 0 ∥ 2 - 1 L k j=1 ∥W k . . . W j+1 ∥ 2 ∥θ j ∥ 2 ∥ hj-1 ∥ 2 ,by ( 12). Finally, by ( 9), the bound on ∥θ j ∥ 2 from the assumptions, and the upper bound we just proved on ∥ hj-1 ∥ 2 , we have∥ hk ∥ 2 ⩾ m∥h 0 ∥ 2 - 1 L k j=1 M m • m 2 4M 2 • 2M ∥h 0 ∥ 2 ⩾ m 2 ∥h 0 ∥ 2 .This concludes the proof.For a twice continuously differentiable function f : R D → R, the largest eigenvalue S of its Hessian H(f ) : R D → R D×D at some x ∈ R D admits the variational characterizationS(x) = lim ξ→0 sup ∥x-x∥⩽ξ ∥∇f (x) -∇f (x)∥ 2 ∥x -x∥ 2 .In our case, the parameters are a set of matrices and the formula above translates intoEQUATIONWe now take W to be an arbitrary minimizer of R L . To obtain the lower bounds, consider for ξ ⩽ 1 Wk (ξ) = W k + ξM k ,where the M k ∈ R d k ×d k-1 are parameters that will be chosen later (depending on the W k but not on ξ). ThenEQUATIONTo alleviate notations, we drop the dependence of Wk on ξ. Recall that, for any parameters W 0 ,EQUATION)Minimizers of the empirical risk also satisfy that X ⊤ Xw prod = X ⊤ Xw ⋆ . Thus, by adding and subtracting differences,EQUATIONHere, and in the remainder of this proof, the notation O is taken with respect to the limit when ξ → 0, everything else being fixed. In particular, inspecting the expression above, we see that ∇R( wprod ) = O(ξ), and therefore, going back to (15), thatEQUATIONBy the inequality of arithmetic and geometric means, and by subadditivity of the operator norm,EQUATION)By definition of ∆k ,∆L • • • ∆1 = (-1) L p∇R 1 ( wprod ) ⊤ w prod ∇R 1 ( wprod ) ⊤ • • • w prod ∇R 1 ( wprod ) ⊤ L-1 times = (-1) L p(∇R 1 ( wprod ) ⊤ w prod ) L-1 ∇R 1 ( wprod ) ⊤ = (-1) L p(∇R 1 ( wprod ) ⊤ w ⋆ ) L-1 ∇R 1 ( wprod ) ⊤ ,where the last identity comes from the formulas ∇R 1 ( wprod ) = 2 n X ⊤ X( wprod -w prod ) andX ⊤ Xw prod = X ⊤ Xw ⋆ . Thus ∥ ∆L • • • ∆1 ∥ 2 = (∇R 1 ( wprod ) ⊤ w ⋆ ) L-1 ∥p∥ 2 ∥∇R 1 ( wprod )∥ 2 ⩾ (∇R 1 ( wprod ) ⊤ w ⋆ ) L-1 ∥p∥ 2 ∇R 1 ( wprod ) ⊤ w ⋆ ∥w ⋆ ∥ 2 = (∇R 1 ( wprod ) ⊤ w ⋆ ) L ∥p∥ 2 ∥w ⋆ ∥ 2 ,and, by (18),EQUATION)Coming back to (16), we havew ⋆⊤ ∇R 1 ( wprod ) = ξ L k=1 w ⋆⊤ 2 n X ⊤ XW ⊤ 1 . . . W ⊤ k-1 M ⊤ k W ⊤ k+1 . . . W ⊤ L p + O(ξ 2 ) .At this point, the computations diverge for the two lower bounds we want to prove. For the first inequality, we take M k = β k W k , where the β k are free parameters to be optimized. We havew ⋆⊤ ∇R 1 ( wprod ) = ξ L k=1 β k w ⋆⊤ 2 n X ⊤ XW ⊤ 1 . . . W ⊤ k-1 W ⊤ k W ⊤ k+1 . . . W ⊤ L p + O(ξ 2 ) = ξ L k=1 β k w ⋆⊤ 2 n X ⊤ Xw prod + O(ξ 2 ) = ξ L k=1 β k w ⋆⊤ 2 n X ⊤ Xw ⋆ + O(ξ 2 ) = 2ξa∥w ⋆ ∥ 2 L j=k β k + O(ξ 2 ) .Therefore, by ( 19),L k=1 ∥∆ k ∥ 2 F ⩾ 4Lξ 2 a 2 ∥w ⋆ ∥ 4-2 L 2 ∥p∥ 2 L 2 L k=1 β k 2 + O(ξ 3 ) . Furthermore, L k=1 ∥M k ∥ 2 F = L k=1 β 2 k ∥W k ∥ 2 F .By ( 14), we getS(W) 2 ⩾ 4La 2 ∥w ⋆ ∥ 4-2 L 2 ∥p∥ 2 L 2 L k=1 β k 2 L k=1 β 2 k ∥W k ∥ 2 F .The first lower bound unfolds by takingβ k = 1/∥W k ∥ F .Moving on to the second lower-bound, we now takeM k = β k u k v ⊤ k ∥u k ∥ 2 ∥v k ∥ 2 , u k = W k-1 . . . W 1 2 n X ⊤ Xw ⋆ , v k = W ⊤ k+1 . . . W ⊤ L p ,where again the β k are free parameters to be optimized. We therefore havew ⋆⊤ ∇R 1 ( wprod ) = L k=1 w ⋆⊤ 2 n X ⊤ XW ⊤ 1 . . . W ⊤ k-1 =u ⊤ k M ⊤ k W ⊤ k+1 . . . W ⊤ L p =v k +O(ξ 2 ) = ξ L k=1 β k ∥u k ∥ 2 ∥v k ∥ 2 + O(ξ 2 ) = ξ L k=1 β k ∥W ⊤ k+1 . . . W ⊤ L p∥ 2 ∥w ⋆⊤ 2 n X ⊤ XW ⊤ 1 . . . W ⊤ k-1 ∥ 2 + O(ξ 2 ) ⩾ ξ L k=1 β k ∥W ⊤ k+1 . . . W ⊤ L pw ⋆⊤ 2 n X ⊤ XW ⊤ 1 . . . W ⊤ k-1 ∥ 2 + O(ξ 2 ) ,where the last line unfolds from subadditivity of the operator norm.We letA k = W ⊤ k+1 . . . W ⊤ L pw ⋆⊤ 2 n X ⊤ XW ⊤ 1 . . . W ⊤ k-1 , and choose β k = ∥A k ∥ 2 . Then w ⊤ prod ∇R 1 ( wprod ) ⩾ ξ L k=1 ∥A k ∥ 2 2 + O(ξ 2 ) .Coming back to (19), we getL k=1 ∥∆ k ∥ 2 F ⩾ Lξ 2 ∥p∥ 2/L 2 ∥w ⋆ ∥ 2/L 2 L k=1 ∥A k ∥ 2 2 2 + O(ξ 3 ) . Furthermore, L k=1 ∥W k -Wk ∥ 2 F = ξ 2 L k=1 β 2 k ∥u k v ⊤ k ∥ 2 F ∥u k ∥ 2 2 ∥v k ∥ 2 2 = ξ 2 L k=1 β 2 k = ξ 2 L k=1 ∥A k ∥ 2 2 .Therefore, we obtain, by ( 14),S(W) 2 ⩾ L∥p∥ 2/L 2 ∥w ⋆ ∥ 2/L 2 L k=1 ∥A k ∥ 2 2 .We can lower bound the sum similarly to (18):EQUATIONA L • • • A 1 = pw ⋆⊤ 2 n X ⊤ X w prod w ⋆⊤ 2 n X ⊤ X • • • w prod w ⋆⊤ 2 n X ⊤ X L-1 times = p w ⋆⊤ 2 n X ⊤ Xw ⋆ L-1 w ⋆⊤ 2 n X ⊤ X .Thus, using the Cauchy-Schwarz inequality, we get∥A L • • • A 1 ∥ 2 = w ⋆⊤ 2 n X ⊤ Xw ⋆ L-1 ∥pw ⋆⊤ 2 n X ⊤ X∥ 2 = w ⋆⊤ 2 n X ⊤ Xw ⋆ L-1 ∥p∥ 2 ∥w ⋆⊤ 2 n X ⊤ X∥ 2 ⩾ ∥p∥ 2 ∥w ⋆ ∥ 2 w ⋆⊤ 2 n X ⊤ Xw ⋆ L = 2 L a L ∥p∥ 2 ∥w ⋆ ∥ 2L-1 2 . Therefore L k=1 ∥A k ∥ 2 2 ⩾ 4La 2 ∥w ⋆ ∥ 4-2 L 2 ∥p∥ 2 L 2 .We finally obtainS(W) 2 ⩾ 4L 2 a 2 ∥w ⋆ ∥ 4-4 L 2 ∥p∥ 4 L2 , which gives the second lower bound.To obtain the upper bound on S min , we construct explicitly a minimizer, and upper bound its sharpness. More precisely, let the W k be rank-one matrices, such that ∥W k ∥ 2 = ∥w ⋆ ∥ 1/L /∥p∥ 1/L , successive matrices have aligned first singular vectors, the first right singular vector of W 1 is aligned with w ⋆ , and the first left singular vector of W L is aligned with p. We then havep ⊤ W L . . . W 1 = w ⋆ ,meaning that the network minimizes the loss. We now upper bound its sharpness using (13), where we recall that the matrix appearing in the numerator is denoted by ∆ k and satisfies ∆ k = ∆k +O(ξ 2 ), where ∆k is given by (17). Contrarily to the proof of the lower bounds where we exhibited a specific direction W , we here seek an upper bound valid for all W . Recall that, by ( 16), we haveEQUATION)By subadditivity of the operator norm,EQUATION)Moving on to bounding the squared Frobenius norm of ∆ k , we observe that ∆k decomposes as a rank-one matrix. We split cases for k = 1 and k > 1. First, for k = 1, we have, by subadditivity of the operator norm and by ( 21),∥∆ 1 ∥ F = ∥W ⊤ 2 . . . W ⊤ L p∥ 2 ∥∇R 1 ( wprod )∥ 2 + O(ξ 2 ) ⩽ ∥W 2 ∥ 2 . . . ∥W L ∥ 2 ∥p∥ 2 ∥∇R 1 ( wprod )∥ 2 + O(ξ 2 ) ⩽ 2Λ∥w ⋆ ∥ 2(L-1)/L 2 ∥p∥ 2/L 2 L k=1 ∥ Wk -W k ∥ 2 + O(ξ 2 ) . Then ∥∆ 1 ∥ 2 F ⩽ 4Λ 2 ∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 L k=1 ∥ Wk -W k ∥ 2 2 + O(ξ 3 ) ⩽ 4LΛ 2 ∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 L k=1 ∥ Wk -W k ∥ 2 2 + O(ξ 3 ) ⩽ 4LΛ 2 ∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 L k=1 ∥ Wk -W k ∥ 2 F + O(ξ 3 ) .For k > 1, we haveEQUATION)Let us now bound ∥∇R 1 ( wprod ) ⊤ W ⊤ 1 ∥ 2 . By (20), separating the first term, we have∥∇R 1 ( wprod ) ⊤ W ⊤ 1 ∥ 2 ⩽ ∥p ⊤ W L . . . W 2 ( W1 -W 1 ) 2 n X ⊤ XW ⊤ 1 ∥ 2 + L k=2 ∥p ⊤ W L . . . W k+1 ( Wk -W k )W k-1 . . . W 1 2 n X ⊤ XW ⊤ 1 ∥ 2 + O(ξ 2 ) ⩽ 2Λ∥p∥ 2 ∥W L ∥ 2 . . . ∥W 2 ∥ 2 ∥ W1 -W 1 ∥ 2 ∥W 1 ∥ 2 + L k=2 ∥p∥ 2 ∥W L ∥ 2 . . . ∥W k+1 ∥ 2 ∥ Wk -W k ∥ 2 ∥W k-1 ∥ 2 . . . ∥W 2 ∥ 2 W 1 2 n X ⊤ XW ⊤ 1 2 + O(ξ 2 ) = 2Λ∥w ⋆ ∥ 2 ∥ W1 -W 1 ∥ 2 + L k=2 ∥w ⋆ ∥ (L-2)/L 2 ∥p∥ 2/L 2 W 1 2 n X ⊤ XW ⊤ 1 2 ∥ Wk -W k ∥ 2 + O(ξ 2 ) .Finally, recall that W 1 is rank-one and its first right singular vector is aligned with w ⋆ . A short computation therefore shows thatW 1 2 n X ⊤ XW ⊤ 1 2 = 2a∥W 1 ∥ 2 = 2a∥w ⋆ ∥ 2/L 2 ∥p∥ -2/L 2 . Thus ∥∇R 1 ( wprod ) ⊤ W ⊤ 1 ∥ 2 ⩽ 2∥w ⋆ ∥ 2 Λ∥ W1 -W 1 ∥ 2 + a L k=2 ∥ Wk -W k ∥ 2 + O(ξ 2 ) .Then, coming back to ( 22), for k > 1,∥∆ k ∥ F ⩽ 2∥w ⋆ ∥ 2(L-1)/L 2 ∥p∥ 2/L 2 Λ∥ W1 -W 1 ∥ 2 + a L k=2 ∥ Wk -W k ∥ 2 + O(ξ 2 ) .Thus∥∆ k ∥ 2 F ⩽ 4∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 Λ∥ W1 -W 1 ∥ 2 + a L k=2 ∥ Wk -W k ∥ 2 2 + O(ξ 3 ) ⩽ 4∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 (Λ 2 + (L -1)a 2 ) L k=1 ∥ Wk -W k ∥ 2 2 + O(ξ 3 ) ⩽ 4∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 (Λ 2 + (L -1)a 2 ) L k=1 ∥ Wk -W k ∥ 2 F + O(ξ 3 ) ,where the second inequality holds by the Cauchy-Schwarz inequality. Thus, by ( 13), putting together the bounds on ∥∆ k ∥ 2 F for k = 1 and k > 1,S(W) 2 = lim ξ→0 sup ∥W k -Wk ∥ F ⩽ξ L k=1 ∥∆ k ∥ 2 F L k=1 ∥W k -Wk ∥ 2 F ⩽ lim ξ→0 4∥w ⋆ ∥ 4(L-1)/L 2 ∥p∥ 4/L 2 (LΛ 2 + (L -1)(Λ 2 + (L -1)a 2 )) + O(ξ) . Therefore S(W) ⩽ 2∥w ⋆ ∥ 2-2 L 2 ∥p∥ 2 L 2 (2L -1)Λ 2 + (L -1) 2 a 2 ,which concludes the proof.This identity can be shown by noting that the identity is trivially true for t = 0, then differentiating on both sides with respect to time, and using (5). We refer, e.g., to Arora et al. (2018) for a detailed proof.Before proving the three statements of the lemma in order, we letε = max 1⩽k⩽L ∥W k (0)∥ 2 F + L-1 k=1 ∥W ⊤ k+1 (0)W k+1 (0) -W k (0)W ⊤ k (0)∥ 2 .Note that ε ⩽ ε.First statement. The claim is true for k = L since W L is a (row) vector. For k ∈ {1, . . . , L -1}, taking the 2-norm in Lemma 1, we have∥W ⊤ k+1 W k+1 ∥ 2 = ∥W k W ⊤ k + W ⊤ k+1 (0)W k+1 (0) -W k (0)W ⊤ k (0)∥ 2 ⩽ ∥W k W ⊤ k ∥ 2 + ∥W ⊤ k+1 (0)W k+1 (0) -W k (0)W ⊤ k (0)∥ 2 .Thus, usingEQUATION)We now take the trace in Lemma 1 to obtain∥W k+1 ∥ 2 F -∥W k+1 (0)∥ 2 F = ∥W k ∥ 2 F -∥W k (0)∥ 2 F .Combining with the inequality above, we have that∥W k ∥ 2 F -∥W k ∥ 2 2 ⩽ ∥W k+1 ∥ 2 F -∥W k+1 ∥ 2 2 + ∥W k (0)∥ 2 F -∥W k+1 (0)∥ 2 F + ∥W ⊤ k+1 (0)W k+1 (0) -W k (0)W ⊤ k (0)∥ 2Summing from k to L -1 and telescoping, we have∥W k ∥ 2 F -∥W k ∥ 2 2 ⩽ ∥W L ∥ 2 F -∥W L ∥ 2 2 + ∥W k (0)∥ 2 F -∥W L (0)∥ 2 F + L-1 k ′ =k ∥W ⊤ k ′ +1 (0)W k ′ +1 (0) -W k ′ (0)W ⊤ k ′ (0)∥ 2 .The first two terms compensate since W L is a vector, and the remainder of the terms is less than ε ⩽ ε by definition. This gives the first statement of the Lemma.Second statement. Assume without loss of generality that j > k. By recurrence over (23), we obtain thatσ 2 j ⩽ σ 2 k + j-1 k ′ =k ∥W ⊤ k ′ +1 (0)W k ′ +1 (0) -W k ′ (0)W ⊤ k ′ (0)∥ 2 ⩽ σ 2 k + ε .The reverse inequality can be shown similarly: by considering again Lemma 1 and taking the 2-norm, we get∥W k W ⊤ k ∥ 2 ⩽ ∥W ⊤ k+1 W k+1 ∥ 2 + ∥W ⊤ k+1 (0)W k+1 (0) -W k (0)W ⊤ k (0)∥ 2 Thus ∥W k ∥ 2 2 ⩽ ∥W k+1 ∥ 2 2 + ∥W ⊤ k+1 (0)W k+1 (0) -W k (0)W ⊤ k (0)∥ 2 .As previously, we get by recurrence that σ 2 k ⩽ σ 2 j + ε . Combined with the reverse bound above, this gives the second statement of the lemma.⊤ k W ⊤ k+1 W k+1 u k .We have on the one hand, by Lemma 1,u ⊤ k W ⊤ k+1 W k+1 u k = u ⊤ k W k W ⊤ k u k -u ⊤ k W k (0)W ⊤ k (0)u k + u ⊤ k W ⊤ k+1 (0)W k+1 (0)u k ⩾ u ⊤ k W k W ⊤ k u k -u ⊤ k W k (0)W ⊤ k (0)u k ⩾ σ 2 k -∥W k (0)∥ 2 2 . On the other hand, u ⊤ k W ⊤ k+1 W k+1 u k = u ⊤ k W ⊤ k+1 W k+1 -v k+1 σ 2 k+1 v ⊤ k+1 u k + u ⊤ k v k+1 σ 2 k+1 v ⊤ k+1 u k ⩽ ∥W ⊤ k+1 W k+1 -v k+1 σ 2 k+1 v ⊤ k+1 ∥ 2 + ⟨v k+1 , u k ⟩ 2 σ 2 k+1The 2-norm above is equal to the second largest eigenvalue of W ⊤ k+1 W k+1 , which is the square of the second largest singular value of W k+1 . In particular, it is lower than ∥W k+1 ∥ 2 F -∥W k+1 ∥ 2 2 which is the sum of the squared singular values of W k+1 except the largest one. We obtainu ⊤ k W ⊤ k+1 W k+1 u k ⩽ ∥W k+1 ∥ 2 F -∥W k+1 ∥ 2 2 + ⟨v k+1 , u k ⟩ 2 σ 2 k+1 ⩽ ε + ⟨v k+1 , u k ⟩ 2 σ 2 k+1by the first statement of the Lemma. Combining the lower and upper bound ofu ⊤ k W ⊤ k+1 W k+1 u k , we get ε + ⟨v k+1 , u k ⟩ 2 σ 2 k+1 ⩾ σ 2 k -∥W k (0)∥ 2 2 . Thus ⟨v k+1 , u k ⟩ 2 ⩾ σ 2 k -∥W k (0)∥ 2 2 - ε σ 2 k+1 ⩾ σ 2 k+1 -ε -∥W k (0)∥ 2 2 - ε σ 2 k+1 ,by the second statement of the Lemma. We finally obtain⟨v k+1 , u k ⟩ 2 ⩾ 1 - 2ε + ∥W k (0)∥ 2 2 σ 2 k+1 ⩾ 1 - ε σ 2 k+1, by definition of ε and ε.We lower bound the first term in the sum of the left-hand side. Recall that, by (5), we haveEQUATION)We show that ∥W L (t) . . . W 2 (t)∥ 2 is large by distinguishing between two cases depending on the magnitude of σ 1 (t) = ∥W 1 (t)∥ 2 . To this aim, let C > √ 2εL.Large spectral norm. We first consider the case whereEQUATION)By Lemma 2, for k ∈ {1, . . . , L -1},σ 2 k ⩾ σ 2 1 -ε > ε ,where the second inequality unfolds from (25). Then, again by Lemma 2, for k ∈ {1, . . . , L -1},⟨v k+1 , u k ⟩ 2 ⩾ 1 - ε σ 2 k+1 > 0 .It is always possible (without loss of generality) to choose the orientation of the u k and v k such that ⟨v k+1 , u k ⟩ ⩾ 0 for any k ∈ {1, . . . L -1}. Making this choice, the equation above implies thatEQUATION)Let us show that this implies a lower bound on ∥W L . . . W 2 ∥ 2 . To do so, let us denote recursivelyx 1 = v 2 , x k+1 = W k+1 x k for k ∈ {1, . . . , L -1} .We then havex L = W L . . . W 2 x 1 , thus ∥W L . . . W 2 ∥ 2 ⩾ ∥x L ∥ 2 ∥x 1 ∥ 2 = ∥x L ∥ 2 .Our goal is thus to lower bound ∥x L ∥ 2 , which entails a lower bound on ∥W L . . . W 2 ∥ 2 . To this aim, first note that, for any k ∈ {1, . . . , L -1},EQUATION)where the last equation stems from (26). Denoteα k = ⟨ x k ∥x k ∥2 , u k ⟩.Then the previous equation shows thatα k+1 ⩾ ∥x k ∥ 2 ∥x k+1 ∥ 2 (σ k+1 α k - √ 2ε) = ∥x k ∥ 2 σ k+1 ∥x k+1 ∥ 2 α k - √ 2ε σ k+1 . Further note that ∥x k+1 ∥ 2 ⩽ σ k+1 ∥x k ∥ 2 , thus α k+1 ⩾ α k - √ 2ε σ k+1.By recurrence,α k ⩾ α 2 - √ 2ε k-1 k ′ =2 1 σ k+1 ⩾ α 2 - √ 2ε k-1 k ′ =2 1 σ 2 1 -ε = α 2 - √ 2ε(k -2) σ 2 1 -ε .Coming back to (27), we have∥x k+1 ∥ 2 ⩾ ⟨x k+1 , u k+1 ⟩ ⩾ ∥x k ∥ 2 (σ k+1 α k - √ 2ε) , thus by recurrence, ∥x L ∥ 2 ⩾ ∥x 2 ∥ 2 L-1 k=2 (σ k+1 α k - √ 2ε) ⩾ ∥x 2 ∥ 2 L-1 k=2 σ 2 1 -ε α 2 - √ 2ε(k -2) σ 2 1 -ε - √ 2ε ⩾ ∥x 2 ∥ 2 L-1 k=2 σ 2 1 -εα 2 - √ 2εk .Finally, by definition,α 2 = ⟨ x2 ∥x2∥2 , u 2 ⟩. Since x 2 = W 2 x 1 = W 2 v 2 = σ 2 u 2, we obtain that α 2 = 1, and thus∥W L . . . W 2 ∥ 2 ⩾ ∥x L ∥ 2 ⩾ σ 2 L-1 k=2 σ 2 1 -ε - √ 2εk ⩾ σ 2 1 -ε σ 2 1 -ε - √ 2ε(L -1) L-2 ⩾ σ 2 1 - √ ε - √ 2ε(L -1) L-1 ⩾ σ 1 - √ 2εL L-1 .We finally getEQUATION) which is a positive quantity by definition of C.Small spectral norm. We now inspect the case where ( 25) is not satisfied, that is, σ 1 (t) ⩽ C. First note that, by the formula (6) for the gradient of R 1 ,∥∇R 1 (w prod )∥ 2 ⩽ 2Λ∥w ⋆ -w prod ∥ 2 ⩽ 2Λ(∥w ⋆ ∥ 2 + ∥w prod ∥ 2 ) .Thus, for ∥w prod ∥ 2 ⩽ ∥w ⋆ ∥ 2 , w prod → R 1 (w prod ) is 4Λ∥w ⋆ ∥ 2 -Lipschitz. Let us use this property to lower bound ∥w prod (t)∥ 2 by a constant independent of t, for t ⩾ 1. Either we have ∥w prod (t)∥ 2 ⩾ ∥w ⋆ ∥ 2 , or ∥w prod ∥ 2 ⩽ ∥w ⋆ ∥ 2 , but then, by the Lipschitzness property,|R 1 (w prod (t)) -R 0 | ⩽ 4Λ∥w ⋆ ∥ 2 ∥w prod (t) -0∥ 2 = 4Λ∥w ⋆ ∥ 2 ∥w prod (t)∥ 2 ,where we recall that R 0 is the risk associated to the null parameters. Furthermore, for t ⩾ 1,R 1 (w prod (t)) = R L (W(t)) ⩽ R L (W(1))(The risk is decreasing along the gradient flow)< R L (W(0)) (∇R L (W(0)) ̸ = 0 by Assumption (A 1 )) ⩽ R 0 .(by Assumption (A 1 ))Thus, for t ⩾ 1, |R 1 (w prod (t)) -R 0 | > R 0 -R L (W(1)) > 0 .To summarize, we proved that, for t ⩾ 1,∥w prod (t)∥ 2 ⩾ min R 0 -R L (W(1)) 4Λ∥w ⋆ ∥ 2 , ∥w ⋆ ∥ 2 > 0 ,where we recall that ∥w ⋆ ∥ 2 > 0 by assumption (see Section 2). Furthermore,∥w prod (t)∥ 2 ⩽ ∥W L (t) . . . W 2 (t)∥ 2 ∥W 1 (t)∥ 2 = ∥W L (t) . . . W 2 (t)∥ 2 σ 1 (t) .Then, for t ⩾ 1,EQUATION)Conclusion. Combining ( 28) and ( 29), we obtain that, for t ⩾ 1,∥W L (t) . . . W 2 (t)∥ 2 ⩾ min C - √ 2εL L-1 , 1 C min R 0 -R L (W(1)) 4Λ∥w ⋆ ∥ 2 , ∥w ⋆ ∥ 2 =: √ µ 1 ,where µ 1 > 0 by the proof above. Then, for t ⩾ 1, by ( 24),∇ 1 R L (W (t)) 2 F ⩾ µ 1 ∥∇R 1 (w prod (t))∥ 2 2 . By Lemma 8, ∥∇R 1 (w prod (t))∥ 2 2 ⩾ 4λ(R L (W(t)) -R min ) . Thus, taking µ = 4µ 1 λ > 0, for t ⩾ 1, L k=1 ∇ k R L (W(t)) 2 F ⩾ ∇ 1 R L (W(t)) 2 F ⩾ µ(R L (W(t)) -R min ) .We first show that Assumption (A 2 ) implies a number of estimates that are useful in the following. Since ∥w ⋆ ∥ 2 ⩾ 1, we haveEQUATION)Let us now prove the Corollary. We first note that Theorem 3 implies exponential convergence of the empirical risk to its minimum by Lemma 7. This also implies the (exponential) convergence of w prod to w ⋆ , since the covariance matrix X ⊤ X is full rank, so w ⋆ is the unique minimizer of R 1 , and, by Lemma 8,R L (W) -R min = 1 n ∥X(w ⋆ -w prod )∥ 2 2 ⩾ λ∥w ⋆ -w prod ∥ 2 2 ,Furthermore,∥w prod ∥ 2 = ∥W ⊤ 1 . . . W ⊤ L ∥ 2 ⩽ ∥W 1 ∥ 2 . . . ∥W L ∥ 2 ⩽ L max k=1 σ k L .Let us show that, for t large enough,EQUATION)If it were not the case, then∥w prod ∥ 2 ⩽ L max k=1 σ k L ⩽ ∥w ⋆ ∥ 2 2 1/L + ε L = ∥w ⋆ ∥ 2 2 1 + ε 2 ∥w ⋆ ∥ 2 1/L L ⩽ ∥w ⋆ ∥ 2 2 1 + 2Lε 2 ∥w ⋆ ∥ 2 1/L, where the last inequality holds by Lemma 4 sinceLε 2 ∥w ⋆ ∥ 2 1/L ⩽ 1 ⇔ Lε ⩽ ∥w ⋆ ∥ 2 2 1/L, which holds by ( 31). Then, we obtain∥w prod ∥ 2 ⩽ 3∥w ⋆ ∥ 2 4 , since 1 + 2Lε 2 ∥w ⋆ ∥ 2 1/L ⩽ 3 2 ⇔ 4Lε ⩽ ∥w ⋆ ∥ 2 2 1/L, which also holds by ( 31). The inequality ∥w prod ∥ 2 ⩽ 3∥w ⋆ ∥2 4 contradicts the fact that w prod converges to w ⋆ , thus proving (33). Then, we haveL max k=1 σ 2 k ⩾ ∥w ⋆ ∥ 2 2 2/L + 2ε ∥w ⋆ ∥ 2 2 1/L ⩾ ∥w ⋆ ∥ 2 2 2/L + ε ,where the last inequality holds by ( 30). Furthermore, by Lemma 2, for all k ∈ {1, . . . , L},σ 2 k ⩾ L max j=1 σ 2 j -ε .This brings the first inequality of the Corollary. We now show the second inequality of the Corollary. First note that, by Lemma 2,∥W k -σ k u k v ⊤ k ∥ 2 F = ∥W k ∥ 2 F -∥W k ∥ 22 ⩽ ε , since both quantities are equal to the sum of the squared singular values of W k except the largest one. Then, adding and substracting,∥W ⊤ 1 . . . W ⊤ L -σ 1 . . . σ L v 1 u ⊤ 1 . . . v L u ⊤ L ∥ 2 ⩽ L k=1 ∥W ⊤ 1 . . . W ⊤ k-1 (W ⊤ k -σ k u k v ⊤ k )σ k+1 . . . σ L v k+1 u ⊤ k+1 . . . v L u ⊤ L ∥ 2 ⩽ L k=1 ∥W 1 ∥ 2 . . . ∥W k-1 ∥ 2 ∥W ⊤ k -σ k u k v ⊤ k ∥ 2 σ k+1 . . . σ L ⩽ L k=1 σ 1 . . . σ k-1 ∥W ⊤ k -σ k u k v ⊤ k ∥ F σ k+1 . . . σ L ⩽ √ ε L k=1 j̸ =k σ j = √ ε L j=1 σ j L k=1 1 σ k ⩽ √ ε L j=1 σ j L min σ k ⩽ L √ ε ∥w ⋆ ∥ 2 1/L L j=1 σ j ,by the first inequality of the Corollary. Moreover, using again the first inequality of the Corollary and Lemma 2, we have that⟨v k+1 , u k ⟩ 2 ⩾ 1 - ε σ 2 k+1 ⩾ 1 - ε ∥w ⋆ ∥2 2 2/L . Since u L = 1, we deduce that ∥σ 1 . . . σ L v 1 u ⊤ 1 . . . v L u ⊤ L ∥ 2 = L k=1 σ k ∥v 1 ∥ 2 L-1 k=1 u ⊤ k v k+1 ⩾ L k=1 σ k 1 - ε ∥w ⋆ ∥2 2 2/L L-1 2 ⩾ L k=1 σ k 1 - (L -1)ε ∥w ⋆ ∥2 2 2/Lwhere in the last step we used Lemma 4 which is valid sinceε ∥w ⋆ ∥2 2 2/L ⩽ 1 2 ⇔ ε ⩽ 1 2 ∥w ⋆ ∥ 2 2 2/L, which holds by (32) . By the triangular inequality, we now haveL k=1 σ k 1 - (L -1)ε ∥w ⋆ ∥2 2 2/L ⩽ ∥σ 1 . . . σ L v 1 u ⊤ 1 . . . v L u ⊤ L ∥ 2 ⩽ ∥W ⊤ 1 . . . W ⊤ L ∥ 2 + ∥W ⊤ 1 . . . W ⊤ L -σ 1 . . . σ L v 1 u ⊤ 1 . . . v L u ⊤ L ∥ 2 ⩽ ∥W ⊤ 1 . . . W ⊤ L ∥ 2 + L √ ε ∥w ⋆ ∥2 2 1/L L k=1 σ k , Thus L k=1 σ k 1 - (L -1)ε ∥w ⋆ ∥2 2 2/L - L √ ε ∥w ⋆ ∥2 2 1/L ⩽ ∥W ⊤ 1 . . . W ⊤ L ∥ 2 .Using ( 31) and ( 32),1 - (L -1)ε ∥w ⋆ ∥2 2 2/L - L √ ε ∥w ⋆ ∥2 2 1/L ⩾ 1 - 1 8 - 1 8 = 3 4 .Moreover, the product of the singular values can be lower-bounded by the smallest one to the power L, so 3(minσ k ) L 4 ⩽ ∥W ⊤ 1 . . . W ⊤ L ∥ 2 .Let us show that, for t large enough,min σ k ⩽ (2∥w ⋆ ∥ 2 ) 1/L -ε .If it were not the case, we would haveEQUATION)where the second inequality holds by Lemma 4 sinceε (2∥w ⋆ ∥ 2 ) 1/L ⩽ 1 2, by ( 31), and the third since2Lε (2∥w ⋆ ∥ 2 ) 1/L ⩽ 1 4 ⇔ Lε ⩽ (2∥w ⋆ ∥ 2 ) 1/L 8 ,also by (31). Since 9/8 > 1, the inequality ( 34) is a contradiction with the fact that ∥W⊤ 1 . . . W ⊤ L ∥ 2 = ∥w prod ∥ 2 → ∥w ⋆ ∥ 2 , showing that, for t large enough, min σ k ⩽ (2∥w ⋆ ∥ 2 ) 1/L -ε . Thus min σ 2 k ⩽ (2∥w ⋆ ∥ 2 ) 2/L + ε 2 -2ε(2∥w ⋆ ∥ 2 ) 1/L . By Lemma 2, for all k ∈ {1, . . . , L}, σ 2 k ⩽ min σ 2 j + ε ⩽ (2∥w ⋆ ∥ 2 ) 2/L + ε 2 -2ε(2∥w ⋆ ∥ 2 ) 1/L + ε ⩽ (2∥w ⋆ ∥ 2 ) 2/L , since ε 2 -2ε(2∥w ⋆ ∥ 2 ) 1/L + ε ⩽ 0 ⇔ ε ⩽ 2(2∥w ⋆ ∥ 2 ) 1/L -1 , which holds true by Assumption (A 2 ) since ε ⩽ 1 ⩽ 2(2∥w ⋆ ∥ 2 ) 1/L -1 since ∥w ⋆ ∥ 2 ⩾ 1.The lower bound unfolds by definition of S min since W SI is a minimizer of the empirical risk. To obtain the upper bound, we proceed similarly to the proof of Theorem 2. For simplicity, we denote W = W SI in the remainder of the proof. We have, as in the proof of Theorem 2,EQUATIONwithEQUATION∇R 1 ( wprod ) = 2 n X ⊤ X L k=1 W ⊤ 1 . . . W ⊤ k-1 ( W ⊤ k -W ⊤ k )W ⊤ k+1 . . . W ⊤ L + O(ξ 2 ) .We recall that, as in the proof of Theorem 2, the notation O is taken with respect to the limit when ξ → 0, everything else being fixed. By subadditivity of the operator norm and by Corollary 1, we haveEQUATION)Moving on to bounding the Frobenius norm of ∆ k , we observe that the dominating term of this matrix decomposes as a rank-one matrix. Thus, again by subadditivity of the operator norm, equation ( 36) and Corollary 1,∥∆ k ∥ F = ∥W ⊤ k+1 . . . W ⊤ L ∥ 2 ∥∇R 1 ( wprod ) ⊤ W ⊤ 1 . . . W ⊤ k-1 ∥ 2 + O(ξ 2 ) ⩽ ∥W k+1 ∥ 2 . . . ∥W L ∥ 2 ∥∇R 1 ( wprod )∥ 2 ∥W 1 ∥ 2 . . . ∥W k-1 ∥ 2 + O(ξ 2 ) ⩽ 2Λ(2∥w ⋆ ∥ 2 ) 2(L-1)/L L k=1 ∥W k -Wk ∥ 2 + O(ξ 2 ) . Then ∥∆ k ∥ 2 F ⩽ 4Λ 2 (2∥w ⋆ ∥ 2 ) 4(L-1)/L L k=1 ∥W k -Wk ∥ 2 2 + O(ξ 3 ) ⩽ 4LΛ 2 (2∥w ⋆ ∥ 2 ) 4(L-1)/L L k=1 ∥W k -Wk ∥ 2 2 + O(ξ 3 ) ⩽ 4LΛ 2 (2∥w ⋆ ∥ 2 ) 4(L-1)/L L k=1 ∥W k -Wk ∥ 2 F + O(ξ 3 ) .Thus, by ( 35),S(W) 2 ⩽ lim ξ→0 4L 2 Λ 2 (2∥w ⋆ ∥ 2 ) 4(L-1)/L + O(ξ) ⩽ 2 6 L 2 Λ 2 ∥w ⋆ ∥ 4-4 L 2 . Therefore S(W) ⩽ 8LΛ∥w ⋆ ∥ 2-2 L 2, which concludes the proof by the second lower bound on S min of Theorem 2, where here ∥p∥ = 1 and a ⩾ λ > 0.To alleviate notations, we omit in this proof to write the explicit dependence of parameters on time. Starting from (5), since the gradient decomposes as a rank-one matrix, we have∥∇ k R L (W)∥ 2 F = ∥W ⊤ k+1 . . . W ⊤ L p∥ 2 2 ∥∇R 1 (w prod ) ⊤ W ⊤ 1 . . . W ⊤ k-1 ∥ 2 2 ⩾ σ 2 min (W ⊤ k+1 . . . W ⊤ L )∥p∥ 2 2 σ 2 min (W ⊤ 1 . . . W ⊤ k-1 )∥∇R 1 (w prod )∥ 2 2 By Lemma 8, ∥∇R 1 (w prod )∥ 2 2 ⩾ 4λ(R L (W) -R min ) . Recall the notation introduced in Section 5 Π L:k := W L . . . W k and Π k:1 := W k . . . W 1 . Then ∥∇ k R L (W)∥ 2 F ⩾ 4λ∥p∥ 2 2 σ 2 min (Π L:k+1 )σ 2 min (Π k-1:1 )(R L (W) -R min ) .(37) We now let u > 0 (whose value will be specified next), andt * = inf t ∈ R + , ∃k ∈ {1, . . . , L}, ∥θ k ∥ F > 1 64 exp(-2s 2 -4u) .Let us lower bound the minimum singular values of Π k: and Π :k uniformly for t ∈ [0, t ⋆ ] by Lemma 3. By definition of t ⋆ , for t ⩽ t ⋆ and for any k ∈ {1, . . . , L}, ∥θ k ∥ 2 ⩽ 1 64 exp(-2s 2 -4u). Therefore, renaming C1 to C4 the constants of Lemma 3, and taking u = C3 , we get that, ifEQUATIONthen, with probability at least1 -8 exp - d C2 3 32s 2 ,it holds for t ⩽ t ⋆ and for all k ∈ {1, . . . , L} thatEQUATIONandσ min (Π k:1 ) ⩾ 1 4 exp - 2s 2 d -C3 .By symmetry, the same statement holds for Π L:k instead of Π k:1 with the same probability. By the union bound, the eventE 1 := 1⩽k⩽L ∥Π k:1 ∥ 2 , ∥Π L:k ∥ 2 ⩽ 4 exp s 2 2 + C3 and σ min (Π k:1 ), σ min (Π L:k ) ⩾ 1 4 exp - 2s 2 d -C3holds with probability at least1 -16 exp - d C2 3 32s 2 . Let now E 2 := R L (W(0)) -R min ⩽ C 4 λ 2 ∥p∥ 2 2 Λand E 3 the event that the gradient flow dynamics converge to a global minimizer W RI of the risk satisfying the statement (8). We show next that, if E 1 and E 2 hold, then E 3 must hold, which shall conclude the proof of the Theorem withEQUATIONUnder E 1 , we haveσ 2 min (Π L:k+1 )σ 2 min (Π k-1:1 ) ⩾ 1 2 8 exp - 8s 2 d -4 C3 ,thus by ( 37)∥∇ k R L (W)∥ 2 F ⩾ 1 2 6 exp - 8s 2 d -4 C3 λ∥p∥ 2 2 (R L (W) -R min ) .Therefore we get the PL condition, fort ⩽ t ⋆ , L k=1 ∇ k R L (W) 2 F ⩾ µ(R L (W) -R min ) , with µ := 1 2 6 exp - 8s 2 d -4 C3 λ∥p∥ 2 L .By Lemma 7, this implies that, forEQUATIONLet us now show that t ⋆ = ∞. We have, since θ(0) = 0 and by definition of the gradient flow,EQUATION)We now upper bound the gradient as follows: starting again from (5),∥∇ k R L (W)∥ F = ∥W ⊤ k+1 . . . W ⊤ L p∥ 2 ∥∇R 1 (w prod ) ⊤ W ⊤ 1 . . . W ⊤ k-1 ∥ 2 ⩽ ∥W ⊤ k+1 . . . W ⊤ L ∥ 2 ∥p∥ 2 ∥W ⊤ 1 . . . W ⊤ k-1 ∥ 2 ∥∇R 1 (w prod )∥ 2 .By Lemma 8, we getEQUATIONBy ( 39) and ( 41), for t ⩽ t ⋆ ,EQUATIONPlugging this into (42), we get, for t ⩽ t ⋆ ,∥θ k ∥ F ⩽ 16 exp(s 2 + 2 C3 ) √ Λ∥p∥ 2 R L (W(0)) -R min L t 0 e -µ 2 τ dτ ⩽ 32 exp(s 2 + 2 C3 ) √ Λ∥p∥ 2 R L (W(0)) -R min L µ = 2 11 exp(s 2 + 8s 2 d + 6 C3 ) √ Λ R L (W(0)) -R min λ∥p∥ 2 ,where the last equality comes from the definition of µ. By E 2 and by definition (40) of C 4 , for t ⩽ t ⋆ , ∥θ k ∥ F ⩽ 2 11 exp(s 2 + 8s 2 d + 6 C3 ) C 3 ⩽ 1 128 exp(-2s 2 -4 C3 ) .If we had t ⋆ < ∞, we would have, by definition of t ⋆ , ∥θ k (t ⋆ )∥ F ⩾ 1 64 exp(-2s 2 -4 C3 ). This contradicts the equation above, showing that t ⋆ = ∞. By (41), this implies convergence of the risk to its minimum. We also see by ( 44) that ∇ k R L (W) is integrable, so W has a limit as t goes to infinity. This limit W RI is a minimizer of the risk and satisfies the condition (8) by definition of t ⋆ = ∞.Extension to multivariate regression. We emphasize that a very similar proof holds in the case of multivariate regression, where the neural network is defined as the linear map from R d to R dx → W L . . . W 1 x , and we now aim at minimizing the mean squared error lossR L (W) = 1 n ∥Y -W L . . . W 1 X ⊤ ∥ 2 2 ,where X, Y ∈ R n×d . In this case, as shown in Zou et al. (2020) ; Sander et al. (2022) , the following bounds on the gradient hold:∥∇ k R L (W)∥ 2 F ⩾ 4λσ 2 min (Π L:k+1 )σ 2 min (Π k-1:1 )(R(W) -R min ) , and ∥∇ k R L (W)∥ 2 F ⩽ 4Λ∥Π L:k+1 ∥ 2 2 ∥Π k-1:1 ∥ 2 2 (R(W) -R min ) .Comparing with (37) and ( 43), the only difference is the absence of ∥p∥ 2 here. From there, the same computations as above hold (taking ∥p∥ 2 = 1), and give the following result.Theorem 5. There exist C 1 , . . . , C 5 > 0 depending only on s such that, if L ⩾ C 1 and d ⩾ C 2 , then, with probability at least1 -16 exp(-C 3 d) , if R L (W(0)) -R min ⩽ C 4 λ 2 ∥p∥ 2 2 Λ ,the gradient flow dynamics (4) converge to a global minimizer W RI of the risk. Furthermore, the minimizer W RI satisfiesW RI k = I + s √ Ld N k + 1 L θ RI k with ∥θ RI k ∥ F ⩽ C 5 , 1 ⩽ k ⩽ L .We begin by proving the result when θ = 0, then explain how to extend to any θ in a ball around 0.Recall that Π k:1 denotes the product of weight matrices up to the k-th. LetM s,u = 2 exp s 2 2 + u , m s,u = 1 2 exp - 2s 2 d -u ,and denote by A the event that for all k, ∥Π k:1 ∥ 2 ⩽ M s,u and B the event that for all k, σ min (Π k:1 ) ⩾ m s,u . We begin by bounding the probability of Ā, then bound the probability of B ∩ A.Useful identities. To this aim, we first introduce some notations and derive some useful identities. We let F k the filtration generated by the random variables N 1 , . . . , N k . For some h 0 ∈ R d , we leth k = Π k:1 h 0 = I + s √ Ld N k . . . I + s √ Ld N 1 h 0In particular, h k+1 and h k are related throughh k+1 = I + s √ Ld N k+1 h k .Taking the squared norm, we get∥h k+1 ∥ 2 2 = ∥h k ∥ 2 2 + s 2 Ld ∥N k+1 k ∥ 2 2 + 2s √ Ld h ⊤ k N k+1 h k .Dividing by ∥h k ∥ 2 2 and taking the logarithm leads toEQUATION) Let Y k,1 = s 2 ∥N k+1 h k ∥ 2 2 Ld∥h k ∥ 2 2 , Y k,2 = 2sh ⊤ k N k+1 h k √ Ld∥h k ∥ 2 2 ,andEQUATIONThis shows in particular that Y k,1 and Y k,2 are independent of N 1 , . . . , N k , and depend only on N k+1 . Then, lettingS k,1 = k-1 j=0 Y j,1 and S k,2 = k-1 j=0 Y j,2 ,both are sums of i.i.d. random variables. In particular,S L,1 ∼ s 2 Ld χ 2 (Ld) and S L,2 ∼ N 0, 4s 2 d .Bound on P( Ā). We first bound the deviation of the norm of h k = Π k:1 h 0 for any fixed h 0 ∈ R d , then conclude on the operator norm of Π k:1 by an ε-net argument. For any (fixed) h 0 ∈ R d and u > 0, we haveP max 1⩽k⩽L ∥h k ∥ 2 2 ∥h 0 ∥ 2 2 ⩾ exp(s 2 + 2u) = P max 1⩽k⩽L ln(∥h k ∥ 2 2 ) -ln(∥h 0 ∥ 2 2 ) ⩾ s 2 + 2u = P max 1⩽k⩽L k-1 j=0 ln 1 + s 2 ∥N j+1 h j ∥ 2 2 Ld∥h j ∥ 2 2 + 2sh ⊤ j N j+1 h j √ Ld∥h j ∥ 2 2 ⩾ s 2 + 2u ⩽ P max 1⩽k⩽L k-1 j=0 s 2 ∥N j+1 h j ∥ 2 2 Ld∥h j ∥ 2 2 + 2sh ⊤ j N j+1 h j √ Ld∥h j ∥ 2 2 ⩾ s 2 + 2u , by using ln(1 + x) ⩽ x. Thus P max 1⩽k⩽L ∥h k ∥ 2 2 ∥h 0 ∥ 2 2 ⩾ exp(s 2 + 2u) ⩽ P max 1⩽k⩽L S k,1 + S k,2 ⩾ s 2 + 2u ⩽ P max 1⩽k⩽L S k,1 ⩾ s 2 + u + P max 1⩽k⩽L S k,2 ⩾ uby the union bound. We now study the two deviation probabilities separately. Beginning by the first one, recall that Ld s 2 S L,1 follows a χ 2 (Ld) distribution. Chi-squared random variables are subexponential, and more precisely satisfy the following property (Ghosh, 2021) : if X ∼ χ 2 (c) and u > 0, thenEQUATIONSince the S k,1 are increasing, we have, for u > 0,P max 1⩽k⩽L S k,1 ⩾ s 2 + u = P(S L,1 ⩾ s 2 + u) = P Ld s 2 S L,1 ⩾ Ld + Ldu s 2 ⩽ exp - L 2 d 2 u 2 4s 4 (Ld + Ldu s 2 ) = exp - Ldu 2 4s 4 + us 2 .Moving on to S k,2 , we have, for u, λ > 0,P max 1⩽k⩽L S k,2 ⩾ u = P max 1⩽k⩽L exp(λS k,2 ) ⩾ exp(λu) .Furthermore, exp(λS k,2 ) is a sub-martingale, since S k,2 is a martingale and by Jensen's inequality:E(exp(λS k+1,2 )|F k ) ⩾ exp(λE(S k+1,2 |F k )) = exp(λS k,2 ) .Thus, by Doob's martingale inequality,P max 1⩽k⩽L S k,2 ⩾ u ⩽ E(exp(λS L,2 )) exp(-λu) . Furthermore, since S L,2 ∼ N (0, 4s 2 d ), E(exp(λS L,2 )) = exp 4s 2 λ 2 d .Therefore,P max 1⩽k⩽L S k,2 ⩾ u ⩽ exp -λu + 4s 2 λ 2 d .This quantity is minimal for λ = du 8s 2 . We get, for all u > 0,EQUATION)Therefore, for any u > 0,EQUATION)To conclude on the operator norm of Π k:1 , consider Σ a 1/2-net of the unit sphere of R d . By Vershynin (2018, Corollary 4.2.13) , it is possible to take such a net of cardinality 5 d . Let us show that, for any u ∈ R,EQUATIONIndeed, assume that there exists k such that ∥Π k:1 2 ⩾ 2 exp( s 2 2 + u). By definition of the operator norm, there exists x on the unit sphere of R d such that∥Π k:1 x∥ = ∥Π k:1 ∥ 2 .But then, by definition of Σ, there existsh 0 ∈ Σ such that ∥x -h 0 ∥ 2 ⩽ 1/2. Then ∥Π k:1 (x -h 0 )∥ 2 ⩽ 1 2 ∥Π k:1 ∥ 2 .Therefore, by the triangular inequality,∥Π k:1 h 0 ∥ 2 ⩾ ∥Π k:1 x∥ 2 -∥Π k:1 (x -h 0 )∥ 2 ⩾ 1 2 ∥Π k:1 ∥ 2 .Then∥Π k:1 h 0 ∥ 2 ⩾ 1 2 ∥Π k:1 ∥ 2 ⩾ exp s 2 2 + u ,which proves (50). By the union bound, we conclude thatP( Ā) = P max 1⩽k⩽L ∥Π k:1 ∥ 2 ⩾ 2 exp s 2 2 + u ⩽ P h0∈Σ max 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩾ exp s 2 2 + u ⩽ |Σ|P max 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩾ exp s 2 2 + u ,where now h 0 denotes any unit-norm fixed vector in R d . Thus, by (49),P( Ā) ⩽ 5 d P max 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 2 ⩾ exp(s 2 + 2u) ⩽ 5 d exp - du 2 16s 2 + exp - Ldu 2 4s 4 + us 2 .This upper bound will be simplified in the conclusion of the proof.Bound on P( B ∩ A). We now move on to proving the lower-bound on σ min (Π k:1 ). We again use an ε-net argument, as follows. Let Σ be an ε-net for ε = ms,u Ms,u . ThenP A ∩ min 1⩽k⩽L σ min (Π k:1 ) ⩽ m s,u ⩽ P h0∈Σ min 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩽ 2m s,u .Indeed, assume that A holds, that is ∥Π k:1 ∥ 2 ⩽ M s,u , and that there exists k such that σ min (Π k:1 ) ⩽ m s,u .By definition of the singular value, there exists x on the unit sphere of R d such that∥Π k:1 x∥ 2 = σ min (Π k:1 ) .But then, by definition of Σ, there existsh 0 ∈ Σ such that ∥x -h 0 ∥ 2 ⩽ ε. Then ∥Π k:1 (x -h 0 )∥ 2 ⩽ ε∥Π k:1 ∥ 2 .Under A, the right-hand side is smaller than εM s,u . Therefore, by the triangular inequality,∥Π k:1 h 0 ∥ 2 ⩽ ∥Π k:1 x∥ 2 + ∥Π k:1 (h 0 -x)∥ 2 ⩽ m s,u + εM s,u = 2m s,u .By the union bound, we conclude thatP A ∩ min 1⩽k⩽L σ min (Π k:1 ) ⩽ m s,u ⩽ P A ∩ h0∈Σ min 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩽ 2m s,u ⩽ P h0∈Σ min 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩽ 2m s,u ⩽ |Σ|P min 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩽ 2m s,u = 2M s,u m s,u + 1 d P min 1⩽k⩽L ∥Π k:1 h 0 ∥ 2 ⩽ 2m s,u ,where the cardinality of the net is given by Vershynin Corollary 4.2.13). We now take a fixed h 0 ∈ R d , and computeP min 1⩽k⩽L ∥h k ∥ 2 ∥h 0 ∥ 2 ⩽ 2m s,u = P min 1⩽k⩽L ∥h k ∥ 2 2 ∥h 0 ∥ 2 2 ⩽ 4m 2 s,u .Denote by E the event1⩽k⩽L Y k,2 ⩾ - 1 2 .We have, by (45),P min 1⩽k⩽L ∥h k ∥ 2 2 ∥h 0 ∥ 2 2 ⩽ 4m 2 s,u = P min 1⩽k⩽L ln(∥h k ∥ 2 2 ) -ln(∥h 0 ∥ 2 2 ) ⩽ ln(4m 2 s,u ) ⩽ P min 1⩽k⩽L k-1 j=0 ln(1 + Y j,2 ) ⩽ ln(4m 2 s,u ) ⩽ P min 1⩽k⩽L k-1 j=0 ln(1 + Y j,2 ) ⩽ ln(4m 2 s,u ) ∩ E + P( Ē) .Using the inequality ln(1 + x) ⩾ x -x 2 for x ⩾ -1 /2, we obtainP ∥h L ∥ 2 2 ∥h 0 ∥ 2 2 ⩽ 4m 2 s,u ⩽ P min 1⩽k⩽L k-1 j=0 (Y j,2 -Y 2 j,2 ) ⩽ ln(4m 2 s,u ) ∩ E + P( Ē) .Thus, by the union bound,P min 1⩽k⩽L ∥h L ∥ 2 2 ∥h 0 ∥ 2 2 ⩽ 4m 2 s,u ⩽ P min 1⩽k⩽L k-1 j=0 (Y j,2 -Y 2 j,2 ) ⩽ ln(4m 2 s,u ) + L-1 k=0 P Y k,2 < - 1 2 =: P 1 + P 2 .We handle both terms separately. Beginning by the second term, we have, for k ∈ {1, . . . L},P Y k,2 < - 1 2 ⩽ exp - Ld 32s 2 ,where we used (46) and the tail bound P(N ⩾ u) ⩽ e -u 2 /2 if N ∼ N (0, 1). Moving on to the first term, we haveP 1 = P min 1⩽k⩽L S k,2 -S k,3 ⩽ ln(4m 2 s,u ) ,where we letS k,3 = k-1 j=0 Y 2 j,2. By definition of m s,u , we haveln(4m 2 s,u ) = - 4s 2 d -2u .Thus we can split the probability into two parts by the union bound:P 1 ⩽ P min 1⩽k⩽L S k,2 ⩽ -u + P max 1⩽k⩽L S k,3 ⩾ 4s 2 d + u .Let us bound each probability separately. We have, for u > 0,P min 1⩽k⩽L S k,2 ⩽ -u = P max 1⩽k⩽L -S k,2 ⩾ u = P max 1⩽k⩽L S k,2 ⩾ u ⩽ exp - du 2 16s 2 ,by symmetry of S k,2 and by (48). Moving on to S k,3 , we have thatY 2 j,2 |F j ∼ 4s 2 Ld N (0, 1) 2 = 4s 2 Ld χ 2 (1) , thus Ld 4s 2 S k,3 ∼ χ 2 (L). By of S k,3 and by (47),P max 1⩽k⩽L S k,3 ⩾ 4s 2 d + u = P S L,3 ⩾ 4s 2 d + u = P Ld 4s 2 S L,3 ⩾ L + Ldu 4s 2 ⩽ exp - L 2 d 2 u 2 64s 4 (L + Ldu 4s 2 ) = exp - Ld 2 u 2 16(4s 2 + du).Putting everything together, we proved thatP( B ∩ A) ⩽ 2M s,u m s,u + 1 d (P 1 + P 2 ) = exp d ln 2M s,u m s,u + 1 (P 1 + P 2 ) ,withP 1 ⩽ exp - du 2 16s 2 + exp - Ld 2 u 2 16(4s 2 + du) and P 2 ⩽ L exp - Ld 32s 2 .Bounding the probability of failure. Putting together the two main bounds we showed, we haveP( Ā ∪ B) = P( Ā) + P( B ∪ A) ⩽ 5 d exp - du 2 16s 2 + exp - Ldu 2 4s 4 + us 2 + exp d ln 2M s,u m s,u + 1 exp - du 2 16s 2 + exp - Ld 2 u 2 16(4s 2 + du) + L exp - Ld 32s 2. This expression is valid for all u, s > 0 and L, d ⩾ 1. We now simplify the expression of our upper bound by algebraic computations using the assumptions of the Theorem. To somewhat alleviate the technicality of the computations, we stop at this point tracking some of the explicit constants and let C denote a positive absolute constant that might vary from equality to equality. We show next that the conditions L ⩾ C, d ⩾ C, u ⩾ C max(s 2 , s -2/3 ), u ⩽ CL 1/4 (51) imply thatP( Ā ∪ B) ⩽ 8 exp - du 2 32s 2 .This shall conclude the proof of the Lemma withC 1 = C, C 2 = C, C 3 = C max(s 2 , s -2/3 ), C 4 = C .First note that the conditions (51) imply that u ⩾ Cs and u ⩾ C .We study the terms of the bound on P( Ā ∪ B) one by one. First, we have by ( 52) that5 d exp - du 2 16s 2 = exp d ln(5) - du 2 16s 2 ⩽ exp - du 2 32s 2 . Next, 5 d exp - Ldu 2 4s 4 + us 2 ⩽ 5 d exp - Ldu 2 8us 2 = 5 d exp - Ldu 8s 2 = exp d ln(5) - Ldu 8s 2 ⩽ exp - Ldu 16s 2 ,where the first inequality uses u ⩾ s 2 by (51), and the last one uses that 51). Thus we can bound the three remaining terms appearing in the bound of P( Ā ∪ B), as follows:exp d ln 2M s,u m s,u + 1 exp - du 2 16s 2 ⩽ exp 3d + 3du - du 2 16s 2 ⩽ exp - du 2 32s 2 , since u 2 ⩾ Cs 2 (1 + u). This is the case by ( 52), which implies thatu 2 = 1 2 u 2 + 1 2 uu ⩾ C 2 2 s 2 + C 2 s 2 u = C 2 + C 2 s 2 (1 + u) .Next, by ( 53 Let us show that this implies thatEQUATION)This statement is true because the three terms appearing inside the exponential of the right-hand side of ( 54) can be bounded by C √ Ldu 2 s 2. More precisely, we haveLd 32s 2 ⩾ C √ Ldu 2 s 2 ⇔ √ L ⩾ Cu 2 ,which holds by (51);3du ⩽ C √ Ldu 2 s 2 ⇔ Cs 2 ⩽ √ Lu ,which is implied by L ⩾ C and u ⩾ Cs 2 by (51); and3d ⩽ C √ Ldu 2 s 2 ⇔ Cs 2 ⩽ √ Lu 2 ,which is implied by L ⩾ C and u ⩾ Cs by ( 51) and ( 52). Finally, we use Lemma 5 to bound the right-hand side of ( 55) by4 exp - du 2 32s 2 .This is possible for L ⩾ C and du 2 32s 2 ⩾ 1, which is implied by d ⩾ C and u ⩾ Cs, by ( 51) and ( 52). Collecting everything, we bounded P( Ā ∪ B) by8 exp - du 2 32s 2 ,which concludes the proof when θ = 0.Summary when θ = 0. In summary, we proved so far the following result: there exist C 1 , . . . , C 4 > 0 depending only on s such that, ifL ⩾ C 1 , d ⩾ C 2 , u ∈ [C 3 , C 4 L 1/4 ] ,then, with probability at least Comparison with the bounds of Zhang et al. (2022) . Theorem 1 of Zhang et al. (2022) gives an upper bound on the singular values of residual networks at initialization, and their Theorem 2 gives a lower bound on the norm of the activations. Comparing with our results, we note two important differences. First, their probability of failure grows quadratically with the depth, whereas ours is independent of depth. This is achieved by a more precise martingale argument making use of Doob's martingale inequality. Second, their lower bound incorrectly assumes that χ 2 random variables are sub-Gaussian (see equation ( 21) of their paper), while in fact they are only sub-exponential (Ghosh, 2021) . Finally, their upper bound holds for the productI + s √ Ld N k + 1 L θ k . . . I + s √ Ld N j + 1 L θ jfor any 1 ⩽ j ⩽ k ⩽ L, which could seem stronger than our result stated for j = 1. In fact, both statements are equivalent, because it is possible to deduce the statement for any j by combining the upper bound and the lower bound for j = 1. The precise argument is given in the beginning of the proof of our Lemma 9.The beginning of the proof is very similar to the one of Corollary 2. Denoting W := W RI , we haveEQUATIONwithEQUATION∇R 1 ( wprod ) = 2 n X ⊤ X L k=1 W ⊤ 1 . . . W ⊤ k-1 ( W ⊤ k -W ⊤ k )W ⊤ k+1 . . . W ⊤ L p + O(ξ 2 ) .At this point, the proofs diverge. We have, by subadditivity of the operator norm,∥∇R 1 ( wprod )∥ 2 ⩽ 2Λ L k=1 ∥W k-1 . . . W 1 ∥ 2 ∥ Wk -W k ∥ 2 ∥W L . . . W k+1 ∥ 2 ∥p∥ 2 + O(ξ 2 ) .Let us now briefly recall the outline of the proof of Theorem 4, which will be useful in bounding the quantity above. The proof shows the existence of C3 depending only on s such that, with high probability (which is exactly the probability in the statement of the Theorem), we have for all t ⩾ 0 and k ∈ {1, . . . , L} that∥W k-1 (t) . . . W 1 (t)∥ 2 ⩽ 4 exp s 2 2 + C3 and ∥W L (t) . . . W k+1 (t)∥ 2 ⩽ 4 exp s 2 2 + C3 ,as well asEQUATIONUnder this high-probability event, the proof of Theorem 4 shows convergence of the gradient flow to W = W RI . In particular, this means that W also verifies the bounds on the operator norm of the matrix products. We therefore obtain∥∇R 1 ( wprod )∥ 2 ⩽ 2Λ L k=1 4 exp s 2 2 + C3 ∥ Wk -W k ∥ 2 4 exp s 2 2 + C3 ∥p∥ 2 + O(ξ 2 ) = 32 exp(s 2 + 2 C3 )Λ∥p∥ 2 L k=1 ∥ Wk -W k ∥ 2 + O(ξ 2 ) .Moving on to bounding the Frobenius norm of ∆ k , we have∥∆ k ∥ F = ∥W ⊤ k+1 . . . W ⊤ L p∥ 2 ∥∇R 1 ( wprod ) ⊤ W ⊤ 1 . . . W ⊤ k-1 ∥ 2 + O(ξ 2 ) ⩽ ∥W L . . . W k+1 ∥ 2 ∥p∥ 2 ∥∇R 1 ( wprod )∥ 2 ∥W k-1 . . . W 1 ∥ 2 + O(ξ 2 ) ⩽ 2 9 exp(2s 2 + 4 C3 )Λ∥p∥ 2 2 L k=1 ∥W k -Wk ∥ 2 + O(ξ 2 ) ,by bounding the three norms by the expressions given above. Then∥∆ k ∥ 2 F ⩽ 2 18 exp(4s 2 + 8 C3 )Λ 2 ∥p∥ 4 2 L k=1 ∥W k -Wk ∥ 2 2 + O(ξ 3 ) ⩽ 2 18 exp(4s 2 + 8 C3 )LΛ 2 ∥p∥ 4 2 L k=1 ∥W k -Wk ∥ 2 2 + O(ξ 3 ) . ⩽ 2 18 exp(4s 2 + 8 C3 )LΛ 2 ∥p∥ 4 2 L k=1 ∥W k -Wk ∥ 2 F + O(ξ 3 ) .Thus, by (56),S(W) 2 ⩽ lim ξ→0 2 18 exp(4s 2 + 8 C3 )L 2 Λ 2 ∥p∥ 4 2 + O(ξ) . Therefore S(W) ⩽ 2 9 exp(2s 2 + 4 C3 )LΛ∥p∥ 2 2 .(58) To conclude, we need to upper bound ∥p∥ 2 by a constant times ∥w ⋆ ∥ 2 . To this aim, we leverage the bound from the assumptions of Theorem 4 on the risk at initialization, to show that ∥p∥ 2 cannot be too far away from ∥w ⋆ ∥ 2 . More precisely, by Lemma 8, since the covariance matrix X ⊤ X is full rank,R L (W(0)) -R min = 1 n ∥X(w prod (0) -w ⋆ )∥ 2 2 = 1 n (w prod (0) -w ⋆ ) ⊤ X ⊤ X(w prod (0) -w ⋆ ) ⩾ λ∥w prod (0) -w ⋆ ∥ 2 2 . Thus R L (W(0)) -R min ⩾ √ λ∥w prod (0) -w ⋆ ∥ 2 .Then, by the triangular inequality,∥w ⋆ ∥ 2 ⩾ ∥w prod (0)∥ 2 -∥w prod (0) -w ⋆ ∥ 2 ⩾ ∥W ⊤ 1 (0) . . . W ⊤ L (0)p∥ 2 - R L (W(0)) -R min λ ⩾ σ min (W L (0) . . . W 1 (0))∥p∥ 2 - R L (W(0)) -R min λ ⩾ 1 4 exp - 2s 2 d -C3 ∥p∥ 2 -C 3 λ Λ ∥p∥ 2 ,by (57) and by the assumption of Theorem 4 on the risk at initialization. We now note that the value of C 3 is given by (40) asC 3 = 2 -36 exp -4s 2 - 16s 2 C2 -20 C3 ,where C2 ⩽ d by (38). ThusC 3 λ Λ ⩽ C 3 = 2 -18 exp(-2s 2 - 8s 2 C2 -10 C3 ) < 1 4 exp - 2s 2 d -C3 . Denoting C ′ = 1 4 exp - 2s 2 d -C3 -C 3 ∈ (0, 1) ,we therefore obtain that ∥w ⋆ ∥ 2 ⩾ C ′ ∥p∥ 2 . Therefore, by (58),S(W) ⩽ 2 9 exp(2s 2 + 4 C3 )(C ′ ) -2+ 2 L LΛ∥w ⋆ ∥ 2-2 L 2 ∥p∥ 2 L2 . Thus, by the second lower bound on S min from Theorem 2, and since a ⩾ λ > 0,S(W) S min ⩽ 2 8 exp(2s 2 + 4 C3 )(C ′ ) -2+ 2 L Λ λ ,which concludes the proof by settingC := 2 8 exp(2s 2 + 4 C3 )(C ′ ) -2 ⩾ 2 8 exp(2s 2 + 4 C3 )(C ′ ) -2+ 2 L .Our code is available at https://github.com/PierreMarion23/implicit-reg-sharpness.Our framework for experiments is JAX (Bradbury et al., 2018) . The experiments take around 3 hours to run on a laptop CPU.Setup. We take n = 50, d = 5, L = 10. The design matrix X is sampled from an isotropic Gaussian distribution. The target y is computed in two steps. First, we compute y 0 = Xw true + ζ, where w true and ζ are standard Gaussian vectors. Then, we compute w ⋆ 0 as the optimal regressor of y 0 on X. Finally, we let y = y 0 /∥w ⋆ 0 ∥ and w ⋆ = w ⋆ 0 /∥w ⋆ 0 ∥. This simplifies the expressions of our bounds by having w ⋆ of unit norm. All Gaussian random variables are independent.Details of Figure 1 . We consider a Gaussian initialization of the weight matrices, where the scale of the initialization (x-axis of some the graphs) is the standard deviation of the entries. All weight matrices are d × d, except the last one which is 1 × d. The square distance to the optimal regressor corresponds to ∥w prod -w ⋆ ∥ 2 2 . The largest eigenvalue of the Hessian is computed by a power iteration method, stopped after 20 iterations. In Figures 1a and 1b , the 95% confidence intervals are plotted. The number of gradient steps and number of independent repetitions depend on the learning rate, and are given below. For large values of the initialization scale, it may happen that the gradient descent diverges. Figure 3 shows the probability of divergence depending on the initialization scale and the learning rate. When the probability of divergence is equal to one, no point is reported in Figure 1 . When it is strictly between 0 and 1, the confidence intervals are computed over non-diverging runs.Figures 1c and 1d show one randomly-chosen run each. The plots are subsampled 5 times for readability, due to the oscillations in Figure 1d .Residual initialization. We now consider the case of a residual initialization as in Section 5. Results are given in Figure 4 . The scale of the initialization now corresponds to the hyperparameter s in (2). The projection vector p ∈ R d is a random isotropic vector of unit norm, which does not change during training. For each learning rate, we use 4, 000 steps of gradient descent, and perform 20 independent repetitions. The plots are similar to the case of Gaussian initialization, apart from the fact that the sharpness at initialization is better conditioned.As previously, for large values of the initialization scale, it may happen that the gradient descent diverges. Figure 5 shows the probability of divergence depending on the initialization scale and the learning rate.Details of Figure 2 . The setup is the same as for the residual initialization. For each learning rate, we use 1, 000 steps of gradient descent, and perform 50 independent repetitions. We take s = 0.25. Underdetermined regression and link to generalization. Although this is not our original motivation, we note that a simple change to our setting allows to make appear the connection between sharpness and generalization. To this aim, we consider the underdetermined case, where the number of data is lower than the dimension (while keeping the rest of the setup identical). Figure 6 shows in this case a correlation between generalization and sharpness. This suggests that the tools developed in the paper could be used in this case to understand the generalization performance of deep (linear) networks, and we leave this analysis for future work. We also qualitatively observe in Figure 6 a similar connection between learning rate, initialization scale and sharpness as in the case of full-rank data (Figure 1b ). We take here n = 15, d = 20, L = 5. The number of gradient steps and number of independent repetitions depend on the learning rate, and are given below. Other technical details are as Figure 1 .Learning (a) Sharpness at initialization and after training, for various learning rates and initialization scales. We qualitatively observe a similar connection between learning rate, initialization scale and sharpness as in the case of a full rank data matrix (Figure 1b of the paper). Non-linear MLPs. As a first attempt to extend our results to non-linear networks, we consider the case of non-linear MLPs. The non-linearity is GELU (Hendrycks and Gimpel, 2016) , a smooth version of ReLU, which we chose because smoothness is necessary in order to compute the sharpness. We qualitatively observe in Figure 7 a similar connection between learning rate, initialization scale and sharpness as for deep linear neural networks (Figure 1b ). For large initialization, the sharpness after training plateaus at 2/η, as in the linear case. For small initialization, the sharpness after training is less that 2/η, and is close to the bounds in the linear case (dotted black lines). For this experiment, we consider noiseless data, meaning that y 0 = Xw true (see paragraph "Setup" above for notations). We perform 20 independent repetitions of each experiment. The number of gradient steps depends on the learning rate, and is given below. Other details are as for Figure 1 . Finally, we also performed the same experiment in the case of noisy data y 0 = Xw true + ζ (no plot reported). We observed that the sharpness of the network reaches 2/η, for every learning rate and initialization scale reported in Figure 7 . We suspect that this is because the network (over)fits the noise in the data, resulting in a high sharpness.Learning rate Number of steps 0.005 160, 000 0.02 40, 000 0.07 16, 000Comments on the connection between gradient flow and gradient descent. In this paper, we show that gradient flow from a small-scale initialization is driven towards low-sharpness regions. This should imply that gradient flow and gradient descent up to a reasonably large learning rate should follow the same trajectory when starting from small-scale initialization, because they do not go in regions of high sharpness where the difference between gradient flow and gradient descent would become significant. This intuition is supported by Figure 1b where we see that, for small initializations, the sharpness after training is independent of the learning rate. We leave further investigation of these questions for future work.Progressive sharpening. Cohen et al. (2021) show that the edge of stability phase is typically preceded by a phase of progressive sharpening, where the sharpness steadily increases until reaching the value of 2/η. Our setting of small-scale initialization presents an example of such a progressive sharpening (although we make no statement on the monotonicity of the increase in sharpness). Other works have proposed analyses of progressive sharpening. Wang et al. (2022) suggest that progressive sharpening is driven by the increase in norm of the output layer. MacDonald et al. (2023) assume from empirical evidence a link between sharpness and the magnitude of the input-output Jacobian, and show that the latter has to be large for the loss to decrease. Finally, Agarwala et al. (2023) propose and analyze a simplified model with quadratic dependence in its parameters, which exhibits a progressive sharpening phenomenon.Connection with deep matrix factorization. Regression with deep linear networks can be seen as an instance of a matrix factorization problem. There is a well-established literature studying the implicit regularization of gradient descent for this class of problem (see, e.g., Gunasekar et al., 2017; Arora et al., 2019b; Li et al., 2021; Yun et al., 2021) . However, this line of work study under-determined settings where there are an infinite number of factorizations reaching zero empirical risk, and study the implicit regularization in function space. On the contrary, we consider an overdetermined setting where there is a single optimal regressor, and study the regularization in parameter space.• The method for calculating the error bars should explained (closed form formula, call to a library function, bootstrap, • 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. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: This is a foundational research paper, not tied to particular applications, let alone deployments. 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.We further note an incorrect use of concentration inequalities for χdistributions that are assumed to be sub-Gaussian while they are in fact only sub-exponential. Details are given in Appendix B.8. 


Decision: The paper proposes an implicit bias of gradient flow for deep linear networks and overdetermined univariate regression that the sharpness of the minimizer is bounded.

The reviewers discussed the paper in detail, most of the concerns (regarding link with generalization, extension to general settings, and tightness of the bounds, etc.) were resolved during the discussion period and thus the reviewers agreed upon positive ratings. Overall I recommend accept.

Reviews: review: paper_summary: paper_summary: The authors show three new results for deep linear networks (DLN).  First, they show that any DLN that implements the optimal linear regressor must have a certain "sharpness".  This amounts to a lower bound on the largest eigenvalue of the Hessian matrix at that set of weights.  They then argue that the weights found by using gradient flow with two different types of initialization find minima that have a sharpness that is within a constant of this lower bound.  They interpret these results as an implicit regularization of gradient flow towards flat minima.
strength_weakness: strength_weakness: The paper kind of has a mixed message, and doesn't really make the connection to generalization power in a way that is intepretable.  Theorem 1 shows that the sharpness has to grow at a rate that is essentially linear in the number of layers.  But then they mention that prior work indicates that flat minima generalize better, which suggests that making the DLN deeper is going to hurt you in this regard.  The subjective interpretation of gradient flow going to flat minima is hard to digest when the authors just proved that there are no flat minima.
questions: questions: Do these results add to our understanding of the generalization power of deep linear networks trained with gradient flow?
limitations: limitations: none
paper_summary: paper_summary: The paper considers a toy non-convex optimization problem, namely overdetermined univariate regression with a deep linear network. Their main contributions are:

1) A lower bound on the sharpness of the minimizers of the empirical risk. In particular, if the step size is chosen too big, then gradient descent will fail to converge.

2) They show that the sharpness of the minimizer found by gradient flow is at most a constant multiple of the above lower bound, for both small-scale and residual initialization. This constant does not depend on the width or depth, but only on the conditioning number of the data. This shows an implicit regularization towards flat minima (the empirical minimizers of the risk can have arbitrarily large sharpness).

On the technical side, they prove convergence of the GF and characterize the solution at convergence in both initialization regimes. Numerical simulations are provided to illustrate and substantiate their claims.
strength_weakness: strength_weakness: The paper considers a simple setting: deep linear network and underdetermined regime. This allows the authors to characterize interesting and non-trivial behavior. However, it is unclear how much these results can extend beyond this simple setting.
questions: questions: Here, it is assumed that the data matrix is full rank. Is it not sufficient to assume that there exists a unique minimizer?
limitations: limitations: Yes
paper_summary: paper_summary: This paper studies the implicit bias of gradient flow on deep linear networks for overdetermined univariate regression problems. A lower bound on the sharpness of any minimizer is first derived for the general data covariance matrix, then it is shown that gradient flow with small initialization finds a minimizer with a sharpness no larger than a constant times the lower bound, where the constant depends on the condition number of the data covariance matrix. A similar result is also derived for GF on linear residual networks.
strength_weakness: strength_weakness: 1. **Tightness of Theorem 2**: In Mulayoff and Michaeli (2020), the lower bound on the sharpness is tight for whitened data: there exists a minimizer that achieves the lower bound. However, Theorem 2 in this paper only provides a lower bound. Is this lower bound improvable? For example, if the true lower bound has $\Lambda$ instead of $\lambda$, then the following results show that the sharpness minimizes found by GF is no larger than a constant times the lower bound, where the constant now does not depend on condition number.
2. **Relevance of the problem setting**:  The sharpness is often considered to be affecting the generalization error of the trained network. However, this paper studies overdetermined linear regression, where every minimizer corresponds to the same input-output map $w^*$, thus having the same generalization error. Why does one care about studying the sharpness of the minimizer in this case, if sharpness does not affect generalization at all?
questions: questions: See Weaknesses
limitations: limitations: See Weaknesses
paper_summary: paper_summary: The paper explores the behavior of deep linear neural networks in the context of overdetermined univariate regression, which is the setting of a univariate response $y$, more samples than input dimensions, and nonsingular data covariance matrix. 

The first result concerns the empirical risk minimizer (this is not equal to the gradient descent or gradient flow estimator)  and shows a lower bound on the sharpness of the ERM.

The next set of results show that gradient flow, a limit of gradient descent with a vanishing learning rate, implicitly regularizes the network towards flat minima, with sharpness close to a constant times the lower bound. This is proven for both small-scale and residual initializations.
strength_weakness: strength_weakness: * I missed a rigorous theoretical connection between the three main results and the behavior of the gradient descent minimizer. If the connection is only empirical in nature, it might be good to highlight this
* The work employs, to my eyes at least, a rather limited definition of sharpness. Are the results really very particular to sharpness being the largest eigenvalue of the Hessian?
* I'm missing some motivation for the two initialization schemes studied. Are they studied simply because they've been studied before and so the results here can borrow from existing work? 
* I find the title misleading. An architecture cannot have a preference for flat minima(?). It must be the architecture in conjunction with how it is learned. In this case, the main results concern very specific configurations of gradient flow. 
* The implication for learning rate design from this theoretical analysis seems rather weak to me. I mean that I cannot envision a practitioner implementing a learning rate design based on the results presented here.
questions: questions: * Apologies for this very naive question, but what does $S(\mathcal W)$, the largest eigenvalue of $\mathcal W=(W_1,\ldots,W_L)$ mean exactly? This $\mathcal W$ is not itself a matrix. 
* Each of your results in Sections 3, 4, and 5 are based on existing proof techniques. Specifically, Section 3 proof is based on Mulayoff and Michaeli 2020, Section 4 on Ji and Telgarsky (2020), and Section 5 on Zou et al 2020, Sander et al 2022. I am not familiar with these papers. Do they attempt to prove similar things to this paper? Do you reach the same conclusions but with different proof techniques?
limitations: limitations: * There are some obvious limitations that are par for the course in papers attempting to prove theoretical results in deep learning. Namely, the architecture here is extremely oversimplified. The data covariance being full rank is also unrealistic. Finally, the main results are about gradient flow and discrete-time stochastic gradient descent is surely very different in nature.

Comments: title: Thanks!
comment: Thank you again for your time and very interesting review.
Thank you for your careful response both here an in the general rebuttal.
Thank you again for your time and very interesting review.
Thank you to the authors for their detailed response. I'd like to maintain my original rating that the paper is a "technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations."
Thank you again for your time and very interesting review.
Thank you for your very thoughtful review and for raising the score. We will add the discussion to the next version of the paper.
I thank the authors for their detailed response. I have no further questions at the moment.
Thank the authors for addressing my concerns. I increase the score from 4 to 6. I suggest the authors adding these discussions to the appendix.