Michael Arbel

Romain Menegaux

Pierre Wolinski

Rethinking Gauss-Newton for learning over-parameterized models

This work studies the global convergence and implicit bias of Gauss Newton's (GN) when optimizing over-parameterized one-hidden layer networks in the mean-field regime. We first establish a global convergence result for GN in the continuous-time limit exhibiting a faster convergence rate compared to GD due to improved conditioning. We then perform an empirical study on a synthetic regression task to investigate the implicit bias of GN's method. While GN is consistently faster than GD in finding a global optimum, the learned model generalizes well on test data when starting from random initial weights with a small variance and using a small step size to slow down convergence. Specifically, our study shows that such a setting results in a hidden learning phenomenon, where the dynamics are able to recover features with good generalization properties despite the model having sub-optimal training and test performances due to an under-optimized linear layer. This study exhibits a trade-off between the convergence speed of GN and the generalization ability of the learned solution.

Introduction

Gauss-Newton (GN) and related methods, like the Natural Gradient (NG) can offer faster convergence rates compared to gradient descent due to the improved dependence on the conditioning of the problem [START_REF] Botev | Practical gauss-newton optimisation for deep learning[END_REF][START_REF] Grosse | A Kronecker-factored Approximate Fisher Matrix for Convolution Layers[END_REF][START_REF] Martens | New insights and perspectives on the natural gradient method[END_REF][START_REF] Martens | Optimizing Neural Networks with Kronecker-factored Approximate Curvature[END_REF]. For these reasons, these methods have attracted attention in machine learning as an alternative to gradient descent when optimizing ill-conditioned problems arising from the use of vastly over-parameterized networks and large training sets [START_REF] Martens | Deep Learning via Hessian-free Optimization[END_REF][START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF]. Unfortunately, GN's high computational cost per iteration, which involves solving an expensive linear system, restricts its applicability in large-scale deep learning optimization. Addressing this challenge has been a primary focus, with extensive efforts dedicated to crafting computationally efficient approximations for GN/NG methods where the aim is to strike a delicate balance between computational cost and optimization speed, thereby enhancing scalability [START_REF] Arbel | Kernelized Wasserstein Natural Gradient[END_REF][START_REF] Botev | Practical gauss-newton optimisation for deep learning[END_REF][START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF][START_REF] Grosse | A Kronecker-factored Approximate Fisher Matrix for Convolution Layers[END_REF][START_REF] Martens | Optimizing Neural Networks with Kronecker-factored Approximate Curvature[END_REF][START_REF] Ren | Efficient subsampled gauss-newton and natural gradient methods for training neural networks[END_REF].

While substantial effort has been put into making GN scalable, little is known about the global convergence of these methods and their generalization properties when optimizing neural networks, notably compared to gradient descent. Understanding these properties is particularly relevant knowing that state-of-the-art machine learning methods rely on over-parameterized networks for which infinitely many solutions interpolating the data exist albeit exhibit varying generalization properties that depend on implicit bias of the optimization procedure [START_REF] Belkin | Does data interpolation contradict statistical optimality?[END_REF][START_REF] Neyshabur | In search of the real inductive bias: On the role of implicit regularization in deep learning[END_REF]. To date, the majority of studies on generalization in over-parameterized networks have primarily focused on analyzing the dynamics of gradient descent and gradient flows -their continuous-time limit -, making important progress in understanding the implicit bias of these optimization procedures. Specifically, [START_REF] Chizat | On lazy training in differentiable programming[END_REF] 1. We provide a global convergence result in Proposition 3 with a linear rate for the continuous-time limit of GN's dynamics in the over-parameterized regime that holds under the mean-field limit.

The achieved rate shows better conditioning than the gradient flow with the same objective, emphasizing the method's advantage from a pure optimization perspective. While similar rates can be achieved under strong non-degeneracy conditions as in [START_REF] Yudin | Flexible modification of gauss-newton method and its stochastic extension[END_REF], these conditions are never met over-parameterized networks. Instead, our result relies on a simple condition on the initial parameters of a similar nature as in [START_REF] Chizat | Sparse Optimization on Measures with Over-parameterized Gradient Descent[END_REF]Thm 3.3]: (i) that the hidden weights are diverse enough so that an interpolating solution can be obtained simply by fitting the linear weights, and (ii) that the linear weights result from a simple fine-tuning procedure to reach near optimality. The result shows that GN can reach the data interpolation regime, often required for studying the implicit bias of optimization (Section 4).

2. The empirical study in Section 5 complements the convergence result in Proposition 3 by further investigating the implicit bias of GN on a student/teacher regression task using synthetic data. We show that, while GN attains an interpolating solution, as predicted by Proposition 3, such a solution generalizes more or less well depending on the regime under which the GN dynamics is running. Specifically, we show that the GN dynamics can exhibit both kernel and feature learning regimes depending on the choice of step size and variance of the weights at initialization. Our experiments also indicate that a hidden learning phenomenon occurs in GN dynamics where hidden features that generalize well are obtained even when both raw training and test performances are sub-optimal due to an under-optimized linear layer. Quite surprisingly, these features are found when using small step sizes and, sometimes, generalize better than those learned by gradient descent. This is contrary to the common prescription of using larger step sizes in gradient descent for better generalization [START_REF] Andriushchenko | SGD with large step sizes learns sparse features[END_REF]. Our results suggest a tradeoff between optimization speed and generalization when using GN which is of practical interest.

Related work

Convergence of Gauss-Newton's methods. A large body of work on inverse problems studied the local convergence of Gauss-Newton and its regularized versions such as the Levenberg-Marquardt method [START_REF] Kaltenbacher | Iterative regularization methods for nonlinear ill-posed problems[END_REF][START_REF] Nesterov | Modified Gauss-Newton scheme with worst case guarantees for global performance[END_REF]. Many such works focus on a continuous-time analysis of the methods as it allows for a simplified study [START_REF] Kaltenbacher | Convergence rates of the continuous regularized Gauss-Newton method[END_REF][START_REF] Ramm | Dynamical systems method for solving operator equations[END_REF]. In this work, we also consider continuous-time dynamics and leave the study of discrete-time algorithms to future work. More recently, increasing attention was given to global convergence analysis. [START_REF] Yudin | Flexible modification of gauss-newton method and its stochastic extension[END_REF] studied the global convergence of GN with an inexact oracle and in a stochastic setting establishing linear convergence. The analysis requires that the smallest singular value of the Jacobian is positive near initialization, a situation that can never occur in the over-parameterized regime. [START_REF] Mishchenko | Regularized Newton Method with Global O(1/k * * 2) Convergence[END_REF] studied the convergence of the Levenberg-Marquardt dynamics, a regularized version of Gauss-Newton's method, in a non-convex setting under a cubic growth assumption and showed global albeit sub-linear convergence rates. It is unclear, however, if such cubic growth assumption holds in our setting of interest. Closest to our work is [START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF] which aims to accelerate optimization by solving a kernel regression in the NTK scaling limit. There, the authors show global convergence of a regularized version of GN in the kernel regime. We aim to investigate the behavior of GD in the feature learning regimes.

The implicit bias of gradient flows. In the over-parameterized setting, the training objective often possesses several global minima that perfectly interpolate the training data. The generalization error thus heavily depends on the solution selected by the optimization procedure. [START_REF] Neyshabur | In search of the real inductive bias: On the role of implicit regularization in deep learning[END_REF] highlighted the importance of understanding the properties of the selected solutions for a given initialization and several works studied this problem in the case of linear neural networks [START_REF] Arora | Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks[END_REF][START_REF] Min | On the explicit role of initialization on the convergence and implicit bias of overparametrized linear networks[END_REF][START_REF] Pesme | Implicit bias of sgd for diagonal linear networks: a provable benefit of stochasticity[END_REF][START_REF] Woodworth | Kernel and rich regimes in overparametrized models[END_REF][START_REF] Yun | A unifying view on implicit bias in training linear neural networks[END_REF]. For nonlinear networks, there is still no general characterization of the implicit bias in the regression setting, although [START_REF] Boursier | Gradient flow dynamics of shallow ReLU networks for square loss and orthogonal inputs[END_REF] recently made important progress in the case of a one-hidden layer network with ReLU activation and orthogonal inputs showing that gradient flows select a minimum norm interpolating solution when the initial weights are close to zero. Recently, the implicit bias of gradient flows was shown to play a crucial role in other optimization problems such as non-convex bi-level optimization [START_REF] Arbel | Non-Convex Bilevel Games with Critical Point Selection Maps[END_REF][START_REF] Vicol | Implicit Regularization in Overparameterized Bilevel Optimization[END_REF], therefore illustrating the ubiquity of such phenomenon in deep learning. The present work empirically studies the implicit bias of Gauss-Newton which is vastly understudied.

3 Setup and preliminaries

Regression using over-parameterized one-hidden layer networks

We consider a regression problem where the goal is to approximate an unknown real-valued function f ⋆ defined over a subset X of R d from a pair of i.i.d. inputs/outputs (x n , y n ) 1≤n≤N where each x n is a sample from some unknown distribution P and y n = f ⋆ (x n ). We assume that f ⋆ belongs to L 2 (P), the set of square-integrable functions w.r.t. P and note that f ⋆ always belongs to L 2 ( P), where P denotes the empirical distribution defined by the samples (x 1 , . . . , x N ). For simplicity, we assume the data are non-degenerate, meaning that x n ̸ =x n ′ whenever n ̸ = n ′ . In this case L 2 ( P) is isomorphic to R N (i.e., L 2 ( P) ∼ = R N ) and we can identify any function f ∈ L 2 ( P) with the evaluation vector (f (x 1 ), . . . , f (x N )). We are interested in approximating f ⋆ using a one-hidden layer network f w with parameter vector w belonging to some, possibly infinite, Hilbert space W. The network's parameters are learned by minimizing an objective of the form L(f w ) where L is a non-negative, L-smooth, and µ-strongly convex function defined over L 2 ( P) and achieving a 0 minimum value at f ⋆ . A typical example that we consider in Section 5 is the mean-squared error over the training set:

min w∈W L(f w ), L(f w ) = 1 2N N n=1 (f w (x n ) -y n ) 2 . ( 1
)
For the convergence results, we will consider both finite-width one-hidden layer networks and their mean-field infinite-width limit. In the experiments, we will restrict to finite-width networks although, we will be using a relatively large number of units M compared to the size of the training data, therefore approximating the mean-field limit.

Finite-width one-hidden layer networks. Given a non-polynomial point-wise activation function γ and some positive integer M , these networks take the form:

f w (x) = 1 M M i=1 v i γ(u ⊤ i x), w = (v, u) ∈ W := R M × R M ×d , (2)
where v ∈ R M is the linear weight vector while u ∈ R M ×d is the hidden weight matrix. Popular choices for γ include ReLU [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] and its smooth approximation SiLU(x) = x(1 + e -βx ) -1 which can be made arbitrary close to ReLU by increasing β [START_REF] Elfwing | Sigmoid-Weighted Linear Units for Neural Network Function Approximation in Reinforcement Learning[END_REF]. The above formulation can also account for a bias term provided the input vector x is augmented with a non-zero constant component.

Mean-field infinite-width limit. We consider some base probability µ measure with full support on X and finite second moment and denote by L 2 (µ) and L 2 (µ, R d ) the set of square integrable functions w.r.t. µ and taking values in R and R d . Given a non-polynomial point-wise non-linearity, we define infinitely-wide one-hidden layer networks f w as:

f w (x) = v(c)γ(u(c) ⊤ x) dµ(c), w := (v, u) ∈ W := L 2 (µ) × L 2 (µ, R d ). (3)
Functions of the form (3) correspond to the mean-field limit of the network defined in [START_REF] Boursier | Gradient flow dynamics of shallow ReLU networks for square loss and orthogonal inputs[END_REF] when the number of units M grows to infinity and appear in the Lagrangian formulation of the Wasserstein gradient flow of infinitely wide networks [START_REF] Wojtowytsch | On the convergence of gradient descent training for two-layer relu-networks in the mean field regime[END_REF]. Existing global convergence results of GN methods were obtained in the context of the NTK limit [START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF] which does not allow feature learning. The meanfield limit we consider here can result in optimization dynamics that exhibit feature learning which is why we consider it in our global convergence result analysis. For completeness, we provide a brief discussion on mean-field vs NTK limits and kernel vs feature learning regimes in Appendix A.

Over-parameterization. We consider an over-parameterized setting where the network has enough parameters to be able to fit the training data exactly, thus achieving a 0 training error. To this end, we define the gram matrix G which is an N × N matrix taking the following forms depending on whether we are using a finite-width network (G F ) or infinite-width network (G I ):

G F n,n ′ (u) = 1 M M i=1 γ(u ⊤ i x n )γ(u ⊤ i x n ′ ), G I n,n ′ (u) = γ(u(c) ⊤ x n )γ(u(c) ⊤ x n ′ ) dµ(c).
We say that a network is over-parameterized if G(u 0 ) is invertible for some hidden-layer weight u 0 :

(A) (Over-parameterization) G(u 0 ) is invertible for an initial parameter w 0 =(v 0 , u 0 )∈W.

When u 0 is sampled randomly from a Gaussian, a common choice to initialize a neural network, Assumption (A) holds for infinitely wide networks as long as the training data are non-degenerate (see [START_REF] Du | Gradient descent provably optimizes overparameterized neural networks[END_REF]Theorem 3.1] for ReLU, and [17, Lemma F.1] for analytic non-polynomial activations). There, the non-degeneracy assumption simply means that the inputs are not parallel to each other (i.e. x i ∦ x j whenever i ̸ = j for any 1 ≤ i, j ≤ N ), a property that holds almost surely if the data distribution P has a density w.r.t. Lebesgue measure. In the case of finite-width networks with M > N , the result still holds with a high probability when u is sampled from a Gaussian [START_REF] Liu | Loss landscapes and optimization in over-parameterized nonlinear systems and neural networks[END_REF].

Generalized Gauss-Newton dynamics

To solve (1), we consider optimization dynamics based on a generalized Gauss-Newton method.

To this end, we denote by J w the Jacobian of (f w (x 1 ), ..., f w (x N )) w.r.t. parameter w which can be viewed as a linear operator from W to R N . Moreover, ∇L(f w ) denotes the vector of size N representing the gradient of L w.r.t. the outputs (f w (x 1 ), ..., f w (x N )). We can then introduce the Gauss-Newton vector field Φ : W → W defined as:

Φ(w) := (J ⊤ w H w J w + ε(w)I) -1 J ⊤ w ∇L(f w ),
where H w is a symmetric positive operator on L 2 ( P) and ε(w) is an optional non-negative (possibly vanishing) damping parameter. Starting from some initial condition w 0 ∈ W and for a given positive step-size γ, the Gauss-Newton updates and their continuous-time limit are given by:

w k+1 = w k -γΦ(w k ), ẇt = -Φ(w t ). (4)
The continuous-time limit is obtained by taking the step-size γ to 0 and rescaling time appropriately. When H w is given by the Hessian of L, the dynamics in (4) recovers the generalized Gauss-Newton dynamics in continuous time when no damping is used [START_REF] Schraudolph | Fast curvature matrix-vector products for second-order gradient descent[END_REF] and recovers the continuoustime Levenberg-Maquardt dynamics when a positive damping term is added [START_REF] Moré | The Levenberg-Marquardt algorithm: implementation and theory[END_REF]. When the matrix H w is the identity H w = I, the resulting dynamics is tightly related to the natural gradient [START_REF] Amari | Natural Gradient Works Efficiently in Learning[END_REF][START_REF] Martens | New insights and perspectives on the natural gradient method[END_REF]. More generally, we only require the following assumption on H w :

(B) H w is continuously differentiable in w with eigenvalues in [µ H , L H ] for positive L H , µ H .

Optional damping. Amongst possible choices for the damping, we focus on ε(w) of the form:

ε(w) = ασ 2 (w) (5)
where α is a non-negative number α ≥ 0, and σ 2 (w) is the smallest eigenvalue of the Neural Tangent Kernel (NTK) A w := J w J ⊤ w which is invertible whenever J w is surjective. This choice allows us to study the impact of scalar damping on the convergence rate of the GN dynamics. While computing σ 2 (w) is challenging in practice, we emphasize that the damping is only optional and that the case when no damping is used, i.e. when α = 0, is covered by our analysis.

Convergence analysis 4.1 Global convergence under no-blow up

We start by showing that the continuous-time dynamics in (4) converge to a global optimum provided that it remains defined at all times. Proposition 1. Let w 0 = (v 0 , u 0 ) ∈ W be an initial condition so that u 0 satisfies Assumption (A). Under Assumption (B) and assuming that the activation γ is twice-continuously differentiable, there exists a unique solution to the continuous-time dynamics in (4) defined up to some maximal, possibly infinite, positive time T . If T is finite, we say that the dynamics blows-up in finite time. If, instead T is infinite, then f wt converges to a global minimizer of L. Moreover, denoting by µ the strong convexity constant of L and defining µ GN := 2µ/(L H (1 + α/µ H )), it holds that:

L(f wt ) ≤ L(f w0 )e -µ GN t . ( 6
)
Proposition 1 is proved in Appendix B.2. It relies on the Cauchy-Lipshitz theorem to show that the dynamics are well-defined up to a positive time T and uses strong convexity of L to obtain the rate in [START_REF] Arora | Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks[END_REF]. The result requires a smooth non-linearity, which excludes ReLU. However, it holds for any smooth approximation to ReLU, such as SiLU. As we discuss in Section 5, this does not make a large difference in practice as the approximation becomes tighter. The rate in (6) is only useful when the dynamics is defined at all times which is far from obvious as the vector field Φ can diverge in a finite time causing the dynamics to explode. Thus, the main challenge is to find conditions ensuring the vector field Φ remains well-behaved which is the object of Section 4.2.

Comparaison with gradient flow. The rate in [START_REF] Arora | Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks[END_REF] only depends on the strong-convexity constant µ of L, the smoothness constant of H, and damping strength α, with the fastest rate achieved when no damping is used, i.e. α = 0. For instance, when choosing H to the identity and α = 0, the linear rate becomes µ GN = 2µ. This is in contrast with a gradient flow of w → L(f w ) for which a linear convergence result, in the kernel regime, follows from [15, Theorem 2.4]:

L(f wt ) ≤ L(f w0 )e µ GD t ,
with µ GF := µσ(w 0 )/4 where σ(w 0 ) is the smallest singular value of the initial NTK matrix A w0 .

In practice, µ GF ≪ µ GN since the linear rate µ GF for a gradient flow is proportional to σ(w 0 ), which can get arbitrarily small as the training sample size N increases thus drastically slowing down convergence [START_REF] Robin | Convergence beyond the over-parameterized regime using Rayleigh quotients[END_REF].

Absence of blow-up for almost-optimal initial linear layer

In this section, we provide a condition on the initialization w 0 = (v 0 , u 0 ) ∈ W ensuring the Gauss-Newton dynamics never blows up. We start with a simple result, with a proof in Appendix B.1, showing the existence of a vector v ⋆ so that w ⋆ := (v ⋆ , u 0 ) ∈ W interpolates the training data whenever the over-parameterization assumption Assumption (A) holds for a given u 0 .

Proposition 2. Let w 0 = (v 0 , u 0 ) ∈ W be an initial condition so that u 0 satisfies Assumption (A). Define Γ(u) := γ(u ⊤ x n ) 1≤n≤N and F ⋆ := (f ⋆ (x 1 ), ..., f ⋆ (x N )) and let v ⋆ be given by:

v ⋆ = Γ(u 0 )G(u 0 )F ⋆ .
Then the vector w ⋆ = (v ⋆ , u 0 ) ∈ W is a global minimizer of w → L(f w ) and a fixed point of the Gauss-Newton dynamics, i.e. Φ(w ⋆ ) = 0.

The above result ensures a global solution can always be found by optimizing the linear weights while keeping the hidden weights fixed. In practice, this can be achieved using any standard optimization algorithm such as gradient descent since the function v → l(v) = L(f (v,u0) ) converges to a global optimum as it satisfies a Polyak-Lojasiewicz inequality (see Proposition 4 in Appendix B.1). However, this procedure is not of interest here as it does not allow us to learn the hidden features. Nonetheless, since w ⋆ = (v ⋆ , u 0 ) is a fixed point of the Gauss-Newton dynamics, we can expect that an initial condition w 0 = (v 0 , u 0 ), where v 0 is close to v ⋆ , results in a well-behaved dynamics.

The following proposition, proven in Appendix B.4, makes the above intuition more precise. Proposition 3. Let w 0 = (v 0 , u 0 ) ∈ W be an initial condition so that Assumption (A) holds, i.e. the gram matrix G(u 0 ) is invertible, and denote by σ 2 0 the smallest eigenvalue of G(u 0 ). Assume that the activation γ is twice-continuously differentiable and that Assumption (B) holds. Let R be any arbitrary positive number and define C R = sup w∈B(w0,R) ∥∂ w J w ∥ op where ∥.∥ op denotes the operator norm. If the linear layer is almost optimal v 0 in the sense that:

∥∇L(f w0 )∥ < ϵ, ϵ := (µµ H µ GN /8LN ) min(R, C -1 R ) min σ 0 , σ 2 0 , (7)
then ( 4) is defined at all times, i.e. T = + ∞, the objective L(f wt ) converges at a linear rate to 0 according to [START_REF] Arora | Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks[END_REF] and the parameters w t remain within a ball of radius R centered around w 0 .

Proposition 3 essentially states that the dynamics never blow up, and hence, converge globally at the linear rate in (6) by Proposition 1, provided the hidden features are diverse enough and the initial linear weights optimize the objective well enough. As discussed in Section 3.1, the first condition on the hidden weights u 0 typically holds for a Gaussian initialization when the data are non-degenerate. Additionally, the near-optimality condition on the linear weights v 0 can always be guaranteed by optimizing the linear weights while keeping the hidden ones fixed.

Remark 1. In the proof of Proposition 3, we show that the occurrence of a finite-time blow-up is tightly related to the Neural Tangent Kernel matrix A wt becoming singular at t increases which causes the vector field Φ to diverge. When the NTK matrix A wt at time t is forced to remain close to the initial one A w0 , as is the case in the NTK limit considered in [START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF], one can expect that the singular values of A wt remain bounded away from 0 since this is already true for A w0 . Therefore, the main technical challenge is the study, performed in Propositions 7 to 9 of Appendix B.3, of the time evolution of the eigenvalues of A wt in the mean-field limit, where A wt is allowed to differ significantly from the initial NTK matrix A w0 . The final step is to deduce that, while the weights w t are allowed to be in a ball centered in w 0 of arbitrarily large radius R, the singular values remain bounded away from 0 provided the initial linear weights v 0 satisfy the condition in [START_REF] Belkin | Does data interpolation contradict statistical optimality?[END_REF].

Remark 2. Comparing Proposition 3 to the convergence results for gradient flows in the kernel regime [15, Theorem 2.4], our result also holds for an infinite dimensional space of parameters such as in the mean-field limit of one-hidden layer network. Moreover, while [15, Theorem 2.4] requires the norm w t to be within a ball of fixed radius R 0 determined by initial smallest singular value σ(w 0 ) of A w0 , our result allows w t to be within arbitrary distance R regardless of the initial σ(w 0 ) provided the condition in (7) holds. Proposition 3 is more similar in flavor to the convergence result of the Wasserstein gradient flow for some probability functionals in [13, Theorem 3.3, Corollary 3.4] which requires the initial objective to be smaller than a given threshold.

While the results of this section guarantee global convergence of the training objective, they do not guarantee learning features that generalize well. Next, we empirically show that GN can exhibit regimes where it yields solutions with good generalization properties.

Empirical study of generalization for a Gauss-Newton algorithm

We perform a comparative study between gradient descent (GD) and Gauss-Newton (GN) method in the context of over-parameterized networks. Since we observed little variability when changing the seed, we have chosen to use a single seed for the sake of clarity in the figures when presenting the results of this section, while deferring the analysis for multiple seeds to Appendix C.1. Additional experiments using MNIST dataset [START_REF] Deng | The mnist database of handwritten digit images for machine learning research[END_REF] are provided in Appendix C.4. The results presented below were obtained by running 720 independent runs, each of which optimizes a network given a specific configuration on a GPU. The total time for all runs was 3600 GPU hours.

Experimental setup

We consider a regression task on a synthetic dataset consisting of N training points

(X n , Y n ) 1≤n≤N .
The objective is to minimize the mean-squared error L(f w ), defined in (1), over the parameters w of a model f w that predicts the target values Y n based on their corresponding input points X n .

Data generation. We generate N i.i.d. samples denoted as (X n ) 1≤n≤N from a standard Gaussian distribution of dimension d = 10. Each corresponding target Y n is obtained by applying a predefined function f ⋆ , referred to as the teacher network, to the input X n (i.e., Y n = f ⋆ (X n )). We choose f ⋆ to be a one-hidden layer network of the form in [START_REF] Boursier | Gradient flow dynamics of shallow ReLU networks for square loss and orthogonal inputs[END_REF] with M ⋆ = 5 hidden units (v ⋆ i , u ⋆ i ) 1≤i≤M ⋆ drawn independently from a standard Gaussian distribution. Furthermore, we consider two nonlinearities γ when defining f ⋆ : the widely used ReLU activation [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] for the main results and its smooth approximation SiLU for additional ablations in Appendix C.3. This choice of target function f ⋆ results in a hidden low-dimensional structure in the regression task as f ⋆ depends only on 5dimensional linear projection of the input data. In most cases, we take the size of the training data to be N = 500 except when studying the effect of the training size. This choice allows us to achieve a balance between conducting an extensive hyper-parameter search and being able to compute precise GN updates on the complete training data within a reasonable timeframe. Finally, we generate 10000 test samples to measure the generalization error.

Model. We consider a well-specified setting where the model w → f w , referred to as the student network, is also a one-hidden layer with the same activation function γ as the teacher network f ⋆ . Importantly, the student network possesses a number M = 5000 of hidden units, denoted as w := (v i , u i ) 1≤i≤M , which is much larger than the teacher's number of units M ⋆ = 5 and allows fitting the training data perfectly (over-parameterized regime). Following [START_REF] Chizat | On lazy training in differentiable programming[END_REF], we also normalize the output of the student f w by M to remain consistent with the mean-field limit as M increases.

Initialization. We initialize the student's hidden units according to a centered Gaussian with standard deviation (std) τ 0 ranging from 10 -3 to 10 3 . Varying the initial std τ 0 allows us to study the optimization dynamics in two regimes: the kernel regime (large values of τ 0 ) and the feature learning regime (small values of τ 0 ). Finally, we initialize the weights of the last layer to be 0. This choice allows us to perform a systematic comparison with the minimum norm least squares solution obtained using random features as described next.

Methods and baselines. We consider three optimization methods for the model's parameters: a regularized version of Gauss-Newton (GN), gradient descent (GD), and optimizing the linear layer parameters alone, which can be viewed as a random features (RF) model. (GN): We use the discrete GN updates in (4) with a constant step-size λ and H w = I. Each update is obtained using Woodbury's matrix identity by writing

Φ(w k )=J ⊤ w k z k with z k the solution of a linear system J w k J ⊤ w k +ϵ(w k )I z k =∇L(f w k ) of size N .
Here, we use the damping defined in (5) with α=1 and ensure it never falls below ϵ 0 = 10 -7 to avoid numerical instabilities. In practice, we found that small values of ϵ 0 had little effect on the results (see Figure 4 (Right) of Appendix C.1). (GD): The model's parameters w are learned using gradient descent with a constant step-size λ. (RF): Instead of optimizing the entire parameter vector w = (v, u) of the model f w , we selectively optimize the parameter v of the linear layer while keeping the weights u of the hidden layer constant. This procedure corresponds to computing a minimal-norm least squares solution v RF for a random features (RF) model, where the features are obtained from the output of the hidden layer at initialization. Specifically, the solution v RF is obtained as

v RF = Γ(u)(Γ(u) ⊤ Γ(u)) † Y, (8)
where

Γ(u) = (γ(u ⊤ X n )) 1≤n≤N ∈ R M ×N
is the feature vector computed over the training data, and Y is a vector of size N consisting of the target values Y n . The choice of the un-regularized solution as a baseline is supported by recent findings [START_REF] Hastie | Surprises in high-dimensional ridgeless least squares interpolation[END_REF], demonstrating its good generalization in the over-parameterized regime.

Stopping criterion. For both (GD) and (GN), we perform as many iterations as needed so that the final training error is at least below 10 -5 . Additionally, we stop the algorithm whenever the training error drops below 10 -7 or when a maximum number of iterations of K GD = 10 6 iterations for (GD) and K GN = 10 5 iterations for (GN) is performed. For (RF) we solve the linear system exactly using a standard linear solver.

Performance metrics

In addition to the training and test losses, we introduce two complementary metrics to assess the quality of the hidden features Γ(u) learned by the student model.

Weighted cosine distance (WCD). Measuring proximity between the student and teacher's hidden weights (u 1 , ..., u M ) and (u ⋆ 1 , ..., u ⋆ M ⋆ ) requires being able to find a correspondance between each student's parameter u i and a teacher's one u ⋆ j . When using a positively homogeneous non-linearity such as ReLU, only the alignment between vectors is relevant. However, since the student's model is vastly over-parameterized, there could be vanishing weights u i = 0 which do not influence on the network's output, while the remaining weights are well aligned with the teacher's weights. We introduce the following weighted distance to account for these specificities:

WCD(u, u ⋆ ) = 2 M i=1 p i 1 -max 1≤j≤M ⋆ u ⊤ i u ⋆ j ∥u i ∥ u ⋆ j , p i = ∥u i ∥ 2 M k=1 ∥u k ∥ 2 . (9)
Equation ( 9) finds a teacher's parameter u ⋆ j that is most aligned with a given student's parameter u i and downweights its cosine distance if u i has a small norm. In practice, we found this measure to be a good indicator for generalization.

Test loss after linear re-fitting (Test-LRfit). To evaluate the relevance of the hidden features Γ(u), we train a new linear model ṽ on those features and measure its generalization error. In practice, we freeze the weights u of the hidden layer and fit the last layer v using the procedure (8) described above for the random features (RF) baseline. In our over-parameterized setting, this new model should always achieve perfect training loss, but its generalization error strongly depends on how the features were learned.

Numerical results

Implicit bias of initialization. Figure 1(Left) shows that GN and GD generalize more or less well compared to a random features model depending on the variance τ 0 of the weights at initialization. First, in the case of the RF model, changing τ 0 barely affects the final training and test error. This is essentially due to the choice of the non-linearity ReLU which is positively homogeneous. Second, for large τ 0 , the final test error of both GN and GD matches that of RF suggesting that the dynamics were running under the kernel regime/lazy regime [START_REF] Chizat | On lazy training in differentiable programming[END_REF][START_REF] Jacot | Neural Tangent Kernel: Convergence and Generalization in Neural Networks[END_REF] where the final layer's weights are learned without changing the hidden features much. Finally, for small values of τ 0 , both GN and GD obtained a better test error than RF which can only be explained by learning suitable hidden features. These results are further confirmed when varying the size of the training set as shown in Appendix C.2. In Appendix C.3, we performed the same experiment using SiLU non-linearity instead of ReLU and observed the same transition between the kernel and feature learning regimes. While prior works such as [START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF] analyzed Gauss-Newton in the kernel learning regime, our results indicate that Gauss-Newton can also exhibit a feature learning regime for small values of τ 0 .

Implicit bias of the step size. Figure 1(Right) shows that increasing the step size in Gauss-Newton results in features that do not generalize well. This is unlike gradient descent where larger step sizes yield better-performing features 3 . This behavior is first observed in the top figures showing an increasing weighted cosine distance (WCD) between the student's hidden weights and the teacher's weights as the step size increases in the case of Gauss-Newton. On the contrary, the distance decreases with larger step sizes in the case of gradient descent, indicating that the algorithm learns better features. This effect is also confirmed when computing the test loss after refitting the linear layer to the training data exactly while keeping the learned hidden weights fixed (see linear re-fitting in Section 5.2). Several works, such as [START_REF] Andriushchenko | SGD with large step sizes learns sparse features[END_REF][START_REF] Mulayoff | The implicit bias of minima stability: A view from function space[END_REF], highlight the importance of taking large step-sizes Step size for gradient descent, to improve generalization, our results show that Gauss-Newton benefits instead of using smaller step-sizes. As shown in Appendix C.3, this behavior persists when using a different non-linearity, such as SiLU with different values for β (1 and 10 6 ). Interestingly, we observed that GN might underperform GD as the parameter β defining SiLU decreases making SiLU nearly linear and less like ReLU. In all cases, while the final test loss remains large for small step sizes (in both GD and GN), the test loss after linear re-fitting is much lower in the case of Gauss-Newton, indicating that the algorithm was able to learn good features while the last layer remained under-optimized.

Hidden feature learning. Figure 2 (Left) shows the evolution of various metrics, including the WCD, the training, and test errors with the number of iterations for both GN and GD. For each method, results for the best step size and variance at initialization are presented ((τ 0 , λ) = (10 -3 , 10 2 ) for GD and (τ 0 , λ) = (10 -3 , 10 -2 ) for GN). The bottom figure clearly indicates that the early iterations of GN essentially optimize internal features first while barely improving the last layer. This can be seen deduced from the training and test losses which remain almost constant at the beginning of optimization, while the test loss after linear re-fitting steadily decreases with the iterations. The last layer is learned only towards the end of training as indicated by a faster decrease in both train and test objectives. Gradient descent displays a different behavior with test loss before and after linear re-fitting remaining close to each other during optimization.

Better generalization of GN comes at the cost of a slower training. Figure 2 (Right) shows the evolution of the same quantities as in the left figures for τ 0 = 10 -3 as a function of time (in seconds) and for the three values of the step-size λ (10 -2 (bright colors) 1 (medium colors) and 10 2 (dark colors)). The figure shows that, while the training loss converges much faster for larger step sizes in the case of GN, it does not result in the best generalization error as measured by the test loss after linear re-fitting. Hence, better generalization comes at the cost of a slower convergence of the training loss. Again, the behavior is quite different in the case of gradient descent for which larger step sizes result in both faster optimization of the training objective and better generalization.

Conclusion

The results illustrate the importance of the initialization and choice of the step size for the generalization performance of GN when optimizing over-parameterized one-hidden layer networks. While Step size 1e2 1e0 1e 2

Step size our theoretical analysis shows that GN can reach a global solution of the training objective when starting from a close-to-optimal initialization, this result does not imply anything about the generalization properties of the obtained solution or the quality of the learned features. Our empirical study instead shows that GN may favor feature learning (thus yielding improved generalization) at the cost of a slower optimization of the training objective due to the use of a smaller step size. Providing a theoretical analysis of these phenomena is an interesting direction for future work. Finally, our study shows that the training or test error may not always be a good indicator of the quality of the learned features due to an under-optimized final layer. Instead, a test/validation error after re-fitting the linear layer may be a better indicator for the quality of the learned features which can be of practical use.

A Background: NTK vs mean-field limits, kernel vs feature learning regimes

We provide a brief discussion of the notions of NTK and mean-field limits as well as the related convergence regimes for a gradient flow. Given a non-polynomial point-wise and differentiable activation function γ and some positive integer M , consider the following one-hidden network:

f w (x) = α(M ) M i=1 v i γ(u ⊤ i x), w = (v, u) ∈ W := R M × R M ×d , (10)
where α(M ) is a scaling factor that depends on the number of neurons M . For simplicity, we assume that the weights (v i , u i ) are i.i.d. samples from a Gaussian. We are interested in the dynamics of a gradient flow when M increases to +∞ and for specific choices of the scaling factor α(M ). First, let's introduce some notations. We denote by F (w) the vector of evaluation of the network on the training data F (w) = (f w (x 1 ), . . . , f w (x N )) and by J w the Jacobian of F (w) w.r.t. parameter w. Moreover, we introduce the NTK matrix A w := J w J ⊤ w which is an N × N symmetric positive semi-definite matrix. Finally, ∇L(f w ) denotes the vector of size N representing the gradient of L w.r.t. F (w). The gradient flow of an objective L(F w ) is given by: ẇt = -J wt ∇ F L(F wt ), Additionally, it is often informative to track the time evolution of the vector F wt , which is given by the following ODE:

∂ t F wt = -J wt J ⊤ wt ∇ F L(F wt ) = -A wt ∇ F L(F wt ). (11)
One can see from [START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF] that the NTK matrix A wt arises naturally when expressing the dynamics of the function F wt . We can even provide an explicit expression of A wt using the NTK kernel, which is a positive semi-definite kernel K w (x, x ′ ) defined by the inner product:

K w (x, x ′ ) = ∂ w f (x)∂ w f (x ′ ) ⊤ = α(M ) 2 M i=1 γ(u ⊤ i x)γ(u ⊤ i x ′ ) + (v i ) 2 γ ′ (u ⊤ i x)γ ′ (u ⊤ i x ′ )x ⊤ x ′ .
The NTK matrix A w is then formed by collecting the values (K w (x i , x j )) 1≤i,j≤N on the training data. It is easy to see that when α(M ) = 1 √ M , then K w (x, x ′ ) converges to an expectation under the distribution of the weights (v i , u i ) 1≤i≤M under mild conditions. In this case, A w depends on the choice of the distribution of the weights (v i , u i ) 1≤i≤M and not on the particular samples. On the other hand, when α(M ) = 1 M , it is the function f w itself that converges to an expectation under the distribution of the weights. This distinction will lead to different qualitative dynamics for the gradient flows as we discuss next.

A.1 The Neural Tangent Kernel limit and the kernel regime (α(M

) = 1 √ M )
When choosing α(M ) = 1 √ M and setting M → +∞, it is shown in [START_REF] Jacot | Neural Tangent Kernel: Convergence and Generalization in Neural Networks[END_REF] that the NTK matrix A wt remains constant in time and depends only on the distribution of the weights w 0 at initialization. Hence, F wt is effectively the solution of the ODE:

∂ t F t = -A w0 ∇ F L(F t ), F 0 = F w0
, which converges at a linear rate to a global minimizer of L whenever A w0 is invertible and L is strongly convex. The NTK limit results in a kernel regime, where the resulting solution is equivalent to the one obtained using a kernel method with the initial NTK kernel K w0 (x, x ′ ) [START_REF] Chizat | On lazy training in differentiable programming[END_REF][START_REF] Jacot | Neural Tangent Kernel: Convergence and Generalization in Neural Networks[END_REF].

A.2 The mean-field limit and the feature learning regime (α(M ) = 1 M )

When choosing instead α(M ) = 1 M , it is possible to view the limiting function as an expectation over weights drawn from some probability distribution ρ f ρ (w) = vγ(u ⊤ x) dρ(v, u).

In [START_REF] Wojtowytsch | On the convergence of gradient descent training for two-layer relu-networks in the mean field regime[END_REF], it is further shown that such a formulation is essentially equivalent to the one considered in [START_REF] Andriushchenko | SGD with large step sizes learns sparse features[END_REF]. As a result, the time evolution in [START_REF] Cai | Gram-gaussnewton method: Learning overparameterized neural networks for regression problems[END_REF] does not simplify as in the NTK limit, since the NTK matrix A wt (and NTK kernel K wt ) is allowed to evolve in time. This is a limit in which feature learning is possible since the dynamics do not necessarily simplify to a kernel method.

B Proofs

Notations. For an operator A between two Hilbert spaces, denote by ∥A∥ op its operator norm and by ∥A∥ F its Frobenius norm and let σ ⋆ (A) be its smallest singular value whenever it is well-defined. Denote by B(w, R) the open ball of radius R centered at an element w in a Hilbert space.

B.1 Over-parameterization implies interpolation.

Proof of Proposition 2. Denote by F ⋆ = (f ⋆ (x 1 ), ..., f ⋆ (x n )). First, it is easy to see that the matrix

∂ v f w ∂ v f ⊤
w is exactly equal to the N × N matrix G(u 0 ) considered in Assumption (A), which is assumed to be invertible. Moreover, ∂ v f w is independent of v since f (v,u0) is linear in v. Hence, we can simply choose:

v ⋆ = (∂ v f (v,u0) ) † F ⋆ = (∂ v f (v,u0) ) ⊤ (∂ v f (v,u0) ∂ v f ⊤ (v,u0) ) † F ⋆ = (∂ v f (v,u0) ) ⊤ G(u 0 ) -1 F ⋆ .
Moreover, using again the linearity of f in v, we have that

f (v ⋆ ,u0) = ∂ v f (v ⋆ ,u0) v ⋆ = G(u 0 )G(u 0 ) -1 F ⋆ = F ⋆ .
We have therefore shown that

w ⋆ = (v ⋆ , u 0 ) is a global minimizer of w → L(f w ). Finally, since ∇L(f w ⋆ ) = ∇L(f ⋆ ) = 0, it follows that ϕ(w ⋆ ) = 0.
Proposition 4. Assume that Assumption (A) holds for some vector u 0 . Then, the function v → l(v) = L(f (v,u0) ) is convex and satisfies a Polyak-Lojasiewicz inequality:

∥∇ u l(u)∥ 2 ≥ 2µσ 2 0 (l(u) -L(f ⋆ ))
. with µ the strong convexity constant of L and σ 2 0 the smallest singular value of the matrix G(u 0 ) appearing in Assumption (A).

Proof. It is easy to see that l(v) is convex as a composition of a strongly convex function L and a linear function v → f (v,u0) . We will show now that l satisfies a Polyak-Lojasiewicz (PL) inequality. To this end, we make the following calculations:

∥∇ u l(u)∥ 2 = ∇ f L(f w ) ⊤ ∂ v f w ∂ v f ⊤ w ∇ f L(f w ). ( 12
)
It is easy to see that the matrix ∂ v f w ∂ v f ⊤ w is exactly equal to the N × N matrix G(u 0 ) considered in Assumption (A), which is assumed to be invertible. Hence, its smallest singular value σ 2 0 is positive and the following inequality holds:

∥∇ u l(u)∥ 2 ≥ σ 2 0 ∥∇ f L(f w )∥ 2 . ( 13
)
On the other hand, since L is strongly convex and differentiable, then it satisfies the following PL inequality:

∥∇ f L(f w )∥ 2 ≥ 2µ(L(f w ) -L(f ⋆ )), ( 14
)
where µ is the strong convexity constant of L. Hence, we get the desired inequality by combining ( 13) and ( 14):

∥∇ u l(u)∥ 2 ≥ 2µσ 2 0 (l(u) -L(f ⋆ )).

B.2 Global convergence under no blow-up: Proof of Proposition 1

In this section we prove the global convergence result under no blow-up stated in Proposition 1. We start by stating the smoothness assumption on the activation function γ as we will be using it in many places.

(C) The non-linearity γ is twice-continuously differentiable.

As a first step, we need to prove a local existence and uniqueness of the solution to the ODE in (4) which results from Cauchy-Lipschitz theorem applied to the vector field Φ. To this end, a key step is to show that Φ is locally Lipschitz near initialization when the Jacobian J w0 is surjective.

Proposition 5 (Regularity of Φ.). Assume (B) and (C) and that w 0 satisfies Assumption (A). Then, there exists a neighborhood of w 0 so that Φ(w) is locally Lipschitz in w.

Proof of Proposition 5. By application of the Woodbury matrix identity, Φ can be equivalently written as:

Φ(w) = J ⊤ w J w J ⊤ w + ε(w)H -1 w -1 H -1 w ∇L(f w ). ( 15
)
Moreover, by the smoothness assumption, J w , H w , ϵ(w) and ∇L(f w ) are all differentiable in w.

Additionally, since J w0 is surjective, it must hold that σ 2 (w) > 0 in a neighborhood of w 0 so that J w J ⊤ w remains invertible. Finally, since H w is always positive, it holds that

J w J ⊤ w + ε(w)H -1 w -1
and H -1 w are differentiable in a neighborhood of w 0 . Therefore, Φ is Locally Lipschitz in a neighborhood of w 0 .

The next proposition establishes the existence and uniqueness of the dynamics defined in (4). Proposition 6 (Local existence). Under the same conditions as Proposition 5, there exists a unique solution (w t ) t≥0 to (4) defined up to a, possibly infinite, maximal time T > 0.

Proposition 6 is a direct consequence of the Cauchy-Lipschitz theorem which holds since the vector field Φ(w) is locally Lipschitz in w in a neighborhood of w 0 by Proposition 5. The above local existence result does not exclude situations where the dynamics blow up in finite time T < +∞ and possibly fail to globally minimize the objective L. We will show, later, that finite-time blow-ups never occur under additional assumptions on the initial condition.

Proof of Proposition 1. Existence and uniqueness of the continuous-time dynamics up to a maximal time T > 0 follows from Proposition 6. It remains to obtain the global convergence rate when T is infinite. For simplicity, we will ignore the dependence on w t in what follows and will write b := ∇L(f wt ). By differentiating the objective in time, we get the following expression:

∂ t [L(f wt )] = -b ⊤ J J ⊤ HJ + εI -1 Jb = -b ⊤ H -1 b + εb ⊤ H -1 JJ ⊤ + εH -1 -1 H -1 b = -b ⊤ H -1 b + εb ⊤ H -1 2 H -1 2 JJ ⊤ + εH -1 -1 H -1 2 H -1 2 b ≤ -b ⊤ H -1 b 1 -ε H -1 2 JJ ⊤ + εH -1 -1 H -1 2 op ≤ -b ⊤ H -1 b 1 - ε ε + σ ⋆ (H)σ ⋆ (JJ ⊤ ) ≥η:=(1+α/µ H ) -1 ≤ -ηb ⊤ H -1 b ≤ -ηL -1 H ∥b∥ 2 ≤ -2ηµL -1
H L(f wt ). We used the Woodbury matrix inversion lemma to get the second line while the third and fourth lines use simple linear algebra. To get the fifth line, we used Lemma 3 to upper-bound the operator norm appearing in the fourth line. The last line follows by definition of ε(w), the fact that H ≤ L H I and the fact that 1 2 ∥∇L(f )∥ 2 ≥ µL(f ) by µ-strong convexity of L. This allows to directly conclude that L(f wt ) ≤ L(f w0 )e -µ GN t , with µ GN := 2ηµL -1 H .

B.3 Control of the singular values of NTK matrix

In this section, we are interested in the evolution of the smallest singular value of the NTK matrix A wt := J wt J ⊤ wt . For any arbitrary positive radius R > 0, define the constant C R and stopping time T R as follows:

C R := sup w∈B(w0,R) ∥∂ w J w ∥ op < +∞ T R := sup {t ≥ 0 |∥w t -w 0 ∥ < R, and σ ⋆ (A wt ) > 0},
The stopping time T R characterizes the smallest-time when the dynamics becomes either degenerate (σ ⋆ (A wt ) = 0) or grows too far away from initialization (w t ̸ = B(w 0 , R)).

To control the singular values of A wt , we find it useful to consider the dynamics of the Frobenius norm of the pseudo-inverse J † w which we denote by a t := J † t F

. Indeed, simple linear algebra provides the following lower-bound on σ ⋆ (A wt ) in terms of a t :

σ ⋆ (A wt ) ≥ a -2 t . (16)
Note that, while the space of parameters can even be infinite-dimensional (mean-field limit), J w has a finite-dimensional range and thus always admits a pseudo-inverse J † w given by:

J † w = J ⊤ w (J w J ⊤ w ) † .
The next proposition provides a differential inequality for a t = J † wt F in terms of ∥ ẇt ∥ that holds for all times smaller than T R . Proposition 7. Assume (B) and (C) and that w 0 satisfies Assumption (A), then for any R > 0, the time T R is positive and smaller than the maximal time T , i.e. 0 < T R ≤ T . Moreover, for any t ∈ [0, T R ), a t satisfies:

| ȧt | ≤ a 2 t C R ∥ ẇt ∥,
Proof. Positive stopping time T R . By construction, T R must be smaller than T as it requires w t to be well-defined. Moreover, J w0 is surjective by assumption and J w is continuous in w, therefore J wt must also be surjective for small enough positive times t > 0. Additionally, it must hold that ∥w t -w 0 ∥ ≤ R for small positive times, by continuity of t → w t . This allows shows that T R > 0.

Dynamics of J † wt . By smoothness of J w it holds that J † wt is differentiable on the open interval (0, T R ) and satisfies the following ODE obtained by direct differentiation:

J † t = -J † t Jt J † t + (I -P t ) J⊤ t J t J ⊤ t -1 . (17)
In [START_REF] Du | Gradient descent finds global minima of deep neural networks[END_REF], we introduced notation J t := J wt and its time derivative Jt for simplicity and denote by P t , the projector:

P t := J ⊤ t J t J ⊤ t -1 J t .
Controlling | ȧt |. Taking the time derivative of a t and recalling the evolution (17) yields:

a t ȧt = -J † t ⊤ J † t Jt J † t .
where we used that J † t ⊤ (I -P t ) = 0 by Lemma 1 in Appendix B.5. Furthermore, taking the absolute value of the above equation, and recalling that the Frobenius norm is sub-multiplicative for the product of operators (Lemma 4 in Appendix B.5), we directly get:

a t | ȧt | ≤ J † t ⊤ J † t Jt J † t ≤ a 3 t Jt F . (18)
By the chain rule, we have that Jt = ∂ w J t ẇt therefore, Jt 18) yields the desired inequality after dividing by a t which remains positive all times t ∈ [0, T R ) since J t does not vanish.

F ≤ ∥∂ w J t ∥ op ∥ ẇt ∥. Moreover, since t < T R , then it must holds that w t ∈ B(w 0 , R) so that ∥∂ w J t ∥ op ≤ C R . Henceforth, Jt F ≤ C T ∥ ẇt ∥ and (
The next proposition provides an estimate of the time derivatives | ȧt | and ∥ ẇt ∥ up to time T R . For simplicity, we introduce the notation ∆ t = ∥∇L(f wt )∥. Proposition 8. Assume (B) and (C) and that w 0 satisfies Assumption (A). Then, for any R > 0, the time T R is positive and smaller than the maximal time T up to which (4) is defined, i.e. 0 < T R ≤ T . Moreover, define C := Lµ -1 µ -1 H ∆ 0 C R . Then a t and ∥ ẇt ∥ satisfy the following differential inequality on [0, T R ):

| ȧt | ≤ Ca 3 t e -µ GN 2 t ∥ ẇt ∥ ≤ CC -1 R a t e -µ GN
Proof of Proposition 8. By application of Proposition 7 to (w t ) 0≤t<T , we directly have that T R is positive and smaller than T and the following estimate holds for any t < T R :

| ȧt | ≤ a 2 t C R ∥ ẇt ∥, It remains to upper-bound ∥ ẇt ∥.
Controlling ∥ ẇt ∥. Ignoring the dependence on w t in the notations for simplicity and setting B := H 1 2 J, we have, by definition of the dynamical system (4):

∥ ẇt ∥ ≤ J ⊤ HJ + εI -1 J ⊤ ∇L(f ) = J ⊤ JJ ⊤ + εH -1 -1 H -1 ∇L(f ) = B ⊤ BB ⊤ + ε -1 H -1 2 ∇L(f ) ≤ BB ⊤ BB ⊤ + εI -2 1 2 F H -2 2 ∇L(f ) ≤ B † F ∥∇L(f )∥ ≤ µ -1 H J † F ∥∇L(f )∥ = µ -1 H a t ∆ t ,
| ȧt | ≤ µ -1 H C R a 3 t ∆ t . Finally, since t ≤ T , it follows from Proposition 1 that ∆ t ≤ L µ ∆ 0 e -µ GN
2 t which allows to conclude.

The next proposition controls a t and ∥w t -w 0 ∥ in terms of T R , C R , and other constants.

T + = -2µ -1 GN log 1 -(ζC R ∆ 0 a 2 0 ) -1 , with the convention that T + = +∞ if 1 ≥ ζC R ∆ 0 a 2 0 .
Then for any time t < T -:= min(T R , T + ) it holds that:

∥w t -w 0 ∥ ≤ ζ∆ 0 a 0 , a t ≤ a 0 1 + ζC R ∆ 0 a 2 0 e -µ GN 2 t -1 -1 2 (19)
Proof of Proposition 9. First, by Proposition 8, we know that a t satisfies the following differential inequality:

| ȧt | ≤ Ca 3 t e -µ GN 2 t , with C := L µ H µ ∆ 0 C R .
We will control a t by application of Grönwall's inequality. To this end, consider the ODE: ḃt = Cb 3 t e -µ GN 2 t , b 0 = a 0 > 0. We know by Lemma 5 in Appendix B.5 that b t is given by:

b t = b 0 1 + 4Cµ -1 GN b 2 0 e -µ GN 2 t -1 -1 2 , = a 0 1 + ζ2C R ∆ 0 µ -1 GN a 2 0 e -µ GN 2 t -1 -1 2
for all times t < T + . Therefore, by Grönwall's inequality, it must hold that a t ≤ b t for all times t < min(T R , T + ). Moreover, recalling now that ∥ ẇt ∥ ≤ Lµ -1 µ -1 H ∆ 0 a t e -µ GN 2 t by Proposition 8, we directly get that:

∥ ẇt ∥ ≤ C C R a 0 e -µ GN 2 t 1 + 4Cµ -1 GN b 2 0 e -µ GN 2 t -1 -1 2 . ( 20
)
Integrating the above inequality in time yields the following estimate on ∥w t -w 0 ∥ for any t < min(T R , T + ):

∥w t -w 0 ∥ ≤ t 0 ∥ ẇs ∥ ds ≤ 1 -1 + 4Cµ -1 GN a 2 0 e -µ GN 2 t -1 1 2 C R a 0 ≤ 4L(µµ H µ GN ) -1 ∆ 0 a 0 = ζ∆ 0 a 0 ,
where the second line follows by explicitly integrating the r.h.s. of (20) while the last line follows by the concavity of the square-root function.

Proposition 9 shows that the Frobenius norm of the J † t remains bounded at all times provided that ∆ 0 a 2 0 is small enough, a quantity that depends only on the initial conditions. Moreover, it also shows that making the product ∆ 0 a 0 small enough at initialization ensures that w t remains in a ball of radius R around w 0 . We exploit these two consequences to prove Proposition 3 in Appendix B.4.

B.4 Absence of blow-up for almost-optimal initial linear layer: proof of Proposition 3

To prove Proposition 3, we rely on the results of Appendix B.3 which allow controlling the evolution of the smallest singular value of the NTK matrix A wt . More precisely, Proposition 3 is a direct consequence of Proposition 9 with a particular choice for the initialization w 0 .

Proof of Proposition 3 . By assumption on w 0 , it holds that

∆ 0 := ∥∇L(f w0 )∥ < ϵ = µµ H µ GN 8LN min(R, C -1 R ) min σ 0 , σ 2 0 .
Moreover, by definition of a 0 := J † 0 F and of σ 2 0 = σ min (G(u 0 )), we know that a 2 0 < N σ -2 0 . Hence, we get

∆ 0 < 1 2ζ min(R, C -1 R ) min a -1 0 N -1 2 , a -2 0 ,
where we introduced ζ := 4L µµ H µ GN . This above inequality directly yields:

ζ∆ 0 a 0 ≤ R 2N 1 2 , ζC R ∆ 0 a 2 0 < 1, (21)
Therefore, we can use the above inequalities in the estimates [START_REF] Elfwing | Sigmoid-Weighted Linear Units for Neural Network Function Approximation in Reinforcement Learning[END_REF] of Proposition 9 to get:

∥w t -w 0 ∥ ≤ R 2N 1 2 , a t ≤ a 0 1 -ζC R ∆ 0 a 2 0 -1 2 ,
for any time t ∈ [0, T R ), where T R is defined as:

T R := sup {t ≥ 0 | ∥w t -w 0 ∥ < R, and σ ⋆ (A wt ) > 0}.
This implies, in particular, that w t never escapes the ball B(w 0 , R). Moreover, using [START_REF] Deng | The mnist database of handwritten digit images for machine learning research[END_REF] we get that the smallest singular value σ ⋆ (A ⊤ wt ) is always lower-bounded by a positive constant σ := a -2 0 1 -ζC R ∆ 0 a 2 0 . Therefore, T R must necessarily be greater or equal to T , the maximal time over which w t is defined. Additionally, if T was finite, then w t must escape any bounded domain. This contradicts the fact that w t ∈ B(w 0 , R). Therefore, T must be infinite and the inequality (6) applies at all times t ≥ 0.

B.5 Auxiliary results

Lemma 1. Let J be a linear operator from a Hilbert space W to a finite-dimensional space H and assume that JJ ⊤ is invertible. Define the projector P = J ⊤ (JJ ⊤ ) -1 J. Then the following relations hold:

(J † ) ⊤ (I -P ) = 0 Lemma 2. Let J be a linear operator from a Hilbert space W to a finite-dimensional space H and assume that JJ ⊤ is invertible. Then for any non-negative number ϵ, it holds that:

JJ ⊤ (JJ ⊤ + ϵI) -2 F ≤ J † 2 F Lemma 3.
Let W be a Hilbert space and R be a finite-dimensional euclidean space. Let J be a linear operator from W to R and H an invertible operator from R to itself. Further, assume that JJ ⊤ is invertible. Then the following holds for any non-negative number ϵ:

H -1 2 JJ ⊤ + ϵH -1 -1 H -1 2 op ≤ 1 ϵ + σ ⋆ (H)σ ⋆ (JJ ⊤ ) .
Lemma 4. Let J and K be two Hilbert-Schmidt operators between two Hilbert spaces. Then, it holds that:

JK ⊤ F ≤ ∥J∥ op ∥K∥ F ≤ ∥J∥ F ∥K∥ F Lemma 5.
Let C and r be two positive constants. For a given initial condition b 0 > 0, consider the following differential equation: ḃt = Cb 3 t e -rt . Then b t is defined for any time t < T + where T + is given by:

T + = r -1 log 1 -r 2Cb 2 0 , 2Cr -1 b 2 0 > 0 +∞ Otherwise.
Moreover, b t is given by:

b t = b 0 1 + 2Cr -1 b 2 0 e -rt -1 -1 2 , t < T + .
Proof. First note that a local solution exists by the Cauchy-Lipschitz theorem. Moreover, it must never vanish since otherwise, it would coincide with the null solution b t = 0 by uniqueness. However, this is impossible, since b 0 > 0. We can therefore divide the ODE by b 3 t and explicitly integrate it which gives: b -2 t = b -2 0 + 2r -1 C e -rt -1 . The solution is defined only for times t so that the r.h.s. is positive, which is exactly equivalent to having t < T + .

C Additional experimental results

C.1 Multi-seed experiments

We present result obtained for 5 independent runs. Each run uses a different initialization for the parameters of both student and teacher networks. Additionally, the training and test data are all generated independently for each run. All these results are obtained using SiLU non-linearity. It is clear from Figure 3, that the results display little variability w.r.t. the seed of the experiment, which justifies the single-seed setting chosen in the main paper.

Effect of over-parameterization Figure 4 (Left) shows the effect of over-parameterization of the network on generalization error. It appears that the generalization error remains stable with an increasing parameterization for both GN and GD, as soon as the network has enough overparameterization to fit the data exactly. On the other hand, the test error of RF improves with increasing over-parameterization which is consistent with [START_REF] Hastie | Surprises in high-dimensional ridgeless least squares interpolation[END_REF].

Evolution of the smallest singular value of A w . Figure 4 (Right) shows that the smallest singular value σ ⋆ (A w ) systematically increases during training, especially in the feature learning regime (bright red colors). Hence, the dynamics are far from blowing up even as the features change substantially. On the other hand, σ ⋆ (A w ) remains nearly constant for large values of τ 0 (darker colors) which indicates that the features barely change in the kernel regime. In both cases, we display the results for two different values for the initial std τ 0 (10 -3 and 1).

Results are reported for the best performing step size for each of GN and GD. The first observation is that the performance gap between random features RF and learned features (using either GN or GD) keeps increasing with the sample size. This observation suggests that the learned features use training samples more efficiently than random features. The second observation is the almost affine relation between the generalization error L test and training sample size in a logarithmic scale with a strong dependence of the slope on the variance at initialization and the optimization method. Such affine relation implies the following upper-bound on L test in terms of sample-size N :

L test ≤ C 1 N α
, for some positive constants C and α. The coefficient α controls the speed at which the estimator approximates f ⋆ as one accesses more training data and usually appears in learning theory to describe the statistical convergence rate of a given estimator [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]. The consistent and strong dependence of α on the optimization methods (GN vs GD) and regime (kernel vs feature learning) further confirms the implicit bias of both initialization and optimization algorithms. All models are optimized until the training objective is smaller than 10 -6 using M = 5000 hidden units. For both GD and GN, results are reported for the best-performing step-size λ selected according to the test loss on a regular logarithmic grid ranging from 10 -3 to 10 3 .

C.3 Choice of the non-linearity

Figure 6 shows the effect of step-size on the final performance of a SiLU network with two choices for the temperature parameter: β = 1 and β = 10 6 . Similarly to Figure 1, it is clear that increasing the step size in Gauss-Newton results in features that do not generalize well, while the opposite is observed in case of gradient descent. On the other hand, depending on the value of β, we observe that GN either outperforms GD (β = 10 6 ) or the opposite (β = 1). This variability suggests that the effectiveness of an optimization method in improving generalization is dependent on the specific characteristics of the problem at hand. Finally, note that for β = 10 6 , we recover almost the same results as those obtained for ReLU in Figure 1 (right). This is due to the fact SiLU becomes a tight smooth approximation to ReLU for large values of β. We also notice that when using a larger std for the initial weights (ex. τ 0 = 1), using either ReLU or SiLU with (β = 1) yields similar results, as shown in Figure 7.

C.4 Experiments on MNIST

We have performed a series of experiments on MNIST to illustrate the behavior of GN on higher dimensional data (notably, with multidimensional targets). To stay as close as possible to the theoretical framework, we modify the original MNIST dataset. Specifically, we construct a student/teacher setup in which the trained network can achieve a zero-loss, which is an assumption made in Section 3.1.

Building datasets based on MNIST by using a teacher. In our theoretical framework, we restricted ourselves to cases where the minimum loss can be achieved on the considered model. To make this condition hold, we build a new dataset by using a teacher trained on the MNIST classification task [START_REF] Deng | The mnist database of handwritten digit images for machine learning research[END_REF]:

1. we construct a training set for the teacher which consists of D tr ∪D ts , where D ts is the original MNIST test dataset and D tr is a balanced subset of size 3000 of the original MNIST training set;

2. we then train a teacher neural network with one hidden layer of size 50 and ReLU activation function on training set D tr ∪ D ts using a cross-entropy objective; Step size Step size Step size Step size Training procedure. Once the pair of training and testing datasets ( Dtr , Dts ) has been built, we use them to train and test a student neural network based on [START_REF] Boursier | Gradient flow dynamics of shallow ReLU networks for square loss and orthogonal inputs[END_REF], with one hidden layer of size 5000 and a ReLU activation function. We use a similar initialization scheme for the student's network and fix the std of the initial weights to τ 0 = 0.001. For training, we use a quadratic loss between the student and teacher's outputs and perform a maximum of 100000 iterations of either GN or GD. The training stops whenever the training loss goes below 10 -6 . We vary the step size λ on a logarithmic range [10 -3 , 10 3 ] for both methods and evaluate the test loss and test loss after linear re-fitting.

23

Figure 1 :

 1 Figure 1: Final values of various metrics vs std of the hidden layer's weights at initialization τ 0 (left) and step size (right) for a ReLU network. (Left figure)The training size is set to N = 500 while τ 0 ranges from 10 -3 to 10 2 . For both GD and GN, results are reported for the best-performing step-size λ selected according to the test loss on a regular logarithmic grid ranging from 10 -3 to 10 3 . (Right figure) The std of the weights at initialization is set to τ 0 = 10 -3 . All models are optimized up to a training error of 10 -6 or until the maximum number of steps is exceeded, (M = 5000 , N = 500).

Figure 2 :

 2 Figure2: Evolution of various metrics during training for both GD and GN. Both the variance of the weights at initialization τ 0 and step-size λ are selected from two sets {10 -3 , 1, 10 2 } and {10 -2 , 1, 10 2 }. In all cases, we report results for the best-performing choice of τ 0 (achieved for τ 0 = 10 -3 ) based on test error after linear re-fitting (Test-LRfit). Both Top-Left and Bottom-Left figures show the evolution of the WCD defined in (9) and the three metrics (the training loss, the test loss and Test-LRfit) for both GD and GN using the step-size λ achieving the lowest Test-LRfit. Top-Right and Bottom-Right figures show the evolution of WCD, train loss and test loss after linear re-fitting for the various choices of step-size (darker colors correspond to larger step-sizes).

  where the second line follows by the Woodbury matrix identity, and the third line follows by simple linear algebra. For the fourth line, we use the properties of the Frobenius norm. The fifth and last lines are direct consequences of Lemmas 2 and 4 in Appendix B.5. Concluding. We can combine the upper-bound on | ȧt | and ∥ ẇt ∥ to get:

Proposition 9 .

 9 Consider the same conditions as in Proposition 8. Let ζ := 4L µµ H µ GN and consider the (possibly infinite) time T + defined as:

Figure 3 :Figure 4 :Figure 5

 345 Figure 3: (Left) Test error vs std of the initial weights τ 0 . All models are trained until the training objective is smaller than 10 -6 using M = 5000 hidden units and N = 500. Confidence interval at 95% estimated from 5 independent runs shown in shaded colors. (Right) Test error vs training sample size. All models are trained until the training objective is smaller than 10 -6 using M = 5000 hidden units and initial std of the weights τ 0 = 1. Confidence interval at 95% estimated from 5 independent runs shown in shaded colors.

Figure 5 :

 5 Figure 5: Final values of various metrics vs training sample size for both ReLU (left) and SiLU (right) networks.Two values for τ 0 are displayed for each method τ 0 ∈ {10 -3 , 1} while training data ranges from N = 100 to N = 1500. All models are optimized until the training objective is smaller than 10 -6 using M = 5000 hidden units. For both GD and GN, results are reported for the best-performing step-size λ selected according to the test loss on a regular logarithmic grid ranging from 10 -3 to 10 3 .

Figure 6 :

 6 Figure 6: Final values of various metrics vs the step size for both SiLU network with β = 1 (left) and β = 10 6 (right). (Right figure) The std of the weights at initialization is set to τ 0 = 10 -3 . All models are optimized up to a training error of 10 -6 or until the maximum number of steps is exceeded, (M = 5000 , N = 500).

Figure 7 :

 7 Figure 7: Final values of various metrics vs the step size for both ReLU (left) and SiLU (right) networks. The std of the weights at initialization is set to τ 0 = 1. All models are optimized up to a training error of 10 -6 or until the maximum number of steps is exceeded, (M = 5000 , N = 500).

3 .

 3 given the trained teacher f T , we build new datasets Dtr and Dts for the student network:Dtr := {(x i , Softmax(f T (x i ))) : (x i , y i ) ∈ D tr }, Dts := {(x i , Softmax(f T (x i ))) : (x i , y i ) ∈ D ts }.

The step size can be taken as large as possible provided the algorithm does not diverge. In our experiments, we found that going beyond a step size of λ = 10 3 for GD results in a divergence of the algorithm

Acknowledgments

The project leading to this work has received funding from the French National Research Agency (ANR-21-JSTM-0001) in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02). This work was granted access to the HPC resources of IDRIS under the allocation 2023-AD011013762R1 made by GENCI.

Step size Results. Figure 8 shows the evolution of the test loss (both with or without re-fitting) as a function of the step size. We observe similar behavior as in the main experiments, with hidden learning occurring for GN when using small step sizes and a degraded generalization error as the step size increases. On the other hand, generalization improves for GD by increasing the step size up to the point where optimization becomes unstable.


Nicolas Comte

Sergi Pujades

Aurélien Courvoisier

Olivier Daniel

Jean-Sébastien Franco

François Faure

Edmond Boyer

3D inference of the scoliotic spine from depth maps of the back

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Adolescent idiopathic scoliosis (AIS) is a progressive disease mostly affecting young people and evolving during the period of growth [START_REF] Weinstein | Adolescent idiopathic scoliosis[END_REF]. The main recommendation is to use braces to stop scoliosis progression with medical monitoring until the end of the adolescence. Diagnosis and follow-up are usually performed with radiographs by measuring the Cobb Angle at the scoliotic curvature on the coronal plane. This method is highly accurate but it involves multiple radiation exposures which increase the risk of cancer [START_REF] Hoffman | Breast Cancer in Women With Scoliosis Exposed to Multiple Diagnostic X Rays[END_REF] [START_REF] Ronckers | Multiple Diagnostic X-rays for Spine Deformities and Risk of Breast Cancer[END_REF]. Therefore, there is a strong interest to examine spines using non-invasive and non-ionising methods, especially with young people [START_REF] Knott | SOSORT 2012 consensus paper: reducing x-ray exposure in pediatric patients with scoliosis[END_REF].

In the recent years, deep learning methods based on Convolutional Neural Networks (CNN) showed promising results into the prediction of the Cobb angle from outer images of the back. Yang et al 2019 [START_REF] Yang | Development and validation of deep learning algorithms for scoliosis screening using back images[END_REF] trained CNN models to classify RGB images of the back according to the scoliosis severity. Kokabu et al 2021 [START_REF] Kokabu | An algorithm for using deep learning convolutional neural networks with three dimensional depth sensor imaging in scoliosis detection[END_REF] extended the Adam's forward bend test by the inference of the severity of scoliosis, i.e. the main Cobb angle, from depth maps. Watanabe et al. 2019 [START_REF] Watanabe | An application of artificial intelligence to diagnostic imaging of spine disease: Estimating spinal alignment from moiré images[END_REF] predicted the coronal spinal alignments from Moiré images and estimate the cobb angles along the predictions. However, these methods are providing a limited 2D characterization of the deformities while scoliosis is inherently a three-dimensional condition [START_REF] Stokes | Three-dimensional terminology of spinal deformity. A report presented to the Scoliosis Research Society by the Scoliosis Research Society Working Group on 3-D terminology of spinal deformity[END_REF].

Our hypothesis is that the information present in a depth map of the back of a person has a very strong correlation with the underlying 3D shape of the spine, i.e. the 3D location of each vertebra. Thus, a model taking as input a depth map could predict the 3D shape of the spine. The goal of this study is to evaluate such a strategy and its accuracy by comparing the predictions made by the proposed CNN model with real-life measurements of spine positions and curvatures.

Data collection and processing

To validate our hypothesis and approach we use two different datasets. One acquired at Grenoble Hospital (GH) and the publicly available NMDID dataset [START_REF] Edgar | New Mexico Decedent Image Database[END_REF].

GH dataset

The study was approved by two Ethical Committees. The GH dataset was collected and published in Courvoisier et al 2019 [START_REF] Courvoisier | Prediction of brace effect in scoliotic patients: blinded evaluation of a novel brace simulator -an observational cross-sectional study[END_REF]: CECIC Rhône-Alpes-Auvergne, Clermont-Ferrand, IRB 5891. The IRB for recent cases is CPP Ile de France 2 on the 07/20/2020: n• ID RCB: 2020-A01071-38. All parents and patients received an information letter.

Two types of data were collected from 32 patients with AIS (7-16 y. old, 81% females): a 3D surface scan of the back or the torso with an optical scan (ScanGogh II Vorum Research Corporation or the Occipital Structure Sensor Mark II), and the corresponding 3D spine reconstruction obtained from biplanar X-rays with the EOS Imaging systems [START_REF] Humbert | 3D reconstruction of the spine from biplanar x-rays using parametric models based on transversal and longitudinal inferences[END_REF].

As the surface scan and the EOS images were not performed simultaneously, the patient pose can vary between both acquisitions. We therefore added an additional step to adjust the skin surface pose to match the EOS pose. We use a 3D kinematic anatomic model in Sofa [START_REF] Dicko | Anatomy Transfer[END_REF] which we personalize to the back surface. We then deform the pose of the avatar to match the visual envelope extracted from the biplanar EOS radiographs. The obtained avatar has thus the skin surface in the same pose as the pose observed during the radiographs.

NMDID dataset.

The GH dataset is composed of AIS patients only. To explore asymptomatic AIS cases we included subjects from the NMDID [START_REF] Edgar | New Mexico Decedent Image Database[END_REF] database, which provides full-body CT-scans of ex vivo subjects. We included 89 cases (15-30 y. old, 32% females, 7% with scoliosis), which had 17 thoracolumbar vertebrae and a Body Mass Index (BMI) under 30kg/m 2 . A selection on BMI was needed since the corpses are in supine position, and the back is flatten during the CT-acquisition. High BMI bodies had completely flat backs with no shape information.

From the selected cases, we automatically segmented the thoracolumbar spine vertebrae using Meng et al. 2023 method [START_REF] Meng | Vertebrae localization, segmentation and identification using a graph optimization and an anatomic consistency cycle[END_REF] and the skin with a binary segmentation approach on each slice of the volumetric image. The results were then turned into meshes using the marching-cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF].

Depth map generation.

From the meshes of the backs of the patients we rendered depth maps using an orthographic camera model. For the NMDID subjects we directly used the segmented skin, whereas for the GH subjects, we used the skin of the repositionned avatar matching the 3D spine shape in the radiographs.

Spine 3D characteristics.

From the 3D vertebra model positions obtained from the EOS Imaging system [START_REF] Humbert | 3D reconstruction of the spine from biplanar x-rays using parametric models based on transversal and longitudinal inferences[END_REF] and the automatic CT segmentation [START_REF] Meng | Vertebrae localization, segmentation and identification using a graph optimization and an anatomic consistency cycle[END_REF], we computed the 3D spine characteristics.

Using the method by Choi et al. [START_REF] Choi | CNN-based spine and Cobb angle estimator using moire images[END_REF] based on cubic B-Splines, we computed the scoliosis severity by selecting the maximum absolute Cobb angle value. Kyphosis and lordosis were obtained with splines fitted in the sagittal plane between T04-T12 and L01-L05 respectively. As we do not model S01 we computed the lordosis Cobb angle using the L05 centroid. From the NMDID and GH datasets we reconstructed 121 pairs of 3D spines and surfaces of the back: 89 from the NMDID and 32 from GH. The input dataset we consider is therefore composed of 45% females aged between 7 and 30 years old. 38 cases with scoliosis are included in our analysis with different type of curvatures (simple, double curvatures, 10-65 • ). Fig. 1 illustrates one case of each dataset with the depth map and the 3D characterization of the associated spine.

Method

Our method takes as input a depth map of the back of a patient and outputs the 3D location of the thoraco-lumbar vertebrae. To ease the spine prediction we use a low dimensional representation of the spine based on Principal Component Analysis (PCA). Fig. 2 presents the overview of the method.

PCA spine representation.

To ensure that the predicted individual vertebrae create a consistent spine shape, we use a low dimensional representation based on PCA (Principal Component Analysis) [START_REF] Pearson | LIII on lines and planes of closest fit to systems of points in space[END_REF]. The learned PCA space can produce a compressed latent representation of the spine data, while preserving most of its information and variability. Two main advantages follow. The prediction of the n θ = 20 PCA coefficients is less complex than the prediction of each of the 17 individual vertebrae 3D locations (17 × 3 = 51). In addition, the compressed latent representation enforces regularization over the estimated spine locations.

Architecture.

We use ResNet-18 [START_REF] He | Deep residual learning for image recognition[END_REF] model as our backbone CNN architecture, which has been shown to be efficient for regression tasks. It encodes the depth image (224×224, 1 channel) into the aforementioned latent PCA representation of the spine, which is in turn easily decoded into the 3D vertebrae positions. We train the network to minimize an L2-loss between the vertebrae coordinates predicted by the neural model, and the annotated coordinates of the training examples. The model is trained using an Adam optimizer with PyTorch [START_REF] Paszke | Pytorch: An imperative style, highperformance deep learning library[END_REF] with 2000 epochs and a learning rate of 1e-4.

Data Processing.

To structure the dataset and ease the learning task we process the data. Following Choi et al. [START_REF] Choi | CNN-based spine and Cobb angle estimator using moire images[END_REF] the depth maps on the spine are centered and cropped according to its length, then resized to a 224 × 224 resolution. Depth pixel observations are further normalized between [-1, 1], and the antero-posterior positions of the vertebrae are defined with respect to a median vertebra (T08). This 3D transformation allows a simpler representation of the spine while keeping the vertebrae alignment information. It differs from [START_REF] Nerot | Estimation of spinal joint centers from external back profile and anatomical landmarks[END_REF] and [START_REF] Caturano | Prediction of the 3d spinal alignment from external shape of the back in ais patients using regression model[END_REF] which are regressing absolute positions.

We augment the dataset applying a random set of different transformations, such as mirroring and rotating the torso in 3D. From the 121 rendered depth maps, the data augmentation raises the amount of images to 10,890.

Evaluation

We start by validating the accuracy of the PCA spine representation. Then we evaluate the vertebrae location accuracy by computing the 3D distances between the ground-truth and the predictions as well as 2D distances in the coronal and sagittal planes. We evaluate the spine curvature accuracy by computing the Mean Absolute Error (MAE) in three angles: the main Cobb angle [START_REF] Choi | CNN-based spine and Cobb angle estimator using moire images[END_REF], kyphosis (T04-T12) and lordosis (L01-L05). As the 3D predicted spine of our method can also be used for spine classification (main Cobb angle ≥ 10 • ) we compute binary classification metrics: sensitivity, specificity, positive predictive value, accuracy and AUC. All our quantitative results are computed using a stratified 20-fold cross-validation (CV) on the full dataset. All metrics are computed by averaging over all test sets.

PCA dimensional reduction.

A first analysis is dedicated to the PCA and how well the resulting representation is able to encode the vertebra positions. For each fold, the PCA model is learned from the training set and evaluated by comparing the PCA reconstructions of the spines of the test set with their true X-rays reconstructions. Our analysis shows that with 20 components, the spines can be reconstructed with a Mean Absolute Error (MAE) in 3D vertebra location of 0.92mm (std: 0.75 mm) and an MAE in scoliosis severity of 2.71deg (std: 2.49). These values are close to the annotation accuracy (< 1 mm and < 5 • ). 1. Our model is able to predict the 3D vertebra positions with an average 3D error below the cm. These predictions allows the characterization of the spinal alignments (scoliosis, kyphosis, lordosis) with automated measurements [START_REF] Choi | CNN-based spine and Cobb angle estimator using moire images[END_REF]. Visuals of the predictions are presented in Fig. 3.

MAE of the vertebrae location predictions are reported in Table

Spine 3D characteristics.

Table 1 reports the angle errors between the predicted spine shapes and the ground truth ones. The last three lines of the Table 3) are provided for the comparison with the literature. Note how our approach is the only one predicting kyphosis and lordosis. Numerical results are provided for informative comparison on the range of values, as the considered datasets are different, both in size and population.

Classification.

Classification results are reported in Table 1. Our method is able to discriminate cases with scoliosis with a sensitivity of 64% and a specificity of 99%. We report a positive predictive value of 95%, an accuracy of 89% and an AUC of 90%. Fig. 4 shows the predicted severity against the ground truth one with the 10 • classification thresholds. Let us note how several miss-classifications arise near the 10°threshold and the underestimation of moderate and severe scoliosis cases. In conclusions we provide leads to improve these predictions.

Table 1: Spine prediction accuracy. Vertebrae errors in mm of the 3D locations (3D) and 2D coronal (Cor.) and sagittal (Sag.) projections. Angle errors in degrees in the scoliosis severity (Sco. sev.), kyphosis (Kyph.) and lordosis (Lord.). Scoliosis classification (≥ 10 • ) sensitivity (Sens.) and specificity (Spec.).

Locations (mm)

Characteristics (

Ablation studies

To further understand the proposed approach we performed two complementary experiments.

Recent work by Klarqvist et al. 2022 [START_REF] Klarqvist | Silhouette images enable estimation of body fat distribution and associated cardiometabolic risk[END_REF] has shown that it is possible to estimate body composition from solely body silhouettes. We consider the analogous case for the prediction of the 3D spine: instead to use a depth map as input we consider the case where only the binary silhouette of the back surface is observed.

In addition, we also experiment by training and testing on the two different datasets (GH and NMDID) to see if the imbalance of scoliotic patients in the datasets has an effect on the predictions. We tested the impact of depth information by training our model to reconstruct the 3D positions on binary images that show only the silhouette of the torso. The results in Table 2 illustrate that model using depth maps achieves overall better performance. Depth is informative for the network.

It is interesting to note the capability of the sihouette network to reconstruct scoliosis curvatures with a severity error of 6.34 • . The biggest difference is obtained in the sagittal plane alignments, which is coherent with the lack of depth information in the binary silhouette image.

Different datasets.

The inclusion of NMDID allows us to validate the approach on asymptomatic cases. However, it could have a negative impact on the prediction of the scoliotic cases by adding bias. We thus evaluate the impact of the NMDID dataset on the network regression performance. We train solely on GH (scoliotic data) and on both GH and NMDID (scoliotic + asymptomatic) data and report the results on Table 3. Including the NMDID cases in the training set does not introduce any bias in the prediction performances (first two lines) on GH subjects and confirms that the inclusion of spines without scoliosis from deceased subjects in supine position does not deteriorate the predictions over living patients acquired with the regular protocol.

Conclusion

In this work, we presented an approach to predict the 3D spine shape of a patient from a depth map of the back. This depth map can be obtained with depth sensors, which are usually inexpensive, portable, making them a promising cost-effective, non-ionizing approach to quantify scoliotic deformities. Our approach provides accurate 3D predictions of the vertebrae locations and allows to compute relevant anatomic curvatures along the spine as well as a fast scoliosis diagnosis.

Our presented work predicts 3D vertebrae locations, but does not consider today their individual 3D orientation. It has been studied that the axial vertebra rotation contains information related to the early detection of scoliosis [START_REF] Courvoisier | Transverse plane 3D analysis of mild scoliosis[END_REF] and its evolution [START_REF] Skalli | Early Detection of Progressive Adolescent Idiopathic Scoliosis: A Severity Index[END_REF]. Our approach could be adapted in the future to also predict the orientations of the individual vertebrae. To that end, the PCA representation should be reconsidered, as a linear approach might not be well suited to capture the 3D rotations.

Another venue of improvement of the presented work would be to add more subjects with moderate and severe scoliosis in the training pipeline. As shown in Fig. 4, our model tends to underestimate severe deviations due to the low number of such cases in our dataset. To make this possible, as well as to foster future research in this direction, we make our trained model and training code available for research purposes: https://gitlab.inria.fr/spine/skin_to_spine.

Aknowledgments

This study received support from the ANRT (Association Nationale de la Recherche et de la Technologie) CIFRE No.2019/1197 and the ANR (Agence Nationale de la Recherche) JCJC SEMBA ANR-19-CE23-0003.

Figure 1 :

 1 Figure 1: Dataset examples. For each subject we have a depth map, the 3D vertebrae locations and the measured angles. Coronal (left) and sagittal (right) views.

Figure 2 :

 2 Figure 2: An input depth map D i is fed into the Resnet that predicts the latent PCA representation θ i . From θ i the 3D vertebrae coordinates Θ i are computed.

Figure 3 :

 3 Figure 3: Selection of results. White: ground-truth; blue: predictions. Top row: coronal view for each case, bottom row: sagittal view. GH: institution; NMDID: New Mexico Decedent Image Database; L: Left; R: Right; P: Posterior; A: Anterior.

4. 2

 2 3D spine prediction accuracy 4.2.1 Vertebra locations.

Figure 4 :

 4 Figure 4: Comparison between predicted scoliosis severity and ground-truth measurements.

Table 2 :

 2 Comparison with binarized depth information, i.e. silhouette (silh). MAE (with std.) on locations in 3D, coronal (Cor.) and sagittal (Sag.) plane and the estimated scoliosis severity. Detection is evaluated sensitivity (Sens). and specificity (Spec.).

	Type	3D (mm)	Cor. (mm) Sag. (mm) Severity (°) Sens. Spec.
	silh	10.02 (6.96) 4.91 (3.34) 8.55 (7.22) 6.34 (7.86)	0.52	0.98
	depth 7.07 (4.69)	4.51 (2.99) 5.60 (4.62) 5.46 (6.19)	0.64	0.99

Table 3 :

 3 MAE in locations in mm and standard-deviations

	Study Train	Test	3D	Coronal	Sagittal
	Ours	GH	GH	8.81 (5.19) 5.87 (3.73) 6.76 (4.97)
	Ours	GH+NMDID GH	8.80 (5.46) 5.90 (3.73) 6.74 (5.26)
	Ours	GH+NMDID NMDID	6.45 (4.21) 4.01 (2.49) 5.19 (4.30)
	Ours	GH+NMDID GH+NMDID 7.07 (4.69) 4.51 (2.99) 5.60 (4.62)
	4.3.1 Predictions from binary images.